Skip to content

[PWGLF] NucleiTask - Change reco condition on MC #17953

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mpuccio merged 1 commit into from
Sep 18, 2026
+132 −22
Merged

[PWGLF] NucleiTask - Change reco condition on MC #17953

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
154 changes: 132 additions & 22 deletions PWGLF/Tasks/Nuspex/lfNucleiBATask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,8 @@ struct lfNucleiBATask {

// Efficiency configurator
std::unordered_set<int> effEvtSet;
bool effEvtSetReady = false;
Configurable<bool> enableEffEvtSet{"enableEffEvtSet", true, "If true, MCGen uses the event-set built by MCReco; if false, MCGen runs stand-alone."};
// bool effEvtSetReady = false;
Configurable<bool> 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};
Expand Down Expand Up @@ -112,6 +112,7 @@ struct lfNucleiBATask {
Configurable<bool> useINELgt0cut{"useINELgt0cut", false, "Apply INEL>0 Event Selection in processData o processMC"};
Configurable<bool> useINELgt1cut{"useINELgt1cut", false, "Apply INEL>1 Event Selection in processData o processMC"};
Configurable<bool> enableGenVzCut{"enableGenVzCut", true, "Apply cut in z-Vertex in the processGen function"};
Configurable<bool> askGenFT0Coincidence{"askGenFT0Coincidence", false, "Require a generated charged-particle coincidence in FT0A and FT0C"};
} evselOptions;

// Set the multiplity event limits
Expand Down Expand Up @@ -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}, ""};
Expand Down Expand Up @@ -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) {
Expand All @@ -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})"};
Expand Down Expand Up @@ -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<TH2>(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<TH1>(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
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -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();
Expand Down Expand Up @@ -6582,21 +6631,53 @@ struct lfNucleiBATask {
////////////

// LOOP OVER GENERATED MC PARTICLES
void processMCGen(soa::Join<aod::McCollisions, aod::McCentFT0Ms>::iterator const& mcCollision,
aod::McParticles const& mcParticles)
// void processMCGen(soa::Join<aod::McCollisions, aod::McCentFT0Ms>::iterator const& mcCollision,
// aod::McParticles const& mcParticles)
void processMCGen(
soa::Join<aod::McCollisions, aod::McCentFT0Ms>::iterator const& mcCollision,
const soa::SmallGroups<soa::Join<EventCandidatesMC, o2::aod::PVMults>>& 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;
Expand All @@ -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;
Expand Down Expand Up @@ -6958,17 +7036,25 @@ struct lfNucleiBATask {
const soa::SmallGroups<soa::Join<EventCandidatesMC, o2::aod::PVMults>>& 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);
if (enableCentrality)
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)
Expand Down Expand Up @@ -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;
Expand Down
Loading