diff --git a/PWGLF/Tasks/Nuspex/lfNucleiBATask.cxx b/PWGLF/Tasks/Nuspex/lfNucleiBATask.cxx index 9cfb1e69a55..27d7abcbddc 100644 --- a/PWGLF/Tasks/Nuspex/lfNucleiBATask.cxx +++ b/PWGLF/Tasks/Nuspex/lfNucleiBATask.cxx @@ -77,8 +77,8 @@ struct lfNucleiBATask { // Efficiency configurator std::unordered_set effEvtSet; - bool effEvtSetReady = false; - Configurable enableEffEvtSet{"enableEffEvtSet", true, "If true, MCGen uses the event-set built by MCReco; if false, MCGen runs stand-alone."}; + // bool effEvtSetReady = false; + Configurable requireRecoEventForGen{"requireRecoEventForGen", true, "Require at least one reco collision passing the event selection in processMCGen"}; HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry spectraGen{"spectraGen", {}, OutputObjHandlingPolicy::AnalysisObject, false, true}; @@ -112,6 +112,7 @@ struct lfNucleiBATask { Configurable useINELgt0cut{"useINELgt0cut", false, "Apply INEL>0 Event Selection in processData o processMC"}; Configurable useINELgt1cut{"useINELgt1cut", false, "Apply INEL>1 Event Selection in processData o processMC"}; Configurable enableGenVzCut{"enableGenVzCut", true, "Apply cut in z-Vertex in the processGen function"}; + Configurable askGenFT0Coincidence{"askGenFT0Coincidence", false, "Require a generated charged-particle coincidence in FT0A and FT0C"}; } evselOptions; // Set the multiplity event limits @@ -167,9 +168,9 @@ struct lfNucleiBATask { // Set the axis used in this task ConfigurableAxis binsPercentile{"binsPercentile", {100, 0, 100}, "Centrality FT0M"}; ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.425, 0.45, 0.475, 0.5, 0.5125, 0.525, 0.5375, 0.55, 0.5625, 0.575, 0.5875, 0.6, 0.6125, 0.625, 0.6375, 0.65, 0.6625, 0.675, 0.6875, 0.7, 0.7125, 0.725, 0.7375, 0.75, 0.7625, 0.775, 0.7875, 0.8, 0.8125, 0.825, 0.8375, 0.85, 0.8625, 0.875, 0.8875, 0.9, 0.9125, 0.925, 0.9375, 0.95, 0.9625, 0.975, 0.9875, 1.0, 1.0125, 1.025, 1.0375, 1.05, 1.0625, 1.075, 1.0875, 1.1, 1.1125, 1.125, 1.1375, 1.15, 1.1625, 1.175, 1.1875, 1.2, 1.2125, 1.225, 1.2375, 1.25, 1.2625, 1.275, 1.2875, 1.3, 1.3125, 1.325, 1.3375, 1.35, 1.3625, 1.375, 1.3875, 1.4, 1.4125, 1.425, 1.4375, 1.45, 1.4625, 1.475, 1.4875, 1.5, 1.5125, 1.525, 1.5375, 1.55, 1.5625, 1.575, 1.5875, 1.6, 1.6125, 1.625, 1.6375, 1.65, 1.6625, 1.675, 1.6875, 1.7, 1.7125, 1.725, 1.7375, 1.75, 1.7625, 1.775, 1.7875, 1.8, 1.8125, 1.825, 1.8375, 1.85, 1.8625, 1.875, 1.8875, 1.9, 1.9125, 1.925, 1.9375, 1.95, 1.9625, 1.975, 1.9875, 2.0, 2.0625, 2.125, 2.1875, 2.25, 2.3125, 2.375, 2.4375, 2.5, 2.625, 2.75, 2.875, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0}, ""}; - ConfigurableAxis binsPtHe{"binsPtHe", {VARIABLE_WIDTH, 1.0, 1.25, 1.50, 1.75, 2.0, 2.25, 2.50, 2.75, 3.0, 3.25, 3.50, 3.75, 4.0, 4.50, 5.0, 6.0, 7.0, 8.0}, ""}; - ConfigurableAxis binsPtZHe{"binsPtZHe", {VARIABLE_WIDTH, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2.0, 2.25, 2.5, 3.0, 3.5, 4.0}, ""}; - ConfigurableAxis binsPtHeMigration{"binsPtHeMigration", {VARIABLE_WIDTH, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.25, 3.50, 3.75, 4.00, 4.50, 5.00, 6.00, 7.00, 8.00}, ""}; + ConfigurableAxis binsPtHe{"binsPtHe", {VARIABLE_WIDTH, 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.50, 1.75, 2.0, 2.25, 2.50, 2.75, 3.0, 3.25, 3.50, 3.75, 4.0, 4.50, 5.0, 6.0, 7.0, 8.0}, ""}; + ConfigurableAxis binsPtZHe{"binsPtZHe", {VARIABLE_WIDTH, 0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0, 1.125, 1.25, 1.375, 1.5, 1.625, 1.75, 1.875, 2.0, 2.25, 2.5, 3.0, 3.5, 4.0}, ""}; + ConfigurableAxis binsPtHeMigration{"binsPtHeMigration", {VARIABLE_WIDTH, 0.0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.25, 3.50, 3.75, 4.00, 4.50, 5.00, 6.00, 7.00, 8.00}, ""}; ConfigurableAxis binsdEdx{"binsdEdx", {600, 0.f, 3000.f}, ""}; ConfigurableAxis binsBeta{"binsBeta", {120, 0.0, 1.2}, ""}; @@ -323,6 +324,38 @@ struct lfNucleiBATask { return 0.5f * std::log((e + pz) / (e - pz)); } + // Added TVX check on MC + bool hasGenTVX(const aod::McParticles& particles, int mcCollisionIndex) + { + bool hasHitFT0A = false; + bool hasHitFT0C = false; + + for (const auto& particle : particles) { + if (particle.mcCollisionId() != mcCollisionIndex) + continue; + + if (!particle.isPhysicalPrimary()) + continue; + + auto* pdgInfo = pdgDB->GetParticle(particle.pdgCode()); + if (!pdgInfo || pdgInfo->Charge() == 0.) + continue; + + const float eta = particle.eta(); + + if (eta > 3.5f && eta < 4.9f) + hasHitFT0A = true; + + if (eta > -3.3f && eta < -2.1f) + hasHitFT0C = true; + + if (hasHitFT0A && hasHitFT0C) + return true; + } + + return false; + } + void init(o2::framework::InitContext& context) { if (initITSPID) { @@ -333,7 +366,7 @@ struct lfNucleiBATask { } effEvtSet.clear(); - effEvtSetReady = false; + // effEvtSetReady = false; const AxisSpec pAxis{binsPt, "#it{p} (GeV/#it{c})"}; const AxisSpec ptAxis{binsPt, "#it{p}_{T} (GeV/#it{c})"}; @@ -2279,6 +2312,21 @@ struct lfNucleiBATask { } LOG(info) << " MC mode "; // MC histograms - all, primary, sec. from weak decay, sec. from material + + if (enableDebug) { + spectraGen.add("hEffEvtSetVsSG", "effEvtSet vs SmallGroups;in effEvtSet;has accepted reco collision", HistType::kTH2D, {{2, -0.5, 1.5}, {2, -0.5, 1.5}}); + auto hEffCheck = spectraGen.get(HIST("hEffEvtSetVsSG")); + hEffCheck->GetXaxis()->SetBinLabel(1, "Not in effEvtSet"); + hEffCheck->GetXaxis()->SetBinLabel(2, "In effEvtSet"); + hEffCheck->GetYaxis()->SetBinLabel(1, "No accepted reco"); + hEffCheck->GetYaxis()->SetBinLabel(2, "Accepted reco"); + + spectraGen.add("hEffEvtSetFalsePositive", "False positives in effEvtSet", HistType::kTH1D, {{2, 0., 2.}}); + auto hFalsePos = spectraGen.get(HIST("hEffEvtSetFalsePositive")); + hFalsePos->GetXaxis()->SetBinLabel(1, "No associated reco"); + hFalsePos->GetXaxis()->SetBinLabel(2, "Reco exists, fails cuts"); + } + if (enableCentrality) spectraGen.add("histGenVetxZ", "PosZ generated events", HistType::kTH2F, {{1500, -15.f, 15.f, "Vertex Z (cm)"}, {binsPercentile, "Centrality FT0M"}}); else @@ -2415,7 +2463,6 @@ struct lfNucleiBATask { } LOG(info) << "MC Histograms defined"; - // LOG(info) << "Histograms of lfNucleiBATask:"; LOG(info) << "Histograms of lfNucleiBATask:"; histos.print(); spectraGen.print(); @@ -6382,12 +6429,13 @@ struct lfNucleiBATask { } float mcCentFT0M = -1.f; + const bool passesRecoVertex = event.posZ() >= cfgVzCutLow && event.posZ() <= cfgVzCutHigh; - if (doRecoSep && event.has_mcCollision()) { + if (doRecoSep && passesRecoVertex && event.has_mcCollision()) { const int mcIdx = event.mcCollisionId(); if (mcIdx >= 0) { effEvtSet.insert(mcIdx); - effEvtSetReady = true; + // effEvtSetReady = true; auto mcColIter = mcCollisions.iteratorAt(mcIdx); mcCentFT0M = mcColIter.centFT0M(); @@ -6423,12 +6471,13 @@ struct lfNucleiBATask { } float mcCentFT0M = -1.f; + const bool passesRecoVertex = event.posZ() >= cfgVzCutLow && event.posZ() <= cfgVzCutHigh; - if (doRecoSep && event.has_mcCollision()) { + if (doRecoSep && passesRecoVertex && event.has_mcCollision()) { const int mcIdx = event.mcCollisionId(); if (mcIdx >= 0) { effEvtSet.insert(mcIdx); - effEvtSetReady = true; + // effEvtSetReady = true; auto mcColIter = mcCollisions.iteratorAt(mcIdx); mcCentFT0M = mcColIter.centFT0M(); @@ -6582,21 +6631,53 @@ struct lfNucleiBATask { //////////// // LOOP OVER GENERATED MC PARTICLES - void processMCGen(soa::Join::iterator const& mcCollision, - aod::McParticles const& mcParticles) + // void processMCGen(soa::Join::iterator const& mcCollision, + // aod::McParticles const& mcParticles) + void processMCGen( + soa::Join::iterator const& mcCollision, + const soa::SmallGroups>& collisions, + aod::McParticles const& mcParticles) { - // Only events that are reconstructed const int mcIdx = mcCollision.globalIndex(); - if (enableEffEvtSet) { - if (!effEvtSetReady) + + // Add GenTVX condition + if (evselOptions.askGenFT0Coincidence) { + if (!hasGenTVX(mcParticles, mcIdx)) return; - if (!effEvtSet.contains(mcIdx)) + } + + // Require Reco conditions (as in signal loss) + if (requireRecoEventForGen) { + bool hasAcceptedRecoCollision = false; + for (const auto& collision : collisions) { + const bool hasTVX = collision.selection_bit(aod::evsel::kIsTriggerTVX); + const bool hasNoTFB = collision.selection_bit(aod::evsel::kNoTimeFrameBorder); + const bool hasNoItsRofFB = collision.selection_bit(aod::evsel::kNoITSROFrameBorder); + if (evselOptions.useSel8 && !collision.sel8()) + continue; + if (evselOptions.useTVXtrigger && !hasTVX) + continue; + if (evselOptions.removeTFBorder && !hasNoTFB) + continue; + if (evselOptions.removeITSROFBorder && !hasNoItsRofFB) + continue; + if (collision.posZ() < cfgVzCutLow || collision.posZ() > cfgVzCutHigh) + continue; + if (evselOptions.useINELgt0cut && !collision.isInelGt0()) + continue; + if (evselOptions.useINELgt1cut && !collision.isInelGt1()) + continue; + hasAcceptedRecoCollision = true; + break; + } + if (!hasAcceptedRecoCollision) return; } if (enableCentrality && (mcCollision.centFT0M() < cfgMultCutLow || mcCollision.centFT0M() > cfgMultCutHigh)) return; + // Vz cut at Gen level if (evselOptions.enableGenVzCut) { if (mcCollision.posZ() < cfgVzCutLow || mcCollision.posZ() > cfgVzCutHigh) return; @@ -6607,9 +6688,6 @@ struct lfNucleiBATask { else spectraGen.fill(HIST("histGenVetxZ"), mcCollision.posZ()); - // const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcIdx, cache); - // const auto particlesInCollision = mcParticles.sliceBy(perMCCol, mcIdx); - for (const auto& mcParticleGen : mcParticles) { if (mcParticleGen.mcCollisionId() != mcIdx) continue; @@ -6958,8 +7036,14 @@ struct lfNucleiBATask { const soa::SmallGroups>& collisions, aod::McParticles const& mcParticles) { + const int mcIdx = mcCollision.globalIndex(); + const bool genFT0Coincidence = hasGenTVX(mcParticles, mcIdx); bool isINELgt0true = pwglf::isINELgtNmc(mcParticles, 0, pdgDB); + // GenFT0 coincidence + if (evselOptions.askGenFT0Coincidence && !genFT0Coincidence) + return; + // EVENT LOSS DENOMINATOR // No cuts histoGen.fill(HIST("events/hMCGen"), 0.5); @@ -6967,8 +7051,10 @@ struct lfNucleiBATask { histoGen.fill(HIST("events/hMCGenVsMult"), 0.5, mcCollision.centFT0M()); // Vtz cut - if (mcCollision.posZ() < cfgVzCutLow || mcCollision.posZ() > cfgVzCutHigh) - return; + if (evselOptions.enableGenVzCut) { + if (mcCollision.posZ() < cfgVzCutLow || mcCollision.posZ() > cfgVzCutHigh) + return; + } histoGen.fill(HIST("events/hMCGen"), 1.5); if (enableCentrality) @@ -7159,6 +7245,30 @@ struct lfNucleiBATask { histoGen.fill(HIST("events/hNRecoCollPerMCCollVsMult_INEL"), recoIdxINEL, mcCollision.centFT0M()); histoGen.fill(HIST("events/hNRecoCollPerMCCollVsMult_INELgt0"), recoIdxINELgt0, mcCollision.centFT0M()); } + // Add checker for EffEvtSet + const int mcIdxCheck = mcCollision.globalIndex(); + const bool inEffEvtSet = effEvtSet.contains(mcIdxCheck); + + if (enableDebug) { + double effEvtSetValue = 0.0; + double inelValue = 0.0; + + if (inEffEvtSet) { + effEvtSetValue = 1.0; + } + if (atLeastOneINEL) { + inelValue = 1.0; + } + + spectraGen.fill(HIST("hEffEvtSetVsSG"), effEvtSetValue, inelValue); + + if (inEffEvtSet && !atLeastOneINEL) { + if (collisions.size() == 0) + spectraGen.fill(HIST("hEffEvtSetFalsePositive"), 0.5); + else + spectraGen.fill(HIST("hEffEvtSetFalsePositive"), 1.5); + } + } if (!atLeastOneINEL) return;