Skip to content

[PWGLF/NuSpEx] Added process function to compute event loss MC #15373

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
maciacco merged 3 commits into from
Mar 17, 2026
+182 −13
Merged

[PWGLF/NuSpEx] Added process function to compute event loss MC #15373

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d5f63b6
Added process function to compute event loss
mapalhares Mar 13, 2026
233bc9a
Please consider the following formatting changes
alibuild Mar 13, 2026
8068253
Merge pull request #62 from alibuild/alibot-cleanup-15373
mapalhares Mar 13, 2026
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
195 changes: 182 additions & 13 deletions PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@
#include "MathUtils/BetheBlochAleph.h"
#include "ReconstructionDataFormats/Track.h"

#include "Math/Vector4D.h"

#include <algorithm>
#include <array>
#include <memory>
Expand All @@ -49,13 +51,18 @@
using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;
using std::array;

using CollBracket = o2::math_utils::Bracket<int>;
using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TOFSignal, aod::TOFEvTime>;
using CollisionsFull = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>;
using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>;

using CollisionsFullWithFlow = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms, aod::FT0Mults, aod::FV0Mults, aod::TPCMults, aod::EPCalibrationTables>;

using McCollisionMults = soa::Join<aod::McCollisions, aod::MultMCExtras>;
using EventCandidatesMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As, aod::Mults>;

namespace
{
constexpr double betheBlochDefault[1][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}};
Expand All @@ -77,6 +84,20 @@ std::shared_ptr<TH1> hH4LMassTracked;
std::shared_ptr<TH1> hDecayChannel;
std::shared_ptr<TH1> hIsMatterGen;
std::shared_ptr<TH1> hIsMatterGenTwoBody;
std::shared_ptr<TH1> hEvtMC;
std::shared_ptr<TH1> hImpactParamGen;
std::shared_ptr<TH1> hImpactParamReco;
std::shared_ptr<TH1> hGen3HLBeforeEvtSel;
std::shared_ptr<TH1> hGen3HLAfterSel;
std::shared_ptr<TH2> hRecoCentralityColvsMultiplicityRecoEta05;
std::shared_ptr<TH2> hRecoCentralityColvsImpactParamReco;
std::shared_ptr<TH2> hGen3HLvsImpactParameterBeforeEvtSel;
std::shared_ptr<TH2> hGen3HLvsImpactParameterAfterSel;
std::shared_ptr<TH2> hGen3HLvsMultiplicityGenEta05BeforeEvtSel;
std::shared_ptr<TH2> hGen3HLvsMultiplicityGenEta05AfterSel;
std::shared_ptr<TH2> hGen3HLvsMultiplicityFT0CBeforeEvtSel;
std::shared_ptr<TH2> hGen3HLvsMultiplicityFT0CAfterSel;

} // namespace

struct hyperCandidate {
Expand Down Expand Up @@ -147,6 +168,7 @@ struct hyperRecoTask {
// PDG codes
Configurable<int> hyperPdg{"hyperPDG", 1010010030, "PDG code of the hyper-mother (could be 3LamH or 4LamH)"};
Configurable<int> heDauPdg{"heDauPDG", 1000020030, "PDG code of the helium (could be 3He or 4He)"};
Configurable<int> piDauPdg{"piDauPdg", 211, "PDG code of pion"};

// Selection criteria
Configurable<double> v0cospacut{"hypcospa", 0.95, "V0 CosPA"};
Expand Down Expand Up @@ -199,6 +221,13 @@ struct hyperRecoTask {
ConfigurableAxis zVtxBins{"zVtxBins", {100, -20.f, 20.f}, "Binning for n sigma"};
ConfigurableAxis centBins{"centBins", {100, 0.f, 100.f}, "Binning for centrality"};

// histogram axes for EvtLossMC
ConfigurableAxis binsImpactPar{"binsImpactPar", {80, 0, 16}, "Binning of the impact parameter axis"};
ConfigurableAxis binsCent{"binsCent", {10, 0.0, 100.0}, "Binning of the centrality axis"};
ConfigurableAxis binsPt{"binsPt", {20, 0, 10}, "Binning of the pt"};
ConfigurableAxis binsFT0CMult{"binsFT0CMult", {500, 0.0f, +500.0f}, "Binning of the FT0C multiplicity"};
ConfigurableAxis binsMult{"binsMult", {500, 0.0f, +500.0f}, ""};

// std vector of candidates
std::vector<hyperCandidate> hyperCandidates;
// vector to keep track of MC mothers already filled
Expand All @@ -219,7 +248,6 @@ struct hyperRecoTask {

void init(InitContext const&)
{

zorroSummary.setObject(zorro.getZorroSummary());

mRunNumber = 0;
Expand Down Expand Up @@ -250,6 +278,11 @@ struct hyperRecoTask {
const AxisSpec nSigma3HeAxis{nSigmaBins, "n_{#sigma}({}^{3}He)"};
const AxisSpec zVtxAxis{zVtxBins, "z_{vtx} (cm)"};
const AxisSpec centAxis{centBins, "Centrality"};
const AxisSpec impactParamAxis{binsImpactPar, "Impact Parameter (b)"};
const AxisSpec centFT0CAxis{binsCent, "Centrality (FT0C %)"};
const AxisSpec binsFT0CMultAxis{binsFT0CMult, "FT0C multiplicity"};
const AxisSpec ptAxis{binsPt, "#it{p}_{T} (GeV/#it{c})"};
const AxisSpec multAxis = {binsMult, "Multiplicity #eta <0.5"};

hNsigma3HeSel = qaRegistry.add<TH2>("hNsigma3HeSel", "; p_{TPC}/z (GeV/#it{c}); n_{#sigma} ({}^{3}He)", HistType::kTH2F, {rigidityAxis, nSigma3HeAxis});
hDeDx3HeSel = qaRegistry.add<TH2>("hDeDx3HeSel", ";p_{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {rigidityAxis, dedxAxis});
Expand All @@ -259,11 +292,12 @@ struct hyperRecoTask {
hH4LMassBefSel = qaRegistry.add<TH1>("hH4LMassBefSel", ";M (GeV/#it{c}^{2}); ", HistType::kTH1D, {{60, 3.76, 3.84}});
hH4LMassTracked = qaRegistry.add<TH1>("hH4LMassTracked", ";M (GeV/#it{c}^{2}); ", HistType::kTH1D, {{60, 3.76, 3.84}});

hEvents = qaRegistry.add<TH1>("hEvents", ";Events; ", HistType::kTH1D, {{4, -0.5, 3.5}});
hEvents = qaRegistry.add<TH1>("hEvents", ";Events; ", HistType::kTH1D, {{5, -0.5, 4.5}});
hEvents->GetXaxis()->SetBinLabel(1, "All");
hEvents->GetXaxis()->SetBinLabel(2, "sel8");
hEvents->GetXaxis()->SetBinLabel(3, "kNoSameBunchPileup");
hEvents->GetXaxis()->SetBinLabel(4, "kIsGoodZvtxFT0vsPV");
hEvents->GetXaxis()->SetBinLabel(3, "z_{vtx}");
hEvents->GetXaxis()->SetBinLabel(4, "kNoSameBunchPileup");
hEvents->GetXaxis()->SetBinLabel(5, "kIsGoodZvtxFT0vsPV");

hEventsZorro = qaRegistry.add<TH1>("hEventsZorro", ";Events; ", HistType::kTH1D, {{2, -0.5, 1.5}});
hEventsZorro->GetXaxis()->SetBinLabel(1, "Zorro before evsel");
Expand All @@ -284,6 +318,29 @@ struct hyperRecoTask {
hCentFT0A = qaRegistry.add<TH1>("hCentFT0A", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}});
hCentFT0C = qaRegistry.add<TH1>("hCentFT0C", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}});
hCentFT0M = qaRegistry.add<TH1>("hCentFT0M", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}});

if (doprocessEventLossMC) {
hEvtMC = qaRegistry.add<TH1>("QAEvent/hEvtMC", ";; ", HistType::kTH1D, {{3, -0.5, 2.5}});
hEvtMC->GetXaxis()->SetBinLabel(1, "All gen evts");
hEvtMC->GetXaxis()->SetBinLabel(2, "Gen evts with al least one reconstructed");
hEvtMC->GetXaxis()->SetBinLabel(3, "Gen evts with no reconstructed collisions");
// Infomation for all generated collisions collisions
hImpactParamGen = qaRegistry.add<TH1>("QAEvent/McColAll/hImpactParamGen", "Impact parameter of generated MC events; Impact Parameter (b); Counts", HistType::kTH1D, {impactParamAxis});
// Infomation for generated collisions collisions with at least one rec. event
hImpactParamReco = qaRegistry.add<TH1>("QAEvent/McColAll/hImpactParamReco", "Impact parameter of generated MC events with at least one rec. evt; Impact Parameter (b); Counts", HistType::kTH1D, {impactParamAxis});
hRecoCentralityColvsMultiplicityRecoEta05 = qaRegistry.add<TH2>("QAEvent/McColAll/hRecoCentralityColvsMultiplicityRecoEta05", "Correlation between FT0C centrality and charged particle multiplicity in generated MC events with at least one rec. evt; Multiplicity #eta <0.5; Counts", HistType::kTH2D, {centFT0CAxis, multAxis});
hRecoCentralityColvsImpactParamReco = qaRegistry.add<TH2>("QAEvent/McColAll/hRecoCentralityColvsImpactParamReco", "Correlation between FT0C centrality and impact parameter in generated MC events with at least one rec. evt; Impact Parameter (b); Counts", HistType::kTH2D, {centFT0CAxis, impactParamAxis});
// Information of generated 3HL in generated events
hGen3HLBeforeEvtSel = qaRegistry.add<TH1>("QAEvent/McCol3HL/hGen3HLBeforeEvtSel", "3HL generated #it{p}_{T} distribution in all gen evt;#it{p}_{T} (GeV/#it{c}); Counts", HistType::kTH1D, {ptAxis});
hGen3HLvsImpactParameterBeforeEvtSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsImpactParameterBeforeEvtSel", "Correlation 3HL generated #it{p}_{T} and impact parameter in all gen evt;#it{p}_{T} (GeV/#it{c}); Impact parameter (b)", HistType::kTH2D, {ptAxis, impactParamAxis});
hGen3HLvsMultiplicityGenEta05BeforeEvtSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityGenEta05BeforeEvtSel", "Correlation 3HL generated #it{p}_{T} and charged particle multiplicity in all gen evt;#it{p}_{T} (GeV/#it{c}); Multiplicity #eta <0.5", HistType::kTH2D, {ptAxis, multAxis});
hGen3HLvsMultiplicityFT0CBeforeEvtSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityFT0CBeforeEvtSel", "Correlation 3HL generated #it{p}_{T} and FT0C multiplicity in all gen evt;#it{p}_{T} (GeV/#it{c}); FT0C Multiplicity", HistType::kTH2D, {ptAxis, binsFT0CMultAxis});
// Information of generated 3HL in generated events with at least one rec. event
hGen3HLAfterSel = qaRegistry.add<TH1>("QAEvent/McCol3HL/hGen3HLAfterSel", "3HL generated #it{p}_{T} distribution in gen. evts with at least one rec. evt; #it{p}_{T} (GeV/#it{c}); Counts", HistType::kTH1D, {ptAxis});
hGen3HLvsImpactParameterAfterSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsImpactParameterAfterSel", "Correlation 3HL generated #it{p}_{T} and impact parameter in gen. evts with at least one rec. evt;#it{p}_{T} (GeV/#it{c}); Impact parameter (b)", HistType::kTH2D, {ptAxis, impactParamAxis});
hGen3HLvsMultiplicityGenEta05AfterSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityGenEta05AfterSel", "Correlation 3HL generated #it{p}_{T} and charged particle multiplicity in gen. evts with at least one rec. evt;#it{p}_{T} (GeV/#it{c}); Multiplicity #eta <0.5", HistType::kTH2D, {ptAxis, multAxis});
hGen3HLvsMultiplicityFT0CAfterSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityFT0CAfterSel", "Correlation 3HL generated #it{p}_{T} and FT0C multiplicity in gen. evts with at least one rec;#it{p}_{T} (GeV/#it{c}); FT0C Multiplicity", HistType::kTH2D, {ptAxis, binsFT0CMultAxis});
}
}

void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
Expand Down Expand Up @@ -366,12 +423,17 @@ struct hyperRecoTask {
}
}

if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > 10) {
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
continue;
}

hEvents->Fill(1.);

if (std::abs(collision.posZ()) > 10) {
hEvents->Fill(2.);
continue;
}

if (zorroSelected) {
hEventsZorro->Fill(1.);
}
Expand All @@ -380,14 +442,13 @@ struct hyperRecoTask {
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
continue;
}
hEvents->Fill(2.);
hEvents->Fill(3.);
}

if (cfgEvSelkIsGoodZvtxFT0vsPV) {
if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) {
continue;
}
hEvents->Fill(3.);
hEvents->Fill(4.);
}

goodCollision[collision.globalIndex()] = true;
Expand All @@ -408,23 +469,27 @@ struct hyperRecoTask {
if (collision.has_mcCollision()) {
recoCollisionIds[collision.mcCollisionId()] = collision.globalIndex();
}
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > 10)
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder))
continue;

hEvents->Fill(1.);

if (std::abs(collision.posZ()) > 10) {
hEvents->Fill(2.);
continue;
}

if (cfgEvSelkNoSameBunchPileup) {
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
continue;
}
hEvents->Fill(2.);
hEvents->Fill(3.);
}

if (cfgEvSelkIsGoodZvtxFT0vsPV) {
if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) {
continue;
}
hEvents->Fill(3.);
hEvents->Fill(4.);
}

if (collision.has_mcCollision()) {
Expand All @@ -437,7 +502,6 @@ struct hyperRecoTask {
hCentFT0M->Fill(collision.centFT0M());
}
}

template <class Ttrack, class Tcolls>
void fillHyperCand(Ttrack& heTrack, Ttrack& piTrack, CollBracket collBracket, const Tcolls& collisions, hyperCandidate& hypCand)
{
Expand Down Expand Up @@ -922,6 +986,111 @@ struct hyperRecoTask {
processMC(collisions, mcCollisions, V0s, tracks, ambiTracks, bcs, trackLabelsMC, particlesMC);
}
PROCESS_SWITCH(hyperRecoTask, processMCTracked, "MC analysis with tracked V0s", false);

template <typename CollType>
bool passEvtSel(const CollType& collision)
{
if (!collision.sel8())
return false;

if ((std::abs(collision.posZ())) > 10)
return false;

if (cfgEvSelkNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup))
return false;

if (cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV))
return false;

return true;
}

void processEventLossMC(McCollisionMults::iterator const& mcCollision, soa::SmallGroups<EventCandidatesMC> const& collisions, aod::McParticles const& GenParticles)
{
if (std::abs(mcCollision.posZ()) > 10) {
return;
}

//////////// Event loss estimation via impact parameter and multiplicity by MCFT0C

// Fill all generated events
hEvtMC->Fill(0);
hImpactParamGen->Fill(mcCollision.impactParameter());

// Fill generated events with no reconstructed collisions
if (collisions.size() == 0) {
hEvtMC->Fill(1);
}

// Define the generated events with at least one reconstructed event
bool atLeastOneRecoEvt = false;
auto centralityFT0C = -999.;

for (auto const& col : collisions) {
if (!passEvtSel(col)) {
continue;
}
centralityFT0C = col.centFT0C();
atLeastOneRecoEvt = true;
}

if (atLeastOneRecoEvt) {
hEvtMC->Fill(2);
hImpactParamReco->Fill(mcCollision.impactParameter());
hRecoCentralityColvsMultiplicityRecoEta05->Fill(centralityFT0C, mcCollision.multMCNParticlesEta05());
hRecoCentralityColvsImpactParamReco->Fill(centralityFT0C, mcCollision.impactParameter());
}
// Construct the H3L 4-vector based on the generated daugthers identification by PDG
ROOT::Math::PxPyPzMVector daugh1, daugh2, mother;

for (const auto& genParticle : GenParticles) {
if (std::abs(genParticle.y()) > 1)
continue;
if (std::abs(genParticle.pdgCode()) != hyperPdg)
continue;

auto daughters = genParticle.daughters_as<aod::McParticles>();

bool dauHe3 = false, dauPiMinus = false, dauAntiHe3 = false, dauPiPos = true;

for (const auto& daughter : daughters) {
if (daughter.pdgCode() == heDauPdg) {
dauHe3 = true;
daugh1 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), he3Mass);
} else if (daughter.pdgCode() == -piDauPdg) {
dauPiMinus = true;
daugh2 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), piMass);
}
if (daughter.pdgCode() == -heDauPdg) {
dauAntiHe3 = true;
daugh1 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), he3Mass);
} else if (daughter.pdgCode() == piDauPdg) {
dauPiPos = true;
daugh2 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), piMass);
}
}
// Check pairs to avoid wrong charge associations
if (!((dauHe3 && dauPiMinus) || !(dauAntiHe3 && dauPiPos)))
continue;

mother = daugh1 + daugh2;

// Fill informations for generated 3HL in all generated events
hGen3HLBeforeEvtSel->Fill(mother.pt());
hGen3HLvsImpactParameterBeforeEvtSel->Fill(mother.pt(), mcCollision.impactParameter());
hGen3HLvsMultiplicityGenEta05BeforeEvtSel->Fill(mother.pt(), mcCollision.multMCNParticlesEta05());
hGen3HLvsMultiplicityFT0CBeforeEvtSel->Fill(mother.pt(), mcCollision.multMCFT0C());

// Fill informations for generated 3HL in generated events with at least one reconstructed event
if (atLeastOneRecoEvt) {
hGen3HLAfterSel->Fill(mother.pt());
hGen3HLvsImpactParameterAfterSel->Fill(mother.pt(), mcCollision.impactParameter());
hGen3HLvsMultiplicityGenEta05AfterSel->Fill(mother.pt(), mcCollision.multMCNParticlesEta05());
hGen3HLvsMultiplicityFT0CAfterSel->Fill(mother.pt(), mcCollision.multMCFT0C());
}
}
}
PROCESS_SWITCH(hyperRecoTask, processEventLossMC, "Event loss analysis", false);
};

WorkflowSpec
Expand Down
Loading