diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index f826d3ff620..eb493b131c0 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -956,6 +956,40 @@ enum DecayChannelDToKKPi { DplusToK0starK // used to describe D+ in MC production for Ds analysis }; +// KF related properties +DECLARE_SOA_COLUMN(KfXPVError, kfXPVError, float); //! error of X coordinate of the event's primary vertex +DECLARE_SOA_COLUMN(KfYPVError, kfYPVError, float); //! error of Y coordinate of the event's primary vertex +DECLARE_SOA_COLUMN(KfZPVError, kfZPVError, float); //! error of Z coordinate of the event's primary vertex +DECLARE_SOA_COLUMN(KfXError, kfXError, float); //! error of candidate's decay point X coordinate from the KFParticle fit +DECLARE_SOA_COLUMN(KfYError, kfYError, float); //! error of candidate's decay point Y coordinate from the KFParticle fit +DECLARE_SOA_COLUMN(KfZError, kfZError, float); //! error of candidate's decay point Z coordinate from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassPKPi, kfMassPKPi, float); //! mass of the PKPi candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassPiKP, kfMassPiKP, float); //! mass of the PiKP candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassPiKPi, kfMassPiKPi, float); //! mass of the PiKPi candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassKKPi, kfMassKKPi, float); //! mass of the KKPi candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassPiKK, kfMassPiKK, float); //! mass of the PiKK candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassKPi, kfMassKPi, float); //! mass of the KPi pair from the KFParticle fit +DECLARE_SOA_COLUMN(KfMassPiK, kfMassPiK, float); //! mass of the PiK pair from the KFParticle fit +DECLARE_SOA_COLUMN(KfPx, kfPx, float); //! Px of the candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfPy, kfPy, float); //! Py of the candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfPz, kfPz, float); //! Pz of the candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfErrorPx, kfErrorPx, float); //! Px error of the candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfErrorPy, kfErrorPy, float); //! Py error of the candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfErrorPz, kfErrorPz, float); //! Pz error of the candidate from the KFParticle fit +DECLARE_SOA_COLUMN(KfChi2PrimProng0, kfChi2PrimProng0, float); //! chi2 primary of the first prong +DECLARE_SOA_COLUMN(KfChi2PrimProng1, kfChi2PrimProng1, float); //! chi2 primary of the second prong +DECLARE_SOA_COLUMN(KfChi2PrimProng2, kfChi2PrimProng2, float); //! chi2 primary of the third prong +DECLARE_SOA_COLUMN(KfDcaProng0Prong1, kfDcaProng0Prong1, float); //! DCA between first and second prongs +DECLARE_SOA_COLUMN(KfDcaProng0Prong2, kfDcaProng0Prong2, float); //! DCA between first and third prongs +DECLARE_SOA_COLUMN(KfDcaProng1Prong2, kfDcaProng1Prong2, float); //! DCA between second and third prongs +DECLARE_SOA_COLUMN(KfChi2GeoProng0Prong1, kfChi2GeoProng0Prong1, float); //! chi2 geo between first and second prongs +DECLARE_SOA_COLUMN(KfChi2GeoProng0Prong2, kfChi2GeoProng0Prong2, float); //! chi2 geo between first and third prongs +DECLARE_SOA_COLUMN(KfChi2GeoProng1Prong2, kfChi2GeoProng1Prong2, float); //! chi2 geo between second and third prongs +DECLARE_SOA_COLUMN(KfChi2Geo, kfChi2Geo, float); //! chi2 geo of the full candidate +DECLARE_SOA_COLUMN(KfChi2Topo, kfChi2Topo, float); //! chi2 topo of the full candidate (chi2prim of candidate to PV) +DECLARE_SOA_COLUMN(KfDecayLength, kfDecayLength, float); //! decay length +DECLARE_SOA_COLUMN(KfDecayLengthError, kfDecayLengthError, float); //! decay length error + } // namespace hf_cand_3prong // 3-prong decay candidate table @@ -1005,6 +1039,17 @@ DECLARE_SOA_EXTENDED_TABLE_USER(HfCand3ProngExt, HfCand3ProngBase, "HFCAND3PEXT" using HfCand3Prong = HfCand3ProngExt; +DECLARE_SOA_TABLE(HfCand3ProngKF, "AOD", "HFCAND3PKF", + hf_cand_3prong::KfXError, hf_cand_3prong::KfYError, hf_cand_3prong::KfZError, + hf_cand_3prong::KfXPVError, hf_cand_3prong::KfYPVError, hf_cand_3prong::KfZPVError, + hf_cand_3prong::KfMassPKPi, hf_cand_3prong::KfMassPiKP, hf_cand_3prong::KfMassPiKPi, hf_cand_3prong::KfMassKKPi, hf_cand_3prong::KfMassPiKK, hf_cand_3prong::KfMassKPi, hf_cand_3prong::KfMassPiK, + hf_cand_3prong::KfPx, hf_cand_3prong::KfPy, hf_cand_3prong::KfPz, + hf_cand_3prong::KfErrorPx, hf_cand_3prong::KfErrorPy, hf_cand_3prong::KfErrorPz, + hf_cand_3prong::KfChi2PrimProng0, hf_cand_3prong::KfChi2PrimProng1, hf_cand_3prong::KfChi2PrimProng2, + hf_cand_3prong::KfDcaProng1Prong2, hf_cand_3prong::KfDcaProng0Prong2, hf_cand_3prong::KfDcaProng0Prong1, + hf_cand_3prong::KfChi2GeoProng1Prong2, hf_cand_3prong::KfChi2GeoProng0Prong2, hf_cand_3prong::KfChi2GeoProng0Prong1, + hf_cand_3prong::KfChi2Geo, hf_cand_3prong::KfDecayLength, hf_cand_3prong::KfDecayLengthError, hf_cand_3prong::KfChi2Topo); + // table with results of reconstruction level MC matching DECLARE_SOA_TABLE(HfCand3ProngMcRec, "AOD", "HFCAND3PMCREC", //! hf_cand_3prong::FlagMcMatchRec, diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index a9d71f86da1..b039c0b6fb7 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -42,7 +42,7 @@ o2physics_add_dpl_workflow(candidate-creator-2prong o2physics_add_dpl_workflow(candidate-creator-3prong SOURCES candidateCreator3Prong.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter O2Physics::EventFilteringUtils + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(candidate-creator-b0 diff --git a/PWGHF/TableProducer/candidateCreator3Prong.cxx b/PWGHF/TableProducer/candidateCreator3Prong.cxx index d59cf2989b2..38329122281 100644 --- a/PWGHF/TableProducer/candidateCreator3Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator3Prong.cxx @@ -15,11 +15,21 @@ /// /// \author Vít Kučera , CERN +#ifndef HomogeneousField +#define HomogeneousField +#endif + #include #include #include #include +#include +#include +#include +#include +#include + #include #include "CommonConstants/PhysicsConstants.h" @@ -31,6 +41,7 @@ #include "ReconstructionDataFormats/DCA.h" #include "Common/Core/trackUtilities.h" +#include "Tools/KFparticle/KFUtilities.h" #include "PWGHF/Core/CentralityEstimation.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" @@ -52,6 +63,7 @@ using namespace o2::framework::expressions; /// Reconstruction of heavy-flavour 3-prong decay candidates struct HfCandidateCreator3Prong { Produces rowCandidateBase; + Produces rowCandidateKF; // vertexing Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; @@ -79,11 +91,20 @@ struct HfCandidateCreator3Prong { int runNumber{0}; float toMicrometers = 10000.; // from cm to µm + double massP{0.}; double massPi{0.}; double massK{0.}; + double massPKPi{0.}; + double massPiKP{0.}; double massPiKPi{0.}; + double massKKPi{0.}; + double massPiKK{0.}; + double massKPi{0.}; + double massPiK{0.}; double bz{0.}; + constexpr static float UndefValueFloat{-999.f}; + using FilteredHf3Prongs = soa::Filtered; using FilteredPvRefitHf3Prongs = soa::Filtered>; @@ -95,10 +116,13 @@ struct HfCandidateCreator3Prong { void init(InitContext const&) { - std::array processes = {doprocessPvRefit, doprocessNoPvRefit, - doprocessPvRefitCentFT0C, doprocessNoPvRefitCentFT0C, - doprocessPvRefitCentFT0M, doprocessNoPvRefitCentFT0M}; - if (std::accumulate(processes.begin(), processes.end(), 0) != 1) { + std::array doprocessDF{doprocessPvRefitWithDCAFitterN, doprocessNoPvRefitWithDCAFitterN, + doprocessPvRefitWithDCAFitterNCentFT0C, doprocessNoPvRefitWithDCAFitterNCentFT0C, + doprocessPvRefitWithDCAFitterNCentFT0M, doprocessNoPvRefitWithDCAFitterNCentFT0M}; + std::array doprocessKF{doprocessPvRefitWithKFParticle, doprocessNoPvRefitWithKFParticle, + doprocessPvRefitWithKFParticleCentFT0C, doprocessNoPvRefitWithKFParticleCentFT0C, + doprocessPvRefitWithKFParticleCentFT0M, doprocessNoPvRefitWithKFParticleCentFT0M}; + if ((std::accumulate(doprocessDF.begin(), doprocessDF.end(), 0) + std::accumulate(doprocessKF.begin(), doprocessKF.end(), 0)) != 1) { LOGP(fatal, "One and only one process function must be enabled at a time."); } @@ -108,13 +132,13 @@ struct HfCandidateCreator3Prong { LOGP(fatal, "At most one process function for collision monitoring can be enabled at a time."); } if (nProcessesCollisions == 1) { - if ((doprocessPvRefit || doprocessNoPvRefit) && !doprocessCollisions) { + if ((doprocessPvRefitWithDCAFitterN || doprocessNoPvRefitWithDCAFitterN || doprocessPvRefitWithKFParticle || doprocessNoPvRefitWithKFParticle) && !doprocessCollisions) { LOGP(fatal, "Process function for collision monitoring not correctly enabled. Did you enable \"processCollisions\"?"); } - if ((doprocessPvRefitCentFT0C || doprocessNoPvRefitCentFT0C) && !doprocessCollisionsCentFT0C) { + if ((doprocessPvRefitWithDCAFitterNCentFT0C || doprocessNoPvRefitWithDCAFitterNCentFT0C || doprocessPvRefitWithKFParticleCentFT0C || doprocessNoPvRefitWithKFParticleCentFT0C) && !doprocessCollisionsCentFT0C) { LOGP(fatal, "Process function for collision monitoring not correctly enabled. Did you enable \"processCollisionsCentFT0C\"?"); } - if ((doprocessPvRefitCentFT0M || doprocessNoPvRefitCentFT0M) && !doprocessCollisionsCentFT0M) { + if ((doprocessPvRefitWithDCAFitterNCentFT0M || doprocessNoPvRefitWithDCAFitterNCentFT0M || doprocessPvRefitWithKFParticleCentFT0M || doprocessNoPvRefitWithKFParticleCentFT0M) && !doprocessCollisionsCentFT0M) { LOGP(fatal, "Process function for collision monitoring not correctly enabled. Did you enable \"processCollisionsCentFT0M\"?"); } } @@ -125,7 +149,13 @@ struct HfCandidateCreator3Prong { } // histograms - registry.add("hMass3", "3-prong candidates;inv. mass (#pi K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1F, {{500, 1.6, 2.1}}}); + registry.add("hMass3PKPi", "3-prong candidates;inv. mass (pK#pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{1200, 1.8, 3.0}}}); + registry.add("hMass3PiKP", "3-prong candidates;inv. mass (#pi Kp) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{1200, 1.8, 3.0}}}); + registry.add("hMass3PiKPi", "3-prong candidates;inv. mass (#pi K#pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{600, 1.6, 2.2}}}); + registry.add("hMass3KKPi", "3-prong candidates;inv. mass (KK #pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{600, 1.7, 2.3}}}); + registry.add("hMass3PiKK", "3-prong candidates;inv. mass (#pi KK) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{600, 1.7, 2.3}}}); + registry.add("hMass2KPi", "2-prong pairs;inv. mass (K#pi) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{1200, 0.8, 2.0}}}); + registry.add("hMass2PiK", "2-prong pairs;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH1D, {{1200, 0.8, 2.0}}}); registry.add("hCovPVXX", "3-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", {HistType::kTH1F, {{100, 0., 1.e-4}}}); registry.add("hCovSVXX", "3-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", {HistType::kTH1F, {{100, 0., 0.2}}}); registry.add("hCovPVYY", "3-prong candidates;YY element of cov. matrix of prim. vtx. position (cm^{2});entries", {HistType::kTH1F, {{100, 0., 1.e-4}}}); @@ -139,6 +169,7 @@ struct HfCandidateCreator3Prong { hCandidates = registry.add("hCandidates", "candidates counter", {HistType::kTH1D, {axisCands}}); hfEvSel.addHistograms(registry); // collision monitoring + massP = MassProton; massPi = MassPiPlus; massK = MassKPlus; @@ -162,10 +193,10 @@ struct HfCandidateCreator3Prong { } template - void runCreator3Prong(Coll const&, - Cand const& rowsTrackIndexProng3, - aod::TracksWCovExtra const&, - aod::BCsWithTimestamps const& /*bcWithTimeStamps*/) + void runCreator3ProngWithDCAFitterN(Coll const&, + Cand const& rowsTrackIndexProng3, + aod::TracksWCovExtra const&, + aod::BCsWithTimestamps const& /*bcWithTimeStamps*/) { // loop over triplets of track indices for (const auto& rowTrackIndexProng3 : rowsTrackIndexProng3) { @@ -308,8 +339,239 @@ struct HfCandidateCreator3Prong { if (fillHistograms) { // calculate invariant mass auto arrayMomenta = std::array{pvec0, pvec1, pvec2}; - massPiKPi = RecoDecay::m(std::move(arrayMomenta), std::array{massPi, massK, massPi}); - registry.fill(HIST("hMass3"), massPiKPi); + massPKPi = RecoDecay::m(arrayMomenta, std::array{massP, massK, massPi}); + massPiKP = RecoDecay::m(arrayMomenta, std::array{massPi, massK, massP}); + massPiKPi = RecoDecay::m(arrayMomenta, std::array{massPi, massK, massPi}); + massKKPi = RecoDecay::m(arrayMomenta, std::array{massK, massK, massPi}); + massPiKK = RecoDecay::m(arrayMomenta, std::array{massPi, massK, massK}); + massKPi = RecoDecay::m(std::array{arrayMomenta.at(1), arrayMomenta.at(2)}, std::array{massK, massPi}); + massPiK = RecoDecay::m(std::array{arrayMomenta.at(0), arrayMomenta.at(1)}, std::array{massPi, massK}); + registry.fill(HIST("hMass3PiKPi"), massPiKPi); + registry.fill(HIST("hMass3PKPi"), massPKPi); + registry.fill(HIST("hMass3PiKP"), massPiKP); + registry.fill(HIST("hMass3KKPi"), massKKPi); + registry.fill(HIST("hMass3PiKK"), massPiKK); + registry.fill(HIST("hMass2KPi"), massKPi); + registry.fill(HIST("hMass2PiK"), massPiK); + } + } + } + + template + void runCreator3ProngWithKFParticle(Coll const&, + Cand const& rowsTrackIndexProng3, + aod::TracksWCovExtra const&, + aod::BCsWithTimestamps const& /*bcWithTimeStamps*/) + { + for (const auto& rowTrackIndexProng3 : rowsTrackIndexProng3) { + /// reject candidates in collisions not satisfying the event selections + auto collision = rowTrackIndexProng3.template collision_as(); + float centrality{-1.f}; + const auto rejectionMask = hfEvSel.getHfCollisionRejectionMask(collision, centrality, ccdb, registry); + if (rejectionMask != 0) { + /// at least one event selection not satisfied --> reject the candidate + continue; + } + + auto track0 = rowTrackIndexProng3.template prong0_as(); + auto track1 = rowTrackIndexProng3.template prong1_as(); + auto track2 = rowTrackIndexProng3.template prong2_as(); + + /// Set the magnetic field from ccdb. + /// The static instance of the propagator was already modified in the HFTrackIndexSkimCreator, + /// but this is not true when running on Run2 data/MC already converted into AO2Ds. + auto bc = collision.template bc_as(); + if (runNumber != bc.runNumber()) { + LOG(info) << ">>>>>>>>>>>> Current run number: " << runNumber; + initCCDB(bc, runNumber, ccdb, isRun2 ? ccdbPathGrp : ccdbPathGrpMag, nullptr, isRun2); + bz = o2::base::Propagator::Instance()->getNominalBz(); + LOG(info) << ">>>>>>>>>>>> Magnetic field: " << bz; + // df.setBz(bz); /// put it outside the 'if'! Otherwise we have a difference wrt bz Configurable (< 1 permille) in Run2 conv. data + // df.print(); + } + float covMatrixPV[6]; + + KFParticle::SetField(bz); + KFPVertex kfpVertex = createKFPVertexFromCollision(collision); + + if constexpr (doPvRefit) { + /// use PV refit + /// Using it in the rowCandidateBase all dynamic columns shall take it into account + // coordinates + kfpVertex.SetXYZ(rowTrackIndexProng3.pvRefitX(), rowTrackIndexProng3.pvRefitY(), rowTrackIndexProng3.pvRefitZ()); + // covariance matrix + kfpVertex.SetCovarianceMatrix(rowTrackIndexProng3.pvRefitSigmaX2(), rowTrackIndexProng3.pvRefitSigmaXY(), rowTrackIndexProng3.pvRefitSigmaY2(), rowTrackIndexProng3.pvRefitSigmaXZ(), rowTrackIndexProng3.pvRefitSigmaYZ(), rowTrackIndexProng3.pvRefitSigmaZ2()); + } + kfpVertex.GetCovarianceMatrix(covMatrixPV); + KFParticle KFPV(kfpVertex); + registry.fill(HIST("hCovPVXX"), covMatrixPV[0]); + registry.fill(HIST("hCovPVYY"), covMatrixPV[2]); + registry.fill(HIST("hCovPVXZ"), covMatrixPV[3]); + registry.fill(HIST("hCovPVZZ"), covMatrixPV[5]); + + KFPTrack kfpTrack0 = createKFPTrackFromTrack(track0); + KFPTrack kfpTrack1 = createKFPTrackFromTrack(track1); + KFPTrack kfpTrack2 = createKFPTrackFromTrack(track2); + + KFParticle kfFirstProton(kfpTrack0, kProton); + KFParticle kfFirstPion(kfpTrack0, kPiPlus); + KFParticle kfFirstKaon(kfpTrack0, kKPlus); + KFParticle kfSecondKaon(kfpTrack1, kKPlus); + KFParticle kfThirdProton(kfpTrack2, kProton); + KFParticle kfThirdPion(kfpTrack2, kPiPlus); + KFParticle kfThirdKaon(kfpTrack2, kKPlus); + + float impactParameter0XY = 0., errImpactParameter0XY = 0., impactParameter1XY = 0., errImpactParameter1XY = 0., impactParameter2XY = 0., errImpactParameter2XY = 0.; + if (!kfFirstProton.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + registry.fill(HIST("hDcaXYProngs"), track0.pt(), impactParameter0XY * toMicrometers); + registry.fill(HIST("hDcaZProngs"), track0.pt(), std::sqrt(kfFirstProton.GetDistanceFromVertex(KFPV) * kfFirstProton.GetDistanceFromVertex(KFPV) - impactParameter0XY * impactParameter0XY) * toMicrometers); + } else { + registry.fill(HIST("hDcaXYProngs"), track0.pt(), UndefValueFloat); + registry.fill(HIST("hDcaZProngs"), track0.pt(), UndefValueFloat); + } + if (!kfSecondKaon.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + registry.fill(HIST("hDcaXYProngs"), track1.pt(), impactParameter1XY * toMicrometers); + registry.fill(HIST("hDcaZProngs"), track1.pt(), std::sqrt(kfSecondKaon.GetDistanceFromVertex(KFPV) * kfSecondKaon.GetDistanceFromVertex(KFPV) - impactParameter1XY * impactParameter1XY) * toMicrometers); + } else { + registry.fill(HIST("hDcaXYProngs"), track1.pt(), UndefValueFloat); + registry.fill(HIST("hDcaZProngs"), track1.pt(), UndefValueFloat); + } + if (!kfThirdProton.GetDistanceFromVertexXY(KFPV, impactParameter2XY, errImpactParameter2XY)) { + registry.fill(HIST("hDcaXYProngs"), track2.pt(), impactParameter2XY * toMicrometers); + registry.fill(HIST("hDcaZProngs"), track2.pt(), std::sqrt(kfThirdProton.GetDistanceFromVertex(KFPV) * kfThirdProton.GetDistanceFromVertex(KFPV) - impactParameter2XY * impactParameter2XY) * toMicrometers); + } else { + registry.fill(HIST("hDcaXYProngs"), track2.pt(), UndefValueFloat); + registry.fill(HIST("hDcaZProngs"), track2.pt(), UndefValueFloat); + } + + auto [impactParameter0Z, errImpactParameter0Z] = kfCalculateImpactParameterZ(kfFirstProton, KFPV); + auto [impactParameter1Z, errImpactParameter1Z] = kfCalculateImpactParameterZ(kfSecondKaon, KFPV); + auto [impactParameter2Z, errImpactParameter2Z] = kfCalculateImpactParameterZ(kfThirdProton, KFPV); + + const float chi2primFirst = kfCalculateChi2ToPrimaryVertex(kfFirstProton, KFPV); + const float chi2primSecond = kfCalculateChi2ToPrimaryVertex(kfSecondKaon, KFPV); + const float chi2primThird = kfCalculateChi2ToPrimaryVertex(kfThirdPion, KFPV); + + const float dcaSecondThird = kfCalculateDistanceBetweenParticles(kfSecondKaon, kfThirdPion); + const float dcaFirstThird = kfCalculateDistanceBetweenParticles(kfFirstProton, kfThirdPion); + const float dcaFirstSecond = kfCalculateDistanceBetweenParticles(kfFirstProton, kfSecondKaon); + + const float chi2geoSecondThird = kfCalculateChi2geoBetweenParticles(kfSecondKaon, kfThirdPion); + const float chi2geoFirstThird = kfCalculateChi2geoBetweenParticles(kfFirstProton, kfThirdPion); + const float chi2geoFirstSecond = kfCalculateChi2geoBetweenParticles(kfFirstProton, kfSecondKaon); + + // Λc± → p± K∓ π±, Ξc± → p± K∓ π± + KFParticle kfCandPKPi; + const KFParticle* kfDaughtersPKPi[3] = {&kfFirstProton, &kfSecondKaon, &kfThirdPion}; + kfCandPKPi.SetConstructMethod(2); + kfCandPKPi.Construct(kfDaughtersPKPi, 3); + KFParticle kfCandPiKP; + const KFParticle* kfDaughtersPiKP[3] = {&kfFirstPion, &kfSecondKaon, &kfThirdProton}; + kfCandPiKP.SetConstructMethod(2); + kfCandPiKP.Construct(kfDaughtersPiKP, 3); + + // D± → π± K∓ π± + KFParticle kfCandPiKPi; + const KFParticle* kfDaughtersPiKPi[3] = {&kfFirstPion, &kfSecondKaon, &kfThirdPion}; + kfCandPiKPi.SetConstructMethod(2); + kfCandPiKPi.Construct(kfDaughtersPiKPi, 3); + + // Ds± → K± K∓ π± + KFParticle kfCandKKPi; + const KFParticle* kfDaughtersKKPi[3] = {&kfFirstKaon, &kfSecondKaon, &kfThirdPion}; + kfCandKKPi.SetConstructMethod(2); + kfCandKKPi.Construct(kfDaughtersKKPi, 3); + KFParticle kfCandPiKK; + const KFParticle* kfDaughtersPiKK[3] = {&kfFirstPion, &kfSecondKaon, &kfThirdKaon}; + kfCandPiKK.SetConstructMethod(2); + kfCandPiKK.Construct(kfDaughtersPiKK, 3); + + KFParticle kfPairKPi; + const KFParticle* kfDaughtersKPi[3] = {&kfSecondKaon, &kfThirdPion}; + kfPairKPi.SetConstructMethod(2); + kfPairKPi.Construct(kfDaughtersKPi, 2); + + KFParticle kfPairPiK; + const KFParticle* kfDaughtersPiK[3] = {&kfFirstPion, &kfSecondKaon}; + kfPairPiK.SetConstructMethod(2); + kfPairPiK.Construct(kfDaughtersPiK, 2); + + const float massPKPi = kfCandPKPi.GetMass(); + const float massPiKP = kfCandPiKP.GetMass(); + const float massPiKPi = kfCandPiKPi.GetMass(); + const float massKKPi = kfCandKKPi.GetMass(); + const float massPiKK = kfCandPiKK.GetMass(); + const float massKPi = kfPairKPi.GetMass(); + const float massPiK = kfPairPiK.GetMass(); + + const float chi2geo = kfCandPKPi.Chi2() / kfCandPKPi.NDF(); + const float chi2topo = kfCalculateChi2ToPrimaryVertex(kfCandPKPi, KFPV); + const std::pair ldl = kfCalculateLdL(kfCandPKPi, KFPV); + + std::array pProng0 = kfCalculateProngMomentumInSecondaryVertex(kfFirstProton, kfCandPiKP); + std::array pProng1 = kfCalculateProngMomentumInSecondaryVertex(kfSecondKaon, kfCandPiKP); + std::array pProng2 = kfCalculateProngMomentumInSecondaryVertex(kfThirdPion, kfCandPiKP); + + registry.fill(HIST("hCovSVXX"), kfCandPKPi.Covariance(0, 0)); + registry.fill(HIST("hCovSVYY"), kfCandPKPi.Covariance(1, 1)); + registry.fill(HIST("hCovSVXZ"), kfCandPKPi.Covariance(2, 0)); + registry.fill(HIST("hCovSVZZ"), kfCandPKPi.Covariance(2, 2)); + + auto covMatrixSV = kfCandPKPi.CovarianceMatrix(); + double phi, theta; + getPointDirection(std::array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, std::array{kfCandPKPi.GetX(), kfCandPKPi.GetY(), kfCandPKPi.GetZ()}, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixSV, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixSV, phi, 0.)); + + auto indexCollision = collision.globalIndex(); + uint8_t bitmapProngsContributorsPV = 0; + if (indexCollision == track0.collisionId() && track0.isPVContributor()) { + SETBIT(bitmapProngsContributorsPV, 0); + } + if (indexCollision == track1.collisionId() && track1.isPVContributor()) { + SETBIT(bitmapProngsContributorsPV, 1); + } + if (indexCollision == track2.collisionId() && track2.isPVContributor()) { + SETBIT(bitmapProngsContributorsPV, 2); + } + uint8_t nProngsContributorsPV = hf_trkcandsel::countOnesInBinary(bitmapProngsContributorsPV); + + // fill candidate table rows + rowCandidateBase(indexCollision, + KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), + kfCandPKPi.GetX(), kfCandPKPi.GetY(), kfCandPKPi.GetZ(), + errorDecayLength, errorDecayLengthXY, + kfCandPKPi.GetChi2() / kfCandPKPi.GetNDF(), + pProng0.at(0), pProng0.at(1), pProng0.at(2), + pProng1.at(0), pProng1.at(1), pProng1.at(2), + pProng2.at(0), pProng2.at(1), pProng2.at(2), + impactParameter0XY, impactParameter1XY, impactParameter2XY, + errImpactParameter0XY, errImpactParameter1XY, errImpactParameter2XY, + impactParameter0Z, impactParameter1Z, impactParameter2Z, + errImpactParameter0Z, errImpactParameter1Z, errImpactParameter2Z, + rowTrackIndexProng3.prong0Id(), rowTrackIndexProng3.prong1Id(), rowTrackIndexProng3.prong2Id(), nProngsContributorsPV, bitmapProngsContributorsPV, + rowTrackIndexProng3.hfflag()); + + // fill KF info + rowCandidateKF(kfCandPKPi.GetErrX(), kfCandPKPi.GetErrY(), kfCandPKPi.GetErrZ(), + std::sqrt(KFPV.Covariance(0, 0)), std::sqrt(KFPV.Covariance(1, 1)), std::sqrt(KFPV.Covariance(2, 2)), + massPKPi, massPiKP, massPiKPi, massKKPi, massPiKK, massKPi, massPiK, + kfCandPKPi.GetPx(), kfCandPKPi.GetPy(), kfCandPKPi.GetPz(), + kfCandPKPi.GetErrPx(), kfCandPKPi.GetErrPy(), kfCandPKPi.GetErrPz(), + chi2primFirst, chi2primSecond, chi2primThird, + dcaSecondThird, dcaFirstThird, dcaFirstSecond, + chi2geoSecondThird, chi2geoFirstThird, chi2geoFirstSecond, + chi2geo, ldl.first, ldl.second, chi2topo); + + // fill histograms + if (fillHistograms) { + registry.fill(HIST("hMass3PiKPi"), massPiKPi); + registry.fill(HIST("hMass3PKPi"), massPKPi); + registry.fill(HIST("hMass3PiKP"), massPiKP); + registry.fill(HIST("hMass3KKPi"), massKKPi); + registry.fill(HIST("hMass3PiKK"), massPiKK); + registry.fill(HIST("hMass2KPi"), massKPi); + registry.fill(HIST("hMass2PiK"), massPiK); } } } @@ -320,25 +582,45 @@ struct HfCandidateCreator3Prong { /// /// /////////////////////////////////// - /// @brief process function w/ PV refit and w/o centrality selections - void processPvRefit(soa::Join const& collisions, - FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, - aod::TracksWCovExtra const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + /// @brief process function using DCA fitter w/ PV refit and w/o centrality selections + void processPvRefitWithDCAFitterN(soa::Join const& collisions, + FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } + PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitWithDCAFitterN, "Run candidate creator using DCA fitter with PV refit and w/o centrality selections", false); + + /// @brief process function using DCA fitter w/o PV refit and w/o centrality selections + void processNoPvRefitWithDCAFitterN(soa::Join const& collisions, + FilteredHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } + PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithDCAFitterN, "Run candidate creator using DCA fitter without PV refit and w/o centrality selections", true); + + /// @brief process function using KFParticle package w/ PV refit and w/o centrality selections + void processPvRefitWithKFParticle(soa::Join const& collisions, + FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3Prong(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + runCreator3ProngWithKFParticle(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefit, "Run candidate creator with PV refit and w/o centrality selections", false); + PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitWithKFParticle, "Run candidate creator using KFParticle package with PV refit and w/o centrality selections", false); - /// @brief process function w/o PV refit and w/o centrality selections - void processNoPvRefit(soa::Join const& collisions, - FilteredHf3Prongs const& rowsTrackIndexProng3, - aod::TracksWCovExtra const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + /// @brief process function using KFParticle package w/o PV refit and w/o centrality selections + void processNoPvRefitWithKFParticle(soa::Join const& collisions, + FilteredHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3Prong(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + runCreator3ProngWithKFParticle(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefit, "Run candidate creator without PV refit and w/o centrality selections", true); + PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithKFParticle, "Run candidate creator using KFParticle package without PV refit and w/o centrality selections", false); ///////////////////////////////////////////// /// /// @@ -346,25 +628,45 @@ struct HfCandidateCreator3Prong { /// /// ///////////////////////////////////////////// - /// @brief process function w/ PV refit and w/ centrality selection on FT0C - void processPvRefitCentFT0C(soa::Join const& collisions, - FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, - aod::TracksWCovExtra const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + /// @brief process function using DCA fitter w/ PV refit and w/ centrality selection on FT0C + void processPvRefitWithDCAFitterNCentFT0C(soa::Join const& collisions, + FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3Prong(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitCentFT0C, "Run candidate creator with PV refit and w/ centrality selection on FT0C", false); + PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitWithDCAFitterNCentFT0C, "Run candidate creator using DCA fitter with PV refit and w/ centrality selection on FT0C", false); - /// @brief process function w/o PV refit and w/ centrality selection on FT0C - void processNoPvRefitCentFT0C(soa::Join const& collisions, - FilteredHf3Prongs const& rowsTrackIndexProng3, - aod::TracksWCovExtra const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + /// @brief process function using DCA fitter w/o PV refit and w/ centrality selection on FT0C + void processNoPvRefitWithDCAFitterNCentFT0C(soa::Join const& collisions, + FilteredHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3Prong(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitCentFT0C, "Run candidate creator without PV refit and w/ centrality selection on FT0C", false); + PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithDCAFitterNCentFT0C, "Run candidate creator using DCA fitter without PV refit and w/ centrality selection on FT0C", false); + + /// @brief process function using KFParticle package w/ PV refit and w/ centrality selection on FT0C + void processPvRefitWithKFParticleCentFT0C(soa::Join const& collisions, + FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator3ProngWithKFParticle(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } + PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitWithKFParticleCentFT0C, "Run candidate creator using KFParticle package with PV refit and w/ centrality selection on FT0C", false); + + /// @brief process function using KFParticle package w/o PV refit and w/ centrality selection on FT0C + void processNoPvRefitWithKFParticleCentFT0C(soa::Join const& collisions, + FilteredHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator3ProngWithKFParticle(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } + PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithKFParticleCentFT0C, "Run candidate creator using KFParticle package without PV refit and w/ centrality selection on FT0C", false); ///////////////////////////////////////////// /// /// @@ -372,25 +674,45 @@ struct HfCandidateCreator3Prong { /// /// ///////////////////////////////////////////// - /// @brief process function w/ PV refit and w/ centrality selection on FT0M - void processPvRefitCentFT0M(soa::Join const& collisions, - FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, - aod::TracksWCovExtra const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + /// @brief process function using DCA fitter w/ PV refit and w/ centrality selection on FT0M + void processPvRefitWithDCAFitterNCentFT0M(soa::Join const& collisions, + FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } + PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitWithDCAFitterNCentFT0M, "Run candidate creator using DCA fitter with PV refit and w/ centrality selection on FT0M", false); + + /// @brief process function using DCA fitter w/o PV refit and w/ centrality selection on FT0M + void processNoPvRefitWithDCAFitterNCentFT0M(soa::Join const& collisions, + FilteredHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator3ProngWithDCAFitterN(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + } + PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithDCAFitterNCentFT0M, "Run candidate creator using DCA fitter without PV refit and w/ centrality selection on FT0M", false); + + /// @brief process function using KFParticle package w/ PV refit and w/ centrality selection on FT0M + void processPvRefitWithKFParticleCentFT0M(soa::Join const& collisions, + FilteredPvRefitHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3Prong(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + runCreator3ProngWithKFParticle(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitCentFT0M, "Run candidate creator with PV refit and w/ centrality selection on FT0M", false); + PROCESS_SWITCH(HfCandidateCreator3Prong, processPvRefitWithKFParticleCentFT0M, "Run candidate creator using KFParticle package with PV refit and w/ centrality selection on FT0M", false); - /// @brief process function w/o PV refit and w/ centrality selection on FT0M - void processNoPvRefitCentFT0M(soa::Join const& collisions, - FilteredHf3Prongs const& rowsTrackIndexProng3, - aod::TracksWCovExtra const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + /// @brief process function using KFParticle package w/o PV refit and w/ centrality selection on FT0M + void processNoPvRefitWithKFParticleCentFT0M(soa::Join const& collisions, + FilteredHf3Prongs const& rowsTrackIndexProng3, + aod::TracksWCovExtra const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator3Prong(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); + runCreator3ProngWithKFParticle(collisions, rowsTrackIndexProng3, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitCentFT0M, "Run candidate creator without PV refit and w/ centrality selection on FT0M", false); + PROCESS_SWITCH(HfCandidateCreator3Prong, processNoPvRefitWithKFParticleCentFT0M, "Run candidate creator using KFParticle package without PV refit and w/ centrality selection on FT0M", false); /////////////////////////////////////////////////////////// /// /// diff --git a/PWGHF/TableProducer/candidateSelectorLc.cxx b/PWGHF/TableProducer/candidateSelectorLc.cxx index 818d65bb269..eacc06a872c 100644 --- a/PWGHF/TableProducer/candidateSelectorLc.cxx +++ b/PWGHF/TableProducer/candidateSelectorLc.cxx @@ -108,7 +108,7 @@ struct HfCandidateSelectorLc { void init(InitContext const&) { - std::array processes = {doprocessNoBayesPid, doprocessBayesPid}; + std::array processes = {doprocessNoBayesPidWithDCAFitterN, doprocessBayesPidWithDCAFitterN, doprocessNoBayesPidWithKFParticle, doprocessBayesPidWithKFParticle}; if (std::accumulate(processes.begin(), processes.end(), 0) != 1) { LOGP(fatal, "One and only one process function must be enabled at a time."); } @@ -232,7 +232,7 @@ struct HfCandidateSelectorLc { /// \param trackPion is the track with the pion hypothesis /// \param trackKaon is the track with the kaon hypothesis /// \return true if candidate passes all cuts for the given Conjugate - template + template bool selectionTopolConjugate(const T1& candidate, const T2& trackProton, const T2& trackKaon, const T2& trackPion) { @@ -248,32 +248,34 @@ struct HfCandidateSelectorLc { } // cut on Lc->pKpi, piKp mass values - if (trackProton.globalIndex() == candidate.prong0Id()) { - if (std::abs(hfHelper.invMassLcToPKPi(candidate) - o2::constants::physics::MassLambdaCPlus) > cuts->get(pTBin, "m")) { - return false; + /// cut on the Kpi pair invariant mass, to study Lc->pK*(->Kpi) + float massLc, massKPi; + if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) { + if (trackProton.globalIndex() == candidate.prong0Id()) { + massLc = hfHelper.invMassLcToPKPi(candidate); + massKPi = hfHelper.invMassKPiPairLcToPKPi(candidate); + } else { + massLc = hfHelper.invMassLcToPiKP(candidate); + massKPi = hfHelper.invMassKPiPairLcToPiKP(candidate); } - } else { - if (std::abs(hfHelper.invMassLcToPiKP(candidate) - o2::constants::physics::MassLambdaCPlus) > cuts->get(pTBin, "m")) { - return false; + } else if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) { + if (trackProton.globalIndex() == candidate.prong0Id()) { + massLc = candidate.kfMassPKPi(); + massKPi = candidate.kfMassKPi(); + } else { + massLc = candidate.kfMassPiKP(); + massKPi = candidate.kfMassPiK(); } } + if (std::abs(massLc - o2::constants::physics::MassLambdaCPlus) > cuts->get(pTBin, "m")) { + return false; + } + /// cut on the Kpi pair invariant mass, to study Lc->pK*(->Kpi) const double cutMassKPi = cuts->get(pTBin, "mass (Kpi)"); - if (cutMassKPi > 0) { - if (trackProton.globalIndex() == candidate.prong0Id()) { - // inspecting the pKpi hypothesis - // K: prong1, pi: prong 2 - if (std::abs(hfHelper.invMassKPiPairLcToPKPi(candidate) - massK0Star892) > cutMassKPi) { - return false; - } - } else { - // inspecting the piKp hypothesis - // K: prong1, pi: prong 0 - if (std::abs(hfHelper.invMassKPiPairLcToPiKP(candidate) - massK0Star892) > cutMassKPi) { - return false; - } - } + if (cutMassKPi > 0 && std::abs(massKPi - massK0Star892) > cutMassKPi) { + return false; } return true; @@ -306,10 +308,11 @@ struct HfCandidateSelectorLc { } /// \brief function to apply Lc selections + /// \param reconstructionType is the reconstruction type (DCAFitterN or KFParticle) /// \param candidates Lc candidate table /// \param tracks track table - template - void runSelectLc(aod::HfCand3Prong const& candidates, TTracks const&) + template + void runSelectLc(CandType const& candidates, TTracks const&) { // looping over 3-prong candidates for (const auto& candidate : candidates) { @@ -364,8 +367,8 @@ struct HfCandidateSelectorLc { } // conjugate-dependent topological selection for Lc - bool topolLcToPKPi = selectionTopolConjugate(candidate, trackPos1, trackNeg, trackPos2); - bool topolLcToPiKP = selectionTopolConjugate(candidate, trackPos2, trackNeg, trackPos1); + bool topolLcToPKPi = selectionTopolConjugate(candidate, trackPos1, trackNeg, trackPos2); + bool topolLcToPiKP = selectionTopolConjugate(candidate, trackPos2, trackNeg, trackPos1); if (!topolLcToPKPi && !topolLcToPiKP) { hfSelLcCandidate(statusLcToPKPi, statusLcToPiKP); @@ -481,25 +484,45 @@ struct HfCandidateSelectorLc { } } - /// \brief process function w/o Bayes PID + /// \brief process function w/o Bayes PID with DCAFitterN + /// \param candidates Lc candidate table + /// \param tracks track table + void processNoBayesPidWithDCAFitterN(aod::HfCand3Prong const& candidates, + TracksSel const& tracks) + { + runSelectLc(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorLc, processNoBayesPidWithDCAFitterN, "Process Lc selection w/o Bayes PID with DCAFitterN", true); + + /// \brief process function with Bayes PID with DCAFitterN + /// \param candidates Lc candidate table + /// \param tracks track table with Bayes PID information + void processBayesPidWithDCAFitterN(aod::HfCand3Prong const& candidates, + TracksSelBayesPid const& tracks) + { + runSelectLc(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorLc, processBayesPidWithDCAFitterN, "Process Lc selection with Bayes PID with DCAFitterN", false); + + /// \brief process function w/o Bayes PID with KFParticle /// \param candidates Lc candidate table /// \param tracks track table - void processNoBayesPid(aod::HfCand3Prong const& candidates, - TracksSel const& tracks) + void processNoBayesPidWithKFParticle(soa::Join const& candidates, + TracksSel const& tracks) { - runSelectLc(candidates, tracks); + runSelectLc(candidates, tracks); } - PROCESS_SWITCH(HfCandidateSelectorLc, processNoBayesPid, "Process Lc selection w/o Bayes PID", true); + PROCESS_SWITCH(HfCandidateSelectorLc, processNoBayesPidWithKFParticle, "Process Lc selection w/o Bayes PID with KFParticle", false); - /// \brief process function with Bayes PID + /// \brief process function with Bayes PID with KFParticle /// \param candidates Lc candidate table /// \param tracks track table with Bayes PID information - void processBayesPid(aod::HfCand3Prong const& candidates, - TracksSelBayesPid const& tracks) + void processBayesPidWithKFParticle(soa::Join const& candidates, + TracksSelBayesPid const& tracks) { - runSelectLc(candidates, tracks); + runSelectLc(candidates, tracks); } - PROCESS_SWITCH(HfCandidateSelectorLc, processBayesPid, "Process Lc selection with Bayes PID", false); + PROCESS_SWITCH(HfCandidateSelectorLc, processBayesPidWithKFParticle, "Process Lc selection with Bayes PID with KFParticle", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx b/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx index 07f869382a1..8dfd9f517ee 100644 --- a/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorLcToPKPi.cxx @@ -30,6 +30,7 @@ using namespace o2; using namespace o2::framework; +using namespace o2::constants::physics; namespace o2::aod { @@ -96,31 +97,87 @@ DECLARE_SOA_COLUMN(CentFDDM, centFDDM, float); DECLARE_SOA_COLUMN(MultZeqNTracksPV, multZeqNTracksPV, float); } // namespace full +namespace kf +{ +DECLARE_SOA_COLUMN(X, x, float); //! decay vertex X coordinate +DECLARE_SOA_COLUMN(Y, y, float); //! decay vertex Y coordinate +DECLARE_SOA_COLUMN(Z, z, float); //! decay vertex Z coordinate +DECLARE_SOA_COLUMN(ErrX, errX, float); //! decay vertex X coordinate error +DECLARE_SOA_COLUMN(ErrY, errY, float); //! decay vertex Y coordinate error +DECLARE_SOA_COLUMN(ErrZ, errZ, float); //! decay vertex Z coordinate error +DECLARE_SOA_COLUMN(ErrPVX, errPVX, float); //! event vertex X coordinate error +DECLARE_SOA_COLUMN(ErrPVY, errPVY, float); //! event vertex Y coordinate error +DECLARE_SOA_COLUMN(ErrPVZ, errPVZ, float); //! event vertex Z coordinate error +DECLARE_SOA_COLUMN(Chi2PrimProton, chi2PrimProton, float); //! Chi2 of prong's approach to the PV +DECLARE_SOA_COLUMN(Chi2PrimKaon, chi2PrimKaon, float); //! Chi2 of prong's approach to the PV +DECLARE_SOA_COLUMN(Chi2PrimPion, chi2PrimPion, float); //! Chi2 of prong's approach to the PV +DECLARE_SOA_COLUMN(DcaProtonKaon, dcaProtonKaon, float); //! Distance of closest approach between 2 prongs, cm +DECLARE_SOA_COLUMN(DcaProtonPion, dcaProtonPion, float); //! Distance of closest approach between 2 prongs, cm +DECLARE_SOA_COLUMN(DcaPionKaon, dcaPionKaon, float); //! Distance of closest approach between 2 prongs, cm +DECLARE_SOA_COLUMN(Chi2GeoProtonKaon, chi2GeoProtonKaon, float); //! Chi2 of two prongs' approach to each other +DECLARE_SOA_COLUMN(Chi2GeoProtonPion, chi2GeoProtonPion, float); //! Chi2 of two prongs' approach to each other +DECLARE_SOA_COLUMN(Chi2GeoPionKaon, chi2GeoPionKaon, float); //! Chi2 of two prongs' approach to each other +DECLARE_SOA_COLUMN(Chi2Geo, chi2Geo, float); //! chi2 geo of the full candidate +DECLARE_SOA_COLUMN(Chi2Topo, chi2Topo, float); //! chi2 topo of the full candidate (chi2prim of candidate to PV) +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(MassInv, massInv, float); //! invariant mass +DECLARE_SOA_COLUMN(P, p, float); //! momentum +DECLARE_SOA_COLUMN(Pt, pt, float); //! transverse momentum +DECLARE_SOA_COLUMN(ErrP, errP, float); //! momentum error +DECLARE_SOA_COLUMN(ErrPt, errPt, float); //! transverse momentum error +DECLARE_SOA_COLUMN(IsSelected, isSelected, int); //! flag whether candidate was selected in candidateSelectorLc task +DECLARE_SOA_COLUMN(SigBgStatus, sigBgStatus, int); //! 0 bg, 1 prompt, 2 non-prompt, 3 wrong order of prongs, -1 default value (impossible, should not be the case), -999 for data +} // namespace kf + +namespace mc_match +{ +DECLARE_SOA_COLUMN(P, p, float); //! Momentum, GeV/c +DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum, GeV/c +DECLARE_SOA_COLUMN(XDecay, xDecay, float); //! Secondary (decay) vertex X coordinate, cm +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(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 +} // namespace mc_match + +DECLARE_SOA_TABLE(HfCandLcMCs, "AOD", "HFCANDLCMC", + mc_match::P, mc_match::Pt, + mc_match::XDecay, mc_match::YDecay, mc_match::ZDecay, mc_match::LDecay, + mc_match::TDecay, + mc_match::XEvent, mc_match::YEvent, mc_match::ZEvent) + +DECLARE_SOA_TABLE(HfCandLcKFs, "AOD", "HFCANDLCKF", + kf::X, kf::Y, kf::Z, kf::ErrX, kf::ErrY, kf::ErrZ, + kf::ErrPVX, kf::ErrPVY, kf::ErrPVZ, + kf::Chi2PrimProton, kf::Chi2PrimKaon, kf::Chi2PrimPion, + kf::DcaProtonKaon, kf::DcaProtonPion, kf::DcaPionKaon, + kf::Chi2GeoProtonKaon, kf::Chi2GeoProtonPion, kf::Chi2GeoPionKaon, + kf::Chi2Geo, kf::Chi2Topo, kf::DecayLength, kf::DecayLengthError, kf::DecayLengthNormalised, kf::T, kf::ErrT, + kf::MassInv, kf::P, kf::Pt, kf::ErrP, kf::ErrPt, + kf::IsSelected, kf::SigBgStatus); + DECLARE_SOA_TABLE(HfCandLcLites, "AOD", "HFCANDLCLITE", collision::PosX, collision::PosY, collision::PosZ, hf_cand::NProngsContributorsPV, hf_cand::BitmapProngsContributorsPV, - // hf_cand::ErrorDecayLength, - // hf_cand::ErrorDecayLengthXY, hf_cand::Chi2PCA, full::DecayLength, full::DecayLengthXY, - // full::DecayLengthNormalised, - // full::DecayLengthXYNormalised, - // full::ImpactParameterNormalised0, full::PtProng0, - // full::ImpactParameterNormalised1, full::PtProng1, - // full::ImpactParameterNormalised2, full::PtProng2, hf_cand::ImpactParameter0, hf_cand::ImpactParameter1, hf_cand::ImpactParameter2, - // hf_cand::ErrorImpactParameter0, - // hf_cand::ErrorImpactParameter1, - // hf_cand::ErrorImpactParameter2, full::NSigTpcPi0, full::NSigTpcPr0, full::NSigTofPi0, @@ -252,34 +309,96 @@ DECLARE_SOA_TABLE(HfCandLcFullPs, "AOD", "HFCANDLCFULLP", full::Y, full::FlagMc, full::OriginMcGen); + } // namespace o2::aod /// Writes the full information in an output TTree struct HfTreeCreatorLcToPKPi { Produces rowCandidateFull; Produces rowCandidateLite; + Produces rowCandidateKF; + Produces rowCandidateMC; Produces rowCollisionId; Produces rowCandidateFullEvents; Produces rowCandidateFullParticles; + Configurable selectionFlagLc{"selectionFlagLc", 1, "Selection Flag for Lc"}; Configurable fillCandidateLiteTable{"fillCandidateLiteTable", false, "Switch to fill lite table with candidate properties"}; Configurable fillCollIdTable{"fillCollIdTable", false, "Fill a single-column table with collision index"}; + Configurable fillCandidateMcTable{"fillCandidateMcTable", false, "Switch to fill a table with MC particles matched to candidates"}; Configurable keepOnlySignalMc{"keepOnlySignalMc", false, "Fill MC tree only with signal candidates"}; Configurable keepOnlyBkg{"keepOnlyBkg", false, "Fill MC tree only with background candidates"}; Configurable downSampleBkgFactor{"downSampleBkgFactor", 1., "Fraction of candidates to store in the tree"}; Configurable downSampleBkgPtMax{"downSampleBkgPtMax", 100.f, "Max. pt for background downsampling"}; + constexpr static float UndefValueFloat = -999.f; + constexpr static int UndefValueInt = -999; + constexpr static float NanoToPico = 1000.f; + HfHelper hfHelper; using TracksWPid = soa::Join; using Cents = soa::Join; + // 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 + Default = -1 // impossible, should not be the case, to catch logical error if any + }; + + /// \brief function which determines if the candidate corresponds to MC-particle or belongs to a combinatorial background + /// \param candidate candidate to be checked for being signal or background + /// \param CandFlag 0 for PKPi hypothesis and 1 for PiKP hypothesis + /// \return SigBgStatus enum with value encoding MC status of the candidate + template + SigBgStatus determineSignalBgStatus(const CandType& candidate, int CandFlag) + { + const int flag = candidate.flagMcMatchRec(); + const int origin = candidate.originMcRec(); + const int swapped = candidate.isCandidateSwapped(); + + SigBgStatus status{Default}; + + if (std::abs(flag) == (1 << o2::aod::hf_cand_3prong::DecayType::LcToPKPi)) { + if (swapped == 0) { + if (CandFlag == 0) { + if (origin == RecoDecay::OriginType::Prompt) + status = Prompt; + else if (origin == RecoDecay::OriginType::NonPrompt) + status = NonPrompt; + } else { + status = WrongOrder; + } + } else { + if (CandFlag == 1) { + if (origin == RecoDecay::OriginType::Prompt) + status = Prompt; + else if (origin == RecoDecay::OriginType::NonPrompt) + status = NonPrompt; + } else { + status = WrongOrder; + } + } + } else { + status = Background; + } + + return status; + } + void init(InitContext const&) { - std::array processes = {doprocessDataNoCentrality, doprocessDataWithCentrality, doprocessMcNoCentrality, doprocessMcWithCentrality}; + std::array processes = {doprocessDataNoCentralityWithDCAFitterN, doprocessDataWithCentralityWithDCAFitterN, doprocessDataNoCentralityWithKFParticle, doprocessDataWithCentralityWithKFParticle, + doprocessMcNoCentralityWithDCAFitterN, doprocessMcWithCentralityWithDCAFitterN, doprocessMcNoCentralityWithKFParticle, doprocessMcWithCentralityWithKFParticle}; if (std::accumulate(processes.begin(), processes.end(), 0) != 1) { LOGP(fatal, "One and only one process function must be enabled at a time."); } + if (std::accumulate(processes.begin(), processes.begin() + 4, 0) && fillCandidateMcTable) { + LOGP(fatal, "fillCandidateMcTable can be activated only in case of MC processing."); + } } /// \brief core function to fill tables in MC @@ -287,12 +406,12 @@ struct HfTreeCreatorLcToPKPi { /// \param mcCollisions MC collision table /// \param candidates Lc->pKpi candidate table /// \param particles Generated particle table - template + template void fillTablesMc(Colls const& collisions, aod::McCollisions const&, - soa::Join const& candidates, + CandType const& candidates, soa::Join const& particles, - TracksWPid const&, aod::BCs const&) + soa::Join const&, aod::BCs const&) { // Filling event properties @@ -330,29 +449,43 @@ struct HfTreeCreatorLcToPKPi { } // Filling candidate properties - if (fillCandidateLiteTable) { - rowCandidateLite.reserve(candidates.size()); + if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) { + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size() * 2); + } else { + rowCandidateFull.reserve(candidates.size() * 2); + } } else { - rowCandidateFull.reserve(candidates.size()); + rowCandidateKF.reserve(candidates.size() * 2); } if (fillCollIdTable) { /// save also candidate collision indices rowCollisionId.reserve(candidates.size()); } + if (fillCandidateMcTable) { + rowCandidateMC.reserve(candidates.size() * 2); + } for (const auto& candidate : candidates) { - auto trackPos1 = candidate.prong0_as(); // positive daughter (negative for the antiparticles) - auto trackNeg = candidate.prong1_as(); // negative daughter (positive for the antiparticles) - auto trackPos2 = candidate.prong2_as(); // positive daughter (negative for the antiparticles) - bool isMcCandidateSignal = std::abs(candidate.flagMcMatchRec()) == (1 << o2::aod::hf_cand_3prong::DecayType::LcToPKPi); - auto fillTable = [&](int CandFlag, - int FunctionSelection, - float FunctionInvMass, - float FunctionCt, - float FunctionY, - float FunctionE, - float FunctionInvMassKPi) { - double pseudoRndm = trackPos1.pt() * 1000. - (int64_t)(trackPos1.pt() * 1000); - if (FunctionSelection >= 1 && (/*keep all*/ (!keepOnlySignalMc && !keepOnlyBkg) || /*keep only signal*/ (keepOnlySignalMc && isMcCandidateSignal) || /*keep only background and downsample it*/ (keepOnlyBkg && !isMcCandidateSignal && (candidate.pt() > downSampleBkgPtMax || (pseudoRndm < downSampleBkgFactor && candidate.pt() < downSampleBkgPtMax))))) { + auto trackPos1 = candidate.template prong0_as>(); // positive daughter (negative for the antiparticles) + auto trackNeg = candidate.template prong1_as>(); // negative daughter (positive for the antiparticles) + auto trackPos2 = candidate.template prong2_as>(); // positive daughter (negative for the antiparticles) + auto fillTable = [&](int CandFlag) { + double pseudoRndm = trackPos1.pt() * 1000. - static_cast(trackPos1.pt() * 1000); + const int FunctionSelection = CandFlag == 0 ? candidate.isSelLcToPKPi() : candidate.isSelLcToPiKP(); + const int sigbgstatus = determineSignalBgStatus(candidate, CandFlag); + bool isMcCandidateSignal = (sigbgstatus == Prompt) || (sigbgstatus == NonPrompt); + if (FunctionSelection >= selectionFlagLc && (/*keep all*/ (!keepOnlySignalMc && !keepOnlyBkg) || /*keep only signal*/ (keepOnlySignalMc && isMcCandidateSignal) || /*keep only background and downsample it*/ (keepOnlyBkg && !isMcCandidateSignal && (candidate.pt() > downSampleBkgPtMax || (pseudoRndm < downSampleBkgFactor && candidate.pt() < downSampleBkgPtMax))))) { + float FunctionInvMass, FunctionInvMassKPi; + if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) { + FunctionInvMass = CandFlag == 0 ? hfHelper.invMassLcToPKPi(candidate) : hfHelper.invMassLcToPiKP(candidate); + FunctionInvMassKPi = CandFlag == 0 ? hfHelper.invMassKPiPairLcToPKPi(candidate) : hfHelper.invMassKPiPairLcToPiKP(candidate); + } else { + FunctionInvMass = CandFlag == 0 ? candidate.kfMassPKPi() : candidate.kfMassPiKP(); + FunctionInvMassKPi = CandFlag == 0 ? candidate.kfMassKPi() : candidate.kfMassPiK(); + } + const float FunctionCt = hfHelper.ctLc(candidate); + const float FunctionY = hfHelper.yLc(candidate); + const float FunctionE = hfHelper.eLc(candidate); if (fillCandidateLiteTable) { rowCandidateLite( candidate.posX(), @@ -360,25 +493,15 @@ struct HfTreeCreatorLcToPKPi { candidate.posZ(), candidate.nProngsContributorsPV(), candidate.bitmapProngsContributorsPV(), - // candidate.errorDecayLength(), - // candidate.errorDecayLengthXY(), candidate.chi2PCA(), candidate.decayLength(), candidate.decayLengthXY(), - // candidate.decayLengthNormalised(), - // candidate.decayLengthXYNormalised(), - // candidate.impactParameterNormalised0(), candidate.ptProng0(), - // candidate.impactParameterNormalised1(), candidate.ptProng1(), - // candidate.impactParameterNormalised2(), candidate.ptProng2(), candidate.impactParameter0(), candidate.impactParameter1(), candidate.impactParameter2(), - // candidate.errorImpactParameter0(), - // candidate.errorImpactParameter1(), - // candidate.errorImpactParameter2(), trackPos1.tpcNSigmaPi(), trackPos1.tpcNSigmaPr(), trackPos1.tofNSigmaPi(), @@ -408,13 +531,11 @@ struct HfTreeCreatorLcToPKPi { candidate.isCandidateSwapped(), candidate.flagMcDecayChanRec(), FunctionInvMassKPi); - // candidate.globalIndex()); if (fillCollIdTable) { /// save also candidate collision indices rowCollisionId(candidate.collisionId()); } - } else { rowCandidateFull( candidate.collisionId(), @@ -491,11 +612,92 @@ struct HfTreeCreatorLcToPKPi { candidate.flagMcDecayChanRec(), FunctionInvMassKPi); } + + if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) { + const float svX = candidate.xSecondaryVertex(); + const float svY = candidate.ySecondaryVertex(); + const float svZ = candidate.zSecondaryVertex(); + const float svErrX = candidate.kfXError(); + const float svErrY = candidate.kfYError(); + const float svErrZ = candidate.kfZError(); + const float pvErrX = candidate.kfXPVError(); + const float pvErrY = candidate.kfYPVError(); + const float pvErrZ = candidate.kfZPVError(); + const float chi2primProton = CandFlag == 0 ? candidate.kfChi2PrimProng0() : candidate.kfChi2PrimProng2(); + const float chi2primKaon = candidate.kfChi2PrimProng1(); + const float chi2primPion = CandFlag == 0 ? candidate.kfChi2PrimProng2() : candidate.kfChi2PrimProng0(); + const float dcaProtonKaon = CandFlag == 0 ? candidate.kfDcaProng0Prong1() : candidate.kfDcaProng1Prong2(); + const float dcaProtonPion = candidate.kfDcaProng0Prong2(); + const float dcaPionKaon = CandFlag == 0 ? candidate.kfDcaProng1Prong2() : candidate.kfDcaProng0Prong1(); + const float chi2GeoProtonKaon = CandFlag == 0 ? candidate.kfChi2GeoProng0Prong1() : candidate.kfChi2GeoProng1Prong2(); + const float chi2GeoProtonPion = candidate.kfChi2GeoProng0Prong2(); + const float chi2GeoPionKaon = CandFlag == 0 ? candidate.kfChi2GeoProng1Prong2() : candidate.kfChi2GeoProng0Prong1(); + const float chi2Geo = candidate.kfChi2Geo(); + const float chi2Topo = candidate.kfChi2Topo(); + const float l = candidate.kfDecayLength(); + const float dl = candidate.kfDecayLengthError(); + const float pt = std::sqrt(candidate.kfPx() * candidate.kfPx() + candidate.kfPy() * candidate.kfPy()); + const float deltaPt = std::sqrt(candidate.kfPx() * candidate.kfPx() * candidate.kfErrorPx() * candidate.kfErrorPx() + + candidate.kfPy() * candidate.kfPy() * candidate.kfErrorPy() * candidate.kfErrorPy()) / + pt; + const float p = std::sqrt(pt * pt + candidate.kfPz() * candidate.kfPz()); + const float deltaP = std::sqrt(pt * pt * deltaPt * deltaPt + + candidate.kfPz() * candidate.kfPz() * candidate.kfErrorPz() * candidate.kfErrorPz()) / + p; + const float t = l * MassLambdaCPlus / LightSpeedCm2PS / p; + const float deltaT = dl * MassLambdaCPlus / LightSpeedCm2PS / p; + const float mass = CandFlag == 0 ? candidate.kfMassPKPi() : candidate.kfMassPiKP(); + rowCandidateKF( + svX, svY, svZ, svErrX, svErrY, svErrZ, + pvErrX, pvErrY, pvErrZ, + chi2primProton, chi2primKaon, chi2primPion, + dcaProtonKaon, dcaProtonPion, dcaPionKaon, + chi2GeoProtonKaon, chi2GeoProtonPion, chi2GeoPionKaon, + chi2Geo, chi2Topo, l, dl, l / dl, t, deltaT, + mass, p, pt, deltaP, deltaPt, + FunctionSelection, sigbgstatus); + } + if (fillCandidateMcTable) { + float p, pt, svX, svY, svZ, pvX, pvY, pvZ, l, t; + if (!isMcCandidateSignal) { + p = UndefValueFloat; + pt = UndefValueFloat; + svX = UndefValueFloat; + svY = UndefValueFloat; + svZ = UndefValueFloat; + pvX = UndefValueFloat; + pvY = UndefValueFloat; + pvZ = UndefValueFloat; + l = UndefValueFloat; + t = UndefValueFloat; + } else { + auto mcParticleProng0 = candidate.template prong0_as>().template mcParticle_as>(); + auto indexMother = RecoDecay::getMother(particles, mcParticleProng0, o2::constants::physics::Pdg::kLambdaCPlus, true); + auto particleMother = particles.rawIteratorAt(indexMother); + auto mcCollision = particleMother.template mcCollision_as(); + p = particleMother.p(); + pt = particleMother.pt(); + const float p2m = p / MassLambdaCPlus; + const float gamma = std::sqrt(1 + p2m * p2m); // mother's particle Lorentz factor + pvX = mcCollision.posX(); + pvY = mcCollision.posY(); + pvZ = mcCollision.posZ(); + svX = mcParticleProng0.vx(); + svY = mcParticleProng0.vy(); + svZ = mcParticleProng0.vz(); + l = std::sqrt((svX - pvX) * (svX - pvX) + (svY - pvY) * (svY - pvY) + (svZ - pvZ) * (svZ - pvZ)); + t = mcParticleProng0.vt() * NanoToPico / gamma; // from ns to ps * from lab time to proper time + } + rowCandidateMC( + p, pt, + svX, svY, svZ, l, t, + pvX, pvY, pvZ); + } } }; - fillTable(0, candidate.isSelLcToPKPi(), hfHelper.invMassLcToPKPi(candidate), hfHelper.ctLc(candidate), hfHelper.yLc(candidate), hfHelper.eLc(candidate), hfHelper.invMassKPiPairLcToPKPi(candidate)); - fillTable(1, candidate.isSelLcToPiKP(), hfHelper.invMassLcToPiKP(candidate), hfHelper.ctLc(candidate), hfHelper.yLc(candidate), hfHelper.eLc(candidate), hfHelper.invMassKPiPairLcToPiKP(candidate)); + fillTable(0); + fillTable(1); } // Filling particle properties @@ -520,15 +722,48 @@ struct HfTreeCreatorLcToPKPi { /// \param particles Generated particle table /// \param tracks Track table /// \param bcs Bunch-crossing table - void processMcNoCentrality(soa::Join const& collisions, - aod::McCollisions const& mcCollisions, - soa::Join const& candidates, - soa::Join const& particles, - TracksWPid const& tracks, aod::BCs const& bcs) + void processMcNoCentralityWithDCAFitterN(soa::Join const& collisions, + aod::McCollisions const& mcCollisions, + soa::Join const& candidates, + soa::Join const& particles, + soa::Join const& tracks, aod::BCs const& bcs) + { + fillTablesMc(collisions, mcCollisions, candidates, particles, tracks, bcs); + } + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMcNoCentralityWithDCAFitterN, "Process MC tree writer w/o centrality with DCAFitterN", false); + + /// \brief process function for MC with centrality + /// \param collisions Collision table with join of the centrality table + /// \param mcCollisions MC collision table + /// \param candidates Lc->pKpi candidate table + /// \param tracks Track table + /// \param bcs Bunch-crossing table + void processMcWithCentralityWithDCAFitterN(soa::Join const& collisions, + aod::McCollisions const& mcCollisions, + soa::Join const& candidates, + soa::Join const& particles, + soa::Join const& tracks, aod::BCs const& bcs) + { + fillTablesMc(collisions, mcCollisions, candidates, particles, tracks, bcs); + } + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMcWithCentralityWithDCAFitterN, "Process MC tree writer with centrality with DCAFitterN", false); + + /// \brief process function for MC w/o centrality + /// \param collisions Collision table w/o join of the centrality table + /// \param mcCollisions MC collision table + /// \param candidates Lc->pKpi candidate table + /// \param particles Generated particle table + /// \param tracks Track table + /// \param bcs Bunch-crossing table + void processMcNoCentralityWithKFParticle(soa::Join const& collisions, + aod::McCollisions const& mcCollisions, + soa::Join const& candidates, + soa::Join const& particles, + soa::Join const& tracks, aod::BCs const& bcs) { - fillTablesMc(collisions, mcCollisions, candidates, particles, tracks, bcs); + fillTablesMc(collisions, mcCollisions, candidates, particles, tracks, bcs); } - PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMcNoCentrality, "Process MC tree writer w/o centrality", false); + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMcNoCentralityWithKFParticle, "Process MC tree writer w/o centrality with KFParticle", false); /// \brief process function for MC with centrality /// \param collisions Collision table with join of the centrality table @@ -536,22 +771,22 @@ struct HfTreeCreatorLcToPKPi { /// \param candidates Lc->pKpi candidate table /// \param tracks Track table /// \param bcs Bunch-crossing table - void processMcWithCentrality(soa::Join const& collisions, - aod::McCollisions const& mcCollisions, - soa::Join const& candidates, - soa::Join const& particles, - TracksWPid const& tracks, aod::BCs const& bcs) + void processMcWithCentralityWithKFParticle(soa::Join const& collisions, + aod::McCollisions const& mcCollisions, + soa::Join const& candidates, + soa::Join const& particles, + soa::Join const& tracks, aod::BCs const& bcs) { - fillTablesMc(collisions, mcCollisions, candidates, particles, tracks, bcs); + fillTablesMc(collisions, mcCollisions, candidates, particles, tracks, bcs); } - PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMcWithCentrality, "Process MC tree writer with centrality", false); + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processMcWithCentralityWithKFParticle, "Process MC tree writer with centrality with KFParticle", false); /// \brief core function to fill tables in data /// \param collisions Collision table /// \param candidates Lc->pKpi candidate table - template + template void fillTablesData(Colls const& collisions, - soa::Join const& candidates, + CandType const& candidates, TracksWPid const&, aod::BCs const&) { @@ -590,28 +825,38 @@ struct HfTreeCreatorLcToPKPi { } // Filling candidate properties - if (fillCandidateLiteTable) { - rowCandidateLite.reserve(candidates.size()); + if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) { + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(candidates.size() * 2); + } else { + rowCandidateFull.reserve(candidates.size() * 2); + } } else { - rowCandidateFull.reserve(candidates.size()); + rowCandidateKF.reserve(candidates.size() * 2); } if (fillCollIdTable) { /// save also candidate collision indices rowCollisionId.reserve(candidates.size()); } for (const auto& candidate : candidates) { - auto trackPos1 = candidate.prong0_as(); // positive daughter (negative for the antiparticles) - auto trackNeg = candidate.prong1_as(); // negative daughter (positive for the antiparticles) - auto trackPos2 = candidate.prong2_as(); // positive daughter (negative for the antiparticles) - auto fillTable = [&](int CandFlag, - int FunctionSelection, - float FunctionInvMass, - float FunctionCt, - float FunctionY, - float FunctionE, - float FunctionInvMassKPi) { - double pseudoRndm = trackPos1.pt() * 1000. - (int64_t)(trackPos1.pt() * 1000); - if (FunctionSelection >= 1 && (candidate.pt() > downSampleBkgPtMax || (pseudoRndm < downSampleBkgFactor && candidate.pt() < downSampleBkgPtMax))) { + auto trackPos1 = candidate.template prong0_as(); // positive daughter (negative for the antiparticles) + auto trackNeg = candidate.template prong1_as(); // negative daughter (positive for the antiparticles) + auto trackPos2 = candidate.template prong2_as(); // positive daughter (negative for the antiparticles) + auto fillTable = [&](int CandFlag) { + double pseudoRndm = trackPos1.pt() * 1000. - static_cast(trackPos1.pt() * 1000); + const int FunctionSelection = CandFlag == 0 ? candidate.isSelLcToPKPi() : candidate.isSelLcToPiKP(); + if (FunctionSelection >= selectionFlagLc && (candidate.pt() > downSampleBkgPtMax || (pseudoRndm < downSampleBkgFactor && candidate.pt() < downSampleBkgPtMax))) { + float FunctionInvMass, FunctionInvMassKPi; + if constexpr (reconstructionType == aod::hf_cand::VertexerType::DCAFitter) { + FunctionInvMass = CandFlag == 0 ? hfHelper.invMassLcToPKPi(candidate) : hfHelper.invMassLcToPiKP(candidate); + FunctionInvMassKPi = CandFlag == 0 ? hfHelper.invMassKPiPairLcToPKPi(candidate) : hfHelper.invMassKPiPairLcToPiKP(candidate); + } else { + FunctionInvMass = CandFlag == 0 ? candidate.kfMassPKPi() : candidate.kfMassPiKP(); + FunctionInvMassKPi = CandFlag == 0 ? candidate.kfMassKPi() : candidate.kfMassPiK(); + } + const float FunctionCt = hfHelper.ctLc(candidate); + const float FunctionY = hfHelper.yLc(candidate); + const float FunctionE = hfHelper.eLc(candidate); if (fillCandidateLiteTable) { rowCandidateLite( candidate.posX(), @@ -619,25 +864,15 @@ struct HfTreeCreatorLcToPKPi { candidate.posZ(), candidate.nProngsContributorsPV(), candidate.bitmapProngsContributorsPV(), - // candidate.errorDecayLength(), - // candidate.errorDecayLengthXY(), candidate.chi2PCA(), candidate.decayLength(), candidate.decayLengthXY(), - // candidate.decayLengthNormalised(), - // candidate.decayLengthXYNormalised(), - // candidate.impactParameterNormalised0(), candidate.ptProng0(), - // candidate.impactParameterNormalised1(), candidate.ptProng1(), - // candidate.impactParameterNormalised2(), candidate.ptProng2(), candidate.impactParameter0(), candidate.impactParameter1(), candidate.impactParameter2(), - // candidate.errorImpactParameter0(), - // candidate.errorImpactParameter1(), - // candidate.errorImpactParameter2(), trackPos1.tpcNSigmaPi(), trackPos1.tpcNSigmaPr(), trackPos1.tofNSigmaPi(), @@ -667,7 +902,6 @@ struct HfTreeCreatorLcToPKPi { 0., -1, FunctionInvMassKPi); - // candidate.globalIndex()); if (fillCollIdTable) { /// save also candidate collision indices @@ -750,11 +984,56 @@ struct HfTreeCreatorLcToPKPi { -1, FunctionInvMassKPi); } + + if constexpr (reconstructionType == aod::hf_cand::VertexerType::KfParticle) { + const float X = candidate.xSecondaryVertex(); + const float Y = candidate.ySecondaryVertex(); + const float Z = candidate.zSecondaryVertex(); + const float ErrX = candidate.kfXError(); + const float ErrY = candidate.kfYError(); + const float ErrZ = candidate.kfZError(); + const float ErrPVX = candidate.kfXPVError(); + const float ErrPVY = candidate.kfYPVError(); + const float ErrPVZ = candidate.kfZPVError(); + const float chi2prim_proton = CandFlag == 0 ? candidate.kfChi2PrimProng0() : candidate.kfChi2PrimProng2(); + const float chi2prim_kaon = candidate.kfChi2PrimProng1(); + const float chi2prim_pion = CandFlag == 0 ? candidate.kfChi2PrimProng2() : candidate.kfChi2PrimProng0(); + const float dca_proton_kaon = CandFlag == 0 ? candidate.kfDcaProng0Prong1() : candidate.kfDcaProng1Prong2(); + const float dca_proton_pion = candidate.kfDcaProng0Prong2(); + const float dca_pion_kaon = CandFlag == 0 ? candidate.kfDcaProng1Prong2() : candidate.kfDcaProng0Prong1(); + const float chi2Geo_proton_kaon = CandFlag == 0 ? candidate.kfChi2GeoProng0Prong1() : candidate.kfChi2GeoProng1Prong2(); + const float chi2Geo_proton_pion = candidate.kfChi2GeoProng0Prong2(); + const float chi2Geo_pion_kaon = CandFlag == 0 ? candidate.kfChi2GeoProng1Prong2() : candidate.kfChi2GeoProng0Prong1(); + const float chi2Geo = candidate.kfChi2Geo(); + const float chi2Topo = candidate.kfChi2Topo(); + const float l = candidate.kfDecayLength(); + const float dl = candidate.kfDecayLengthError(); + const float pt = std::sqrt(candidate.kfPx() * candidate.kfPx() + candidate.kfPy() * candidate.kfPy()); + const float deltaPt = std::sqrt(candidate.kfPx() * candidate.kfPx() * candidate.kfErrorPx() * candidate.kfErrorPx() + + candidate.kfPy() * candidate.kfPy() * candidate.kfErrorPy() * candidate.kfErrorPy()) / + pt; + const float p = std::sqrt(pt * pt + candidate.kfPz() * candidate.kfPz()); + const float deltaP = std::sqrt(pt * pt * deltaPt * deltaPt + + candidate.kfPz() * candidate.kfPz() * candidate.kfErrorPz() * candidate.kfErrorPz()) / + p; + const float T = l * MassLambdaCPlus / LightSpeedCm2PS / p; + const float deltaT = dl * MassLambdaCPlus / LightSpeedCm2PS / p; + const float mass = CandFlag == 0 ? candidate.kfMassPKPi() : candidate.kfMassPiKP(); + rowCandidateKF( + X, Y, Z, ErrX, ErrY, ErrZ, + ErrPVX, ErrPVY, ErrPVZ, + chi2prim_proton, chi2prim_kaon, chi2prim_pion, + dca_proton_kaon, dca_proton_pion, dca_pion_kaon, + chi2Geo_proton_kaon, chi2Geo_proton_pion, chi2Geo_pion_kaon, + chi2Geo, chi2Topo, l, dl, l / dl, T, deltaT, + mass, p, pt, deltaP, deltaPt, + FunctionSelection, UndefValueInt); + } } }; - fillTable(0, candidate.isSelLcToPKPi(), hfHelper.invMassLcToPKPi(candidate), hfHelper.ctLc(candidate), hfHelper.yLc(candidate), hfHelper.eLc(candidate), hfHelper.invMassKPiPairLcToPKPi(candidate)); - fillTable(1, candidate.isSelLcToPiKP(), hfHelper.invMassLcToPiKP(candidate), hfHelper.ctLc(candidate), hfHelper.yLc(candidate), hfHelper.eLc(candidate), hfHelper.invMassKPiPairLcToPiKP(candidate)); + fillTable(0); + fillTable(1); } } @@ -763,26 +1042,52 @@ struct HfTreeCreatorLcToPKPi { /// \param candidates Lc->pKpi candidate table /// \param tracks Track table /// \param bcs Bunch-crossing table - void processDataNoCentrality(soa::Join const& collisions, - soa::Join const& candidates, - TracksWPid const& tracks, aod::BCs const& bcs) + void processDataNoCentralityWithDCAFitterN(soa::Join const& collisions, + soa::Join const& candidates, + TracksWPid const& tracks, aod::BCs const& bcs) + { + fillTablesData(collisions, candidates, tracks, bcs); + } + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processDataNoCentralityWithDCAFitterN, "Process data tree writer w/o centrality with DCAFitterN", false); + + /// \brief process function for data with centrality + /// \param collisions Collision table with join of the centrality table + /// \param candidates Lc->pKpi candidate table + /// \param tracks Track table + /// \param bcs Bunch-crossing table + void processDataWithCentralityWithDCAFitterN(soa::Join const& collisions, + soa::Join const& candidates, + TracksWPid const& tracks, aod::BCs const& bcs) + { + fillTablesData(collisions, candidates, tracks, bcs); + } + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processDataWithCentralityWithDCAFitterN, "Process data tree writer with centrality with DCAFitterN", true); + + /// \brief process function for data w/o centrality + /// \param collisions Collision table w/o join of the centrality table + /// \param candidates Lc->pKpi candidate table + /// \param tracks Track table + /// \param bcs Bunch-crossing table + void processDataNoCentralityWithKFParticle(soa::Join const& collisions, + soa::Join const& candidates, + TracksWPid const& tracks, aod::BCs const& bcs) { - fillTablesData(collisions, candidates, tracks, bcs); + fillTablesData(collisions, candidates, tracks, bcs); } - PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processDataNoCentrality, "Process data tree writer w/o centrality", false); + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processDataNoCentralityWithKFParticle, "Process data tree writer w/o centrality with KFParticle", false); /// \brief process function for data with centrality /// \param collisions Collision table with join of the centrality table /// \param candidates Lc->pKpi candidate table /// \param tracks Track table /// \param bcs Bunch-crossing table - void processDataWithCentrality(soa::Join const& collisions, - soa::Join const& candidates, - TracksWPid const& tracks, aod::BCs const& bcs) + void processDataWithCentralityWithKFParticle(soa::Join const& collisions, + soa::Join const& candidates, + TracksWPid const& tracks, aod::BCs const& bcs) { - fillTablesData(collisions, candidates, tracks, bcs); + fillTablesData(collisions, candidates, tracks, bcs); } - PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processDataWithCentrality, "Process data tree writer with centrality", true); + PROCESS_SWITCH(HfTreeCreatorLcToPKPi, processDataWithCentralityWithKFParticle, "Process data tree writer with centrality with KFParticle", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/Tools/KFparticle/KFUtilities.h b/Tools/KFparticle/KFUtilities.h index 7f03cb1c4af..e27b2340866 100644 --- a/Tools/KFparticle/KFUtilities.h +++ b/Tools/KFparticle/KFUtilities.h @@ -22,6 +22,8 @@ #define HomogeneousField #endif +#include + #include // FIXME #include "KFParticle.h" @@ -313,4 +315,91 @@ float ldlXYFromKF(KFParticle kfpParticle, KFParticle PV) return l_particle / dl_particle; } +/// @brief squared distance between track and primary vertex normalised by its uncertainty evaluated in matrix form +/// @param track KFParticle track (must be passed as a copy) +/// @param vtx KFParticle primary vertex +/// @return chi2 to primary vertex +float kfCalculateChi2ToPrimaryVertex(KFParticle track, const KFParticle& vtx) +{ + const float PvPoint[3] = {vtx.X(), vtx.Y(), vtx.Z()}; + + track.TransportToPoint(PvPoint); + return track.GetDeviationFromVertex(vtx); +} + +/// @brief prong's momentum in the secondary (decay) vertex +/// @param track KFParticle track (must be passed as a copy) +/// @param vtx KFParticle secondary vertex +/// @return array with components of prong's momentum in the secondary (decay) vertex +std::array kfCalculateProngMomentumInSecondaryVertex(KFParticle track, const KFParticle& vtx) +{ + const float SvPoint[3] = {vtx.X(), vtx.Y(), vtx.Z()}; + + track.TransportToPoint(SvPoint); + return {track.GetPx(), track.GetPy(), track.GetPz()}; +} + +/// @brief distance of closest approach between two tracks, cm +/// @param track1 KFParticle first track (must be passed as a copy) +/// @param track2 KFParticle second track (must be passed as a copy) +/// @return DCA [cm] in the PCA +float kfCalculateDistanceBetweenParticles(KFParticle track1, KFParticle track2) +{ + float dS[2]; + float dsdr[4][6]; + float params1[8], params2[8]; + float covs1[36], covs2[36]; + track1.GetDStoParticle(track2, dS, dsdr); + track1.Transport(dS[0], dsdr[0], params1, covs1); + track2.Transport(dS[1], dsdr[3], params2, covs2); + const float dx = params1[0] - params2[0]; + const float dy = params1[1] - params2[1]; + const float dz = params1[2] - params2[2]; + return std::sqrt(dx * dx + dy * dy + dz * dz); +} + +/// @brief squared distance between two tracks normalised by its uncertainty evaluated in matrix form +/// @param track1 KFParticle first track (must be passed as a copy) +/// @param track2 KFParticle second track (must be passed as a copy) +/// @return chi2 in PCA +float kfCalculateChi2geoBetweenParticles(KFParticle track1, KFParticle track2) +{ + KFParticle kfPair; + const KFParticle* kfDaughters[3] = {&track1, &track2}; + kfPair.SetConstructMethod(2); + kfPair.Construct(kfDaughters, 2); + + return kfPair.Chi2() / kfPair.NDF(); +} + +/// @brief signed distance between primary and secondary vertex and its uncertainty, cm +/// @param candidate KFParticle decay candidate (must be passed as a copy) +/// @param vtx KFParticle primary vertex +/// @return pair of l and delta l +std::pair kfCalculateLdL(KFParticle candidate, const KFParticle& vtx) +{ + float l, dl; + candidate.SetProductionVertex(vtx); + candidate.KFParticleBase::GetDecayLength(l, dl); + + return std::make_pair(l, dl); +} + +/// @brief Z projection of the impact parameter from the track to the primary vertex, cm +/// @param candidate KFParticle prong +/// @param vtx KFParticle primary vertex +/// @return pair of impact parameter and its error +std::pair kfCalculateImpactParameterZ(const KFParticle& candidate, const KFParticle& vtx) +{ + float distanceToVertexXY, errDistanceToVertexXY; + candidate.GetDistanceFromVertexXY(vtx, distanceToVertexXY, errDistanceToVertexXY); + const float distanceToVertex = candidate.GetDistanceFromVertex(vtx); + const float chi2ToVertex = candidate.GetDeviationFromVertex(vtx); + const float distanceToVertexZ2 = distanceToVertex * distanceToVertex - distanceToVertexXY * distanceToVertexXY; + const float distanceToVertexZ = distanceToVertexZ2 > 0 ? std::sqrt(distanceToVertexZ2) : -std::sqrt(-distanceToVertexZ2); + const float errDistanceToVertexZ2 = (distanceToVertex * distanceToVertex * distanceToVertex * distanceToVertex / chi2ToVertex - distanceToVertexXY * distanceToVertexXY * errDistanceToVertexXY * errDistanceToVertexXY) / distanceToVertexZ2; + const float errDistanceToVertexZ = errDistanceToVertexZ2 > 0 ? std::sqrt(errDistanceToVertexZ2) : -std::sqrt(-errDistanceToVertexZ2); + return std::make_pair(distanceToVertexZ, errDistanceToVertexZ); +} + #endif // TOOLS_KFPARTICLE_KFUTILITIES_H_