From 0c135b432b76d6f2b76ae2278d4e110bc4787387 Mon Sep 17 00:00:00 2001 From: Wooseok Ham Date: Tue, 29 Sep 2026 13:31:05 +0900 Subject: [PATCH 1/2] PWGJE: Add hasColl and collision-environment QA to jet cross-section efficiency --- PWGJE/Tasks/jetCrossSectionEfficiency.cxx | 734 +++++++++++++++++++++- 1 file changed, 708 insertions(+), 26 deletions(-) diff --git a/PWGJE/Tasks/jetCrossSectionEfficiency.cxx b/PWGJE/Tasks/jetCrossSectionEfficiency.cxx index 9bcc1195454..874abdea4fb 100644 --- a/PWGJE/Tasks/jetCrossSectionEfficiency.cxx +++ b/PWGJE/Tasks/jetCrossSectionEfficiency.cxx @@ -43,13 +43,19 @@ #include #include #include +#include +#include #include #include #include #include #include +#include #include +#include +#include +#include #include using namespace o2; @@ -101,6 +107,7 @@ struct JetCrossSectionEfficiency { 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"}; + Configurable tvxEnvironmentRequireOtherBcBaseline{"tvxEnvironmentRequireOtherBcBaseline", false, "require RCT, TVX, NoTFB and NoITSROFB for other BCs in the TVX-BC environment QA"}; o2::aod::rctsel::RCTFlagsChecker rctChecker; uint64_t rctMask = 0; @@ -143,6 +150,15 @@ struct JetCrossSectionEfficiency { static constexpr float ConfigSwitchHigh = 9998.0f; static constexpr float BrokenPtHardSentinel = 1.0f; static constexpr int MinITSClustersForOccupancy = 5; + static constexpr int HasCollQaMaxCount = 500; + static constexpr int HasCollQaMaxDeltaGlobalBC = 100000; + static constexpr int64_t NoOtherTruthCollision = -1; + + static std::vector getNearestAbsDeltaGlobalBCEdges() + { + // Keep the no-other-BC sentinel separate, then resolve the near-BC region most relevant for reconstruction losses. + return {-1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 7.5, 15.5, 31.5, 63.5, 127.5, 255.5, 511.5, 1000.5, 2000.5, 4000.5, 8000.5, 16000.5, 32000.5, 64000.5, static_cast(HasCollQaMaxDeltaGlobalBC) + 0.5}; + } struct TruthPbPbSelections { bool valid = false; @@ -150,12 +166,44 @@ struct JetCrossSectionEfficiency { bool noCollInRofStandard = false; }; + struct TruthHasCollEnvironment { + bool valid = false; + int nOtherSameTruthRof = 0; + int64_t nearestAbsDeltaGlobalBC = NoOtherTruthCollision; + }; + + struct VisibleBCEnvironment { + bool valid = false; + int nOtherTVXSameRof = 0; + int64_t nearestAbsDeltaTVXGlobalBC = NoOtherTruthCollision; + int64_t originalBCIndex = -1; + }; + + struct VisibleBCEnvironmentCache { + std::map byBCIndex; + std::map, VisibleBCEnvironment> byRunAndGlobalBC; + }; + using JetCollisionsMCDWithParent = soa::Join; + using JMcCollisionsWithParent = soa::Join; + using JetCollisionsWithParent = soa::Join; using CollisionsWithEvSels = soa::Join; + using OriginalBCsWithSels = soa::Join; + using OriginalCollisionsWithMcLabels = soa::Join; using FullTracksIU = soa::Join; + using FullTracksIUWithMcLabels = soa::Join; + using ChargedMCDJetsWithConstituents = soa::Join; + using ChargedMCPJetsWithConstituents = soa::Join; Partition pvTracks = ((aod::track::flags & static_cast(o2::aod::track::PVContributor)) == static_cast(o2::aod::track::PVContributor)); Preslice pvTracksPerCollision = aod::track::collisionId; + PresliceUnsorted collisionsPerFoundBC = aod::evsel::foundBCId; + PresliceUnsorted jetCollisionsPerMcCollision = aod::jmccollisionlb::mcCollisionId; + Preslice mcpJetsPerMcCollision = aod::jet::mcCollisionId; + Preslice mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId; + PresliceUnsorted recoTracksPerMcParticle = aod::mctracklabel::mcParticleId; + PresliceUnsorted jetCollisionsPerOriginalCollision = aod::jcollision::collisionId; + Preslice detectorJetsPerJetCollision = aod::jet::collisionId; enum EventSelectionPreset { PresetSelTvx = 0, @@ -242,7 +290,24 @@ struct JetCrossSectionEfficiency { break; } - if (doprocessCrossSectionEfficiencyBcBitsFirst && isSel8FullPbPb && + const bool isDataTVXBCEnvironmentQAEnabled = doprocessDataTVXBCEnvironmentQA; + const bool isMCQAEnabled = doprocessMCTruthMultFT0CvsRecoAmplitude || + doprocessMCTruthMultFT0CvsRecoITS567 || + doprocessMCEnvironmentQA || + doprocessMCHasCollTrackFateQA || + doprocessMCHasCollJetSurvivalQA || + doprocessMCHasCollJetEnvironmentQA || + doprocessCrossSectionEfficiency || + doprocessCrossSectionEfficiencyBcBitsFirst; + if (isDataTVXBCEnvironmentQAEnabled && isMCQAEnabled) { + LOGF(fatal, "processDataTVXBCEnvironmentQA cannot run with MC-only QA or cross-section-efficiency processes; enable the data process alone for data AO2Ds"); + } + + const bool needsRofParameters = (doprocessCrossSectionEfficiencyBcBitsFirst && isSel8FullPbPb) || + doprocessMCEnvironmentQA || + doprocessMCHasCollJetEnvironmentQA || + doprocessDataTVXBCEnvironmentQA; + if (needsRofParameters && (truthRofOffsetInBC < 0 || truthRofLengthInBC <= 0)) { ccdb->setURL(static_cast(truthRofCcdbUrl)); ccdb->setCaching(true); @@ -327,7 +392,7 @@ struct JetCrossSectionEfficiency { } } - if (doprocessTruthMultFT0CvsRecoAmplitude) { + if (doprocessMCTruthMultFT0CvsRecoAmplitude) { 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", @@ -335,7 +400,7 @@ struct JetCrossSectionEfficiency { {HistType::kTH2F, {truthMultFT0CAxis, recoFT0CAmplitudeAxis}}); } - if (doprocessTruthMultFT0CvsRecoITS567) { + if (doprocessMCTruthMultFT0CvsRecoITS567) { 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"}; @@ -352,6 +417,170 @@ struct JetCrossSectionEfficiency { binaryActivity->GetYaxis()->SetBinLabel(1, "N_{ITS567} = 0"); binaryActivity->GetYaxis()->SetBinLabel(2, "N_{ITS567} > 0"); } + + if (doprocessMCEnvironmentQA) { + AxisSpec nOtherSameTruthRofAxis = {HasCollQaMaxCount + 1, -0.5, static_cast(HasCollQaMaxCount) + 0.5, "N_{other same truth ROF}"}; + AxisSpec nearestAbsDeltaGlobalBCAxis = {getNearestAbsDeltaGlobalBCEdges(), "nearest |#Delta global BC|"}; + AxisSpec targetMultFT0CAxis = {500, -0.5, 499.5, "target truth multFT0C"}; + AxisSpec hasCollAxis = {2, -0.5, 1.5, "JetDerived hasColl"}; + + registry.add("h2_n_other_same_truth_rof_vs_has_coll", + "other MC collisions in the same truth ITS ROF vs hasColl;N_{other same truth ROF};JetDerived hasColl;weighted counts", + {HistType::kTH2F, {nOtherSameTruthRofAxis, hasCollAxis}}); + registry.add("h2_nearest_abs_delta_global_bc_vs_has_coll", + "nearest other MC collision |#Delta global BC| vs hasColl;nearest |#Delta global BC|;JetDerived hasColl;weighted counts", + {HistType::kTH2F, {nearestAbsDeltaGlobalBCAxis, hasCollAxis}}); + registry.add("h2_target_mult_ft0c_vs_has_coll", + "target truth multFT0C vs hasColl;target truth multFT0C;JetDerived hasColl;weighted counts", + {HistType::kTH2F, {targetMultFT0CAxis, hasCollAxis}}); + registry.add("h2_target_mult_ft0c_vs_n_other_same_truth_rof_all", + "target truth multFT0C vs other MC collisions in the same truth ITS ROF;target truth multFT0C;N_{other same truth ROF};weighted counts", + {HistType::kTH2F, {targetMultFT0CAxis, nOtherSameTruthRofAxis}}); + registry.add("h2_target_mult_ft0c_vs_n_other_same_truth_rof_has_coll", + "target truth multFT0C vs other MC collisions in the same truth ITS ROF for hasColl;target truth multFT0C;N_{other same truth ROF};weighted counts", + {HistType::kTH2F, {targetMultFT0CAxis, nOtherSameTruthRofAxis}}); + registry.add("h2_target_mult_ft0c_vs_nearest_abs_delta_global_bc_all", + "target truth multFT0C vs nearest other MC collision distance;target truth multFT0C;nearest |#Delta global BC|;weighted counts", + {HistType::kTH2F, {targetMultFT0CAxis, nearestAbsDeltaGlobalBCAxis}}); + registry.add("h2_target_mult_ft0c_vs_nearest_abs_delta_global_bc_has_coll", + "target truth multFT0C vs nearest other MC collision distance for hasColl;target truth multFT0C;nearest |#Delta global BC|;weighted counts", + {HistType::kTH2F, {targetMultFT0CAxis, nearestAbsDeltaGlobalBCAxis}}); + + auto nOtherSameTruthRof = registry.get(HIST("h2_n_other_same_truth_rof_vs_has_coll")); + nOtherSameTruthRof->GetYaxis()->SetBinLabel(1, "noColl"); + nOtherSameTruthRof->GetYaxis()->SetBinLabel(2, "hasColl"); + auto nearestAbsDeltaGlobalBC = registry.get(HIST("h2_nearest_abs_delta_global_bc_vs_has_coll")); + nearestAbsDeltaGlobalBC->GetXaxis()->SetBinLabel(1, "no other valid MC collision"); + nearestAbsDeltaGlobalBC->GetYaxis()->SetBinLabel(1, "noColl"); + nearestAbsDeltaGlobalBC->GetYaxis()->SetBinLabel(2, "hasColl"); + auto targetMultFT0C = registry.get(HIST("h2_target_mult_ft0c_vs_has_coll")); + targetMultFT0C->GetYaxis()->SetBinLabel(1, "noColl"); + targetMultFT0C->GetYaxis()->SetBinLabel(2, "hasColl"); + registry.get(HIST("h2_target_mult_ft0c_vs_nearest_abs_delta_global_bc_all"))->GetYaxis()->SetBinLabel(1, "no other valid MC collision"); + registry.get(HIST("h2_target_mult_ft0c_vs_nearest_abs_delta_global_bc_has_coll"))->GetYaxis()->SetBinLabel(1, "no other valid MC collision"); + } + + if (doprocessDataTVXBCEnvironmentQA || doprocessMCEnvironmentQA) { + AxisSpec nOtherTVXSameRofAxis = {HasCollQaMaxCount + 1, -0.5, static_cast(HasCollQaMaxCount) + 0.5, "N_{other TVX BCs in same ITS ROF}"}; + AxisSpec nearestAbsDeltaTVXGlobalBCEnvironmentAxis = {getNearestAbsDeltaGlobalBCEdges(), "nearest other TVX |#Delta global BC|"}; + AxisSpec hasRecoCollisionAtTargetBCAxis = {2, -0.5, 1.5, "reconstructed collision at target TVX BC"}; + AxisSpec targetFT0CAmplitudeAxis = {1000, 0., 20000., "A_{FT0C} at target TVX BC"}; + + registry.add("h2_n_other_tvx_same_rof_vs_nearest_abs_delta_global_bc_tvx", + "TVX-BC environment;N_{other TVX BCs in same ITS ROF};nearest other TVX |#Delta global BC|;BCs", + {HistType::kTH2F, {nOtherTVXSameRofAxis, nearestAbsDeltaTVXGlobalBCEnvironmentAxis}}); + registry.add("h2_n_other_tvx_same_rof_vs_nearest_abs_delta_global_bc_tvx_has_reco", + "TVX-BC environment for target BCs with a reconstructed collision;N_{other TVX BCs in same ITS ROF};nearest other TVX |#Delta global BC|;BCs", + {HistType::kTH2F, {nOtherTVXSameRofAxis, nearestAbsDeltaTVXGlobalBCEnvironmentAxis}}); + registry.add("h2_n_other_tvx_same_rof_vs_has_reco_collision", + "TVX-BC environment vs reconstructed-collision proxy;N_{other TVX BCs in same ITS ROF};reconstructed collision at target TVX BC;BCs", + {HistType::kTH2F, {nOtherTVXSameRofAxis, hasRecoCollisionAtTargetBCAxis}}); + registry.add("h2_nearest_abs_delta_global_bc_tvx_vs_has_reco_collision", + "TVX-BC environment vs reconstructed-collision proxy;nearest other TVX |#Delta global BC|;reconstructed collision at target TVX BC;BCs", + {HistType::kTH2F, {nearestAbsDeltaTVXGlobalBCEnvironmentAxis, hasRecoCollisionAtTargetBCAxis}}); + registry.add("h2_target_ft0c_amplitude_vs_n_other_tvx_same_rof", + "FT0C amplitude vs TVX-BC environment;A_{FT0C} at target TVX BC;N_{other TVX BCs in same ITS ROF};BCs", + {HistType::kTH2F, {targetFT0CAmplitudeAxis, nOtherTVXSameRofAxis}}); + + auto nOtherTVXSameRof = registry.get(HIST("h2_n_other_tvx_same_rof_vs_has_reco_collision")); + nOtherTVXSameRof->GetYaxis()->SetBinLabel(1, "no reco collision"); + nOtherTVXSameRof->GetYaxis()->SetBinLabel(2, "reco collision exists"); + auto nearestAbsDeltaTVXGlobalBC = registry.get(HIST("h2_nearest_abs_delta_global_bc_tvx_vs_has_reco_collision")); + nearestAbsDeltaTVXGlobalBC->GetXaxis()->SetBinLabel(1, "no other TVX BC in run"); + nearestAbsDeltaTVXGlobalBC->GetYaxis()->SetBinLabel(1, "no reco collision"); + nearestAbsDeltaTVXGlobalBC->GetYaxis()->SetBinLabel(2, "reco collision exists"); + registry.get(HIST("h2_n_other_tvx_same_rof_vs_nearest_abs_delta_global_bc_tvx"))->GetYaxis()->SetBinLabel(1, "no other TVX BC in run"); + registry.get(HIST("h2_n_other_tvx_same_rof_vs_nearest_abs_delta_global_bc_tvx_has_reco"))->GetYaxis()->SetBinLabel(1, "no other TVX BC in run"); + } + + if (doprocessMCEnvironmentQA) { + AxisSpec nOtherSameTruthRofAxis = {HasCollQaMaxCount + 1, -0.5, static_cast(HasCollQaMaxCount) + 0.5, "N_{other MC collisions in same truth ITS ROF}"}; + AxisSpec nOtherTVXSameRofAxis = {HasCollQaMaxCount + 1, -0.5, static_cast(HasCollQaMaxCount) + 0.5, "N_{other TVX BCs in same ITS ROF}"}; + AxisSpec nearestAbsDeltaGlobalBCEnvironmentAxis = {getNearestAbsDeltaGlobalBCEdges(), "nearest |#Delta global BC|"}; + AxisSpec hasRecoCollisionAtTargetBCAxis = {2, -0.5, 1.5, "reconstructed collision at target TVX BC"}; + AxisSpec hasCollAxis = {2, -0.5, 1.5, "JetDerived hasColl"}; + + registry.add("h2_n_other_truth_same_rof_vs_n_other_tvx_same_rof", + "truth and TVX-BC same-ROF environment;N_{other MC collisions in same truth ITS ROF};N_{other TVX BCs in same ITS ROF};weighted MC collisions", + {HistType::kTH2F, {nOtherSameTruthRofAxis, nOtherTVXSameRofAxis}}); + registry.add("h2_nearest_abs_delta_global_bc_truth_vs_nearest_abs_delta_global_bc_tvx", + "truth and TVX-BC nearest collision environment;truth nearest |#Delta global BC|;TVX nearest |#Delta global BC|;weighted MC collisions", + {HistType::kTH2F, {nearestAbsDeltaGlobalBCEnvironmentAxis, nearestAbsDeltaGlobalBCEnvironmentAxis}}); + registry.add("h2_reco_collision_proxy_vs_has_coll", + "reconstructed-collision proxy vs JetDerived hasColl;reconstructed collision at target TVX BC;JetDerived hasColl;weighted MC collisions", + {HistType::kTH2F, {hasRecoCollisionAtTargetBCAxis, hasCollAxis}}); + auto recoCollisionProxy = registry.get(HIST("h2_reco_collision_proxy_vs_has_coll")); + recoCollisionProxy->GetXaxis()->SetBinLabel(1, "no reco collision"); + recoCollisionProxy->GetXaxis()->SetBinLabel(2, "reco collision exists"); + recoCollisionProxy->GetYaxis()->SetBinLabel(1, "noColl"); + recoCollisionProxy->GetYaxis()->SetBinLabel(2, "hasColl"); + } + + if (doprocessMCHasCollJetEnvironmentQA) { + AxisSpec nOtherTVXSameRofAxis = {HasCollQaMaxCount + 1, -0.5, static_cast(HasCollQaMaxCount) + 0.5, "N_{other TVX BCs in same ITS ROF}"}; + AxisSpec nearestAbsDeltaTVXGlobalBCEnvironmentAxis = {getNearestAbsDeltaGlobalBCEdges(), "nearest other TVX |#Delta global BC|"}; + AxisSpec hasCollAxis = {2, -0.5, 1.5, "JetDerived hasColl"}; + registry.add("hs_mcp_jet_pt_tvx_environment_has_coll", + "selected MCP jets vs TVX-BC environment and hasColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});N_{other TVX BCs in same ITS ROF};nearest other TVX |#Delta global BC|;JetDerived hasColl;weighted jets", + {HistType::kTHnSparseF, {jetPtAxis, nOtherTVXSameRofAxis, nearestAbsDeltaTVXGlobalBCEnvironmentAxis, hasCollAxis}}, true); + auto jetEnvironment = registry.get(HIST("hs_mcp_jet_pt_tvx_environment_has_coll")); + jetEnvironment->GetAxis(3)->SetBinLabel(1, "noColl"); + jetEnvironment->GetAxis(3)->SetBinLabel(2, "hasColl"); + jetEnvironment->GetAxis(2)->SetBinLabel(1, "no other TVX BC in run"); + } + + if (doprocessMCHasCollTrackFateQA) { + AxisSpec trackFateCountAxis = {HasCollQaMaxCount + 1, -0.5, static_cast(HasCollQaMaxCount) + 0.5, "N_{target-origin reconstructed tracks}"}; + AxisSpec hasCollAxis = {2, -0.5, 1.5, "JetDerived hasColl"}; + registry.add("h2_n_reco_tracks_from_target_vs_has_coll", + "reconstructed tracks from target MC collision vs hasColl;N_{target-origin reconstructed tracks};JetDerived hasColl;weighted counts", + {HistType::kTH2F, {trackFateCountAxis, hasCollAxis}}); + registry.add("h2_n_unassigned_tracks_from_target_vs_has_coll", + "unassigned tracks from target MC collision vs hasColl;N_{target-origin unassigned tracks};JetDerived hasColl;weighted counts", + {HistType::kTH2F, {trackFateCountAxis, hasCollAxis}}); + registry.add("h2_n_target_tracks_attached_to_other_labeled_collision_vs_has_coll", + "target-origin tracks attached to other-labeled collision vs hasColl;N_{tracks attached to other-labeled collision};JetDerived hasColl;weighted counts", + {HistType::kTH2F, {trackFateCountAxis, hasCollAxis}}); + + auto setHasCollAxisLabels = [](const auto& histogram) { + histogram->GetYaxis()->SetBinLabel(1, "noColl"); + histogram->GetYaxis()->SetBinLabel(2, "hasColl"); + }; + setHasCollAxisLabels(registry.get(HIST("h2_n_reco_tracks_from_target_vs_has_coll"))); + setHasCollAxisLabels(registry.get(HIST("h2_n_unassigned_tracks_from_target_vs_has_coll"))); + setHasCollAxisLabels(registry.get(HIST("h2_n_target_tracks_attached_to_other_labeled_collision_vs_has_coll"))); + } + + if (doprocessMCHasCollJetSurvivalQA) { + AxisSpec detectorJetExistsAxis = {2, -0.5, 1.5, "detector jet with target-origin reco track"}; + AxisSpec constituentPtFractionAxis = {300, 0., 3., "p_{T}^{reco target tracks} / #Sigma p_{T}^{truth constituents}"}; + + registry.add("h2_mcp_jet_pt_vs_has_detector_jet_hasColl", + "MCP jet pT vs detector-level jet survival for hasColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});detector jet with target-origin reco track;weighted counts", + {HistType::kTH2F, {jetPtAxis, detectorJetExistsAxis}}); + registry.add("h2_mcp_jet_pt_vs_has_detector_jet_noColl", + "MCP jet pT vs detector-level jet survival for noColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});detector jet with target-origin reco track;weighted counts", + {HistType::kTH2F, {jetPtAxis, detectorJetExistsAxis}}); + registry.add("h2_mcp_jet_pt_vs_reco_constituent_pt_fraction_hasColl", + "MCP jet pT vs target-origin reconstructed constituent pT fraction for hasColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});#Sigma p_{T}^{reco target tracks} / #Sigma p_{T}^{truth constituents};weighted counts", + {HistType::kTH2F, {jetPtAxis, constituentPtFractionAxis}}); + registry.add("h2_mcp_jet_pt_vs_reco_constituent_pt_fraction_noColl", + "MCP jet pT vs target-origin reconstructed constituent pT fraction for noColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});#Sigma p_{T}^{reco target tracks} / #Sigma p_{T}^{truth constituents};weighted counts", + {HistType::kTH2F, {jetPtAxis, constituentPtFractionAxis}}); + registry.add("h2_mcp_jet_pt_vs_best_detector_jet_capture_fraction_hasColl", + "MCP jet pT vs best detector-jet target-track capture fraction for hasColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});best #Sigma p_{T}^{reco target tracks} / #Sigma p_{T}^{truth constituents};weighted counts", + {HistType::kTH2F, {jetPtAxis, constituentPtFractionAxis}}); + registry.add("h2_mcp_jet_pt_vs_best_detector_jet_capture_fraction_noColl", + "MCP jet pT vs best detector-jet target-track capture fraction for noColl;#it{p}_{T,jet}^{MCP} (GeV/#it{c});best #Sigma p_{T}^{reco target tracks} / #Sigma p_{T}^{truth constituents};weighted counts", + {HistType::kTH2F, {jetPtAxis, constituentPtFractionAxis}}); + + auto setDetectorJetExistsLabels = [](const auto& histogram) { + histogram->GetYaxis()->SetBinLabel(1, "no detector jet"); + histogram->GetYaxis()->SetBinLabel(2, "detector jet exists"); + }; + setDetectorJetExistsLabels(registry.get(HIST("h2_mcp_jet_pt_vs_has_detector_jet_hasColl"))); + setDetectorJetExistsLabels(registry.get(HIST("h2_mcp_jet_pt_vs_has_detector_jet_noColl"))); + } } template @@ -369,7 +598,7 @@ struct JetCrossSectionEfficiency { } template - bool configureTruthRofParameters(TBC const& truthBC) + bool configureRofParameters(TBC const& bc) { if (truthRofOffsetInBC >= 0 && truthRofLengthInBC > 0) { cachedTruthRofOffsetInBC = truthRofOffsetInBC; @@ -377,28 +606,28 @@ struct JetCrossSectionEfficiency { return true; } - if (cachedTruthRofRun == truthBC.runNumber() && cachedTruthRofLengthInBC > 0) { + if (cachedTruthRofRun == bc.runNumber() && cachedTruthRofLengthInBC > 0) { return true; } - auto alppar = ccdb->getForTimeStamp>("ITS/Config/AlpideParam", truthBC.timestamp()); + auto alppar = ccdb->getForTimeStamp>("ITS/Config/AlpideParam", bc.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())); + LOGF(fatal, "Could not retrieve ITS/Config/AlpideParam for ITS ROF calculations (run %d, timestamp %" PRIu64 ")", bc.runNumber(), static_cast(bc.timestamp())); return false; } - cachedTruthRofRun = truthBC.runNumber(); + cachedTruthRofRun = bc.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)); + LOGF(fatal, "Invalid ITS ROF length %" PRId64 " BC", 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); + LOGF(info, "ITS ROF calculations use 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 + int64_t rofId(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; @@ -426,11 +655,11 @@ struct JetCrossSectionEfficiency { requireValidTruthFT0CActivity(mccollision); auto truthBC = mccollision.template bc_as(); const uint64_t currentGlobalBC = truthBC.globalBC(); - if (currentGlobalBC == std::numeric_limits::max() || !configureTruthRofParameters(truthBC)) { + if (currentGlobalBC == std::numeric_limits::max() || !configureRofParameters(truthBC)) { return selections; } - const int64_t currentRof = truthRofId(currentGlobalBC); + const int64_t currentRof = rofId(currentGlobalBC); bool hasNarrowActivity = false; bool hasHighActivityInTimeRange = false; bool hasHighActivityInSameRof = false; @@ -455,7 +684,7 @@ struct JetCrossSectionEfficiency { 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; + const bool isSameRof = rofId(otherGlobalBC) == currentRof; if (!inNarrowWindow && !inStandardTimeWindow && !isSameRof) { continue; } @@ -488,6 +717,177 @@ struct JetCrossSectionEfficiency { return selections; } + template + bool passesTVXBCEnvironmentBaseline(TBC const& bc) + { + const bool passesRct = !static_cast(applyRCT) || (bc.rct_raw() & rctMask) == 0; + return passesRct && + bc.selection_bit(aod::evsel::kIsTriggerTVX) && + bc.selection_bit(aod::evsel::kNoTimeFrameBorder) && + bc.selection_bit(aod::evsel::kNoITSROFrameBorder); + } + + template + bool isOtherTVXBCForEnvironment(TBC const& bc) + { + if (!bc.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + return !static_cast(tvxEnvironmentRequireOtherBcBaseline) || passesTVXBCEnvironmentBaseline(bc); + } + + template + VisibleBCEnvironmentCache buildVisibleBCEnvironmentCache(TBCs const& bcs) + { + struct TVXBCInfo { + int64_t globalIndex = -1; + int runNumber = 0; + uint64_t globalBC = 0; + int64_t rof = -1; + bool passesTargetBaseline = false; + }; + + std::map> tvxBCsPerRun; + for (const auto& bc : bcs) { + if (!isOtherTVXBCForEnvironment(bc) || bc.globalBC() == std::numeric_limits::max()) { + continue; + } + if (!configureRofParameters(bc)) { + continue; + } + tvxBCsPerRun[bc.runNumber()].push_back({static_cast(bc.globalIndex()), bc.runNumber(), bc.globalBC(), rofId(bc.globalBC()), passesTVXBCEnvironmentBaseline(bc)}); + } + + VisibleBCEnvironmentCache cache; + for (auto& [runNumber, tvxBCs] : tvxBCsPerRun) { + std::sort(tvxBCs.begin(), tvxBCs.end(), [](const TVXBCInfo& lhs, const TVXBCInfo& rhs) { + return lhs.globalBC < rhs.globalBC; + }); + + std::unordered_map nTVXBCsPerRof; + for (const auto& tvxBC : tvxBCs) { + ++nTVXBCsPerRof[tvxBC.rof]; + } + + for (size_t i = 0; i < tvxBCs.size(); ++i) { + const auto& tvxBC = tvxBCs[i]; + if (!tvxBC.passesTargetBaseline) { + continue; + } + + uint64_t nearestAbsDeltaGlobalBC = std::numeric_limits::max(); + if (i > 0) { + nearestAbsDeltaGlobalBC = tvxBC.globalBC - tvxBCs[i - 1].globalBC; + } + if (i + 1 < tvxBCs.size()) { + const uint64_t nextAbsDeltaGlobalBC = tvxBCs[i + 1].globalBC - tvxBC.globalBC; + nearestAbsDeltaGlobalBC = std::min(nearestAbsDeltaGlobalBC, nextAbsDeltaGlobalBC); + } + + VisibleBCEnvironment environment; + environment.valid = true; + environment.nOtherTVXSameRof = nTVXBCsPerRof[tvxBC.rof] - 1; + environment.originalBCIndex = tvxBC.globalIndex; + if (nearestAbsDeltaGlobalBC != std::numeric_limits::max()) { + environment.nearestAbsDeltaTVXGlobalBC = static_cast(nearestAbsDeltaGlobalBC); + } + cache.byBCIndex[tvxBC.globalIndex] = environment; + cache.byRunAndGlobalBC[{runNumber, tvxBC.globalBC}] = environment; + } + } + return cache; + } + + template + std::map buildTruthHasCollEnvironmentCache(TMcCollisions const& mcCollisions) + { + struct TruthBCInfo { + int64_t mcCollisionIndex = -1; + int runNumber = 0; + uint64_t globalBC = 0; + int64_t rof = -1; + }; + + std::map> truthBCsPerRun; + for (const auto& mccollision : mcCollisions) { + const auto truthBC = mccollision.template bc_as(); + if (truthBC.globalBC() == std::numeric_limits::max() || !configureRofParameters(truthBC)) { + continue; + } + truthBCsPerRun[truthBC.runNumber()].push_back({static_cast(mccollision.globalIndex()), truthBC.runNumber(), truthBC.globalBC(), rofId(truthBC.globalBC())}); + } + + std::map cache; + for (auto& [runNumber, truthBCs] : truthBCsPerRun) { + (void)runNumber; + std::sort(truthBCs.begin(), truthBCs.end(), [](const TruthBCInfo& lhs, const TruthBCInfo& rhs) { + return lhs.globalBC < rhs.globalBC; + }); + + std::unordered_map nTruthBCsPerRof; + for (const auto& truthBC : truthBCs) { + ++nTruthBCsPerRof[truthBC.rof]; + } + + for (size_t i = 0; i < truthBCs.size(); ++i) { + const auto& truthBC = truthBCs[i]; + uint64_t nearestAbsDeltaGlobalBC = std::numeric_limits::max(); + if (i > 0) { + nearestAbsDeltaGlobalBC = truthBC.globalBC - truthBCs[i - 1].globalBC; + } + if (i + 1 < truthBCs.size()) { + const uint64_t nextAbsDeltaGlobalBC = truthBCs[i + 1].globalBC - truthBC.globalBC; + nearestAbsDeltaGlobalBC = std::min(nearestAbsDeltaGlobalBC, nextAbsDeltaGlobalBC); + } + + TruthHasCollEnvironment environment; + environment.valid = true; + environment.nOtherSameTruthRof = nTruthBCsPerRof[truthBC.rof] - 1; + if (nearestAbsDeltaGlobalBC != std::numeric_limits::max()) { + environment.nearestAbsDeltaGlobalBC = static_cast(nearestAbsDeltaGlobalBC); + } + cache[truthBC.mcCollisionIndex] = environment; + } + } + return cache; + } + + template + void fillTVXBCEnvironmentHistogramsFromCache(VisibleBCEnvironmentCache const& environments, + TBCs const& bcs, + CollisionsWithEvSels const& collisions, + aod::FT0s const& ft0s) + { + for (const auto& bc : bcs) { + const auto environment = environments.byBCIndex.find(static_cast(bc.globalIndex())); + if (environment == environments.byBCIndex.end()) { + continue; + } + + // This is a data-level TVX-BC to reconstructed-collision proxy, not the MC truth hasColl definition. + const auto collisionsAtTargetBC = collisions.sliceBy(collisionsPerFoundBC, bc.globalIndex()); + const bool hasRecoCollisionAtTargetBC = collisionsAtTargetBC.size() >= 1; + const auto& visibleEnvironment = environment->second; + registry.fill(HIST("h2_n_other_tvx_same_rof_vs_nearest_abs_delta_global_bc_tvx"), visibleEnvironment.nOtherTVXSameRof, visibleEnvironment.nearestAbsDeltaTVXGlobalBC); + if (hasRecoCollisionAtTargetBC) { + registry.fill(HIST("h2_n_other_tvx_same_rof_vs_nearest_abs_delta_global_bc_tvx_has_reco"), visibleEnvironment.nOtherTVXSameRof, visibleEnvironment.nearestAbsDeltaTVXGlobalBC); + } + registry.fill(HIST("h2_n_other_tvx_same_rof_vs_has_reco_collision"), visibleEnvironment.nOtherTVXSameRof, hasRecoCollisionAtTargetBC); + registry.fill(HIST("h2_nearest_abs_delta_global_bc_tvx_vs_has_reco_collision"), visibleEnvironment.nearestAbsDeltaTVXGlobalBC, hasRecoCollisionAtTargetBC); + if (bc.has_foundFT0()) { + const auto foundFT0 = ft0s.rawIteratorAt(bc.foundFT0Id()); + registry.fill(HIST("h2_target_ft0c_amplitude_vs_n_other_tvx_same_rof"), foundFT0.sumAmpC(), visibleEnvironment.nOtherTVXSameRof); + } + } + } + + template + void fillTVXBCEnvironmentHistograms(TBCs const& bcs, CollisionsWithEvSels const& collisions, aod::FT0s const& ft0s) + { + const auto environments = buildVisibleBCEnvironmentCache(bcs); + fillTVXBCEnvironmentHistogramsFromCache(environments, bcs, collisions, ft0s); + } + template bool isAcceptedJet(TJets const& jet) { @@ -519,10 +919,10 @@ struct JetCrossSectionEfficiency { return true; } - void processTruthMultFT0CvsRecoAmplitude(aod::JetMcCollisions::iterator const& mccollision, - soa::SmallGroups const& collisions, - CollisionsWithEvSels const&, - aod::FT0s const& ft0s) + void processMCTruthMultFT0CvsRecoAmplitude(aod::JetMcCollisions::iterator const& mccollision, + soa::SmallGroups const& collisions, + CollisionsWithEvSels const&, + aod::FT0s const& ft0s) { if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { return; @@ -543,13 +943,13 @@ struct JetCrossSectionEfficiency { 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); + PROCESS_SWITCH(JetCrossSectionEfficiency, processMCTruthMultFT0CvsRecoAmplitude, + "[MC ONLY] 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&) + void processMCTruthMultFT0CvsRecoITS567(aod::JetMcCollisions::iterator const& mccollision, + soa::SmallGroups const& collisions, + CollisionsWithEvSels const&, + FullTracksIU const&) { if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { return; @@ -576,8 +976,290 @@ struct JetCrossSectionEfficiency { 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); + PROCESS_SWITCH(JetCrossSectionEfficiency, processMCTruthMultFT0CvsRecoITS567, + "[MC ONLY] truth multFT0C activity vs reconstructed PV-contributor ITS567 activity for one-to-one MC/reco collision associations", false); + + void processMCEnvironmentQA(aod::JetMcCollisions const& mcCollisions, + aod::JetCollisionsMCD const& jetCollisions, + aod::JBCs const&, + OriginalBCsWithSels const& bcs, + CollisionsWithEvSels const& collisions, + aod::FT0s const& ft0s) + { + const auto visibleEnvironments = buildVisibleBCEnvironmentCache(bcs); + const auto truthEnvironments = buildTruthHasCollEnvironmentCache(mcCollisions); + // The visible-BC QA is deliberately unweighted so its definition remains directly comparable to data. + fillTVXBCEnvironmentHistogramsFromCache(visibleEnvironments, bcs, collisions, ft0s); + for (const auto& mccollision : mcCollisions) { + if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + continue; + } + + const auto truthBC = mccollision.template bc_as(); + if (!passesTVXBCEnvironmentBaseline(truthBC)) { + continue; + } + + const auto environment = truthEnvironments.find(static_cast(mccollision.globalIndex())); + if (environment == truthEnvironments.end() || !environment->second.valid) { + continue; + } + + const bool hasRecoColl = jetCollisions.sliceBy(jetCollisionsPerMcCollision, mccollision.globalIndex()).size() >= 1; + const float weight = mccollision.weight(); + registry.fill(HIST("h2_n_other_same_truth_rof_vs_has_coll"), environment->second.nOtherSameTruthRof, hasRecoColl, weight); + registry.fill(HIST("h2_nearest_abs_delta_global_bc_vs_has_coll"), environment->second.nearestAbsDeltaGlobalBC, hasRecoColl, weight); + + const auto visibleEnvironment = visibleEnvironments.byRunAndGlobalBC.find({truthBC.runNumber(), truthBC.globalBC()}); + if (visibleEnvironment != visibleEnvironments.byRunAndGlobalBC.end() && visibleEnvironment->second.valid) { + const auto collisionsAtTargetBC = collisions.sliceBy(collisionsPerFoundBC, visibleEnvironment->second.originalBCIndex); + const bool hasRecoCollisionAtTargetBC = collisionsAtTargetBC.size() >= 1; + registry.fill(HIST("h2_n_other_truth_same_rof_vs_n_other_tvx_same_rof"), environment->second.nOtherSameTruthRof, visibleEnvironment->second.nOtherTVXSameRof, weight); + registry.fill(HIST("h2_nearest_abs_delta_global_bc_truth_vs_nearest_abs_delta_global_bc_tvx"), environment->second.nearestAbsDeltaGlobalBC, visibleEnvironment->second.nearestAbsDeltaTVXGlobalBC, weight); + registry.fill(HIST("h2_reco_collision_proxy_vs_has_coll"), hasRecoCollisionAtTargetBC, hasRecoColl, weight); + } + + if (mccollision.multFT0C() < 0.0f) { + continue; + } + + registry.fill(HIST("h2_target_mult_ft0c_vs_has_coll"), mccollision.multFT0C(), hasRecoColl, weight); + registry.fill(HIST("h2_target_mult_ft0c_vs_n_other_same_truth_rof_all"), mccollision.multFT0C(), environment->second.nOtherSameTruthRof, weight); + registry.fill(HIST("h2_target_mult_ft0c_vs_nearest_abs_delta_global_bc_all"), mccollision.multFT0C(), environment->second.nearestAbsDeltaGlobalBC, weight); + if (hasRecoColl) { + registry.fill(HIST("h2_target_mult_ft0c_vs_n_other_same_truth_rof_has_coll"), mccollision.multFT0C(), environment->second.nOtherSameTruthRof, weight); + registry.fill(HIST("h2_target_mult_ft0c_vs_nearest_abs_delta_global_bc_has_coll"), mccollision.multFT0C(), environment->second.nearestAbsDeltaGlobalBC, weight); + } + } + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processMCEnvironmentQA, + "[MC ONLY] truth and visible TVX collision-environment QA and JetDerived hasColl validation", false); + + void processMCHasCollTrackFateQA(JMcCollisionsWithParent::iterator const& mccollision, + soa::SmallGroups const& collisions, + OriginalCollisionsWithMcLabels const&, + FullTracksIUWithMcLabels const& recoTracks, + aod::McCollisions const&, + aod::McParticles const& mcParticles, + aod::JBCs const&) + { + if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + return; + } + + auto truthBC = mccollision.bc_as(); + const bool passesRct = !applyRCT || (truthBC.rct_raw() & rctMask) == 0; + const bool passesTVXTruth = truthBC.selection_bit(aod::evsel::kIsTriggerTVX); + const bool passesNoTFBTruth = truthBC.selection_bit(aod::evsel::kNoTimeFrameBorder); + const bool passesNoITSROFBTruth = truthBC.selection_bit(aod::evsel::kNoITSROFrameBorder); + if (!passesRct || !passesTVXTruth || !passesNoTFBTruth || !passesNoITSROFBTruth) { + return; + } + + const auto originalMcCollision = mccollision.template mcCollision_as(); + const auto targetMcParticles = mcParticles.sliceBy(mcParticlesPerMcCollision, originalMcCollision.globalIndex()); + int nRecoTracksFromTarget = 0; + int nUnassignedTracksFromTarget = 0; + int nTracksAttachedToOtherLabeledRecoCollision = 0; + + for (const auto& mcParticle : targetMcParticles) { + const auto targetRecoTracks = recoTracks.sliceBy(recoTracksPerMcParticle, mcParticle.globalIndex()); + for (const auto& recoTrack : targetRecoTracks) { + if (!recoTrack.has_mcParticle()) { + continue; + } + + ++nRecoTracksFromTarget; + if (!recoTrack.has_collision()) { + ++nUnassignedTracksFromTarget; + continue; + } + + const auto recoCollision = recoTrack.template collision_as(); + if (recoCollision.has_mcCollision() && recoCollision.mcCollisionId() != originalMcCollision.globalIndex()) { + ++nTracksAttachedToOtherLabeledRecoCollision; + } + } + } + + const bool hasRecoColl = collisions.size() >= 1; + const float weight = mccollision.weight(); + registry.fill(HIST("h2_n_reco_tracks_from_target_vs_has_coll"), nRecoTracksFromTarget, hasRecoColl, weight); + registry.fill(HIST("h2_n_unassigned_tracks_from_target_vs_has_coll"), nUnassignedTracksFromTarget, hasRecoColl, weight); + registry.fill(HIST("h2_n_target_tracks_attached_to_other_labeled_collision_vs_has_coll"), nTracksAttachedToOtherLabeledRecoCollision, hasRecoColl, weight); + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processMCHasCollTrackFateQA, + "[MC ONLY] track fate for target MC collisions before JetDerived hasColl", false); + + void processMCHasCollJetSurvivalQA(JMcCollisionsWithParent::iterator const& mccollision, + soa::SmallGroups const& collisions, + ChargedMCPJetsWithConstituents const& mcpJets, + ChargedMCDJetsWithConstituents const& detectorJets, + JetCollisionsWithParent const& jetCollisions, + FullTracksIUWithMcLabels const& recoTracks, + aod::JetTracks const&, + aod::JTrackPIs const& jTrackParents, + aod::JetParticles const&, + aod::JMcParticlePIs const& jMcParticleParents, + aod::JBCs const&) + { + if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + return; + } + + auto truthBC = mccollision.bc_as(); + const bool passesRct = !applyRCT || (truthBC.rct_raw() & rctMask) == 0; + const bool passesTVXTruth = truthBC.selection_bit(aod::evsel::kIsTriggerTVX); + const bool passesNoTFBTruth = truthBC.selection_bit(aod::evsel::kNoTimeFrameBorder); + const bool passesNoITSROFBTruth = truthBC.selection_bit(aod::evsel::kNoITSROFrameBorder); + if (!passesRct || !passesTVXTruth || !passesNoTFBTruth || !passesNoITSROFBTruth) { + return; + } + + const float pTHat = computePtHat(mccollision); + if (pTHat < pTHatAbsoluteMin) { + return; + } + + const bool hasRecoColl = collisions.size() >= 1; + const float weight = mccollision.weight(); + for (const auto& mcpJet : mcpJets) { + if (!jetfindingutilities::isInEtaAcceptance(mcpJet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || + mcpJet.pt() < jetPtMin || mcpJet.pt() > pTHatMaxMCP * pTHat || + !isAcceptedJet(mcpJet)) { + continue; + } + + float sumTruthConstituentPt = 0.0f; + float sumRecoTrackPtFromTruthJet = 0.0f; + std::unordered_map recoTrackPtByOriginalTrackId; + std::unordered_set jetCollisionIdsReceivingTargetTracks; + for (const auto& constituent : mcpJet.template tracks_as()) { + sumTruthConstituentPt += constituent.pt(); + const auto originalMcParticle = jMcParticleParents.rawIteratorAt(constituent.globalIndex()); + const auto targetRecoTracks = recoTracks.sliceBy(recoTracksPerMcParticle, originalMcParticle.mcParticleId()); + for (const auto& recoTrack : targetRecoTracks) { + if (!recoTrack.has_mcParticle()) { + continue; + } + + const int64_t originalTrackId = recoTrack.globalIndex(); + if (!recoTrackPtByOriginalTrackId.emplace(originalTrackId, recoTrack.pt()).second) { + continue; + } + sumRecoTrackPtFromTruthJet += recoTrack.pt(); + + if (!recoTrack.has_collision()) { + continue; + } + const auto associatedJetCollisions = jetCollisions.sliceBy(jetCollisionsPerOriginalCollision, recoTrack.collisionId()); + for (const auto& jetCollision : associatedJetCollisions) { + jetCollisionIdsReceivingTargetTracks.insert(jetCollision.globalIndex()); + } + } + } + + if (sumTruthConstituentPt <= 0.0f) { + continue; + } + + std::unordered_map capturedRecoPtByDetectorJetId; + for (const int64_t jetCollisionId : jetCollisionIdsReceivingTargetTracks) { + const auto associatedDetectorJets = detectorJets.sliceBy(detectorJetsPerJetCollision, jetCollisionId); + for (const auto& detectorJet : associatedDetectorJets) { + float capturedRecoPt = 0.0f; + for (const auto& detectorConstituent : detectorJet.template tracks_as()) { + const auto originalTrack = jTrackParents.rawIteratorAt(detectorConstituent.globalIndex()); + const auto recoTrack = recoTrackPtByOriginalTrackId.find(originalTrack.trackId()); + if (recoTrack != recoTrackPtByOriginalTrackId.end()) { + capturedRecoPt += recoTrack->second; + } + } + if (capturedRecoPt > 0.0f) { + capturedRecoPtByDetectorJetId[detectorJet.globalIndex()] += capturedRecoPt; + } + } + } + + float bestDetectorJetCapturedTargetRecoPt = 0.0f; + for (const auto& [detectorJetId, capturedRecoPt] : capturedRecoPtByDetectorJetId) { + (void)detectorJetId; + if (capturedRecoPt > bestDetectorJetCapturedTargetRecoPt) { + bestDetectorJetCapturedTargetRecoPt = capturedRecoPt; + } + } + + const bool hasDetectorJetWithTargetTrack = !capturedRecoPtByDetectorJetId.empty(); + const float recoConstituentPtFraction = sumRecoTrackPtFromTruthJet / sumTruthConstituentPt; + const float bestDetectorJetCaptureFraction = bestDetectorJetCapturedTargetRecoPt / sumTruthConstituentPt; + if (hasRecoColl) { + registry.fill(HIST("h2_mcp_jet_pt_vs_has_detector_jet_hasColl"), mcpJet.pt(), hasDetectorJetWithTargetTrack, weight); + registry.fill(HIST("h2_mcp_jet_pt_vs_reco_constituent_pt_fraction_hasColl"), mcpJet.pt(), recoConstituentPtFraction, weight); + registry.fill(HIST("h2_mcp_jet_pt_vs_best_detector_jet_capture_fraction_hasColl"), mcpJet.pt(), bestDetectorJetCaptureFraction, weight); + } else { + registry.fill(HIST("h2_mcp_jet_pt_vs_has_detector_jet_noColl"), mcpJet.pt(), hasDetectorJetWithTargetTrack, weight); + registry.fill(HIST("h2_mcp_jet_pt_vs_reco_constituent_pt_fraction_noColl"), mcpJet.pt(), recoConstituentPtFraction, weight); + registry.fill(HIST("h2_mcp_jet_pt_vs_best_detector_jet_capture_fraction_noColl"), mcpJet.pt(), bestDetectorJetCaptureFraction, weight); + } + } + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processMCHasCollJetSurvivalQA, + "[MC ONLY] truth-constituent overlap QA for detector-level charged jet survival before JetDerived hasColl", false); + + void processMCHasCollJetEnvironmentQA(aod::JetMcCollisions const& mcCollisions, + aod::JetCollisionsMCD const& jetCollisions, + ChargedMCPJetsWithConstituents const& mcpJets, + aod::JetParticles const&, + aod::JBCs const&, + OriginalBCsWithSels const& bcs) + { + const auto visibleEnvironments = buildVisibleBCEnvironmentCache(bcs); + for (const auto& mccollision : mcCollisions) { + if (skipMBGapEvents && mccollision.getSubGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) { + continue; + } + + const auto truthBC = mccollision.template bc_as(); + if (!passesTVXBCEnvironmentBaseline(truthBC)) { + continue; + } + + const auto visibleEnvironment = visibleEnvironments.byRunAndGlobalBC.find({truthBC.runNumber(), truthBC.globalBC()}); + if (visibleEnvironment == visibleEnvironments.byRunAndGlobalBC.end() || !visibleEnvironment->second.valid) { + continue; + } + + const float pTHat = computePtHat(mccollision); + if (pTHat < pTHatAbsoluteMin) { + continue; + } + + const bool hasRecoColl = jetCollisions.sliceBy(jetCollisionsPerMcCollision, mccollision.globalIndex()).size() >= 1; + const float weight = mccollision.weight(); + const auto targetMcpJets = mcpJets.sliceBy(mcpJetsPerMcCollision, mccollision.globalIndex()); + for (const auto& mcpJet : targetMcpJets) { + if (!jetfindingutilities::isInEtaAcceptance(mcpJet, jetEtaMin, jetEtaMax, trackEtaMin, trackEtaMax) || + mcpJet.pt() < jetPtMin || mcpJet.pt() > pTHatMaxMCP * pTHat || + !isAcceptedJet(mcpJet)) { + continue; + } + + // hasRecoColl remains a histogram axis rather than a selection, so noColl MCP jets are part of the denominator. + registry.fill(HIST("hs_mcp_jet_pt_tvx_environment_has_coll"), mcpJet.pt(), visibleEnvironment->second.nOtherTVXSameRof, visibleEnvironment->second.nearestAbsDeltaTVXGlobalBC, hasRecoColl, weight); + } + } + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processMCHasCollJetEnvironmentQA, + "[MC ONLY] selected MCP jets vs original-TVX-BC environment and JetDerived hasColl", false); + + void processDataTVXBCEnvironmentQA(OriginalBCsWithSels const& bcs, + CollisionsWithEvSels const& collisions, + aod::FT0s const& ft0s) + { + fillTVXBCEnvironmentHistograms(bcs, collisions, ft0s); + } + PROCESS_SWITCH(JetCrossSectionEfficiency, processDataTVXBCEnvironmentQA, + "[DATA ONLY] TVX-BC collision-environment QA using original BC and reconstructed-collision tables", false); void processCrossSectionEfficiency(aod::JetMcCollisions::iterator const& mccollision, soa::SmallGroups const& collisions, From ec27efea75b434474b1a8133cce72c3b513d7445 Mon Sep 17 00:00:00 2001 From: Wooseok Ham Date: Tue, 29 Sep 2026 18:11:12 +0900 Subject: [PATCH 2/2] PWGJE: fix the data type error --- PWGJE/Tasks/jetCrossSectionEfficiency.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGJE/Tasks/jetCrossSectionEfficiency.cxx b/PWGJE/Tasks/jetCrossSectionEfficiency.cxx index 874abdea4fb..cfc984b2c50 100644 --- a/PWGJE/Tasks/jetCrossSectionEfficiency.cxx +++ b/PWGJE/Tasks/jetCrossSectionEfficiency.cxx @@ -35,6 +35,7 @@ #include #include #include +#include #include #include #include @@ -1000,7 +1001,7 @@ struct JetCrossSectionEfficiency { continue; } - const auto environment = truthEnvironments.find(static_cast(mccollision.globalIndex())); + const auto environment = truthEnvironments.find(mccollision.globalIndex()); if (environment == truthEnvironments.end() || !environment->second.valid) { continue; }