From b57555c0fab4392c3bf11218621d185c8a7d8223 Mon Sep 17 00:00:00 2001 From: Changhwan Choi Date: Wed, 19 Aug 2026 16:59:26 +0900 Subject: [PATCH 1/4] Updated GNN bjet: UE subtraction, SV mass template fit, refactoring using templates --- PWGJE/Tasks/bjetTaggingGnn.cxx | 1371 +++++++++++++++++++++++-------- PWGJE/Tasks/bjetTreeCreator.cxx | 19 +- 2 files changed, 1033 insertions(+), 357 deletions(-) diff --git a/PWGJE/Tasks/bjetTaggingGnn.cxx b/PWGJE/Tasks/bjetTaggingGnn.cxx index fd47112f547..92aef55088e 100644 --- a/PWGJE/Tasks/bjetTaggingGnn.cxx +++ b/PWGJE/Tasks/bjetTaggingGnn.cxx @@ -39,11 +39,16 @@ #include #include +#include +#include +#include #include #include #include +#include #include +#include #include #include @@ -61,9 +66,11 @@ enum class EvtSelFlag : uint8_t { kTVX = 1 << 2, kNoTFB = 1 << 3, kNoITSROFB = 1 << 4, - kZvtx = 1 << 5, - kINELgt0 = 1 << 6, - kINELgt0rec = 1 << 7, + kNoPileup = 1 << 5, + kIsGoodZvtx = 1 << 6, + kZvtx = 1 << 7, + // kINELgt0 = 1 << 8, + // kINELgt0rec = 1 << 9, INEL = kINEL, INELZvtx = kINEL | kZvtx, @@ -75,8 +82,12 @@ enum class EvtSelFlag : uint8_t { SelMCZvtx = kINEL | kColl | kTVX | kNoTFB | kZvtx, Sel8 = kINEL | kColl | kTVX | kNoTFB | kNoITSROFB, Sel8Zvtx = kINEL | kColl | kTVX | kNoTFB | kNoITSROFB | kZvtx, - INELgt0 = kINEL | kZvtx | kINELgt0, - INELgt0rec = kINEL | kZvtx | kColl | kTVX | kNoTFB | kNoITSROFB | kINELgt0rec + Sel8Full = kINEL | kColl | kTVX | kNoTFB | kNoITSROFB | kNoPileup, + Sel8FullZvtx = kINEL | kColl | kTVX | kNoTFB | kNoITSROFB | kNoPileup | kZvtx, + Sel8FullGood = kINEL | kColl | kTVX | kNoTFB | kNoITSROFB | kNoPileup | kIsGoodZvtx, + Sel8FullGoodZvtx = kINEL | kColl | kTVX | kNoTFB | kNoITSROFB | kNoPileup | kIsGoodZvtx | kZvtx, + // INELgt0 = kINEL | kZvtx | kINELgt0 + // INELgt0rec = kINEL | kZvtx | kColl | kTVX | kNoTFB | kNoITSROFB | kINELgt0rec }; constexpr EvtSelFlag operator|(EvtSelFlag a, EvtSelFlag b) { @@ -112,9 +123,54 @@ enum class EvtSel { SelMCZvtx, Sel8, Sel8Zvtx, - INELgt0, - INELgt0rec + Sel8Full, + Sel8FullZvtx, + Sel8FullGood, + Sel8FullGoodZvtx + // INELgt0, + // INELgt0rec }; +// Ordered list of event-selection stages shared by all "evtsel"-differential histograms below; +// bin index in the histogram == index in this table (all "evtsel" histograms span the full table, +// so a given event-selection stage always lands on the same bin at particle- and reco-level, matching +// hCollCounter/hMcCollCounter). Reco-level histograms only ever get bins filled for stages +// [kEvtSelStageRecoFirst, kEvtSelStageRecoLast]; bins before kEvtSelStageRecoFirst are unfillable +// placeholders (see setEvtSelStageAxisLabels()), mirroring hCollCounter's "_1"/"_2" bins. +// MC-truth (MCP) ones fill the full range. +struct EvtSelStage { + EvtSelFlag flag; + const char* label; +}; +constexpr std::array kEvtSelStages{{{EvtSelFlag::INEL, "INEL"}, + {EvtSelFlag::INELZvtx, "INEL+Zvtx"}, + {EvtSelFlag::Coll, "Coll"}, + {EvtSelFlag::CollZvtx, "Coll+Zvtx"}, + {EvtSelFlag::TVX, "TVX"}, + {EvtSelFlag::TVXZvtx, "TVX+Zvtx"}, + {EvtSelFlag::SelMC, "SelMC"}, + {EvtSelFlag::SelMCZvtx, "SelMC+Zvtx"}, + {EvtSelFlag::Sel8, "Sel8"}, + {EvtSelFlag::Sel8Zvtx, "Sel8+Zvtx"}, + {EvtSelFlag::Sel8Full, "Sel8Full"}, + {EvtSelFlag::Sel8FullZvtx, "Sel8Full+Zvtx"}, + {EvtSelFlag::Sel8FullGood, "Sel8FullGood"}, + {EvtSelFlag::Sel8FullGoodZvtx, "Sel8FullGood+Zvtx"}}}; +constexpr int kEvtSelStageRecoFirst = 2; // Coll +constexpr int kEvtSelStageRecoLast = 13; // Sel8FullGood+Zvtx + +// Labels a full-range [0, kEvtSelStages.size()-1] evtsel axis. Bins below `first` (unreachable at +// reco level) get numbered placeholder labels ("_1", "_2", ...) instead of a stage label, matching +// hCollCounter's placeholder bins for the same reason: keeps every "evtsel" histogram's bin index +// aligned to the same event-selection stage across particle- and reco-level histograms. +void setEvtSelStageAxisLabels(TAxis* axis, int first = 0, int last = kEvtSelStages.size() - 1) +{ + for (int i = 0; i < first; ++i) { + axis->SetBinLabel(i + 1, ("_" + std::to_string(i + 1)).c_str()); + } + for (int i = first; i <= last; ++i) { + axis->SetBinLabel(i + 1, kEvtSelStages[i].label); + } +} }; // namespace bjet_tagging_gnn_evtsel using namespace bjet_tagging_gnn_evtsel; @@ -133,7 +189,7 @@ struct BjetTaggingGnn { Configurable trackPtMax{"trackPtMax", 1000.0, "maximum track pT"}; Configurable trackEtaMin{"trackEtaMin", -0.9, "minimum track eta"}; Configurable trackEtaMax{"trackEtaMax", 0.9, "maximum track eta"}; - Configurable trackPtMinGnn{"trackPtMinGnn", 0.5, "minimum track pT for GNN inputs"}; + Configurable trackPtMinGnn{"trackPtMinGnn", 0.15, "minimum track pT for GNN inputs"}; Configurable maxIPxy{"maxIPxy", 10, "maximum track DCA in xy plane"}; Configurable maxIPz{"maxIPz", 10, "maximum track DCA in z direction"}; @@ -141,22 +197,23 @@ struct BjetTaggingGnn { Configurable trackNppCrit{"trackNppCrit", 0.95, "track not physical primary ratio"}; // jet level configurables - Configurable jetPtMin{"jetPtMin", 5.0, "minimum jet pT"}; + Configurable jetPtMin{"jetPtMin", 0.0, "minimum jet pT"}; Configurable jetPtMax{"jetPtMax", 1000.0, "maximum jet pT"}; - Configurable jetEtaMin{"jetEtaMin", -99.0, "minimum jet pseudorapidity"}; - Configurable jetEtaMax{"jetEtaMax", 99.0, "maximum jet pseudorapidity"}; + Configurable jetEtaMin{"jetEtaMin", -0.9, "minimum jet pseudorapidity"}; + Configurable jetEtaMax{"jetEtaMax", 0.9, "maximum jet pseudorapidity"}; Configurable leadingConstituentPtMin{"leadingConstituentPtMin", -99.0, "minimum pT selection on jet constituent"}; Configurable leadingConstituentPtMax{"leadingConstituentPtMax", 9999.0, "maximum pT selection on jet constituent"}; Configurable jetAreaFractionMin{"jetAreaFractionMin", -99.0, "used to make a cut on the jet areas"}; Configurable> jetRadii{"jetRadii", std::vector{0.4}, "jet resolution parameters"}; - Configurable dbMin{"dbMin", -10., "minimum GNN Db"}; - Configurable dbMax{"dbMax", 20., "maximum GNN Db"}; + Configurable dbMin{"dbMin", -15., "minimum GNN Db"}; + Configurable dbMax{"dbMax", 15., "maximum GNN Db"}; Configurable dbNbins{"dbNbins", 3000, "number of bins in axisDbFine"}; Configurable doDataDriven{"doDataDriven", false, "Flag whether to use fill THnSpase for data driven methods"}; Configurable doDataDrivenExtra{"doDataDrivenExtra", false, "Flag whether to add extra axes to THnSparses"}; + Configurable doDataDrivenSV{"doDataDrivenSV", false, "Flag whether to subscribe SV tables and to use fill THnSparse for data driven methods for SV"}; Configurable callSumw2{"callSumw2", false, "Flag whether to call THnSparse::Sumw2() for error calculation"}; Configurable trainingDatasetRatioParam{"trainingDatasetRatioParam", 0, "Parameter for splitting training/evaluation datasets by collisionId"}; @@ -177,6 +234,8 @@ struct BjetTaggingGnn { std::vector eventSelectionBitsTVX; std::vector eventSelectionBitsSelMC; std::vector eventSelectionBitsSel8; + std::vector eventSelectionBitsSel8Full; + std::vector eventSelectionBitsSel8FullGood; int trackSelectionBits; @@ -190,6 +249,8 @@ struct BjetTaggingGnn { eventSelectionBitsTVX = jetderiveddatautilities::initialiseEventSelectionBits("TVX"); eventSelectionBitsSel8 = jetderiveddatautilities::initialiseEventSelectionBits("sel8"); eventSelectionBitsSelMC = jetderiveddatautilities::initialiseEventSelectionBits("selMC"); + eventSelectionBitsSel8Full = jetderiveddatautilities::initialiseEventSelectionBits("sel8Full"); + eventSelectionBitsSel8FullGood = jetderiveddatautilities::initialiseEventSelectionBits("sel8Full+IsGoodZvtxFT0vsPV"); trackSelectionBits = jetderiveddatautilities::initialiseTrackSelection(static_cast(trackSelections)); @@ -200,7 +261,7 @@ struct BjetTaggingGnn { registry.add("h_event_counter", ";analysis collision counter", {HistType::kTH1F, {{1, 0.0, 1.0}}}, callSumw2); registry.add("h_event_counter_mcp", ";analysis collision matched MC collision counter", {HistType::kTH1F, {{1, 0.0, 1.0}}}, callSumw2); registry.add("h_vertexZ", "Vertex Z;#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); - registry.add("hCollCounter", ";collision counter", {HistType::kTH1F, {{12, 1.0, 13.0}}}, callSumw2); + registry.add("hCollCounter", ";collision counter", {HistType::kTH1F, {{14, 1.0, 15.0}}}, callSumw2); auto hCollCounter = registry.get(HIST("hCollCounter")); hCollCounter->GetXaxis()->SetBinLabel(1, "_1"); hCollCounter->GetXaxis()->SetBinLabel(2, "_2"); @@ -212,9 +273,13 @@ struct BjetTaggingGnn { hCollCounter->GetXaxis()->SetBinLabel(8, "Coll+TVX+NoTFB+Zvtx"); hCollCounter->GetXaxis()->SetBinLabel(9, "Coll+TVX+NoTFB+NoITSROFB"); // sel8 hCollCounter->GetXaxis()->SetBinLabel(10, "Coll+TVX+NoTFB+NoITSROFB+Zvtx"); - hCollCounter->GetXaxis()->SetBinLabel(11, "_11"); - hCollCounter->GetXaxis()->SetBinLabel(12, "INELgt0+Zvtx(rec)"); // sel8 - registry.add("hMcCollCounter", ";MC collision counter", {HistType::kTH1F, {{12, 1.0, 13.0}}}, callSumw2); + hCollCounter->GetXaxis()->SetBinLabel(11, "Coll+TVX+NoTFB+NoITSROFB+NoPileup"); // sel8Full + hCollCounter->GetXaxis()->SetBinLabel(12, "Coll+TVX+NoTFB+NoITSROFB+NoPileup+Zvtx"); + hCollCounter->GetXaxis()->SetBinLabel(13, "Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx"); // sel8FullGood + hCollCounter->GetXaxis()->SetBinLabel(14, "Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx+Zvtx"); + // hCollCounter->GetXaxis()->SetBinLabel(11, "_11"); + // hCollCounter->GetXaxis()->SetBinLabel(12, "INELgt0+Zvtx(rec)"); // sel8 + registry.add("hMcCollCounter", ";MC collision counter", {HistType::kTH1F, {{14, 1.0, 15.0}}}, callSumw2); auto hMcCollCounter = registry.get(HIST("hMcCollCounter")); hMcCollCounter->GetXaxis()->SetBinLabel(1, "McColl(INEL)"); hMcCollCounter->GetXaxis()->SetBinLabel(2, "McColl+Zvtx"); @@ -226,8 +291,31 @@ struct BjetTaggingGnn { hMcCollCounter->GetXaxis()->SetBinLabel(8, "McColl+Zvtx(-> Coll+TVX+NoTFB+Zvtx)"); hMcCollCounter->GetXaxis()->SetBinLabel(9, "McColl(-> Coll+TVX+NoTFB+NoITSROFB)"); // sel8 hMcCollCounter->GetXaxis()->SetBinLabel(10, "McColl+Zvtx(-> Coll+TVX+NoTFB+NoITSROFB+Zvtx)"); - hMcCollCounter->GetXaxis()->SetBinLabel(11, "INELgt0+Zvtx"); - hMcCollCounter->GetXaxis()->SetBinLabel(12, "INELgt0+Zvtx(-> INELgt0+Zvtx(rec))"); // sel8 + hMcCollCounter->GetXaxis()->SetBinLabel(11, "McColl(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup)"); // sel8Full + hMcCollCounter->GetXaxis()->SetBinLabel(12, "McColl+Zvtx(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup+Zvtx)"); + hMcCollCounter->GetXaxis()->SetBinLabel(13, "McColl(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx)"); // sel8FullGood + hMcCollCounter->GetXaxis()->SetBinLabel(14, "McColl+Zvtx(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx+Zvtx)"); + // hMcCollCounter->GetXaxis()->SetBinLabel(11, "INELgt0+Zvtx"); + // hMcCollCounter->GetXaxis()->SetBinLabel(12, "INELgt0+Zvtx(-> INELgt0+Zvtx(rec))"); // sel8 + // Same 14-bin layout as hCollCounter/hMcCollCounter (see fillMcCollCounterBC()), but the ladder is + // evaluated on the MC-truth BC's raw evsel bits instead of a matched reco collision - bins 1-2 + // (INEL/INEL+Zvtx) are placeholders here too since that's already covered by hMcCollCounter. + registry.add("hMcCollCounterBC", ";MC collision counter (BC-based)", {HistType::kTH1F, {{14, 1.0, 15.0}}}, callSumw2); + auto hMcCollCounterBC = registry.get(HIST("hMcCollCounterBC")); + hMcCollCounterBC->GetXaxis()->SetBinLabel(1, "_1"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(2, "_2"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(3, "McColl(-> BC)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(4, "McColl+Zvtx(-> BC)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(5, "McColl(-> BC+TVX)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(6, "McColl+Zvtx(-> BC+TVX)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(7, "McColl(-> BC+TVX+NoTFB)"); // selMC + hMcCollCounterBC->GetXaxis()->SetBinLabel(8, "McColl+Zvtx(-> BC+TVX+NoTFB)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(9, "McColl(-> BC+TVX+NoTFB+NoITSROFB)"); // sel8 + hMcCollCounterBC->GetXaxis()->SetBinLabel(10, "McColl+Zvtx(-> BC+TVX+NoTFB+NoITSROFB)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(11, "McColl(-> BC+TVX+NoTFB+NoITSROFB+NoPileup)"); // sel8Full + hMcCollCounterBC->GetXaxis()->SetBinLabel(12, "McColl+Zvtx(-> BC+TVX+NoTFB+NoITSROFB+NoPileup)"); + hMcCollCounterBC->GetXaxis()->SetBinLabel(13, "McColl(-> BC+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx)"); // sel8FullGood + hMcCollCounterBC->GetXaxis()->SetBinLabel(14, "McColl+Zvtx(-> BC+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx)"); registry.add("hBCCounter", "", {HistType::kTH1F, {{10, 0.0, 10.}}}, callSumw2); auto hBCCounter = registry.get(HIST("hBCCounter")); hBCCounter->GetXaxis()->SetBinLabel(1, "BC"); @@ -243,22 +331,78 @@ struct BjetTaggingGnn { const AxisSpec axisTrackpT{200, 0., 200., "#it{p}_{T} (GeV/#it{c})"}; const AxisSpec axisTrackpTFine{1000, 0., 10., "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec axisJetpT{250, 0., 250., "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisJetpT{250, 0., 250., "#it{p}_{T,~ch~jet} (GeV/#it{c})"}; + // Used in place of axisJetpT for every "_sub"-suffixed histogram (UE-subtracted jet pT can go negative). + const AxisSpec axisJetpTSub{300, -50., 250., "#it{p}_{T,~ch~jet}^{sub} (GeV/#it{c})"}; const AxisSpec axisJetEta{200, -0.8, 0.8, "#it{#eta}_{jet}"}; const AxisSpec axisDb{200, dbMin, dbMax, "#it{D}_{b}"}; 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}"}; + // Jet-flavour category axis shared by every "_flavor"-suffixed histogram below (see getJetFlavorCat()/setJetFlavorCatAxisLabels()). + const AxisSpec axisJetFlavorCat{3, -0.5, 2.5, "flavour"}; + const AxisSpec axisLfMatchStatus{2, -0.5, 1.5, "lf match status"}; + // Event-selection-stage axis shared by every "_evtsel"-suffixed histogram below (reco and particle-level + // alike span the full table so a given stage always lands on the same bin, see kEvtSelStages/fillEvtSelStages()). + const AxisSpec axisEvtSelStage{static_cast(kEvtSelStages.size()), -0.5, static_cast(kEvtSelStages.size()) - 0.5, "evt. sel. stage"}; + const AxisSpec axisSVMass{200, 0., 10., "#it{m}_{SV} (GeV/#it{c}^{2})"}; + const AxisSpec axisSVLxyS{200, 0., 100., "L_{xy}/#sigma"}; + const AxisSpec axisSVDispersion{100, 0., 0.5, "SV dispersion"}; + + // `suffix` ("" or "_sub") is appended to jet-pT-dependent histogram names below. Histograms whose value + // doesn't depend on jet pT (h_SVMass, h_jetPhi, ...) are booked once, unsuffixed, outside these lambdas. + if (doprocessDataJetsSV || doprocessMCDJetsSV || doprocessDataJetsSVSub || doprocessMCDJetsSVSub) { + registry.add("h_SVMass", "", {HistType::kTH1F, {axisSVMass}}, callSumw2); + } + auto addSVHistograms = [&](const AxisSpec& axisJetpT, const std::string& suffix) { + // Best SV per jet (highest decay length significance); see fillDataJetHistogramsSV()/fillMCDJetHistogramsSV(). + registry.add("h2_SVMass_jetpT" + suffix, "", HistType::kTH2F, {axisJetpT, axisSVMass}, callSumw2); + registry.add("h2_SVLxyS_jetpT" + suffix, "", HistType::kTH2F, {axisJetpT, axisSVLxyS}, callSumw2); + registry.add("h2_SVDispersion_jetpT" + suffix, "", HistType::kTH2F, {axisJetpT, axisSVDispersion}, callSumw2); + }; + if (doprocessMCDJetsSV || doprocessMCDJetsSVSub) { + registry.add("h2_SVMass_flavor", "", {HistType::kTH2F, {axisSVMass, axisJetFlavorCat}}, callSumw2); + } + auto addSVFlavorHistograms = [&](const AxisSpec& axisJetpT, const std::string& suffix) { + auto hMass = registry.add("h3_SVMass_jetpT_flavor" + suffix, "", HistType::kTH3F, {axisJetpT, axisSVMass, axisJetFlavorCat}, callSumw2); + auto hLxyS = registry.add("h3_SVLxyS_jetpT_flavor" + suffix, "", HistType::kTH3F, {axisJetpT, axisSVLxyS, axisJetFlavorCat}, callSumw2); + auto hDisp = registry.add("h3_SVDispersion_jetpT_flavor" + suffix, "", HistType::kTH3F, {axisJetpT, axisSVDispersion, axisJetFlavorCat}, callSumw2); + setJetFlavorCatAxisLabels(hMass->GetZaxis()); + setJetFlavorCatAxisLabels(hLxyS->GetZaxis()); + setJetFlavorCatAxisLabels(hDisp->GetZaxis()); + }; + auto addCoreJetHistograms = [&](const AxisSpec& axisJetpT, const std::string& suffix) { + registry.add("h_jetpT" + suffix, "", HistType::kTH1F, {axisJetpT}, callSumw2); + registry.add("h2_jetpT_Db" + suffix, "", HistType::kTH2F, {axisJetpT, axisDb}); + }; + + if (doprocessDataJetsSV || doprocessMCDJetsSV) { + addSVHistograms(axisJetpT, ""); + } + if (doprocessDataJetsSVSub || doprocessMCDJetsSVSub) { + addSVHistograms(axisJetpTSub, "_sub"); + } + if (doprocessMCDJetsSV) { + addSVFlavorHistograms(axisJetpT, ""); + } + if (doprocessMCDJetsSVSub) { + addSVFlavorHistograms(axisJetpTSub, "_sub"); + } - registry.add("h_jetpT", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); + // h_jetEta/h_jetPhi/h_jetMass/h_Db/h2_nTracks_Db don't depend on jet pT - one shared copy regardless of + // `withSub` (see addCoreJetHistograms()'s comment above). registry.add("h_jetEta", "", {HistType::kTH1F, {axisJetEta}}, callSumw2); registry.add("h_jetPhi", "", {HistType::kTH1F, {{200, 0., 2. * M_PI, "#it{phi}_{jet}"}}}); registry.add("h_jetMass", "", {HistType::kTH1F, {axisJetMass}}); registry.add("h_Db", "", {HistType::kTH1F, {axisDbFine}}); - registry.add("h2_jetpT_Db", "", {HistType::kTH2F, {axisJetpT, axisDb}}); registry.add("h2_nTracks_Db", "", {HistType::kTH2F, {axisNTracks, axisDb}}); + addCoreJetHistograms(axisJetpT, ""); + if (doprocessDataJetsSub || doprocessDataJetsSVSub || doprocessMCDJetsSub || doprocessMCDJetsSVSub || doprocessDataJetsTrigSub) { + addCoreJetHistograms(axisJetpTSub, "_sub"); + } + registry.add("h_gnnfeat_trackpT", "", {HistType::kTH1F, {{200, 0., 100., "#it{p}_{T} (GeV/#it{c})"}}}); registry.add("h_gnnfeat_trackPhi", "", {HistType::kTH1F, {{200, 0., 2. * M_PI, "#it{#phi}"}}}); registry.add("h_gnnfeat_trackEta", "", {HistType::kTH1F, {{200, -0.9, 0.9, "#it{#eta}"}}}); @@ -281,7 +425,7 @@ struct BjetTaggingGnn { registry.add("h_trackphi", "", {HistType::kTH1F, {{100, 0.0, 2.0 * M_PI, "#it{#phi}"}}}, callSumw2); registry.add("h_dcaXY", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); registry.add("h_dcaZ", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); } if (doprocessMCDTracks) { @@ -308,177 +452,212 @@ struct BjetTaggingGnn { 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); - registry.add("h3_dca_pt_coll_fake", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_fake", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_coll_matched", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_coll_matched_b", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_coll_matched_c", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_coll_matched_lf", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_coll_mismatched", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_npp", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); - registry.add("h3_dca_pt_npp_mismatched", "", {HistType::kTH3F, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_coll_fake", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_fake", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_coll_matched", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_coll_matched_b", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_coll_matched_c", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_coll_matched_lf", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_coll_mismatched", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_npp", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); + registry.add("hSparse_dca_pt_npp_mismatched", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); } if (doprocessDataJetsSel || doprocessMCDJetsSel) { - registry.add("h_jetpT_sel8", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - 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_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); - registry.add("h_jetpT_coll_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - } - - if (doprocessMCDJets) { - registry.add("h_jetpT_b", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_c", "c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_lf", "lf-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_Db_b", "b-jet", {HistType::kTH1F, {axisDbFine}}); - registry.add("h_Db_c", "c-jet", {HistType::kTH1F, {axisDbFine}}); - registry.add("h_Db_lf", "lf-jet", {HistType::kTH1F, {axisDbFine}}); - registry.add("h2_jetpT_Db_b", "b-jet", {HistType::kTH2F, {axisJetpT, axisDb}}); - registry.add("h2_jetpT_Db_c", "c-jet", {HistType::kTH2F, {axisJetpT, axisDb}}); - registry.add("h2_jetpT_Db_lf", "lf-jet", {HistType::kTH2F, {axisJetpT, axisDb}}); - registry.add("h2_nTracks_Db_b", "b-jet", {HistType::kTH2F, {axisNTracks, axisDb}}); - registry.add("h2_nTracks_Db_c", "c-jet", {HistType::kTH2F, {axisNTracks, axisDb}}); - registry.add("h2_nTracks_Db_lf", "lf-jet", {HistType::kTH2F, {axisNTracks, axisDb}}); - registry.add("h2_Response_DetjetpT_PartjetpT", "", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_b", "b-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_c", "c-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_lf", "lf-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}, callSumw2}); - registry.add("h2_jetpT_Db_lf_none", "lf-jet (none)", {HistType::kTH2F, {axisJetpT, axisDb}}, callSumw2); - registry.add("h2_jetpT_Db_lf_matched", "lf-jet (matched)", {HistType::kTH2F, {axisJetpT, axisDb}}, callSumw2); - registry.add("h2_jetpT_Db_npp", "NotPhysPrim", {HistType::kTH2F, {axisJetpT, axisDb}}); - registry.add("h2_jetpT_Db_npp_b", "NotPhysPrim b-jet", {HistType::kTH2F, {axisJetpT, axisDb}}); - registry.add("h2_jetpT_Db_npp_c", "NotPhysPrim c-jet", {HistType::kTH2F, {axisJetpT, axisDb}}); - registry.add("h2_jetpT_Db_npp_lf", "NotPhysPrim lf-jet", {HistType::kTH2F, {axisJetpT, axisDb}}); + auto hEvtsel = registry.add("h2_jetpT_evtsel", "", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + setEvtSelStageAxisLabels(hEvtsel->GetYaxis(), kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + } + if (doprocessDataJetsSelSub || doprocessMCDJetsSelSub) { + auto hEvtselSub = registry.add("h2_jetpT_evtsel_sub", "", HistType::kTH2F, {axisJetpTSub, axisEvtSelStage}, callSumw2); + setEvtSelStageAxisLabels(hEvtselSub->GetYaxis(), kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + } + + // b/c/lf categories collapsed onto axisJetFlavorCat instead of one histogram per category. Also includes + // the geometric matched-response and pTHat-sparse histograms, which get UE-subtracted variants too, not + // just the core Db/flavor spectra. Histograms whose value doesn't depend on jetpT are booked once, above. + if (doprocessMCDJets || doprocessMCDJetsSV || doprocessMCDJetsSub || doprocessMCDJetsSVSub) { + registry.add("h2_Db_flavor", "", {HistType::kTH2F, {axisDbFine, axisJetFlavorCat}}); + registry.add("h3_nTracks_Db_flavor", "", {HistType::kTH3F, {axisNTracks, axisDb, axisJetFlavorCat}}); + registry.add("h2_Db_lfmatch", "lf-jet", {HistType::kTH2F, {axisDbFine, axisLfMatchStatus}}); registry.add("h_Db_npp", "NotPhysPrim", {HistType::kTH1F, {axisDbFine}}); - registry.add("h_Db_npp_b", "NotPhysPrim b-jet", {HistType::kTH1F, {axisDbFine}}); - registry.add("h_Db_npp_c", "NotPhysPrim c-jet", {HistType::kTH1F, {axisDbFine}}); - registry.add("h_Db_npp_lf", "NotPhysPrim lf-jet", {HistType::kTH1F, {axisDbFine}}); - registry.add("h_jetpT_matched", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_matched", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_matched", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b_matched", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h2_Db_npp_flavor", "NotPhysPrim", {HistType::kTH2F, {axisDbFine, axisJetFlavorCat}}); + registry.add("h2_Db_npp_lfmatch", "NotPhysPrim lf-jet", {HistType::kTH2F, {axisDbFine, axisLfMatchStatus}}); + setLfMatchStatusAxisLabels(registry.get(HIST("h2_Db_lfmatch"))->GetYaxis()); + setLfMatchStatusAxisLabels(registry.get(HIST("h2_Db_npp_lfmatch"))->GetYaxis()); + setJetFlavorCatAxisLabels(registry.get(HIST("h2_Db_flavor"))->GetYaxis()); + setJetFlavorCatAxisLabels(registry.get(HIST("h3_nTracks_Db_flavor"))->GetZaxis()); + setJetFlavorCatAxisLabels(registry.get(HIST("h2_Db_npp_flavor"))->GetYaxis()); + } + auto addMCDFlavorAndMatchedHistograms = [&](const AxisSpec& axisJetpT, const std::string& suffix) { + auto hJetpTFlavor = registry.add("h2_jetpT_flavor" + suffix, "", HistType::kTH2F, {axisJetpT, axisJetFlavorCat}, callSumw2); + auto hJetpTDbFlavor = registry.add("h3_jetpT_Db_flavor" + suffix, "", HistType::kTH3F, {axisJetpT, axisDb, axisJetFlavorCat}); + auto hJetpTLfmatch = registry.add("h2_jetpT_lfmatch" + suffix, "lf-jet", HistType::kTH2F, {axisJetpT, axisLfMatchStatus}, callSumw2); + auto hJetpTDbLfmatch = registry.add("h3_jetpT_Db_lfmatch" + suffix, "lf-jet", HistType::kTH3F, {axisJetpT, axisDb, axisLfMatchStatus}, callSumw2); + setLfMatchStatusAxisLabels(hJetpTLfmatch->GetYaxis()); + setLfMatchStatusAxisLabels(hJetpTDbLfmatch->GetZaxis()); + auto hResponse = registry.add("h2_Response_DetjetpT_PartjetpT" + suffix, "", HistType::kTH2F, {axisJetpT, axisJetpT}, callSumw2); + auto hResponseFlavor = registry.add("h3_Response_DetjetpT_PartjetpT_flavor" + suffix, "", HistType::kTH3F, {axisJetpT, axisJetpT, axisJetFlavorCat}, callSumw2); + registry.add("h2_jetpT_Db_npp" + suffix, "NotPhysPrim", HistType::kTH2F, {axisJetpT, axisDb}); + auto hJetpTDbNppFlavor = registry.add("h3_jetpT_Db_npp_flavor" + suffix, "NotPhysPrim", HistType::kTH3F, {axisJetpT, axisDb, axisJetFlavorCat}); + setJetFlavorCatAxisLabels(hJetpTFlavor->GetYaxis()); + setJetFlavorCatAxisLabels(hJetpTDbFlavor->GetZaxis()); + setJetFlavorCatAxisLabels(hResponseFlavor->GetZaxis()); + setJetFlavorCatAxisLabels(hJetpTDbNppFlavor->GetZaxis()); + (void)hResponse; + + registry.add("h_jetpT_matched" + suffix, "", HistType::kTH1F, {axisJetpT}, callSumw2); + registry.add("h_jetpT_particle_matched" + suffix, "", HistType::kTH1F, {axisJetpT}, callSumw2); + auto hJetpTMatchedFlavor = registry.add("h2_jetpT_matched_flavor" + suffix, "", HistType::kTH2F, {axisJetpT, axisJetFlavorCat}, callSumw2); + auto hJetpTParticleMatchedFlavor = registry.add("h2_jetpT_particle_matched_flavor" + suffix, "", HistType::kTH2F, {axisJetpT, axisJetFlavorCat}, callSumw2); + setJetFlavorCatAxisLabels(hJetpTMatchedFlavor->GetYaxis()); + setJetFlavorCatAxisLabels(hJetpTParticleMatchedFlavor->GetYaxis()); + // MCP jets + registry.add("h_jetpT_particle" + suffix, "", HistType::kTH1F, {axisJetpT}, callSumw2); + auto hJetpTParticleFlavor = registry.add("h2_jetpT_particle_flavor" + suffix, "", HistType::kTH2F, {axisJetpT, axisJetFlavorCat}, callSumw2); + setJetFlavorCatAxisLabels(hJetpTParticleFlavor->GetYaxis()); // pTHat study - registry.add("h3_pthat_jetpT", "", {HistType::kTH3F, {{300, 0., 300., "#hat{#it{p}}_{T} (GeV/#it{c})"}, axisJetpT, axisJetpT}}, callSumw2); - registry.add("h3_pthat_jetpT_b", "b-jet", {HistType::kTH3F, {{300, 0., 300., "#hat{#it{p}}_{T} (GeV/#it{c})"}, axisJetpT, axisJetpT}}, callSumw2); + registry.add("hSparse_pthat_jetpT" + suffix, "", HistType::kTHnSparseF, {{300, 0., 300., "#hat{#it{p}}_{T} (GeV/#it{c})"}, axisJetpT, axisJetpT}, callSumw2); + registry.add("hSparse_pthat_jetpT_b" + suffix, "b-jet", HistType::kTHnSparseF, {{300, 0., 300., "#hat{#it{p}}_{T} (GeV/#it{c})"}, axisJetpT, axisJetpT}, callSumw2); + registry.add("hSparse_pthat_jetpT_c" + suffix, "c-jet", HistType::kTHnSparseF, {{300, 0., 300., "#hat{#it{p}}_{T} (GeV/#it{c})"}, axisJetpT, axisJetpT}, callSumw2); + }; + + if (doprocessMCDJets || doprocessMCDJetsSV) { + addMCDFlavorAndMatchedHistograms(axisJetpT, ""); + } + if (doprocessMCDJetsSub || doprocessMCDJetsSVSub) { + addMCDFlavorAndMatchedHistograms(axisJetpTSub, "_sub"); } if (doprocessMCDJetsSel) { - registry.add("h_jetpT_b_sel8", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_sel8_zvtx", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_sel8", "", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_b_sel8", "b-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_selmc", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_selmc_zvtx", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_selmc", "", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_b_selmc", "b-jet", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_tvx", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_b_tvx_zvtx", "b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h2_Response_DetjetpT_PartjetpT_tvx", "", {HistType::kTH2F, {axisJetpT, axisJetpT}}, callSumw2); - 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("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); + auto hEvtselB = registry.add("h2_jetpT_evtsel_b", "b-jet", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + auto hEvtselC = registry.add("h2_jetpT_evtsel_c", "c-jet", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + setEvtSelStageAxisLabels(hEvtselB->GetYaxis(), kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + setEvtSelStageAxisLabels(hEvtselC->GetYaxis(), kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + } + if (doprocessMCDJetsSelSub) { + auto hEvtselBSub = registry.add("h2_jetpT_evtsel_b_sub", "b-jet", HistType::kTH2F, {axisJetpTSub, axisEvtSelStage}, callSumw2); + auto hEvtselCSub = registry.add("h2_jetpT_evtsel_c_sub", "c-jet", HistType::kTH2F, {axisJetpTSub, axisEvtSelStage}, callSumw2); + setEvtSelStageAxisLabels(hEvtselBSub->GetYaxis(), kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + setEvtSelStageAxisLabels(hEvtselCSub->GetYaxis(), kEvtSelStageRecoFirst, kEvtSelStageRecoLast); } - if (doprocessMCPJets) { + // Used by both processMCPJets ("", axisJetpT) and processMCPJetsSub ("_sub", axisJetpTSub). + auto addMCPEvtSelHistograms = [&](const AxisSpec& axisJetpT, const std::string& suffix) { + // All INEL/Coll/TVX/SelMC/Sel8(+Zvtx) variants collapsed onto axisEvtSelStage. + auto hEvtsel = registry.add("h2_jetpT_particle_evtsel" + suffix, "", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + auto hEvtselB = registry.add("h2_jetpT_particle_evtsel_b" + suffix, "particle b-jet", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + auto hEvtselC = registry.add("h2_jetpT_particle_evtsel_c" + suffix, "particle c-jet", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + setEvtSelStageAxisLabels(hEvtsel->GetYaxis()); + setEvtSelStageAxisLabels(hEvtselB->GetYaxis()); + setEvtSelStageAxisLabels(hEvtselC->GetYaxis()); + + // Same as above, but the event-selection stage is evaluated on the MC-truth BC directly + // (see fillMcCollCounterBC()) instead of via a matched reco collision. + auto hEvtselbc = registry.add("h2_jetpT_particle_evtselbc" + suffix, "", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + auto hEvtselbcB = registry.add("h2_jetpT_particle_evtselbc_b" + suffix, "particle b-jet", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + auto hEvtselbcC = registry.add("h2_jetpT_particle_evtselbc_c" + suffix, "particle c-jet", HistType::kTH2F, {axisJetpT, axisEvtSelStage}, callSumw2); + setEvtSelStageAxisLabels(hEvtselbc->GetYaxis()); + setEvtSelStageAxisLabels(hEvtselbcB->GetYaxis()); + setEvtSelStageAxisLabels(hEvtselbcC->GetYaxis()); + }; + + // processMCPJetsCommon() fills these unconditionally regardless of `withSub` (see fillMcCollCounterBC()/ + // hMcCollCounter's convention) - must be booked whenever either the raw or Sub process is on, not just + // processMCPJets (processMCPJetsSub alone would otherwise fill an unbooked histogram and throw at runtime). + if (doprocessMCPJets || doprocessMCPJetsSub) { registry.add("h_vertexZ_truth", "Vertex Z truth;#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); registry.add("h_vertexZ_truth_coll", "Vertex Z truth (Coll);#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); registry.add("h_vertexZ_truth_tvx", "Vertex Z truth (TVX);#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); - - registry.add("h_jetpT_particle", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_b", "particle b-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_c", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); - 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_selmc", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_tvx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_coll", "", {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_selmc_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_tvx_zvtx", "", {HistType::kTH1F, {axisJetpT}}, callSumw2); - registry.add("h_jetpT_particle_coll_zvtx", "", {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_c_inelgt0", "particle c-jet", {HistType::kTH1F, {axisJetpT}}, callSumw2); + registry.add("h_vertexZ_truth_bctvx", "Vertex Z truth (BC+TVX);#it{Z} (cm)", {HistType::kTH1F, {{100, -20.0, 20.0}}}, callSumw2); + } + if (doprocessMCPJets) { + addMCPEvtSelHistograms(axisJetpT, ""); + } + if (doprocessMCPJetsSub) { + addMCPEvtSelHistograms(axisJetpTSub, "_sub"); } - if (doDataDriven) { - if (doDataDrivenExtra) { - registry.add("hSparse_Incljets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}}, callSumw2); - if (doprocessMCDJets) { - registry.add("hSparse_bjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}}, callSumw2); - registry.add("hSparse_cjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}}, callSumw2); - registry.add("hSparse_lfjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}}, callSumw2); - registry.add("hSparse_lfjets_none", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}}, callSumw2); - registry.add("hSparse_lfjets_matched", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}}, callSumw2); + auto addDataDrivenHistograms = [&](const AxisSpec& axisJetpT, const std::string& suffix, bool mcdJetsEnabled) { + // hSparse_Incljets' last axis is the best-SV mass instead of the jet mass when doDataDrivenSV is on + // (only meaningful together with processDataJetsSV/processMCDJetsSV, see fill*JetHistogramsSV()). + if (doDataDrivenSV) { + registry.add("hSparse_Incljets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisSVMass}, callSumw2); + } else if (doDataDrivenExtra) { + registry.add("hSparse_Incljets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}, callSumw2); + } else { + registry.add("hSparse_Incljets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}, callSumw2); + } + if (doDataDrivenSV) { + if (mcdJetsEnabled) { + registry.add("hSparse_bjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisSVMass}, callSumw2); + registry.add("hSparse_cjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisSVMass}, callSumw2); + registry.add("hSparse_lfjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisSVMass}, callSumw2); + registry.add("hSparse_lfjets_none" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisSVMass}, callSumw2); + registry.add("hSparse_lfjets_matched" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisSVMass}, callSumw2); + } + } else if (doDataDrivenExtra) { + if (mcdJetsEnabled) { + registry.add("hSparse_bjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}, callSumw2); + registry.add("hSparse_cjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}, callSumw2); + registry.add("hSparse_lfjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}, callSumw2); + registry.add("hSparse_lfjets_none" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}, callSumw2); + registry.add("hSparse_lfjets_matched" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks, axisJetMass}, callSumw2); } } else { - registry.add("hSparse_Incljets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); - if (doprocessMCDJets) { - registry.add("hSparse_bjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); - registry.add("hSparse_cjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); - registry.add("hSparse_lfjets", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); - registry.add("hSparse_lfjets_none", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); - registry.add("hSparse_lfjets_matched", "", {HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}}, callSumw2); + if (mcdJetsEnabled) { + registry.add("hSparse_bjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}, callSumw2); + registry.add("hSparse_cjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}, callSumw2); + registry.add("hSparse_lfjets" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}, callSumw2); + registry.add("hSparse_lfjets_none" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}, callSumw2); + registry.add("hSparse_lfjets_matched" + suffix, "", HistType::kTHnSparseF, {axisJetpT, axisDbFine, axisNTracks}, callSumw2); } } + }; + + if (doDataDriven) { + addDataDrivenHistograms(axisJetpT, "", doprocessMCDJets || doprocessMCDJetsSV); + if (doprocessDataJetsSub || doprocessDataJetsSVSub || doprocessMCDJetsSub || doprocessMCDJetsSVSub || doprocessDataJetsTrigSub) { + addDataDrivenHistograms(axisJetpTSub, "_sub", doprocessMCDJetsSub || doprocessMCDJetsSVSub); + } } } + // Filters Filter collisionFilter = nabs(aod::jcollision::posZ) < vertexZCut; Filter trackFilter = (aod::jtrack::pt >= trackPtMin && aod::jtrack::pt < trackPtMax && aod::jtrack::eta > trackEtaMin && aod::jtrack::eta < trackEtaMax); Filter jetFilter = (aod::jet::pt >= jetPtMin && aod::jet::pt < jetPtMax && aod::jet::eta < jetEtaMax - aod::jet::r / 100.f && aod::jet::eta > jetEtaMin + aod::jet::r / 100.f); - using AnalysisCollisions = soa::Join; + // Data + using AnalysisCollisions = soa::Join; using FilteredCollisions = soa::Filtered; - using AnalysisCollisionsTriggered = soa::Join; + using AnalysisCollisionsTriggered = soa::Join; using FilteredCollisionsTriggered = soa::Filtered; // using OrigCollisions = soa::Join; + using DataJets = soa::Join; using FilteredDataJets = soa::Filtered; + using DataJetsSV = soa::Join; + using FilteredDataJetsSV = soa::Filtered; + using AnalysisTracks = soa::Join; using FilteredTracks = soa::Filtered; // using OriginalTracks = soa::Join; + // MCD + using AnalysisCollisionsMCD = soa::Join; + using FilteredCollisionsMCD = soa::Filtered; + using MCDJets = soa::Join; using FilteredMCDJets = soa::Filtered; + using MCDJetsSV = soa::Join; + using FilteredMCDJetsSV = soa::Filtered; + using AnalysisTracksMCD = soa::Join; using FilteredTracksMCD = soa::Filtered; - using AnalysisCollisionsMCD = soa::Join; - using FilteredCollisionsMCD = soa::Filtered; + // MCP + using AnalysisCollisionsMCP = soa::Join; - using AnalysisCollisionsMCP = soa::Join; using MCPJets = soa::Join; using FilteredMCPJets = soa::Filtered; @@ -555,33 +734,169 @@ struct BjetTaggingGnn { return true; } - template - void fillDataJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& allTracks) + // Jet-flavour categories collapsed onto a single axis (instead of one histogram per category): b, c, lf. + static double getJetFlavorCat(int8_t jetFlavor) + { + if (jetFlavor == JetTaggingSpecies::beauty) { + return 0.; + } + if (jetFlavor == JetTaggingSpecies::charm) { + return 1.; + } + return 2.; // lf + } + + static void setJetFlavorCatAxisLabels(TAxis* axis) + { + axis->SetBinLabel(1, "b"); + axis->SetBinLabel(2, "c"); + axis->SetBinLabel(3, "lf"); + } + + // Fills one bin per satisfied event-selection stage in [first, last] (see kEvtSelStages). Bin index + // == stage index i (not offset by `first`), since the histogram's axis always spans the full table. + // Callers pass the "_sub"-suffixed HIST name directly for a UE-subtracted or MC-particle-context variant + // of the same histogram (single `registry`, see its declaration). + void fillEvtSelStages(auto const& histName, double xVal, EvtSelFlag evtselCode, double weight, int first, int last) + { + for (int i = first; i <= last; ++i) { + if (hasAll(evtselCode, kEvtSelStages[i].flag)) { + registry.fill(histName, xVal, static_cast(i), weight); + } + } + } + + // lf-jets further split by MC-particle-jet match status (see h3_jetpT_Db_lfmatch). + static double getLfMatchStatus(int8_t jetFlavor) + { + return (jetFlavor == JetTaggingSpecies::none) ? 1. : 0.; + } + + static void setLfMatchStatusAxisLabels(TAxis* axis) + { + axis->SetBinLabel(1, "matched"); + axis->SetBinLabel(2, "none"); + } + + // Selects the SV with the highest decay length significance out of all SVs matched to the jet, + // shared by fillDataJetHistogramsSV()/fillMCDJetHistogramsSV(). Fills the inclusive SV histograms + // (plus the flavour-split variants when flavorCat >= 0, i.e. for MCD jets) and returns the + // selected SV's mass (-1 if the jet has no SV). `jetpT` is passed in already raw-or-UE-subtracted by + // the caller (per `withSub`) rather than recomputed here, since this function doesn't otherwise need rho. + template + float fillSVHistograms(AnalysisJet const& analysisJet, SecondaryVertices const& /*allSVs*/, float jetpT, double weightEvt = 1.0, double flavorCat = -1.0) + { + auto svs = analysisJet.template secondaryVertices_as(); + if (svs.size() == 0) { + // Not jet-pT-dependent (no-SV sentinel) - one shared histogram regardless of `withSub`. + registry.fill(HIST("h_SVMass"), -1.f, weightEvt); + return -1.f; + } + const auto& sv = *std::max_element(svs.begin(), svs.end(), [](auto const& svA, auto const& svB) { + return (svA.decayLengthXY() / svA.errorDecayLengthXY()) < (svB.decayLengthXY() / svB.errorDecayLengthXY()); + }); + float massSV = sv.m(); + float lxyS = sv.decayLengthXY() / sv.errorDecayLengthXY(); + // h_SVMass/h2_SVMass_flavor don't use jetpT as a fill value - shared regardless of `withSub`. The + // jetpT-axis histograms below get the "_sub" name instead of a second registry (see `registry`'s + // declaration). + registry.fill(HIST("h_SVMass"), massSV, weightEvt); + if (flavorCat >= 0.) { + registry.fill(HIST("h2_SVMass_flavor"), massSV, flavorCat, weightEvt); + } + if constexpr (withSub) { + registry.fill(HIST("h2_SVMass_jetpT_sub"), jetpT, massSV, weightEvt); + registry.fill(HIST("h2_SVLxyS_jetpT_sub"), jetpT, lxyS, weightEvt); + registry.fill(HIST("h2_SVDispersion_jetpT_sub"), jetpT, sv.dispersion(), weightEvt); + if (flavorCat >= 0.) { + registry.fill(HIST("h3_SVMass_jetpT_flavor_sub"), jetpT, massSV, flavorCat, weightEvt); + registry.fill(HIST("h3_SVLxyS_jetpT_flavor_sub"), jetpT, lxyS, flavorCat, weightEvt); + registry.fill(HIST("h3_SVDispersion_jetpT_flavor_sub"), jetpT, sv.dispersion(), flavorCat, weightEvt); + } + } else { + registry.fill(HIST("h2_SVMass_jetpT"), jetpT, massSV, weightEvt); + registry.fill(HIST("h2_SVLxyS_jetpT"), jetpT, lxyS, weightEvt); + registry.fill(HIST("h2_SVDispersion_jetpT"), jetpT, sv.dispersion(), weightEvt); + if (flavorCat >= 0.) { + registry.fill(HIST("h3_SVMass_jetpT_flavor"), jetpT, massSV, flavorCat, weightEvt); + registry.fill(HIST("h3_SVLxyS_jetpT_flavor"), jetpT, lxyS, flavorCat, weightEvt); + registry.fill(HIST("h3_SVDispersion_jetpT_flavor"), jetpT, sv.dispersion(), flavorCat, weightEvt); + } + } + return massSV; + } + + // `withSub` selects UE(rho*area)-subtracted jet pT + the "_sub"-suffixed histogram set in the single + // `registry` instead of raw pT + the unsuffixed names - used by processDataJets(SV)/processDataJets(SV)Sub. + template + int fillDataJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& allTracks, float rho = 0.0) { int nTracks = analyzeJetTrackInfo(analysisJet, allTracks); - registry.fill(HIST("h_jetpT"), analysisJet.pt()); + float jetpT = withSub ? (analysisJet.pt() - rho * analysisJet.area()) : analysisJet.pt(); + // h_jetPhi/h_jetEta/h_jetMass/h_Db/h2_nTracks_Db don't use jetpT as a fill value - one shared histogram + // regardless of `withSub` (see addCoreJetHistograms()'s init() comment). Only h_jetpT/h2_jetpT_Db/ + // hSparse_Incljets actually differ between the raw and Sub call, so only those get the "_sub" name. registry.fill(HIST("h_jetPhi"), analysisJet.phi()); registry.fill(HIST("h_jetEta"), analysisJet.eta()); registry.fill(HIST("h_jetMass"), analysisJet.mass()); registry.fill(HIST("h_Db"), analysisJet.scoreML()); - registry.fill(HIST("h2_jetpT_Db"), analysisJet.pt(), analysisJet.scoreML()); registry.fill(HIST("h2_nTracks_Db"), nTracks, analysisJet.scoreML()); - if (doDataDriven) { - if (doDataDrivenExtra) { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass()); + if constexpr (withSub) { + registry.fill(HIST("h_jetpT_sub"), jetpT); + registry.fill(HIST("h2_jetpT_Db_sub"), jetpT, analysisJet.scoreML()); + + if (doDataDriven && !doDataDrivenSV) { + if (doDataDrivenExtra) { + registry.fill(HIST("hSparse_Incljets_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass()); + } else { + registry.fill(HIST("hSparse_Incljets_sub"), jetpT, analysisJet.scoreML(), nTracks); + } + } + } else { + registry.fill(HIST("h_jetpT"), jetpT); + registry.fill(HIST("h2_jetpT_Db"), jetpT, analysisJet.scoreML()); + + if (doDataDriven && !doDataDrivenSV) { + if (doDataDrivenExtra) { + registry.fill(HIST("hSparse_Incljets"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass()); + } else { + registry.fill(HIST("hSparse_Incljets"), jetpT, analysisJet.scoreML(), nTracks); + } + } + } + return nTracks; + } + + template + void fillDataJetHistogramsSV(AnalysisJet const& analysisJet, AnyTracks const& allTracks, AnySVs const& allSVs, float rho = 0.0) + { + int nTracks = fillDataJetHistograms(analysisJet, allTracks, rho); + float jetpT = withSub ? (analysisJet.pt() - rho * analysisJet.area()) : analysisJet.pt(); + float massSV = fillSVHistograms(analysisJet, allSVs, jetpT); + + if (doDataDriven && doDataDrivenSV) { + if constexpr (withSub) { + registry.fill(HIST("hSparse_Incljets_sub"), jetpT, analysisJet.scoreML(), nTracks, massSV); } else { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), nTracks); + registry.fill(HIST("hSparse_Incljets"), jetpT, analysisJet.scoreML(), nTracks, massSV); } } } + // jetFlavor: JetTaggingSpecies of the jet; nTracks: number of GNN-input constituents (reused by fillMCDJetHistogramsSV()). + struct JetHistFillResult { + int8_t jetFlavor; + int nTracks; + }; + + // Constituent loop factored out of fillMCDJetHistograms(): the h_gnnfeat_* track-QA histograms it fills + // don't depend on the jet's own (possibly UE-subtracted) pT, so they always use the unsuffixed names + // regardless of `withSub` - no "_sub" duplicates needed. template - int8_t fillMCDJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, double weightEvt = 1.0) + std::pair analyzeMCDJetTrackInfo(AnalysisJet const& analysisJet, AnyTracks const& /*allTracks*/, double weightEvt = 1.0) { - int8_t jetFlavor = analysisJet.origin(); - int nTracks = 0; int nNppTracks = 0; for (const auto& constituent : analysisJet.template tracks_as()) { @@ -614,85 +929,152 @@ struct BjetTaggingGnn { ++nTracks; } + return {nTracks, nNppTracks}; + } + + // `withSub` selects UE(rho*area)-subtracted jet pT + the "_sub"-suffixed histogram set, same convention as + // fillDataJetHistograms(). Histograms whose value doesn't depend on jetpT are shared, unsuffixed. + template + JetHistFillResult fillMCDJetHistograms(AnalysisJet const& analysisJet, AnyTracks const& allTracks, float rho = 0.f, double weightEvt = 1.0) + { + int8_t jetFlavor = analysisJet.origin(); + auto [nTracks, nNppTracks] = analyzeMCDJetTrackInfo(analysisJet, allTracks, weightEvt); + + float jetpT = withSub ? (analysisJet.pt() - rho * analysisJet.area()) : analysisJet.pt(); + double flavorCat = getJetFlavorCat(jetFlavor); - registry.fill(HIST("h_jetpT"), analysisJet.pt(), weightEvt); registry.fill(HIST("h_jetPhi"), analysisJet.phi(), weightEvt); registry.fill(HIST("h_jetEta"), analysisJet.eta(), weightEvt); registry.fill(HIST("h_jetMass"), analysisJet.mass(), weightEvt); registry.fill(HIST("h_Db"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); registry.fill(HIST("h2_nTracks_Db"), nTracks, analysisJet.scoreML(), weightEvt); - - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h_jetpT_b"), analysisJet.pt(), weightEvt); - registry.fill(HIST("h_Db_b"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_b"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_nTracks_Db_b"), nTracks, analysisJet.scoreML(), weightEvt); - } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("h_jetpT_c"), analysisJet.pt(), weightEvt); - registry.fill(HIST("h_Db_c"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_c"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_nTracks_Db_c"), nTracks, analysisJet.scoreML(), weightEvt); - } else { - registry.fill(HIST("h_jetpT_lf"), analysisJet.pt(), weightEvt); - registry.fill(HIST("h_Db_lf"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_lf"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_nTracks_Db_lf"), nTracks, analysisJet.scoreML(), weightEvt); - if (jetFlavor == JetTaggingSpecies::none) { - registry.fill(HIST("h2_jetpT_Db_lf_none"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - } else { - registry.fill(HIST("h2_jetpT_Db_lf_matched"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - } + registry.fill(HIST("h2_Db_flavor"), analysisJet.scoreML(), flavorCat, weightEvt); + registry.fill(HIST("h3_nTracks_Db_flavor"), nTracks, analysisJet.scoreML(), flavorCat, weightEvt); + bool isLf = jetFlavor != JetTaggingSpecies::beauty && jetFlavor != JetTaggingSpecies::charm; + if (isLf) { + registry.fill(HIST("h2_Db_lfmatch"), analysisJet.scoreML(), getLfMatchStatus(jetFlavor), weightEvt); } - - // Inspection for jets with predominant non-physical primary tracks - if (static_cast(nNppTracks) / nTracks > trackNppCrit) { + bool isNpp = static_cast(nNppTracks) / nTracks > trackNppCrit; + if (isNpp) { registry.fill(HIST("h_Db_npp"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_npp"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h_Db_npp_b"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_npp_b"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("h_Db_npp_c"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_npp_c"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); - } else { - registry.fill(HIST("h_Db_npp_lf"), analysisJet.scoreML(), weightEvt); - registry.fill(HIST("h2_jetpT_Db_npp_lf"), analysisJet.pt(), analysisJet.scoreML(), weightEvt); + registry.fill(HIST("h2_Db_npp_flavor"), analysisJet.scoreML(), flavorCat, weightEvt); + if (isLf) { + registry.fill(HIST("h2_Db_npp_lfmatch"), analysisJet.scoreML(), getLfMatchStatus(jetFlavor), weightEvt); } } - if (doDataDriven) { - if (doDataDrivenExtra) { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("hSparse_bjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); - } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("hSparse_cjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + if constexpr (withSub) { + registry.fill(HIST("h_jetpT_sub"), jetpT, weightEvt); + registry.fill(HIST("h2_jetpT_Db_sub"), jetpT, analysisJet.scoreML(), weightEvt); + registry.fill(HIST("h2_jetpT_flavor_sub"), jetpT, flavorCat, weightEvt); + registry.fill(HIST("h3_jetpT_Db_flavor_sub"), jetpT, analysisJet.scoreML(), flavorCat, weightEvt); + if (isLf) { + registry.fill(HIST("h2_jetpT_lfmatch_sub"), jetpT, getLfMatchStatus(jetFlavor), weightEvt); + registry.fill(HIST("h3_jetpT_Db_lfmatch_sub"), jetpT, analysisJet.scoreML(), getLfMatchStatus(jetFlavor), weightEvt); + } + if (isNpp) { + registry.fill(HIST("h2_jetpT_Db_npp_sub"), jetpT, analysisJet.scoreML(), weightEvt); + registry.fill(HIST("h3_jetpT_Db_npp_flavor_sub"), jetpT, analysisJet.scoreML(), flavorCat, weightEvt); + } + + if (doDataDriven && !doDataDrivenSV) { + if (doDataDrivenExtra) { + registry.fill(HIST("hSparse_Incljets_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_bjets_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_cjets_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } else { + registry.fill(HIST("hSparse_lfjets_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + if (jetFlavor == JetTaggingSpecies::none) { + registry.fill(HIST("hSparse_lfjets_none_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } else { + registry.fill(HIST("hSparse_lfjets_matched_sub"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } + } } else { - registry.fill(HIST("hSparse_lfjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); - if (jetFlavor == JetTaggingSpecies::none) { - registry.fill(HIST("hSparse_lfjets_none"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + registry.fill(HIST("hSparse_Incljets_sub"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_bjets_sub"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_cjets_sub"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); } else { - registry.fill(HIST("hSparse_lfjets_matched"), analysisJet.pt(), analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + registry.fill(HIST("hSparse_lfjets_sub"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + if (jetFlavor == JetTaggingSpecies::none) { + registry.fill(HIST("hSparse_lfjets_none_sub"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + } else { + registry.fill(HIST("hSparse_lfjets_matched_sub"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + } } } - } else { - registry.fill(HIST("hSparse_Incljets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weightEvt); - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("hSparse_bjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weightEvt); - } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("hSparse_cjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weightEvt); + } + } else { + registry.fill(HIST("h_jetpT"), jetpT, weightEvt); + registry.fill(HIST("h2_jetpT_Db"), jetpT, analysisJet.scoreML(), weightEvt); + registry.fill(HIST("h2_jetpT_flavor"), jetpT, flavorCat, weightEvt); + registry.fill(HIST("h3_jetpT_Db_flavor"), jetpT, analysisJet.scoreML(), flavorCat, weightEvt); + if (isLf) { + registry.fill(HIST("h2_jetpT_lfmatch"), jetpT, getLfMatchStatus(jetFlavor), weightEvt); + registry.fill(HIST("h3_jetpT_Db_lfmatch"), jetpT, analysisJet.scoreML(), getLfMatchStatus(jetFlavor), weightEvt); + } + if (isNpp) { + registry.fill(HIST("h2_jetpT_Db_npp"), jetpT, analysisJet.scoreML(), weightEvt); + registry.fill(HIST("h3_jetpT_Db_npp_flavor"), jetpT, analysisJet.scoreML(), flavorCat, weightEvt); + } + + if (doDataDriven && !doDataDrivenSV) { + if (doDataDrivenExtra) { + registry.fill(HIST("hSparse_Incljets"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_bjets"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_cjets"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } else { + registry.fill(HIST("hSparse_lfjets"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + if (jetFlavor == JetTaggingSpecies::none) { + registry.fill(HIST("hSparse_lfjets_none"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } else { + registry.fill(HIST("hSparse_lfjets_matched"), jetpT, analysisJet.scoreML(), nTracks, analysisJet.mass(), weightEvt); + } + } } else { - registry.fill(HIST("hSparse_lfjets"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weightEvt); - if (jetFlavor == JetTaggingSpecies::none) { - registry.fill(HIST("hSparse_lfjets_none"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weightEvt); + registry.fill(HIST("hSparse_Incljets"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_bjets"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_cjets"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); } else { - registry.fill(HIST("hSparse_lfjets_matched"), analysisJet.pt(), analysisJet.scoreML(), nTracks, weightEvt); + registry.fill(HIST("hSparse_lfjets"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + if (jetFlavor == JetTaggingSpecies::none) { + registry.fill(HIST("hSparse_lfjets_none"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + } else { + registry.fill(HIST("hSparse_lfjets_matched"), jetpT, analysisJet.scoreML(), nTracks, weightEvt); + } } } } } + return {jetFlavor, nTracks}; + } + + // Same as fillMCDJetHistograms(), plus SV histograms (see fillSVHistograms()) and the SV-mass variant + // of hSparse_Incljets (only when doDataDrivenSV is on); used by processMCDJetsSV. + template + int8_t fillMCDJetHistogramsSV(AnalysisJet const& analysisJet, AnyTracks const& allTracks, AnySVs const& allSVs, float rho = 0.f, double weightEvt = 1.0) + { + auto [jetFlavor, nTracks] = fillMCDJetHistograms(analysisJet, allTracks, rho, weightEvt); + float jetpT = withSub ? (analysisJet.pt() - rho * analysisJet.area()) : analysisJet.pt(); + float massSV = fillSVHistograms(analysisJet, allSVs, jetpT, weightEvt, getJetFlavorCat(jetFlavor)); + + if (doDataDriven && doDataDrivenSV) { + if constexpr (withSub) { + registry.fill(HIST("hSparse_Incljets_sub"), jetpT, analysisJet.scoreML(), nTracks, massSV, weightEvt); + } else { + registry.fill(HIST("hSparse_Incljets"), jetpT, analysisJet.scoreML(), nTracks, massSV, weightEvt); + } + } + return jetFlavor; } @@ -715,6 +1097,28 @@ struct BjetTaggingGnn { return false; } + // Shared entry checks for MCD collisions, used by processMCDJets(SV) and processMCDTracks: + // event selection, MC-outlier rejection (own and matched mcCollision's), and eval/train split. + template + bool isRejectedMCDCollision(AnalysisCollision const& collision) + { + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return true; + } + // Reject outlier MC collisions + if (collision.isOutlier()) { + return true; + } + if (collision.has_mcCollision() && collision.template mcCollision_as().isOutlier()) { + return true; + } + // Uses only collisionId % trainingDatasetRaioParam != 0 for evaluation dataset + if (trainingDatasetRatioParam && collision.collisionId() % trainingDatasetRatioParam == 0) { + return true; + } + return false; + } + template EvtSelFlag fillCollCounter(AnalysisCollision const& collision, float weightEvt = 1.f) { @@ -749,6 +1153,77 @@ struct BjetTaggingGnn { // registry.fill(HIST("hCollCounter"), static_cast(EvtSel::INELgt0rec), weightEvt); // INELgt0+Zvtx(rec) // } } + if (jetderiveddatautilities::selectCollision(collision, eventSelectionBitsSel8Full)) { + evtselCode |= EvtSelFlag::kNoPileup; + registry.fill(HIST("hCollCounter"), static_cast(EvtSel::Sel8Full), weightEvt); // Coll+TVX+NoTFB+NoITSROFB+NoPileup + if (zvtx) { + registry.fill(HIST("hCollCounter"), static_cast(EvtSel::Sel8FullZvtx), weightEvt); // Coll+TVX+NoTFB+NoITSROFB+NoPileup+Zvtx + } + if (jetderiveddatautilities::selectCollision(collision, eventSelectionBitsSel8FullGood)) { + evtselCode |= EvtSelFlag::kIsGoodZvtx; + registry.fill(HIST("hCollCounter"), static_cast(EvtSel::Sel8FullGood), weightEvt); // Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx + if (zvtx) { + registry.fill(HIST("hCollCounter"), static_cast(EvtSel::Sel8FullGoodZvtx), weightEvt); // Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx+Zvtx + } + } + } + } + } + } + + return evtselCode; + } + + // Same stage ladder as fillCollCounter(), but evaluated directly on the MC-truth BC's raw evsel bits + // (mcCollision.bc_as(), a mandatory 1:1 index -> always resolves, no "has BC" gate needed) + // instead of a matched reco collision's packed EventSel word. Reuses EvtSelFlag::kColl/EvtSel::Coll(+Zvtx) + // to mean "BC" here (a BC always exists for an McCollision, so that stage is trivially satisfied) purely + // so this cascade lines up bin-for-bin with hCollCounter/hMcCollCounter. + template + EvtSelFlag fillMcCollCounterBC(AnalysisBC const& bc, float posZ, float weightEvt = 1.f) + { + EvtSelFlag evtselCode = EvtSelFlag::Coll; // "Coll" == BC here, see above + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Coll), weightEvt); // BC + + bool zvtx = std::fabs(posZ) < vertexZCut; + if (zvtx) { + evtselCode |= EvtSelFlag::kZvtx; + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::CollZvtx), weightEvt); // BC+Zvtx + } + + if (bc.selection_bit(o2::aod::evsel::kIsTriggerTVX)) { + evtselCode |= EvtSelFlag::kTVX; + registry.fill(HIST("h_vertexZ_truth_bctvx"), posZ, weightEvt); + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::TVX), weightEvt); // BC+TVX + if (zvtx) { + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::TVXZvtx), weightEvt); // BC+TVX, +Zvtx + } + if (bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { + evtselCode |= EvtSelFlag::kNoTFB; + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::SelMC), weightEvt); // BC+TVX+NoTFB + if (zvtx) { + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::SelMCZvtx), weightEvt); // BC+TVX+NoTFB, +Zvtx + } + if (bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + evtselCode |= EvtSelFlag::kNoITSROFB; + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Sel8), weightEvt); // BC+TVX+NoTFB+NoITSROFB + if (zvtx) { + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Sel8Zvtx), weightEvt); // BC+TVX+NoTFB+NoITSROFB, +Zvtx + } + if (bc.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + evtselCode |= EvtSelFlag::kNoPileup; + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Sel8Full), weightEvt); // BC+TVX+NoTFB+NoITSROFB+NoPileup + if (zvtx) { + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Sel8FullZvtx), weightEvt); // BC+TVX+NoTFB+NoITSROFB+NoPileup, +Zvtx + } + if (bc.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + evtselCode |= EvtSelFlag::kIsGoodZvtx; + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Sel8FullGood), weightEvt); // BC+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx + if (zvtx) { + registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Sel8FullGoodZvtx), weightEvt); // BC+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx, +Zvtx + } + } + } } } } @@ -771,25 +1246,60 @@ 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) + // Shared body of processDataJets(Sub)()/processDataJetsSV(Sub)(); see processMCDJetsCommon() for the + // withSV/withSub convention. + template + void processDataJetsCommon(FilteredCollisions::iterator const& collision, AnalysisJets const& alljets, FilteredTracks const& allTracks, SecondaryVertices const* allSVs = nullptr) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { return; } + // Keeps hCollCounter populated consistently across every process[Data,MCD]Jets(SV)(Sub) variant, not + // just the Sel ones - these should all behave identically except for which jet histograms get filled. + fillCollCounter(collision); registry.fill(HIST("h_event_counter"), 0.0); // Coll+TVX+NoTFB+NoITSROFB+Zvtx registry.fill(HIST("h_vertexZ"), collision.posZ()); + float rho = collision.rho(); + for (const auto& analysisJet : alljets) { if (!isAcceptedJet(analysisJet)) { continue; } - fillDataJetHistograms(analysisJet, allTracks); + if constexpr (withSV) { + fillDataJetHistogramsSV(analysisJet, allTracks, *allSVs, rho); + } else { + fillDataJetHistograms(analysisJet, allTracks, rho); + } } } + + void processDataJets(FilteredCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& allTracks) + { + processDataJetsCommon(collision, alljets, allTracks); + } PROCESS_SWITCH(BjetTaggingGnn, processDataJets, "jet information in Data", false); + void processDataJetsSub(FilteredCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& allTracks) + { + processDataJetsCommon(collision, alljets, allTracks); + } + PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSub, "jet information in Data, UE(rho*area)-subtracted jet pT", false); + + void processDataJetsSV(FilteredCollisions::iterator const& collision, FilteredDataJetsSV const& alljets, aod::DataSecondaryVertex3Prongs const& allSVs, FilteredTracks const& allTracks) + { + processDataJetsCommon(collision, alljets, allTracks, &allSVs); + } + PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSV, "jet information in Data with secondary vertex info", false); + + void processDataJetsSVSub(FilteredCollisions::iterator const& collision, FilteredDataJetsSV const& alljets, aod::DataSecondaryVertex3Prongs const& allSVs, FilteredTracks const& allTracks) + { + processDataJetsCommon(collision, alljets, allTracks, &allSVs); + } + PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSVSub, "jet information in Data with secondary vertex info, UE(rho*area)-subtracted jet pT", false); + 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 @@ -808,16 +1318,48 @@ struct BjetTaggingGnn { registry.fill(HIST("h_vertexZ"), collision.posZ()); + float rho = collision.rho(); + for (const auto& analysisJet : alljets) { if (!isAcceptedJet(analysisJet)) { continue; } - fillDataJetHistograms(analysisJet, allTracks); + fillDataJetHistograms(analysisJet, allTracks, rho); } } PROCESS_SWITCH(BjetTaggingGnn, processDataJetsTrig, "jet information in software triggered Data", false); + void processDataJetsTrigSub(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 + const auto& bc = collision.bc_as(); + // Initialize CCDB objects using the BC info + initCCDB(bc); + // If SoftwareTriggerSelection (i.e. skimming) is enabled, skip this event unless it passes Zorro selection + if (doSoftwareTriggerSelection && !zorro.isSelected(bc.globalBC())) { + return; + } + + if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { + return; + } + registry.fill(HIST("h_event_counter"), 0.0); // Coll+TVX+NoTFB+NoITSROFB+Zvtx + + registry.fill(HIST("h_vertexZ"), collision.posZ()); + + float rho = collision.rho(); + + for (const auto& analysisJet : alljets) { + if (!isAcceptedJet(analysisJet)) { + continue; + } + + fillDataJetHistograms(analysisJet, allTracks, rho); + } + } + PROCESS_SWITCH(BjetTaggingGnn, processDataJetsTrigSub, "jet information in software triggered Data, UE(rho*area)-subtracted jet pT", false); + void processDataJetsSel(AnalysisCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& /*allTracks*/) { EvtSelFlag evtselCode = fillCollCounter(collision); @@ -827,18 +1369,28 @@ struct BjetTaggingGnn { continue; } - registry.fill(HIST("h_jetpT_coll"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::Coll)); - registry.fill(HIST("h_jetpT_coll_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::CollZvtx)); - registry.fill(HIST("h_jetpT_tvx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::TVX)); - registry.fill(HIST("h_jetpT_tvx_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx)); - registry.fill(HIST("h_jetpT_selmc"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC)); - 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)); + fillEvtSelStages(HIST("h2_jetpT_evtsel"), analysisJet.pt(), evtselCode, 1.0, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); } } PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSel, "jet information in Data (event selection)", false); + void processDataJetsSelSub(AnalysisCollisions::iterator const& collision, FilteredDataJets const& alljets, FilteredTracks const& /*allTracks*/) + { + EvtSelFlag evtselCode = fillCollCounter(collision); + + float rho = collision.rho(); + + for (const auto& analysisJet : alljets) { + if (!isAcceptedJet(analysisJet)) { + continue; + } + + float jetptSub = analysisJet.pt() - rho * analysisJet.area(); + fillEvtSelStages(HIST("h2_jetpT_evtsel_sub"), jetptSub, evtselCode, 1.0, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + } + } + PROCESS_SWITCH(BjetTaggingGnn, processDataJetsSelSub, "jet information in Data (event selection), UE(rho*area)-subtracted jet pT", false); + void processDataTracks(FilteredCollisions::iterator const& collision, AnalysisTracks const& tracks) { if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { @@ -855,7 +1407,7 @@ struct BjetTaggingGnn { if (track.pt() >= trackPtMin) { registry.fill(HIST("h_dcaXY"), std::fabs(track.dcaXY())); registry.fill(HIST("h_dcaZ"), std::fabs(track.dcaZ())); - registry.fill(HIST("h3_dca_pt"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY())); + registry.fill(HIST("hSparse_dca_pt"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY())); } } } @@ -863,30 +1415,34 @@ struct BjetTaggingGnn { Preslice mcpjetsPerMCPCollision = aod::jet::mcCollisionId; - void processMCDJets(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& allTracks, FilteredMCPJets const& MCPjets, aod::JetParticles const& /*mcParticles*/, AnalysisCollisionsMCP const& /*mcCollisions*/) + // Shared body of processMCDJets(Sub)()/processMCDJetsSV(Sub)(): identical except for whether each jet's + // histograms are filled via fillMCDJetHistograms() or fillMCDJetHistogramsSV(), and whether raw or + // UE(rho*area)-subtracted jet pT (`withSub`) + the "_sub"-suffixed names are used throughout, including the + // matched geometric-response and pTHat-sparse histograms below, not just the core Db/flavor spectra. + // `withSV` selects the branch at compile time via `if constexpr` (bool-template-flag + defaulted-pointer + // pattern, matches PWGHF/D2H/Tasks/taskCharmPolarisation.cxx's `WithEp`/`QVecs*`). + template + void processMCDJetsCommon(FilteredCollisionsMCD::iterator const& collision, AnalysisJets const& MCDjets, FilteredTracksMCD const& allTracks, FilteredMCPJets const& MCPjets, SecondaryVertices const* allSVs = nullptr) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { - return; - } - // Reject outlier MC collisions - if (collision.isOutlier()) { + if (isRejectedMCDCollision(collision)) { return; } bool matchedMcColl = collision.has_mcCollision(); - if (matchedMcColl && collision.template mcCollision_as().isOutlier()) { - return; - } - // Uses only collisionId % trainingDatasetRaioParam != 0 for evaluation dataset - if (trainingDatasetRatioParam && collision.collisionId() % trainingDatasetRatioParam == 0) { - return; - } + float rho = collision.rho(); float weightEvt = useEventWeight ? collision.weight() : 1.f; + // Keeps hCollCounter populated consistently across every process[Data,MCD]Jets(SV)(Sub) variant, not + // just the Sel ones - these should all behave identically except for which jet histograms get filled. + fillCollCounter(collision, weightEvt); registry.fill(HIST("h_event_counter"), 0.0, weightEvt); registry.fill(HIST("h_vertexZ"), collision.posZ(), weightEvt); + // Only valid once matchedMcColl is confirmed true below; computed eagerly here (guarded by the ternary) + // since it's needed inside the MCDjets loop's matched-mcpjet block, not just after it. + float rhoMc = matchedMcColl ? collision.template mcCollision_as().rho() : 0.f; + // Store matched particle jet indices to avoid double-counting in mcpjets loop std::unordered_set matchedMcpJetIndices; @@ -895,12 +1451,20 @@ struct BjetTaggingGnn { continue; } - int8_t jetFlavor = fillMCDJetHistograms(analysisJet, allTracks, weightEvt); + int8_t jetFlavor; + if constexpr (withSV) { + jetFlavor = fillMCDJetHistogramsSV(analysisJet, allTracks, *allSVs, rho, weightEvt); + } else { + jetFlavor = fillMCDJetHistograms(analysisJet, allTracks, rho, weightEvt).jetFlavor; + } if (!matchedMcColl) { continue; } + float jetpT = withSub ? (analysisJet.pt() - rho * analysisJet.area()) : analysisJet.pt(); bool matchedJet = false; + // Everything below is jetpT/mcpJetpT-dependent, so the whole block is duplicated once per `withSub` + // with the "_sub"-suffixed HIST names. for (const auto& mcpjet : analysisJet.template matchedJetGeo_as()) { // matchedJetGeo_as is not Filtered. if (mcpjet.pt() < jetPtMin || mcpjet.pt() >= jetPtMax || mcpjet.eta() >= jetEtaMax - mcpjet.r() / 100.f || mcpjet.eta() <= jetEtaMin + mcpjet.r() / 100.f) { @@ -908,25 +1472,46 @@ struct BjetTaggingGnn { } matchedJet = true; matchedMcpJetIndices.insert(mcpjet.globalIndex()); - 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); - registry.fill(HIST("h3_pthat_jetpT"), collision.template mcCollision_as().ptHard(), analysisJet.pt(), mcpjet.pt(), weightEvt); // Matched jets - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_b"), analysisJet.pt(), mcpjet.pt(), weightEvt); - registry.fill(HIST("h_jetpT_b_matched"), analysisJet.pt(), weightEvt); - registry.fill(HIST("h_jetpT_particle_b_matched"), mcpjet.pt(), weightEvt); - registry.fill(HIST("h3_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), analysisJet.pt(), mcpjet.pt(), weightEvt); // Matched b-jets - } else if (jetFlavor == JetTaggingSpecies::charm) { - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_c"), analysisJet.pt(), mcpjet.pt(), weightEvt); + float mcpJetpT = withSub ? (mcpjet.pt() - rhoMc * mcpjet.area()) : mcpjet.pt(); + if constexpr (withSub) { + registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_sub"), jetpT, mcpJetpT, weightEvt); + registry.fill(HIST("h_jetpT_matched_sub"), jetpT, weightEvt); + registry.fill(HIST("h_jetpT_particle_matched_sub"), mcpJetpT, weightEvt); + registry.fill(HIST("hSparse_pthat_jetpT_sub"), collision.template mcCollision_as().ptHard(), jetpT, mcpJetpT, weightEvt); // Matched jets + registry.fill(HIST("h2_jetpT_matched_flavor_sub"), jetpT, getJetFlavorCat(jetFlavor), weightEvt); + registry.fill(HIST("h2_jetpT_particle_matched_flavor_sub"), mcpJetpT, getJetFlavorCat(jetFlavor), weightEvt); + registry.fill(HIST("h3_Response_DetjetpT_PartjetpT_flavor_sub"), jetpT, mcpJetpT, getJetFlavorCat(jetFlavor), weightEvt); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_pthat_jetpT_b_sub"), collision.template mcCollision_as().ptHard(), jetpT, mcpJetpT, weightEvt); // Matched b-jets + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_pthat_jetpT_c_sub"), collision.template mcCollision_as().ptHard(), jetpT, mcpJetpT, weightEvt); // Matched c-jets + } } else { - registry.fill(HIST("h2_Response_DetjetpT_PartjetpT_lf"), analysisJet.pt(), mcpjet.pt(), weightEvt); + registry.fill(HIST("h2_Response_DetjetpT_PartjetpT"), jetpT, mcpJetpT, weightEvt); + registry.fill(HIST("h_jetpT_matched"), jetpT, weightEvt); + registry.fill(HIST("h_jetpT_particle_matched"), mcpJetpT, weightEvt); + registry.fill(HIST("hSparse_pthat_jetpT"), collision.template mcCollision_as().ptHard(), jetpT, mcpJetpT, weightEvt); // Matched jets + registry.fill(HIST("h2_jetpT_matched_flavor"), jetpT, getJetFlavorCat(jetFlavor), weightEvt); + registry.fill(HIST("h2_jetpT_particle_matched_flavor"), mcpJetpT, getJetFlavorCat(jetFlavor), weightEvt); + registry.fill(HIST("h3_Response_DetjetpT_PartjetpT_flavor"), jetpT, mcpJetpT, getJetFlavorCat(jetFlavor), weightEvt); + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), jetpT, mcpJetpT, weightEvt); // Matched b-jets + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_pthat_jetpT_c"), collision.template mcCollision_as().ptHard(), jetpT, mcpJetpT, weightEvt); // Matched c-jets + } } } if (!matchedJet) { - registry.fill(HIST("h3_pthat_jetpT"), collision.template mcCollision_as().ptHard(), analysisJet.pt(), -1.f, weightEvt); // Fake jets, overflow-pTpart jets - if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h3_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), analysisJet.pt(), -1.f, weightEvt); // Overflow-pTpart b-jets + if constexpr (withSub) { + registry.fill(HIST("hSparse_pthat_jetpT_sub"), collision.template mcCollision_as().ptHard(), jetpT, -1.f, weightEvt); // Fake jets, overflow-pTpart jets + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_pthat_jetpT_b_sub"), collision.template mcCollision_as().ptHard(), jetpT, -1.f, weightEvt); // Overflow-pTpart b-jets + } + } else { + registry.fill(HIST("hSparse_pthat_jetpT"), collision.template mcCollision_as().ptHard(), jetpT, -1.f, weightEvt); // Fake jets, overflow-pTpart jets + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), jetpT, -1.f, weightEvt); // Overflow-pTpart b-jets + } } } } @@ -938,30 +1523,91 @@ struct BjetTaggingGnn { // Fill histograms for jets matched to the analysis event selection const auto& mcpjetspermcpcollision = MCPjets.sliceBy(mcpjetsPerMCPCollision, collision.mcCollisionId()); for (const auto& mcpjet : mcpjetspermcpcollision) { - registry.fill(HIST("h_jetpT_particle"), mcpjet.pt(), weightEvt); - - // Fill h3_pthat_jetpT only for unmatched particle jets (reco pT = -1) - if (matchedMcpJetIndices.find(mcpjet.globalIndex()) == matchedMcpJetIndices.end()) { - registry.fill(HIST("h3_pthat_jetpT"), collision.template mcCollision_as().ptHard(), -1.f, mcpjet.pt(), weightEvt); // Missing jets, overflow-pTreco jets - if (mcpjet.origin() == JetTaggingSpecies::beauty) { - registry.fill(HIST("h3_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), -1.f, mcpjet.pt(), weightEvt); // Missing b-jets, overflow-pTpart b-jets + float mcpJetpT = withSub ? (mcpjet.pt() - rhoMc * mcpjet.area()) : mcpjet.pt(); + int8_t jetFlavor = mcpjet.origin(); + if constexpr (withSub) { + registry.fill(HIST("h_jetpT_particle_sub"), mcpJetpT, weightEvt); + // Fill hSparse_pthat_jetpT only for unmatched particle jets (reco pT = -1) + if (matchedMcpJetIndices.find(mcpjet.globalIndex()) == matchedMcpJetIndices.end()) { + registry.fill(HIST("hSparse_pthat_jetpT_sub"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing jets, overflow-pTreco jets + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_pthat_jetpT_b_sub"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing b-jets, overflow-pTpart b-jets + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_pthat_jetpT_c_sub"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing c-jets, overflow-pTpart c-jets + } } + registry.fill(HIST("h2_jetpT_particle_flavor_sub"), mcpJetpT, getJetFlavorCat(jetFlavor), weightEvt); + } else { + registry.fill(HIST("h_jetpT_particle"), mcpJetpT, weightEvt); + // Fill hSparse_pthat_jetpT only for unmatched particle jets (reco pT = -1) + if (matchedMcpJetIndices.find(mcpjet.globalIndex()) == matchedMcpJetIndices.end()) { + registry.fill(HIST("hSparse_pthat_jetpT"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing jets, overflow-pTreco jets + if (jetFlavor == JetTaggingSpecies::beauty) { + registry.fill(HIST("hSparse_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing b-jets, overflow-pTpart b-jets + } else if (jetFlavor == JetTaggingSpecies::charm) { + registry.fill(HIST("hSparse_pthat_jetpT_c"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing c-jets, overflow-pTpart c-jets + } + } + registry.fill(HIST("h2_jetpT_particle_flavor"), mcpJetpT, getJetFlavorCat(jetFlavor), weightEvt); } + } + } - int8_t jetFlavor = mcpjet.origin(); + void processMCDJets(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& allTracks, FilteredMCPJets const& MCPjets, aod::JetParticles const& /*mcParticles*/, AnalysisCollisionsMCP const& /*mcCollisions*/) + { + processMCDJetsCommon(collision, MCDjets, allTracks, MCPjets); + } + PROCESS_SWITCH(BjetTaggingGnn, processMCDJets, "jet information in MC", false); + + void processMCDJetsSub(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& allTracks, FilteredMCPJets const& MCPjets, aod::JetParticles const& /*mcParticles*/, AnalysisCollisionsMCP const& /*mcCollisions*/) + { + processMCDJetsCommon(collision, MCDjets, allTracks, MCPjets); + } + PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSub, "jet information in MC, UE(rho*area)-subtracted jet pT", false); + + void processMCDJetsSV(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJetsSV const& MCDjets, aod::MCDSecondaryVertex3Prongs const& allSVs, FilteredTracksMCD const& allTracks, FilteredMCPJets const& MCPjets, aod::JetParticles const& /*mcParticles*/, AnalysisCollisionsMCP const& /*mcCollisions*/) + { + processMCDJetsCommon(collision, MCDjets, allTracks, MCPjets, &allSVs); + } + PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSV, "jet information in MC with secondary vertex info", false); + + void processMCDJetsSVSub(FilteredCollisionsMCD::iterator const& collision, FilteredMCDJetsSV const& MCDjets, aod::MCDSecondaryVertex3Prongs const& allSVs, FilteredTracksMCD const& allTracks, FilteredMCPJets const& MCPjets, aod::JetParticles const& /*mcParticles*/, AnalysisCollisionsMCP const& /*mcCollisions*/) + { + processMCDJetsCommon(collision, MCDjets, allTracks, MCPjets, &allSVs); + } + PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSVSub, "jet information in MC with secondary vertex info, UE(rho*area)-subtracted jet pT", false); + + void processMCDJetsSel(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& /*allTracks*/, FilteredMCPJets const& /*MCPjets*/, AnalysisCollisionsMCP 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); + for (const auto& analysisJet : MCDjets) { + if (!isAcceptedJet(analysisJet)) { + continue; + } + + int8_t jetFlavor = analysisJet.origin(); + + fillEvtSelStages(HIST("h2_jetpT_evtsel"), analysisJet.pt(), evtselCode, weightEvt, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); if (jetFlavor == JetTaggingSpecies::beauty) { - registry.fill(HIST("h_jetpT_particle_b"), mcpjet.pt(), weightEvt); + fillEvtSelStages(HIST("h2_jetpT_evtsel_b"), analysisJet.pt(), evtselCode, weightEvt, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); } 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); + fillEvtSelStages(HIST("h2_jetpT_evtsel_c"), analysisJet.pt(), evtselCode, weightEvt, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); } } } - PROCESS_SWITCH(BjetTaggingGnn, processMCDJets, "jet information in MC", false); + PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSel, "jet information in MC (event selection)", false); - void processMCDJetsSel(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& /*allTracks*/, FilteredMCPJets const& /*MCPjets*/, AnalysisCollisionsMCP const& /*mcCollisions*/) + void processMCDJetsSelSub(AnalysisCollisionsMCD::iterator const& collision, FilteredMCDJets const& MCDjets, FilteredTracksMCD const& /*allTracks*/, FilteredMCPJets const& /*MCPjets*/, AnalysisCollisionsMCP const& /*mcCollisions*/) { // Reject outlier MC collisions if (collision.isOutlier()) { @@ -974,38 +1620,33 @@ struct BjetTaggingGnn { float weightEvt = useEventWeight ? collision.weight() : 1.f; EvtSelFlag evtselCode = fillCollCounter(collision, weightEvt); + float rho = collision.rho(); + for (const auto& analysisJet : MCDjets) { if (!isAcceptedJet(analysisJet)) { continue; } int8_t jetFlavor = analysisJet.origin(); + float jetptSub = analysisJet.pt() - rho * analysisJet.area(); - 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); - registry.fill(HIST("h_jetpT_tvx_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_selmc"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); - 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); + fillEvtSelStages(HIST("h2_jetpT_evtsel_sub"), jetptSub, evtselCode, weightEvt, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); 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); - registry.fill(HIST("h_jetpT_b_tvx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::TVX) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_b_tvx_zvtx"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::TVXZvtx) ? weightEvt : 0.0); - registry.fill(HIST("h_jetpT_b_selmc"), analysisJet.pt(), hasAll(evtselCode, EvtSelFlag::SelMC) ? weightEvt : 0.0); - 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); + fillEvtSelStages(HIST("h2_jetpT_evtsel_b_sub"), jetptSub, evtselCode, weightEvt, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); + } else if (jetFlavor == JetTaggingSpecies::charm) { + fillEvtSelStages(HIST("h2_jetpT_evtsel_c_sub"), jetptSub, evtselCode, weightEvt, kEvtSelStageRecoFirst, kEvtSelStageRecoLast); } } } - PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSel, "jet information in MC (event selection)", false); + PROCESS_SWITCH(BjetTaggingGnn, processMCDJetsSelSub, "jet information in MC (event selection), UE(rho*area)-subtracted jet pT", false); PresliceUnsorted collisionsPerMCPCollision = aod::jmccollisionlb::mcCollisionId; - void processMCPJets(AnalysisCollisionsMCP const& mcCollisions, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, aod::JetParticles const& mcParticles) + // `withSub` selects UE(rho*area)-subtracted mcp-jet pT + the "_sub"-suffixed evtsel histogram names below + // instead of raw pT + the unsuffixed names. hMcCollCounter/h_vertexZ_truth*/h_event_counter_mcp aren't + // jet-pT-dependent and always use the plain name. + template + void processMCPJetsCommon(AnalysisCollisionsMCP const& mcCollisions, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, aod::JBCs const&) { // Subscribing AnalysisCollisionsMCP::iterator causes an issue related to unsorted JMcCollisionLbs index. for (const auto& mcCollision : mcCollisions) { @@ -1018,25 +1659,32 @@ struct BjetTaggingGnn { continue; } } + // mcCollision -> BC is a mandatory 1:1 index (jmccollision::JBCId), unlike the reco-collision match + // above (which can be 0, 1, or split into several matchedCollisions) - no slicing/size check needed. + const auto& matchedBC = mcCollision.bc_as(); + + float rho = mcCollision.rho(); float weightEvt = useEventWeight ? mcCollision.weight() : 1.f; EvtSelFlag evtselCode = EvtSelFlag::INEL; + EvtSelFlag evtselBCCode = 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 isTrueINELgt0 = isTrueINEL0(mcCollision, mcParticles); + // 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 - } + // if (isTrueINELgt0) { + // evtselCode |= EvtSelFlag::kINELgt0; + // registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELgt0), weightEvt); // INELgt0 + // } } + // Coll if (matchedCollisions.size() >= 1) { zvtxMatched = std::fabs(matchedCollisions.begin().posZ()) < vertexZCut; evtselCode |= EvtSelFlag::kColl; @@ -1068,6 +1716,20 @@ struct BjetTaggingGnn { // registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::INELgt0rec), weightEvt); // INELgt0+Zvtx(rec) // } } + if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsSel8Full)) { + evtselCode |= EvtSelFlag::kNoPileup; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8Full), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup) + if (zvtxMatched) { + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8FullZvtx), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup+Zvtx) + } + if (jetderiveddatautilities::selectCollision(matchedCollisions.begin(), eventSelectionBitsSel8FullGood)) { + evtselCode |= EvtSelFlag::kIsGoodZvtx; + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8FullGood), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx) + if (zvtxMatched) { + registry.fill(HIST("hMcCollCounter"), static_cast(EvtSel::Sel8FullGoodZvtx), weightEvt); // McColl(-> Coll+TVX+NoTFB+NoITSROFB+NoPileup+GoodZvtx+Zvtx) + } + } + } } } } @@ -1076,71 +1738,68 @@ struct BjetTaggingGnn { } } + // BC (mcCollision's own truth-level BC, always resolved - see fillMcCollCounterBC()). + evtselBCCode |= fillMcCollCounterBC(matchedBC, mcCollision.posZ(), weightEvt); + const auto& mcpjetspermcpcollision = mcpjets.sliceBy(mcpjetsPerMCPCollision, mcCollision.mcCollisionId()); for (const auto& mcpjet : mcpjetspermcpcollision) { int8_t jetFlavor = mcpjet.origin(); + float mcpJetpT = withSub ? (mcpjet.pt() - rho * mcpjet.area()) : mcpjet.pt(); + + if constexpr (withSub) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtsel_sub"), mcpJetpT, evtselCode, weightEvt, 0, kEvtSelStages.size() - 1); + if (jetFlavor == JetTaggingSpecies::beauty) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtsel_b_sub"), mcpJetpT, evtselCode, weightEvt, 0, kEvtSelStages.size() - 1); + } else if (jetFlavor == JetTaggingSpecies::charm) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtsel_c_sub"), mcpJetpT, evtselCode, weightEvt, 0, kEvtSelStages.size() - 1); + } - 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_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_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); + fillEvtSelStages(HIST("h2_jetpT_particle_evtselbc_sub"), mcpJetpT, evtselBCCode, weightEvt, 0, kEvtSelStages.size() - 1); + if (jetFlavor == JetTaggingSpecies::beauty) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtselbc_b_sub"), mcpJetpT, evtselBCCode, weightEvt, 0, kEvtSelStages.size() - 1); + } else if (jetFlavor == JetTaggingSpecies::charm) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtselbc_c_sub"), mcpJetpT, evtselBCCode, weightEvt, 0, kEvtSelStages.size() - 1); + } + } else { + fillEvtSelStages(HIST("h2_jetpT_particle_evtsel"), mcpJetpT, evtselCode, weightEvt, 0, kEvtSelStages.size() - 1); + if (jetFlavor == JetTaggingSpecies::beauty) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtsel_b"), mcpJetpT, evtselCode, weightEvt, 0, kEvtSelStages.size() - 1); + } else if (jetFlavor == JetTaggingSpecies::charm) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtsel_c"), mcpJetpT, evtselCode, weightEvt, 0, kEvtSelStages.size() - 1); + } + + fillEvtSelStages(HIST("h2_jetpT_particle_evtselbc"), mcpJetpT, evtselBCCode, weightEvt, 0, kEvtSelStages.size() - 1); + if (jetFlavor == JetTaggingSpecies::beauty) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtselbc_b"), mcpJetpT, evtselBCCode, weightEvt, 0, kEvtSelStages.size() - 1); + } else if (jetFlavor == JetTaggingSpecies::charm) { + fillEvtSelStages(HIST("h2_jetpT_particle_evtselbc_c"), mcpJetpT, evtselBCCode, weightEvt, 0, kEvtSelStages.size() - 1); + } } } } } + + void processMCPJets(AnalysisCollisionsMCP const& mcCollisions, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, aod::JBCs const& bcs) + { + processMCPJetsCommon(mcCollisions, mcpjets, collisions, bcs); + } PROCESS_SWITCH(BjetTaggingGnn, processMCPJets, "mc collision information", false); + void processMCPJetsSub(AnalysisCollisionsMCP const& mcCollisions, FilteredMCPJets const& mcpjets, AnalysisCollisionsMCD const& collisions, aod::JBCs const& bcs) + { + processMCPJetsCommon(mcCollisions, mcpjets, collisions, bcs); + } + PROCESS_SWITCH(BjetTaggingGnn, processMCPJetsSub, "mc collision information, UE(rho*area)-subtracted jet pT", false); + Preslice mcparticlesPerMCPCollision = aod::jmcparticle::mcCollisionId; void processMCDTracks(FilteredCollisionsMCD::iterator const& collision, AnalysisTracksMCD const& tracks, AnalysisCollisionsMCP const& /*mcCollisions*/, aod::JetParticles const& allParticles) { - if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) { - return; - } - // Reject outlier MC collisions - if (collision.isOutlier()) { + if (isRejectedMCDCollision(collision)) { return; } bool matchedMcColl = collision.has_mcCollision(); - if (matchedMcColl && 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; @@ -1154,14 +1813,14 @@ struct BjetTaggingGnn { if (track.pt() >= trackPtMin) { registry.fill(HIST("h_dcaXY"), std::fabs(track.dcaXY()), weightEvt); registry.fill(HIST("h_dcaZ"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } if (!matchedMcColl) { if (track.pt() >= trackPtMin) { registry.fill(HIST("h_dcaXY_coll_fake"), std::fabs(track.dcaXY()), weightEvt); registry.fill(HIST("h_dcaZ_coll_fake"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_coll_fake"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_coll_fake"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } continue; } @@ -1169,7 +1828,7 @@ struct BjetTaggingGnn { if (track.pt() >= trackPtMin) { registry.fill(HIST("h_dcaXY_fake"), std::fabs(track.dcaXY()), weightEvt); registry.fill(HIST("h_dcaZ_fake"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_fake"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_fake"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } continue; } @@ -1185,35 +1844,35 @@ struct BjetTaggingGnn { if (particle.mcCollisionId() == collision.mcCollisionId()) { registry.fill(HIST("h_dcaXY_coll_matched"), std::fabs(track.dcaXY()), weightEvt); // Matched to particle from the same MC collision registry.fill(HIST("h_dcaZ_coll_matched"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_coll_matched"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_coll_matched"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); int origin = RecoDecay::getParticleOrigin(allParticles, particle, false); if (origin == RecoDecay::OriginType::NonPrompt) { registry.fill(HIST("h_dcaXY_coll_matched_b"), std::fabs(track.dcaXY()), weightEvt); registry.fill(HIST("h_dcaZ_coll_matched_b"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_coll_matched_b"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_coll_matched_b"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } else if (origin == RecoDecay::OriginType::Prompt) { registry.fill(HIST("h_dcaXY_coll_matched_c"), std::fabs(track.dcaXY()), weightEvt); registry.fill(HIST("h_dcaZ_coll_matched_c"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_coll_matched_c"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_coll_matched_c"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } else { registry.fill(HIST("h_dcaXY_coll_matched_lf"), std::fabs(track.dcaXY()), weightEvt); registry.fill(HIST("h_dcaZ_coll_matched_lf"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_coll_matched_lf"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_coll_matched_lf"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } } else { registry.fill(HIST("h_dcaXY_coll_mismatched"), std::fabs(track.dcaXY()), weightEvt); // Matched to particle from a different MC collision registry.fill(HIST("h_dcaZ_coll_mismatched"), std::fabs(track.dcaZ()), weightEvt); - registry.fill(HIST("h3_dca_pt_coll_mismatched"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_coll_mismatched"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } } else { 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); - registry.fill(HIST("h3_dca_pt_npp"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_npp"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), 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); - registry.fill(HIST("h3_dca_pt_npp_mismatched"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); + registry.fill(HIST("hSparse_dca_pt_npp_mismatched"), track.pt(), std::fabs(track.dcaZ()), std::fabs(track.dcaXY()), weightEvt); } } } diff --git a/PWGJE/Tasks/bjetTreeCreator.cxx b/PWGJE/Tasks/bjetTreeCreator.cxx index 4adb20dc780..2b370f21d73 100644 --- a/PWGJE/Tasks/bjetTreeCreator.cxx +++ b/PWGJE/Tasks/bjetTreeCreator.cxx @@ -112,6 +112,8 @@ DECLARE_SOA_COLUMN(TrackTPCChi2NCl, tracktpcchi2ncl, float); //! The t DECLARE_SOA_COLUMN(TrackITSNCls, trackitsncls, float); //! The track ITS NCls DECLARE_SOA_COLUMN(TrackTPCNCls, tracktpcncls, float); //! The track TPC NCls (Found) DECLARE_SOA_COLUMN(TrackTPCNCrossedRows, tracktpcncrossedrows, float); //! The track TPC NCrossedRows +// DECLARE_SOA_COLUMN(TrackTPCNSigmaPi, tracktpcnsigmapi, float); //! The track TPC nSigma Pi +// DECLARE_SOA_COLUMN(TrackTOFNSigmaPi, tracktofnsigmapi, float); //! The track TOF nSigma Pi DECLARE_SOA_COLUMN(TrackOrigin, trk_origin, int); //! The track origin label for GNN track origin predictions DECLARE_SOA_COLUMN(TrackVtxIndex, trk_vtx_index, int); //! The track vertex index for GNN vertex predictions // DECLARE_SOA_COLUMN(DCATrackJet, dcatrackjet, float); //! The distance between track and jet, unfortunately it cannot be calculated in O2 @@ -158,6 +160,18 @@ DECLARE_SOA_TABLE(bjetTracksParamsExtrb, "AOD", "BJETTRACKSEXTRB", using bjetTracksParamExtrb = bjetTracksParamsExtrb::iterator; +// // PID information +// DECLARE_SOA_TABLE(bjetTracksParamsExtrc, "AOD", "BJETTRACKSEXTRC", +// // o2::soa::Index<>, +// trackInfo::TrackPhi, +// trackInfo::TrackCharge, +// trackInfo::TrackTPCNSigmaPi, +// trackInfo::TrackTOFNSigmaPi, +// trackInfo::TrackOrigin, +// trackInfo::TrackVtxIndex); + +// using bjetTracksParamExtrc = bjetTracksParamsExtrc::iterator; + namespace SVInfo { DECLARE_SOA_INDEX_COLUMN(bjetParam, jetindex); //! The jet index @@ -554,7 +568,10 @@ struct BJetTreeCreator { trkIdx++; - if (constituent.pt() < trackPtMin || !jettaggingutilities::trackAcceptanceWithDca(constituent, maxIPxy, maxIPz)) { + // No DCA-acceptance cut here: matches BjetTaggingGnn::fillMCDJetHistograms(), which only applies + // a track pT cut (trackPtMinGnn) for GNN input tracks, so the training tree stays consistent with + // what is evaluated at inference time. + if (constituent.pt() < trackPtMin) { continue; } From 2d5605e459024000c6abfc1f50bf2f80a5ec1f7b Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 19 Aug 2026 08:09:07 +0000 Subject: [PATCH 2/4] Please consider the following formatting changes --- PWGJE/Tasks/bjetTaggingGnn.cxx | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/PWGJE/Tasks/bjetTaggingGnn.cxx b/PWGJE/Tasks/bjetTaggingGnn.cxx index 92aef55088e..f601a8c54cb 100644 --- a/PWGJE/Tasks/bjetTaggingGnn.cxx +++ b/PWGJE/Tasks/bjetTaggingGnn.cxx @@ -142,21 +142,21 @@ struct EvtSelStage { const char* label; }; constexpr std::array kEvtSelStages{{{EvtSelFlag::INEL, "INEL"}, - {EvtSelFlag::INELZvtx, "INEL+Zvtx"}, - {EvtSelFlag::Coll, "Coll"}, - {EvtSelFlag::CollZvtx, "Coll+Zvtx"}, - {EvtSelFlag::TVX, "TVX"}, - {EvtSelFlag::TVXZvtx, "TVX+Zvtx"}, - {EvtSelFlag::SelMC, "SelMC"}, - {EvtSelFlag::SelMCZvtx, "SelMC+Zvtx"}, - {EvtSelFlag::Sel8, "Sel8"}, - {EvtSelFlag::Sel8Zvtx, "Sel8+Zvtx"}, - {EvtSelFlag::Sel8Full, "Sel8Full"}, - {EvtSelFlag::Sel8FullZvtx, "Sel8Full+Zvtx"}, - {EvtSelFlag::Sel8FullGood, "Sel8FullGood"}, - {EvtSelFlag::Sel8FullGoodZvtx, "Sel8FullGood+Zvtx"}}}; -constexpr int kEvtSelStageRecoFirst = 2; // Coll -constexpr int kEvtSelStageRecoLast = 13; // Sel8FullGood+Zvtx + {EvtSelFlag::INELZvtx, "INEL+Zvtx"}, + {EvtSelFlag::Coll, "Coll"}, + {EvtSelFlag::CollZvtx, "Coll+Zvtx"}, + {EvtSelFlag::TVX, "TVX"}, + {EvtSelFlag::TVXZvtx, "TVX+Zvtx"}, + {EvtSelFlag::SelMC, "SelMC"}, + {EvtSelFlag::SelMCZvtx, "SelMC+Zvtx"}, + {EvtSelFlag::Sel8, "Sel8"}, + {EvtSelFlag::Sel8Zvtx, "Sel8+Zvtx"}, + {EvtSelFlag::Sel8Full, "Sel8Full"}, + {EvtSelFlag::Sel8FullZvtx, "Sel8Full+Zvtx"}, + {EvtSelFlag::Sel8FullGood, "Sel8FullGood"}, + {EvtSelFlag::Sel8FullGoodZvtx, "Sel8FullGood+Zvtx"}}}; +constexpr int kEvtSelStageRecoFirst = 2; // Coll +constexpr int kEvtSelStageRecoLast = 13; // Sel8FullGood+Zvtx // Labels a full-range [0, kEvtSelStages.size()-1] evtsel axis. Bins below `first` (unreachable at // reco level) get numbered placeholder labels ("_1", "_2", ...) instead of a stage label, matching @@ -1182,7 +1182,7 @@ struct BjetTaggingGnn { template EvtSelFlag fillMcCollCounterBC(AnalysisBC const& bc, float posZ, float weightEvt = 1.f) { - EvtSelFlag evtselCode = EvtSelFlag::Coll; // "Coll" == BC here, see above + EvtSelFlag evtselCode = EvtSelFlag::Coll; // "Coll" == BC here, see above registry.fill(HIST("hMcCollCounterBC"), static_cast(EvtSel::Coll), weightEvt); // BC bool zvtx = std::fabs(posZ) < vertexZCut; From 07ae072e2c55173db8c2b2b3454eb63849ada43a Mon Sep 17 00:00:00 2001 From: Changhwan Choi Date: Wed, 19 Aug 2026 17:53:41 +0900 Subject: [PATCH 3/4] Resolved MegaLinter --- PWGJE/Tasks/bjetTaggingGnn.cxx | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/PWGJE/Tasks/bjetTaggingGnn.cxx b/PWGJE/Tasks/bjetTaggingGnn.cxx index f601a8c54cb..6c87c920b2e 100644 --- a/PWGJE/Tasks/bjetTaggingGnn.cxx +++ b/PWGJE/Tasks/bjetTaggingGnn.cxx @@ -27,6 +27,7 @@ #include "Common/DataModel/EventSelection.h" #include +#include #include #include #include @@ -237,7 +238,7 @@ struct BjetTaggingGnn { std::vector eventSelectionBitsSel8Full; std::vector eventSelectionBitsSel8FullGood; - int trackSelectionBits; + int trackSelectionBits{}; std::vector jetRadiiValues; @@ -393,7 +394,7 @@ struct BjetTaggingGnn { // h_jetEta/h_jetPhi/h_jetMass/h_Db/h2_nTracks_Db don't depend on jet pT - one shared copy regardless of // `withSub` (see addCoreJetHistograms()'s comment above). registry.add("h_jetEta", "", {HistType::kTH1F, {axisJetEta}}, callSumw2); - registry.add("h_jetPhi", "", {HistType::kTH1F, {{200, 0., 2. * M_PI, "#it{phi}_{jet}"}}}); + registry.add("h_jetPhi", "", {HistType::kTH1F, {{200, 0., o2::constants::math::TwoPI, "#it{phi}_{jet}"}}}); registry.add("h_jetMass", "", {HistType::kTH1F, {axisJetMass}}); registry.add("h_Db", "", {HistType::kTH1F, {axisDbFine}}); registry.add("h2_nTracks_Db", "", {HistType::kTH2F, {axisNTracks, axisDb}}); @@ -404,7 +405,7 @@ struct BjetTaggingGnn { } registry.add("h_gnnfeat_trackpT", "", {HistType::kTH1F, {{200, 0., 100., "#it{p}_{T} (GeV/#it{c})"}}}); - registry.add("h_gnnfeat_trackPhi", "", {HistType::kTH1F, {{200, 0., 2. * M_PI, "#it{#phi}"}}}); + registry.add("h_gnnfeat_trackPhi", "", {HistType::kTH1F, {{200, 0., o2::constants::math::TwoPI, "#it{#phi}"}}}); registry.add("h_gnnfeat_trackEta", "", {HistType::kTH1F, {{200, -0.9, 0.9, "#it{#eta}"}}}); registry.add("h_gnnfeat_trackCharge", "", {HistType::kTH1F, {{3, -1., 2., "#it{q}"}}}); registry.add("h_gnnfeat_trackDCAxy", "", {HistType::kTH1F, {{200, -5., 5., "DCA_{#it{xy}} (cm)"}}}); @@ -422,7 +423,7 @@ struct BjetTaggingGnn { if (doprocessDataTracks || doprocessMCDTracks) { registry.add("h_trackpT", "", {HistType::kTH1F, {axisTrackpT}}, callSumw2); registry.add("h_tracketa", "", {HistType::kTH1F, {{100, trackEtaMin, trackEtaMax, "#it{#eta}"}}}, callSumw2); - registry.add("h_trackphi", "", {HistType::kTH1F, {{100, 0.0, 2.0 * M_PI, "#it{#phi}"}}}, callSumw2); + registry.add("h_trackphi", "", {HistType::kTH1F, {{100, 0.0, o2::constants::math::TwoPI, "#it{#phi}"}}}, callSumw2); registry.add("h_dcaXY", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); registry.add("h_dcaZ", "", {HistType::kTH1F, {{200, 0., 4., "|DCA_{#it{z}}| (cm)"}}}, callSumw2); registry.add("hSparse_dca_pt", "", {HistType::kTHnSparseF, {{1000, 0., 100., "#it{p}_{T} (GeV/#it{c})"}, {200, 0., 4., "|DCA_{#it{z}}| (cm)"}, {200, 0., 4., "|DCA_{#it{xy}}| (cm)"}}}, callSumw2); @@ -703,7 +704,7 @@ struct BjetTaggingGnn { bool isAcceptedJet(AnalysisJet const& jet) { if (jetAreaFractionMin > largeNegativeNumber) { - if (jet.area() < jetAreaFractionMin * M_PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { + if (jet.area() < jetAreaFractionMin * o2::constants::math::PI * (jet.r() / 100.0) * (jet.r() / 100.0)) { return false; } } @@ -888,6 +889,7 @@ struct BjetTaggingGnn { // jetFlavor: JetTaggingSpecies of the jet; nTracks: number of GNN-input constituents (reused by fillMCDJetHistogramsSV()). struct JetHistFillResult { int8_t jetFlavor; + // cppcheck-suppress unusedStructMember int nTracks; }; From 3302581bd26c48389890d1d698ff7fdaa3ca64a7 Mon Sep 17 00:00:00 2001 From: Changhwan Choi Date: Fri, 21 Aug 2026 00:25:19 +0900 Subject: [PATCH 4/4] Resolved code checker issues --- PWGJE/Tasks/bjetTaggingGnn.cxx | 55 +++++++++++++++++---------------- PWGJE/Tasks/bjetTreeCreator.cxx | 4 +-- 2 files changed, 29 insertions(+), 30 deletions(-) diff --git a/PWGJE/Tasks/bjetTaggingGnn.cxx b/PWGJE/Tasks/bjetTaggingGnn.cxx index 6c87c920b2e..dec0b8acdb8 100644 --- a/PWGJE/Tasks/bjetTaggingGnn.cxx +++ b/PWGJE/Tasks/bjetTaggingGnn.cxx @@ -18,6 +18,7 @@ #include "PWGJE/Core/JetTaggingUtilities.h" #include "PWGJE/DataModel/Jet.h" #include "PWGJE/DataModel/JetReducedData.h" +#include "PWGJE/DataModel/JetSubtraction.h" #include "PWGJE/DataModel/JetTagging.h" #include "Common/CCDB/EventSelectionParams.h" @@ -39,6 +40,7 @@ #include #include +#include #include #include #include @@ -47,13 +49,10 @@ #include #include #include -#include #include #include #include -#include - using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -142,20 +141,20 @@ struct EvtSelStage { EvtSelFlag flag; const char* label; }; -constexpr std::array kEvtSelStages{{{EvtSelFlag::INEL, "INEL"}, - {EvtSelFlag::INELZvtx, "INEL+Zvtx"}, - {EvtSelFlag::Coll, "Coll"}, - {EvtSelFlag::CollZvtx, "Coll+Zvtx"}, - {EvtSelFlag::TVX, "TVX"}, - {EvtSelFlag::TVXZvtx, "TVX+Zvtx"}, - {EvtSelFlag::SelMC, "SelMC"}, - {EvtSelFlag::SelMCZvtx, "SelMC+Zvtx"}, - {EvtSelFlag::Sel8, "Sel8"}, - {EvtSelFlag::Sel8Zvtx, "Sel8+Zvtx"}, - {EvtSelFlag::Sel8Full, "Sel8Full"}, - {EvtSelFlag::Sel8FullZvtx, "Sel8Full+Zvtx"}, - {EvtSelFlag::Sel8FullGood, "Sel8FullGood"}, - {EvtSelFlag::Sel8FullGoodZvtx, "Sel8FullGood+Zvtx"}}}; +constexpr std::array kEvtSelStages{{{.flag = EvtSelFlag::INEL, .label = "INEL"}, + {.flag = EvtSelFlag::INELZvtx, .label = "INEL+Zvtx"}, + {.flag = EvtSelFlag::Coll, .label = "Coll"}, + {.flag = EvtSelFlag::CollZvtx, .label = "Coll+Zvtx"}, + {.flag = EvtSelFlag::TVX, .label = "TVX"}, + {.flag = EvtSelFlag::TVXZvtx, .label = "TVX+Zvtx"}, + {.flag = EvtSelFlag::SelMC, .label = "SelMC"}, + {.flag = EvtSelFlag::SelMCZvtx, .label = "SelMC+Zvtx"}, + {.flag = EvtSelFlag::Sel8, .label = "Sel8"}, + {.flag = EvtSelFlag::Sel8Zvtx, .label = "Sel8+Zvtx"}, + {.flag = EvtSelFlag::Sel8Full, .label = "Sel8Full"}, + {.flag = EvtSelFlag::Sel8FullZvtx, .label = "Sel8Full+Zvtx"}, + {.flag = EvtSelFlag::Sel8FullGood, .label = "Sel8FullGood"}, + {.flag = EvtSelFlag::Sel8FullGoodZvtx, .label = "Sel8FullGood+Zvtx"}}}; constexpr int kEvtSelStageRecoFirst = 2; // Coll constexpr int kEvtSelStageRecoLast = 13; // Sel8FullGood+Zvtx @@ -227,8 +226,8 @@ struct BjetTaggingGnn { OutputObj zorroSummary{"zorroSummary"}; // Service - Service pdg; - Service ccdb; + Service pdg{}; + Service ccdb{}; // Event selection bits std::vector eventSelectionBits; @@ -697,8 +696,8 @@ struct BjetTaggingGnn { return nTracks; } - const float largeNegativeNumber = -98.0f; - const float largePositiveNumber = 9999.0f; + static constexpr float largeNegativeNumber = -98.0f; + static constexpr float largePositiveNumber = 9999.0f; template bool isAcceptedJet(AnalysisJet const& jet) @@ -1081,18 +1080,20 @@ struct BjetTaggingGnn { } // Check if the collision is INEL>0 - const int nPartInel0 = 3; + static constexpr int nPartInel0 = 3; template bool isTrueINEL0(MCColl const& /*mccoll*/, MCPart const& mcparts) { for (const auto& mcparticle : mcparts) { - if (!mcparticle.isPhysicalPrimary()) + if (!mcparticle.isPhysicalPrimary()) { continue; + } const auto p = pdg->GetParticle(mcparticle.pdgCode()); if (p != nullptr) { if (std::abs(p->Charge()) >= nPartInel0) { - if (std::abs(mcparticle.eta()) < 1) + if (std::abs(mcparticle.eta()) < 1) { return true; + } } } } @@ -1453,7 +1454,7 @@ struct BjetTaggingGnn { continue; } - int8_t jetFlavor; + int8_t jetFlavor = 0; if constexpr (withSV) { jetFlavor = fillMCDJetHistogramsSV(analysisJet, allTracks, *allSVs, rho, weightEvt); } else { @@ -1530,7 +1531,7 @@ struct BjetTaggingGnn { if constexpr (withSub) { registry.fill(HIST("h_jetpT_particle_sub"), mcpJetpT, weightEvt); // Fill hSparse_pthat_jetpT only for unmatched particle jets (reco pT = -1) - if (matchedMcpJetIndices.find(mcpjet.globalIndex()) == matchedMcpJetIndices.end()) { + if (!matchedMcpJetIndices.contains(mcpjet.globalIndex())) { registry.fill(HIST("hSparse_pthat_jetpT_sub"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing jets, overflow-pTreco jets if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("hSparse_pthat_jetpT_b_sub"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing b-jets, overflow-pTpart b-jets @@ -1542,7 +1543,7 @@ struct BjetTaggingGnn { } else { registry.fill(HIST("h_jetpT_particle"), mcpJetpT, weightEvt); // Fill hSparse_pthat_jetpT only for unmatched particle jets (reco pT = -1) - if (matchedMcpJetIndices.find(mcpjet.globalIndex()) == matchedMcpJetIndices.end()) { + if (!matchedMcpJetIndices.contains(mcpjet.globalIndex())) { registry.fill(HIST("hSparse_pthat_jetpT"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing jets, overflow-pTreco jets if (jetFlavor == JetTaggingSpecies::beauty) { registry.fill(HIST("hSparse_pthat_jetpT_b"), collision.template mcCollision_as().ptHard(), -1.f, mcpJetpT, weightEvt); // Missing b-jets, overflow-pTpart b-jets diff --git a/PWGJE/Tasks/bjetTreeCreator.cxx b/PWGJE/Tasks/bjetTreeCreator.cxx index 2b370f21d73..52e203d91af 100644 --- a/PWGJE/Tasks/bjetTreeCreator.cxx +++ b/PWGJE/Tasks/bjetTreeCreator.cxx @@ -49,8 +49,6 @@ #include #include -#include - using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -811,7 +809,7 @@ struct BJetTreeCreator { } // Uses only collisionId % trainingDatasetRaioParam == 0 for training dataset - if (trainingDatasetRatioParam && collision.collisionId() % trainingDatasetRatioParam != 0) { + if (trainingDatasetRatioParam != 0 && collision.collisionId() % trainingDatasetRatioParam != 0) { return; }