diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index a87293c9cca..992ecd11714 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -436,7 +436,7 @@ if(FastJet_FOUND) o2physics_add_dpl_workflow(jet-cross-section-efficiency SOURCES jetCrossSectionEfficiency.cxx REUSE_FROM JetSubstructureHFPCH - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2::Framework O2::DataFormatsITSMFT O2Physics::PWGJECore O2Physics::AnalysisCore O2Physics::AnalysisCCDB COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(gamma-jet-tree-producer SOURCES gammaJetTreeProducer.cxx diff --git a/PWGJE/Tasks/jetCrossSectionEfficiency.cxx b/PWGJE/Tasks/jetCrossSectionEfficiency.cxx index 21efb950459..9bcc1195454 100644 --- a/PWGJE/Tasks/jetCrossSectionEfficiency.cxx +++ b/PWGJE/Tasks/jetCrossSectionEfficiency.cxx @@ -12,9 +12,10 @@ // Step-by-step event selection efficiency for MC particle-level jets. // Two ordering variants exposed as separate process functions: // - CollRecoFirst: reco collision required first, BC bits read from the reco-coll EvSel bitmask. -// - BcBitsFirst: BC bits read from the MC truth BC, then truth-side SBP, then reco at the end. +// - BcBitsFirst: BC bits read from the MC truth BC, then truth-side selections, then reco at the end. // /// \author Joonsuk Bae +/// \author Wooseok Ham #include "PWGJE/Core/JetDerivedDataUtilities.h" #include "PWGJE/Core/JetFindingUtilities.h" @@ -23,9 +24,15 @@ #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/DataModel/EventSelection.h" +#include +#include #include +#include #include +#include +#include #include #include #include @@ -33,12 +40,17 @@ #include #include +#include #include #include +#include #include +#include #include +#include #include +#include using namespace o2; using namespace o2::framework; @@ -48,7 +60,7 @@ struct JetCrossSectionEfficiency { HistogramRegistry registry; - Configurable eventSelections{"eventSelections", "sel8", "selTVX | selMC | selMCFull | sel8 | sel8Full"}; + Configurable eventSelections{"eventSelections", "sel8", "selTVX | selMC | selMCFull | sel8 | sel8Full | sel8FullPbPb"}; Configurable skipMBGapEvents{"skipMBGapEvents", false, "reject min-bias gap events from hybrid MB+JJ MC productions"}; Configurable vertexZCut{"vertexZCut", 10.0f, "Accepted z-vertex range"}; Configurable centralityMin{"centralityMin", -999.0f, "minimum centrality"}; @@ -77,12 +89,36 @@ struct JetCrossSectionEfficiency { Configurable applyRCT{"applyRCT", true, "apply RCT_pass step in the cascade (false: force the RCT step to pass)"}; Configurable rctSelectionsLabel{"rctSelectionsLabel", "CBT_hadronPID", "RCT selection preset name (see RCTSelectionFlags)"}; + // The reconstructed occupancy bits cannot be evaluated before hasColl. + // These configurables define the truth-only approximation used by sel8FullPbPb. + Configurable truthRofCcdbUrl{"truthRofCcdbUrl", "http://alice-ccdb.cern.ch", "CCDB URL used to retrieve ITS ROF parameters for sel8FullPbPb truth selections"}; + Configurable truthRofOffsetInBC{"truthRofOffsetInBC", -1, "ITS ROF bias in BC for sel8FullPbPb truth selections; -1 retrieves DPLAlpideParam from CCDB"}; + Configurable truthRofLengthInBC{"truthRofLengthInBC", -1, "ITS ROF length in BC for sel8FullPbPb truth selections; -1 retrieves DPLAlpideParam from CCDB"}; + Configurable truthTimeRangeStandardMinUs{"truthTimeRangeStandardMinUs", -4.0f, "truth lower time range in us; follows kNoCollInTimeRangeStandard"}; + Configurable truthTimeRangeStandardMaxUs{"truthTimeRangeStandardMaxUs", 2.0f, "truth upper time range in us; follows kNoCollInTimeRangeStandard"}; + Configurable truthTimeRangeNarrowUs{"truthTimeRangeNarrowUs", 0.25f, "truth narrow time range in us; follows kNoCollInTimeRangeStandard"}; + Configurable truthFT0CActivityMinForTrackProxy{"truthFT0CActivityMinForTrackProxy", 0.0f, "minimum multFT0C truth activity used as an ITS-track-presence proxy"}; + Configurable truthFT0CActivityThresholdTimeRange{"truthFT0CActivityThresholdTimeRange", 8000.0f, "multFT0C truth proxy threshold for NoCollInTimeRangeStandard; reco uses FT0C amplitude, so validate or tune this value"}; + Configurable truthFT0CActivityThresholdROF{"truthFT0CActivityThresholdROF", 5000.0f, "multFT0C truth proxy threshold for NoCollInRofStandard; reco uses FT0C amplitude, so validate or tune this value"}; + Configurable truthRofCloseVzMax{"truthRofCloseVzMax", 0.3f, "maximum truth |delta z| in cm for the NoCollInRofStandard track-presence proxy"}; + o2::aod::rctsel::RCTFlagsChecker rctChecker; uint64_t rctMask = 0; + Service ccdb{}; + + int cachedTruthRofRun = std::numeric_limits::min(); + int64_t cachedTruthRofOffsetInBC = -1; + int64_t cachedTruthRofLengthInBC = -1; bool applyTFB = true; bool applyROFB = true; bool applySBP = true; + bool applyNoCollInTimeRangeStandard = false; + bool applyNoCollInRofStandard = false; + bool isSel8FullPbPb = false; + + std::vector collRecoFirstLabels; + std::vector bcBitsFirstLabels; enum AcceptSplitCollisionsOptions { NonSplitOnly = 0, @@ -90,38 +126,36 @@ struct JetCrossSectionEfficiency { SplitOkCheckFirstAssocCollOnly // 2 }; + enum BinPbPbTruthSelectionOnly { + PbPbTruthSelectionOnlyInel = 1, + PbPbTruthSelectionOnlyRct = 2, + PbPbTruthSelectionOnlyTvx = 3, + PbPbTruthSelectionOnlyNoTimeFrameBorder = 4, + PbPbTruthSelectionOnlyNoItsRofBorder = 5, + PbPbTruthSelectionOnlySelection = 6, + PbPbTruthSelectionOnlyHasCollision = 7, + PbPbTruthSelectionOnlyVertexZ = 8, + PbPbTruthSelectionOnlyNoSplit = 9, + PbPbTruthSelectionOnlyNBins = 9 + }; + static constexpr float ConfigSwitchLow = -98.0f; static constexpr float ConfigSwitchHigh = 9998.0f; static constexpr float BrokenPtHardSentinel = 1.0f; + static constexpr int MinITSClustersForOccupancy = 5; - // CollRecoFirst: reco collision required first; BC bits read from the reco-coll EvSels. - enum BinCollRecoFirst { - CollRecoFirstInel = 1, - CollRecoFirstRct, - CollRecoFirstHasCollision, - CollRecoFirstVertexZ, - CollRecoFirstNoSplit, - CollRecoFirstTvx, - CollRecoFirstNoTimeFrameBorder, - CollRecoFirstNoItsRofBorder, - CollRecoFirstNoSameBunchPileup, - CollRecoFirstNBins = CollRecoFirstNoSameBunchPileup + struct TruthPbPbSelections { + bool valid = false; + bool noCollInTimeRangeStandard = false; + bool noCollInRofStandard = false; }; - // BcBitsFirst: BC bits read from the MC truth BC; SBP from a Preslice count - // (exactly one MC collision per truth BC) so it works before requiring reco. - enum BinBcBitsFirst { - BcBitsFirstInel = 1, - BcBitsFirstRct, - BcBitsFirstTvx, - BcBitsFirstNoTimeFrameBorder, - BcBitsFirstNoItsRofBorder, - BcBitsFirstNoSameBunchPileupTruth, - BcBitsFirstHasCollision, - BcBitsFirstVertexZ, - BcBitsFirstNoSplit, - BcBitsFirstNBins = BcBitsFirstNoSplit - }; + using JetCollisionsMCDWithParent = soa::Join; + using CollisionsWithEvSels = soa::Join; + using FullTracksIU = soa::Join; + + Partition pvTracks = ((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)); + Preslice pvTracksPerCollision = aod::track::collisionId; enum EventSelectionPreset { PresetSelTvx = 0, @@ -129,6 +163,7 @@ struct JetCrossSectionEfficiency { PresetSelMcFull, PresetSel8, PresetSel8Full, + PresetSel8FullPbPb, PresetInvalid }; @@ -149,6 +184,9 @@ struct JetCrossSectionEfficiency { if (preset == "sel8Full") { return PresetSel8Full; } + if (preset == "sel8FullPbPb") { + return PresetSel8FullPbPb; + } return PresetInvalid; } @@ -156,7 +194,9 @@ struct JetCrossSectionEfficiency { void init(InitContext&) { - if (!(acceptSplitCollisions == NonSplitOnly || acceptSplitCollisions == SplitOkCheckAnyAssocColl || acceptSplitCollisions == SplitOkCheckFirstAssocCollOnly)) { + if (acceptSplitCollisions != NonSplitOnly && + acceptSplitCollisions != SplitOkCheckAnyAssocColl && + acceptSplitCollisions != SplitOkCheckFirstAssocCollOnly) { LOGF(fatal, "Configurable acceptSplitCollisions has wrong input value; stopping workflow"); } @@ -189,73 +229,128 @@ struct JetCrossSectionEfficiency { applyROFB = true; applySBP = true; break; + case PresetSel8FullPbPb: + applyTFB = true; + applyROFB = true; + applySBP = false; + applyNoCollInTimeRangeStandard = true; + applyNoCollInRofStandard = true; + isSel8FullPbPb = true; + break; default: - LOGF(fatal, "Configurable eventSelections=%s not supported; use selTVX, selMC, selMCFull, sel8, or sel8Full", static_cast(eventSelections).c_str()); + LOGF(fatal, "Configurable eventSelections=%s not supported; use selTVX, selMC, selMCFull, sel8, sel8Full, or sel8FullPbPb", static_cast(eventSelections).c_str()); break; } + if (doprocessCrossSectionEfficiencyBcBitsFirst && isSel8FullPbPb && + (truthRofOffsetInBC < 0 || truthRofLengthInBC <= 0)) { + ccdb->setURL(static_cast(truthRofCcdbUrl)); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + } + + // Preserve the upstream binning and order for existing presets. + // sel8FullPbPb replaces the SBP stage with the two PbPb occupancy selections. + collRecoFirstLabels = {"INEL", "+RCT_pass", "+hasRecoColl", "+|zReco|<10", "+noSplit", "+kTVX", "+kNoTFB", "+kNoITSROFB"}; + bcBitsFirstLabels = {"INEL", "+RCT_pass", "+kTVX(truth)", "+kNoTFB(truth)", "+kNoITSROFB(truth)"}; + if (isSel8FullPbPb) { + collRecoFirstLabels.emplace_back("+kNoCollInTimeRangeStandard"); + collRecoFirstLabels.emplace_back("+kNoCollInRofStandard"); + bcBitsFirstLabels.emplace_back("+kNoCollInTimeRangeStandard(truth)"); + bcBitsFirstLabels.emplace_back("+kNoCollInRofStandard(truth)"); + bcBitsFirstLabels.emplace_back("+hasColl"); + bcBitsFirstLabels.emplace_back("+|zReco|<10"); + bcBitsFirstLabels.emplace_back("+noSplit"); + } else { + collRecoFirstLabels.emplace_back("+kNoSBP"); + bcBitsFirstLabels.emplace_back("+kNoSBP(truth)"); + bcBitsFirstLabels.emplace_back("+hasColl"); + bcBitsFirstLabels.emplace_back("+|zReco|<10"); + bcBitsFirstLabels.emplace_back("+noSplit"); + } + + auto setAxisLabels = [](TAxis* axis, const std::vector& labels) { + for (size_t i = 0; i < labels.size(); ++i) { + axis->SetBinLabel(static_cast(i + 1), labels[i].c_str()); + } + }; + AxisSpec jetPtAxis = {200, 0., jetPtMax, "#it{p}_{T} (GeV/#it{c})"}; if (doprocessCrossSectionEfficiency) { - AxisSpec axisSelectionCollRecoFirst = {CollRecoFirstNBins, 0.5, static_cast(CollRecoFirstNBins) + 0.5, "event selection (CollRecoFirst)"}; + AxisSpec axisSelectionCollRecoFirst = {static_cast(collRecoFirstLabels.size()), 0.5, static_cast(collRecoFirstLabels.size()) + 0.5, "event selection (CollRecoFirst)"}; registry.add("h2_jet_pt_part_eventselection_collRecoFirst", "part jet pT vs event selection (CollRecoFirst);#it{p}_{T,jet}^{part} (GeV/#it{c});event selection;counts", {HistType::kTH2F, {jetPtAxis, axisSelectionCollRecoFirst}}); - auto hJetPtCollRecoFirst = registry.get(HIST("h2_jet_pt_part_eventselection_collRecoFirst")); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstInel, "INEL"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstRct, "+RCT_pass"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstHasCollision, "+hasRecoColl"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstVertexZ, "+|zReco|<10"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstNoSplit, "+noSplit"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstTvx, "+kTVX"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstNoTimeFrameBorder, "+kNoTFB"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstNoItsRofBorder, "+kNoITSROFB"); - hJetPtCollRecoFirst->GetYaxis()->SetBinLabel(CollRecoFirstNoSameBunchPileup, "+kNoSBP"); + setAxisLabels(registry.get(HIST("h2_jet_pt_part_eventselection_collRecoFirst"))->GetYaxis(), collRecoFirstLabels); registry.add("h_mccollisions_eventselection_collRecoFirst", "number of mc events vs event selection (CollRecoFirst);event selection;entries", - {HistType::kTH1F, {{CollRecoFirstNBins, 0.5, static_cast(CollRecoFirstNBins) + 0.5}}}); - auto hMcCollisionsCollRecoFirst = registry.get(HIST("h_mccollisions_eventselection_collRecoFirst")); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstInel, "INEL"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstRct, "+RCT_pass"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstHasCollision, "+hasRecoColl"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstVertexZ, "+|zReco|<10"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstNoSplit, "+noSplit"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstTvx, "+kTVX"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstNoTimeFrameBorder, "+kNoTFB"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstNoItsRofBorder, "+kNoITSROFB"); - hMcCollisionsCollRecoFirst->GetXaxis()->SetBinLabel(CollRecoFirstNoSameBunchPileup, "+kNoSBP"); + {HistType::kTH1F, {axisSelectionCollRecoFirst}}); + setAxisLabels(registry.get(HIST("h_mccollisions_eventselection_collRecoFirst"))->GetXaxis(), collRecoFirstLabels); } if (doprocessCrossSectionEfficiencyBcBitsFirst) { - AxisSpec axisSelectionBcBitsFirst = {BcBitsFirstNBins, 0.5, static_cast(BcBitsFirstNBins) + 0.5, "event selection (BcBitsFirst)"}; + AxisSpec axisSelectionBcBitsFirst = {static_cast(bcBitsFirstLabels.size()), 0.5, static_cast(bcBitsFirstLabels.size()) + 0.5, "event selection (BcBitsFirst)"}; registry.add("h2_jet_pt_part_eventselection_bcBitsFirst", "part jet pT vs event selection (BcBitsFirst);#it{p}_{T,jet}^{part} (GeV/#it{c});event selection;counts", {HistType::kTH2F, {jetPtAxis, axisSelectionBcBitsFirst}}); - auto hJetPtBcBitsFirst = registry.get(HIST("h2_jet_pt_part_eventselection_bcBitsFirst")); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstInel, "INEL"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstRct, "+RCT_pass"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstTvx, "+kTVX(truth)"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstNoTimeFrameBorder, "+kNoTFB(truth)"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstNoItsRofBorder, "+kNoITSROFB(truth)"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstNoSameBunchPileupTruth, "+kNoSBP(truth)"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstHasCollision, "+hasColl"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstVertexZ, "+|zReco|<10"); - hJetPtBcBitsFirst->GetYaxis()->SetBinLabel(BcBitsFirstNoSplit, "+noSplit"); + setAxisLabels(registry.get(HIST("h2_jet_pt_part_eventselection_bcBitsFirst"))->GetYaxis(), bcBitsFirstLabels); registry.add("h_mccollisions_eventselection_bcBitsFirst", "number of mc events vs event selection (BcBitsFirst);event selection;entries", - {HistType::kTH1F, {{BcBitsFirstNBins, 0.5, static_cast(BcBitsFirstNBins) + 0.5}}}); - auto hMcCollisionsBcBitsFirst = registry.get(HIST("h_mccollisions_eventselection_bcBitsFirst")); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstInel, "INEL"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstRct, "+RCT_pass"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstTvx, "+kTVX(truth)"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstNoTimeFrameBorder, "+kNoTFB(truth)"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstNoItsRofBorder, "+kNoITSROFB(truth)"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstNoSameBunchPileupTruth, "+kNoSBP(truth)"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstHasCollision, "+hasColl"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstVertexZ, "+|zReco|<10"); - hMcCollisionsBcBitsFirst->GetXaxis()->SetBinLabel(BcBitsFirstNoSplit, "+noSplit"); + {HistType::kTH1F, {axisSelectionBcBitsFirst}}); + setAxisLabels(registry.get(HIST("h_mccollisions_eventselection_bcBitsFirst"))->GetXaxis(), bcBitsFirstLabels); + + if (isSel8FullPbPb) { + AxisSpec axisPbPbTruthSelectionOnly = {PbPbTruthSelectionOnlyNBins, 0.5, static_cast(PbPbTruthSelectionOnlyNBins) + 0.5, "event selection (PbPb truth-only)"}; + std::vector timeRangeOnlyLabels = {"INEL", "+RCT_pass", "+kTVX(truth)", "+kNoTFB(truth)", "+kNoITSROFB(truth)", "+kNoCollInTimeRangeStandard(truth)", "+hasColl", "+|zReco|<10", "+noSplit"}; + std::vector rofOnlyLabels = {"INEL", "+RCT_pass", "+kTVX(truth)", "+kNoTFB(truth)", "+kNoITSROFB(truth)", "+kNoCollInRofStandard(truth)", "+hasColl", "+|zReco|<10", "+noSplit"}; + + registry.add("h2_jet_pt_part_eventselection_bcBitsFirst_timeRangeOnly_truth", + "part jet pT vs TimeRange-only truth selection;#it{p}_{T,jet}^{part} (GeV/#it{c});event selection;weighted counts", + {HistType::kTH2F, {jetPtAxis, axisPbPbTruthSelectionOnly}}); + registry.add("h_mccollisions_eventselection_bcBitsFirst_timeRangeOnly_truth", + "weighted MC events vs TimeRange-only truth selection;event selection;weighted events", + {HistType::kTH1F, {axisPbPbTruthSelectionOnly}}); + setAxisLabels(registry.get(HIST("h2_jet_pt_part_eventselection_bcBitsFirst_timeRangeOnly_truth"))->GetYaxis(), timeRangeOnlyLabels); + setAxisLabels(registry.get(HIST("h_mccollisions_eventselection_bcBitsFirst_timeRangeOnly_truth"))->GetXaxis(), timeRangeOnlyLabels); + + registry.add("h2_jet_pt_part_eventselection_bcBitsFirst_rofOnly_truth", + "part jet pT vs ROF-only truth selection;#it{p}_{T,jet}^{part} (GeV/#it{c});event selection;weighted counts", + {HistType::kTH2F, {jetPtAxis, axisPbPbTruthSelectionOnly}}); + registry.add("h_mccollisions_eventselection_bcBitsFirst_rofOnly_truth", + "weighted MC events vs ROF-only truth selection;event selection;weighted events", + {HistType::kTH1F, {axisPbPbTruthSelectionOnly}}); + setAxisLabels(registry.get(HIST("h2_jet_pt_part_eventselection_bcBitsFirst_rofOnly_truth"))->GetYaxis(), rofOnlyLabels); + setAxisLabels(registry.get(HIST("h_mccollisions_eventselection_bcBitsFirst_rofOnly_truth"))->GetXaxis(), rofOnlyLabels); + } + } + + if (doprocessTruthMultFT0CvsRecoAmplitude) { + AxisSpec truthMultFT0CAxis = {500, -0.5, 499.5, "truth multFT0C"}; + AxisSpec recoFT0CAmplitudeAxis = {1000, 0., 20000., "A_{FT0C}^{reco}"}; + registry.add("h2_mccollision_mult_ft0c_found_ft0_sum_amp_c", + "truth multFT0C vs reconstructed A_{FT0C};truth multFT0C;A_{FT0C}^{reco};weighted counts", + {HistType::kTH2F, {truthMultFT0CAxis, recoFT0CAmplitudeAxis}}); + } + + if (doprocessTruthMultFT0CvsRecoITS567) { + AxisSpec truthMultFT0CAxis = {500, -0.5, 499.5, "truth multFT0C"}; + AxisSpec recoNITS567Axis = {1000, -0.5, 999.5, "N_{ITS567}^{reco}"}; + AxisSpec truthMultFT0CPositiveAxis = {2, -0.5, 1.5, "truth multFT0C > 0"}; + AxisSpec recoNITS567PositiveAxis = {2, -0.5, 1.5, "N_{ITS567}^{reco} > 0"}; + registry.add("h2_mccollision_mult_ft0c_reco_n_its567", + "truth multFT0C vs reconstructed N_{ITS567};truth multFT0C;N_{ITS567}^{reco};weighted counts", + {HistType::kTH2F, {truthMultFT0CAxis, recoNITS567Axis}}); + registry.add("h2_mccollision_mult_ft0c_positive_vs_reco_its567_positive", + "truth multFT0C activity vs reconstructed ITS567 activity;truth multFT0C > 0;N_{ITS567}^{reco} > 0;weighted counts", + {HistType::kTH2F, {truthMultFT0CPositiveAxis, recoNITS567PositiveAxis}}); + auto binaryActivity = registry.get(HIST("h2_mccollision_mult_ft0c_positive_vs_reco_its567_positive")); + binaryActivity->GetXaxis()->SetBinLabel(1, "multFT0C = 0"); + binaryActivity->GetXaxis()->SetBinLabel(2, "multFT0C > 0"); + binaryActivity->GetYaxis()->SetBinLabel(1, "N_{ITS567} = 0"); + binaryActivity->GetYaxis()->SetBinLabel(2, "N_{ITS567} > 0"); } } @@ -273,6 +368,126 @@ struct JetCrossSectionEfficiency { : simPtRef / std::pow(weight, 1.0f / pTHatExponent); } + template + bool configureTruthRofParameters(TBC const& truthBC) + { + if (truthRofOffsetInBC >= 0 && truthRofLengthInBC > 0) { + cachedTruthRofOffsetInBC = truthRofOffsetInBC; + cachedTruthRofLengthInBC = truthRofLengthInBC; + return true; + } + + if (cachedTruthRofRun == truthBC.runNumber() && cachedTruthRofLengthInBC > 0) { + return true; + } + + auto alppar = ccdb->getForTimeStamp>("ITS/Config/AlpideParam", truthBC.timestamp()); + if (alppar == nullptr) { + LOGF(fatal, "Could not retrieve ITS/Config/AlpideParam for sel8FullPbPb truth selections (run %d, timestamp %" PRIu64 ")", truthBC.runNumber(), static_cast(truthBC.timestamp())); + return false; + } + + cachedTruthRofRun = truthBC.runNumber(); + cachedTruthRofOffsetInBC = truthRofOffsetInBC >= 0 ? truthRofOffsetInBC : alppar->roFrameBiasInBC; + cachedTruthRofLengthInBC = truthRofLengthInBC > 0 ? truthRofLengthInBC : alppar->roFrameLengthInBC; + if (cachedTruthRofLengthInBC <= 0) { + LOGF(fatal, "Invalid ITS ROF length %" PRId64 " BC for sel8FullPbPb truth selections", static_cast(cachedTruthRofLengthInBC)); + return false; + } + LOGF(info, "sel8FullPbPb truth selections use ITS ROF offset %" PRId64 " and length %" PRId64 " BC for run %d", static_cast(cachedTruthRofOffsetInBC), static_cast(cachedTruthRofLengthInBC), cachedTruthRofRun); + return true; + } + + int64_t truthRofId(uint64_t globalBC) const + { + // Match EventSelectionModule: use the ITS ROF bias and length from DPLAlpideParam, with one orbit added to avoid a negative numerator. + return (static_cast(globalBC) + o2::constants::lhc::LHCMaxBunches - cachedTruthRofOffsetInBC) / cachedTruthRofLengthInBC; + } + + template + void requireValidTruthFT0CActivity(TMcCollision const& mccollision) const + { + if (mccollision.multFT0C() < 0.0f) { + LOGF(fatal, "sel8FullPbPb truth selection requires valid MC multFT0C, but multFT0C < 0 was found. Enable MC multiplicity information in the derived-data production."); + } + } + + template + float truthFT0CActivity(TMcCollision const& mccollision) const + { + requireValidTruthFT0CActivity(mccollision); + return mccollision.multFT0C(); + } + + template + TruthPbPbSelections evaluateTruthPbPbSelections(TMcCollision const& mccollision, TAllMcCollisions const& allMcCollisions) + { + TruthPbPbSelections selections; + requireValidTruthFT0CActivity(mccollision); + auto truthBC = mccollision.template bc_as(); + const uint64_t currentGlobalBC = truthBC.globalBC(); + if (currentGlobalBC == std::numeric_limits::max() || !configureTruthRofParameters(truthBC)) { + return selections; + } + + const int64_t currentRof = truthRofId(currentGlobalBC); + bool hasNarrowActivity = false; + bool hasHighActivityInTimeRange = false; + bool hasHighActivityInSameRof = false; + bool hasCloseVzActivityInSameRof = false; + + // EventSelectionModule uses reconstructed foundGlobalBC, FT0C digit amplitudes and ITS layer 5-7 tracks. + // Those are unavailable for hasColl=false MC collisions. + // This truth-only approximation uses nominal MC BCs and multFT0C; its thresholds are intentionally configurable. + // The scan is limited to the current input chunk because JMcCollisions has no truth-side found-BC or time-frame association. + for (auto const& otherMcCollision : allMcCollisions) { + if (otherMcCollision.globalIndex() == mccollision.globalIndex()) { + continue; + } + + auto otherTruthBC = otherMcCollision.template bc_as(); + const uint64_t otherGlobalBC = otherTruthBC.globalBC(); + if (otherTruthBC.runNumber() != truthBC.runNumber() || otherGlobalBC == std::numeric_limits::max()) { + continue; + } + + const int64_t deltaGlobalBC = static_cast(otherGlobalBC) - static_cast(currentGlobalBC); + const float deltaTimeUs = static_cast(deltaGlobalBC) * o2::constants::lhc::LHCBunchSpacingNS / 1000.0f; + const bool inNarrowWindow = std::abs(deltaTimeUs) < truthTimeRangeNarrowUs; + const bool inStandardTimeWindow = deltaTimeUs > truthTimeRangeStandardMinUs && deltaTimeUs < truthTimeRangeStandardMaxUs; + const bool isSameRof = truthRofId(otherGlobalBC) == currentRof; + if (!inNarrowWindow && !inStandardTimeWindow && !isSameRof) { + continue; + } + + const float otherActivity = truthFT0CActivity(otherMcCollision); + const bool hasTrackPresenceProxy = otherActivity > truthFT0CActivityMinForTrackProxy; + + if (inNarrowWindow && hasTrackPresenceProxy) { + hasNarrowActivity = true; + } + if (inStandardTimeWindow && otherActivity > truthFT0CActivityThresholdTimeRange) { + hasHighActivityInTimeRange = true; + } + + if (isSameRof) { + if (otherActivity > truthFT0CActivityThresholdROF) { + hasHighActivityInSameRof = true; + } + // The reco bit uses ITS layer 5-7 tracks and reconstructed PV z. + // The truth approximation uses multFT0C and MC collision z instead. + if (hasTrackPresenceProxy && std::abs(otherMcCollision.posZ() - mccollision.posZ()) < truthRofCloseVzMax) { + hasCloseVzActivityInSameRof = true; + } + } + } + + selections.valid = true; + selections.noCollInTimeRangeStandard = !hasNarrowActivity && !hasHighActivityInTimeRange; + selections.noCollInRofStandard = !hasHighActivityInSameRof && !hasCloseVzActivityInSameRof; + return selections; + } + template bool isAcceptedJet(TJets const& jet) { @@ -304,6 +519,66 @@ struct JetCrossSectionEfficiency { return true; } + void processTruthMultFT0CvsRecoAmplitude(aod::JetMcCollisions::iterator const& mccollision, + soa::SmallGroups const& collisions, + CollisionsWithEvSels const&, + aod::FT0s const& ft0s) + { + if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + return; + } + if (mccollision.multFT0C() < 0.0f) { + return; + } + // A single associated JCollision gives a one-to-one MC-collision to reconstructed-collision association. + // The parent index restores the original EvSels relation used by EventSelectionModule. + if (collisions.size() != 1) { + return; + } + auto const& collision = collisions.begin(); + auto originalCollision = collision.template collision_as(); + if (!originalCollision.has_foundFT0()) { + return; + } + auto foundFT0 = ft0s.rawIteratorAt(originalCollision.foundFT0Id()); + registry.fill(HIST("h2_mccollision_mult_ft0c_found_ft0_sum_amp_c"), mccollision.multFT0C(), foundFT0.sumAmpC(), mccollision.weight()); + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processTruthMultFT0CvsRecoAmplitude, + "truth multFT0C vs reconstructed foundFT0 A_FT0C for one-to-one MC/reco collision associations", false); + + void processTruthMultFT0CvsRecoITS567(aod::JetMcCollisions::iterator const& mccollision, + soa::SmallGroups const& collisions, + CollisionsWithEvSels const&, + FullTracksIU const&) + { + if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + return; + } + if (mccollision.multFT0C() < 0.0f) { + return; + } + // This matches EventSelectionModule::vTracksITS567perColl: only PV-contributor tracks with at least five ITS clusters are counted. + if (collisions.size() != 1) { + return; + } + auto const& collision = collisions.begin(); + auto originalCollision = collision.template collision_as(); + auto const& collisionPvTracks = pvTracks.sliceBy(pvTracksPerCollision, originalCollision.globalIndex()); + int nITS567 = 0; + for (const auto& track : collisionPvTracks) { + if (track.itsNCls() >= MinITSClustersForOccupancy) { + ++nITS567; + } + } + + const float weight = mccollision.weight(); + registry.fill(HIST("h2_mccollision_mult_ft0c_reco_n_its567"), mccollision.multFT0C(), nITS567, weight); + registry.fill(HIST("h2_mccollision_mult_ft0c_positive_vs_reco_its567_positive"), + mccollision.multFT0C() > 0.0f, nITS567 > 0, weight); + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processTruthMultFT0CvsRecoITS567, + "truth multFT0C activity vs reconstructed PV-contributor ITS567 activity for one-to-one MC/reco collision associations", false); + void processCrossSectionEfficiency(aod::JetMcCollisions::iterator const& mccollision, soa::SmallGroups const& collisions, soa::Join const& jets, @@ -331,7 +606,12 @@ struct JetCrossSectionEfficiency { } bool noSplitPass = (acceptSplitCollisions == NonSplitOnly) ? (collisions.size() == 1) : true; - bool passesTVX = false, passesNoTFB = false, passesNoITSROFB = false, passesNoSBP = false; + bool passesTVX = false; + bool passesNoTFB = false; + bool passesNoITSROFB = false; + bool passesNoSBP = false; + bool passesNoCollInTimeRangeStandard = false; + bool passesNoCollInRofStandard = false; if (hasRecoColl) { auto const& col = collisions.begin(); auto evSel = col.eventSel(); @@ -339,24 +619,33 @@ struct JetCrossSectionEfficiency { passesNoTFB = (evSel & (1u << jetderiveddatautilities::JCollisionSel::selNoTimeFrameBorder)) != 0u; passesNoITSROFB = (evSel & (1u << jetderiveddatautilities::JCollisionSel::selNoITSROFrameBorder)) != 0u; passesNoSBP = (evSel & (1u << jetderiveddatautilities::JCollisionSel::selNoSameBunchPileup)) != 0u; + passesNoCollInTimeRangeStandard = (evSel & (1u << jetderiveddatautilities::JCollisionSel::selNoCollInTimeRangeStandard)) != 0u; + passesNoCollInRofStandard = (evSel & (1u << jetderiveddatautilities::JCollisionSel::selNoCollInRofStandard)) != 0u; } bool passesRct = applyRCT ? (mccollision.bc_as().rct_raw() & rctMask) == 0 : true; - bool pass[CollRecoFirstNBins + 1] = {false, true, passesRct, hasRecoColl, passesZvtxCutReco, - hasRecoColl && noSplitPass, passesTVX, - applyTFB ? passesNoTFB : true, - applyROFB ? passesNoITSROFB : true, - applySBP ? passesNoSBP : true}; + std::vector pass = {true, passesRct, hasRecoColl, passesZvtxCutReco, + hasRecoColl && noSplitPass, passesTVX, + applyTFB ? passesNoTFB : true, + applyROFB ? passesNoITSROFB : true}; + if (applyNoCollInTimeRangeStandard) { + pass.push_back(passesNoCollInTimeRangeStandard); + } + if (applyNoCollInRofStandard) { + pass.push_back(passesNoCollInRofStandard); + } + if (!isSel8FullPbPb) { + pass.push_back(applySBP ? passesNoSBP : true); + } - // Unified weight handling: MB MC -> weight=1 (no-op), JJ MC -> per-event sigma fraction. float weight = mccollision.weight(); - int sMax = 0; - for (int s = CollRecoFirstInel; s <= CollRecoFirstNBins; ++s) { - if (!pass[s]) + for (size_t s = 0; s < pass.size(); ++s) { + if (!pass[s]) { break; - registry.fill(HIST("h_mccollisions_eventselection_collRecoFirst"), static_cast(s), weight); - sMax = s; + } + registry.fill(HIST("h_mccollisions_eventselection_collRecoFirst"), static_cast(s + 1), weight); + sMax = static_cast(s + 1); } if (sMax == 0) { return; @@ -373,7 +662,7 @@ struct JetCrossSectionEfficiency { !isAcceptedJet(jet)) { continue; } - for (int s = CollRecoFirstInel; s <= sMax; ++s) { + for (int s = 1; s <= sMax; ++s) { registry.fill(HIST("h2_jet_pt_part_eventselection_collRecoFirst"), jet.pt(), static_cast(s), weight); } } @@ -397,8 +686,11 @@ struct JetCrossSectionEfficiency { bool passesNoTFBTruth = truthBC.selection_bit(aod::evsel::kNoTimeFrameBorder); bool passesNoITSROFBTruth = truthBC.selection_bit(aod::evsel::kNoITSROFrameBorder); - auto sameBC = allMcCollisions.sliceBy(mcCollsPerBC, mccollision.bcId()); - bool truthNoSBP = (sameBC.size() == 1); + bool truthNoSBP = true; + if (applySBP) { + auto sameBC = allMcCollisions.sliceBy(mcCollsPerBC, mccollision.bcId()); + truthNoSBP = (sameBC.size() == 1); + } bool hasRecoColl = (collisions.size() >= 1); bool passesZvtxCutReco = false; @@ -417,21 +709,64 @@ struct JetCrossSectionEfficiency { } bool noSplitPass = (acceptSplitCollisions == NonSplitOnly) ? (collisions.size() == 1) : true; + TruthPbPbSelections truthPbPbSelections; + if (isSel8FullPbPb) { + truthPbPbSelections = evaluateTruthPbPbSelections(mccollision, allMcCollisions); + } + bool passesRct = applyRCT ? (truthBC.rct_raw() & rctMask) == 0 : true; - bool pass[BcBitsFirstNBins + 1] = {false, true, passesRct, passesTVXTruth, - applyTFB ? passesNoTFBTruth : true, - applyROFB ? passesNoITSROFBTruth : true, - applySBP ? truthNoSBP : true, - hasRecoColl, passesZvtxCutReco, hasRecoColl && noSplitPass}; + std::vector pass = {true, passesRct, passesTVXTruth, + applyTFB ? passesNoTFBTruth : true, + applyROFB ? passesNoITSROFBTruth : true}; + if (applyNoCollInTimeRangeStandard) { + pass.push_back(truthPbPbSelections.valid && truthPbPbSelections.noCollInTimeRangeStandard); + } + if (applyNoCollInRofStandard) { + pass.push_back(truthPbPbSelections.valid && truthPbPbSelections.noCollInRofStandard); + } + if (!isSel8FullPbPb) { + pass.push_back(applySBP ? truthNoSBP : true); + } + pass.push_back(hasRecoColl); + pass.push_back(passesZvtxCutReco); + pass.push_back(hasRecoColl && noSplitPass); float weight = mccollision.weight(); + int sMaxTimeRangeOnly = 0; + int sMaxRofOnly = 0; + if (isSel8FullPbPb) { + std::vector timeRangeOnlyPass = {true, passesRct, passesTVXTruth, + passesNoTFBTruth, passesNoITSROFBTruth, + truthPbPbSelections.valid && truthPbPbSelections.noCollInTimeRangeStandard, + hasRecoColl, passesZvtxCutReco, hasRecoColl && noSplitPass}; + std::vector rofOnlyPass = {true, passesRct, passesTVXTruth, + passesNoTFBTruth, passesNoITSROFBTruth, + truthPbPbSelections.valid && truthPbPbSelections.noCollInRofStandard, + hasRecoColl, passesZvtxCutReco, hasRecoColl && noSplitPass}; + for (size_t s = 0; s < timeRangeOnlyPass.size(); ++s) { + if (!timeRangeOnlyPass[s]) { + break; + } + registry.fill(HIST("h_mccollisions_eventselection_bcBitsFirst_timeRangeOnly_truth"), static_cast(s + 1), weight); + sMaxTimeRangeOnly = static_cast(s + 1); + } + for (size_t s = 0; s < rofOnlyPass.size(); ++s) { + if (!rofOnlyPass[s]) { + break; + } + registry.fill(HIST("h_mccollisions_eventselection_bcBitsFirst_rofOnly_truth"), static_cast(s + 1), weight); + sMaxRofOnly = static_cast(s + 1); + } + } + int sMax = 0; - for (int s = BcBitsFirstInel; s <= BcBitsFirstNBins; ++s) { - if (!pass[s]) + for (size_t s = 0; s < pass.size(); ++s) { + if (!pass[s]) { break; - registry.fill(HIST("h_mccollisions_eventselection_bcBitsFirst"), static_cast(s), weight); - sMax = s; + } + registry.fill(HIST("h_mccollisions_eventselection_bcBitsFirst"), static_cast(s + 1), weight); + sMax = static_cast(s + 1); } if (sMax == 0) { return; @@ -448,9 +783,17 @@ struct JetCrossSectionEfficiency { !isAcceptedJet(jet)) { continue; } - for (int s = BcBitsFirstInel; s <= sMax; ++s) { + for (int s = 1; s <= sMax; ++s) { registry.fill(HIST("h2_jet_pt_part_eventselection_bcBitsFirst"), jet.pt(), static_cast(s), weight); } + if (isSel8FullPbPb) { + for (int s = 1; s <= sMaxTimeRangeOnly; ++s) { + registry.fill(HIST("h2_jet_pt_part_eventselection_bcBitsFirst_timeRangeOnly_truth"), jet.pt(), static_cast(s), weight); + } + for (int s = 1; s <= sMaxRofOnly; ++s) { + registry.fill(HIST("h2_jet_pt_part_eventselection_bcBitsFirst_rofOnly_truth"), jet.pt(), static_cast(s), weight); + } + } } } PROCESS_SWITCH(JetCrossSectionEfficiency, processCrossSectionEfficiencyBcBitsFirst,