diff --git a/PWGHF/D2H/Tasks/taskLc.cxx b/PWGHF/D2H/Tasks/taskLc.cxx index 0334f01ada9..6c7ce834e12 100644 --- a/PWGHF/D2H/Tasks/taskLc.cxx +++ b/PWGHF/D2H/Tasks/taskLc.cxx @@ -63,6 +63,7 @@ using namespace o2::framework::expressions; using namespace o2::hf_centrality; using namespace o2::hf_occupancy; using namespace o2::hf_evsel; +using namespace o2::constants::physics; /// Λc± → p± K∓ π± analysis task struct HfTaskLc { @@ -74,7 +75,7 @@ struct HfTaskLc { Configurable fillTHn{"fillTHn", false, "fill THn"}; Configurable storeOccupancy{"storeOccupancy", true, "Flag to store occupancy information"}; Configurable occEstimator{"occEstimator", 2, "Occupancy estimation (None: 0, ITS: 1, FT0C: 2)"}; - Configurable storeProperLifetime{"storeProperLifetime", false, "Flag to store proper lifetime"}; + Configurable storeProperDecayTime{"storeProperDecayTime", false, "Flag to store proper decay time"}; // CCDB configuration Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "Path of the grp file (Run 2)"}; @@ -117,16 +118,9 @@ struct HfTaskLc { ConfigurableAxis thnConfigAxisGenPtB{"thnConfigAxisGenPtB", {1000, 0, 100}, "Gen Pt B"}; ConfigurableAxis thnConfigAxisNumPvContr{"thnConfigAxisNumPvContr", {200, -0.5, 199.5}, "Number of PV contributors"}; ConfigurableAxis thnConfigAxisOccupancy{"thnConfigAxisOccupancy", {14, 0, 14000}, "axis for centrality"}; - ConfigurableAxis thnConfigAxisProperLifetime{"thnConfigAxisProperLifetime", {200, 0, 2}, "Proper lifetime, ps"}; + ConfigurableAxis thnConfigAxisProperDecayTime{"thnConfigAxisProperDecayTime", {200, 0, 2}, "Proper decay time, ps"}; HistogramRegistry registry{"registry", {}}; - // Factors for conversion between units - constexpr static float CtToProperLifetimePs = 1.f / o2::constants::physics::LightSpeedCm2PS; - constexpr static float NanoToPico = 1000.f; - // Names of folders and suffixes for MC signal histograms - constexpr static std::string_view SignalFolders[] = {"signal", "prompt", "nonprompt"}; - constexpr static std::string_view SignalSuffixes[] = {"", "Prompt", "NonPrompt"}; - enum MlClasses : int { MlClassBackground = 0, MlClassPrompt, @@ -137,9 +131,16 @@ struct HfTaskLc { enum SignalClasses : int { Signal = 0, Prompt, - NonPrompt + NonPrompt, + NumberOfSignalClasses }; + // Factor for conversion between units + constexpr static float CtToProperDecayTimePs = 1.f / o2::constants::physics::LightSpeedCm2PS; + // Names of folders and suffixes for MC signal histograms + constexpr static std::array SignalFolders = {"signal", "prompt", "nonprompt"}; + constexpr static std::array SignalSuffixes = {"", "Prompt", "NonPrompt"}; + void init(InitContext&) { const std::array doprocess{doprocessDataStd, doprocessDataStdWithFT0C, doprocessDataStdWithFT0M, doprocessDataWithMl, doprocessDataWithMlWithFT0C, doprocessDataWithMlWithFT0M, doprocessMcStd, doprocessMcStdWithFT0C, doprocessMcStdWithFT0M, doprocessMcWithMl, doprocessMcWithMlWithFT0C, doprocessMcWithMlWithFT0M}; @@ -186,8 +187,8 @@ struct HfTaskLc { addHistogramsRec("hDecLength", "decay length (cm)", "entries", {HistType::kTH1F, {{400, 0., 1.}}}); /// decay length xy candidate addHistogramsRec("hDecLengthxy", "decay length xy (cm)", "entries", {HistType::kTH1F, {{400, 0., 1.}}}); - /// proper lifetime - addHistogramsRec("hCt", "proper lifetime (#Lambda_{c}) * #it{c} (cm)", "entries", {HistType::kTH1F, {{100, 0., 0.2}}}); + /// proper decay time + addHistogramsRec("hCt", "proper decay time (#Lambda_{c}) * #it{c} (cm)", "entries", {HistType::kTH1F, {{100, 0., 0.2}}}); /// cosine of pointing angle addHistogramsRec("hCPA", "cosine of pointing angle", "entries", {HistType::kTH1F, {{110, -1.1, 1.1}}}); /// cosine of pointing angle xy @@ -219,8 +220,8 @@ struct HfTaskLc { /// decay length xy candidate addHistogramsRec("hDecLengthxyVsPt", "decay length xy (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{400, 0., 1.}, {vbins}}}); - /// proper lifetime - addHistogramsRec("hCtVsPt", "proper lifetime (#Lambda_{c}) * #it{c} (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{100, 0., 0.2}, {vbins}}}); + /// proper decay time + addHistogramsRec("hCtVsPt", "proper decay time (#Lambda_{c}) * #it{c} (cm)", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{100, 0., 0.2}, {vbins}}}); /// cosine of pointing angle addHistogramsRec("hCPAVsPt", "cosine of pointing angle", "#it{p}_{T} (GeV/#it{c})", {HistType::kTH2F, {{110, -1.1, 1.1}, {vbins}}}); @@ -268,7 +269,8 @@ struct HfTaskLc { const AxisSpec thnAxisPtB{thnConfigAxisGenPtB, "#it{p}_{T}^{B} (GeV/#it{c})"}; const AxisSpec thnAxisTracklets{thnConfigAxisNumPvContr, "Number of PV contributors"}; const AxisSpec thnAxisOccupancy{thnConfigAxisOccupancy, "Occupancy"}; - const AxisSpec thnAxisProperLifetime{thnConfigAxisProperLifetime, "T_{proper} (ps)"}; + const AxisSpec thnAxisProperDecayTime{thnConfigAxisProperDecayTime, "#it{t}_{proper} (ps)"}; + const AxisSpec thnAxisProperDecayTimeGen{thnConfigAxisProperDecayTime, "#it{t}_{proper, gen} (ps)"}; bool const isDataWithMl = doprocessDataWithMl || doprocessDataWithMlWithFT0C || doprocessDataWithMlWithFT0M; bool const isMcWithMl = doprocessMcWithMl || doprocessMcWithMlWithFT0C || doprocessMcWithMlWithFT0M; @@ -300,10 +302,13 @@ struct HfTaskLc { } } } - if (storeProperLifetime) { + if (storeProperDecayTime) { for (const auto& axes : std::array*, 3>{&axesWithBdt, &axesStd, &axesGen}) { if (!axes->empty()) { - axes->push_back(thnAxisProperLifetime); + axes->push_back(thnAxisProperDecayTime); + if (!isData && axes != &axesGen) { + axes->push_back(thnAxisProperDecayTimeGen); + } } } } @@ -334,13 +339,34 @@ struct HfTaskLc { return o2::hf_centrality::getCentralityColl(collision); } + /// Evaluate decay time of generated particle + /// \param mcParticleProng0 one of generated particle's daughters + /// \param motherParticle generated particle + /// \return decay time in picoseconds. For nonprompt particles it is evaluated as if it was prompt (as it is calculated for data) + float evaluateMcGenDecayTime(const McParticles3ProngMatched::iterator& mcParticleProng0, const McParticles3ProngMatched::iterator& motherParticle) + { + const auto mcCollision = motherParticle.template mcCollision_as(); + const float pMother = motherParticle.p(); + const float pvX = mcCollision.posX(); + const float pvY = mcCollision.posY(); + const float pvZ = mcCollision.posZ(); + const float svX = mcParticleProng0.vx(); + const float svY = mcParticleProng0.vy(); + const float svZ = mcParticleProng0.vz(); + + const float decayLength = static_cast(RecoDecay::distance(std::array{svX, svY, svZ}, std::array{pvX, pvY, pvZ})); + const float properDecayTime = decayLength * static_cast(MassLambdaCPlus) / LightSpeedCm2PS / pMother; + + return properDecayTime; + } + /// Helper function for filling MC reconstructed histograms for prompt, nonpromt and common (signal) /// \param candidate is a reconstructed candidate /// \tparam SignalType is an enum defining which histogram in which folder (signal, prompt or nonpromt) to fill template void fillHistogramsRecSig(CandidateType const& candidate) { - const auto& mcParticleProng0 = candidate.template prong0_as().template mcParticle_as>(); + const auto& mcParticleProng0 = candidate.template prong0_as().template mcParticle_as(); const auto pdgCodeProng0 = std::abs(mcParticleProng0.pdgCode()); if ((candidate.isSelLcToPKPi() >= selectionFlagLc) && pdgCodeProng0 == kProton) { registry.fill(HIST("MC/reconstructed/") + HIST(SignalFolders[SignalType]) + HIST("/hMassRecSig") + HIST(SignalSuffixes[SignalType]), HfHelper::invMassLcToPKPi(candidate)); @@ -385,8 +411,8 @@ struct HfTaskLc { /// Fill MC histograms at reconstruction level /// \tparam FillMl switch to fill ML histograms - template - void fillHistosMcRec(CollType const& collision, CandLcMcRec const& candidates, CandLcMcGen const& mcParticles) + template + void fillHistosMcRec(CollType const& collision, CandLcMcRec const& candidates, McParticles3ProngMatched const& mcParticles) { const auto thisCollId = collision.globalIndex(); const auto& groupedLcCandidates = candidates.sliceBy(candLcPerCollision, thisCollId); @@ -403,7 +429,7 @@ struct HfTaskLc { if (std::abs(candidate.flagMcMatchRec()) == hf_decay::hf_cand_3prong::DecayChannelMain::LcToPKPi) { // Get the corresponding MC particle. - const auto& mcParticleProng0 = candidate.template prong0_as().template mcParticle_as>(); + const auto& mcParticleProng0 = candidate.template prong0_as().template mcParticle_as(); const auto pdgCodeProng0 = std::abs(mcParticleProng0.pdgCode()); const auto indexMother = RecoDecay::getMother(mcParticles, mcParticleProng0, o2::constants::physics::Pdg::kLambdaCPlus, true); const auto particleMother = mcParticles.rawIteratorAt(indexMother); @@ -420,6 +446,8 @@ struct HfTaskLc { const auto numPvContributors = collision.numContrib(); const auto ptRecB = candidate.ptBhadMotherPart(); + const float properDecayTimeGen = evaluateMcGenDecayTime(mcParticleProng0, particleMother); + /// MC reconstructed signal fillHistogramsRecSig(candidate); @@ -438,7 +466,7 @@ struct HfTaskLc { occ = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); } double outputBkg(-1), outputPrompt(-1), outputFD(-1); - const float properLifetime = HfHelper::ctLc(candidate) * CtToProperLifetimePs; + const float properDecayTime = HfHelper::ctLc(candidate) * CtToProperDecayTimePs; auto fillTHnRecSig = [&](bool isPKPi) { const auto massLc = isPKPi ? HfHelper::invMassLcToPKPi(candidate) : HfHelper::invMassLcToPiKP(candidate); @@ -461,8 +489,9 @@ struct HfTaskLc { if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { valuesToFill.push_back(occ); } - if (storeProperLifetime) { - valuesToFill.push_back(properLifetime); + if (storeProperDecayTime) { + valuesToFill.push_back(properDecayTime); + valuesToFill.push_back(properDecayTimeGen); } if constexpr (FillMl) { registry.get(HIST("hnLcVarsWithBdt"))->Fill(valuesToFill.data()); @@ -485,8 +514,8 @@ struct HfTaskLc { /// Helper function for filling MC generated histograms for prompt, nonpromt and common (signal) /// \param particle is a generated particle /// \tparam SignalType is an enum defining which histogram in which folder (signal, prompt or nonpromt) to fill - template - void fillHistogramsGen(ParticleType const& particle) + template + void fillHistogramsGen(McParticles3ProngMatched::iterator const& particle) { registry.fill(HIST("MC/generated/") + HIST(SignalFolders[SignalType]) + HIST("/hPtGen") + HIST(SignalSuffixes[SignalType]), particle.pt()); registry.fill(HIST("MC/generated/") + HIST(SignalFolders[SignalType]) + HIST("/hEtaGen") + HIST(SignalSuffixes[SignalType]), particle.eta()); @@ -498,8 +527,8 @@ struct HfTaskLc { } /// Fill MC histograms at generated level - template - void fillHistosMcGen(CandLcMcGen const& mcParticles, Coll const& recoCollisions) + template + void fillHistosMcGen(McParticles3ProngMatched const& mcParticles, Coll const& recoCollisions) { // MC gen. for (const auto& particle : mcParticles) { @@ -522,10 +551,9 @@ struct HfTaskLc { occ = o2::hf_occupancy::getOccupancyGenColl(recoCollsPerMcColl, occEstimator); } - const auto& mcDaughter0 = particle.template daughters_as>().begin(); - const float p2m = particle.p() / o2::constants::physics::MassLambdaCPlus; - const float gamma = std::sqrt(1 + p2m * p2m); // mother's particle Lorentz factor - const float properLifetime = mcDaughter0.vt() * NanoToPico / gamma; // from ns to ps * from lab time to proper time + const auto mcDaughter0 = particle.template daughters_as().begin(); + + const float properDecayTime = evaluateMcGenDecayTime(mcDaughter0, particle); fillHistogramsGen(particle); @@ -537,8 +565,8 @@ struct HfTaskLc { if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { valuesToFill.push_back(occ); } - if (storeProperLifetime) { - valuesToFill.push_back(properLifetime); + if (storeProperDecayTime) { + valuesToFill.push_back(properDecayTime); } registry.get(HIST("hnLcVarsGen"))->Fill(valuesToFill.data()); } @@ -631,7 +659,7 @@ struct HfTaskLc { occ = o2::hf_occupancy::getOccupancyColl(collision, occEstimator); } double outputBkg(-1), outputPrompt(-1), outputFD(-1); - const float properLifetime = HfHelper::ctLc(candidate) * CtToProperLifetimePs; + const float properDecayTime = HfHelper::ctLc(candidate) * CtToProperDecayTimePs; auto fillTHnData = [&](bool isPKPi) { const auto massLc = isPKPi ? HfHelper::invMassLcToPKPi(candidate) : HfHelper::invMassLcToPiKP(candidate); @@ -654,8 +682,8 @@ struct HfTaskLc { if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) { valuesToFill.push_back(occ); } - if (storeProperLifetime) { - valuesToFill.push_back(properLifetime); + if (storeProperDecayTime) { + valuesToFill.push_back(properDecayTime); } if constexpr (FillMl) { registry.get(HIST("hnLcVarsWithBdt"))->Fill(valuesToFill.data()); @@ -687,10 +715,10 @@ struct HfTaskLc { /// Run the analysis on MC data /// \tparam FillMl switch to fill ML histograms - template + template void runAnalysisPerCollisionMc(CollType const& collisions, CandType const& candidates, - CandLcMcGen const& mcParticles) + McParticles3ProngMatched const& mcParticles) { for (const auto& collision : collisions) { // MC Rec. diff --git a/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx b/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx index 6bf2d001097..7b0eb1f6a24 100644 --- a/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx @@ -144,8 +144,8 @@ DECLARE_SOA_COLUMN(Chi2Topo, chi2Topo, float); //! chi DECLARE_SOA_COLUMN(DecayLength, decayLength, float); //! decay length, cm DECLARE_SOA_COLUMN(DecayLengthError, decayLengthError, float); //! decay length error DECLARE_SOA_COLUMN(DecayLengthNormalised, decayLengthNormalised, float); //! decay length over its error -DECLARE_SOA_COLUMN(T, t, float); //! proper lifetime, ps -DECLARE_SOA_COLUMN(ErrT, errT, float); //! lifetime error +DECLARE_SOA_COLUMN(T, t, float); //! proper decay time, ps +DECLARE_SOA_COLUMN(ErrT, errT, float); //! decay time error DECLARE_SOA_COLUMN(MassInv, massInv, float); //! invariant mass DECLARE_SOA_COLUMN(P, p, float); //! momentum DECLARE_SOA_COLUMN(Pt, pt, float); //! transverse momentum @@ -183,7 +183,7 @@ DECLARE_SOA_COLUMN(XDecay, xDecay, float); //! Secondary (decay) vertex X coordi DECLARE_SOA_COLUMN(YDecay, yDecay, float); //! Secondary (decay) vertex Y coordinate, cm DECLARE_SOA_COLUMN(ZDecay, zDecay, float); //! Secondary (decay) vertex Z coordinate, cm DECLARE_SOA_COLUMN(LDecay, lDecay, float); //! Decay length, cm (distance between PV and SV, curvature is neglected) -DECLARE_SOA_COLUMN(TDecay, tDecay, float); //! Proper lifetime, ps +DECLARE_SOA_COLUMN(TDecay, tDecay, float); //! Proper decay time, ps DECLARE_SOA_COLUMN(XEvent, xEvent, float); //! Primary (event) vertex X coordinate, cm DECLARE_SOA_COLUMN(YEvent, yEvent, float); //! Primary (event) vertex Y coordinate, cm DECLARE_SOA_COLUMN(ZEvent, zEvent, float); //! Primary (event) vertex Z coordinate, cm @@ -412,7 +412,6 @@ struct HfTreeCreatorLcToPKPi { constexpr static float UndefValueFloat = -999.f; constexpr static int UndefValueInt = -999; - constexpr static float NanoToPico = 1000.f; using TracksWPid = soa::Join; using Cents = soa::Join; @@ -420,9 +419,9 @@ struct HfTreeCreatorLcToPKPi { // number showing MC status of the candidate (signal or background, prompt or non-prompt etc.) enum SigBgStatus : int { Background = 0, // combinatorial background, at least one of the prongs do not originate from the Lc decay - Prompt, // signal with Lc produced directly in the event - NonPrompt, // signal with Lc produced aftewards the event, e.g. during decay of beauty particle - WrongOrder, // all the prongs are from Lc decay, but proton and pion hypothesis are swapped + Prompt = 1, // signal with Lc produced directly in the event + NonPrompt = 2, // signal with Lc produced aftewards the event, e.g. during decay of beauty particle + WrongOrder = 3, // all the prongs are from Lc decay, but proton and pion hypothesis are swapped Default = -1 // impossible, should not be the case, to catch logical error if any }; @@ -911,7 +910,7 @@ struct HfTreeCreatorLcToPKPi { const float deltaP = std::sqrt(pt * pt * deltaPt * deltaPt + candidate.kfPz() * candidate.kfPz() * candidate.kfErrorPz() * candidate.kfErrorPz()) / p; - const float lifetime = decayLength * static_cast(MassLambdaCPlus) / LightSpeedCm2PS / p; + const float decayTime = decayLength * static_cast(MassLambdaCPlus) / LightSpeedCm2PS / p; const float deltaT = dl * static_cast(MassLambdaCPlus) / LightSpeedCm2PS / p; rowCandidateKF( svX, svY, svZ, svErrX, svErrY, svErrZ, @@ -919,7 +918,7 @@ struct HfTreeCreatorLcToPKPi { chi2primProton, chi2primKaon, chi2primPion, dcaProtonKaon, dcaProtonPion, dcaPionKaon, chi2GeoProtonKaon, chi2GeoProtonPion, chi2GeoPionKaon, - chi2Geo, chi2Topo, decayLength, dl, decayLength / dl, lifetime, deltaT, + chi2Geo, chi2Topo, decayLength, dl, decayLength / dl, decayTime, deltaT, mass, p, pt, deltaP, deltaPt, functionSelection, sigbgstatus, collision.multNTracksPV(), @@ -953,7 +952,7 @@ struct HfTreeCreatorLcToPKPi { fillEventProperties(collisions); - const int64_t candidatesSize = static_cast(candidates.size()); + const auto candidatesSize = static_cast(candidates.size()); reserveTables(candidatesSize, IsMc); int iCand{0}; @@ -983,7 +982,7 @@ struct HfTreeCreatorLcToPKPi { fillKFTable(candidate, collision, candFlag, functionSelection, sigbgstatus); } if (fillCandidateMcTable) { - float p{}, pt{}, svX{}, svY{}, svZ{}, pvX{}, pvY{}, pvZ{}, decayLength{}, lifetime{}; + float p{}, pt{}, svX{}, svY{}, svZ{}, pvX{}, pvY{}, pvZ{}, decayLength{}, decayTime{}; if (!isMcCandidateSignal) { p = UndefValueFloat; pt = UndefValueFloat; @@ -994,7 +993,7 @@ struct HfTreeCreatorLcToPKPi { pvY = UndefValueFloat; pvZ = UndefValueFloat; decayLength = UndefValueFloat; - lifetime = UndefValueFloat; + decayTime = UndefValueFloat; } else { const auto mcParticleProng0 = candidate.template prong0_as>().template mcParticle_as>(); const auto indexMother = RecoDecay::getMother(particles, mcParticleProng0, o2::constants::physics::Pdg::kLambdaCPlus, true); @@ -1002,8 +1001,6 @@ struct HfTreeCreatorLcToPKPi { const auto mcCollision = particleMother.template mcCollision_as(); p = particleMother.p(); pt = particleMother.pt(); - const float p2m = p / static_cast(MassLambdaCPlus); - const float gamma = std::sqrt(1 + p2m * p2m); // mother's particle Lorentz factor pvX = mcCollision.posX(); pvY = mcCollision.posY(); pvZ = mcCollision.posZ(); @@ -1011,11 +1008,11 @@ struct HfTreeCreatorLcToPKPi { svY = mcParticleProng0.vy(); svZ = mcParticleProng0.vz(); decayLength = static_cast(RecoDecay::distance(std::array{svX, svY, svZ}, std::array{pvX, pvY, pvZ})); - lifetime = mcParticleProng0.vt() * NanoToPico / gamma; // from ns to ps * from lab time to proper time + decayTime = decayLength * static_cast(MassLambdaCPlus) / LightSpeedCm2PS / p; } rowCandidateMC( p, pt, - svX, svY, svZ, decayLength, lifetime, + svX, svY, svZ, decayLength, decayTime, pvX, pvY, pvZ); } } @@ -1032,16 +1029,14 @@ struct HfTreeCreatorLcToPKPi { const auto mcDaughter0 = particle.template daughters_as>().begin(); const auto mcCollision = particle.template mcCollision_as(); const auto p = particle.p(); - const float p2m = p / static_cast(MassLambdaCPlus); - const float gamma = std::sqrt(1 + p2m * p2m); // mother's particle Lorentz factor const float pvX = mcCollision.posX(); const float pvY = mcCollision.posY(); const float pvZ = mcCollision.posZ(); const float svX = mcDaughter0.vx(); const float svY = mcDaughter0.vy(); const float svZ = mcDaughter0.vz(); - const float l = static_cast(RecoDecay::distance(std::array{svX, svY, svZ}, std::array{pvX, pvY, pvZ})); - const float t = mcDaughter0.vt() * NanoToPico / gamma; // from ns to ps * from lab time to proper time + const float decayLength = static_cast(RecoDecay::distance(std::array{svX, svY, svZ}, std::array{pvX, pvY, pvZ})); + const float decayTime = decayLength * static_cast(MassLambdaCPlus) / LightSpeedCm2PS / p; rowCandidateFullParticles( particle.pt(), particle.eta(), @@ -1050,7 +1045,7 @@ struct HfTreeCreatorLcToPKPi { particle.flagMcMatchGen(), particle.originMcGen(), p, - svX, svY, svZ, l, t, + svX, svY, svZ, decayLength, decayTime, pvX, pvY, pvZ); } }