diff --git a/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx b/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx index 096b427617e..e200e66ae57 100644 --- a/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx @@ -40,6 +40,8 @@ #include "MathUtils/BetheBlochAleph.h" #include "ReconstructionDataFormats/Track.h" +#include "Math/Vector4D.h" + #include #include #include @@ -49,6 +51,8 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using std::array; + using CollBracket = o2::math_utils::Bracket; using TracksFull = soa::Join; using CollisionsFull = soa::Join; @@ -56,6 +60,9 @@ using CollisionsFullMC = soa::Join; +using McCollisionMults = soa::Join; +using EventCandidatesMC = soa::Join; + namespace { constexpr double betheBlochDefault[1][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}}; @@ -77,6 +84,20 @@ std::shared_ptr hH4LMassTracked; std::shared_ptr hDecayChannel; std::shared_ptr hIsMatterGen; std::shared_ptr hIsMatterGenTwoBody; +std::shared_ptr hEvtMC; +std::shared_ptr hImpactParamGen; +std::shared_ptr hImpactParamReco; +std::shared_ptr hGen3HLBeforeEvtSel; +std::shared_ptr hGen3HLAfterSel; +std::shared_ptr hRecoCentralityColvsMultiplicityRecoEta05; +std::shared_ptr hRecoCentralityColvsImpactParamReco; +std::shared_ptr hGen3HLvsImpactParameterBeforeEvtSel; +std::shared_ptr hGen3HLvsImpactParameterAfterSel; +std::shared_ptr hGen3HLvsMultiplicityGenEta05BeforeEvtSel; +std::shared_ptr hGen3HLvsMultiplicityGenEta05AfterSel; +std::shared_ptr hGen3HLvsMultiplicityFT0CBeforeEvtSel; +std::shared_ptr hGen3HLvsMultiplicityFT0CAfterSel; + } // namespace struct hyperCandidate { @@ -147,6 +168,7 @@ struct hyperRecoTask { // PDG codes Configurable hyperPdg{"hyperPDG", 1010010030, "PDG code of the hyper-mother (could be 3LamH or 4LamH)"}; Configurable heDauPdg{"heDauPDG", 1000020030, "PDG code of the helium (could be 3He or 4He)"}; + Configurable piDauPdg{"piDauPdg", 211, "PDG code of pion"}; // Selection criteria Configurable v0cospacut{"hypcospa", 0.95, "V0 CosPA"}; @@ -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 hyperCandidates; // vector to keep track of MC mothers already filled @@ -219,7 +248,6 @@ struct hyperRecoTask { void init(InitContext const&) { - zorroSummary.setObject(zorro.getZorroSummary()); mRunNumber = 0; @@ -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("hNsigma3HeSel", "; p_{TPC}/z (GeV/#it{c}); n_{#sigma} ({}^{3}He)", HistType::kTH2F, {rigidityAxis, nSigma3HeAxis}); hDeDx3HeSel = qaRegistry.add("hDeDx3HeSel", ";p_{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {rigidityAxis, dedxAxis}); @@ -259,11 +292,12 @@ struct hyperRecoTask { hH4LMassBefSel = qaRegistry.add("hH4LMassBefSel", ";M (GeV/#it{c}^{2}); ", HistType::kTH1D, {{60, 3.76, 3.84}}); hH4LMassTracked = qaRegistry.add("hH4LMassTracked", ";M (GeV/#it{c}^{2}); ", HistType::kTH1D, {{60, 3.76, 3.84}}); - hEvents = qaRegistry.add("hEvents", ";Events; ", HistType::kTH1D, {{4, -0.5, 3.5}}); + hEvents = qaRegistry.add("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("hEventsZorro", ";Events; ", HistType::kTH1D, {{2, -0.5, 1.5}}); hEventsZorro->GetXaxis()->SetBinLabel(1, "Zorro before evsel"); @@ -284,6 +318,29 @@ struct hyperRecoTask { hCentFT0A = qaRegistry.add("hCentFT0A", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); hCentFT0C = qaRegistry.add("hCentFT0C", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); hCentFT0M = qaRegistry.add("hCentFT0M", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); + + if (doprocessEventLossMC) { + hEvtMC = qaRegistry.add("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("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("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("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("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("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("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("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("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("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("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("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("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) @@ -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.); } @@ -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; @@ -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()) { @@ -437,7 +502,6 @@ struct hyperRecoTask { hCentFT0M->Fill(collision.centFT0M()); } } - template void fillHyperCand(Ttrack& heTrack, Ttrack& piTrack, CollBracket collBracket, const Tcolls& collisions, hyperCandidate& hypCand) { @@ -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 + 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 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(); + + 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