From 83167c069a30f30e8ccc080cd9072e14250cc93d Mon Sep 17 00:00:00 2001 From: Maria Paula Martins Palhares Date: Wed, 2 Sep 2026 15:20:04 +0200 Subject: [PATCH 1/2] Added SigEvtLoss process function + nTPCCLus for candidates selection studies + minor changes --- PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx | 503 +++++++++++++++------ 1 file changed, 371 insertions(+), 132 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx b/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx index a9a823d364c..5e3eaa7c85c 100644 --- a/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx @@ -19,6 +19,7 @@ #include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -41,8 +42,11 @@ #include #include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include #include #include +#include #include #include @@ -59,16 +63,20 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using std::array; + using TracksFull = soa::Join; using TracksFullMC = soa::Join; -using CollisionsFull = soa::Join; -using CollisionsFullMC = soa::Join; +using CollisionsFull = soa::Join; +using CollisionsFullMC = soa::Join; using CollisionCandidates = o2::soa::Join; using CollisionCandidatesMC = o2::soa::Join; using TrackCandidates = o2::soa::Join; using TrackCandidatesMC = o2::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}}; @@ -83,6 +91,7 @@ std::shared_ptr hCentFT0M; std::shared_ptr hCentFV0A; std::shared_ptr hNsigma3HSel; std::shared_ptr hNsigma3HSelTOF; +std::shared_ptr h2Nsigma3HSelTPCTOF; std::shared_ptr hdEdx3HSel; std::shared_ptr hdEdx3HPosTrack; std::shared_ptr hdEdx3HNegTrack; @@ -92,6 +101,8 @@ std::shared_ptr h3HSignalPtTOF; std::shared_ptr hIsMatterGen; std::shared_ptr hLnnCandLoss; std::shared_ptr hDecayChannel; +// QA centrality scaling +std::shared_ptr h2GlobalTrackMultVsCent; // QA histograms before track selections std::shared_ptr h2FT0CnClusTPCtoTrBfSel; std::shared_ptr h2FT0CnClusTPCtoPiBfSel; @@ -106,20 +117,42 @@ std::shared_ptr h2FT0CptPiBfSelItsTpc; // QA MC tracks from the lnn decay std::shared_ptr hPtGeneratedTrStr; std::shared_ptr hPtGeneratedPiStr; +std::shared_ptr hPtGeneratedLnnStr; +std::shared_ptr hPtTpcTrStr; std::shared_ptr hPtItsTpcTrStr; +std::shared_ptr hPtTrkTpcTrStr; std::shared_ptr hPtTrkItsTpcTrStr; std::shared_ptr hPtItsTpcTofTrStr; std::shared_ptr hPtTrkItsTpcTofTrStr; -std::shared_ptr hPtItsTpcPiStr; -std::shared_ptr hPtTrkItsTpcPiStr; + // QA reco generated candidate and daugher particles from secondary vertex std::shared_ptr h2FT0CPtGenColRecCandMC; std::shared_ptr h2FT0CPtGenColRecTrStrMC; -std::shared_ptr h2FT0CPtGenColRecPiStrMC; // QA signal generated candidate and daugher particles from secondary vertex std::shared_ptr h2FT0CPtGenColGenCandMC; std::shared_ptr h2FT0CPtGenColGenTrStrMC; -std::shared_ptr h2FT0CPtGenColGenPiStrMC; +// Event and signal loss +std::shared_ptr hEvtMC; +std::shared_ptr hImpactParamGen; +std::shared_ptr hImpactParamGenOneReco; +std::shared_ptr hImpactParamReco; +std::shared_ptr hGenLnnBeforeEvtSel; +std::shared_ptr hGenLnnAfterSel; +std::shared_ptr hGenOneRecoCentrality; +std::shared_ptr hRecoCentrality; +std::shared_ptr hGenEventsNchEta08; +std::shared_ptr hGenCentralityColvsMultiplicityGenEta08; +std::shared_ptr hGenCentralityColvsImpactParamGen; +std::shared_ptr hGenCentralityColvsFT0Cmultiplicity; +std::shared_ptr hRecoCentralityColvsMultiplicityRecoEta08; +std::shared_ptr hRecoCentralityColvsImpactParamReco; +std::shared_ptr hRecoCentralityColvsFT0Cmultiplicity; +std::shared_ptr hGenLnnvsImpactParameterBeforeEvtSel; +std::shared_ptr hGenLnnvsImpactParameterAfterSel; +std::shared_ptr hGenLnnvsMultiplicityGenEta08BeforeEvtSel; +std::shared_ptr hGenLnnvsMultiplicityGenEta08AfterSel; +std::shared_ptr hGenLnnvsMultiplicityFT0CBeforeEvtSel; +std::shared_ptr hGenLnnvsMultiplicityFT0CAfterSel; float alphaAP(std::array const& momB, std::array const& momC) { @@ -149,10 +182,6 @@ struct LnnCandidate { float genPhi3H() const { return std::atan2(gMom3H[1], gMom3H[0]); } float genEta3H() const { return std::asinh(gMom3H[2] / genPt3H()); } - float genPtPi() const { return std::hypot(gMomPi[0], gMomPi[1]); } - float genPhiPi() const { return std::atan2(gMomPi[1], gMomPi[0]); } - float genEtaPi() const { return std::asinh(gMomPi[2] / genPtPi()); } - int posTrackID; int negTrackID; float dcaV0dau = -10; @@ -172,7 +201,6 @@ struct LnnCandidate { std::array decVtx; std::array gMom; std::array gMom3H; - std::array gMomPi; std::array gDecVtx; uint16_t tpcSignal3H = 0u; uint16_t tpcSignalPi = 0u; @@ -191,7 +219,7 @@ struct LnnCandidate { uint8_t flags = 0u; // flags for dughter particles }; -struct lnnRecoTask { +struct LnnRecoTask { Produces outputDataTable; Produces outputMCTable; @@ -246,7 +274,7 @@ struct lnnRecoTask { Configurable piDauPdg{"piDauPdg", 211, "PDG Pi"}; // PDG Pi Configurable lnnPdg{"lnnPdg", 1010000030, "PDG Lnn"}; // PDG Lnn - // Histogram configuration QA lnn-task + // histograms axes ConfigurableAxis rigidityBins{"rigidityBins", {200, -10.f, 10.f}, "Binning for rigidity"}; ConfigurableAxis dEdxBins{"dEdxBins", {5000, 0.f, 1000.f}, "Binning for dE/dx"}; ConfigurableAxis nSigmaBins{"nSigmaBins", {200, -5.f, 5.f}, "Binning for n#sigma_{TPC}"}; @@ -263,9 +291,17 @@ struct lnnRecoTask { ConfigurableAxis tpcNClsCrossedRowsBins{"tpcNClsCrossedRowsBins", {260, 30, 165}, "Binning for TPCNClsCrossedRows"}; ConfigurableAxis tpcChi2NClusBins{"tpcChi2NClusBins", {20, 0.5, 10}, "Binning for chi2NClusTPC"}; ConfigurableAxis itsChi2NClusBins{"itsChi2NClusBins", {72, 0, 36}, "Binning for chi2NClusTPC"}; - ConfigurableAxis candPtBins{"candPtBins", {160, 0, 8}, "Binning for lnn cand pt"}; + ConfigurableAxis candPtBins{"candPtBins", {32, -8, 8}, "Binning for lnn cand pt"}; ConfigurableAxis candEtaBins{"candEtaBins", {160, -0.8, 0.8}, "Binning for eta"}; + // 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}, ""}; + ConfigurableAxis binsglobalTrackMult{"binsglobalTrackMult", {100, 0.f, 1000.f}, "Binning for global track multiplicity"}; + // std vector of candidates std::vector lnnCandidates; // vector to keep track of MC mothers already filled @@ -277,11 +313,15 @@ struct lnnRecoTask { // vector to armazenade h3Track Preslice perCollision = o2::aod::v0::collisionId; + + Preslice mcParticlesPerMcCollision = aod::mcparticle::mcCollisionId; + Preslice tracksPerCollision = aod::track::collisionId; + HistogramRegistry qaRegistry{"QA", {}, OutputObjHandlingPolicy::AnalysisObject}; - int mRunNumber; - float d_bz; - std::array mBBparams3H; + int mRunNumber{-1}; + float dBz{0.f}; + std::array mBBparams3H{}; static constexpr float KallEvents = 0.; static constexpr float KevAfterSel8 = 1.; @@ -292,7 +332,7 @@ struct lnnRecoTask { void init(InitContext const&) { mRunNumber = 0; - d_bz = 0; + dBz = 0; ccdb->setURL(ccdburl); ccdb->setCaching(true); @@ -311,7 +351,8 @@ struct lnnRecoTask { const AxisSpec rigidityAxis{rigidityBins, "#it{p}^{TPC}/#it{z}"}; const AxisSpec dEdxAxis{dEdxBins, "d#it{E}/d#it{x}"}; - const AxisSpec nSigma3HAxis{nSigmaBins, "n_{#sigma}({}^{3}H)"}; + const AxisSpec nSigma3HtpcAxis{nSigmaBins, "n_{#sigma TPC}({}^{3}H)"}; + const AxisSpec nSigma3HtofAxis{nSigmaBins, "n_{#sigma TOF}({}^{3}H)"}; const AxisSpec zVtxAxis{zVtxBins, "z_{vtx} (cm)"}; const AxisSpec centAxis{centBins, "FT0C (%)"}; const AxisSpec mTOFAxis{mTOFBins, "#frac{m^{2}}{z^{2}}"}; @@ -326,8 +367,17 @@ struct lnnRecoTask { const AxisSpec candPtAxis(candPtBins, "#it{p}_{T} (Gev/#it{c})"); const AxisSpec candEtaAxis(candEtaBins, "#eta"); - hNsigma3HSel = qaRegistry.add("PID/hNsigma3HSel", ";#it{p}^{TPC}/z (GeV/#it{c}); n_{#sigma} ({}^{3}H)", HistType::kTH2F, {rigidityAxis, nSigma3HAxis}); - hNsigma3HSelTOF = qaRegistry.add("PID/hNsigma3HSelTOF", ";#it{p}_{T} (GeV/#it{c}); n_{#sigma} ({}^{3}H)", HistType::kTH2F, {tPtAxis, nSigma3HAxis}); + 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"}; + const AxisSpec globalTrackMultAxis = {binsglobalTrackMult, "N_{global track mult.}"}; + + + hNsigma3HSel = qaRegistry.add("PID/hNsigma3HSel", ";#it{p}^{TPC}/z (GeV/#it{c}); n_{#sigma} ({}^{3}H)", HistType::kTH2F, {rigidityAxis, nSigma3HtpcAxis}); + hNsigma3HSelTOF = qaRegistry.add("PID/hNsigma3HSelTOF", ";#it{p}_{T} (GeV/#it{c}); n_{#sigma} ({}^{3}H)", HistType::kTH2F, {tPtAxis, nSigma3HtpcAxis}); + h2Nsigma3HSelTPCTOF = qaRegistry.add("PID/h2Nsigma3HSelTPCTOF", ";n_{#sigma TPC} ({}^{3}H); n_{#sigma TOF} ({}^{3}H)", HistType::kTH2F, {nSigma3HtpcAxis, nSigma3HtofAxis}); hdEdx3HSel = qaRegistry.add("hdEdx3HSel", ";#it{p}^{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {rigidityAxis, dEdxAxis}); hdEdx3HPosTrack = qaRegistry.add("PID/hdEdx3HPosTrack", "; #it{p}^{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {tPosRigidityAxis, dEdxAxis}); hdEdx3HNegTrack = qaRegistry.add("PID/hdEdx3HNegTrack", "; #it{p}^{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {tPNegRigidityAxis, dEdxAxis}); @@ -335,27 +385,16 @@ struct lnnRecoTask { h3HMassPtTOF = qaRegistry.add("PID/hTrMassPtTOF", "; #it{p}_{T} ({}^{3}H) (GeV/#it{c}); #frac{m^{2}}{z^{2}} (GeV^{2}/#it{c}^{4})", HistType::kTH2F, {tPtAxis, mTOFAxis}); h3HSignalPtTOF = qaRegistry.add("PID/h3HSignalPtTOF", "; #it{p}_{T}({}^{3}H) (GeV/#it{c}); #beta_{TOF}", HistType::kTH2F, {tPtAxis, betaAxis}); hEvents = qaRegistry.add("hEvents", ";Events; ", HistType::kTH1D, {{4, -0.5, 3.5}}); - hLnnCandLoss = qaRegistry.add("CandCounts/hLnnCandLoss", ";CandLoss; ", HistType::kTH1D, {{7, -0.5, 6.5}}); - // QA its-tpc and its-tpc-tof tracks before selection - h2FT0CnClusTPCtoTrBfSel = qaRegistry.add("QATracks/h2FT0CnClusTPCtoTrBfSel", ";FT0C (%);N_{clus}^{TPC}", HistType::kTH2F, {centAxis, tpcNClusAxis}); - h2FT0CnClusTPCtoPiBfSel = qaRegistry.add("QATracks/h2FT0CnClusTPCtoPiBfSel", ";FT0C (%);N_{clus}^{TPC}", HistType::kTH2F, {centAxis, tpcNClusAxis}); - h2FT0Cchi2NClTPCtoTrBfSel = qaRegistry.add("QATracks/h2FT0Cchi2NClTPCtoTrBfSel", ";FT0C (%);{#Chi}^{2}/N_{clus}^{TPC} ", HistType::kTH2F, {centAxis, tpcChi2NClusAxis}); - h2FT0Cchi2NClITStoTrBfSel = qaRegistry.add("QATracks/h2FT0Cchi2NClITStoTrBfSel", ";FT0C (%);{#Chi}^{2}/N_{clus}^{ITS}", HistType::kTH2F, {centAxis, itsChi2NClusAxis}); - h2FT0CnTPCNClsCrossedRows3HBfSel = qaRegistry.add("QATracks/h2FT0CnTPCNClsCrossedRows3H", ";FT0C (%);N_{TPC} crossed rows", HistType::kTH2F, {centAxis, tpcNClsCrossedRowsAxis}); - h2FT0CnTPCNClsCrossedRowsPiBfSel = qaRegistry.add("QATracks/h2FT0CnTPCNClsCrossedRowsPi", ";FT0C (%);N_{TPC} crossed rows", HistType::kTH2F, {centAxis, tpcNClsCrossedRowsAxis}); - h2FT0CptTrBfSelItsTpc = qaRegistry.add("QATracks/h2FT0CptTrBfSelItsTpc", ";FT0C (%);#it{p}_{T}", HistType::kTH2F, {centAxis, candPtBins}); - h2FT0CptTrBfSelItsTpcTof = qaRegistry.add("QATracks/h2FT0CptTrBfSelItsTpcTof", ";FT0C (%);#it{p}_{T}", HistType::kTH2F, {centAxis, candPtBins}); - h2FT0CptPiBfSelItsTpc = qaRegistry.add("QATracks/h2FT0CptPiBfSelItsTpc", ";FT0C (%);#it{p}_{T}", HistType::kTH2F, {centAxis, candPtBins}); - // QA its-tpc, its-tpc-tof, and generated MC - hPtItsTpcTrStr = qaRegistry.add("MC/McTracks/hPtItsTpcTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtTrkItsTpcTrStr = qaRegistry.add("MC/McTracks/hPtTrkItsTpcTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtItsTpcTofTrStr = qaRegistry.add("MC/McTracks/hPtItsTpcTofTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtTrkItsTpcTofTrStr = qaRegistry.add("MC/McTracks/hPtTrkItsTpcTofTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtItsTpcPiStr = qaRegistry.add("MC/McTracks/hPtItsTpcPiStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtTrkItsTpcPiStr = qaRegistry.add("MC/McTracks/hPtTrkItsTpcPiStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtGeneratedTrStr = qaRegistry.add("MC/McGen/hPtGeneratedTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); - hPtGeneratedPiStr = qaRegistry.add("MC/McGen/hPtGeneratedPiStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {160, 0, 8}}); + h2GlobalTrackMultVsCent = qaRegistry.add("QATracks/h2GlobalTrackMultVsCent", ";FT0C (%); N_{global track mult.}", HistType::kTH2F, {centAxis, globalTrackMultAxis}); + + hZvtx = qaRegistry.add("hZvtx", ";z_{vtx} (cm); ", HistType::kTH1D, {{100, -20, 20}}); + 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}}); + hCentFV0A = qaRegistry.add("hCentFV0A", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); + + hLnnCandLoss = qaRegistry.add("CandCounts/hLnnCandLoss", ";CandLoss; ", HistType::kTH1D, {{7, -0.5, 6.5}}); hEvents->GetXaxis()->SetBinLabel(1, "All"); hEvents->GetXaxis()->SetBinLabel(2, "sel8"); hEvents->GetXaxis()->SetBinLabel(3, "z_{vtx}"); @@ -369,6 +408,27 @@ struct lnnRecoTask { hLnnCandLoss->GetXaxis()->SetBinLabel(5, "!isLnnMass"); hLnnCandLoss->GetXaxis()->SetBinLabel(6, "DCA #it{V}_{0} daughter"); hLnnCandLoss->GetXaxis()->SetBinLabel(7, "cosPA"); + + // QA its-tpc and its-tpc-tof tracks before selection + h2FT0CnClusTPCtoTrBfSel = qaRegistry.add("QATracks/h2FT0CnClusTPCtoTrBfSel", ";FT0C (%);N_{clus}^{TPC}", HistType::kTH2F, {centAxis, tpcNClusAxis}); + h2FT0CnClusTPCtoPiBfSel = qaRegistry.add("QATracks/h2FT0CnClusTPCtoPiBfSel", ";FT0C (%);N_{clus}^{TPC}", HistType::kTH2F, {centAxis, tpcNClusAxis}); + h2FT0Cchi2NClTPCtoTrBfSel = qaRegistry.add("QATracks/h2FT0Cchi2NClTPCtoTrBfSel", ";FT0C (%);{#Chi}^{2}/N_{clus}^{TPC} ", HistType::kTH2F, {centAxis, tpcChi2NClusAxis}); + h2FT0Cchi2NClITStoTrBfSel = qaRegistry.add("QATracks/h2FT0Cchi2NClITStoTrBfSel", ";FT0C (%);{#Chi}^{2}/N_{clus}^{ITS}", HistType::kTH2F, {centAxis, itsChi2NClusAxis}); + h2FT0CnTPCNClsCrossedRows3HBfSel = qaRegistry.add("QATracks/h2FT0CnTPCNClsCrossedRows3H", ";FT0C (%);N_{TPC} crossed rows", HistType::kTH2F, {centAxis, tpcNClsCrossedRowsAxis}); + h2FT0CnTPCNClsCrossedRowsPiBfSel = qaRegistry.add("QATracks/h2FT0CnTPCNClsCrossedRowsPi", ";FT0C (%);N_{TPC} crossed rows", HistType::kTH2F, {centAxis, tpcNClsCrossedRowsAxis}); + h2FT0CptTrBfSelItsTpc = qaRegistry.add("QATracks/h2FT0CptTrBfSelItsTpc", ";FT0C (%);#it{p}_{T}", HistType::kTH2F, {centAxis, candPtAxis}); + h2FT0CptTrBfSelItsTpcTof = qaRegistry.add("QATracks/h2FT0CptTrBfSelItsTpcTof", ";FT0C (%);#it{p}_{T}", HistType::kTH2F, {centAxis, candPtAxis}); + h2FT0CptPiBfSelItsTpc = qaRegistry.add("QATracks/h2FT0CptPiBfSelItsTpc", ";FT0C (%);#it{p}_{T}", HistType::kTH2F, {centAxis, candPtAxis}); + // QA its-tpc, its-tpc-tof, and generated MC + hPtTpcTrStr = qaRegistry.add("MC/McTracks/hPtTpcTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtTrkTpcTrStr = qaRegistry.add("MC/McTracks/hPtTrkTpcTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtItsTpcTrStr = qaRegistry.add("MC/McTracks/hPtItsTpcTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtTrkItsTpcTrStr = qaRegistry.add("MC/McTracks/hPtTrkItsTpcTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtItsTpcTofTrStr = qaRegistry.add("MC/McTracks/hPtItsTpcTofTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtTrkItsTpcTofTrStr = qaRegistry.add("MC/McTracks/hPtTrkItsTpcTofTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtGeneratedTrStr = qaRegistry.add("MC/McGen/hPtGeneratedTrStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + hPtGeneratedLnnStr = qaRegistry.add("MC/McGen/hPtGeneratedLnnStr", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{10, 0, 100}, {320, -8, 8}}); + if (doprocessMC) { hDecayChannel = qaRegistry.add("MC/hDecayChannel", ";Decay channel; ", HistType::kTH1D, {{2, -0.5, 1.5}}); hDecayChannel->GetXaxis()->SetBinLabel(1, "All"); @@ -379,17 +439,44 @@ struct lnnRecoTask { // QA for generated mother candidate and daughter particles h2FT0CPtGenColRecCandMC = qaRegistry.add("MC/RecMcCol/h2FT0CPtGenColRecCandMC", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {centAxis, candPtAxis}); h2FT0CPtGenColRecTrStrMC = qaRegistry.add("MC/RecMcCol/h2FT0CPtGenColRecTrStrMC", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {centAxis, candPtAxis}); - h2FT0CPtGenColRecPiStrMC = qaRegistry.add("MC/RecMcCol/h2FT0CPtGenColRecPiStrMC", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {centAxis, candPtAxis}); // QA signal generated candidate and daugher particles from secondary vertex h2FT0CPtGenColGenCandMC = qaRegistry.add("MC/QASignalGenSV/h2FT0CPtGenColGenCandMC", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {centAxis, candPtAxis}); h2FT0CPtGenColGenTrStrMC = qaRegistry.add("MC/QASignalGenSV/h2FT0CPtGenColGenTrStrMC", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {centAxis, candPtAxis}); - h2FT0CPtGenColGenPiStrMC = qaRegistry.add("MC/QASignalGenSV/h2FT0CPtGenColGenPiStrMC", ";FT0C (%);#it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {centAxis, candPtAxis}); } - hZvtx = qaRegistry.add("hZvtx", ";z_{vtx} (cm); ", HistType::kTH1D, {{100, -20, 20}}); - 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}}); - hCentFV0A = qaRegistry.add("hCentFV0A", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); + + if (doprocessSigEvtLossMC) { + // Infomation for all generated collisions + 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"); + hGenEventsNchEta08 = qaRegistry.add("QAEvent/hGenEventsNchEta08", ";;", HistType::kTH2D, {{multAxis}, {2, -0.5f, +1.5f}}); + hGenEventsNchEta08->GetYaxis()->SetBinLabel(1, "All gen. events"); + hGenEventsNchEta08->GetYaxis()->SetBinLabel(2, "Gen evts with at least 1 rec. collisions"); + hImpactParamGen = qaRegistry.add("QAEvent/McColAll/hImpactParamGen", "Impact parameter of generated MC events; Impact Parameter (b); Counts", HistType::kTH1D, {impactParamAxis}); + // Infomation for generated collisions which at least one rec. collision and passed the event selection + hImpactParamGenOneReco = qaRegistry.add("QAEvent/McColPassedEvSel/hImpactParamGenOneReco", "Impact parameter of generated MC events with at least one rec. evt and passed the event selection; Impact Parameter (b); Counts", HistType::kTH1D, {impactParamAxis}); + hGenOneRecoCentrality = qaRegistry.add("QAEvent/McColPassedEvSel/hGenOneRecoCentrality", "Centrality distribution of generated MC events with at least one rec. evt and passed the event selection; Centrality (FT0C %); Counts", HistType::kTH1D, {centFT0CAxis}); + hGenCentralityColvsMultiplicityGenEta08 = qaRegistry.add("QAEvent/McColPassedEvSel/hGenCentralityColvsMultiplicityGenEta08", "Correlation between FT0C centrality and charged particle multiplicity in generated MC events with at least one rec. evt and passed the event selection; Multiplicity #eta <0.8; Counts", HistType::kTH2D, {centFT0CAxis, multAxis}); + hGenCentralityColvsImpactParamGen = qaRegistry.add("QAEvent/McColPassedEvSel/hGenCentralityColvsImpactParamGen", "Correlation between FT0C centrality and impact parameter in generated MC events with at least one rec. evt and passed the event selection; Multiplicity #eta <0.8; Counts", HistType::kTH2D, {centFT0CAxis, impactParamAxis}); + hGenCentralityColvsFT0Cmultiplicity = qaRegistry.add("QAEvent/McColPassedEvSel/hGenCentralityColvsFT0Cmultiplicity", "Correlation between FT0C centrality and FT0C multiplicity in generated MC events with at least one rec. evt and passed the event selection; FT0c multiplicity", HistType::kTH2D, {centFT0CAxis, binsFT0CMultAxis}); + // Infomation for all reconstructed collisions passed the event selection (for Event Splitting) + 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}); + hRecoCentrality = qaRegistry.add("QAEvent/McColAll/hRecoCentrality", "Centrality distribution of reconstructed MC events passed the event selection; Centrality (FT0C %); Counts", HistType::kTH1D, {centFT0CAxis}); + hRecoCentralityColvsMultiplicityRecoEta08 = qaRegistry.add("QAEvent/McColAll/hRecoCentralityColvsMultiplicityRecoEta08", "Correlation between FT0C centrality and charged particle multiplicity in reconstructed MC events passed the event selection; Multiplicity #eta <0.8; Counts", HistType::kTH2D, {centFT0CAxis, multAxis}); + hRecoCentralityColvsImpactParamReco = qaRegistry.add("QAEvent/McColAll/hRecoCentralityColvsImpactParamReco", "Correlation between FT0C centrality and impact parameter in reconstructed MC events passed the event selection; Impact Parameter (b); Counts", HistType::kTH2D, {centFT0CAxis, impactParamAxis}); + hRecoCentralityColvsFT0Cmultiplicity = qaRegistry.add("QAEvent/McColAll/hRecoCentralityColvsFT0Cmultiplicity", "Correlation between FT0C centrality and FT0C multiplicity in reconstructed MC events passed the event selection; FT0C (%); FT0c multiplicity", HistType::kTH2D, {centFT0CAxis, binsFT0CMultAxis}); + // Information of generated Lnn in generated events + hGenLnnBeforeEvtSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnBeforeEvtSel", "Lnn generated #it{p}_{T} distribution in all gen evt;#it{p}_{T} (GeV/#it{c}); Counts", HistType::kTH1D, {ptAxis}); + hGenLnnvsImpactParameterBeforeEvtSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnvsImpactParameterBeforeEvtSel", "Correlation Lnn generated #it{p}_{T} and impact parameter in all gen evt;#it{p}_{T} (GeV/#it{c}); Impact parameter (b)", HistType::kTH2D, {ptAxis, impactParamAxis}); + hGenLnnvsMultiplicityGenEta08BeforeEvtSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnvsMultiplicityGenEta08BeforeEvtSel", "Correlation Lnn generated #it{p}_{T} and charged particle multiplicity in all gen evt;#it{p}_{T} (GeV/#it{c}); Multiplicity #eta <0.8", HistType::kTH2D, {ptAxis, multAxis}); + hGenLnnvsMultiplicityFT0CBeforeEvtSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnvsMultiplicityFT0CBeforeEvtSel", "Correlation Lnn 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 Lnn in generated events with at least one rec. event and passed the event selection + hGenLnnAfterSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnAfterSel", "Lnn generated #it{p}_{T} distribution in gen. evts with at least one rec. evt; #it{p}_{T} (GeV/#it{c}); Counts", HistType::kTH1D, {ptAxis}); + hGenLnnvsImpactParameterAfterSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnvsImpactParameterAfterSel", "Correlation Lnn 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}); + hGenLnnvsMultiplicityGenEta08AfterSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnvsMultiplicityGenEta08AfterSel", "Correlation Lnn 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.8", HistType::kTH2D, {ptAxis, multAxis}); + hGenLnnvsMultiplicityFT0CAfterSel = qaRegistry.add("QAEvent/McColLnn/hGenLnnvsMultiplicityFT0CAfterSel", "Correlation Lnn 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}); + } } // group BCs @@ -398,37 +485,37 @@ struct lnnRecoTask { if (mRunNumber == bc.runNumber()) { return; } - auto run3grp_timestamp = bc.timestamp(); + auto run3grptimestamp = bc.timestamp(); static const double kBzAutoThreshold = -990.; - o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); + o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(grpPath, run3grptimestamp); o2::parameters::GRPMagField* grpmag = 0x0; if (grpo) { o2::base::Propagator::initFieldFromGRP(grpo); if (dBzInput < kBzAutoThreshold) { // Fetch magnetic field from ccdb for current collision - d_bz = grpo->getNominalL3Field(); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; + dBz = grpo->getNominalL3Field(); + LOG(info) << "Retrieved GRP for timestamp " << run3grptimestamp << " with magnetic field of " << dBz << " kZG"; } else { - d_bz = dBzInput; + dBz = dBzInput; } } else { - grpmag = ccdb->getForTimeStamp(grpmagPath, run3grp_timestamp); + grpmag = ccdb->getForTimeStamp(grpmagPath, run3grptimestamp); if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp; + LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grptimestamp; } o2::base::Propagator::initFieldFromGRP(grpmag); if (dBzInput < kBzAutoThreshold) { // Fetch magnetic field from ccdb for current collision - d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; + dBz = std::lround(5.f * grpmag->getL3Current() / 30000.f); + LOG(info) << "Retrieved GRP for timestamp " << run3grptimestamp << " with magnetic field of " << dBz << " kZG"; } else { - d_bz = dBzInput; + dBz = dBzInput; } } if (!pidPath.value.empty()) { - auto h3pid = ccdb->getForTimeStamp>(pidPath.value + "_3H", run3grp_timestamp); + auto h3pid = ccdb->getForTimeStamp>(pidPath.value + "_3H", run3grptimestamp); std::copy(h3pid->begin(), h3pid->end(), mBBparams3H.begin()); } else { int kNBetheBlochParams = 5; @@ -438,7 +525,7 @@ struct lnnRecoTask { } mBBparams3H[kResolutionIndex] = cfgBetheBlochParams->get("3H", "resolution"); } - fitter.setBz(d_bz); + fitter.setBz(dBz); mRunNumber = bc.runNumber(); } @@ -502,15 +589,16 @@ struct lnnRecoTask { auto& h3Rigidity = lnnCand.isMatter ? posRigidity : negRigidity; auto& piRigidity = lnnCand.isMatter ? negRigidity : posRigidity; - // fill QA track histogram studies to check track signal before selections - h2FT0CnClusTPCtoTrBfSel->Fill(collision.centFT0C(), h3track.tpcNClsFound()); - h2FT0CnClusTPCtoPiBfSel->Fill(collision.centFT0C(), pitrack.tpcNClsFound()); - h2FT0CnTPCNClsCrossedRows3HBfSel->Fill(collision.centFT0C(), h3track.tpcNClsCrossedRows()); - h2FT0CnTPCNClsCrossedRowsPiBfSel->Fill(collision.centFT0C(), pitrack.tpcNClsCrossedRows()); - h2FT0Cchi2NClTPCtoTrBfSel->Fill(collision.centFT0C(), h3track.tpcChi2NCl()); - h2FT0Cchi2NClITStoTrBfSel->Fill(collision.centFT0C(), h3track.itsChi2NCl()); - if (doTrackQA) { + // fill QA track histogram studies to check track signal before selections + + h2FT0CnClusTPCtoTrBfSel->Fill(collision.centFT0C(), h3track.tpcNClsFound()); + h2FT0CnClusTPCtoPiBfSel->Fill(collision.centFT0C(), pitrack.tpcNClsFound()); + h2FT0CnTPCNClsCrossedRows3HBfSel->Fill(collision.centFT0C(), h3track.tpcNClsCrossedRows()); + h2FT0CnTPCNClsCrossedRowsPiBfSel->Fill(collision.centFT0C(), pitrack.tpcNClsCrossedRows()); + h2FT0Cchi2NClTPCtoTrBfSel->Fill(collision.centFT0C(), h3track.tpcChi2NCl()); + h2FT0Cchi2NClITStoTrBfSel->Fill(collision.centFT0C(), h3track.itsChi2NCl()); + bool passedTrTrackITS = h3track.hasITS(); bool passedTrTrackTOF = h3track.hasTOF(); bool passedPiTrackITS = pitrack.hasITS(); @@ -526,11 +614,13 @@ struct lnnRecoTask { } if (h3Rigidity < tpcRigidityMin3H || - h3track.tpcNClsCrossedRows() < nTPCNClsCrossedRows3H || h3track.tpcChi2NCl() < chi2nClusTPCMin || h3track.tpcChi2NCl() > chi2nClusTPCMax || h3track.itsChi2NCl() > chi2nClusITS || - pitrack.tpcNClsCrossedRows() < nTPCNClsCrossedRowsPi) { + h3track.tpcNClsCrossedRows() < nTPCNClsCrossedRows3H || + pitrack.tpcNClsCrossedRows() < nTPCNClsCrossedRowsPi || + h3track.tpcNClsFound() < nTPCClusMin3H || + pitrack.tpcNClsFound() < nTPCClusMinPi) { continue; } @@ -610,10 +700,10 @@ struct lnnRecoTask { continue; } // Definition of lnn mass - float mLNN_HypHI = 3.00; // , but 2993.7 MeV/c**2 + float mLnnHypHi = 3.00; // , but 2993.7 MeV/c**2 float massLNNL = std::sqrt(h3lE * h3lE - lnnMom[0] * lnnMom[0] - lnnMom[1] * lnnMom[1] - lnnMom[2] * lnnMom[2]); bool isLNNMass = false; - if (massLNNL > mLNN_HypHI - masswidth && massLNNL < mLNN_HypHI + masswidth) { + if (massLNNL > mLnnHypHi - masswidth && massLNNL < mLnnHypHi + masswidth) { isLNNMass = true; } if (!isLNNMass) { @@ -658,14 +748,16 @@ struct lnnRecoTask { // Fill 2D map after all selections hdEdx3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, h3track.tpcSignal()); - hdEdx3HPosTrack->Fill(lnnCand.mom3HTPC, h3track.tpcSignal()); - if (!lnnCand.isMatter) { - hdEdx3HNegTrack->Fill(-lnnCand.mom3HTPC, h3track.tpcSignal()); - } + if (lnnCand.isMatter) { + hdEdx3HPosTrack->Fill(lnnCand.mom3HTPC, h3track.tpcSignal()); + } else { + hdEdx3HNegTrack->Fill(-lnnCand.mom3HTPC, h3track.tpcSignal());} + hNsigma3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, lnnCand.nSigma3H); if (h3track.hasTOF()) { h3HSignalPtTOF->Fill(chargeFactor * h3track.pt(), beta); - hNsigma3HSelTOF->Fill(chargeFactor * h3track.p(), h3track.tofNSigmaTr()); + hNsigma3HSelTOF->Fill(chargeFactor * h3track.pt(), h3track.tofNSigmaTr()); + h2Nsigma3HSelTPCTOF->Fill(lnnCand.nSigma3H, h3track.tofNSigmaTr()); h3HMassPtTOF->Fill(chargeFactor * h3track.pt(), lnnCand.mass2TrTOF); } } @@ -674,7 +766,8 @@ struct lnnRecoTask { // Monte Carlo information void fillMCinfo(aod::McTrackLabels const& trackLabels, aod::McParticles const&) { - for (auto& lnnCand : lnnCandidates) { + for(size_t iCand = 0; iCand < lnnCandidates.size(); ++iCand) { + auto& lnnCand = lnnCandidates[iCand]; auto mcLabPos = trackLabels.rawIteratorAt(lnnCand.posTrackID); auto mcLabNeg = trackLabels.rawIteratorAt(lnnCand.negTrackID); @@ -705,8 +798,6 @@ struct lnnRecoTask { lnnCand.gMom3H = isTrTrack ? mcTrackPos.pVector() : mcTrackNeg.pVector(); - lnnCand.gMomPi = isTrTrack ? mcTrackNeg.pVector() : mcTrackPos.pVector(); - int kNDimGen = 3; for (int i = 0; i < kNDimGen; i++) { lnnCand.gDecVtx[i] = secVtx[i] - posPrimVtx[i]; @@ -744,36 +835,47 @@ struct lnnRecoTask { return motherIsAccepted; } + // Generated MC particles template - void fillMcHistograms(TracksFullMC::iterator const& mcTrack, aod::McParticles::iterator const& mc, Tcoll const& collision) + void fillGeneratedMcHistograms(aod::McParticles::iterator const& mc, Tcoll const& collision) { + const int pdg = std::abs(mc.pdgCode()); + int chargeFactor = -1 + 2 * (mc.pdgCode() > 0); + if (pdg == lnnPdg) { + hPtGeneratedLnnStr->Fill(collision.centFT0C(), chargeFactor * mc.pt()); + } + if (pdg == h3DauPdg) { + if (isLnnDecay(mc)) { + hPtGeneratedTrStr->Fill(collision.centFT0C(), chargeFactor * mc.pt()); + } + } + } + // Reconstructed tracks associated to MC generated particles + template + void fillMcHistograms(TracksFullMC::iterator const& mcTrack, aod::McParticles::iterator const& mc, Tcoll const& collision) + { bool passedTrackITS = mcTrack.hasITS(); bool passedTrackTPC = mcTrack.hasTPC(); bool passedTrackTOF = mcTrack.hasTOF(); int pdg = std::abs(mc.pdgCode()); - + int chargeFactor = -1 + 2 * (mc.pdgCode() > 0); + int chargeTrackFactor = mcTrack.sign(); if (pdg == h3DauPdg) { - hPtGeneratedTrStr->Fill(collision.centFT0C(), mc.pt()); - - if (passedTrackITS && passedTrackTPC) { - hPtItsTpcTrStr->Fill(collision.centFT0C(), mc.pt()); - hPtTrkItsTpcTrStr->Fill(collision.centFT0C(), mcTrack.pt()); - if (passedTrackTOF) { - hPtItsTpcTofTrStr->Fill(collision.centFT0C(), mc.pt()); - hPtTrkItsTpcTofTrStr->Fill(collision.centFT0C(), mcTrack.pt()); + if (passedTrackITS) { + hPtTpcTrStr->Fill(collision.centFT0C(), chargeFactor * mc.pt()); + hPtTrkTpcTrStr->Fill(collision.centFT0C(), chargeTrackFactor * mcTrack.pt()); + if (passedTrackTPC) { + hPtItsTpcTrStr->Fill(collision.centFT0C(), chargeFactor * mc.pt()); + hPtTrkItsTpcTrStr->Fill(collision.centFT0C(), chargeTrackFactor * mcTrack.pt()); + if (passedTrackTOF) { + hPtItsTpcTofTrStr->Fill(collision.centFT0C(), chargeFactor * mc.pt()); + hPtTrkItsTpcTofTrStr->Fill(collision.centFT0C(), chargeTrackFactor * mcTrack.pt()); + } } } } - - if (pdg == piDauPdg) { - hPtGeneratedPiStr->Fill(collision.centFT0C(), mc.pt()); - if (passedTrackITS && passedTrackTPC) { - hPtItsTpcPiStr->Fill(collision.centFT0C(), mc.pt()); - hPtTrkItsTpcPiStr->Fill(collision.centFT0C(), mcTrack.pt()); - } - } } void processData(CollisionsFull const& collisions, aod::V0s const& V0s, TracksFull const& tracks, aod::BCsWithTimestamps const&) @@ -798,6 +900,8 @@ struct lnnRecoTask { } hEvents->Fill(KevPileupCut); + h2GlobalTrackMultVsCent->Fill(collision.centFT0C(), collision.multNTracksGlobal()); + hZvtx->Fill(collision.posZ()); hCentFT0A->Fill(collision.centFT0A()); hCentFT0C->Fill(collision.centFT0C()); @@ -805,10 +909,10 @@ struct lnnRecoTask { hCentFV0A->Fill(collision.centFV0A()); const uint64_t collIdx = collision.globalIndex(); - auto v0Table_thisCollision = V0s.sliceBy(perCollision, collIdx); - v0Table_thisCollision.bindExternalIndices(&tracks); + auto v0TablethisCollision = V0s.sliceBy(perCollision, collIdx); + v0TablethisCollision.bindExternalIndices(&tracks); - fillCandidateData(collision, v0Table_thisCollision); + fillCandidateData(collision, v0TablethisCollision); for (const auto& lnnCand : lnnCandidates) { outputDataTable(collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), @@ -825,32 +929,21 @@ struct lnnRecoTask { } } } - PROCESS_SWITCH(lnnRecoTask, processData, "Data analysis", true); + PROCESS_SWITCH(LnnRecoTask, processData, "Data analysis", true); // MC process void processMC(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, aod::V0s const& V0s, aod::BCsWithTimestamps const&, TracksFull const& tracks, TracksFullMC const& tracksMC, aod::McTrackLabels const& trackLabelsMC, aod::McParticles const& particlesMC) { filledMothers.clear(); isGoodCollision.clear(); - isGoodCollision.resize(mcCollisions.size(), false); collisionFT0Ccent.clear(); - collisionFT0Ccent.resize(mcCollisions.size(), -1.f); recoCollisionIds.clear(); + + isGoodCollision.resize(mcCollisions.size(), false); + collisionFT0Ccent.resize(mcCollisions.size(), -1.f); recoCollisionIds.resize(mcCollisions.size(), -1); for (const auto& collision : collisions) { - for (auto const& trackMC : tracksMC) { - - if (!trackMC.has_mcParticle()) - continue; - - auto mc = trackMC.mcParticle(); - - if (!isLnnDecay(mc)) - continue; - - fillMcHistograms(trackMC, mc, collision); - } lnnCandidates.clear(); auto bc = collision.bc_as(); @@ -874,6 +967,8 @@ struct lnnRecoTask { } hEvents->Fill(KevPileupCut); + h2GlobalTrackMultVsCent->Fill(collision.centFT0C(), collision.multNTracksGlobal()); + hZvtx->Fill(collision.posZ()); hCentFT0A->Fill(collision.centFT0A()); hCentFT0C->Fill(collision.centFT0C()); @@ -886,12 +981,33 @@ struct lnnRecoTask { } const uint64_t collIdx = collision.globalIndex(); - auto v0Table_thisCollision = V0s.sliceBy(perCollision, collIdx); - v0Table_thisCollision.bindExternalIndices(&tracks); + auto v0TablethisCollision = V0s.sliceBy(perCollision, collIdx); + v0TablethisCollision.bindExternalIndices(&tracks); - fillCandidateData(collision, v0Table_thisCollision); + fillCandidateData(collision, v0TablethisCollision); fillMCinfo(trackLabelsMC, particlesMC); + // Fill MC gen. particles producing to the current MC collision + const auto mcCollId = collision.mcCollisionId(); + auto particlesThisMcCollision = particlesMC.sliceBy(mcParticlesPerMcCollision, mcCollId); + for (const auto& mc : particlesThisMcCollision) { + fillGeneratedMcHistograms(mc, collision); + } + + // Fill only tracks producing to the current reconstructed collision + auto tracksThisCollision = tracksMC.sliceBy( tracksPerCollision, collIdx); + for (const auto& trackMC : tracksThisCollision) { + if (!trackMC.has_mcParticle()) { + continue; + } + auto mc = trackMC.mcParticle(); + // Select only daughters from the Λnn decay + if (!isLnnDecay(mc)) { + continue; + } + fillMcHistograms(trackMC, mc,collision); + } + for (const auto& lnnCand : lnnCandidates) { if (!lnnCand.isSignal && mcSignalOnly) { continue; @@ -900,7 +1016,6 @@ struct lnnRecoTask { if (lnnCand.recoMcColl && lnnCand.survEvSelection) { h2FT0CPtGenColRecCandMC->Fill(collision.centFT0C(), chargeFactor * lnnCand.genPt()); h2FT0CPtGenColRecTrStrMC->Fill(collision.centFT0C(), chargeFactor * lnnCand.genPt3H()); - h2FT0CPtGenColRecPiStrMC->Fill(collision.centFT0C(), chargeFactor * lnnCand.genPtPi()); } outputMCTable(collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), collision.posX(), collision.posY(), collision.posZ(), @@ -926,17 +1041,19 @@ struct lnnRecoTask { } float cent = collisionFT0Ccent[mcPart.mcCollisionId()]; - h2FT0CPtGenColGenCandMC->Fill(cent, mcPart.pt()); + int chargeFactorGen = mcPart.pdgCode() > 0 ? 1 : -1; + + h2FT0CPtGenColGenCandMC->Fill(cent, chargeFactorGen * mcPart.pt()); constexpr std::size_t KVtxDim = 3; std::array secVtx; std::array primVtx = {mcPart.vx(), mcPart.vy(), mcPart.vz()}; - constexpr std::size_t kArrayDim = 3; - std::array momMother = mcPart.pVector(); + constexpr std::size_t KArrayDim = 3; + std::array momMother = mcPart.pVector(); + + std::array mom3H; - std::array mom3H; - std::array momPi; bool is3HFound = false; for (const auto& mcDaught : mcPart.daughters_as()) { @@ -945,13 +1062,13 @@ struct lnnRecoTask { if (pdg == h3DauPdg) { secVtx = {mcDaught.vx(), mcDaught.vy(), mcDaught.vz()}; mom3H = mcDaught.pVector(); - h2FT0CPtGenColGenTrStrMC->Fill(cent, mcDaught.pt()); + h2FT0CPtGenColGenTrStrMC->Fill(cent, chargeFactorGen * mcDaught.pt()); is3HFound = true; break; } - if (pdg == piDauPdg) { - momPi = mcDaught.pVector(); - h2FT0CPtGenColGenPiStrMC->Fill(cent, mcDaught.pt()); + + if (is3HFound) { + break; } } @@ -973,6 +1090,7 @@ struct lnnRecoTask { LnnCandidate lnnCand; lnnCand.pdgCode = mcPart.pdgCode(); lnnCand.survEvSelection = isGoodCollision[mcPart.mcCollisionId()]; + int chargeFactor = -1 + 2 * (lnnCand.pdgCode > 0); int kDimGen = 3; @@ -980,7 +1098,6 @@ struct lnnRecoTask { lnnCand.gDecVtx[i] = secVtx[i] - primVtx[i]; lnnCand.gMom[i] = momMother[i]; lnnCand.gMom3H[i] = mom3H[i]; - lnnCand.gMomPi[i] = momPi[i]; } lnnCand.posTrackID = -1; @@ -1007,12 +1124,134 @@ struct lnnRecoTask { lnnCand.gDecVtx[0], lnnCand.gDecVtx[1], lnnCand.gDecVtx[2], lnnCand.isReco, lnnCand.isSignal, lnnCand.recoMcColl, lnnCand.survEvSelection); } } - PROCESS_SWITCH(lnnRecoTask, processMC, "MC analysis", false); + PROCESS_SWITCH(LnnRecoTask, processMC, "MC analysis", false); + + template + bool passEvtSel(const CollType& collision) + { + if (!collision.sel8()) + return false; + + if ((std::abs(collision.posZ())) > vtxZCut) + return false; + + if (useNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) + return false; + + return true; + } + + void processSigEvtLossMC(McCollisionMults::iterator const& mcCollision, soa::SmallGroups const& collisions, aod::McParticles const& GenParticles) + { + if (std::abs(mcCollision.posZ()) > vtxZCut) { + return; + } + //////////// Event loss estimation via impact parameter and multiplicity by FT0C + + // Fill all generated events + hEvtMC->Fill(0); + hImpactParamGen->Fill(mcCollision.impactParameter()); + hGenEventsNchEta08->Fill(mcCollision.multMCNParticlesEta08(), 0); + + // 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.; + int biggestNContribs = -1; + + for (auto const& col : collisions) { + if (!passEvtSel(col)) { + continue; + } + // In case of multiple reconstructed collisions associated to the same generated one, only consider the one with the biggest number of contributors + if (biggestNContribs < col.numContrib()) { + biggestNContribs = col.numContrib(); + centralityFT0C = col.centFT0C(); + } + atLeastOneRecoEvt = true; + hImpactParamReco->Fill(mcCollision.impactParameter()); + hRecoCentrality->Fill(col.centFT0C()); + hRecoCentralityColvsMultiplicityRecoEta08->Fill(col.centFT0C(), mcCollision.multMCNParticlesEta08()); + hRecoCentralityColvsImpactParamReco->Fill(col.centFT0C(), mcCollision.impactParameter()); + hRecoCentralityColvsFT0Cmultiplicity->Fill(col.centFT0C(), mcCollision.multMCFT0C()); + } + + if (atLeastOneRecoEvt) { + hEvtMC->Fill(2); + hGenEventsNchEta08->Fill(mcCollision.multMCNParticlesEta08(), 1); + hImpactParamGenOneReco->Fill(mcCollision.impactParameter()); + hGenOneRecoCentrality->Fill(centralityFT0C); + hGenCentralityColvsMultiplicityGenEta08->Fill(centralityFT0C, mcCollision.multMCNParticlesEta08()); + hGenCentralityColvsImpactParamGen->Fill(centralityFT0C, mcCollision.impactParameter()); + hGenCentralityColvsFT0Cmultiplicity->Fill(centralityFT0C, mcCollision.multMCFT0C()); + } + // 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()) != lnnPdg) + continue; + + auto daughters = genParticle.daughters_as(); + + bool dauTr = false; + bool dauPi = false; + + int trSign = 0; + int piSign = 0; + + for (const auto& daughter : daughters) { + + if (std::abs(daughter.pdgCode()) == h3DauPdg) { + + dauTr = true; + trSign = daughter.pdgCode() > 0 ? 1 : -1; + + daugh1 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), h3Mass); + + } + else if (std::abs(daughter.pdgCode()) == piDauPdg) { + + dauPi = true; + piSign = daughter.pdgCode() > 0 ? 1 : -1; + + daugh2 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), piMass); + } + } + + if (!dauTr || !dauPi || (trSign * piSign) > 0) { + continue; + } + + mother = daugh1 + daugh2; + + // Fill informations for generated 3HL in all generated events + hGenLnnBeforeEvtSel->Fill(mother.pt()); + hGenLnnvsImpactParameterBeforeEvtSel->Fill(mother.pt(), mcCollision.impactParameter()); + hGenLnnvsMultiplicityGenEta08BeforeEvtSel->Fill(mother.pt(), mcCollision.multMCNParticlesEta08()); + hGenLnnvsMultiplicityFT0CBeforeEvtSel->Fill(mother.pt(), mcCollision.multMCFT0C()); + + // Fill informations for generated 3HL in generated events with at least one reconstructed event + if (atLeastOneRecoEvt) { + hGenLnnAfterSel->Fill(mother.pt()); + hGenLnnvsImpactParameterAfterSel->Fill(mother.pt(), mcCollision.impactParameter()); + hGenLnnvsMultiplicityGenEta08AfterSel->Fill(mother.pt(), mcCollision.multMCNParticlesEta08()); + hGenLnnvsMultiplicityFT0CAfterSel->Fill(mother.pt(), mcCollision.multMCFT0C()); + } + } + } + PROCESS_SWITCH(LnnRecoTask, processSigEvtLossMC, "Signal and event loss analysis", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc)}; } From 0bc8f549d0e8963835aae5b17a3740663ff791ce Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 2 Sep 2026 13:29:10 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx | 29 +++++++++++----------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx b/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx index 5e3eaa7c85c..4c1a56ce899 100644 --- a/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx +++ b/PWGLF/TableProducer/Nuspex/lnnRecoTask.cxx @@ -19,8 +19,8 @@ #include "Common/Core/RecoDecay.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" -#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -67,7 +67,7 @@ using std::array; using TracksFull = soa::Join; using TracksFullMC = soa::Join; using CollisionsFull = soa::Join; -using CollisionsFullMC = soa::Join; +using CollisionsFullMC = soa::Join; using CollisionCandidates = o2::soa::Join; using CollisionCandidatesMC = o2::soa::Join; @@ -131,7 +131,7 @@ std::shared_ptr h2FT0CPtGenColRecTrStrMC; // QA signal generated candidate and daugher particles from secondary vertex std::shared_ptr h2FT0CPtGenColGenCandMC; std::shared_ptr h2FT0CPtGenColGenTrStrMC; -// Event and signal loss +// Event and signal loss std::shared_ptr hEvtMC; std::shared_ptr hImpactParamGen; std::shared_ptr hImpactParamGenOneReco; @@ -374,7 +374,6 @@ struct LnnRecoTask { const AxisSpec multAxis = {binsMult, "Multiplicity #eta <0.5"}; const AxisSpec globalTrackMultAxis = {binsglobalTrackMult, "N_{global track mult.}"}; - hNsigma3HSel = qaRegistry.add("PID/hNsigma3HSel", ";#it{p}^{TPC}/z (GeV/#it{c}); n_{#sigma} ({}^{3}H)", HistType::kTH2F, {rigidityAxis, nSigma3HtpcAxis}); hNsigma3HSelTOF = qaRegistry.add("PID/hNsigma3HSelTOF", ";#it{p}_{T} (GeV/#it{c}); n_{#sigma} ({}^{3}H)", HistType::kTH2F, {tPtAxis, nSigma3HtpcAxis}); h2Nsigma3HSelTPCTOF = qaRegistry.add("PID/h2Nsigma3HSelTPCTOF", ";n_{#sigma TPC} ({}^{3}H); n_{#sigma TOF} ({}^{3}H)", HistType::kTH2F, {nSigma3HtpcAxis, nSigma3HtofAxis}); @@ -393,7 +392,7 @@ struct LnnRecoTask { hCentFT0C = qaRegistry.add("hCentFT0C", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); hCentFT0M = qaRegistry.add("hCentFT0M", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); hCentFV0A = qaRegistry.add("hCentFV0A", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}}); - + hLnnCandLoss = qaRegistry.add("CandCounts/hLnnCandLoss", ";CandLoss; ", HistType::kTH1D, {{7, -0.5, 6.5}}); hEvents->GetXaxis()->SetBinLabel(1, "All"); hEvents->GetXaxis()->SetBinLabel(2, "sel8"); @@ -751,8 +750,9 @@ struct LnnRecoTask { if (lnnCand.isMatter) { hdEdx3HPosTrack->Fill(lnnCand.mom3HTPC, h3track.tpcSignal()); } else { - hdEdx3HNegTrack->Fill(-lnnCand.mom3HTPC, h3track.tpcSignal());} - + hdEdx3HNegTrack->Fill(-lnnCand.mom3HTPC, h3track.tpcSignal()); + } + hNsigma3HSel->Fill(chargeFactor * lnnCand.mom3HTPC, lnnCand.nSigma3H); if (h3track.hasTOF()) { h3HSignalPtTOF->Fill(chargeFactor * h3track.pt(), beta); @@ -766,7 +766,7 @@ struct LnnRecoTask { // Monte Carlo information void fillMCinfo(aod::McTrackLabels const& trackLabels, aod::McParticles const&) { - for(size_t iCand = 0; iCand < lnnCandidates.size(); ++iCand) { + for (size_t iCand = 0; iCand < lnnCandidates.size(); ++iCand) { auto& lnnCand = lnnCandidates[iCand]; auto mcLabPos = trackLabels.rawIteratorAt(lnnCand.posTrackID); auto mcLabNeg = trackLabels.rawIteratorAt(lnnCand.negTrackID); @@ -995,7 +995,7 @@ struct LnnRecoTask { } // Fill only tracks producing to the current reconstructed collision - auto tracksThisCollision = tracksMC.sliceBy( tracksPerCollision, collIdx); + auto tracksThisCollision = tracksMC.sliceBy(tracksPerCollision, collIdx); for (const auto& trackMC : tracksThisCollision) { if (!trackMC.has_mcParticle()) { continue; @@ -1005,7 +1005,7 @@ struct LnnRecoTask { if (!isLnnDecay(mc)) { continue; } - fillMcHistograms(trackMC, mc,collision); + fillMcHistograms(trackMC, mc, collision); } for (const auto& lnnCand : lnnCandidates) { @@ -1090,7 +1090,7 @@ struct LnnRecoTask { LnnCandidate lnnCand; lnnCand.pdgCode = mcPart.pdgCode(); lnnCand.survEvSelection = isGoodCollision[mcPart.mcCollisionId()]; - + int chargeFactor = -1 + 2 * (lnnCand.pdgCode > 0); int kDimGen = 3; @@ -1214,9 +1214,8 @@ struct LnnRecoTask { trSign = daughter.pdgCode() > 0 ? 1 : -1; daugh1 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), h3Mass); - - } - else if (std::abs(daughter.pdgCode()) == piDauPdg) { + + } else if (std::abs(daughter.pdgCode()) == piDauPdg) { dauPi = true; piSign = daughter.pdgCode() > 0 ? 1 : -1; @@ -1226,7 +1225,7 @@ struct LnnRecoTask { } if (!dauTr || !dauPi || (trSign * piSign) > 0) { - continue; + continue; } mother = daugh1 + daugh2;