Skip to content

[PWGJE] MC process functions updated #9809

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
nzardosh merged 7 commits into from
Feb 8, 2025
Merged

[PWGJE] MC process functions updated #9809

Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
308828f
Change particle level. Added pure det lvl process function. Added cen…
joachimckh Jan 29, 2025
02627a1
clang-format
joachimckh Jan 29, 2025
6fa1deb
modification to MC process functions.. centrality to MCP
joachimckh Feb 6, 2025
2e49ee6
readable event histograms
joachimckh Feb 7, 2025
09c8049
Merge branch 'master' into jet-unf
joachimckh Feb 7, 2025
9c69cd3
format
joachimckh Feb 7, 2025
ba7817d
fix err
joachimckh Feb 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 62 additions & 20 deletions PWGJE/Tasks/jetSpectraEseTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ struct JetSpectraEseTask {
ConfigurableAxis binESE{"binESE", {100, 0, 100}, ""};
ConfigurableAxis binCos{"binCos", {100, -1.05, 1.05}, ""};
ConfigurableAxis binOccupancy{"binOccupancy", {5000, 0, 25000}, ""};
ConfigurableAxis binQVec{"binQVec", {100, -3, 3}, ""};
ConfigurableAxis binQVec{"binQVec", {500, -3, 3}, ""};

Configurable<float> jetPtMin{"jetPtMin", 5.0, "minimum jet pT cut"};
Configurable<float> jetR{"jetR", 0.2, "jet resolution parameter"};
Expand All @@ -60,7 +60,8 @@ struct JetSpectraEseTask {
Configurable<std::string> eventSelections{"eventSelections", "sel8FullPbPb", "choose event selection"};
Configurable<std::string> trackSelections{"trackSelections", "globalTracks", "set track selections"};

Configurable<bool> cfgEvSelOccupancy{"cfgEvSelOccupancy", false, "Flag for occupancy cut"};
Configurable<bool> cfgEvSelOccupancy{"cfgEvSelOccupancy", true, "Flag for occupancy cut"};

Configurable<std::vector<int>> cfgCutOccupancy{"cfgCutOccupancy", {0, 1000}, "Occupancy cut"};
Configurable<std::vector<float>> cfgOccupancyPtCut{"cfgOccupancyPtCut", {0, 100}, "pT cut"};

Expand All @@ -83,6 +84,9 @@ struct JetSpectraEseTask {
std::vector<int> eventSelectionBits;
int trackSelection{-1};

Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f) && nabs(aod::jet::eta) < 0.9f - jetR;
Filter colFilter = nabs(aod::jcollision::posZ) < vertexZCut;
Filter mcCollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut;
using ChargedMCDJets = soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>>;
Preslice<ChargedMCDJets> mcdjetsPerJCollision = o2::aod::jet::collisionId;

Expand All @@ -104,7 +108,8 @@ struct JetSpectraEseTask {
if (doprocessESEDataCharged) {
LOGF(info, "JetSpectraEseTask::init() - processESEDataCharged");
registry.add("hEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
registry.add("hCentralityMult", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
registry.add("hCentralitySel", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
registry.add("hCentralityAnalyzed", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
registry.add("hJetPt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
registry.add("hJetPt_bkgsub", "jet pT background sub;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
registry.add("hJetEta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
Expand All @@ -128,30 +133,50 @@ struct JetSpectraEseTask {
registry.addClone("hPsi2FT0C", "hEPUncorV2");
registry.addClone("hPsi2FT0C", "hEPRectrV2");
registry.addClone("hPsi2FT0C", "hEPTwistV2");

registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(2, "Event selection");
registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(3, "Occupancy cut");
registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(4, "ESE available");
registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(5, "Lead track");
registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(6, "Jet loop");
registry.get<TH1>(HIST("hEventCounter"))->GetXaxis()->SetBinLabel(7, "Centrality analyzed");
}
if (doprocessMCParticleLevel) {
LOGF(info, "JetSpectraEseTask::init() - processMCParticleLevel");
registry.add("hMCPartEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
registry.add("hPartCentralitySel", ";centr;entries", {HistType::kTH1F, {{100, 0.0, 100.0}}});
registry.add("hPartJetPt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
registry.add("hPartJetPtSubBkg", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
registry.add("hPartJetSparse", ";Centrality;#it{p}_{T,jet part} (GeV/#it{c})", {HistType::kTHnSparseF, {{110, -10, 100}, {jetPtAxis}, {100, -1.0, 1.0}, {80, -1.0, 7}}});
registry.add("hPartJetEta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
registry.add("hPartJetPhi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}});

registry.get<TH1>(HIST("hMCPartEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
registry.get<TH1>(HIST("hMCPartEventCounter"))->GetXaxis()->SetBinLabel(2, "Collision size < 1");
registry.get<TH1>(HIST("hMCPartEventCounter"))->GetXaxis()->SetBinLabel(3, "MCD size != 1");
}
if (doprocessMCDetectorLevel) {
LOGF(info, "JetSpectraEseTask::init() - processMCDetectorLevel");
registry.add("hMCDetEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
registry.add("hDetCentralitySel", ";centr;entries", {HistType::kTH1F, {{100, 0.0, 100.0}}});
registry.add("hDetJetPt", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
registry.add("hDetJetPtSubBkg", "particle level jet pT;#it{p}_{T,jet part} (GeV/#it{c});entries", {HistType::kTH1F, {{jetPtAxis}}});
registry.add("hDetJetSparse", ";Centr;#it{p}_{T,jet part} (GeV/#it{c})", {HistType::kTHnSparseF, {{110, -10, 100}, {jetPtAxis}, {100, -1.0, 1.0}, {80, -1.0, 7}}});
registry.add("hDetJetEta", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
registry.add("hDetJetPhi", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}});

registry.get<TH1>(HIST("hMCDetEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
registry.get<TH1>(HIST("hMCDetEventCounter"))->GetXaxis()->SetBinLabel(2, "Event eelection");
registry.get<TH1>(HIST("hMCDetEventCounter"))->GetXaxis()->SetBinLabel(3, "Occupancy cut");
}
if (doprocessMCChargedMatched) {
LOGF(info, "JetSpectraEseTask::init() - processMCChargedMatched");
registry.add("hMCEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
registry.add("hMCDMatchedEventCounter", "event status;event status;entries", {HistType::kTH1F, {{10, 0.0, 10.0}}});
registry.add("hCentralityMult", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
registry.add("hCentralityAnalyzed", ";Centrality;entries", {HistType::kTH1F, {{100, 0, 100}}});
registry.add("hPartJetPtMatch", ";Centrality;#it{p}_{T,jet part} (GeV/#it{c})", {HistType::kTH2F, {{100, 0, 100}, {jetPtAxis}}});
registry.add("hPartJetPtMatchSubBkg", ";Centrality;#it{p}_{T,jet part} (GeV/#it{c})", {HistType::kTH2F, {{100, 0, 100}, {jetPtAxis}}});
registry.add("hPartJetMatchSparse", ";Centrality;#it{p}_{T,jet part} (GeV/#it{c})", {HistType::kTHnSparseF, {{110, -10, 100}, {jetPtAxis}, {100, -1.0, 1.0}, {80, -1.0, 7}}});
registry.add("hPartJetEtaMatch", "particle level jet #eta;#eta_{jet part};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}});
registry.add("hPartJetPhiMatch", "particle level jet #phi;#phi_{jet part};entries", {HistType::kTH1F, {{80, -1.0, 7.}}});
registry.add("hDetectorJetPt", ";Centrality;#it{p}_{T,jet det} (GeV/#it{c})", {HistType::kTH2F, {{100, 0, 100}, {jetPtAxis}}});
Expand All @@ -165,6 +190,17 @@ struct JetSpectraEseTask {
registry.add("hMatchedJetsPhiDelta", "#phi_{jet part}; det - part", {HistType::kTH2F, {{80, -1.0, 7.}, {200, -10.0, 10.}}});
registry.add("hRespMcDMcPMatch", ";Centrality,#it{p}_{T, jet det}; #it{p}_{T, jet part}", HistType::kTHnSparseF, {{100, 0, 100}, jetPtAxis, jetPtAxis});
registry.add("hRespMcDMcPMatchSubBkg", ";Centrality,#it{p}_{T, jet det}; #it{p}_{T, jet part}", HistType::kTHnSparseF, {{100, 0, 100}, jetPtAxis, jetPtAxis});

registry.get<TH1>(HIST("hMCEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
registry.get<TH1>(HIST("hMCEventCounter"))->GetXaxis()->SetBinLabel(2, "Collision size < 1");
registry.get<TH1>(HIST("hMCEventCounter"))->GetXaxis()->SetBinLabel(3, "Vertex cut");
registry.get<TH1>(HIST("hMCEventCounter"))->GetXaxis()->SetBinLabel(4, "After analysis");
registry.get<TH1>(HIST("hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(1, "Input event");
registry.get<TH1>(HIST("hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(2, "Vertex cut");
registry.get<TH1>(HIST("hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(3, "Event selection");
registry.get<TH1>(HIST("hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(4, "Occupancy cut");
registry.get<TH1>(HIST("hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(5, "Centrality cut1:cut2");
registry.get<TH1>(HIST("hMCDMatchedEventCounter"))->GetXaxis()->SetBinLabel(6, "After analysis");
}
if (doprocessESEOccupancy) {
LOGF(info, "JetSpectraEseTask::init() - processESEOccupancy");
Expand All @@ -177,10 +213,6 @@ struct JetSpectraEseTask {
}
}

Filter jetCuts = aod::jet::pt > jetPtMin&& aod::jet::r == nround(jetR.node() * 100.0f) && nabs(aod::jet::eta) < 0.9f - jetR;
Filter colFilter = nabs(aod::jcollision::posZ) < vertexZCut;
Filter mcCollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut;

void processESEDataCharged(soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::Qvectors, aod::QPercentileFT0Cs>::iterator const& collision,
soa::Filtered<aod::ChargedJets> const& jets,
aod::JetTracks const& tracks)
Expand All @@ -194,6 +226,7 @@ struct JetSpectraEseTask {
if (cfgEvSelOccupancy && !isOccupancyWithin(collision))
return;
registry.fill(HIST("hEventCounter"), counter++);
registry.fill(HIST("hCentralitySel"), collision.centrality());

const auto vPsi2{procEP<true>(collision)};
const auto qPerc{collision.qPERCFT0C()};
Expand All @@ -206,7 +239,7 @@ struct JetSpectraEseTask {

registry.fill(HIST("hEventCounter"), counter++);
registry.fill(HIST("hRho"), collision.rho());
registry.fill(HIST("hCentralityMult"), collision.centrality());
registry.fill(HIST("hCentralityAnalyzed"), collision.centrality());
for (auto const& jet : jets) {
float jetpTbkgsub = jet.pt() - (collision.rho() * jet.area());
registry.fill(HIST("hJetPt"), jet.pt());
Expand Down Expand Up @@ -256,25 +289,33 @@ struct JetSpectraEseTask {
}
PROCESS_SWITCH(JetSpectraEseTask, processESEOccupancy, "process occupancy", false);

void processMCParticleLevel(soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>::iterator const& collision,
void processMCParticleLevel(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCollision,
soa::SmallGroups<aod::JetCollisionsMCD> const& collisions,
soa::Filtered<aod::ChargedMCParticleLevelJets> const& jets)
{
float counter{0.5f};
registry.fill(HIST("hMCPartEventCounter"), counter++);
if (collision.size() < 1) {
if (mcCollision.size() < 1) {
return;
}
registry.fill(HIST("hMCPartEventCounter"), counter++);
if (!(std::abs(collision.posZ()) < vertexZCut)) {
if (collisions.size() != 1) {
return;
}
registry.fill(HIST("hMCPartEventCounter"), counter++);
auto centrality{-1};
for (const auto& col : collisions) {
centrality = col.centrality();
}

registry.fill(HIST("hPartCentralitySel"), centrality);
for (const auto& jet : jets) {
const auto mcPt{jet.pt() - (mcCollision.rho() * jet.area())};
registry.fill(HIST("hPartJetPt"), jet.pt());
registry.fill(HIST("hPartJetPtSubBkg"), jet.pt() - (collision.rho() * jet.area()));
registry.fill(HIST("hPartJetPtSubBkg"), mcPt);
registry.fill(HIST("hPartJetEta"), jet.eta());
registry.fill(HIST("hPartJetPhi"), jet.phi());
registry.fill(HIST("hPartJetSparse"), centrality, mcPt, jet.eta(), jet.phi());
}
}
PROCESS_SWITCH(JetSpectraEseTask, processMCParticleLevel, "jets on particle level MC", false);
Expand All @@ -285,27 +326,26 @@ struct JetSpectraEseTask {
aod::JetParticles const&)
{
float counter{0.5f};
registry.fill(HIST("hMCDetEventCounter"), counter++);
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
return;
registry.fill(HIST("hMCDetEventCounter"), counter++);

if (cfgEvSelOccupancy && !isOccupancyWithin(collision))
return;
registry.fill(HIST("hMCDetEventCounter"), counter++);

registry.fill(HIST("hDetCentralitySel"), collision.centrality());
for (const auto& mcdjet : mcdjets) {
auto mcdetPt{mcdjet.pt() - (collision.rho() * mcdjet.area())};
registry.fill(HIST("hDetJetPt"), mcdjet.pt());
registry.fill(HIST("hDetJetPtSubBkg"), mcdetPt);
registry.fill(HIST("hDetJetEta"), mcdjet.eta());
registry.fill(HIST("hDetJetPhi"), mcdjet.phi());
registry.fill(HIST("hDetJetSparse"), collision.centrality(), mcdetPt, mcdjet.eta(), mcdjet.phi());
}
}
PROCESS_SWITCH(JetSpectraEseTask, processMCDetectorLevel, "jets on detector level", false);

using JetMCPTable = soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents, aod::ChargedMCParticleLevelJetsMatchedToChargedMCDetectorLevelJets>>;
void processMCChargedMatched(/*soa::Filtered<aod::JetCollisionsMCD>::iterator const& collision*/
soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>::iterator const& mcCol,
void processMCChargedMatched(soa::Filtered<soa::Join<aod::JetMcCollisions, aod::BkgChargedMcRhos>>::iterator const& mcCol,
soa::SmallGroups<soa::Join<aod::JetCollisionsMCD, aod::BkgChargedRhos>> const& collisions,
ChargedMCDJets const& mcdjets,
JetMCPTable const&,
Expand All @@ -329,6 +369,7 @@ struct JetSpectraEseTask {
if (!(std::abs(collision.posZ()) < vertexZCut)) {
return;
}
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);

if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits))
return;
Expand All @@ -341,7 +382,7 @@ struct JetSpectraEseTask {
if (collision.centrality() < centRange->at(0) || collision.centrality() > centRange->at(1))
registry.fill(HIST("hMCDMatchedEventCounter"), secCount++);

registry.fill(HIST("hCentralityMult"), collision.centrality());
registry.fill(HIST("hCentralityAnalyzed"), collision.centrality());
auto collmcdJets{mcdjets.sliceBy(mcdjetsPerJCollision, collision.globalIndex())};
for (const auto& mcdjet : collmcdJets) {
auto mcdPt{mcdjet.pt() - (collision.rho() * mcdjet.area())};
Expand All @@ -355,6 +396,7 @@ struct JetSpectraEseTask {
registry.fill(HIST("hPartJetPtMatchSubBkg"), collision.centrality(), mcpPt);
registry.fill(HIST("hPartJetEtaMatch"), mcpjet.eta());
registry.fill(HIST("hPartJetPhiMatch"), mcpjet.phi());
registry.fill(HIST("hPartJetMatchSparse"), collision.centrality(), mcpPt, mcpjet.eta(), mcpjet.phi());
registry.fill(HIST("hMatchedJetsPtDelta"), mcpjet.pt(), mcdjet.pt() - mcpjet.pt());
registry.fill(HIST("hGenMatchedJetsPtDeltadPt"), collision.centrality(), mcpPt, (mcdPt - mcpPt) / mcpPt);
registry.fill(HIST("hRecoMatchedJetsPtDeltadPt"), collision.centrality(), mcdPt, (mcdPt - mcpPt) / mcdPt);
Expand Down
Loading