From ad7b5f4472d9344fd7811e988fe0db206e35a789 Mon Sep 17 00:00:00 2001 From: Changhwan Choi Date: Mon, 11 May 2026 18:49:35 +0900 Subject: [PATCH] Modified GNN b-jet analysis codes to utilise JE outlier rejector task, minimised access to original AOD tables --- PWGJE/Tasks/bjetTaggingGnn.cxx | 421 ++++++++++++++++---------------- PWGJE/Tasks/bjetTreeCreator.cxx | 67 +++-- 2 files changed, 262 insertions(+), 226 deletions(-) diff --git a/PWGJE/Tasks/bjetTaggingGnn.cxx b/PWGJE/Tasks/bjetTaggingGnn.cxx index 8a48eab4af5..e72fba2ccef 100644 --- a/PWGJE/Tasks/bjetTaggingGnn.cxx +++ b/PWGJE/Tasks/bjetTaggingGnn.cxx @@ -128,10 +128,6 @@ struct BjetTaggingGnn { Configurable eventSelections{"eventSelections", "sel8", "choose event selection"}; Configurable useEventWeight{"useEventWeight", true, "Flag whether to scale histograms with the event weight"}; - Configurable pTHatMaxMCD{"pTHatMaxMCD", 999.0, "maximum fraction of hard scattering for jet acceptance in detector MC"}; - Configurable pTHatMaxMCP{"pTHatMaxMCP", 999.0, "maximum fraction of hard scattering for jet acceptance in particle MC"}; - Configurable pTHatExponent{"pTHatExponent", 6.0, "exponent of the event weight for the calculation of pTHat"}; - // track level configurables Configurable trackSelections{"trackSelections", "QualityTracks", "set track selections"}; Configurable trackPtMin{"trackPtMin", 0.15, "minimum track pT"}; @@ -254,7 +250,7 @@ struct BjetTaggingGnn { const AxisSpec axisDbFine{dbNbins, dbMin, dbMax, "#it{D}_{b}"}; const AxisSpec axisJetMass{200, 0., 50., "#it{m}_{jet} (GeV/#it{c}^{2})"}; const AxisSpec axisJetProb{200, 0., 40., "-ln(JP)"}; - const AxisSpec axisNTracks{42, 0, 42, "#it{n}_{tracks}"}; + const AxisSpec axisNTracks{42, 0, 42, "#it{N}_{tracks}"}; registry.add("h_jetpT", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetEta", "", {HistType::kTH1F, {axisJetEta}}, callSumw2); @@ -301,7 +297,8 @@ struct BjetTaggingGnn { registry.add("h_dcaXY_coll_matched_c", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); // tracks from matched collisions, c hadron decay registry.add("h_dcaXY_coll_matched_lf", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); // tracks from matched collisions, others registry.add("h_dcaXY_coll_mismatched", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); // tracks from mismatched collisions - registry.add("h_dcaXY_npp", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); // non-physical primary tracks (GenStatusCode=-1) + registry.add("h_dcaXY_npp", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); // non-physical primary tracks (GenStatusCode=-1) from matched collisions + registry.add("h_dcaXY_npp_mismatched", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); // non-physical primary tracks from mismatched collisions registry.add("h_dcaZ_coll_fake", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); registry.add("h_dcaZ_fake", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); registry.add("h_dcaZ_coll_matched", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); @@ -310,6 +307,7 @@ struct BjetTaggingGnn { registry.add("h_dcaZ_coll_matched_lf", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); registry.add("h_dcaZ_coll_mismatched", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); registry.add("h_dcaZ_npp", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); + registry.add("h_dcaZ_npp_mismatched", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); } if (doprocessDataJetsSel || doprocessMCDJetsSel) { @@ -317,7 +315,6 @@ struct BjetTaggingGnn { registry.add("h_jetpT_selmc", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_tvx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_coll", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_inelgt0rec", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_sel8_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_selmc_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_tvx_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); @@ -372,7 +369,6 @@ struct BjetTaggingGnn { registry.add("h2_Response_DetjetpT_PartjetpT_b_tvx", "b-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); registry.add("h_jetpT_b_coll", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_b_coll_zvtx", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_inelgt0rec", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h2_Response_DetjetpT_PartjetpT_inelgt0", "", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); registry.add("h2_Response_DetjetpT_PartjetpT_b_inelgt0", "b-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); } @@ -388,31 +384,40 @@ struct BjetTaggingGnn { registry.add("h_jetpT_particle_lf", "particle lf-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_sel8", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_sel8", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_selmc", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_selmc", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_tvx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_tvx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_coll", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_coll", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_inel", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_sel8", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_selmc", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_tvx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_coll", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_b_inel", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_sel8", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_selmc", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_tvx", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_coll", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_inel", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_sel8_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_sel8_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_selmc_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_selmc_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_tvx_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_tvx_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_coll_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_coll_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_inel_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_sel8_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_selmc_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_tvx_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_b_coll_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_b_inel_zvtx", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_sel8_zvtx", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_selmc_zvtx", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_tvx_zvtx", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_coll_zvtx", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_inel_zvtx", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_inelgt0", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); registry.add("h_jetpT_particle_b_inelgt0", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_inelgt0rec", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_inelgt0rec", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_jetpT_particle_c_inelgt0", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); } if (doDataDriven) { @@ -446,28 +451,29 @@ struct BjetTaggingGnn { using FilteredCollisions = soa::Filtered; using AnalysisCollisionsTriggered = soa::Join; using FilteredCollisionsTriggered = soa::Filtered; - using OrigCollisions = soa::Join; + // using OrigCollisions = soa::Join; using DataJets = soa::Join; using FilteredDataJets = soa::Filtered; using AnalysisTracks = soa::Join; using FilteredTracks = soa::Filtered; - using OriginalTracks = soa::Join; + // using OriginalTracks = soa::Join; using MCDJets = soa::Join; using FilteredMCDJets = soa::Filtered; using AnalysisTracksMCD = soa::Join; using FilteredTracksMCD = soa::Filtered; - using AnalysisCollisionsMCD = soa::Join; + using AnalysisCollisionsMCD = soa::Join; using FilteredCollisionsMCD = soa::Filtered; Filter mccollisionFilter = nabs(aod::jmccollision::posZ) < vertexZCut; - using FilteredCollisionsMCP = soa::Filtered; + using AnalysisCollisionsMCP = soa::Join; + using FilteredCollisionsMCP = soa::Filtered; using MCPJets = soa::Join; using FilteredMCPJets = soa::Filtered; - template - int analyzeJetTrackInfo(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, AnyOriginalTracks const /*origTracks*/ /*, int8_t jetFlavor = 0, double weight = 1.0*/) + template + int analyzeJetTrackInfo(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/ /*, int8_t jetFlavor = 0, double weight = 1.0*/) { int nTracks = 0; for (const auto& constituent : analysisJet.template tracks_as()) { @@ -477,21 +483,21 @@ struct BjetTaggingGnn { } int sign = jettaggingutilities::getGeoSign(analysisJet, constituent); - auto origConstit = constituent.template track_as(); + // auto origConstit = constituent.template track_as(); registry.fill(HIST("h_gnnfeat_trackpT"), constituent.pt()); - registry.fill(HIST("h_gnnfeat_trackPhi"), origConstit.phi()); + registry.fill(HIST("h_gnnfeat_trackPhi"), constituent.phi()); registry.fill(HIST("h_gnnfeat_trackEta"), constituent.eta()); registry.fill(HIST("h_gnnfeat_trackCharge"), static_cast(constituent.sign())); registry.fill(HIST("h_gnnfeat_trackDCAxy"), std::abs(constituent.dcaXY()) * sign); registry.fill(HIST("h_gnnfeat_trackSigmaDCAxy"), constituent.sigmadcaXY()); registry.fill(HIST("h_gnnfeat_trackDCAz"), std::abs(constituent.dcaZ()) * sign); registry.fill(HIST("h_gnnfeat_trackSigmaDCAz"), constituent.sigmadcaZ()); - registry.fill(HIST("h_gnnfeat_trackITSNCls"), static_cast(origConstit.itsNCls())); - registry.fill(HIST("h_gnnfeat_trackTPCNCls"), static_cast(origConstit.tpcNClsFound())); - registry.fill(HIST("h_gnnfeat_trackTPCNCrossedRows"), static_cast(origConstit.tpcNClsCrossedRows())); - registry.fill(HIST("h_gnnfeat_trackITSChi2NCl"), origConstit.itsChi2NCl()); - registry.fill(HIST("h_gnnfeat_trackTPCChi2NCl"), origConstit.tpcChi2NCl()); + // registry.fill(HIST("h_gnnfeat_trackITSNCls"), static_cast(origConstit.itsNCls())); + // registry.fill(HIST("h_gnnfeat_trackTPCNCls"), static_cast(origConstit.tpcNClsFound())); + // registry.fill(HIST("h_gnnfeat_trackTPCNCrossedRows"), static_cast(origConstit.tpcNClsCrossedRows())); + // registry.fill(HIST("h_gnnfeat_trackITSChi2NCl"), origConstit.itsChi2NCl()); + // registry.fill(HIST("h_gnnfeat_trackTPCChi2NCl"), origConstit.tpcChi2NCl()); registry.fill(HIST("h_gnnfeat_tracksIPxy"), std::abs(constituent.dcaXY()) * sign / constituent.sigmadcaXY()); registry.fill(HIST("h_gnnfeat_tracksIPz"), std::abs(constituent.dcaZ()) * sign / constituent.sigmadcaZ()); @@ -504,7 +510,6 @@ struct BjetTaggingGnn { template bool isAcceptedJet(AnalysisJet const& jet) { - if (jetAreaFractionMin > -98.0) { if (jet.area() < jetAreaFractionMin * M_PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { return false; @@ -537,10 +542,10 @@ struct BjetTaggingGnn { return true; } - template - void fillDataJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& allTracks, AnyOriginalTracks const& origTracks) + template + void fillDataJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& allTracks) { - int nTracks = analyzeJetTrackInfo(analysisJet, allTracks, origTracks); + int nTracks = analyzeJetTrackInfo(analysisJet, allTracks); registry.fill(HIST("h_jetpT"), analysisJet.pt()); registry.fill(HIST("h_jetPhi"), analysisJet.phi()); @@ -559,8 +564,8 @@ struct BjetTaggingGnn { } } - template - int8_t fillMCDJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, AnyOriginalTracks const& /*origTracks*/, double weightEvt = 1.0) + template + int8_t fillMCDJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, double weightEvt = 1.0) { int8_t jetFlavor = analysisJet.origin(); @@ -575,24 +580,24 @@ struct BjetTaggingGnn { } int sign = jettaggingutilities::getGeoSign(analysisJet, constituent); - auto origConstit = constituent.template track_as(); - - registry.fill(HIST("h_gnnfeat_trackpT"), constituent.pt()); - registry.fill(HIST("h_gnnfeat_trackPhi"), origConstit.phi()); - registry.fill(HIST("h_gnnfeat_trackEta"), constituent.eta()); - registry.fill(HIST("h_gnnfeat_trackCharge"), static_cast(constituent.sign())); - registry.fill(HIST("h_gnnfeat_trackDCAxy"), std::abs(constituent.dcaXY()) * sign); - registry.fill(HIST("h_gnnfeat_trackSigmaDCAxy"), constituent.sigmadcaXY()); - registry.fill(HIST("h_gnnfeat_trackDCAz"), std::abs(constituent.dcaZ()) * sign); - registry.fill(HIST("h_gnnfeat_trackSigmaDCAz"), constituent.sigmadcaZ()); - registry.fill(HIST("h_gnnfeat_trackITSNCls"), static_cast(origConstit.itsNCls())); - registry.fill(HIST("h_gnnfeat_trackTPCNCls"), static_cast(origConstit.tpcNClsFound())); - registry.fill(HIST("h_gnnfeat_trackTPCNCrossedRows"), static_cast(origConstit.tpcNClsCrossedRows())); - registry.fill(HIST("h_gnnfeat_trackITSChi2NCl"), origConstit.itsChi2NCl()); - registry.fill(HIST("h_gnnfeat_trackTPCChi2NCl"), origConstit.tpcChi2NCl()); - - registry.fill(HIST("h_gnnfeat_tracksIPxy"), std::abs(constituent.dcaXY()) * sign / constituent.sigmadcaXY()); - registry.fill(HIST("h_gnnfeat_tracksIPz"), std::abs(constituent.dcaZ()) * sign / constituent.sigmadcaZ()); + // auto origConstit = constituent.template track_as(); + + registry.fill(HIST("h_gnnfeat_trackpT"), constituent.pt(), weightEvt); + registry.fill(HIST("h_gnnfeat_trackPhi"), constituent.phi(), weightEvt); + registry.fill(HIST("h_gnnfeat_trackEta"), constituent.eta(), weightEvt); + registry.fill(HIST("h_gnnfeat_trackCharge"), static_cast(constituent.sign()), weightEvt); + registry.fill(HIST("h_gnnfeat_trackDCAxy"), std::abs(constituent.dcaXY()) * sign, weightEvt); + registry.fill(HIST("h_gnnfeat_trackSigmaDCAxy"), constituent.sigmadcaXY(), weightEvt); + registry.fill(HIST("h_gnnfeat_trackDCAz"), std::abs(constituent.dcaZ()) * sign, weightEvt); + registry.fill(HIST("h_gnnfeat_trackSigmaDCAz"), constituent.sigmadcaZ(), weightEvt); + // registry.fill(HIST("h_gnnfeat_trackITSNCls"), static_cast(origConstit.itsNCls()), weightEvt); + // registry.fill(HIST("h_gnnfeat_trackTPCNCls"), static_cast(origConstit.tpcNClsFound()), weightEvt); + // registry.fill(HIST("h_gnnfeat_trackTPCNCrossedRows"), static_cast(origConstit.tpcNClsCrossedRows()), weightEvt); + // registry.fill(HIST("h_gnnfeat_trackITSChi2NCl"), origConstit.itsChi2NCl(), weightEvt); + // registry.fill(HIST("h_gnnfeat_trackTPCChi2NCl"), origConstit.tpcChi2NCl(), weightEvt); + + registry.fill(HIST("h_gnnfeat_tracksIPxy"), std::abs(constituent.dcaXY()) * sign / constituent.sigmadcaXY(), weightEvt); + registry.fill(HIST("h_gnnfeat_tracksIPz"), std::abs(constituent.dcaZ()) * sign / constituent.sigmadcaZ(), weightEvt); ++nTracks; } @@ -725,10 +730,10 @@ struct BjetTaggingGnn { registry.fill(HIST("hCollCounter"), static_cast(EvtSel::Sel8), weightEvt); // Coll+TVX+NoTFB+NoITSROFB if (zvtx) { registry.fill(HIST("hCollCounter"), static_cast(EvtSel::Sel8Zvtx), weightEvt); // Coll+TVX+NoTFB+NoITSROFB+Zvtx - if (collision.template collision_as().isInelGt0()) { - evtselCode |= EvtSelFlag::kINELgt0rec; - registry.fill(HIST("hCollCounter"), static_cast(EvtSel::INELgt0rec), weightEvt); // INELgt0+Zvtx(rec) - } + // if (collision.template collision_as().isInelGt0()) { + // evtselCode |= EvtSelFlag::kINELgt0rec; + // registry.fill(HIST("hCollCounter"), static_cast(EvtSel::INELgt0rec), weightEvt); // INELgt0+Zvtx(rec) + // } } } } @@ -752,7 +757,7 @@ struct BjetTaggingGnn { } PROCESS_SWITCH(BjetTaggingGnn, processDummy, "Dummy process function turned on by default", true); - void processDataJets(FilteredCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& allTracks, OriginalTracks const& origTracks) + void processDataJets(FilteredCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& allTracks) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; @@ -766,12 +771,12 @@ struct BjetTaggingGnn { continue; } - fillDataJetHistograms(analysisJet, allTracks, origTracks); + fillDataJetHistograms(analysisJet, allTracks); } } PROCESS_SWITCH(BjetTaggingGnn, processDataJets, "jet information in Data", false); - void processDataJetsTrig(FilteredCollisionsTriggered::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& allTracks, OriginalTracks const& origTracks, aod::JBCs const& /*bcInfo*/) + void processDataJetsTrig(FilteredCollisionsTriggered::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& allTracks, aod::JBCs const& /*bcInfo*/) { // Get BC info associated with the collision before applying any event selections auto bc = collision.bc_as(); @@ -794,12 +799,12 @@ struct BjetTaggingGnn { continue; } - fillDataJetHistograms(analysisJet, allTracks, origTracks); + fillDataJetHistograms(analysisJet, allTracks); } } PROCESS_SWITCH(BjetTaggingGnn, processDataJetsTrig, "jet information in software triggered Data", false); - void processDataJetsSel(AnalysisCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& /*allTracks*/, OrigCollisions const& /*origCollisions*/) + void processDataJetsSel(AnalysisCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& /*allTracks*/) { EvtSelFlag evtselCode = fillCollCounter(collision); @@ -816,7 +821,6 @@ struct BjetTaggingGnn { registry.fill(HIST("h_jetpT_selmc_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx)); registry.fill(HIST("h_jetpT_sel8"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8)); registry.fill(HIST("h_jetpT_sel8_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx)); - registry.fill(HIST("h_jetpT_inelgt0rec"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0rec)); } } PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSel, "jet information in Data (event selection)", false); @@ -842,11 +846,18 @@ struct BjetTaggingGnn { } PROCESS_SWITCH(BjetTaggingGnn, processDataTracks, "track information in Data", false); - void processMCDJets(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& allTracks, OriginalTracks const& origTracks, FilteredMCPJets const& /*MCPjets*/, aod::JetParticles const& /*mcParticles*/, FilteredCollisionsMCP const& /*mcCollisions*/) + void processMCDJets(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& allTracks, FilteredMCPJets const& /*MCPjets*/, aod::JetParticles const& /*mcParticles*/, FilteredCollisionsMCP const& /*mcCollisions*/) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } + // Reject outlier MC collisions + if (collision.isOutlier()) { + return; + } + if (collision.has_mcCollision() && collision.template mcCollision_as().isOutlier()) { + return; + } // Uses only collisionId % trainingDatasetRaioParam != 0 for evaluation dataset if (trainingDatasetRatioParam && collision.collisionId() % trainingDatasetRatioParam == 0) { return; @@ -864,21 +875,12 @@ struct BjetTaggingGnn { continue; } - float pTHat = 10. / (std::pow(weightEvt, 1.0 / pTHatExponent)); - if (useEventWeight && analysisJet.pt() > pTHatMaxMCD * pTHat) { - continue; - } - - int8_t jetFlavor = fillMCDJetHistograms(analysisJet, allTracks, origTracks, weightEvt); + int8_t jetFlavor = fillMCDJetHistograms(analysisJet, allTracks, weightEvt); if (!matchedMcColl) { continue; } for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { - if (mcpjet.pt() > pTHatMaxMCP * pTHat) { - continue; - } - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT"), analysisJet.pt(), mcpjet.pt(), weightEvt); registry.fill(HIST("h_jetpT_matched"), analysisJet.pt(), weightEvt); registry.fill(HIST("h_jetpT_particle_matched"), mcpjet.pt(), weightEvt); @@ -896,34 +898,26 @@ struct BjetTaggingGnn { } PROCESS_SWITCH(BjetTaggingGnn, processMCDJets, "jet information in MC", false); - void processMCDJetsSel(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& /*allTracks*/, FilteredMCPJets const& /*MCPjets*/, OrigCollisions const& /*origCollisions*/, FilteredCollisionsMCP const& /*mcCollisions*/, aod::JetParticles const& mcParticles) + void processMCDJetsSel(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& /*allTracks*/, FilteredMCPJets const& /*MCPjets*/, FilteredCollisionsMCP const& /*mcCollisions*/) { + // Reject outlier MC collisions + if (collision.isOutlier()) { + return; + } + if (collision.has_mcCollision() && collision.template mcCollision_as().isOutlier()) { + return; + } + float weightEvt = useEventWeight ? collision.weight() : 1.f; EvtSelFlag evtselCode = fillCollCounter(collision, weightEvt); - bool isTrueINELgt0 = collision.has_mcCollision() && isTrueINEL0(collision.template mcCollision_as(), mcParticles); for (const auto& analysisJet : MCDjets) { if (!isAcceptedJet(analysisJet)) { continue; } - float pTHat = 10. / (std::pow(weightEvt, 1.0 / pTHatExponent)); - if (useEventWeight && analysisJet.pt() > pTHatMaxMCD * pTHat) { - continue; - } - int8_t jetFlavor = analysisJet.origin(); - // Get matched particle-level jet pT - double mcpjetpT = -1.0; - for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { - if (useEventWeight && mcpjet.pt() > pTHatMaxMCP * pTHat) { - continue; - } - mcpjetpT = mcpjet.pt(); - } - bool isMatched = mcpjetpT > 0.0; - registry.fill(HIST("h_jetpT_coll"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_coll_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_tvx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); @@ -932,12 +926,6 @@ struct BjetTaggingGnn { registry.fill(HIST("h_jetpT_selmc_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_sel8"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_sel8_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_inelgt0rec"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0rec) ? weightEvt : 0.0); - if (isMatched) { - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_selmc"), analysisJet.pt(), mcpjetpT, hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_sel8"), analysisJet.pt(), mcpjetpT, hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_inelgt0"), analysisJet.pt(), mcpjetpT, isTrueINELgt0 && (hasAll(evtselCode, EvtSelFlag::INELgt0rec)) ? weightEvt : 0.0); - } if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("h_jetpT_b_coll"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_b_coll_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); @@ -947,142 +935,147 @@ struct BjetTaggingGnn { registry.fill(HIST("h_jetpT_b_selmc_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_b_sel8"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); registry.fill(HIST("h_jetpT_b_sel8_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_b_inelgt0rec"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0rec) ? weightEvt : 0.0); - if (isMatched) { - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_b_selmc"), analysisJet.pt(), mcpjetpT, hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_b_sel8"), analysisJet.pt(), mcpjetpT, hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_b_inelgt0"), analysisJet.pt(), mcpjetpT, isTrueINELgt0 && (hasAll(evtselCode, EvtSelFlag::INELgt0rec)) ? weightEvt : 0.0); - } } } } PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSel, "jet information in MC (event selection)", false); PresliceUnsorted collisionsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; - Preslice mcpjetsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; + Preslice mcpjetsPerMCPCollision = aod::jet::mcCollisionId; - void processMCPJets(aod::McCollisions::iterator const& mcCollision, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, OrigCollisions const& /*origCollisions*/, aod::JetParticles const& mcParticles) + void processMCPJets(AnalysisCollisionsMCP const& mcCollisions, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, aod::JetParticles const& mcParticles) { - float weightEvt = useEventWeight ? mcCollision.weight() : 1.f; - auto matchedCollisions = collisions.sliceBy(collisionsPerMCPCollision, mcCollision.globalIndex()); + // Subscribing AnalysisCollisionsMCP::iterator causes an issue related to unsorted JMcCollisionLbs index. + for (const auto& mcCollision : mcCollisions) { + if (mcCollision.isOutlier()) { + continue; + } + float weightEvt = useEventWeight ? mcCollision.weight() : 1.f; + auto matchedCollisions = collisions.sliceBy(collisionsPerMCPCollision, mcCollision.mcCollisionId()); - EvtSelFlag evtselCode = EvtSelFlag::INEL; - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INEL), weightEvt); // INEL - registry.fill(HIST("h_vertexZ_truth"), mcCollision.posZ(), weightEvt); + EvtSelFlag evtselCode = EvtSelFlag::INEL; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INEL), weightEvt); // INEL + registry.fill(HIST("h_vertexZ_truth"), mcCollision.posZ(), weightEvt); - bool zvtx = std::fabs(mcCollision.posZ()) < vertexZCut; - bool zvtxMatched = false; + bool zvtx = std::fabs(mcCollision.posZ()) < vertexZCut; + bool zvtxMatched = false; - bool isTrueINELgt0 = isTrueINEL0(mcCollision, mcParticles); - if (zvtx) { - evtselCode |= EvtSelFlag::kZvtx; - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELZvtx), weightEvt); // INEL+Zvtx - if (isTrueINELgt0) { - evtselCode |= EvtSelFlag::kINELgt0; - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELgt0), weightEvt); // INELgt0 + bool isTrueINELgt0 = isTrueINEL0(mcCollision, mcParticles); + if (zvtx) { + evtselCode |= EvtSelFlag::kZvtx; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELZvtx), weightEvt); // INEL+Zvtx + if (isTrueINELgt0) { + evtselCode |= EvtSelFlag::kINELgt0; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELgt0), weightEvt); // INELgt0 + } } - } - bool isMatchedToAnalysisSelection = false; + bool isMatchedToAnalysisSelection = false; - if (matchedCollisions.size() >= 1) { - zvtxMatched = std::fabs(matchedCollisions.begin().posZ()) < vertexZCut; - evtselCode |= EvtSelFlag::kColl; - registry.fill(HIST("h_vertexZ_truth_coll"), mcCollision.posZ(), weightEvt); - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Coll), weightEvt); // McColl(-> Coll) - if (zvtxMatched) { - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::CollZvtx), weightEvt); // McColl(-> Coll+Zvtx) - } - if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsTVX)) { - evtselCode |= EvtSelFlag::kTVX; - registry.fill(HIST("h_vertexZ_truth_tvx"), mcCollision.posZ(), weightEvt); - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::TVX), weightEvt); // McColl(-> Coll+TVX) + if (matchedCollisions.size() >= 1) { + zvtxMatched = std::fabs(matchedCollisions.begin().posZ()) < vertexZCut; + evtselCode |= EvtSelFlag::kColl; + registry.fill(HIST("h_vertexZ_truth_coll"), mcCollision.posZ(), weightEvt); + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Coll), weightEvt); // McColl(-> Coll) if (zvtxMatched) { - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::TVXZvtx), weightEvt); // McColl(-> Coll+TVX+Zvtx) + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::CollZvtx), weightEvt); // McColl(-> Coll+Zvtx) } - if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsSelMC)) { - evtselCode |= EvtSelFlag::kNoTFB; - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::SelMC), weightEvt); // McColl(-> Coll+TVX+NoTFB) + if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsTVX)) { + evtselCode |= EvtSelFlag::kTVX; + registry.fill(HIST("h_vertexZ_truth_tvx"), mcCollision.posZ(), weightEvt); + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::TVX), weightEvt); // McColl(-> Coll+TVX) if (zvtxMatched) { - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::SelMCZvtx), weightEvt); // McColl(-> Coll+TVX+NoTFB+Zvtx) + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::TVXZvtx), weightEvt); // McColl(-> Coll+TVX+Zvtx) } - if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsSel8)) { - evtselCode |= EvtSelFlag::kNoITSROFB; - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB) + if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsSelMC)) { + evtselCode |= EvtSelFlag::kNoTFB; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::SelMC), weightEvt); // McColl(-> Coll+TVX+NoTFB) if (zvtxMatched) { - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8Zvtx), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB+Zvtx) - if (matchedCollisions.begin().template collision_as().isInelGt0()) { - evtselCode |= EvtSelFlag::kINELgt0rec; - registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELgt0rec), weightEvt); // INELgt0+Zvtx(rec) + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::SelMCZvtx), weightEvt); // McColl(-> Coll+TVX+NoTFB+Zvtx) + } + if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsSel8)) { + evtselCode |= EvtSelFlag::kNoITSROFB; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB) + if (zvtxMatched) { + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8Zvtx), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB+Zvtx) + // if (matchedCollisions.begin().template collision_as().isInelGt0()) { + // evtselCode |= EvtSelFlag::kINELgt0rec; + // registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELgt0rec), weightEvt); // INELgt0+Zvtx(rec) + // } } } } } - } - if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBits) && zvtxMatched) { - isMatchedToAnalysisSelection = true; - registry.fill(HIST("h_event_counter_mcp"), 0.0, weightEvt); - } - } - - auto mcpjetspermcpcollision = mcpjets.sliceBy(mcpjetsPerMCPCollision, mcCollision.globalIndex()); - for (const auto& mcpjet : mcpjetspermcpcollision) { - bool jetIncluded = false; - for (const auto& jetR : jetRadiiValues) { - if (mcpjet.r() == static_cast(jetR * 100)) { - jetIncluded = true; - break; + if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBits) && zvtxMatched) { + isMatchedToAnalysisSelection = true; + registry.fill(HIST("h_event_counter_mcp"), 0.0, weightEvt); } } - if (!jetIncluded) { - continue; - } - - float pTHat = 10. / (std::pow(weightEvt, 1.0 / pTHatExponent)); - if (useEventWeight && mcpjet.pt() > pTHatMaxMCP * pTHat) { - continue; - } - - int8_t jetFlavor = mcpjet.origin(); - - registry.fill(HIST("h_jetpT_particle_inel"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INEL) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_inel_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_coll"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_coll_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_tvx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_tvx_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_selmc"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_selmc_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_sel8"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_sel8_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_inelgt0"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_inelgt0rec"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0rec) ? weightEvt : 0.0); - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h_jetpT_particle_b_inel"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INEL) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_inel_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_coll"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_coll_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_tvx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_tvx_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_selmc"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_selmc_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_sel8"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_sel8_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_inelgt0"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_particle_b_inelgt0rec"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0rec) ? weightEvt : 0.0); - } + auto mcpjetspermcpcollision = mcpjets.sliceBy(mcpjetsPerMCPCollision, mcCollision.mcCollisionId()); + for (const auto& mcpjet : mcpjetspermcpcollision) { + bool jetIncluded = false; + for (const auto& jetR : jetRadiiValues) { + if (mcpjet.r() == static_cast(jetR * 100)) { + jetIncluded = true; + break; + } + } - // Fill histograms for jets matched to the analysis event selection - if (isMatchedToAnalysisSelection) { - registry.fill(HIST("h_jetpT_particle"), mcpjet.pt(), weightEvt); + if (!jetIncluded) { + continue; + } + int8_t jetFlavor = mcpjet.origin(); + + registry.fill(HIST("h_jetpT_particle_inel"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INEL) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_inel_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_coll"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_coll_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_tvx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_tvx_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_selmc"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_selmc_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_sel8"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_sel8_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_inelgt0"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0) ? weightEvt : 0.0); if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h_jetpT_particle_b"), mcpjet.pt(), weightEvt); + registry.fill(HIST("h_jetpT_particle_b_inel"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INEL) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_inel_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_coll"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_coll_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_tvx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_tvx_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_selmc"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_selmc_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_sel8"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_sel8_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_b_inelgt0"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0) ? weightEvt : 0.0); } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("h_jetpT_particle_c"), mcpjet.pt(), weightEvt); - } else { - registry.fill(HIST("h_jetpT_particle_lf"), mcpjet.pt(), weightEvt); + registry.fill(HIST("h_jetpT_particle_c_inel"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INEL) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_inel_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_coll"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Coll) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_coll_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_tvx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_tvx_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_selmc"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_selmc_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::SelMCZvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_sel8"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_sel8_zvtx"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::Sel8Zvtx) ? weightEvt : 0.0); + registry.fill(HIST("h_jetpT_particle_c_inelgt0"), mcpjet.pt(), hasAll(evtselCode, EvtSelFlag::INELgt0) ? weightEvt : 0.0); + } + + // Fill histograms for jets matched to the analysis event selection + if (isMatchedToAnalysisSelection) { + registry.fill(HIST("h_jetpT_particle"), mcpjet.pt(), weightEvt); + + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("h_jetpT_particle_b"), mcpjet.pt(), weightEvt); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("h_jetpT_particle_c"), mcpjet.pt(), weightEvt); + } else { + registry.fill(HIST("h_jetpT_particle_lf"), mcpjet.pt(), weightEvt); + } } } } @@ -1093,15 +1086,22 @@ struct BjetTaggingGnn { void processMCDTracks(FilteredCollisionsMCD::iterator const& collision, AnalysisTracksMCD const& tracks, FilteredCollisionsMCP const& /*mcCollisions*/, aod::JetParticles const& allParticles) { - float weightEvt = useEventWeight ? collision.weight() : 1.f; if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } + // Reject outlier MC collisions + if (collision.isOutlier()) { + return; + } + if (collision.has_mcCollision() && collision.template mcCollision_as().isOutlier()) { + return; + } // Uses only collisionId % trainingDatasetRaioParam != 0 for evaluation dataset if (trainingDatasetRatioParam && collision.collisionId() % trainingDatasetRatioParam == 0) { return; } + float weightEvt = useEventWeight ? collision.weight() : 1.f; bool matchedMcColl = collision.has_mcCollision() && std::fabs(collision.template mcCollision_as().posZ()) < vertexZCut; for (const auto& track : tracks) { @@ -1158,8 +1158,13 @@ struct BjetTaggingGnn { registry.fill(HIST("h_dcaZ_coll_mismatched"), std::fabs(track.dcaZ()), weightEvt); } } else { - registry.fill(HIST("h_dcaXY_npp"), std::fabs(track.dcaXY()), weightEvt); - registry.fill(HIST("h_dcaZ_npp"), std::fabs(track.dcaZ()), weightEvt); + if (particle.mcCollisionId() == collision.mcCollisionId()) { + registry.fill(HIST("h_dcaXY_npp"), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("h_dcaZ_npp"), std::fabs(track.dcaZ()), weightEvt); + } else { + registry.fill(HIST("h_dcaXY_npp_mismatched"), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("h_dcaZ_npp_mismatched"), std::fabs(track.dcaZ()), weightEvt); + } } } } diff --git a/PWGJE/Tasks/bjetTreeCreator.cxx b/PWGJE/Tasks/bjetTreeCreator.cxx index a5c85560f26..4adb20dc780 100644 --- a/PWGJE/Tasks/bjetTreeCreator.cxx +++ b/PWGJE/Tasks/bjetTreeCreator.cxx @@ -136,6 +136,15 @@ DECLARE_SOA_TABLE(bjetTracksParams, "AOD", "BJETTRACKSPARAM", using bjetTracksParam = bjetTracksParams::iterator; DECLARE_SOA_TABLE(bjetTracksParamsExtra, "AOD", "BJETTRACKSEXTRA", + // o2::soa::Index<>, + trackInfo::TrackPhi, + trackInfo::TrackCharge, + trackInfo::TrackOrigin, + trackInfo::TrackVtxIndex); + +using bjetTracksParamExtra = bjetTracksParamsExtra::iterator; + +DECLARE_SOA_TABLE(bjetTracksParamsExtrb, "AOD", "BJETTRACKSEXTRB", // o2::soa::Index<>, trackInfo::TrackPhi, trackInfo::TrackCharge, @@ -147,7 +156,7 @@ DECLARE_SOA_TABLE(bjetTracksParamsExtra, "AOD", "BJETTRACKSEXTRA", trackInfo::TrackOrigin, trackInfo::TrackVtxIndex); -using bjetTracksParamExtra = bjetTracksParamsExtra::iterator; +using bjetTracksParamExtrb = bjetTracksParamsExtrb::iterator; namespace SVInfo { @@ -204,6 +213,7 @@ struct BJetTreeCreator { Produces bjetParamsExtraTable; Produces bjetTracksParamsTable; Produces bjetTracksExtraTable; + Produces bjetTracksExtrbTable; Produces bjetSVParamsTable; Produces bjetConstituentsTable; @@ -257,6 +267,7 @@ struct BJetTreeCreator { Configurable vtxRes{"vtxRes", 0.01, "Vertex position resolution (cluster size) for GNN vertex predictions (cm)"}; Configurable trainingDatasetRatioParam{"trainingDatasetRatioParam", 0, "Parameter for splitting training/evaluation datasets by collisionId"}; + Configurable produceExtrbTableGnn{"produceExtrbTableGnn", false, "Flag whether to produce the extra table with track quality variables for GNN training"}; std::vector eventSelectionBits; @@ -345,16 +356,18 @@ struct BJetTreeCreator { registry.add("h_trk_dcaz", "trk_dcaxyz;#it{DCA}_{z} (cm);Entries", {HistType::kTH1F, {{200, -0.1, 0.1}}}); registry.add("h_trk_sigmadcaxy", "trk_sigmadcaxy;#it{#sigma}_{#it{DCA}_{xy}} (cm);Entries", {HistType::kTH1F, {{200, 0., 0.1}}}); registry.add("h_trk_sigmadcaz", "trk_sigmadcaxyz;#it{#sigma}_{#it{DCA}_{z}} (cm);Entries", {HistType::kTH1F, {{200, 0., 0.1}}}); - registry.add("h_trk_itsncls", "trk_itsncls;ITS NCls;Entries", {HistType::kTH1F, {{10, 0., 10.}}}); - registry.add("h_trk_tpcncls", "trk_tpcncls;TPC NCls (Found);Entries", {HistType::kTH1F, {{200, 0., 200.}}}); - registry.add("h_trk_tpcncrs", "trk_tpcncrs;TPC NCrossedRows;Entries", {HistType::kTH1F, {{200, 0., 200.}}}); - registry.add("h_trk_itschi2ncl", "trk_itschi2ncl;ITS #it{#chi}^{2}/ndf;Entries", {HistType::kTH1F, {{200, 0., 20.}}}); - registry.add("h_trk_tpcchi2ncl", "trk_tpcchi2ncl;TPC #it{#chi}^{2}/ndf;Entries", {HistType::kTH1F, {{200, 0., 10.}}}); - registry.add("h2_trk_jtrackpt_vs_origtrackpt", "JTracks::pt vs Tracks::pt", {HistType::kTH2F, {{200, 0., 100.}, {200, 0., 100.}}}); + if (produceExtrbTableGnn) { + registry.add("h_trk_itsncls", "trk_itsncls;ITS NCls;Entries", {HistType::kTH1F, {{10, 0., 10.}}}); + registry.add("h_trk_tpcncls", "trk_tpcncls;TPC NCls (Found);Entries", {HistType::kTH1F, {{200, 0., 200.}}}); + registry.add("h_trk_tpcncrs", "trk_tpcncrs;TPC NCrossedRows;Entries", {HistType::kTH1F, {{200, 0., 200.}}}); + registry.add("h_trk_itschi2ncl", "trk_itschi2ncl;ITS #it{#chi}^{2}/ndf;Entries", {HistType::kTH1F, {{200, 0., 20.}}}); + registry.add("h_trk_tpcchi2ncl", "trk_tpcchi2ncl;TPC #it{#chi}^{2}/ndf;Entries", {HistType::kTH1F, {{200, 0., 10.}}}); + registry.add("h2_trk_jtrackpt_vs_origtrackpt", "JTracks::pt vs Tracks::pt", {HistType::kTH2F, {{200, 0., 100.}, {200, 0., 100.}}}); + } registry.add("h_trk_vtx_index", "trk_vtx_index;Vertex index;Entries", {HistType::kTH1F, {{20, 0., 20.}}}); registry.add("h_trk_origin", "trk_origin;Track origin;Entries", {HistType::kTH1F, {{5, 0., 5.}}}); auto hTrackOrigin = registry.get(HIST("h_trk_origin")); - hTrackOrigin->GetXaxis()->SetBinLabel(1, "NotPhysPrim"); + hTrackOrigin->GetXaxis()->SetBinLabel(1, "Background"); hTrackOrigin->GetXaxis()->SetBinLabel(2, "Charm"); hTrackOrigin->GetXaxis()->SetBinLabel(3, "Beauty"); hTrackOrigin->GetXaxis()->SetBinLabel(4, "Primary"); @@ -577,23 +590,29 @@ struct BJetTreeCreator { //+trk registry.fill(HIST("h_trk_pt"), constituent.pt(), eventweight); registry.fill(HIST("h_trk_eta"), constituent.eta(), eventweight); - registry.fill(HIST("h_trk_phi"), origConstit.phi(), eventweight); + registry.fill(HIST("h_trk_phi"), constituent.phi(), eventweight); registry.fill(HIST("h_trk_charge"), constituent.sign(), eventweight); registry.fill(HIST("h_trk_dcaxy"), std::abs(constituent.dcaXY()) * sign, eventweight); registry.fill(HIST("h_trk_dcaz"), std::abs(constituent.dcaZ()) * sign, eventweight); registry.fill(HIST("h_trk_sigmadcaxy"), constituent.sigmadcaXY(), eventweight); registry.fill(HIST("h_trk_sigmadcaz"), constituent.sigmadcaZ(), eventweight); - registry.fill(HIST("h_trk_itsncls"), origConstit.itsNCls(), eventweight); - registry.fill(HIST("h_trk_tpcncls"), origConstit.tpcNClsFound(), eventweight); - registry.fill(HIST("h_trk_tpcncrs"), origConstit.tpcNClsCrossedRows(), eventweight); - registry.fill(HIST("h_trk_itschi2ncl"), origConstit.itsChi2NCl(), eventweight); - registry.fill(HIST("h_trk_tpcchi2ncl"), origConstit.tpcChi2NCl(), eventweight); - registry.fill(HIST("h2_trk_jtrackpt_vs_origtrackpt"), constituent.pt(), origConstit.pt(), eventweight); // jtrack & new extra table are well joined (linear correlation) + if (produceExtrbTableGnn) { + registry.fill(HIST("h_trk_itsncls"), origConstit.itsNCls(), eventweight); + registry.fill(HIST("h_trk_tpcncls"), origConstit.tpcNClsFound(), eventweight); + registry.fill(HIST("h_trk_tpcncrs"), origConstit.tpcNClsCrossedRows(), eventweight); + registry.fill(HIST("h_trk_itschi2ncl"), origConstit.itsChi2NCl(), eventweight); + registry.fill(HIST("h_trk_tpcchi2ncl"), origConstit.tpcChi2NCl(), eventweight); + registry.fill(HIST("h2_trk_jtrackpt_vs_origtrackpt"), constituent.pt(), origConstit.pt(), eventweight); // jtrack & new extra table are well joined (linear correlation) + } registry.fill(HIST("h_trk_vtx_index"), trkVtxIndex); registry.fill(HIST("h_trk_origin"), trkOrigin); if (produceTree) { - bjetTracksExtraTable(/*bjetParamsTable.lastIndex() + 1, */ origConstit.phi(), constituent.sign(), origConstit.itsChi2NCl(), origConstit.tpcChi2NCl(), origConstit.itsNCls(), origConstit.tpcNClsFound(), origConstit.tpcNClsCrossedRows(), trkOrigin, trkVtxIndex); //+ + if (produceExtrbTableGnn) { + bjetTracksExtrbTable(/*bjetParamsTable.lastIndex() + 1, */ constituent.phi(), constituent.sign(), origConstit.itsChi2NCl(), origConstit.tpcChi2NCl(), origConstit.itsNCls(), origConstit.tpcNClsFound(), origConstit.tpcNClsCrossedRows(), trkOrigin, trkVtxIndex); //+ + } else { + bjetTracksExtraTable(/*bjetParamsTable.lastIndex() + 1, */ constituent.phi(), constituent.sign(), trkOrigin, trkVtxIndex); //+ + } } } @@ -762,10 +781,13 @@ struct BJetTreeCreator { } PROCESS_SWITCH(BJetTreeCreator, processMCJets, "jet information in MC", false); + using FilteredCollisionMCDOutlier = soa::Filtered>; using MCDJetTableNoSV = soa::Filtered>; + using AnalysisCollisionsMCPOutlier = soa::Join; using JetParticleswID = soa::Join; - void processMCJetsForGNN(FilteredCollisionMCD::iterator const& collision, aod::JMcCollisions const&, MCDJetTableNoSV const& MCDjets, MCPJetTable const& MCPjets, JetTracksMCDwID const& allTracks, JetParticleswID const& MCParticles, OriginalTracks const& origTracks, aod::McParticles const& origParticles) + // aod::McParticles::vx(), vy(), vz() are necessary for jettaggingutilities::vertexClustering() + void processMCJetsForGNN(FilteredCollisionMCDOutlier::iterator const& collision, MCDJetTableNoSV const& MCDjets, MCPJetTable const& MCPjets, JetTracksMCDwID const& allTracks, JetParticleswID const& MCParticles, AnalysisCollisionsMCPOutlier const&, OriginalTracks const& origTracks, aod::McParticles const& origParticles) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits) || (static_cast(std::rand()) / RAND_MAX < eventReductionFactor)) { return; @@ -776,6 +798,14 @@ struct BJetTreeCreator { return; } + // Reject outlier MC collisions + if (collision.isOutlier()) { + return; + } + if (collision.has_mcCollision() && collision.template mcCollision_as().isOutlier()) { + return; + } + float eventWeight = collision.weight(); registry.fill(HIST("h_vertexZ"), collision.posZ(), eventWeight); @@ -799,7 +829,7 @@ struct BJetTreeCreator { //+ TrackLabelMap trkLabels{{"trkVtxIndex", {}}, {"trkOrigin", {}}}; - int nVertices = jettaggingutilities::vertexClustering(collision.template mcCollision_as(), analysisJet, allTracks, MCParticles, origParticles, trkLabels, true, vtxRes, trackPtMin); + int nVertices = jettaggingutilities::vertexClustering(collision.template mcCollision_as(), analysisJet, allTracks, MCParticles, origParticles, trkLabels, true, vtxRes, trackPtMin); analyzeJetTrackInfoForGNN(collision, analysisJet, allTracks, origTracks, indicesTracks, jetFlavor, eventWeight, &trkLabels); registry.fill(HIST("h2_jetMass_jetpT"), analysisJet.pt(), analysisJet.mass(), eventWeight); @@ -828,6 +858,7 @@ struct BJetTreeCreator { if (produceTree) { bjetConstituentsTable(bjetParamsTable.lastIndex() + 1, indicesTracks, indicesSVs); bjetParamsTable(analysisJet.pt(), analysisJet.eta(), analysisJet.phi(), indicesTracks.size(), nVertices, analysisJet.mass(), jetFlavor, analysisJet.r()); + bjetParamsExtraTable(eventWeight); } } }