From 5a30cebc20c87798479322c00340dbd9bb39586b Mon Sep 17 00:00:00 2001 From: Maximiliano Puccio Date: Sat, 8 Apr 2023 03:42:55 +0200 Subject: [PATCH 01/24] Skip DF with 0 collisions (#2351) --- EventFiltering/selectBCRange.cxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/EventFiltering/selectBCRange.cxx b/EventFiltering/selectBCRange.cxx index 87fadf094c1..ff288657716 100644 --- a/EventFiltering/selectBCRange.cxx +++ b/EventFiltering/selectBCRange.cxx @@ -74,6 +74,11 @@ struct BCRangeSelector { if (cols.size() != fdecs.size()) { throw std::runtime_error("Collision table and CefpDecision do not have the same number of rows! "); } + if (cols.size() == 0) { + LOGF(warning, "No collisions found!"); + pc.outputs().snapshot({"PPF", "IFRAMES", 0, Lifetime::Timeframe}, res); + return; + } // 1. loop over collisions auto filt = fdecs.begin(); @@ -149,7 +154,7 @@ struct BCRangeSelector { return a.first < b.first; }); if (bcRanges.empty()) { - LOGF(error, "No BCs selected! This should not happen! Adding a bogus BC range to avoid crashes."); + LOGF(warning, "No BCs selected! This should not happen! Adding a bogus BC range to avoid crashes."); bcRanges.push_back(std::make_pair(0, 0)); } std::vector> bcRangesMerged(1, bcRanges[0]); From f8140e6632ad35a49bf8166136b834d77e103c25 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 12 Jul 2023 11:49:27 +0200 Subject: [PATCH 02/24] PWGHF: Add KFParticle to D0 creator --- PWGHF/TableProducer/CMakeLists.txt | 2 +- .../TableProducer/candidateCreator2Prong.cxx | 307 +++++++++++++----- 2 files changed, 220 insertions(+), 89 deletions(-) diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index 749c6381277..9140dcbecbc 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -25,7 +25,7 @@ o2physics_add_dpl_workflow(refit-pv-dummy o2physics_add_dpl_workflow(candidate-creator-2prong SOURCES candidateCreator2Prong.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(candidate-creator-3prong diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 9d32f800eaa..1ded149473a 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -25,6 +25,14 @@ #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/Utils/utilsBfieldCCDB.h" +/// includes KFParticle +#include "Tools/KFparticle/KFUtilities.h" +#include "KFParticle.h" +#include "KFPTrack.h" +#include "KFPVertex.h" +#include "KFParticleBase.h" +#include "KFVertex.h" + using namespace o2; using namespace o2::framework; using namespace o2::aod::hf_cand; @@ -35,6 +43,8 @@ struct HfCandidateCreator2Prong { Produces rowCandidateBase; // vertexing + Configurable useDCAFitterN{"useDCAFitterN", true, "calculate vertex with DCAFitterN"}; + Configurable useKFParticle{"useKFParticle", false, "calculate vertex with KFParticle"}; Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; @@ -74,6 +84,8 @@ struct HfCandidateCreator2Prong { OutputObj hCovSVZZ{TH1F("hCovSVZZ", "2-prong candidates;ZZ element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hDcaXYProngs{TH2F("hDcaXYProngs", "DCAxy of 2-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{xy}) (#mum);entries", 100, 0., 20., 200, -500., 500.)}; OutputObj hDcaZProngs{TH2F("hDcaZProngs", "DCAz of 2-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{z}) (#mum);entries", 100, 0., 20., 200, -500., 500.)}; + OutputObj hUseKForDCAFitter{TH1F("hUseKForDCAFitter", ";Use KF or DCAFitterN;entries", 2, -0.5, 1.5)}; // 0 for KF, 1 for DCAFitterN + OutputObj hKFMassFailed{TH1F("hKFMassFailed", ";KF geometrical fitting failed;entries", 3, -0.5, 2.5)}; // 0 for D0 failed, 1 for D0bar failed void init(InitContext const&) { @@ -95,14 +107,26 @@ struct HfCandidateCreator2Prong { { // 2-prong vertex fitter o2::vertexing::DCAFitterN<2> df; - // df.setBz(bz); - df.setPropagateToPCA(propagateToPCA); - df.setMaxR(maxR); - df.setMaxDZIni(maxDZIni); - df.setMinParamChange(minParamChange); - df.setMinRelChi2Change(minRelChi2Change); - df.setUseAbsDCA(useAbsDCA); - df.setWeightedFinalPCA(useWeightedFinalPCA); + if (useDCAFitterN) { + // df.setBz(bz); + df.setPropagateToPCA(propagateToPCA); + df.setMaxR(maxR); + df.setMaxDZIni(maxDZIni); + df.setMinParamChange(minParamChange); + df.setMinRelChi2Change(minRelChi2Change); + df.setUseAbsDCA(useAbsDCA); + df.setWeightedFinalPCA(useWeightedFinalPCA); + hUseKForDCAFitter->Fill(1.0); + } + if (useKFParticle) { + hUseKForDCAFitter->Fill(0.0); + } + if (!useDCAFitterN && !useKFParticle) { + LOGP(fatal, "At least one package between useDCAFitterN and useDCAFitterN should be enabled at a time."); + } + if (useDCAFitterN && useKFParticle) { + LOGP(fatal, "Only one package between useDCAFitterN and useDCAFitterN can be enabled at a time."); + } // loop over pairs of track indices for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { @@ -124,88 +148,195 @@ struct HfCandidateCreator2Prong { // df.setBz(bz); /// put it outside the 'if'! Otherwise we have a difference wrt bz Configurable (< 1 permille) in Run2 conv. data // df.print(); } - df.setBz(bz); - // reconstruct the 2-prong secondary vertex - if (df.process(trackParVarPos1, trackParVarNeg1) == 0) { - continue; - } - const auto& secondaryVertex = df.getPCACandidate(); - auto chi2PCA = df.getChi2AtPCACandidate(); - auto covMatrixPCA = df.calcPCACovMatrixFlat(); - hCovSVXX->Fill(covMatrixPCA[0]); // FIXME: Calculation of errorDecayLength(XY) gives wrong values without this line. - hCovSVYY->Fill(covMatrixPCA[2]); - hCovSVXZ->Fill(covMatrixPCA[3]); - hCovSVZZ->Fill(covMatrixPCA[5]); - auto trackParVar0 = df.getTrack(0); - auto trackParVar1 = df.getTrack(1); - - // get track momenta - array pvec0; - array pvec1; - trackParVar0.getPxPyPzGlo(pvec0); - trackParVar1.getPxPyPzGlo(pvec1); - - // get track impact parameters - // This modifies track momenta! - auto primaryVertex = getPrimaryVertex(collision); - auto covMatrixPV = primaryVertex.getCov(); - if constexpr (doPvRefit) { - /// use PV refit - /// Using it in the rowCandidateBase all dynamic columns shall take it into account - // coordinates - primaryVertex.setX(rowTrackIndexProng2.pvRefitX()); - primaryVertex.setY(rowTrackIndexProng2.pvRefitY()); - primaryVertex.setZ(rowTrackIndexProng2.pvRefitZ()); - // covariance matrix - primaryVertex.setSigmaX2(rowTrackIndexProng2.pvRefitSigmaX2()); - primaryVertex.setSigmaXY(rowTrackIndexProng2.pvRefitSigmaXY()); - primaryVertex.setSigmaY2(rowTrackIndexProng2.pvRefitSigmaY2()); - primaryVertex.setSigmaXZ(rowTrackIndexProng2.pvRefitSigmaXZ()); - primaryVertex.setSigmaYZ(rowTrackIndexProng2.pvRefitSigmaYZ()); - primaryVertex.setSigmaZ2(rowTrackIndexProng2.pvRefitSigmaZ2()); - covMatrixPV = primaryVertex.getCov(); + if (useDCAFitterN) { + df.setBz(bz); + + // reconstruct the 2-prong secondary vertex + if (df.process(trackParVarPos1, trackParVarNeg1) == 0) { + continue; + } + const auto& secondaryVertex = df.getPCACandidate(); + auto chi2PCA = df.getChi2AtPCACandidate(); + auto covMatrixPCA = df.calcPCACovMatrixFlat(); + hCovSVXX->Fill(covMatrixPCA[0]); // FIXME: Calculation of errorDecayLength(XY) gives wrong values without this line. + hCovSVYY->Fill(covMatrixPCA[2]); + hCovSVXZ->Fill(covMatrixPCA[3]); + hCovSVZZ->Fill(covMatrixPCA[5]); + auto trackParVar0 = df.getTrack(0); + auto trackParVar1 = df.getTrack(1); + + // get track momenta + array pvec0; + array pvec1; + trackParVar0.getPxPyPzGlo(pvec0); + trackParVar1.getPxPyPzGlo(pvec1); + + // get track impact parameters + // This modifies track momenta! + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + if constexpr (doPvRefit) { + /// use PV refit + /// Using it in the rowCandidateBase all dynamic columns shall take it into account + // coordinates + primaryVertex.setX(rowTrackIndexProng2.pvRefitX()); + primaryVertex.setY(rowTrackIndexProng2.pvRefitY()); + primaryVertex.setZ(rowTrackIndexProng2.pvRefitZ()); + // covariance matrix + primaryVertex.setSigmaX2(rowTrackIndexProng2.pvRefitSigmaX2()); + primaryVertex.setSigmaXY(rowTrackIndexProng2.pvRefitSigmaXY()); + primaryVertex.setSigmaY2(rowTrackIndexProng2.pvRefitSigmaY2()); + primaryVertex.setSigmaXZ(rowTrackIndexProng2.pvRefitSigmaXZ()); + primaryVertex.setSigmaYZ(rowTrackIndexProng2.pvRefitSigmaYZ()); + primaryVertex.setSigmaZ2(rowTrackIndexProng2.pvRefitSigmaZ2()); + covMatrixPV = primaryVertex.getCov(); + } + hCovPVXX->Fill(covMatrixPV[0]); + hCovPVYY->Fill(covMatrixPV[2]); + hCovPVXZ->Fill(covMatrixPV[3]); + hCovPVZZ->Fill(covMatrixPV[5]); + o2::dataformats::DCA impactParameter0; + o2::dataformats::DCA impactParameter1; + trackParVar0.propagateToDCA(primaryVertex, bz, &impactParameter0); + trackParVar1.propagateToDCA(primaryVertex, bz, &impactParameter1); + hDcaXYProngs->Fill(track0.pt(), impactParameter0.getY() * toMicrometers); + hDcaXYProngs->Fill(track1.pt(), impactParameter1.getY() * toMicrometers); + hDcaZProngs->Fill(track0.pt(), impactParameter0.getZ() * toMicrometers); + hDcaZProngs->Fill(track1.pt(), impactParameter1.getZ() * toMicrometers); + + // get uncertainty of the decay length + double phi, theta; + getPointDirection(array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + // fill candidate table rows + rowCandidateBase(collision.globalIndex(), + primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pvec0[0], pvec0[1], pvec0[2], + pvec1[0], pvec1[1], pvec1[2], + impactParameter0.getY(), impactParameter1.getY(), + std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), + rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), + rowTrackIndexProng2.hfflag()); + + // fill histograms + if (fillHistograms) { + // calculate invariant masses + auto arrayMomenta = array{pvec0, pvec1}; + massPiK = RecoDecay::m(arrayMomenta, array{massPi, massK}); + massKPi = RecoDecay::m(arrayMomenta, array{massK, massPi}); + hMass2->Fill(massPiK); + hMass2->Fill(massKPi); + } } - hCovPVXX->Fill(covMatrixPV[0]); - hCovPVYY->Fill(covMatrixPV[2]); - hCovPVXZ->Fill(covMatrixPV[3]); - hCovPVZZ->Fill(covMatrixPV[5]); - o2::dataformats::DCA impactParameter0; - o2::dataformats::DCA impactParameter1; - trackParVar0.propagateToDCA(primaryVertex, bz, &impactParameter0); - trackParVar1.propagateToDCA(primaryVertex, bz, &impactParameter1); - hDcaXYProngs->Fill(track0.pt(), impactParameter0.getY() * toMicrometers); - hDcaXYProngs->Fill(track1.pt(), impactParameter1.getY() * toMicrometers); - hDcaZProngs->Fill(track0.pt(), impactParameter0.getZ() * toMicrometers); - hDcaZProngs->Fill(track1.pt(), impactParameter1.getZ() * toMicrometers); - - // get uncertainty of the decay length - double phi, theta; - getPointDirection(array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); - auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); - auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - - // fill candidate table rows - rowCandidateBase(collision.globalIndex(), - primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), - secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], - errorDecayLength, errorDecayLengthXY, - chi2PCA, - pvec0[0], pvec0[1], pvec0[2], - pvec1[0], pvec1[1], pvec1[2], - impactParameter0.getY(), impactParameter1.getY(), - std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), - rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), - rowTrackIndexProng2.hfflag()); - - // fill histograms - if (fillHistograms) { - // calculate invariant masses - auto arrayMomenta = array{pvec0, pvec1}; - massPiK = RecoDecay::m(arrayMomenta, array{massPi, massK}); - massKPi = RecoDecay::m(arrayMomenta, array{massK, massPi}); - hMass2->Fill(massPiK); - hMass2->Fill(massKPi); + + if (useKFParticle) { + Float_t 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(rowTrackIndexProng2.pvRefitX(), rowTrackIndexProng2.pvRefitY(), rowTrackIndexProng2.pvRefitZ()); + // covariance matrix + kfpVertex.SetCovarianceMatrix(rowTrackIndexProng2.pvRefitSigmaX2(), rowTrackIndexProng2.pvRefitSigmaXY(), rowTrackIndexProng2.pvRefitSigmaY2(), rowTrackIndexProng2.pvRefitSigmaXZ(), rowTrackIndexProng2.pvRefitSigmaYZ(), rowTrackIndexProng2.pvRefitSigmaZ2()); + } + kfpVertex.GetCovarianceMatrix(covMatrixPV); + KFParticle KFPV(kfpVertex); + hCovPVXX->Fill(covMatrixPV[0]); + hCovPVYY->Fill(covMatrixPV[2]); + hCovPVXZ->Fill(covMatrixPV[3]); + hCovPVZZ->Fill(covMatrixPV[5]); + + KFPTrack kfpTrackPosPi; + KFPTrack kfpTrackNegPi; + KFPTrack kfpTrackPosKa; + KFPTrack kfpTrackNegKa; + + if (track0.sign() == 1 && track1.sign() == -1) { + kfpTrackPosPi = createKFPTrackFromTrack(track0); + kfpTrackNegKa = createKFPTrackFromTrack(track1); + kfpTrackPosKa = createKFPTrackFromTrack(track0); + kfpTrackNegPi = createKFPTrackFromTrack(track1); + } + + KFParticle KFPosPion(kfpTrackPosPi, 211); + KFParticle KFNegPion(kfpTrackNegPi, 211); + KFParticle KFPosKaon(kfpTrackPosKa, 321); + KFParticle KFNegKaon(kfpTrackNegKa, 321); + + float impactParameter0XY = 0., errImpactParameter0XY = 0., impactParameter1XY = 0., errImpactParameter1XY = 0.; + if (!KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + hDcaXYProngs->Fill(track0.pt(), impactParameter0XY * toMicrometers); + hDcaZProngs->Fill(track0.pt(), sqrt(KFPosPion.GetDistanceFromVertex(KFPV) * KFPosPion.GetDistanceFromVertex(KFPV) - impactParameter0XY * impactParameter0XY) * toMicrometers); + } + if (KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + hDcaXYProngs->Fill(track0.pt(), -999.); + hDcaZProngs->Fill(track0.pt(), -999.); + } + if (!KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + hDcaXYProngs->Fill(track1.pt(), impactParameter1XY * toMicrometers); + hDcaZProngs->Fill(track1.pt(), sqrt(KFNegPion.GetDistanceFromVertex(KFPV) * KFNegPion.GetDistanceFromVertex(KFPV) - impactParameter1XY * impactParameter1XY) * toMicrometers); + } + if (KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + hDcaXYProngs->Fill(track1.pt(), -999.); + hDcaZProngs->Fill(track1.pt(), -999.); + } + + KFParticle KFDZero; + const KFParticle* D0Daughters[2] = {&KFPosPion, &KFNegKaon}; + KFDZero.SetConstructMethod(2); + KFDZero.Construct(D0Daughters, 2); + KFParticle KFDZeroBar; + const KFParticle* D0BarDaughters[2] = {&KFNegPion, &KFPosKaon}; + KFDZeroBar.SetConstructMethod(2); + KFDZeroBar.Construct(D0BarDaughters, 2); + + hCovSVXX->Fill(KFDZero.Covariance(0, 0)); + hCovSVYY->Fill(KFDZero.Covariance(1, 1)); + hCovSVXZ->Fill(KFDZero.Covariance(2, 0)); + hCovSVZZ->Fill(KFDZero.Covariance(2, 2)); + auto covMatrixSV = KFDZero.CovarianceMatrix(); + + double phi, theta; + getPointDirection(array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, array{KFDZero.GetX(), KFDZero.GetY(), KFDZero.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.)); + + // fill candidate table rows + rowCandidateBase(collision.globalIndex(), + KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), + KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ(), + errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one + KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one + KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), + KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), + impactParameter0XY, impactParameter1XY, + errImpactParameter0XY, errImpactParameter1XY, + rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), + rowTrackIndexProng2.hfflag()); + + // fill histograms + if (fillHistograms) { + Float_t massDZero = 0., errMassDZero = 0.; + Float_t massDZeroBar = 0., errMassDZeroBar = 0.; + if (!KFDZero.GetMass(massDZero, errMassDZero)) + hMass2->Fill(massDZero); + if (KFDZero.GetMass(massDZero, errMassDZero)) + hKFMassFailed->Fill(0); + if (!KFDZeroBar.GetMass(massDZeroBar, errMassDZeroBar)) + hMass2->Fill(massDZeroBar); + if (KFDZeroBar.GetMass(massDZeroBar, errMassDZeroBar)) + hKFMassFailed->Fill(1); + } } } } From 5bf3cb1cedb336a4cb58c594433a363f611e8620 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Sun, 23 Jul 2023 22:57:24 +0200 Subject: [PATCH 03/24] PWGHF: Added topological chi2 to HfCand2Prong table --- .../DataModel/CandidateReconstructionTables.h | 4 ++ .../TableProducer/candidateCreator2Prong.cxx | 44 ++++++++++++++----- PWGHF/TableProducer/candidateSelectorD0.cxx | 3 ++ 3 files changed, 41 insertions(+), 10 deletions(-) diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index b08dc4e1298..63771869163 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -253,6 +253,8 @@ DECLARE_SOA_COLUMN(ZSecondaryVertex, zSecondaryVertex, float); //! DECLARE_SOA_DYNAMIC_COLUMN(RSecondaryVertex, rSecondaryVertex, //! [](float xVtxS, float yVtxS) -> float { return RecoDecay::sqrtSumOfSquares(xVtxS, yVtxS); }); DECLARE_SOA_COLUMN(Chi2PCA, chi2PCA, float); //! sum of (non-weighted) distances of the secondary vertex to its prongs +DECLARE_SOA_COLUMN(KFTopChi2OverNDF_DZero, kfTopChi2OverNDF_DZero, float); //! +DECLARE_SOA_COLUMN(KFTopChi2OverNDF_DZeroBar, kfTopChi2OverNDF_DZeroBar, float); //! // prong properties DECLARE_SOA_COLUMN(PxProng0, pxProng0, float); //! DECLARE_SOA_COLUMN(PyProng0, pyProng0, float); //! @@ -475,6 +477,8 @@ DECLARE_SOA_TABLE(HfCand2ProngBase, "AOD", "HFCAND2PBASE", //! o2::soa::Index<>, // general columns HFCAND_COLUMNS, + hf_cand::KFTopChi2OverNDF_DZero, + hf_cand::KFTopChi2OverNDF_DZeroBar, // 2-prong specific columns hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 1ded149473a..5ffc1e8f0cd 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -45,6 +45,7 @@ struct HfCandidateCreator2Prong { // vertexing Configurable useDCAFitterN{"useDCAFitterN", true, "calculate vertex with DCAFitterN"}; Configurable useKFParticle{"useKFParticle", false, "calculate vertex with KFParticle"}; + Configurable useKFTop2PV{"useKFTop2PV", false, "constraint KFParticle to PV, only if useKFParticle is true"}; Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; @@ -74,6 +75,7 @@ struct HfCandidateCreator2Prong { double bz = 0.; OutputObj hMass2{TH1F("hMass2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; + OutputObj hKFMassFailed{TH2F("hKFMassFailed", ";KF fit failed;entries", 500, 0., 5., 3, -0.5, 2.5)}; // 0 for D0 failed, 1 for D0bar failed OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hCovPVYY{TH1F("hCovPVYY", "2-prong candidates;YY element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; @@ -217,6 +219,7 @@ struct HfCandidateCreator2Prong { secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], errorDecayLength, errorDecayLengthXY, chi2PCA, + -999., -999., pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2], impactParameter0.getY(), impactParameter1.getY(), @@ -311,12 +314,33 @@ struct HfCandidateCreator2Prong { auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixSV, phi, theta)); auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixSV, phi, 0.)); + Float_t topChi2PerNDF_DZero = -999.; + Float_t topChi2PerNDF_DZeroBar = -999.; + + KFParticle KFDZeroTop2PV, KFDZeroBarTop2PV; + KFParticle KFDZeroTop2PV_DV, KFDZeroBarTop2PV_DV; + if (useKFTop2PV) { + KFDZeroTop2PV = KFDZero; + KFDZeroTop2PV.SetProductionVertex(KFPV); + KFDZeroTop2PV_DV = KFDZeroTop2PV; + KFDZeroTop2PV_DV.TransportToDecayVertex(); + + KFDZeroBarTop2PV = KFDZeroBar; + KFDZeroBarTop2PV.SetProductionVertex(KFPV); + KFDZeroBarTop2PV_DV = KFDZeroBarTop2PV; + KFDZeroBarTop2PV_DV.TransportToDecayVertex(); + + topChi2PerNDF_DZero = KFDZeroTop2PV.GetChi2() / KFDZeroTop2PV.GetNDF(); + topChi2PerNDF_DZeroBar = KFDZeroBarTop2PV.GetChi2() / KFDZeroBarTop2PV.GetNDF(); + } + // fill candidate table rows rowCandidateBase(collision.globalIndex(), KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ(), errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one + topChi2PerNDF_DZero, topChi2PerNDF_DZeroBar, KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), impactParameter0XY, impactParameter1XY, @@ -326,16 +350,16 @@ struct HfCandidateCreator2Prong { // fill histograms if (fillHistograms) { - Float_t massDZero = 0., errMassDZero = 0.; - Float_t massDZeroBar = 0., errMassDZeroBar = 0.; - if (!KFDZero.GetMass(massDZero, errMassDZero)) - hMass2->Fill(massDZero); - if (KFDZero.GetMass(massDZero, errMassDZero)) - hKFMassFailed->Fill(0); - if (!KFDZeroBar.GetMass(massDZeroBar, errMassDZeroBar)) - hMass2->Fill(massDZeroBar); - if (KFDZeroBar.GetMass(massDZeroBar, errMassDZeroBar)) - hKFMassFailed->Fill(1); + hMass2->Fill(KFDZero.GetMass()); + hMass2->Fill(KFDZeroBar.GetMass()); + if (useKFTop2PV) { + Float_t massDZeroTop = 0., errMassDZeroTop = 0.; + Float_t massDZeroBarTop = 0., errMassDZeroBarTop = 0.; + if (KFDZeroTop2PV.GetMass(massDZeroTop, errMassDZeroTop)) + hKFMassFailed->Fill(KFDZero.GetMass(), 0); // Set topological constraint to DZero failed + if (KFDZeroBarTop2PV.GetMass(massDZeroBarTop, errMassDZeroBarTop)) + hKFMassFailed->Fill(KFDZeroBar.GetMass(), 1); // Set topological constraint to DZeroBar failed + } } } } diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index d6dfbd28d7b..b58698c8e7d 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -102,6 +102,9 @@ struct HfCandidateSelectorD0 { } // candidate DCA // if (candidate.chi2PCA() > cuts[pTBin][1]) return false; + // candidate chi2 + // if (candidate.kfTopChi2OverNDF_DZero() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; + // if (candidate.kfTopChi2OverNDF_DZeroBar() > cuts->get(pTBin, "topological chi2overndf as D0bar")) return false; // decay exponentail law, with tau = beta*gamma*ctau // decay length > ctau retains (1-1/e) From 594ee2d6f575b7a4c4878758d957be9bb77a1ad6 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Mon, 24 Jul 2023 00:21:03 +0200 Subject: [PATCH 04/24] PWGHF: Add topological chi2 to treeCreatorD0ToKPi --- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index 341333b1e0d..224a7eafcdc 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -84,6 +84,8 @@ DECLARE_SOA_TABLE(HfCand2ProngFull, "AOD", "HFCAND2PFull", hf_cand::ErrorDecayLength, hf_cand::ErrorDecayLengthXY, hf_cand::Chi2PCA, + hf_cand::KFTopChi2OverNDF_DZero, + hf_cand::KFTopChi2OverNDF_DZeroBar, full::RSecondaryVertex, full::DecayLength, full::DecayLengthXY, @@ -190,6 +192,8 @@ struct HfTreeCreatorD0ToKPi { candidate.errorDecayLength(), candidate.errorDecayLengthXY(), candidate.chi2PCA(), + candidate.kfTopChi2OverNDF_DZero(), + candidate.kfTopChi2OverNDF_DZeroBar(), candidate.rSecondaryVertex(), candidate.decayLength(), candidate.decayLengthXY(), From 76117d472ef0d43ca8e17e41d42430b382a62d48 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Mon, 24 Jul 2023 11:05:07 +0200 Subject: [PATCH 05/24] PWGHF: Added histograms for D0 and D0bar after topological constraint --- PWGHF/TableProducer/candidateCreator2Prong.cxx | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 5ffc1e8f0cd..278c913d1f1 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -75,7 +75,8 @@ struct HfCandidateCreator2Prong { double bz = 0.; OutputObj hMass2{TH1F("hMass2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; - OutputObj hKFMassFailed{TH2F("hKFMassFailed", ";KF fit failed;entries", 500, 0., 5., 3, -0.5, 2.5)}; // 0 for D0 failed, 1 for D0bar failed + OutputObj hKFMassFailed{TH2F("hKFMassFailed", ";inv. mass with topo fit failed;entries", 500, 0., 5., 3, -0.5, 2.5)}; // 0 for D0 failed, 1 for D0bar failed + OutputObj hMassTop2{TH1F("hMassTop2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 600, -1.0, 5.)}; OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hCovPVYY{TH1F("hCovPVYY", "2-prong candidates;YY element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; @@ -86,8 +87,7 @@ struct HfCandidateCreator2Prong { OutputObj hCovSVZZ{TH1F("hCovSVZZ", "2-prong candidates;ZZ element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hDcaXYProngs{TH2F("hDcaXYProngs", "DCAxy of 2-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{xy}) (#mum);entries", 100, 0., 20., 200, -500., 500.)}; OutputObj hDcaZProngs{TH2F("hDcaZProngs", "DCAz of 2-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{z}) (#mum);entries", 100, 0., 20., 200, -500., 500.)}; - OutputObj hUseKForDCAFitter{TH1F("hUseKForDCAFitter", ";Use KF or DCAFitterN;entries", 2, -0.5, 1.5)}; // 0 for KF, 1 for DCAFitterN - OutputObj hKFMassFailed{TH1F("hKFMassFailed", ";KF geometrical fitting failed;entries", 3, -0.5, 2.5)}; // 0 for D0 failed, 1 for D0bar failed + OutputObj hUseKForDCAFitter{TH1F("hUseKForDCAFitter", ";Use KF or DCAFitterN;entries", 2, -0.5, 1.5)}; // 0 for KF, 1 for DCAFitterN void init(InitContext const&) { @@ -355,6 +355,8 @@ struct HfCandidateCreator2Prong { if (useKFTop2PV) { Float_t massDZeroTop = 0., errMassDZeroTop = 0.; Float_t massDZeroBarTop = 0., errMassDZeroBarTop = 0.; + hMassTop2->Fill(KFDZeroTop2PV.GetMass()); + hMassTop2->Fill(KFDZeroBarTop2PV.GetMass()); if (KFDZeroTop2PV.GetMass(massDZeroTop, errMassDZeroTop)) hKFMassFailed->Fill(KFDZero.GetMass(), 0); // Set topological constraint to DZero failed if (KFDZeroBarTop2PV.GetMass(massDZeroBarTop, errMassDZeroBarTop)) From bbd85d8defdbe69d7293e76b00f33f3451e35ead Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 2 Aug 2023 11:27:43 +0200 Subject: [PATCH 06/24] PWGHF: Add geometrical mass to the output table --- PWGHF/D2H/Tasks/taskD0.cxx | 4 ++++ PWGHF/DataModel/CandidateReconstructionTables.h | 3 +++ PWGHF/TableProducer/candidateCreator2Prong.cxx | 2 ++ PWGHF/TableProducer/candidateSelectorD0.cxx | 4 ++-- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 4 ++++ 5 files changed, 15 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index ed9095d7b36..0174805de61 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -173,6 +173,10 @@ struct HfTaskD0 { auto massD0 = invMassD0ToPiK(candidate); auto massD0bar = invMassD0barToKPi(candidate); + if (candidate.kfGeoMass_DZero() != -999. && candidate.kfGeoMass_DZeroBar() != -999.) { + massD0 = candidate.kfGeoMass_DZero(); + massD0bar = candidate.kfGeoMass_DZeroBar(); + } auto ptCandidate = candidate.pt(); if (candidate.isSelD0() >= selectionFlagD0) { diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 63771869163..a69fa76c2d3 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -255,6 +255,8 @@ DECLARE_SOA_DYNAMIC_COLUMN(RSecondaryVertex, rSecondaryVertex, //! DECLARE_SOA_COLUMN(Chi2PCA, chi2PCA, float); //! sum of (non-weighted) distances of the secondary vertex to its prongs DECLARE_SOA_COLUMN(KFTopChi2OverNDF_DZero, kfTopChi2OverNDF_DZero, float); //! DECLARE_SOA_COLUMN(KFTopChi2OverNDF_DZeroBar, kfTopChi2OverNDF_DZeroBar, float); //! +DECLARE_SOA_COLUMN(KFGeoMass_DZero, kfGeoMass_DZero, float); //! +DECLARE_SOA_COLUMN(KFGeoMass_DZeroBar, kfGeoMass_DZeroBar, float); //! // prong properties DECLARE_SOA_COLUMN(PxProng0, pxProng0, float); //! DECLARE_SOA_COLUMN(PyProng0, pyProng0, float); //! @@ -479,6 +481,7 @@ DECLARE_SOA_TABLE(HfCand2ProngBase, "AOD", "HFCAND2PBASE", //! HFCAND_COLUMNS, hf_cand::KFTopChi2OverNDF_DZero, hf_cand::KFTopChi2OverNDF_DZeroBar, + hf_cand::KFGeoMass_DZero, hf_cand::KFGeoMass_DZeroBar, // 2-prong specific columns hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 278c913d1f1..555f06f3609 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -220,6 +220,7 @@ struct HfCandidateCreator2Prong { errorDecayLength, errorDecayLengthXY, chi2PCA, -999., -999., + -999., -999., pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2], impactParameter0.getY(), impactParameter1.getY(), @@ -341,6 +342,7 @@ struct HfCandidateCreator2Prong { errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one topChi2PerNDF_DZero, topChi2PerNDF_DZeroBar, + KFDZero.GetMass(), KFDZeroBar.GetMass(), KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), impactParameter0XY, impactParameter1XY, diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index b58698c8e7d..12d55ea9225 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -103,8 +103,8 @@ struct HfCandidateSelectorD0 { // candidate DCA // if (candidate.chi2PCA() > cuts[pTBin][1]) return false; // candidate chi2 - // if (candidate.kfTopChi2OverNDF_DZero() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; - // if (candidate.kfTopChi2OverNDF_DZeroBar() > cuts->get(pTBin, "topological chi2overndf as D0bar")) return false; + // if (candidate.kfTopChi2OverNDF_DZero() !=-999. && candidate.kfTopChi2OverNDF_DZero() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; + // if (candidate.kfTopChi2OverNDF_DZeroBar() !=-999. && candidate.kfTopChi2OverNDF_DZeroBar() > cuts->get(pTBin, "topological chi2overndf as D0bar")) return false; // decay exponentail law, with tau = beta*gamma*ctau // decay length > ctau retains (1-1/e) diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index 08bd16d927c..0875522dbe9 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -121,6 +121,8 @@ DECLARE_SOA_TABLE(HfCandD0Fulls, "AOD", "HFCANDD0FULL", hf_cand::Chi2PCA, hf_cand::KFTopChi2OverNDF_DZero, hf_cand::KFTopChi2OverNDF_DZeroBar, + hf_cand::KFGeoMass_DZero, + hf_cand::KFGeoMass_DZeroBar, full::RSecondaryVertex, full::DecayLength, full::DecayLengthXY, @@ -280,6 +282,8 @@ struct HfTreeCreatorD0ToKPi { candidate.chi2PCA(), candidate.kfTopChi2OverNDF_DZero(), candidate.kfTopChi2OverNDF_DZeroBar(), + candidate.kfGeoMass_DZero(), + candidate.kfGeoMass_DZeroBar(), candidate.rSecondaryVertex(), candidate.decayLength(), candidate.decayLengthXY(), From 94279bb60bfb132d6880fcb270be4c57029b4427 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Tue, 8 Aug 2023 22:39:19 +0200 Subject: [PATCH 07/24] PWGHF: Use the extra table to store KFParticle topological chi2 and geometrical mass --- PWGHF/D2H/Tasks/taskD0.cxx | 63 +++++- .../DataModel/CandidateReconstructionTables.h | 14 +- .../TableProducer/candidateCreator2Prong.cxx | 124 ++++++----- PWGHF/TableProducer/candidateSelectorD0.cxx | 49 ++++- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 193 +++++++++++------- 5 files changed, 289 insertions(+), 154 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index 5d89e734450..b8dc38e4f8c 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -28,6 +28,9 @@ using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_2prong; using namespace o2::analysis::hf_cuts_d0_to_pi_k; +constexpr static int useDCAFitterN = 0; +constexpr static int useKFParticle = 1; + /// D0 analysis task struct HfTaskD0 { Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; @@ -41,7 +44,9 @@ struct HfTaskD0 { Configurable> binsPt{"binsPt", std::vector{hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits"}; Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + Partition> selectedD0CandidatesKF = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; Partition> recoFlag2Prong = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; + Partition> recoFlag2ProngKF = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; HistogramRegistry registry{ "registry", @@ -161,9 +166,10 @@ struct HfTaskD0 { registry.add("hDecLengthxyVsPtSig", "2-prong candidates;decay length xy (cm) vs #it{p}_{T} for signal;entries", {HistType::kTH2F, {{800, 0., 4.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } - void process(soa::Join& candidates) + template + void processData(THfCand2Prong& candidates, THfCand2ProngSel& selectedD0CandidatesSets) { - for (auto& candidate : selectedD0Candidates) { + for (auto& candidate : selectedD0CandidatesSets) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -171,9 +177,12 @@ struct HfTaskD0 { continue; } - auto massD0 = invMassD0ToPiK(candidate); - auto massD0bar = invMassD0barToKPi(candidate); - if (candidate.kfGeoMass_DZero() != -999. && candidate.kfGeoMass_DZeroBar() != -999.) { + float massD0, massD0bar; + if constexpr (ReconstructionType == useDCAFitterN) { + massD0 = invMassD0ToPiK(candidate); + massD0bar = invMassD0barToKPi(candidate); + } + if constexpr (ReconstructionType == useKFParticle) { massD0 = candidate.kfGeoMass_DZero(); massD0bar = candidate.kfGeoMass_DZeroBar(); } @@ -219,24 +228,43 @@ struct HfTaskD0 { registry.fill(HIST("hCPAXYFinerBinning"), candidate.cpaXY(), ptCandidate); } } + void processWithDCAFitterN(soa::Join const& candidates) + { + processData<0>(candidates, selectedD0Candidates); + } + void processWithKFParticle(soa::Join const& candidates) + { + processData<1>(candidates, selectedD0CandidatesKF); + } + + PROCESS_SWITCH(HfTaskD0, processWithDCAFitterN, "process taskD0 with DCAFitterN", true); + PROCESS_SWITCH(HfTaskD0, processWithKFParticle, "process taskD0 with KFParticle", false); - void processMc(soa::Join& candidates, + template + void processMc(THfCand2Prong& candidates, THfCand2ProngFlag& recoFlag2ProngSets, soa::Join const& particlesMC, aod::TracksWMc const& tracks) { // MC rec. // Printf("MC Candidates: %d", candidates.size()); - for (auto& candidate : recoFlag2Prong) { + for (auto& candidate : recoFlag2ProngSets) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } if (yCandRecoMax >= 0. && std::abs(yD0(candidate)) > yCandRecoMax) { continue; } - auto massD0 = invMassD0ToPiK(candidate); - auto massD0bar = invMassD0barToKPi(candidate); + float massD0, massD0bar; + if constexpr (ReconstructionType == useDCAFitterN) { + massD0 = invMassD0ToPiK(candidate); + massD0bar = invMassD0barToKPi(candidate); + } + if constexpr (ReconstructionType == useKFParticle) { + massD0 = candidate.kfGeoMass_DZero(); + massD0bar = candidate.kfGeoMass_DZeroBar(); + } if (std::abs(candidate.flagMcMatchRec()) == 1 << DecayType::D0ToPiK) { // Get the corresponding MC particle. - auto indexMother = RecoDecay::getMother(particlesMC, candidate.prong0_as().mcParticle_as>(), pdg::Code::kD0, true); + auto indexMother = RecoDecay::getMother(particlesMC, candidate.template prong0_as().template mcParticle_as>(), pdg::Code::kD0, true); auto particleMother = particlesMC.rawIteratorAt(indexMother); auto ptGen = particleMother.pt(); // gen. level pT auto yGen = RecoDecay::y(array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::getMassPDG(particleMother.pdgCode())); // gen. level y @@ -407,7 +435,20 @@ struct HfTaskD0 { } } - PROCESS_SWITCH(HfTaskD0, processMc, "Process MC", false); + void processMcWithDCAFitterN(soa::Join& candidates, + soa::Join const& particlesMC, aod::TracksWMc const& tracks) + { + processMc<0>(candidates, recoFlag2Prong, particlesMC, tracks); + } + + void processMcWithKFParticle(soa::Join& candidates, + soa::Join const& particlesMC, aod::TracksWMc const& tracks) + { + processMc<1>(candidates, recoFlag2ProngKF, particlesMC, tracks); + } + + PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); + PROCESS_SWITCH(HfTaskD0, processMcWithKFParticle, "Process MC with KFParticle", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 4a89bf4e882..bf73d64903f 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -261,10 +261,9 @@ DECLARE_SOA_COLUMN(ZSecondaryVertex, zSecondaryVertex, float); //! DECLARE_SOA_DYNAMIC_COLUMN(RSecondaryVertex, rSecondaryVertex, //! [](float xVtxS, float yVtxS) -> float { return RecoDecay::sqrtSumOfSquares(xVtxS, yVtxS); }); DECLARE_SOA_COLUMN(Chi2PCA, chi2PCA, float); //! sum of (non-weighted) distances of the secondary vertex to its prongs -DECLARE_SOA_COLUMN(KFTopChi2OverNDF_DZero, kfTopChi2OverNDF_DZero, float); //! -DECLARE_SOA_COLUMN(KFTopChi2OverNDF_DZeroBar, kfTopChi2OverNDF_DZeroBar, float); //! -DECLARE_SOA_COLUMN(KFGeoMass_DZero, kfGeoMass_DZero, float); //! -DECLARE_SOA_COLUMN(KFGeoMass_DZeroBar, kfGeoMass_DZeroBar, float); //! +DECLARE_SOA_COLUMN(KFTopChi2OverNDF, kfTopChi2OverNDF, float); //! +DECLARE_SOA_COLUMN(KFGeoMass_DZero, kfGeoMass_DZero, float); //! +DECLARE_SOA_COLUMN(KFGeoMass_DZeroBar, kfGeoMass_DZeroBar, float); //! // prong properties DECLARE_SOA_COLUMN(PxProng0, pxProng0, float); //! DECLARE_SOA_COLUMN(PyProng0, pyProng0, float); //! @@ -487,9 +486,6 @@ DECLARE_SOA_TABLE(HfCand2ProngBase, "AOD", "HFCAND2PBASE", //! o2::soa::Index<>, // general columns HFCAND_COLUMNS, - hf_cand::KFTopChi2OverNDF_DZero, - hf_cand::KFTopChi2OverNDF_DZeroBar, - hf_cand::KFGeoMass_DZero, hf_cand::KFGeoMass_DZeroBar, // 2-prong specific columns hf_cand::PxProng0, hf_cand::PyProng0, hf_cand::PzProng0, hf_cand::PxProng1, hf_cand::PyProng1, hf_cand::PzProng1, @@ -526,6 +522,10 @@ DECLARE_SOA_EXTENDED_TABLE_USER(HfCand2ProngExt, HfCand2ProngBase, "HFCAND2PEXT" using HfCand2Prong = HfCand2ProngExt; +DECLARE_SOA_TABLE(HfCand2ProngKF, "AOD", "HFCAND2PKF", + hf_cand::KFTopChi2OverNDF, + hf_cand::KFGeoMass_DZero, hf_cand::KFGeoMass_DZeroBar); + // table with results of reconstruction level MC matching DECLARE_SOA_TABLE(HfCand2ProngMcRec, "AOD", "HFCAND2PMCREC", //! hf_cand_2prong::FlagMcMatchRec, diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 61ac3fc5b05..c5883ed947e 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -41,11 +41,10 @@ using namespace o2::aod::hf_cand_2prong; /// Reconstruction of heavy-flavour 2-prong decay candidates struct HfCandidateCreator2Prong { Produces rowCandidateBase; + Produces rowCandidateKF; // vertexing - Configurable useDCAFitterN{"useDCAFitterN", true, "calculate vertex with DCAFitterN"}; - Configurable useKFParticle{"useKFParticle", false, "calculate vertex with KFParticle"}; - Configurable useKFTop2PV{"useKFTop2PV", false, "constraint KFParticle to PV, only if useKFParticle is true"}; + Configurable useKFTop2PV{"useKFTop2PV", true, "constraint KFParticle to PV, only if useKFParticle is true"}; Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; @@ -91,9 +90,12 @@ struct HfCandidateCreator2Prong { void init(InitContext const&) { - if (doprocessPvRefit && doprocessNoPvRefit) { + if ((doprocessPvRefit && doprocessNoPvRefit) || (doprocessPvRefitWithKF && doprocessNoPvRefitWithKF)) { LOGP(fatal, "Only one process function between processPvRefit and processNoPvRefit can be enabled at a time."); } + if ((doprocessPvRefit && (doprocessPvRefitWithKF || doprocessNoPvRefitWithKF)) || (doprocessNoPvRefit && (doprocessPvRefitWithKF || doprocessNoPvRefitWithKF))) { + LOGP(fatal, "Only one process function between DCAFitterN and KFParticle can be enabled at a time."); + } ccdb->setURL(ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -101,39 +103,28 @@ struct HfCandidateCreator2Prong { runNumber = 0; } - template + template void runCreator2Prong(aod::Collisions const& collisions, Cand const& rowsTrackIndexProng2, - aod::TracksWCov const& tracks, + TTracks const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { // 2-prong vertex fitter o2::vertexing::DCAFitterN<2> df; - if (useDCAFitterN) { - // df.setBz(bz); - df.setPropagateToPCA(propagateToPCA); - df.setMaxR(maxR); - df.setMaxDZIni(maxDZIni); - df.setMinParamChange(minParamChange); - df.setMinRelChi2Change(minRelChi2Change); - df.setUseAbsDCA(useAbsDCA); - df.setWeightedFinalPCA(useWeightedFinalPCA); - hUseKForDCAFitter->Fill(1.0); - } - if (useKFParticle) { - hUseKForDCAFitter->Fill(0.0); - } - if (!useDCAFitterN && !useKFParticle) { - LOGP(fatal, "At least one package between useDCAFitterN and useDCAFitterN should be enabled at a time."); - } - if (useDCAFitterN && useKFParticle) { - LOGP(fatal, "Only one package between useDCAFitterN and useDCAFitterN can be enabled at a time."); - } + // df.setBz(bz); + df.setPropagateToPCA(propagateToPCA); + df.setMaxR(maxR); + df.setMaxDZIni(maxDZIni); + df.setMinParamChange(minParamChange); + df.setMinRelChi2Change(minRelChi2Change); + df.setUseAbsDCA(useAbsDCA); + df.setWeightedFinalPCA(useWeightedFinalPCA); + hUseKForDCAFitter->Fill(1.0); // loop over pairs of track indices for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { - auto track0 = rowTrackIndexProng2.template prong0_as(); - auto track1 = rowTrackIndexProng2.template prong1_as(); + auto track0 = rowTrackIndexProng2.template prong0_as(); + auto track1 = rowTrackIndexProng2.template prong1_as(); auto trackParVarPos1 = getTrackParCov(track0); auto trackParVarNeg1 = getTrackParCov(track1); auto collision = rowTrackIndexProng2.collision(); @@ -151,7 +142,6 @@ struct HfCandidateCreator2Prong { // df.print(); } - if (useDCAFitterN) { df.setBz(bz); // reconstruct the 2-prong secondary vertex @@ -219,8 +209,6 @@ struct HfCandidateCreator2Prong { secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], errorDecayLength, errorDecayLengthXY, chi2PCA, - -999., -999., - -999., -999., pvec0[0], pvec0[1], pvec0[2], pvec1[0], pvec1[1], pvec1[2], impactParameter0.getY(), impactParameter1.getY(), @@ -238,8 +226,35 @@ struct HfCandidateCreator2Prong { hMass2->Fill(massKPi); } } + } + + template + void runCreator2ProngWithKF(aod::Collisions const& collisions, + Cand const& rowsTrackIndexProng2, + TTracks const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + hUseKForDCAFitter->Fill(0.0); - if (useKFParticle) { + for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { + auto track0 = rowTrackIndexProng2.template prong0_as(); + auto track1 = rowTrackIndexProng2.template prong1_as(); + auto trackParVarPos1 = getTrackParCov(track0); + auto trackParVarNeg1 = getTrackParCov(track1); + auto collision = rowTrackIndexProng2.collision(); + + /// 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, lut, 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_t covMatrixPV[6]; KFParticle::SetField(bz); @@ -316,23 +331,11 @@ struct HfCandidateCreator2Prong { auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixSV, phi, 0.)); Float_t topChi2PerNDF_DZero = -999.; - Float_t topChi2PerNDF_DZeroBar = -999.; - - KFParticle KFDZeroTop2PV, KFDZeroBarTop2PV; - KFParticle KFDZeroTop2PV_DV, KFDZeroBarTop2PV_DV; + KFParticle KFDZeroTop2PV; if (useKFTop2PV) { KFDZeroTop2PV = KFDZero; KFDZeroTop2PV.SetProductionVertex(KFPV); - KFDZeroTop2PV_DV = KFDZeroTop2PV; - KFDZeroTop2PV_DV.TransportToDecayVertex(); - - KFDZeroBarTop2PV = KFDZeroBar; - KFDZeroBarTop2PV.SetProductionVertex(KFPV); - KFDZeroBarTop2PV_DV = KFDZeroBarTop2PV; - KFDZeroBarTop2PV_DV.TransportToDecayVertex(); - topChi2PerNDF_DZero = KFDZeroTop2PV.GetChi2() / KFDZeroTop2PV.GetNDF(); - topChi2PerNDF_DZeroBar = KFDZeroBarTop2PV.GetChi2() / KFDZeroBarTop2PV.GetNDF(); } // fill candidate table rows @@ -341,14 +344,14 @@ struct HfCandidateCreator2Prong { KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ(), errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one - topChi2PerNDF_DZero, topChi2PerNDF_DZeroBar, - KFDZero.GetMass(), KFDZeroBar.GetMass(), KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), impactParameter0XY, impactParameter1XY, errImpactParameter0XY, errImpactParameter1XY, rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), rowTrackIndexProng2.hfflag()); + rowCandidateKF(topChi2PerNDF_DZero, + KFDZero.GetMass(), KFDZeroBar.GetMass()); // fill histograms if (fillHistograms) { @@ -356,17 +359,12 @@ struct HfCandidateCreator2Prong { hMass2->Fill(KFDZeroBar.GetMass()); if (useKFTop2PV) { Float_t massDZeroTop = 0., errMassDZeroTop = 0.; - Float_t massDZeroBarTop = 0., errMassDZeroBarTop = 0.; hMassTop2->Fill(KFDZeroTop2PV.GetMass()); - hMassTop2->Fill(KFDZeroBarTop2PV.GetMass()); if (KFDZeroTop2PV.GetMass(massDZeroTop, errMassDZeroTop)) hKFMassFailed->Fill(KFDZero.GetMass(), 0); // Set topological constraint to DZero failed - if (KFDZeroBarTop2PV.GetMass(massDZeroBarTop, errMassDZeroBarTop)) - hKFMassFailed->Fill(KFDZeroBar.GetMass(), 1); // Set topological constraint to DZeroBar failed - } } } - } + } } void processPvRefit(aod::Collisions const& collisions, @@ -388,6 +386,26 @@ struct HfCandidateCreator2Prong { } PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefit, "Run candidate creator without PV refit", true); + + void processPvRefitWithKF(aod::Collisions const& collisions, + soa::Join const& rowsTrackIndexProng2, + soa::Join const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator2ProngWithKF(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + } + + PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefitWithKF, "Run candidate creator with PV refit", false); + + void processNoPvRefitWithKF(aod::Collisions const& collisions, + aod::Hf2Prongs const& rowsTrackIndexProng2, + soa::Join const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) + { + runCreator2ProngWithKF(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + } + + PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithKF, "Run candidate creator without PV refit", false); }; /// Extends the base table with expression columns. diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index 661cb376d00..ec4c558bb87 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -28,6 +28,9 @@ using namespace o2::framework; using namespace o2::aod::hf_cand_2prong; using namespace o2::analysis::hf_cuts_d0_to_pi_k; +constexpr static int useDCAFitterN = 0; +constexpr static int useKFParticle = 1; + /// Struct for applying D0 selection cuts struct HfCandidateSelectorD0 { Produces hfSelD0Candidate; @@ -72,7 +75,7 @@ struct HfCandidateSelectorD0 { /// Conjugate-independent topological cuts /// \param candidate is candidate /// \return true if candidate passes all cuts - template + template bool selectionTopol(const T& candidate) { auto candpT = candidate.pt(); @@ -103,9 +106,13 @@ struct HfCandidateSelectorD0 { } // candidate DCA // if (candidate.chi2PCA() > cuts[pTBin][1]) return false; + // candidate chi2 - // if (candidate.kfTopChi2OverNDF_DZero() !=-999. && candidate.kfTopChi2OverNDF_DZero() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; - // if (candidate.kfTopChi2OverNDF_DZeroBar() !=-999. && candidate.kfTopChi2OverNDF_DZeroBar() > cuts->get(pTBin, "topological chi2overndf as D0bar")) return false; + // if constexpr (ReconstructionType == useKFParticle) { + + // if (candidate.kfTopChi2OverNDF() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; + // return false; + // } // decay exponentail law, with tau = beta*gamma*ctau // decay length > ctau retains (1-1/e) @@ -134,7 +141,7 @@ struct HfCandidateSelectorD0 { /// \param trackKaon is the track with the kaon hypothesis /// \note trackPion = positive and trackKaon = negative for D0 selection and inverse for D0bar /// \return true if candidate passes all cuts for the given Conjugate - template + template bool selectionTopolConjugate(const T1& candidate, const T2& trackPion, const T2& trackKaon) { auto candpT = candidate.pt(); @@ -144,12 +151,21 @@ struct HfCandidateSelectorD0 { } // invariant-mass cut + float massD0, massD0bar; + if constexpr (ReconstructionType == useDCAFitterN) { + massD0 = invMassD0ToPiK(candidate); + massD0bar = invMassD0barToKPi(candidate); + } + if constexpr (ReconstructionType == useKFParticle) { + massD0 = candidate.kfGeoMass_DZero(); + massD0bar = candidate.kfGeoMass_DZeroBar(); + } if (trackPion.sign() > 0) { - if (std::abs(invMassD0ToPiK(candidate) - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { + if (std::abs(massD0 - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { return false; } } else { - if (std::abs(invMassD0barToKPi(candidate) - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { + if (std::abs(massD0bar - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { return false; } } @@ -190,8 +206,8 @@ struct HfCandidateSelectorD0 { return true; } - - void process(aod::HfCand2Prong const& candidates, TracksSel const&) + template + void processSel(THfCand2Prong const& candidates, TracksSel const&) { // looping over 2-prong candidates for (auto& candidate : candidates) { @@ -210,11 +226,11 @@ struct HfCandidateSelectorD0 { } statusHFFlag = 1; - auto trackPos = candidate.prong0_as(); // positive daughter - auto trackNeg = candidate.prong1_as(); // negative daughter + auto trackPos = candidate.template prong0_as(); // positive daughter + auto trackNeg = candidate.template prong1_as(); // negative daughter // conjugate-independent topological selection - if (!selectionTopol(candidate)) { + if (!selectionTopol(candidate)) { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); continue; } @@ -288,6 +304,17 @@ struct HfCandidateSelectorD0 { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); } } + void processWithDCAFitterN(aod::HfCand2Prong const& candidates, TracksSel const& tracks) + { + processSel<0>(candidates, tracks); + } + void processWithKFParticle(soa::Join const& candidates, TracksSel const& tracks) + { + processSel<1>(candidates, tracks); + } + + PROCESS_SWITCH(HfCandidateSelectorD0, processWithDCAFitterN, "process candidates selection with DCAFitterN", false); + PROCESS_SWITCH(HfCandidateSelectorD0, processWithKFParticle, "process candidates selection with KFParticle", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index b0a324e6f37..9831aae66bf 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -28,6 +28,9 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_2prong; +constexpr static int useDCAFitterN = 0; +constexpr static int useKFParticle = 1; + namespace o2::aod { namespace full @@ -119,10 +122,7 @@ DECLARE_SOA_TABLE(HfCandD0Fulls, "AOD", "HFCANDD0FULL", hf_cand::ErrorDecayLength, hf_cand::ErrorDecayLengthXY, hf_cand::Chi2PCA, - hf_cand::KFTopChi2OverNDF_DZero, - hf_cand::KFTopChi2OverNDF_DZeroBar, - hf_cand::KFGeoMass_DZero, - hf_cand::KFGeoMass_DZeroBar, + hf_cand::KFTopChi2OverNDF, full::RSecondaryVertex, full::DecayLength, full::DecayLengthXY, @@ -207,6 +207,7 @@ struct HfTreeCreatorD0ToKPi { using TracksWPid = soa::Join; using SelectedCandidatesMc = soa::Filtered>; + using SelectedCandidatesMcKF = soa::Filtered>; using MatchedGenCandidatesMc = soa::Filtered>; Filter filterSelectCandidates = aod::hf_sel_candidate_d0::isSelD0 >= 1 || aod::hf_sel_candidate_d0::isSelD0bar >= 1; @@ -214,6 +215,8 @@ struct HfTreeCreatorD0ToKPi { Partition reconstructedCandSig = nabs(aod::hf_cand_2prong::flagMcMatchRec) == static_cast(BIT(DecayType::D0ToPiK)); Partition reconstructedCandBkg = nabs(aod::hf_cand_2prong::flagMcMatchRec) != static_cast(BIT(DecayType::D0ToPiK)); + Partition reconstructedCandSigKF = nabs(aod::hf_cand_2prong::flagMcMatchRec) == static_cast(BIT(DecayType::D0ToPiK)); + Partition reconstructedCandBkgKF = nabs(aod::hf_cand_2prong::flagMcMatchRec) != static_cast(BIT(DecayType::D0ToPiK)); void init(InitContext const&) { @@ -233,7 +236,7 @@ struct HfTreeCreatorD0ToKPi { } template - auto fillTable(const T& candidate, const U& prong0, const U& prong1, int candFlag, double invMass, double cosThetaStar, + auto fillTable(const T& candidate, const U& prong0, const U& prong1, int candFlag, double invMass, double cosThetaStar, double topoChi2, double ct, double y, double e, int8_t flagMc, int8_t origin) { if (fillCandidateLiteTable) { @@ -281,10 +284,7 @@ struct HfTreeCreatorD0ToKPi { candidate.errorDecayLength(), candidate.errorDecayLengthXY(), candidate.chi2PCA(), - candidate.kfTopChi2OverNDF_DZero(), - candidate.kfTopChi2OverNDF_DZeroBar(), - candidate.kfGeoMass_DZero(), - candidate.kfGeoMass_DZeroBar(), + topoChi2, candidate.rSecondaryVertex(), candidate.decayLength(), candidate.decayLengthXY(), @@ -334,8 +334,9 @@ struct HfTreeCreatorD0ToKPi { } } + template void processData(aod::Collisions const& collisions, - soa::Filtered> const& candidates, + THfCand2Prong const& candidates, TracksWPid const&, aod::BCs const&) { // Filling event properties @@ -357,25 +358,50 @@ struct HfTreeCreatorD0ToKPi { continue; } } - auto prong0 = candidate.prong0_as(); - auto prong1 = candidate.prong1_as(); + auto prong0 = candidate.template prong0_as(); + auto prong1 = candidate.template prong1_as(); double yD = yD0(candidate); double eD = eD0(candidate); double ctD = ctD0(candidate); - if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), ctD, yD, eD, 0, 0); - } - if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), ctD, yD, eD, 0, 0); + if constexpr (ReconstructionType == useDCAFitterN) { + if (candidate.isSelD0()) { + fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), -999., ctD, yD, eD, 0, 0); + } + if (candidate.isSelD0bar()) { + fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), -999., ctD, yD, eD, 0, 0); + } + } else if constexpr (ReconstructionType == useKFParticle) { + if (candidate.isSelD0()) { + fillTable(candidate, prong0, prong1, 0, candidate.kfGeoMass_DZero(), cosThetaStarD0(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, 0, 0); + } + if (candidate.isSelD0bar()) { + fillTable(candidate, prong0, prong1, 1, candidate.kfGeoMass_DZeroBar(), cosThetaStarD0bar(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, 0, 0); + } } } } - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processData, "Process data", true); + void processDataWithDCAFitterN(aod::Collisions const& collisions, + soa::Filtered> const& candidates, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processData<0>(collisions, candidates, tracks, BCs); + } + void processDataWithKFParticle(aod::Collisions const& collisions, + soa::Filtered> const& candidates, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processData<1>(collisions, candidates, tracks, BCs); + } + + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithDCAFitterN, "Process data with DCAFitterN", true); + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithKFParticle, "Process data with KFParticle", false); + template void processMc(aod::Collisions const& collisions, aod::McCollisions const&, - SelectedCandidatesMc const& candidates, + THfCand2Prong const& candidates, + THfCand2ProngSel const& reconstructedCandSets, MatchedGenCandidatesMc const& particles, TracksWPid const&, aod::BCs const&) { @@ -386,67 +412,36 @@ struct HfTreeCreatorD0ToKPi { } // Filling candidate properties - if (fillOnlySignal) { - if (fillCandidateLiteTable) { - rowCandidateLite.reserve(reconstructedCandSig.size()); - } else { - rowCandidateFull.reserve(reconstructedCandSig.size()); - } - for (const auto& candidate : reconstructedCandSig) { - auto prong0 = candidate.prong0_as(); - auto prong1 = candidate.prong0_as(); - double yD = yD0(candidate); - double eD = eD0(candidate); - double ctD = ctD0(candidate); - if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); - } - if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); + if (fillCandidateLiteTable) { + rowCandidateLite.reserve(reconstructedCandSets.size()); + } else { + rowCandidateFull.reserve(reconstructedCandSets.size()); + } + for (const auto& candidate : reconstructedCandSets) { + if (onlyBkg && downSampleBkgFactor < 1.) { + float pseudoRndm = candidate.ptProng0() * 1000. - (int64_t)(candidate.ptProng0() * 1000); + if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { + continue; } } - } else if (fillOnlyBackground) { - if (fillCandidateLiteTable) { - rowCandidateLite.reserve(reconstructedCandBkg.size()); - } else { - rowCandidateFull.reserve(reconstructedCandBkg.size()); - } - for (const auto& candidate : reconstructedCandBkg) { - if (downSampleBkgFactor < 1.) { - float pseudoRndm = candidate.ptProng0() * 1000. - (int64_t)(candidate.ptProng0() * 1000); - if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { - continue; - } - } - auto prong0 = candidate.prong0_as(); - auto prong1 = candidate.prong0_as(); - double yD = yD0(candidate); - double eD = eD0(candidate); - double ctD = ctD0(candidate); + auto prong0 = candidate.template prong0_as(); + auto prong1 = candidate.template prong1_as(); + double yD = yD0(candidate); + double eD = eD0(candidate); + double ctD = ctD0(candidate); + if constexpr (ReconstructionType == useDCAFitterN) { if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); + fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), -999., ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); } if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); + fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), -999., ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); } - } - } else { - if (fillCandidateLiteTable) { - rowCandidateLite.reserve(candidates.size()); - } else { - rowCandidateFull.reserve(candidates.size()); - } - for (const auto& candidate : candidates) { - auto prong0 = candidate.prong0_as(); - auto prong1 = candidate.prong0_as(); - double yD = yD0(candidate); - double eD = eD0(candidate); - double ctD = ctD0(candidate); + } else if constexpr (ReconstructionType == useKFParticle) { if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); + fillTable(candidate, prong0, prong1, 0, candidate.kfGeoMass_DZero(), cosThetaStarD0(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); } if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); + fillTable(candidate, prong0, prong1, 1, candidate.kfGeoMass_DZeroBar(), cosThetaStarD0bar(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); } } } @@ -468,7 +463,61 @@ struct HfTreeCreatorD0ToKPi { } } - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMc, "Process MC", false); + void processMcWithDCAFitterOnlySig(aod::Collisions const& collisions, + aod::McCollisions const& McCollisions, + SelectedCandidatesMc const& candidates, + MatchedGenCandidatesMc const& particles, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processMc<0>(collisions, McCollisions, candidates, reconstructedCandSig, particles, tracks, BCs); + } + void processMcWithDCAFitterOnlyBkg(aod::Collisions const& collisions, + aod::McCollisions const& McCollisions, + SelectedCandidatesMc const& candidates, + MatchedGenCandidatesMc const& particles, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processMc<0, true>(collisions, McCollisions, candidates, reconstructedCandBkg, particles, tracks, BCs); + } + void processMcWithDCAFitterAll(aod::Collisions const& collisions, + aod::McCollisions const& McCollisions, + SelectedCandidatesMc const& candidates, + MatchedGenCandidatesMc const& particles, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processMc<0>(collisions, McCollisions, candidates, candidates, particles, tracks, BCs); + } + void processMcWithKFParticleOnlySig(aod::Collisions const& collisions, + aod::McCollisions const& McCollisions, + SelectedCandidatesMcKF const& candidates, + MatchedGenCandidatesMc const& particles, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processMc<1>(collisions, McCollisions, candidates, reconstructedCandSigKF, particles, tracks, BCs); + } + void processMcWithKFParticleOnlyBkg(aod::Collisions const& collisions, + aod::McCollisions const& McCollisions, + SelectedCandidatesMcKF const& candidates, + MatchedGenCandidatesMc const& particles, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processMc<1, true>(collisions, McCollisions, candidates, reconstructedCandBkgKF, particles, tracks, BCs); + } + void processMcWithKFParticleAll(aod::Collisions const& collisions, + aod::McCollisions const& McCollisions, + SelectedCandidatesMcKF const& candidates, + MatchedGenCandidatesMc const& particles, + TracksWPid const& tracks, aod::BCs const& BCs) + { + processMc<1>(collisions, McCollisions, candidates, candidates, particles, tracks, BCs); + } + + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlySig, "Process MC with DCAFitterN only for signals", false); + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlyBkg, "Process MC with DCAFitterN only for background", false); + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterAll, "Process MC with DCAFitterN", false); + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlySig, "Process MC with KFParticle only for signals", false); + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlyBkg, "Process MC with KFParticle only for background", false); + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleAll, "Process MC with KFParticle", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From ee6315d0cdd8e62156186cebed9068f882eedb23 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Tue, 8 Aug 2023 23:36:43 +0200 Subject: [PATCH 08/24] PWGHF: git-cland-format for candidateCreator2Prong --- PWGHF/TableProducer/candidateCreator2Prong.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index c5883ed947e..ecd9e275542 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -514,4 +514,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc, TaskName{"hf-candidate-creator-2prong"}), adaptAnalysisTask(cfgc, TaskName{"hf-candidate-creator-2prong-expressions"})}; -} +} \ No newline at end of file From d6214c00bdaf6da2fbd41d1af9cb2c324a5a85d4 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Tue, 8 Aug 2023 23:56:40 +0200 Subject: [PATCH 09/24] PWGHF:git-clang-format --- PWGHF/TableProducer/candidateCreator2Prong.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index ecd9e275542..c5883ed947e 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -514,4 +514,4 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) return WorkflowSpec{ adaptAnalysisTask(cfgc, TaskName{"hf-candidate-creator-2prong"}), adaptAnalysisTask(cfgc, TaskName{"hf-candidate-creator-2prong-expressions"})}; -} \ No newline at end of file +} From f2068328b5cbafc5751fc428eed2744d1df7177d Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Fri, 11 Aug 2023 10:26:12 +0200 Subject: [PATCH 10/24] PWGHF: candidateCreator2Prong.cxx with git-clang-format --- .../TableProducer/candidateCreator2Prong.cxx | 414 +++++++++--------- 1 file changed, 207 insertions(+), 207 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index c5883ed947e..8a87d946d2f 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -142,90 +142,90 @@ struct HfCandidateCreator2Prong { // df.print(); } - df.setBz(bz); + df.setBz(bz); - // reconstruct the 2-prong secondary vertex - if (df.process(trackParVarPos1, trackParVarNeg1) == 0) { - continue; - } - const auto& secondaryVertex = df.getPCACandidate(); - auto chi2PCA = df.getChi2AtPCACandidate(); - auto covMatrixPCA = df.calcPCACovMatrixFlat(); - hCovSVXX->Fill(covMatrixPCA[0]); // FIXME: Calculation of errorDecayLength(XY) gives wrong values without this line. - hCovSVYY->Fill(covMatrixPCA[2]); - hCovSVXZ->Fill(covMatrixPCA[3]); - hCovSVZZ->Fill(covMatrixPCA[5]); - auto trackParVar0 = df.getTrack(0); - auto trackParVar1 = df.getTrack(1); - - // get track momenta - array pvec0; - array pvec1; - trackParVar0.getPxPyPzGlo(pvec0); - trackParVar1.getPxPyPzGlo(pvec1); - - // get track impact parameters - // This modifies track momenta! - auto primaryVertex = getPrimaryVertex(collision); - auto covMatrixPV = primaryVertex.getCov(); - if constexpr (doPvRefit) { - /// use PV refit - /// Using it in the rowCandidateBase all dynamic columns shall take it into account - // coordinates - primaryVertex.setX(rowTrackIndexProng2.pvRefitX()); - primaryVertex.setY(rowTrackIndexProng2.pvRefitY()); - primaryVertex.setZ(rowTrackIndexProng2.pvRefitZ()); - // covariance matrix - primaryVertex.setSigmaX2(rowTrackIndexProng2.pvRefitSigmaX2()); - primaryVertex.setSigmaXY(rowTrackIndexProng2.pvRefitSigmaXY()); - primaryVertex.setSigmaY2(rowTrackIndexProng2.pvRefitSigmaY2()); - primaryVertex.setSigmaXZ(rowTrackIndexProng2.pvRefitSigmaXZ()); - primaryVertex.setSigmaYZ(rowTrackIndexProng2.pvRefitSigmaYZ()); - primaryVertex.setSigmaZ2(rowTrackIndexProng2.pvRefitSigmaZ2()); - covMatrixPV = primaryVertex.getCov(); - } - hCovPVXX->Fill(covMatrixPV[0]); - hCovPVYY->Fill(covMatrixPV[2]); - hCovPVXZ->Fill(covMatrixPV[3]); - hCovPVZZ->Fill(covMatrixPV[5]); - o2::dataformats::DCA impactParameter0; - o2::dataformats::DCA impactParameter1; - trackParVar0.propagateToDCA(primaryVertex, bz, &impactParameter0); - trackParVar1.propagateToDCA(primaryVertex, bz, &impactParameter1); - hDcaXYProngs->Fill(track0.pt(), impactParameter0.getY() * toMicrometers); - hDcaXYProngs->Fill(track1.pt(), impactParameter1.getY() * toMicrometers); - hDcaZProngs->Fill(track0.pt(), impactParameter0.getZ() * toMicrometers); - hDcaZProngs->Fill(track1.pt(), impactParameter1.getZ() * toMicrometers); - - // get uncertainty of the decay length - double phi, theta; - getPointDirection(array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); - auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); - auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); - - // fill candidate table rows - rowCandidateBase(collision.globalIndex(), - primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), - secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], - errorDecayLength, errorDecayLengthXY, - chi2PCA, - pvec0[0], pvec0[1], pvec0[2], - pvec1[0], pvec1[1], pvec1[2], - impactParameter0.getY(), impactParameter1.getY(), - std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), - rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), - rowTrackIndexProng2.hfflag()); - - // fill histograms - if (fillHistograms) { - // calculate invariant masses - auto arrayMomenta = array{pvec0, pvec1}; - massPiK = RecoDecay::m(arrayMomenta, array{massPi, massK}); - massKPi = RecoDecay::m(arrayMomenta, array{massK, massPi}); - hMass2->Fill(massPiK); - hMass2->Fill(massKPi); - } + // reconstruct the 2-prong secondary vertex + if (df.process(trackParVarPos1, trackParVarNeg1) == 0) { + continue; } + const auto& secondaryVertex = df.getPCACandidate(); + auto chi2PCA = df.getChi2AtPCACandidate(); + auto covMatrixPCA = df.calcPCACovMatrixFlat(); + hCovSVXX->Fill(covMatrixPCA[0]); // FIXME: Calculation of errorDecayLength(XY) gives wrong values without this line. + hCovSVYY->Fill(covMatrixPCA[2]); + hCovSVXZ->Fill(covMatrixPCA[3]); + hCovSVZZ->Fill(covMatrixPCA[5]); + auto trackParVar0 = df.getTrack(0); + auto trackParVar1 = df.getTrack(1); + + // get track momenta + array pvec0; + array pvec1; + trackParVar0.getPxPyPzGlo(pvec0); + trackParVar1.getPxPyPzGlo(pvec1); + + // get track impact parameters + // This modifies track momenta! + auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + if constexpr (doPvRefit) { + /// use PV refit + /// Using it in the rowCandidateBase all dynamic columns shall take it into account + // coordinates + primaryVertex.setX(rowTrackIndexProng2.pvRefitX()); + primaryVertex.setY(rowTrackIndexProng2.pvRefitY()); + primaryVertex.setZ(rowTrackIndexProng2.pvRefitZ()); + // covariance matrix + primaryVertex.setSigmaX2(rowTrackIndexProng2.pvRefitSigmaX2()); + primaryVertex.setSigmaXY(rowTrackIndexProng2.pvRefitSigmaXY()); + primaryVertex.setSigmaY2(rowTrackIndexProng2.pvRefitSigmaY2()); + primaryVertex.setSigmaXZ(rowTrackIndexProng2.pvRefitSigmaXZ()); + primaryVertex.setSigmaYZ(rowTrackIndexProng2.pvRefitSigmaYZ()); + primaryVertex.setSigmaZ2(rowTrackIndexProng2.pvRefitSigmaZ2()); + covMatrixPV = primaryVertex.getCov(); + } + hCovPVXX->Fill(covMatrixPV[0]); + hCovPVYY->Fill(covMatrixPV[2]); + hCovPVXZ->Fill(covMatrixPV[3]); + hCovPVZZ->Fill(covMatrixPV[5]); + o2::dataformats::DCA impactParameter0; + o2::dataformats::DCA impactParameter1; + trackParVar0.propagateToDCA(primaryVertex, bz, &impactParameter0); + trackParVar1.propagateToDCA(primaryVertex, bz, &impactParameter1); + hDcaXYProngs->Fill(track0.pt(), impactParameter0.getY() * toMicrometers); + hDcaXYProngs->Fill(track1.pt(), impactParameter1.getY() * toMicrometers); + hDcaZProngs->Fill(track0.pt(), impactParameter0.getZ() * toMicrometers); + hDcaZProngs->Fill(track1.pt(), impactParameter1.getZ() * toMicrometers); + + // get uncertainty of the decay length + double phi, theta; + getPointDirection(array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + + // fill candidate table rows + rowCandidateBase(collision.globalIndex(), + primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ(), + secondaryVertex[0], secondaryVertex[1], secondaryVertex[2], + errorDecayLength, errorDecayLengthXY, + chi2PCA, + pvec0[0], pvec0[1], pvec0[2], + pvec1[0], pvec1[1], pvec1[2], + impactParameter0.getY(), impactParameter1.getY(), + std::sqrt(impactParameter0.getSigmaY2()), std::sqrt(impactParameter1.getSigmaY2()), + rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), + rowTrackIndexProng2.hfflag()); + + // fill histograms + if (fillHistograms) { + // calculate invariant masses + auto arrayMomenta = array{pvec0, pvec1}; + massPiK = RecoDecay::m(arrayMomenta, array{massPi, massK}); + massKPi = RecoDecay::m(arrayMomenta, array{massK, massPi}); + hMass2->Fill(massPiK); + hMass2->Fill(massKPi); + } + } } template @@ -234,137 +234,137 @@ struct HfCandidateCreator2Prong { TTracks const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { - hUseKForDCAFitter->Fill(0.0); - - for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { - auto track0 = rowTrackIndexProng2.template prong0_as(); - auto track1 = rowTrackIndexProng2.template prong1_as(); - auto trackParVarPos1 = getTrackParCov(track0); - auto trackParVarNeg1 = getTrackParCov(track1); - auto collision = rowTrackIndexProng2.collision(); - - /// 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, lut, 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_t 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(rowTrackIndexProng2.pvRefitX(), rowTrackIndexProng2.pvRefitY(), rowTrackIndexProng2.pvRefitZ()); - // covariance matrix - kfpVertex.SetCovarianceMatrix(rowTrackIndexProng2.pvRefitSigmaX2(), rowTrackIndexProng2.pvRefitSigmaXY(), rowTrackIndexProng2.pvRefitSigmaY2(), rowTrackIndexProng2.pvRefitSigmaXZ(), rowTrackIndexProng2.pvRefitSigmaYZ(), rowTrackIndexProng2.pvRefitSigmaZ2()); - } - kfpVertex.GetCovarianceMatrix(covMatrixPV); - KFParticle KFPV(kfpVertex); - hCovPVXX->Fill(covMatrixPV[0]); - hCovPVYY->Fill(covMatrixPV[2]); - hCovPVXZ->Fill(covMatrixPV[3]); - hCovPVZZ->Fill(covMatrixPV[5]); - - KFPTrack kfpTrackPosPi; - KFPTrack kfpTrackNegPi; - KFPTrack kfpTrackPosKa; - KFPTrack kfpTrackNegKa; - - if (track0.sign() == 1 && track1.sign() == -1) { - kfpTrackPosPi = createKFPTrackFromTrack(track0); - kfpTrackNegKa = createKFPTrackFromTrack(track1); - kfpTrackPosKa = createKFPTrackFromTrack(track0); - kfpTrackNegPi = createKFPTrackFromTrack(track1); - } + hUseKForDCAFitter->Fill(0.0); - KFParticle KFPosPion(kfpTrackPosPi, 211); - KFParticle KFNegPion(kfpTrackNegPi, 211); - KFParticle KFPosKaon(kfpTrackPosKa, 321); - KFParticle KFNegKaon(kfpTrackNegKa, 321); + for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { + auto track0 = rowTrackIndexProng2.template prong0_as(); + auto track1 = rowTrackIndexProng2.template prong1_as(); + auto trackParVarPos1 = getTrackParCov(track0); + auto trackParVarNeg1 = getTrackParCov(track1); + auto collision = rowTrackIndexProng2.collision(); - float impactParameter0XY = 0., errImpactParameter0XY = 0., impactParameter1XY = 0., errImpactParameter1XY = 0.; - if (!KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { - hDcaXYProngs->Fill(track0.pt(), impactParameter0XY * toMicrometers); - hDcaZProngs->Fill(track0.pt(), sqrt(KFPosPion.GetDistanceFromVertex(KFPV) * KFPosPion.GetDistanceFromVertex(KFPV) - impactParameter0XY * impactParameter0XY) * toMicrometers); - } - if (KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { - hDcaXYProngs->Fill(track0.pt(), -999.); - hDcaZProngs->Fill(track0.pt(), -999.); - } - if (!KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { - hDcaXYProngs->Fill(track1.pt(), impactParameter1XY * toMicrometers); - hDcaZProngs->Fill(track1.pt(), sqrt(KFNegPion.GetDistanceFromVertex(KFPV) * KFNegPion.GetDistanceFromVertex(KFPV) - impactParameter1XY * impactParameter1XY) * toMicrometers); - } - if (KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { - hDcaXYProngs->Fill(track1.pt(), -999.); - hDcaZProngs->Fill(track1.pt(), -999.); - } + /// 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, lut, 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_t 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(rowTrackIndexProng2.pvRefitX(), rowTrackIndexProng2.pvRefitY(), rowTrackIndexProng2.pvRefitZ()); + // covariance matrix + kfpVertex.SetCovarianceMatrix(rowTrackIndexProng2.pvRefitSigmaX2(), rowTrackIndexProng2.pvRefitSigmaXY(), rowTrackIndexProng2.pvRefitSigmaY2(), rowTrackIndexProng2.pvRefitSigmaXZ(), rowTrackIndexProng2.pvRefitSigmaYZ(), rowTrackIndexProng2.pvRefitSigmaZ2()); + } + kfpVertex.GetCovarianceMatrix(covMatrixPV); + KFParticle KFPV(kfpVertex); + hCovPVXX->Fill(covMatrixPV[0]); + hCovPVYY->Fill(covMatrixPV[2]); + hCovPVXZ->Fill(covMatrixPV[3]); + hCovPVZZ->Fill(covMatrixPV[5]); + + KFPTrack kfpTrackPosPi; + KFPTrack kfpTrackNegPi; + KFPTrack kfpTrackPosKa; + KFPTrack kfpTrackNegKa; + + if (track0.sign() == 1 && track1.sign() == -1) { + kfpTrackPosPi = createKFPTrackFromTrack(track0); + kfpTrackNegKa = createKFPTrackFromTrack(track1); + kfpTrackPosKa = createKFPTrackFromTrack(track0); + kfpTrackNegPi = createKFPTrackFromTrack(track1); + } - KFParticle KFDZero; - const KFParticle* D0Daughters[2] = {&KFPosPion, &KFNegKaon}; - KFDZero.SetConstructMethod(2); - KFDZero.Construct(D0Daughters, 2); - KFParticle KFDZeroBar; - const KFParticle* D0BarDaughters[2] = {&KFNegPion, &KFPosKaon}; - KFDZeroBar.SetConstructMethod(2); - KFDZeroBar.Construct(D0BarDaughters, 2); - - hCovSVXX->Fill(KFDZero.Covariance(0, 0)); - hCovSVYY->Fill(KFDZero.Covariance(1, 1)); - hCovSVXZ->Fill(KFDZero.Covariance(2, 0)); - hCovSVZZ->Fill(KFDZero.Covariance(2, 2)); - auto covMatrixSV = KFDZero.CovarianceMatrix(); - - double phi, theta; - getPointDirection(array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, array{KFDZero.GetX(), KFDZero.GetY(), KFDZero.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.)); - - Float_t topChi2PerNDF_DZero = -999.; - KFParticle KFDZeroTop2PV; - if (useKFTop2PV) { - KFDZeroTop2PV = KFDZero; - KFDZeroTop2PV.SetProductionVertex(KFPV); - topChi2PerNDF_DZero = KFDZeroTop2PV.GetChi2() / KFDZeroTop2PV.GetNDF(); - } + KFParticle KFPosPion(kfpTrackPosPi, 211); + KFParticle KFNegPion(kfpTrackNegPi, 211); + KFParticle KFPosKaon(kfpTrackPosKa, 321); + KFParticle KFNegKaon(kfpTrackNegKa, 321); - // fill candidate table rows - rowCandidateBase(collision.globalIndex(), - KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), - KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ(), - errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one - KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one - KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), - KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), - impactParameter0XY, impactParameter1XY, - errImpactParameter0XY, errImpactParameter1XY, - rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), - rowTrackIndexProng2.hfflag()); - rowCandidateKF(topChi2PerNDF_DZero, - KFDZero.GetMass(), KFDZeroBar.GetMass()); - - // fill histograms - if (fillHistograms) { - hMass2->Fill(KFDZero.GetMass()); - hMass2->Fill(KFDZeroBar.GetMass()); - if (useKFTop2PV) { - Float_t massDZeroTop = 0., errMassDZeroTop = 0.; - hMassTop2->Fill(KFDZeroTop2PV.GetMass()); - if (KFDZeroTop2PV.GetMass(massDZeroTop, errMassDZeroTop)) - hKFMassFailed->Fill(KFDZero.GetMass(), 0); // Set topological constraint to DZero failed - } + float impactParameter0XY = 0., errImpactParameter0XY = 0., impactParameter1XY = 0., errImpactParameter1XY = 0.; + if (!KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + hDcaXYProngs->Fill(track0.pt(), impactParameter0XY * toMicrometers); + hDcaZProngs->Fill(track0.pt(), sqrt(KFPosPion.GetDistanceFromVertex(KFPV) * KFPosPion.GetDistanceFromVertex(KFPV) - impactParameter0XY * impactParameter0XY) * toMicrometers); + } + if (KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + hDcaXYProngs->Fill(track0.pt(), -999.); + hDcaZProngs->Fill(track0.pt(), -999.); } + if (!KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + hDcaXYProngs->Fill(track1.pt(), impactParameter1XY * toMicrometers); + hDcaZProngs->Fill(track1.pt(), sqrt(KFNegPion.GetDistanceFromVertex(KFPV) * KFNegPion.GetDistanceFromVertex(KFPV) - impactParameter1XY * impactParameter1XY) * toMicrometers); + } + if (KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + hDcaXYProngs->Fill(track1.pt(), -999.); + hDcaZProngs->Fill(track1.pt(), -999.); + } + + KFParticle KFDZero; + const KFParticle* D0Daughters[2] = {&KFPosPion, &KFNegKaon}; + KFDZero.SetConstructMethod(2); + KFDZero.Construct(D0Daughters, 2); + KFParticle KFDZeroBar; + const KFParticle* D0BarDaughters[2] = {&KFNegPion, &KFPosKaon}; + KFDZeroBar.SetConstructMethod(2); + KFDZeroBar.Construct(D0BarDaughters, 2); + + hCovSVXX->Fill(KFDZero.Covariance(0, 0)); + hCovSVYY->Fill(KFDZero.Covariance(1, 1)); + hCovSVXZ->Fill(KFDZero.Covariance(2, 0)); + hCovSVZZ->Fill(KFDZero.Covariance(2, 2)); + auto covMatrixSV = KFDZero.CovarianceMatrix(); + + double phi, theta; + getPointDirection(array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, array{KFDZero.GetX(), KFDZero.GetY(), KFDZero.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.)); + + Float_t topChi2PerNDF_DZero = -999.; + KFParticle KFDZeroTop2PV; + if (useKFTop2PV) { + KFDZeroTop2PV = KFDZero; + KFDZeroTop2PV.SetProductionVertex(KFPV); + topChi2PerNDF_DZero = KFDZeroTop2PV.GetChi2() / KFDZeroTop2PV.GetNDF(); } + + // fill candidate table rows + rowCandidateBase(collision.globalIndex(), + KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), + KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ(), + errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one + KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one + KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), + KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), + impactParameter0XY, impactParameter1XY, + errImpactParameter0XY, errImpactParameter1XY, + rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), + rowTrackIndexProng2.hfflag()); + rowCandidateKF(topChi2PerNDF_DZero, + KFDZero.GetMass(), KFDZeroBar.GetMass()); + + // fill histograms + if (fillHistograms) { + hMass2->Fill(KFDZero.GetMass()); + hMass2->Fill(KFDZeroBar.GetMass()); + if (useKFTop2PV) { + Float_t massDZeroTop = 0., errMassDZeroTop = 0.; + hMassTop2->Fill(KFDZeroTop2PV.GetMass()); + if (KFDZeroTop2PV.GetMass(massDZeroTop, errMassDZeroTop)) + hKFMassFailed->Fill(KFDZero.GetMass(), 0); // Set topological constraint to DZero failed + } + } + } } void processPvRefit(aod::Collisions const& collisions, From adf51db21f45a0e0d8d472695cc787fd42b9b3a3 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 16 Aug 2023 09:42:04 +0200 Subject: [PATCH 11/24] PWGHF: optimize the code to add Kf to D0 --- PWGHF/D2H/Tasks/taskD0.cxx | 54 +++---- .../DataModel/CandidateReconstructionTables.h | 14 +- PWGHF/TableProducer/CMakeLists.txt | 2 +- .../TableProducer/candidateCreator2Prong.cxx | 144 +++++++++--------- PWGHF/TableProducer/candidateSelectorD0.cxx | 26 ++-- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 117 +++++++------- 6 files changed, 171 insertions(+), 186 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index b8dc38e4f8c..b6ed8eb22af 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -28,9 +28,6 @@ using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_2prong; using namespace o2::analysis::hf_cuts_d0_to_pi_k; -constexpr static int useDCAFitterN = 0; -constexpr static int useKFParticle = 1; - /// D0 analysis task struct HfTaskD0 { Configurable selectionFlagD0{"selectionFlagD0", 1, "Selection Flag for D0"}; @@ -166,8 +163,8 @@ struct HfTaskD0 { registry.add("hDecLengthxyVsPtSig", "2-prong candidates;decay length xy (cm) vs #it{p}_{T} for signal;entries", {HistType::kTH2F, {{800, 0., 4.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } - template - void processData(THfCand2Prong& candidates, THfCand2ProngSel& selectedD0CandidatesSets) + template + void processData(THfCand2ProngSel& selectedD0CandidatesSets) { for (auto& candidate : selectedD0CandidatesSets) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { @@ -178,14 +175,13 @@ struct HfTaskD0 { } float massD0, massD0bar; - if constexpr (ReconstructionType == useDCAFitterN) { + if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + massD0 = candidate.kfGeoMassD0(); + massD0bar = candidate.kfGeoMassD0bar(); + } else { massD0 = invMassD0ToPiK(candidate); massD0bar = invMassD0barToKPi(candidate); } - if constexpr (ReconstructionType == useKFParticle) { - massD0 = candidate.kfGeoMass_DZero(); - massD0bar = candidate.kfGeoMass_DZeroBar(); - } auto ptCandidate = candidate.pt(); if (candidate.isSelD0() >= selectionFlagD0) { @@ -228,20 +224,20 @@ struct HfTaskD0 { registry.fill(HIST("hCPAXYFinerBinning"), candidate.cpaXY(), ptCandidate); } } - void processWithDCAFitterN(soa::Join const& candidates) + void processWithDCAFitterN(soa::Join const&) { - processData<0>(candidates, selectedD0Candidates); + processData<0>(selectedD0Candidates); } - void processWithKFParticle(soa::Join const& candidates) + PROCESS_SWITCH(HfTaskD0, processWithDCAFitterN, "process taskD0 with DCAFitterN", true); + + void processWithKFParticle(soa::Join const&) { - processData<1>(candidates, selectedD0CandidatesKF); + processData<1>(selectedD0CandidatesKF); } - - PROCESS_SWITCH(HfTaskD0, processWithDCAFitterN, "process taskD0 with DCAFitterN", true); PROCESS_SWITCH(HfTaskD0, processWithKFParticle, "process taskD0 with KFParticle", false); - template - void processMc(THfCand2Prong& candidates, THfCand2ProngFlag& recoFlag2ProngSets, + template + void processMc(THfCand2ProngFlag& recoFlag2ProngSets, soa::Join const& particlesMC, aod::TracksWMc const& tracks) { // MC rec. @@ -254,14 +250,13 @@ struct HfTaskD0 { continue; } float massD0, massD0bar; - if constexpr (ReconstructionType == useDCAFitterN) { + if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + massD0 = candidate.kfGeoMassD0(); + massD0bar = candidate.kfGeoMassD0bar(); + } else { massD0 = invMassD0ToPiK(candidate); massD0bar = invMassD0barToKPi(candidate); } - if constexpr (ReconstructionType == useKFParticle) { - massD0 = candidate.kfGeoMass_DZero(); - massD0bar = candidate.kfGeoMass_DZeroBar(); - } if (std::abs(candidate.flagMcMatchRec()) == 1 << DecayType::D0ToPiK) { // Get the corresponding MC particle. auto indexMother = RecoDecay::getMother(particlesMC, candidate.template prong0_as().template mcParticle_as>(), pdg::Code::kD0, true); @@ -435,19 +430,16 @@ struct HfTaskD0 { } } - void processMcWithDCAFitterN(soa::Join& candidates, - soa::Join const& particlesMC, aod::TracksWMc const& tracks) + void processMcWithDCAFitterN(soa::Join const& particlesMC, aod::TracksWMc const& tracks) { - processMc<0>(candidates, recoFlag2Prong, particlesMC, tracks); + processMc<0>(recoFlag2Prong, particlesMC, tracks); } + PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); - void processMcWithKFParticle(soa::Join& candidates, - soa::Join const& particlesMC, aod::TracksWMc const& tracks) + void processMcWithKFParticle(soa::Join const& particlesMC, aod::TracksWMc const& tracks) { - processMc<1>(candidates, recoFlag2ProngKF, particlesMC, tracks); + processMc<1>(recoFlag2ProngKF, particlesMC, tracks); } - - PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); PROCESS_SWITCH(HfTaskD0, processMcWithKFParticle, "Process MC with KFParticle", false); }; diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index bf73d64903f..a4fa9044ed1 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -261,9 +261,6 @@ DECLARE_SOA_COLUMN(ZSecondaryVertex, zSecondaryVertex, float); //! DECLARE_SOA_DYNAMIC_COLUMN(RSecondaryVertex, rSecondaryVertex, //! [](float xVtxS, float yVtxS) -> float { return RecoDecay::sqrtSumOfSquares(xVtxS, yVtxS); }); DECLARE_SOA_COLUMN(Chi2PCA, chi2PCA, float); //! sum of (non-weighted) distances of the secondary vertex to its prongs -DECLARE_SOA_COLUMN(KFTopChi2OverNDF, kfTopChi2OverNDF, float); //! -DECLARE_SOA_COLUMN(KFGeoMass_DZero, kfGeoMass_DZero, float); //! -DECLARE_SOA_COLUMN(KFGeoMass_DZeroBar, kfGeoMass_DZeroBar, float); //! // prong properties DECLARE_SOA_COLUMN(PxProng0, pxProng0, float); //! DECLARE_SOA_COLUMN(PyProng0, pyProng0, float); //! @@ -343,6 +340,9 @@ DECLARE_SOA_DYNAMIC_COLUMN(Ct, ct, //! [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz, double m) -> float { return RecoDecay::ct(array{px, py, pz}, RecoDecay::distance(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}), m); }); DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, //! [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) -> float { return RecoDecay::impParXY(array{xVtxP, yVtxP, zVtxP}, array{xVtxS, yVtxS, zVtxS}, array{px, py, pz}); }); + +constexpr static int useDCAFitterN = 0; +constexpr static int useKFParticle = 1; } // namespace hf_cand // specific 2-prong decay properties @@ -371,6 +371,10 @@ DECLARE_SOA_COLUMN(FlagMcMatchRec, flagMcMatchRec, int8_t); //! reconstruction l DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); //! generator level DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); //! particle origin, reconstruction level DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); //! particle origin, generator level +// KF related properties +DECLARE_SOA_COLUMN(KfTopolChi2OverNdf, kfTopolChi2OverNdf, float); //! chi2overndf of the KFParticle topological constraint +DECLARE_SOA_COLUMN(KfGeoMassD0, kfGeoMassD0, float); //! mass of the D0 candidate from the KFParticle geometric fit +DECLARE_SOA_COLUMN(KfGeoMassD0bar, kfGeoMassD0bar, float); //! mass of the D0bar candidate from the KFParticle geometric fit // mapping of decay types enum DecayType { D0ToPiK = 0, @@ -523,8 +527,8 @@ DECLARE_SOA_EXTENDED_TABLE_USER(HfCand2ProngExt, HfCand2ProngBase, "HFCAND2PEXT" using HfCand2Prong = HfCand2ProngExt; DECLARE_SOA_TABLE(HfCand2ProngKF, "AOD", "HFCAND2PKF", - hf_cand::KFTopChi2OverNDF, - hf_cand::KFGeoMass_DZero, hf_cand::KFGeoMass_DZeroBar); + hf_cand_2prong::KfTopolChi2OverNdf, + hf_cand_2prong::KfGeoMassD0, hf_cand_2prong::KfGeoMassD0bar); // table with results of reconstruction level MC matching DECLARE_SOA_TABLE(HfCand2ProngMcRec, "AOD", "HFCAND2PMCREC", //! diff --git a/PWGHF/TableProducer/CMakeLists.txt b/PWGHF/TableProducer/CMakeLists.txt index 86c0db0fe19..23f99038990 100644 --- a/PWGHF/TableProducer/CMakeLists.txt +++ b/PWGHF/TableProducer/CMakeLists.txt @@ -25,7 +25,7 @@ o2physics_add_dpl_workflow(refit-pv-dummy o2physics_add_dpl_workflow(candidate-creator-2prong SOURCES candidateCreator2Prong.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DCAFitter KFParticle::KFParticle COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(candidate-creator-3prong diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 8a87d946d2f..81db0b603f6 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -14,6 +14,7 @@ /// /// \author Gian Michele Innocenti , CERN /// \author Vít Kučera , CERN +/// \author Pengzhong Lu , GSI Darmstadt, USTC #include "DCAFitter/DCAFitterN.h" #include "Framework/AnalysisTask.h" @@ -27,11 +28,11 @@ /// includes KFParticle #include "Tools/KFparticle/KFUtilities.h" -#include "KFParticle.h" -#include "KFPTrack.h" -#include "KFPVertex.h" -#include "KFParticleBase.h" -#include "KFVertex.h" +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; @@ -44,7 +45,7 @@ struct HfCandidateCreator2Prong { Produces rowCandidateKF; // vertexing - Configurable useKFTop2PV{"useKFTop2PV", true, "constraint KFParticle to PV, only if useKFParticle is true"}; + Configurable constrainKfToPv{"constrainKfToPv", true, "constraint KFParticle to PV"}; Configurable propagateToPCA{"propagateToPCA", true, "create tracks version propagated to PCA"}; Configurable useAbsDCA{"useAbsDCA", false, "Minimise abs. distance rather than chi2"}; Configurable useWeightedFinalPCA{"useWeightedFinalPCA", false, "Recalculate vertex position using track covariances, effective only if useAbsDCA is true"}; @@ -74,8 +75,8 @@ struct HfCandidateCreator2Prong { double bz = 0.; OutputObj hMass2{TH1F("hMass2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; - OutputObj hKFMassFailed{TH2F("hKFMassFailed", ";inv. mass with topo fit failed;entries", 500, 0., 5., 3, -0.5, 2.5)}; // 0 for D0 failed, 1 for D0bar failed - OutputObj hMassTop2{TH1F("hMassTop2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 600, -1.0, 5.)}; + OutputObj hKFMassFailed{TH1F("hKFMassFailed", ";inv. mass with topol fit failed;entries", 500, 0., 5.)}; // 0 for D0 failed, 1 for D0bar failed + OutputObj hMass2Topol{TH1F("hMass2Topol", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 600, -1.0, 5.)}; OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hCovPVYY{TH1F("hCovPVYY", "2-prong candidates;YY element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; @@ -90,12 +91,12 @@ struct HfCandidateCreator2Prong { void init(InitContext const&) { - if ((doprocessPvRefit && doprocessNoPvRefit) || (doprocessPvRefitWithKF && doprocessNoPvRefitWithKF)) { - LOGP(fatal, "Only one process function between processPvRefit and processNoPvRefit can be enabled at a time."); - } - if ((doprocessPvRefit && (doprocessPvRefitWithKF || doprocessNoPvRefitWithKF)) || (doprocessNoPvRefit && (doprocessPvRefitWithKF || doprocessNoPvRefitWithKF))) { - LOGP(fatal, "Only one process function between DCAFitterN and KFParticle can be enabled at a time."); + std::array doprocessDF{doprocessPvRefit, doprocessNoPvRefit}; + std::array doprocessKF{doprocessPvRefitWithKF, doprocessNoPvRefitWithKF}; + if ((std::accumulate(doprocessDF.begin(), doprocessDF.end(), 0) + std::accumulate(doprocessKF.begin(), doprocessKF.end(), 0)) != 1) { + LOGP(fatal, "Only one process function can be enabled at a time."); } + ccdb->setURL(ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -199,7 +200,7 @@ struct HfCandidateCreator2Prong { // get uncertainty of the decay length double phi, theta; - getPointDirection(array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); + getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); auto errorDecayLength = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); auto errorDecayLengthXY = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); @@ -255,7 +256,7 @@ struct HfCandidateCreator2Prong { // 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_t covMatrixPV[6]; + float covMatrixPV[6]; KFParticle::SetField(bz); KFPVertex kfpVertex = createKFPVertexFromCollision(collision); @@ -275,93 +276,88 @@ struct HfCandidateCreator2Prong { hCovPVXZ->Fill(covMatrixPV[3]); hCovPVZZ->Fill(covMatrixPV[5]); - KFPTrack kfpTrackPosPi; - KFPTrack kfpTrackNegPi; - KFPTrack kfpTrackPosKa; - KFPTrack kfpTrackNegKa; + KFPTrack kfpTrack0; + KFPTrack kfpTrack1; - if (track0.sign() == 1 && track1.sign() == -1) { - kfpTrackPosPi = createKFPTrackFromTrack(track0); - kfpTrackNegKa = createKFPTrackFromTrack(track1); - kfpTrackPosKa = createKFPTrackFromTrack(track0); - kfpTrackNegPi = createKFPTrackFromTrack(track1); - } + kfpTrack0 = createKFPTrackFromTrack(track0); + kfpTrack1 = createKFPTrackFromTrack(track1); - KFParticle KFPosPion(kfpTrackPosPi, 211); - KFParticle KFNegPion(kfpTrackNegPi, 211); - KFParticle KFPosKaon(kfpTrackPosKa, 321); - KFParticle KFNegKaon(kfpTrackNegKa, 321); + KFParticle kfPosPion(kfpTrack0, kPiPlus); + KFParticle kfNegPion(kfpTrack1, kPiPlus); + KFParticle kfPosKaon(kfpTrack0, kKPlus); + KFParticle kfNegKaon(kfpTrack1, kKPlus); float impactParameter0XY = 0., errImpactParameter0XY = 0., impactParameter1XY = 0., errImpactParameter1XY = 0.; - if (!KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + if (!kfPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { hDcaXYProngs->Fill(track0.pt(), impactParameter0XY * toMicrometers); - hDcaZProngs->Fill(track0.pt(), sqrt(KFPosPion.GetDistanceFromVertex(KFPV) * KFPosPion.GetDistanceFromVertex(KFPV) - impactParameter0XY * impactParameter0XY) * toMicrometers); - } - if (KFPosPion.GetDistanceFromVertexXY(KFPV, impactParameter0XY, errImpactParameter0XY)) { + hDcaZProngs->Fill(track0.pt(), std::sqrt(kfPosPion.GetDistanceFromVertex(KFPV) * kfPosPion.GetDistanceFromVertex(KFPV) - impactParameter0XY * impactParameter0XY) * toMicrometers); + } else { hDcaXYProngs->Fill(track0.pt(), -999.); hDcaZProngs->Fill(track0.pt(), -999.); } - if (!KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + if (!kfNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { hDcaXYProngs->Fill(track1.pt(), impactParameter1XY * toMicrometers); - hDcaZProngs->Fill(track1.pt(), sqrt(KFNegPion.GetDistanceFromVertex(KFPV) * KFNegPion.GetDistanceFromVertex(KFPV) - impactParameter1XY * impactParameter1XY) * toMicrometers); - } - if (KFNegPion.GetDistanceFromVertexXY(KFPV, impactParameter1XY, errImpactParameter1XY)) { + hDcaZProngs->Fill(track1.pt(), std::sqrt(kfNegPion.GetDistanceFromVertex(KFPV) * kfNegPion.GetDistanceFromVertex(KFPV) - impactParameter1XY * impactParameter1XY) * toMicrometers); + } else { hDcaXYProngs->Fill(track1.pt(), -999.); hDcaZProngs->Fill(track1.pt(), -999.); } - KFParticle KFDZero; - const KFParticle* D0Daughters[2] = {&KFPosPion, &KFNegKaon}; - KFDZero.SetConstructMethod(2); - KFDZero.Construct(D0Daughters, 2); - KFParticle KFDZeroBar; - const KFParticle* D0BarDaughters[2] = {&KFNegPion, &KFPosKaon}; - KFDZeroBar.SetConstructMethod(2); - KFDZeroBar.Construct(D0BarDaughters, 2); - - hCovSVXX->Fill(KFDZero.Covariance(0, 0)); - hCovSVYY->Fill(KFDZero.Covariance(1, 1)); - hCovSVXZ->Fill(KFDZero.Covariance(2, 0)); - hCovSVZZ->Fill(KFDZero.Covariance(2, 2)); - auto covMatrixSV = KFDZero.CovarianceMatrix(); + KFParticle kfCandD0; + const KFParticle* kfDaughtersD0[2] = {&kfPosPion, &kfNegKaon}; + kfCandD0.SetConstructMethod(2); + kfCandD0.Construct(kfDaughtersD0, 2); + KFParticle kfCandD0bar; + const KFParticle* kfDaughtersD0bar[2] = {&kfNegPion, &kfPosKaon}; + kfCandD0bar.SetConstructMethod(2); + kfCandD0bar.Construct(kfDaughtersD0bar, 2); + + auto massD0 = kfCandD0.GetMass(); + auto massD0bar = kfCandD0bar.GetMass(); + + hCovSVXX->Fill(kfCandD0.Covariance(0, 0)); + hCovSVYY->Fill(kfCandD0.Covariance(1, 1)); + hCovSVXZ->Fill(kfCandD0.Covariance(2, 0)); + hCovSVZZ->Fill(kfCandD0.Covariance(2, 2)); + auto covMatrixSV = kfCandD0.CovarianceMatrix(); double phi, theta; - getPointDirection(array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, array{KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ()}, phi, theta); + getPointDirection(std::array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, array{kfCandD0.GetX(), kfCandD0.GetY(), kfCandD0.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.)); - Float_t topChi2PerNDF_DZero = -999.; - KFParticle KFDZeroTop2PV; - if (useKFTop2PV) { - KFDZeroTop2PV = KFDZero; - KFDZeroTop2PV.SetProductionVertex(KFPV); - topChi2PerNDF_DZero = KFDZeroTop2PV.GetChi2() / KFDZeroTop2PV.GetNDF(); + float topolChi2PerNdfD0 = -999.; + KFParticle kfCandD0Topol2PV; + if (constrainKfToPv) { + kfCandD0Topol2PV = kfCandD0; + kfCandD0Topol2PV.SetProductionVertex(KFPV); + topolChi2PerNdfD0 = kfCandD0Topol2PV.GetChi2() / kfCandD0Topol2PV.GetNDF(); } // fill candidate table rows rowCandidateBase(collision.globalIndex(), KFPV.GetX(), KFPV.GetY(), KFPV.GetZ(), - KFDZero.GetX(), KFDZero.GetY(), KFDZero.GetZ(), - errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one - KFDZero.GetChi2() / KFDZero.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one - KFPosPion.GetPx(), KFPosPion.GetPy(), KFPosPion.GetPz(), - KFNegKaon.GetPx(), KFNegKaon.GetPy(), KFNegKaon.GetPz(), + kfCandD0.GetX(), kfCandD0.GetY(), kfCandD0.GetZ(), + errorDecayLength, errorDecayLengthXY, // TODO: much different from the DCAFitterN one + kfCandD0.GetChi2() / kfCandD0.GetNDF(), // TODO: to make sure it should be chi2 only or chi2/ndf, much different from the DCAFitterN one + kfPosPion.GetPx(), kfPosPion.GetPy(), kfPosPion.GetPz(), + kfNegKaon.GetPx(), kfNegKaon.GetPy(), kfNegKaon.GetPz(), impactParameter0XY, impactParameter1XY, errImpactParameter0XY, errImpactParameter1XY, rowTrackIndexProng2.prong0Id(), rowTrackIndexProng2.prong1Id(), rowTrackIndexProng2.hfflag()); - rowCandidateKF(topChi2PerNDF_DZero, - KFDZero.GetMass(), KFDZeroBar.GetMass()); + rowCandidateKF(topolChi2PerNdfD0, + massD0, massD0bar); // fill histograms if (fillHistograms) { - hMass2->Fill(KFDZero.GetMass()); - hMass2->Fill(KFDZeroBar.GetMass()); - if (useKFTop2PV) { - Float_t massDZeroTop = 0., errMassDZeroTop = 0.; - hMassTop2->Fill(KFDZeroTop2PV.GetMass()); - if (KFDZeroTop2PV.GetMass(massDZeroTop, errMassDZeroTop)) - hKFMassFailed->Fill(KFDZero.GetMass(), 0); // Set topological constraint to DZero failed + hMass2->Fill(massD0); + hMass2->Fill(massD0bar); + if (constrainKfToPv) { + float massD0Topol = 0., errMassD0Topol = 0.; + if (kfCandD0Topol2PV.GetMass(massD0Topol, errMassD0Topol)) + hKFMassFailed->Fill(massD0); // Set up D0 topological constraints to PV failed + hMass2Topol->Fill(massD0Topol); } } } @@ -374,7 +370,6 @@ struct HfCandidateCreator2Prong { { runCreator2Prong(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefit, "Run candidate creator with PV refit", false); void processNoPvRefit(aod::Collisions const& collisions, @@ -384,7 +379,6 @@ struct HfCandidateCreator2Prong { { runCreator2Prong(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefit, "Run candidate creator without PV refit", true); void processPvRefitWithKF(aod::Collisions const& collisions, @@ -394,7 +388,6 @@ struct HfCandidateCreator2Prong { { runCreator2ProngWithKF(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefitWithKF, "Run candidate creator with PV refit", false); void processNoPvRefitWithKF(aod::Collisions const& collisions, @@ -404,7 +397,6 @@ struct HfCandidateCreator2Prong { { runCreator2ProngWithKF(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithKF, "Run candidate creator without PV refit", false); }; diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index ec4c558bb87..a4f350469b5 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -28,9 +28,6 @@ using namespace o2::framework; using namespace o2::aod::hf_cand_2prong; using namespace o2::analysis::hf_cuts_d0_to_pi_k; -constexpr static int useDCAFitterN = 0; -constexpr static int useKFParticle = 1; - /// Struct for applying D0 selection cuts struct HfCandidateSelectorD0 { Produces hfSelD0Candidate; @@ -59,6 +56,7 @@ struct HfCandidateSelectorD0 { TrackSelectorPi selectorPion; TrackSelectorKa selectorKaon; + using cand2ProngKF = soa::Join; using TracksSel = soa::Join; void init(InitContext& initContext) @@ -109,8 +107,7 @@ struct HfCandidateSelectorD0 { // candidate chi2 // if constexpr (ReconstructionType == useKFParticle) { - - // if (candidate.kfTopChi2OverNDF() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; + // if (candidate.kfTopolChi2OverNdf() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; // return false; // } @@ -152,14 +149,13 @@ struct HfCandidateSelectorD0 { // invariant-mass cut float massD0, massD0bar; - if constexpr (ReconstructionType == useDCAFitterN) { + if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + massD0 = candidate.kfGeoMassD0(); + massD0bar = candidate.kfGeoMassD0bar(); + } else { massD0 = invMassD0ToPiK(candidate); massD0bar = invMassD0barToKPi(candidate); } - if constexpr (ReconstructionType == useKFParticle) { - massD0 = candidate.kfGeoMass_DZero(); - massD0bar = candidate.kfGeoMass_DZeroBar(); - } if (trackPion.sign() > 0) { if (std::abs(massD0 - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { return false; @@ -240,9 +236,9 @@ struct HfCandidateSelectorD0 { // need to add special cuts (additional cuts on decay length and d0 norm) // conjugate-dependent topological selection for D0 - bool topolD0 = selectionTopolConjugate(candidate, trackPos, trackNeg); + bool topolD0 = selectionTopolConjugate(candidate, trackPos, trackNeg); // conjugate-dependent topological selection for D0bar - bool topolD0bar = selectionTopolConjugate(candidate, trackNeg, trackPos); + bool topolD0bar = selectionTopolConjugate(candidate, trackNeg, trackPos); if (!topolD0 && !topolD0bar) { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); @@ -308,12 +304,12 @@ struct HfCandidateSelectorD0 { { processSel<0>(candidates, tracks); } - void processWithKFParticle(soa::Join const& candidates, TracksSel const& tracks) + PROCESS_SWITCH(HfCandidateSelectorD0, processWithDCAFitterN, "process candidates selection with DCAFitterN", true); + + void processWithKFParticle(cand2ProngKF const& candidates, TracksSel const& tracks) { processSel<1>(candidates, tracks); } - - PROCESS_SWITCH(HfCandidateSelectorD0, processWithDCAFitterN, "process candidates selection with DCAFitterN", false); PROCESS_SWITCH(HfCandidateSelectorD0, processWithKFParticle, "process candidates selection with KFParticle", false); }; diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index 9831aae66bf..541af3b535b 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -28,9 +28,6 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_2prong; -constexpr static int useDCAFitterN = 0; -constexpr static int useKFParticle = 1; - namespace o2::aod { namespace full @@ -122,7 +119,7 @@ DECLARE_SOA_TABLE(HfCandD0Fulls, "AOD", "HFCANDD0FULL", hf_cand::ErrorDecayLength, hf_cand::ErrorDecayLengthXY, hf_cand::Chi2PCA, - hf_cand::KFTopChi2OverNDF, + hf_cand_2prong::KfTopolChi2OverNdf, full::RSecondaryVertex, full::DecayLength, full::DecayLengthXY, @@ -200,14 +197,12 @@ struct HfTreeCreatorD0ToKPi { Configurable fillCandidateLiteTable{"fillCandidateLiteTable", false, "Switch to fill lite table with candidate properties"}; // parameters for production of training samples - Configurable fillOnlySignal{"fillOnlySignal", false, "Flag to fill derived tables with signal for ML trainings"}; - Configurable fillOnlyBackground{"fillOnlyBackground", false, "Flag to fill derived tables with background for ML trainings"}; Configurable downSampleBkgFactor{"downSampleBkgFactor", 1., "Fraction of background candidates to keep for ML trainings"}; Configurable ptMaxForDownSample{"ptMaxForDownSample", 10., "Maximum pt for the application of the downsampling factor"}; using TracksWPid = soa::Join; using SelectedCandidatesMc = soa::Filtered>; - using SelectedCandidatesMcKF = soa::Filtered>; + using SelectedCandidatesMcKf = soa::Filtered>; using MatchedGenCandidatesMc = soa::Filtered>; Filter filterSelectCandidates = aod::hf_sel_candidate_d0::isSelD0 >= 1 || aod::hf_sel_candidate_d0::isSelD0bar >= 1; @@ -215,11 +210,16 @@ struct HfTreeCreatorD0ToKPi { Partition reconstructedCandSig = nabs(aod::hf_cand_2prong::flagMcMatchRec) == static_cast(BIT(DecayType::D0ToPiK)); Partition reconstructedCandBkg = nabs(aod::hf_cand_2prong::flagMcMatchRec) != static_cast(BIT(DecayType::D0ToPiK)); - Partition reconstructedCandSigKF = nabs(aod::hf_cand_2prong::flagMcMatchRec) == static_cast(BIT(DecayType::D0ToPiK)); - Partition reconstructedCandBkgKF = nabs(aod::hf_cand_2prong::flagMcMatchRec) != static_cast(BIT(DecayType::D0ToPiK)); + Partition reconstructedCandSigKF = nabs(aod::hf_cand_2prong::flagMcMatchRec) == static_cast(BIT(DecayType::D0ToPiK)); + Partition reconstructedCandBkgKF = nabs(aod::hf_cand_2prong::flagMcMatchRec) != static_cast(BIT(DecayType::D0ToPiK)); void init(InitContext const&) { + std::array doprocessData{doprocessDataWithDCAFitterN, doprocessDataWithKFParticle}; + std::array doprocessMc{doprocessMcWithDCAFitterOnlySig, doprocessMcWithDCAFitterOnlyBkg, doprocessMcWithDCAFitterAll, doprocessMcWithKFParticleOnlySig, doprocessMcWithKFParticleOnlyBkg, doprocessMcWithKFParticleAll}; + if ((std::accumulate(doprocessData.begin(), doprocessData.end(), 0) + std::accumulate(doprocessMc.begin(), doprocessMc.end(), 0)) != 1) { + LOGP(fatal, "Only one process function can be enabled at a time."); + } } template @@ -363,20 +363,21 @@ struct HfTreeCreatorD0ToKPi { double yD = yD0(candidate); double eD = eD0(candidate); double ctD = ctD0(candidate); - if constexpr (ReconstructionType == useDCAFitterN) { - if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), -999., ctD, yD, eD, 0, 0); - } - if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), -999., ctD, yD, eD, 0, 0); - } - } else if constexpr (ReconstructionType == useKFParticle) { - if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, candidate.kfGeoMass_DZero(), cosThetaStarD0(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, 0, 0); - } - if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, candidate.kfGeoMass_DZeroBar(), cosThetaStarD0bar(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, 0, 0); - } + float massD0, massD0bar; + float topolChi2PerNdf = -999.; + if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + massD0 = candidate.kfGeoMassD0(); + massD0bar = candidate.kfGeoMassD0bar(); + topolChi2PerNdf = candidate.kfTopolChi2OverNdf(); + } else { + massD0 = invMassD0ToPiK(candidate); + massD0bar = invMassD0barToKPi(candidate); + } + if (candidate.isSelD0()) { + fillTable(candidate, prong0, prong1, 0, massD0, cosThetaStarD0(candidate), topolChi2PerNdf, ctD, yD, eD, 0, 0); + } + if (candidate.isSelD0bar()) { + fillTable(candidate, prong0, prong1, 1, massD0bar, cosThetaStarD0bar(candidate), topolChi2PerNdf, ctD, yD, eD, 0, 0); } } } @@ -387,20 +388,19 @@ struct HfTreeCreatorD0ToKPi { { processData<0>(collisions, candidates, tracks, BCs); } + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithDCAFitterN, "Process data with DCAFitterN", true); + void processDataWithKFParticle(aod::Collisions const& collisions, soa::Filtered> const& candidates, TracksWPid const& tracks, aod::BCs const& BCs) { processData<1>(collisions, candidates, tracks, BCs); } - - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithDCAFitterN, "Process data with DCAFitterN", true); PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithKFParticle, "Process data with KFParticle", false); - template + template void processMc(aod::Collisions const& collisions, aod::McCollisions const&, - THfCand2Prong const& candidates, THfCand2ProngSel const& reconstructedCandSets, MatchedGenCandidatesMc const& particles, TracksWPid const&, aod::BCs const&) @@ -429,20 +429,21 @@ struct HfTreeCreatorD0ToKPi { double yD = yD0(candidate); double eD = eD0(candidate); double ctD = ctD0(candidate); - if constexpr (ReconstructionType == useDCAFitterN) { - if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, invMassD0ToPiK(candidate), cosThetaStarD0(candidate), -999., ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); - } - if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, invMassD0barToKPi(candidate), cosThetaStarD0bar(candidate), -999., ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); - } - } else if constexpr (ReconstructionType == useKFParticle) { - if (candidate.isSelD0()) { - fillTable(candidate, prong0, prong1, 0, candidate.kfGeoMass_DZero(), cosThetaStarD0(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); - } - if (candidate.isSelD0bar()) { - fillTable(candidate, prong0, prong1, 1, candidate.kfGeoMass_DZeroBar(), cosThetaStarD0bar(candidate), candidate.kfTopChi2OverNDF(), ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); - } + float massD0, massD0bar; + float topolChi2PerNdf = -999.; + if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + massD0 = candidate.kfGeoMassD0(); + massD0bar = candidate.kfGeoMassD0bar(); + topolChi2PerNdf = candidate.kfTopolChi2OverNdf(); + } else { + massD0 = invMassD0ToPiK(candidate); + massD0bar = invMassD0barToKPi(candidate); + } + if (candidate.isSelD0()) { + fillTable(candidate, prong0, prong1, 0, massD0, cosThetaStarD0(candidate), topolChi2PerNdf, ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); + } + if (candidate.isSelD0bar()) { + fillTable(candidate, prong0, prong1, 1, massD0bar, cosThetaStarD0bar(candidate), topolChi2PerNdf, ctD, yD, eD, candidate.flagMcMatchRec(), candidate.originMcRec()); } } @@ -465,58 +466,58 @@ struct HfTreeCreatorD0ToKPi { void processMcWithDCAFitterOnlySig(aod::Collisions const& collisions, aod::McCollisions const& McCollisions, - SelectedCandidatesMc const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<0>(collisions, McCollisions, candidates, reconstructedCandSig, particles, tracks, BCs); + processMc<0>(collisions, McCollisions, reconstructedCandSig, particles, tracks, BCs); } + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlySig, "Process MC with DCAFitterN only for signals", false); + void processMcWithDCAFitterOnlyBkg(aod::Collisions const& collisions, aod::McCollisions const& McCollisions, - SelectedCandidatesMc const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<0, true>(collisions, McCollisions, candidates, reconstructedCandBkg, particles, tracks, BCs); + processMc<0, true>(collisions, McCollisions, reconstructedCandBkg, particles, tracks, BCs); } + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlyBkg, "Process MC with DCAFitterN only for background", false); + void processMcWithDCAFitterAll(aod::Collisions const& collisions, aod::McCollisions const& McCollisions, SelectedCandidatesMc const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<0>(collisions, McCollisions, candidates, candidates, particles, tracks, BCs); + processMc<0>(collisions, McCollisions, candidates, particles, tracks, BCs); } + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterAll, "Process MC with DCAFitterN", false); + void processMcWithKFParticleOnlySig(aod::Collisions const& collisions, aod::McCollisions const& McCollisions, - SelectedCandidatesMcKF const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<1>(collisions, McCollisions, candidates, reconstructedCandSigKF, particles, tracks, BCs); + processMc<1>(collisions, McCollisions, reconstructedCandSigKF, particles, tracks, BCs); } + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlySig, "Process MC with KFParticle only for signals", false); + void processMcWithKFParticleOnlyBkg(aod::Collisions const& collisions, aod::McCollisions const& McCollisions, - SelectedCandidatesMcKF const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<1, true>(collisions, McCollisions, candidates, reconstructedCandBkgKF, particles, tracks, BCs); + processMc<1, true>(collisions, McCollisions, reconstructedCandBkgKF, particles, tracks, BCs); } + PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlyBkg, "Process MC with KFParticle only for background", false); + void processMcWithKFParticleAll(aod::Collisions const& collisions, aod::McCollisions const& McCollisions, - SelectedCandidatesMcKF const& candidates, + SelectedCandidatesMcKf const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<1>(collisions, McCollisions, candidates, candidates, particles, tracks, BCs); + processMc<1>(collisions, McCollisions, candidates, particles, tracks, BCs); } - - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlySig, "Process MC with DCAFitterN only for signals", false); - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlyBkg, "Process MC with DCAFitterN only for background", false); - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterAll, "Process MC with DCAFitterN", false); - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlySig, "Process MC with KFParticle only for signals", false); - PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlyBkg, "Process MC with KFParticle only for background", false); PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleAll, "Process MC with KFParticle", false); }; From 79f7ef21eed9036a2d6f208e57157bee7962a0be Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 16 Aug 2023 10:00:30 +0200 Subject: [PATCH 12/24] PWGHF: candidateSelectorD0.cxx conflict fix --- PWGHF/TableProducer/candidateSelectorD0.cxx | 43 +++++++++++++++++---- 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index 43a6e37578c..5061b9b52de 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -76,6 +76,7 @@ struct HfCandidateSelectorD0 { TrackSelectorPi selectorPion; TrackSelectorKa selectorKaon; + using cand2ProngKF = soa::Join; using TracksSel = soa::Join; void init(InitContext& initContext) @@ -102,7 +103,7 @@ struct HfCandidateSelectorD0 { /// Conjugate-independent topological cuts /// \param candidate is candidate /// \return true if candidate passes all cuts - template + template bool selectionTopol(const T& candidate) { auto candpT = candidate.pt(); @@ -134,6 +135,12 @@ struct HfCandidateSelectorD0 { // candidate DCA // if (candidate.chi2PCA() > cuts[pTBin][1]) return false; + // candidate chi2 + // if constexpr (ReconstructionType == useKFParticle) { + // if (candidate.kfTopolChi2OverNdf() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; + // return false; + // } + // decay exponentail law, with tau = beta*gamma*ctau // decay length > ctau retains (1-1/e) if (std::abs(candidate.impactParameterNormalised0()) < 0.5 || std::abs(candidate.impactParameterNormalised1()) < 0.5) { @@ -161,7 +168,7 @@ struct HfCandidateSelectorD0 { /// \param trackKaon is the track with the kaon hypothesis /// \note trackPion = positive and trackKaon = negative for D0 selection and inverse for D0bar /// \return true if candidate passes all cuts for the given Conjugate - template + template bool selectionTopolConjugate(const T1& candidate, const T2& trackPion, const T2& trackKaon) { auto candpT = candidate.pt(); @@ -171,12 +178,20 @@ struct HfCandidateSelectorD0 { } // invariant-mass cut + float massD0, massD0bar; + if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + massD0 = candidate.kfGeoMassD0(); + massD0bar = candidate.kfGeoMassD0bar(); + } else { + massD0 = invMassD0ToPiK(candidate); + massD0bar = invMassD0barToKPi(candidate); + } if (trackPion.sign() > 0) { - if (std::abs(invMassD0ToPiK(candidate) - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { + if (std::abs(massD0 - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { return false; } } else { - if (std::abs(invMassD0barToKPi(candidate) - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { + if (std::abs(massD0bar - RecoDecay::getMassPDG(pdg::Code::kD0)) > cuts->get(pTBin, "m")) { return false; } } @@ -217,8 +232,8 @@ struct HfCandidateSelectorD0 { return true; } - - void process(aod::HfCand2Prong const& candidates, TracksSel const&) + template + void processSel(THfCand2Prong const& candidates, TracksSel const&) { // looping over 2-prong candidates for (auto& candidate : candidates) { @@ -241,8 +256,8 @@ struct HfCandidateSelectorD0 { statusHFFlag = 1; auto ptCand = candidate.pt(); - auto trackPos = candidate.prong0_as(); // positive daughter - auto trackNeg = candidate.prong1_as(); // negative daughter + auto trackPos = candidate.template prong0_as(); // positive daughter + auto trackNeg = candidate.template prong1_as(); // negative daughter // conjugate-independent topological selection if (!selectionTopol(candidate)) { @@ -349,6 +364,18 @@ struct HfCandidateSelectorD0 { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); } } + + void processWithDCAFitterN(aod::HfCand2Prong const& candidates, TracksSel const& tracks) + { + processSel<0>(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorD0, processWithDCAFitterN, "process candidates selection with DCAFitterN", true); + + void processWithKFParticle(cand2ProngKF const& candidates, TracksSel const& tracks) + { + processSel<1>(candidates, tracks); + } + PROCESS_SWITCH(HfCandidateSelectorD0, processWithKFParticle, "process candidates selection with KFParticle", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From a4c00a0a1e8c8d4349bb9dec80b8f73c6f3d0447 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 13 Sep 2023 15:29:22 +0200 Subject: [PATCH 13/24] Merge commit '76acd5f6ff94b2fa82fd234526cf84f09daf68a9', fix conflicts in PWGHF/TableProducer/candidateCreatorD0.cxx --- PWGHF/TableProducer/candidateSelectorD0.cxx | 5 ----- 1 file changed, 5 deletions(-) diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index 9d76358c34c..805d049054e 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -246,14 +246,9 @@ struct HfCandidateSelectorD0 { return true; } -<<<<<<< HEAD - template - void processSel(THfCand2Prong const& candidates, TracksSel const&) -======= void process(aod::HfCand2Prong const& candidates, TracksSel const&) ->>>>>>> 76acd5f6ff94b2fa82fd234526cf84f09daf68a9 { // looping over 2-prong candidates for (const auto& candidate : candidates) { From 11be4ea02ecfaa74c9a683203d1a8dceffc71a7e Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 13 Sep 2023 16:23:45 +0200 Subject: [PATCH 14/24] Updates for adding KF to HF-D0 --- PWGHF/D2H/Tasks/taskD0.cxx | 29 ++++-- .../DataModel/CandidateReconstructionTables.h | 8 +- .../TableProducer/candidateCreator2Prong.cxx | 96 +++++++++---------- PWGHF/TableProducer/candidateSelectorD0.cxx | 24 ++--- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 38 ++++---- 5 files changed, 101 insertions(+), 94 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index 679278859df..1c998901201 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -26,6 +26,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_2prong; +using namespace o2::aod::hf_cand; using namespace o2::analysis::hf_cuts_d0_to_pi_k; /// D0 analysis task @@ -163,9 +164,10 @@ struct HfTaskD0 { registry.add("hDecLengthxyVsPtSig", "2-prong candidates;decay length xy (cm) vs #it{p}_{T} for signal;entries", {HistType::kTH2F, {{800, 0., 4.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } - void process(soa::Join const& candidates) + template + void processData(THfCand2ProngSel& selectedD0CandidatesSets) { - for (const auto& candidate : selectedD0Candidates) { + for (auto& candidate : selectedD0CandidatesSets) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -174,7 +176,7 @@ struct HfTaskD0 { } float massD0, massD0bar; - if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + if constexpr (ReconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); } else { @@ -225,17 +227,24 @@ struct HfTaskD0 { } void processWithDCAFitterN(soa::Join const&) { - processData<0>(selectedD0Candidates); + processData(selectedD0Candidates); } PROCESS_SWITCH(HfTaskD0, processWithDCAFitterN, "process taskD0 with DCAFitterN", true); - void processMc(soa::Join const& candidates, + void processWithKFParticle(soa::Join const&) + { + processData(selectedD0CandidatesKF); + } + PROCESS_SWITCH(HfTaskD0, processWithKFParticle, "process taskD0 with KFParticle", false); + + template + void processMc(THfCand2Prong const& candidates, soa::Join const& mcParticles, aod::TracksWMc const& tracks) { // MC rec. // Printf("MC Candidates: %d", candidates.size()); - for (const auto& candidate : recoFlag2Prong) { + for (auto& candidate : candidates) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -243,7 +252,7 @@ struct HfTaskD0 { continue; } float massD0, massD0bar; - if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + if constexpr (ReconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); } else { @@ -252,7 +261,7 @@ struct HfTaskD0 { } if (std::abs(candidate.flagMcMatchRec()) == 1 << DecayType::D0ToPiK) { // Get the corresponding MC particle. - auto indexMother = RecoDecay::getMother(mcParticles, candidate.prong0_as().mcParticle_as>(), pdg::Code::kD0, true); + auto indexMother = RecoDecay::getMother(mcParticles, candidate.template prong0_as().template mcParticle_as>(), pdg::Code::kD0, true); auto particleMother = mcParticles.rawIteratorAt(indexMother); auto ptGen = particleMother.pt(); // gen. level pT auto yGen = RecoDecay::y(std::array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::getMassPDG(particleMother.pdgCode())); // gen. level y @@ -425,13 +434,13 @@ struct HfTaskD0 { void processMcWithDCAFitterN(soa::Join const& particlesMC, aod::TracksWMc const& tracks) { - processMc<0>(recoFlag2Prong, particlesMC, tracks); + processMc(recoFlag2Prong, particlesMC, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); void processMcWithKFParticle(soa::Join const& particlesMC, aod::TracksWMc const& tracks) { - processMc<1>(recoFlag2ProngKF, particlesMC, tracks); + processMc(recoFlag2ProngKF, particlesMC, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithKFParticle, "Process MC with KFParticle", false); }; diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 9b6ffc450a5..81182f5ea9a 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -353,6 +353,11 @@ DECLARE_SOA_DYNAMIC_COLUMN(Ct, ct, //! [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz, double m) -> float { return RecoDecay::ct(std::array{px, py, pz}, RecoDecay::distance(std::array{xVtxP, yVtxP, zVtxP}, std::array{xVtxS, yVtxS, zVtxS}), m); }); DECLARE_SOA_DYNAMIC_COLUMN(ImpactParameterXY, impactParameterXY, //! [](float xVtxP, float yVtxP, float zVtxP, float xVtxS, float yVtxS, float zVtxS, float px, float py, float pz) -> float { return RecoDecay::impParXY(std::array{xVtxP, yVtxP, zVtxP}, std::array{xVtxS, yVtxS, zVtxS}, std::array{px, py, pz}); }); +DECLARE_SOA_COLUMN(KfTopolChi2OverNdf, kfTopolChi2OverNdf, float); //! chi2overndf of the KFParticle topological constraint + +// method of secondary-vertex reconstruction +enum VertexerType { DCAFitter = 0, + KfParticle }; } // namespace hf_cand // specific 2-prong decay properties @@ -382,7 +387,6 @@ DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); //! generator level DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); //! particle origin, reconstruction level DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); //! particle origin, generator level // KF related properties -DECLARE_SOA_COLUMN(KfTopolChi2OverNdf, kfTopolChi2OverNdf, float); //! chi2overndf of the KFParticle topological constraint DECLARE_SOA_COLUMN(KfGeoMassD0, kfGeoMassD0, float); //! mass of the D0 candidate from the KFParticle geometric fit DECLARE_SOA_COLUMN(KfGeoMassD0bar, kfGeoMassD0bar, float); //! mass of the D0bar candidate from the KFParticle geometric fit @@ -537,7 +541,7 @@ DECLARE_SOA_EXTENDED_TABLE_USER(HfCand2ProngExt, HfCand2ProngBase, "HFCAND2PEXT" using HfCand2Prong = HfCand2ProngExt; DECLARE_SOA_TABLE(HfCand2ProngKF, "AOD", "HFCAND2PKF", - hf_cand_2prong::KfTopolChi2OverNdf, + hf_cand::KfTopolChi2OverNdf, hf_cand_2prong::KfGeoMassD0, hf_cand_2prong::KfGeoMassD0bar); // table with results of reconstruction level MC matching diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 3136946ca77..002e2590750 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -75,8 +75,6 @@ struct HfCandidateCreator2Prong { double bz = 0.; OutputObj hMass2{TH1F("hMass2", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 500, 0., 5.)}; - OutputObj hKFMassFailed{TH1F("hKFMassFailed", ";inv. mass with topol fit failed;entries", 500, 0., 5.)}; // 0 for D0 failed, 1 for D0bar failed - OutputObj hMass2Topol{TH1F("hMass2Topol", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", 600, -1.0, 5.)}; OutputObj hCovPVXX{TH1F("hCovPVXX", "2-prong candidates;XX element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; OutputObj hCovSVXX{TH1F("hCovSVXX", "2-prong candidates;XX element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hCovPVYY{TH1F("hCovPVYY", "2-prong candidates;YY element of cov. matrix of prim. vtx. position (cm^{2});entries", 100, 0., 1.e-4)}; @@ -91,11 +89,17 @@ struct HfCandidateCreator2Prong { void init(InitContext const&) { - std::array doprocessDF{doprocessPvRefit, doprocessNoPvRefit}; - std::array doprocessKF{doprocessPvRefitWithKF, doprocessNoPvRefitWithKF}; + std::array doprocessDF{doprocessPvRefitWithDCAFitterN, doprocessNoPvRefitWithDCAFitterN}; + std::array doprocessKF{doprocessPvRefitWithKFParticle, doprocessNoPvRefitWithKFParticle}; if ((std::accumulate(doprocessDF.begin(), doprocessDF.end(), 0) + std::accumulate(doprocessKF.begin(), doprocessKF.end(), 0)) != 1) { LOGP(fatal, "Only one process function can be enabled at a time."); } + if (std::accumulate(doprocessDF.begin(), doprocessDF.end(), 0) == 1) { + hUseKForDCAFitter->Fill(1.0); + } + if (std::accumulate(doprocessKF.begin(), doprocessKF.end(), 0) == 1) { + hUseKForDCAFitter->Fill(0.0); + } ccdb->setURL(ccdbUrl); ccdb->setCaching(true); @@ -105,10 +109,10 @@ struct HfCandidateCreator2Prong { } template - void runCreator2Prong(aod::Collisions const& collisions, - Cand const& rowsTrackIndexProng2, - TTracks const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + void runCreator2ProngWithDCAFitterN(aod::Collisions const& collisions, + Cand const& rowsTrackIndexProng2, + TTracks const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { // 2-prong vertex fitter o2::vertexing::DCAFitterN<2> df; @@ -120,7 +124,6 @@ struct HfCandidateCreator2Prong { df.setMinRelChi2Change(minRelChi2Change); df.setUseAbsDCA(useAbsDCA); df.setWeightedFinalPCA(useWeightedFinalPCA); - hUseKForDCAFitter->Fill(1.0); // loop over pairs of track indices for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { @@ -142,7 +145,6 @@ struct HfCandidateCreator2Prong { // df.setBz(bz); /// put it outside the 'if'! Otherwise we have a difference wrt bz Configurable (< 1 permille) in Run2 conv. data // df.print(); } - df.setBz(bz); // reconstruct the 2-prong secondary vertex @@ -230,12 +232,11 @@ struct HfCandidateCreator2Prong { } template - void runCreator2ProngWithKF(aod::Collisions const& collisions, - Cand const& rowsTrackIndexProng2, - TTracks const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + void runCreator2ProngWithKFParticle(aod::Collisions const& collisions, + Cand const& rowsTrackIndexProng2, + TTracks const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - hUseKForDCAFitter->Fill(0.0); for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { auto track0 = rowTrackIndexProng2.template prong0_as(); @@ -276,11 +277,8 @@ struct HfCandidateCreator2Prong { hCovPVXZ->Fill(covMatrixPV[3]); hCovPVZZ->Fill(covMatrixPV[5]); - KFPTrack kfpTrack0; - KFPTrack kfpTrack1; - - kfpTrack0 = createKFPTrackFromTrack(track0); - kfpTrack1 = createKFPTrackFromTrack(track1); + KFPTrack kfpTrack0 = createKFPTrackFromTrack(track0); + KFPTrack kfpTrack1 = createKFPTrackFromTrack(track1); KFParticle kfPosPion(kfpTrack0, kPiPlus); KFParticle kfNegPion(kfpTrack1, kPiPlus); @@ -322,7 +320,7 @@ struct HfCandidateCreator2Prong { auto covMatrixSV = kfCandD0.CovarianceMatrix(); double phi, theta; - getPointDirection(std::array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, array{kfCandD0.GetX(), kfCandD0.GetY(), kfCandD0.GetZ()}, phi, theta); + getPointDirection(std::array{KFPV.GetX(), KFPV.GetY(), KFPV.GetZ()}, std::array{kfCandD0.GetX(), kfCandD0.GetY(), kfCandD0.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.)); @@ -353,51 +351,47 @@ struct HfCandidateCreator2Prong { if (fillHistograms) { hMass2->Fill(massD0); hMass2->Fill(massD0bar); - if (constrainKfToPv) { - float massD0Topol = 0., errMassD0Topol = 0.; - if (kfCandD0Topol2PV.GetMass(massD0Topol, errMassD0Topol)) - hKFMassFailed->Fill(massD0); // Set up D0 topological constraints to PV failed - hMass2Topol->Fill(massD0Topol); - } } } } - void processPvRefit(aod::Collisions const& collisions, - soa::Join const& rowsTrackIndexProng2, - aod::TracksWCov const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + void processPvRefitWithDCAFitterN(aod::Collisions const& collisions, + soa::Join const& rowsTrackIndexProng2, + aod::TracksWCov const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator2Prong(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + runCreator2ProngWithDCAFitterN(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefit, "Run candidate creator with PV refit", false); - void processNoPvRefit(aod::Collisions const& collisions, - aod::Hf2Prongs const& rowsTrackIndexProng2, - aod::TracksWCov const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefitWithDCAFitterN, "Run candidate creator with PV refit", false); + + void processNoPvRefitWithDCAFitterN(aod::Collisions const& collisions, + aod::Hf2Prongs const& rowsTrackIndexProng2, + aod::TracksWCov const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator2Prong(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + runCreator2ProngWithDCAFitterN(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefit, "Run candidate creator without PV refit", true); - void processPvRefitWithKF(aod::Collisions const& collisions, - soa::Join const& rowsTrackIndexProng2, - soa::Join const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithDCAFitterN, "Run candidate creator without PV refit", true); + + void processPvRefitWithKFParticle(aod::Collisions const& collisions, + soa::Join const& rowsTrackIndexProng2, + soa::Join const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator2ProngWithKF(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + runCreator2ProngWithKFParticle(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefitWithKF, "Run candidate creator with PV refit", false); + PROCESS_SWITCH(HfCandidateCreator2Prong, processPvRefitWithKFParticle, "Run candidate creator with PV refit", false); - void processNoPvRefitWithKF(aod::Collisions const& collisions, - aod::Hf2Prongs const& rowsTrackIndexProng2, - soa::Join const& tracks, - aod::BCsWithTimestamps const& bcWithTimeStamps) + void processNoPvRefitWithKFParticle(aod::Collisions const& collisions, + aod::Hf2Prongs const& rowsTrackIndexProng2, + soa::Join const& tracks, + aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator2ProngWithKF(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + runCreator2ProngWithKFParticle(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } - PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithKF, "Run candidate creator without PV refit", false); + PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithKFParticle, "Run candidate creator without PV refit", false); }; /// Extends the base table with expression columns. diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index 805d049054e..15173b45b7f 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -27,6 +27,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::aod::hf_cand_2prong; +using namespace o2::aod::hf_cand; using namespace o2::analysis::hf_cuts_d0_to_pi_k; /// Struct for applying D0 selection cuts @@ -149,10 +150,9 @@ struct HfCandidateSelectorD0 { // candidate DCA // if (candidate.chi2PCA() > cuts[pTBin][1]) return false; - // candidate chi2 - // if constexpr (ReconstructionType == useKFParticle) { + // candidate topological chi2 over ndf when using KFParticle, need to add this selection to the SelectorCuts.h + // if constexpr (ReconstructionType == VertexerType::KfParticle) { // if (candidate.kfTopolChi2OverNdf() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; - // return false; // } // decay exponentail law, with tau = beta*gamma*ctau @@ -193,7 +193,7 @@ struct HfCandidateSelectorD0 { // invariant-mass cut float massD0, massD0bar; - if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + if constexpr (ReconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); } else { @@ -246,9 +246,9 @@ struct HfCandidateSelectorD0 { return true; } - - void process(aod::HfCand2Prong const& candidates, - TracksSel const&) + template + void processSel(THfCand2Prong const& candidates, + TracksSel const&) { // looping over 2-prong candidates for (const auto& candidate : candidates) { @@ -275,7 +275,7 @@ struct HfCandidateSelectorD0 { auto trackNeg = candidate.template prong1_as(); // negative daughter // conjugate-independent topological selection - if (!selectionTopol(candidate)) { + if (!selectionTopol(candidate)) { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); if (applyMl) { hfMlD0Candidate(outputMl); @@ -288,9 +288,9 @@ struct HfCandidateSelectorD0 { // need to add special cuts (additional cuts on decay length and d0 norm) // conjugate-dependent topological selection for D0 - bool topolD0 = selectionTopolConjugate(candidate, trackPos, trackNeg); + bool topolD0 = selectionTopolConjugate(candidate, trackPos, trackNeg); // conjugate-dependent topological selection for D0bar - bool topolD0bar = selectionTopolConjugate(candidate, trackNeg, trackPos); + bool topolD0bar = selectionTopolConjugate(candidate, trackNeg, trackPos); if (!topolD0 && !topolD0bar) { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); @@ -390,13 +390,13 @@ struct HfCandidateSelectorD0 { void processWithDCAFitterN(aod::HfCand2Prong const& candidates, TracksSel const& tracks) { - processSel<0>(candidates, tracks); + processSel(candidates, tracks); } PROCESS_SWITCH(HfCandidateSelectorD0, processWithDCAFitterN, "process candidates selection with DCAFitterN", true); void processWithKFParticle(cand2ProngKF const& candidates, TracksSel const& tracks) { - processSel<1>(candidates, tracks); + processSel(candidates, tracks); } PROCESS_SWITCH(HfCandidateSelectorD0, processWithKFParticle, "process candidates selection with KFParticle", false); }; diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index 79df70a887d..e16e6471c44 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -27,6 +27,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::aod::hf_cand_2prong; +using namespace o2::aod::hf_cand; namespace o2::aod { @@ -119,7 +120,7 @@ DECLARE_SOA_TABLE(HfCandD0Fulls, "AOD", "HFCANDD0FULL", hf_cand::ErrorDecayLength, hf_cand::ErrorDecayLengthXY, hf_cand::Chi2PCA, - hf_cand_2prong::KfTopolChi2OverNdf, + hf_cand::KfTopolChi2OverNdf, full::RSecondaryVertex, full::DecayLength, full::DecayLengthXY, @@ -215,9 +216,8 @@ struct HfTreeCreatorD0ToKPi { void init(InitContext const&) { - std::array doprocessData{doprocessDataWithDCAFitterN, doprocessDataWithKFParticle}; - std::array doprocessMc{doprocessMcWithDCAFitterOnlySig, doprocessMcWithDCAFitterOnlyBkg, doprocessMcWithDCAFitterAll, doprocessMcWithKFParticleOnlySig, doprocessMcWithKFParticleOnlyBkg, doprocessMcWithKFParticleAll}; - if ((std::accumulate(doprocessData.begin(), doprocessData.end(), 0) + std::accumulate(doprocessMc.begin(), doprocessMc.end(), 0)) != 1) { + std::array doprocess{doprocessDataWithDCAFitterN, doprocessDataWithKFParticle, doprocessMcWithDCAFitterOnlySig, doprocessMcWithDCAFitterOnlyBkg, doprocessMcWithDCAFitterAll, doprocessMcWithKFParticleOnlySig, doprocessMcWithKFParticleOnlyBkg, doprocessMcWithKFParticleAll}; + if (std::accumulate(doprocess.begin(), doprocess.end(), 0) != 1) { LOGP(fatal, "Only one process function can be enabled at a time."); } } @@ -365,7 +365,7 @@ struct HfTreeCreatorD0ToKPi { double ctD = ctD0(candidate); float massD0, massD0bar; float topolChi2PerNdf = -999.; - if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + if constexpr (ReconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); topolChi2PerNdf = candidate.kfTopolChi2OverNdf(); @@ -386,7 +386,7 @@ struct HfTreeCreatorD0ToKPi { soa::Filtered> const& candidates, TracksWPid const& tracks, aod::BCs const& BCs) { - processData<0>(collisions, candidates, tracks, BCs); + processData(collisions, candidates, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithDCAFitterN, "Process data with DCAFitterN", true); @@ -394,14 +394,14 @@ struct HfTreeCreatorD0ToKPi { soa::Filtered> const& candidates, TracksWPid const& tracks, aod::BCs const& BCs) { - processData<1>(collisions, candidates, tracks, BCs); + processData(collisions, candidates, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithKFParticle, "Process data with KFParticle", false); - template + template void processMc(aod::Collisions const& collisions, aod::McCollisions const&, - THfCand2ProngSel const& reconstructedCandSets, + THfCand2Prong const& candidates, MatchedGenCandidatesMc const& particles, TracksWPid const&, aod::BCs const&) { @@ -413,11 +413,11 @@ struct HfTreeCreatorD0ToKPi { // Filling candidate properties if (fillCandidateLiteTable) { - rowCandidateLite.reserve(reconstructedCandSets.size()); + rowCandidateLite.reserve(candidates.size()); } else { - rowCandidateFull.reserve(reconstructedCandSets.size()); + rowCandidateFull.reserve(candidates.size()); } - for (const auto& candidate : reconstructedCandSets) { + for (const auto& candidate : candidates) { if (onlyBkg && downSampleBkgFactor < 1.) { float pseudoRndm = candidate.ptProng0() * 1000. - (int64_t)(candidate.ptProng0() * 1000); if (candidate.pt() < ptMaxForDownSample && pseudoRndm >= downSampleBkgFactor) { @@ -431,7 +431,7 @@ struct HfTreeCreatorD0ToKPi { double ctD = ctD0(candidate); float massD0, massD0bar; float topolChi2PerNdf = -999.; - if constexpr (ReconstructionType == o2::aod::hf_cand::useKFParticle) { + if constexpr (ReconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); topolChi2PerNdf = candidate.kfTopolChi2OverNdf(); @@ -469,7 +469,7 @@ struct HfTreeCreatorD0ToKPi { MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<0>(collisions, McCollisions, reconstructedCandSig, particles, tracks, BCs); + processMc(collisions, McCollisions, reconstructedCandSig, particles, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlySig, "Process MC with DCAFitterN only for signals", false); @@ -478,7 +478,7 @@ struct HfTreeCreatorD0ToKPi { MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<0, true>(collisions, McCollisions, reconstructedCandBkg, particles, tracks, BCs); + processMc(collisions, McCollisions, reconstructedCandBkg, particles, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlyBkg, "Process MC with DCAFitterN only for background", false); @@ -488,7 +488,7 @@ struct HfTreeCreatorD0ToKPi { MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<0>(collisions, McCollisions, candidates, particles, tracks, BCs); + processMc(collisions, McCollisions, candidates, particles, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterAll, "Process MC with DCAFitterN", false); @@ -497,7 +497,7 @@ struct HfTreeCreatorD0ToKPi { MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<1>(collisions, McCollisions, reconstructedCandSigKF, particles, tracks, BCs); + processMc(collisions, McCollisions, reconstructedCandSigKF, particles, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlySig, "Process MC with KFParticle only for signals", false); @@ -506,7 +506,7 @@ struct HfTreeCreatorD0ToKPi { MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<1, true>(collisions, McCollisions, reconstructedCandBkgKF, particles, tracks, BCs); + processMc(collisions, McCollisions, reconstructedCandBkgKF, particles, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlyBkg, "Process MC with KFParticle only for background", false); @@ -516,7 +516,7 @@ struct HfTreeCreatorD0ToKPi { MatchedGenCandidatesMc const& particles, TracksWPid const& tracks, aod::BCs const& BCs) { - processMc<1>(collisions, McCollisions, candidates, particles, tracks, BCs); + processMc(collisions, McCollisions, candidates, particles, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleAll, "Process MC with KFParticle", false); }; From 57101d3a1a6ce7f1fbd9445e78a7a2dab21517af Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 13 Sep 2023 16:56:34 +0200 Subject: [PATCH 15/24] git-clang-format for "PWGHF/DataModel/CandidateReconstructionTables.h" --- PWGHF/DataModel/CandidateReconstructionTables.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGHF/DataModel/CandidateReconstructionTables.h b/PWGHF/DataModel/CandidateReconstructionTables.h index 81182f5ea9a..653dcd57240 100644 --- a/PWGHF/DataModel/CandidateReconstructionTables.h +++ b/PWGHF/DataModel/CandidateReconstructionTables.h @@ -387,8 +387,8 @@ DECLARE_SOA_COLUMN(FlagMcMatchGen, flagMcMatchGen, int8_t); //! generator level DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int8_t); //! particle origin, reconstruction level DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int8_t); //! particle origin, generator level // KF related properties -DECLARE_SOA_COLUMN(KfGeoMassD0, kfGeoMassD0, float); //! mass of the D0 candidate from the KFParticle geometric fit -DECLARE_SOA_COLUMN(KfGeoMassD0bar, kfGeoMassD0bar, float); //! mass of the D0bar candidate from the KFParticle geometric fit +DECLARE_SOA_COLUMN(KfGeoMassD0, kfGeoMassD0, float); //! mass of the D0 candidate from the KFParticle geometric fit +DECLARE_SOA_COLUMN(KfGeoMassD0bar, kfGeoMassD0bar, float); //! mass of the D0bar candidate from the KFParticle geometric fit // mapping of decay types enum DecayType { D0ToPiK = 0, From f2ba845891496f0cf89281d1209de63f1679e287 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 13 Sep 2023 23:46:08 +0200 Subject: [PATCH 16/24] PWGHF: Try to respect the grouping and sorting of the head files --- .../TableProducer/candidateCreator2Prong.cxx | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index 002e2590750..02b681b0d5c 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -16,24 +16,27 @@ /// \author Vít Kučera , CERN /// \author Pengzhong Lu , GSI Darmstadt, USTC +#ifndef HomogeneousField +#define HomogeneousField +#endif + +#include +#include +#include +#include +#include + #include "DCAFitter/DCAFitterN.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/DCA.h" #include "Common/Core/trackUtilities.h" +#include "Tools/KFparticle/KFUtilities.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/Utils/utilsBfieldCCDB.h" -/// includes KFParticle -#include "Tools/KFparticle/KFUtilities.h" -#include -#include -#include -#include -#include - using namespace o2; using namespace o2::framework; using namespace o2::aod::hf_cand; From 004d84ec45c07dce5ed5d31ce79bdda543b24d3d Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Fri, 15 Sep 2023 12:57:44 +0200 Subject: [PATCH 17/24] Update PR --- PWGHF/D2H/Tasks/taskD0.cxx | 30 +++++++++---------- .../TableProducer/candidateCreator2Prong.cxx | 14 ++++----- PWGHF/TableProducer/candidateSelectorD0.cxx | 20 +++++++------ PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 14 +++++---- 4 files changed, 41 insertions(+), 37 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index b3aa7d41c65..b5976a337ef 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -43,8 +43,8 @@ struct HfTaskD0 { Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; Partition> selectedD0CandidatesKF = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> recoFlag2Prong = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; - Partition> recoFlag2ProngKF = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; + Partition> selectedD0CandidatesMc = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; + Partition> selectedD0CandidatesMcKF = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; HistogramRegistry registry{ "registry", @@ -164,10 +164,10 @@ struct HfTaskD0 { registry.add("hDecLengthxyVsPtSig", "2-prong candidates;decay length xy (cm) vs #it{p}_{T} for signal;entries", {HistType::kTH2F, {{800, 0., 4.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); } - template - void processData(THfCand2ProngSel& selectedD0CandidatesSets) + template + void processData(CandType& candidates) { - for (auto& candidate : selectedD0CandidatesSets) { + for (const auto& candidate : candidates) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { continue; } @@ -176,7 +176,7 @@ struct HfTaskD0 { } float massD0, massD0bar; - if constexpr (ReconstructionType == VertexerType::KfParticle) { + if constexpr (reconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); } else { @@ -225,20 +225,20 @@ struct HfTaskD0 { registry.fill(HIST("hCPAXYFinerBinning"), candidate.cpaXY(), ptCandidate); } } - void processWithDCAFitterN(soa::Join const&) + void processDataWithDCAFitterN(soa::Join const&) { processData(selectedD0Candidates); } - PROCESS_SWITCH(HfTaskD0, processWithDCAFitterN, "process taskD0 with DCAFitterN", true); + PROCESS_SWITCH(HfTaskD0, processDataWithDCAFitterN, "process taskD0 with DCAFitterN", true); - void processWithKFParticle(soa::Join const&) + void processDataWithKFParticle(soa::Join const&) { processData(selectedD0CandidatesKF); } - PROCESS_SWITCH(HfTaskD0, processWithKFParticle, "process taskD0 with KFParticle", false); + PROCESS_SWITCH(HfTaskD0, processDataWithKFParticle, "process taskD0 with KFParticle", false); - template - void processMc(THfCand2Prong const& candidates, + template + void processMc(CandType const& candidates, soa::Join const& mcParticles, aod::TracksWMc const& tracks) { @@ -251,7 +251,7 @@ struct HfTaskD0 { continue; } float massD0, massD0bar; - if constexpr (ReconstructionType == VertexerType::KfParticle) { + if constexpr (reconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); } else { @@ -432,13 +432,13 @@ struct HfTaskD0 { void processMcWithDCAFitterN(soa::Join const& particlesMC, aod::TracksWMc const& tracks) { - processMc(recoFlag2Prong, particlesMC, tracks); + processMc(selectedD0CandidatesMc, particlesMC, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); void processMcWithKFParticle(soa::Join const& particlesMC, aod::TracksWMc const& tracks) { - processMc(recoFlag2ProngKF, particlesMC, tracks); + processMc(selectedD0CandidatesMcKF, particlesMC, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithKFParticle, "Process MC with KFParticle", false); }; diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index f35b82b846d..c10d3382107 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -88,7 +88,7 @@ struct HfCandidateCreator2Prong { OutputObj hCovSVZZ{TH1F("hCovSVZZ", "2-prong candidates;ZZ element of cov. matrix of sec. vtx. position (cm^{2});entries", 100, 0., 0.2)}; OutputObj hDcaXYProngs{TH2F("hDcaXYProngs", "DCAxy of 2-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{xy}) (#mum);entries", 100, 0., 20., 200, -500., 500.)}; OutputObj hDcaZProngs{TH2F("hDcaZProngs", "DCAz of 2-prong candidates;#it{p}_{T} (GeV/#it{c};#it{d}_{z}) (#mum);entries", 100, 0., 20., 200, -500., 500.)}; - OutputObj hUseKForDCAFitter{TH1F("hUseKForDCAFitter", ";Use KF or DCAFitterN;entries", 2, -0.5, 1.5)}; // 0 for KF, 1 for DCAFitterN + OutputObj hVertexerType{TH1F("hVertexerType", "Use KF or DCAFitterN;Vertexer type;entries", 2, -0.5, 1.5)}; // See o2::aod::hf_cand::VertexerType void init(InitContext const&) { @@ -98,10 +98,10 @@ struct HfCandidateCreator2Prong { LOGP(fatal, "Only one process function can be enabled at a time."); } if (std::accumulate(doprocessDF.begin(), doprocessDF.end(), 0) == 1) { - hUseKForDCAFitter->Fill(1.0); + hVertexerType->Fill(VertexerType::DCAFitter); } if (std::accumulate(doprocessKF.begin(), doprocessKF.end(), 0) == 1) { - hUseKForDCAFitter->Fill(0.0); + hVertexerType->Fill(VertexerType::KfParticle); } ccdb->setURL(ccdbUrl); @@ -111,7 +111,7 @@ struct HfCandidateCreator2Prong { runNumber = 0; } - template + template void runCreator2ProngWithDCAFitterN(aod::Collisions const& collisions, Cand const& rowsTrackIndexProng2, TTracks const& tracks, @@ -234,7 +234,7 @@ struct HfCandidateCreator2Prong { } } - template + template void runCreator2ProngWithKFParticle(aod::Collisions const& collisions, Cand const& rowsTrackIndexProng2, TTracks const& tracks, @@ -373,7 +373,7 @@ struct HfCandidateCreator2Prong { aod::TracksWCov const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator2ProngWithDCAFitterN(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + runCreator2ProngWithDCAFitterN(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithDCAFitterN, "Run candidate creator without PV refit", true); @@ -392,7 +392,7 @@ struct HfCandidateCreator2Prong { soa::Join const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { - runCreator2ProngWithKFParticle(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); + runCreator2ProngWithKFParticle(collisions, rowsTrackIndexProng2, tracks, bcWithTimeStamps); } PROCESS_SWITCH(HfCandidateCreator2Prong, processNoPvRefitWithKFParticle, "Run candidate creator without PV refit", false); }; diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index 15173b45b7f..9e0cb43a523 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -116,9 +116,10 @@ struct HfCandidateSelectorD0 { } /// Conjugate-independent topological cuts + /// \param reconstructionType is the reconstruction type (DCAFitterN or KFParticle) /// \param candidate is candidate /// \return true if candidate passes all cuts - template + template bool selectionTopol(const T& candidate) { auto candpT = candidate.pt(); @@ -151,7 +152,7 @@ struct HfCandidateSelectorD0 { // if (candidate.chi2PCA() > cuts[pTBin][1]) return false; // candidate topological chi2 over ndf when using KFParticle, need to add this selection to the SelectorCuts.h - // if constexpr (ReconstructionType == VertexerType::KfParticle) { + // if constexpr (reconstructionType == VertexerType::KfParticle) { // if (candidate.kfTopolChi2OverNdf() > cuts->get(pTBin, "topological chi2overndf as D0")) return false; // } @@ -177,12 +178,13 @@ struct HfCandidateSelectorD0 { } /// Conjugate-dependent topological cuts + /// \param reconstructionType is the reconstruction type (DCAFitterN or KFParticle) /// \param candidate is candidate /// \param trackPion is the track with the pion hypothesis /// \param trackKaon is the track with the kaon hypothesis /// \note trackPion = positive and trackKaon = negative for D0 selection and inverse for D0bar /// \return true if candidate passes all cuts for the given Conjugate - template + template bool selectionTopolConjugate(const T1& candidate, const T2& trackPion, const T2& trackKaon) { auto candpT = candidate.pt(); @@ -193,7 +195,7 @@ struct HfCandidateSelectorD0 { // invariant-mass cut float massD0, massD0bar; - if constexpr (ReconstructionType == VertexerType::KfParticle) { + if constexpr (reconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); } else { @@ -246,7 +248,7 @@ struct HfCandidateSelectorD0 { return true; } - template + template void processSel(THfCand2Prong const& candidates, TracksSel const&) { @@ -275,7 +277,7 @@ struct HfCandidateSelectorD0 { auto trackNeg = candidate.template prong1_as(); // negative daughter // conjugate-independent topological selection - if (!selectionTopol(candidate)) { + if (!selectionTopol(candidate)) { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); if (applyMl) { hfMlD0Candidate(outputMl); @@ -288,9 +290,9 @@ struct HfCandidateSelectorD0 { // need to add special cuts (additional cuts on decay length and d0 norm) // conjugate-dependent topological selection for D0 - bool topolD0 = selectionTopolConjugate(candidate, trackPos, trackNeg); + bool topolD0 = selectionTopolConjugate(candidate, trackPos, trackNeg); // conjugate-dependent topological selection for D0bar - bool topolD0bar = selectionTopolConjugate(candidate, trackNeg, trackPos); + bool topolD0bar = selectionTopolConjugate(candidate, trackNeg, trackPos); if (!topolD0 && !topolD0bar) { hfSelD0Candidate(statusD0, statusD0bar, statusHFFlag, statusTopol, statusCand, statusPID); @@ -394,7 +396,7 @@ struct HfCandidateSelectorD0 { } PROCESS_SWITCH(HfCandidateSelectorD0, processWithDCAFitterN, "process candidates selection with DCAFitterN", true); - void processWithKFParticle(cand2ProngKF const& candidates, TracksSel const& tracks) + void processWithKFParticle(soa::Join const& candidates, TracksSel const& tracks) { processSel(candidates, tracks); } diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index e16e6471c44..51a02e36a02 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -334,7 +334,7 @@ struct HfTreeCreatorD0ToKPi { } } - template + template void processData(aod::Collisions const& collisions, THfCand2Prong const& candidates, TracksWPid const&, aod::BCs const&) @@ -365,7 +365,7 @@ struct HfTreeCreatorD0ToKPi { double ctD = ctD0(candidate); float massD0, massD0bar; float topolChi2PerNdf = -999.; - if constexpr (ReconstructionType == VertexerType::KfParticle) { + if constexpr (reconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); topolChi2PerNdf = candidate.kfTopolChi2OverNdf(); @@ -384,7 +384,8 @@ struct HfTreeCreatorD0ToKPi { void processDataWithDCAFitterN(aod::Collisions const& collisions, soa::Filtered> const& candidates, - TracksWPid const& tracks, aod::BCs const& BCs) + TracksWPid const& tracks, + aod::BCs const& BCs) { processData(collisions, candidates, tracks, BCs); } @@ -392,13 +393,14 @@ struct HfTreeCreatorD0ToKPi { void processDataWithKFParticle(aod::Collisions const& collisions, soa::Filtered> const& candidates, - TracksWPid const& tracks, aod::BCs const& BCs) + TracksWPid const& tracks, + aod::BCs const& BCs) { processData(collisions, candidates, tracks, BCs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithKFParticle, "Process data with KFParticle", false); - template + template void processMc(aod::Collisions const& collisions, aod::McCollisions const&, THfCand2Prong const& candidates, @@ -431,7 +433,7 @@ struct HfTreeCreatorD0ToKPi { double ctD = ctD0(candidate); float massD0, massD0bar; float topolChi2PerNdf = -999.; - if constexpr (ReconstructionType == VertexerType::KfParticle) { + if constexpr (reconstructionType == VertexerType::KfParticle) { massD0 = candidate.kfGeoMassD0(); massD0bar = candidate.kfGeoMassD0bar(); topolChi2PerNdf = candidate.kfTopolChi2OverNdf(); From 502bbf99cb7eeada2f32974358096e3ef9dbca89 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Mon, 18 Sep 2023 00:14:44 +0200 Subject: [PATCH 18/24] Update PR --- .../TableProducer/candidateCreator2Prong.cxx | 8 +- PWGHF/TableProducer/candidateSelectorD0.cxx | 4 +- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 79 ++++++++++--------- 3 files changed, 49 insertions(+), 42 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index c10d3382107..d576a483a1f 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -111,9 +111,9 @@ struct HfCandidateCreator2Prong { runNumber = 0; } - template + template void runCreator2ProngWithDCAFitterN(aod::Collisions const& collisions, - Cand const& rowsTrackIndexProng2, + CandType const& rowsTrackIndexProng2, TTracks const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { @@ -234,9 +234,9 @@ struct HfCandidateCreator2Prong { } } - template + template void runCreator2ProngWithKFParticle(aod::Collisions const& collisions, - Cand const& rowsTrackIndexProng2, + CandType const& rowsTrackIndexProng2, TTracks const& tracks, aod::BCsWithTimestamps const& bcWithTimeStamps) { diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index 9e0cb43a523..f9abd12a136 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -248,8 +248,8 @@ struct HfCandidateSelectorD0 { return true; } - template - void processSel(THfCand2Prong const& candidates, + template + void processSel(CandType const& candidates, TracksSel const&) { // looping over 2-prong candidates diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index 51a02e36a02..3c78e0f69ca 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -334,9 +334,9 @@ struct HfTreeCreatorD0ToKPi { } } - template + template void processData(aod::Collisions const& collisions, - THfCand2Prong const& candidates, + CandType const& candidates, TracksWPid const&, aod::BCs const&) { // Filling event properties @@ -385,27 +385,28 @@ struct HfTreeCreatorD0ToKPi { void processDataWithDCAFitterN(aod::Collisions const& collisions, soa::Filtered> const& candidates, TracksWPid const& tracks, - aod::BCs const& BCs) + aod::BCs const& bcs) { - processData(collisions, candidates, tracks, BCs); + processData(collisions, candidates, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithDCAFitterN, "Process data with DCAFitterN", true); void processDataWithKFParticle(aod::Collisions const& collisions, soa::Filtered> const& candidates, TracksWPid const& tracks, - aod::BCs const& BCs) + aod::BCs const& bcs) { - processData(collisions, candidates, tracks, BCs); + processData(collisions, candidates, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processDataWithKFParticle, "Process data with KFParticle", false); - template + template void processMc(aod::Collisions const& collisions, aod::McCollisions const&, - THfCand2Prong const& candidates, - MatchedGenCandidatesMc const& particles, - TracksWPid const&, aod::BCs const&) + CandType const& candidates, + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const&, + aod::BCs const&) { // Filling event properties rowCandidateFullEvents.reserve(collisions.size()); @@ -450,8 +451,8 @@ struct HfTreeCreatorD0ToKPi { } // Filling particle properties - rowCandidateFullParticles.reserve(particles.size()); - for (const auto& particle : particles) { + rowCandidateFullParticles.reserve(mcParticles.size()); + for (const auto& particle : mcParticles) { if (std::abs(particle.flagMcMatchGen()) == 1 << DecayType::D0ToPiK) { rowCandidateFullParticles( particle.mcCollisionId(), @@ -467,58 +468,64 @@ struct HfTreeCreatorD0ToKPi { } void processMcWithDCAFitterOnlySig(aod::Collisions const& collisions, - aod::McCollisions const& McCollisions, - MatchedGenCandidatesMc const& particles, - TracksWPid const& tracks, aod::BCs const& BCs) + aod::McCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const& tracks, + aod::BCs const& bcs) { - processMc(collisions, McCollisions, reconstructedCandSig, particles, tracks, BCs); + processMc(collisions, mcCollisions, reconstructedCandSig, mcParticles, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlySig, "Process MC with DCAFitterN only for signals", false); void processMcWithDCAFitterOnlyBkg(aod::Collisions const& collisions, - aod::McCollisions const& McCollisions, - MatchedGenCandidatesMc const& particles, - TracksWPid const& tracks, aod::BCs const& BCs) + aod::McCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const& tracks, + aod::BCs const& bcs) { - processMc(collisions, McCollisions, reconstructedCandBkg, particles, tracks, BCs); + processMc(collisions, mcCollisions, reconstructedCandBkg, mcParticles, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterOnlyBkg, "Process MC with DCAFitterN only for background", false); void processMcWithDCAFitterAll(aod::Collisions const& collisions, - aod::McCollisions const& McCollisions, + aod::McCollisions const& mcCollisions, SelectedCandidatesMc const& candidates, - MatchedGenCandidatesMc const& particles, - TracksWPid const& tracks, aod::BCs const& BCs) + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const& tracks, + aod::BCs const& bcs) { - processMc(collisions, McCollisions, candidates, particles, tracks, BCs); + processMc(collisions, mcCollisions, candidates, mcParticles, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithDCAFitterAll, "Process MC with DCAFitterN", false); void processMcWithKFParticleOnlySig(aod::Collisions const& collisions, - aod::McCollisions const& McCollisions, - MatchedGenCandidatesMc const& particles, - TracksWPid const& tracks, aod::BCs const& BCs) + aod::McCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const& tracks, + aod::BCs const& bcs) { - processMc(collisions, McCollisions, reconstructedCandSigKF, particles, tracks, BCs); + processMc(collisions, mcCollisions, reconstructedCandSigKF, mcParticles, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlySig, "Process MC with KFParticle only for signals", false); void processMcWithKFParticleOnlyBkg(aod::Collisions const& collisions, - aod::McCollisions const& McCollisions, - MatchedGenCandidatesMc const& particles, - TracksWPid const& tracks, aod::BCs const& BCs) + aod::McCollisions const& mcCollisions, + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const& tracks, + aod::BCs const& bcs) { - processMc(collisions, McCollisions, reconstructedCandBkgKF, particles, tracks, BCs); + processMc(collisions, mcCollisions, reconstructedCandBkgKF, mcParticles, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleOnlyBkg, "Process MC with KFParticle only for background", false); void processMcWithKFParticleAll(aod::Collisions const& collisions, - aod::McCollisions const& McCollisions, + aod::McCollisions const& mcCollisions, SelectedCandidatesMcKf const& candidates, - MatchedGenCandidatesMc const& particles, - TracksWPid const& tracks, aod::BCs const& BCs) + MatchedGenCandidatesMc const& mcParticles, + TracksWPid const& tracks, + aod::BCs const& bcs) { - processMc(collisions, McCollisions, candidates, particles, tracks, BCs); + processMc(collisions, mcCollisions, candidates, mcParticles, tracks, bcs); } PROCESS_SWITCH(HfTreeCreatorD0ToKPi, processMcWithKFParticleAll, "Process MC with KFParticle", false); }; From da7f9515715fd6a9832fcb8dcbb8db41b8460811 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Sun, 24 Sep 2023 22:37:04 +0200 Subject: [PATCH 19/24] PWGHF: update PR --- PWGHF/D2H/Tasks/taskD0.cxx | 15 ++++++++++++--- PWGHF/TableProducer/candidateSelectorD0.cxx | 6 +++++- PWGHF/TableProducer/treeCreatorD0ToKPi.cxx | 4 ++++ 3 files changed, 21 insertions(+), 4 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index b5976a337ef..66472febdd8 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -124,6 +124,11 @@ struct HfTaskD0 { void init(InitContext&) { + std::array doprocess{doprocessDataWithDCAFitterN, doprocessDataWithKFParticle, doprocessMcWithDCAFitterN, doprocessMcWithKFParticle}; + if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) { + LOGP(fatal, "Only one process function can be enabled at a time."); + } + auto vbins = (std::vector)binsPt; registry.add("hMass", "2-prong candidates;inv. mass (#pi K) (GeV/#it{c}^{2});entries", {HistType::kTH2F, {{500, 0., 5.}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}}); registry.add("hMassVsPhi", "2-prong candidates vs phi;inv. mass (#pi K) (GeV/#it{c}^{2});phi (rad);entries", {HistType::kTH3F, {{120, 1.5848, 2.1848}, {vbins, "#it{p}_{T} (GeV/#it{c})"}, {32, 0, o2::constants::math::TwoPI}}}); @@ -165,7 +170,7 @@ struct HfTaskD0 { } template - void processData(CandType& candidates) + void processData(CandType const& candidates) { for (const auto& candidate : candidates) { if (!(candidate.hfflag() & 1 << DecayType::D0ToPiK)) { @@ -430,13 +435,17 @@ struct HfTaskD0 { } } - void processMcWithDCAFitterN(soa::Join const& particlesMC, aod::TracksWMc const& tracks) + void processMcWithDCAFitterN(soa::Join const& particlesMC, + aod::TracksWMc const& tracks) { processMc(selectedD0CandidatesMc, particlesMC, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); - void processMcWithKFParticle(soa::Join const& particlesMC, aod::TracksWMc const& tracks) + void processMcWithKFParticle(soa::Join const& particlesMC, + aod::TracksWMc const& tracks) { processMc(selectedD0CandidatesMcKF, particlesMC, tracks); } diff --git a/PWGHF/TableProducer/candidateSelectorD0.cxx b/PWGHF/TableProducer/candidateSelectorD0.cxx index f9abd12a136..cc3569fb86f 100644 --- a/PWGHF/TableProducer/candidateSelectorD0.cxx +++ b/PWGHF/TableProducer/candidateSelectorD0.cxx @@ -76,7 +76,6 @@ struct HfCandidateSelectorD0 { TrackSelectorPi selectorPion; TrackSelectorKa selectorKaon; - using cand2ProngKF = soa::Join; using TracksSel = soa::Join; // Define histograms @@ -87,6 +86,11 @@ struct HfCandidateSelectorD0 { void init(InitContext& initContext) { + std::array doprocess{doprocessWithDCAFitterN, doprocessWithKFParticle}; + if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) { + LOGP(fatal, "Only one process function can be enabled at a time."); + } + if (applyMl) { registry.add("DebugBdt/hBdtScore1VsStatus", ";BDT score;status", {HistType::kTH2F, {axisBdtScore, axisSelStatus}}); registry.add("DebugBdt/hBdtScore2VsStatus", ";BDT score;status", {HistType::kTH2F, {axisBdtScore, axisSelStatus}}); diff --git a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx index 3c78e0f69ca..9c90c2170a7 100644 --- a/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx +++ b/PWGHF/TableProducer/treeCreatorD0ToKPi.cxx @@ -469,6 +469,7 @@ struct HfTreeCreatorD0ToKPi { void processMcWithDCAFitterOnlySig(aod::Collisions const& collisions, aod::McCollisions const& mcCollisions, + SelectedCandidatesMc const&, MatchedGenCandidatesMc const& mcParticles, TracksWPid const& tracks, aod::BCs const& bcs) @@ -479,6 +480,7 @@ struct HfTreeCreatorD0ToKPi { void processMcWithDCAFitterOnlyBkg(aod::Collisions const& collisions, aod::McCollisions const& mcCollisions, + SelectedCandidatesMc const&, MatchedGenCandidatesMc const& mcParticles, TracksWPid const& tracks, aod::BCs const& bcs) @@ -500,6 +502,7 @@ struct HfTreeCreatorD0ToKPi { void processMcWithKFParticleOnlySig(aod::Collisions const& collisions, aod::McCollisions const& mcCollisions, + SelectedCandidatesMcKf const&, MatchedGenCandidatesMc const& mcParticles, TracksWPid const& tracks, aod::BCs const& bcs) @@ -510,6 +513,7 @@ struct HfTreeCreatorD0ToKPi { void processMcWithKFParticleOnlyBkg(aod::Collisions const& collisions, aod::McCollisions const& mcCollisions, + SelectedCandidatesMcKf const&, MatchedGenCandidatesMc const& mcParticles, TracksWPid const& tracks, aod::BCs const& bcs) From 5dc95fbc666673fd72b8c43d73046f851fd7f83c Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Sun, 24 Sep 2023 22:41:51 +0200 Subject: [PATCH 20/24] PWGHF: update git-clang-format --- PWGHF/D2H/Tasks/taskD0.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index 66472febdd8..dde163f933e 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -267,9 +267,9 @@ struct HfTaskD0 { // Get the corresponding MC particle. auto indexMother = RecoDecay::getMother(mcParticles, candidate.template prong0_as().template mcParticle_as>(), pdg::Code::kD0, true); auto particleMother = mcParticles.rawIteratorAt(indexMother); - auto ptGen = particleMother.pt(); // gen. level pT + auto ptGen = particleMother.pt(); // gen. level pT auto yGen = RecoDecay::y(std::array{particleMother.px(), particleMother.py(), particleMother.pz()}, RecoDecay::getMassPDG(particleMother.pdgCode())); // gen. level y - registry.fill(HIST("hPtGenSig"), ptGen); // gen. level pT + registry.fill(HIST("hPtGenSig"), ptGen); // gen. level pT auto ptRec = candidate.pt(); auto yRec = yD0(candidate); if (candidate.isRecoHfFlag() >= selectionFlagHf) { @@ -436,7 +436,7 @@ struct HfTaskD0 { } void processMcWithDCAFitterN(soa::Join const& particlesMC, + aod::HfCand2ProngMcGen> const& particlesMC, aod::TracksWMc const& tracks) { processMc(selectedD0CandidatesMc, particlesMC, tracks); @@ -444,7 +444,7 @@ struct HfTaskD0 { PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); void processMcWithKFParticle(soa::Join const& particlesMC, + aod::HfCand2ProngMcGen> const& particlesMC, aod::TracksWMc const& tracks) { processMc(selectedD0CandidatesMcKF, particlesMC, tracks); From b03deb6b86564bf88a373f0e23e5c530e940e371 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 27 Sep 2023 15:47:10 +0200 Subject: [PATCH 21/24] PWGHF: update PR --- PWGHF/D2H/Tasks/taskD0.cxx | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index dde163f933e..dd3b5139201 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -125,8 +125,11 @@ struct HfTaskD0 { void init(InitContext&) { std::array doprocess{doprocessDataWithDCAFitterN, doprocessDataWithKFParticle, doprocessMcWithDCAFitterN, doprocessMcWithKFParticle}; - if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) != 1) { - LOGP(fatal, "Only one process function can be enabled at a time."); + if ((std::accumulate(doprocess.begin(), doprocess.end(), 0)) == 0) { + LOGP(fatal, "At least one process function should be enabled at a time."); + } + if ((doprocessDataWithDCAFitterN || doprocessMcWithDCAFitterN) && (doprocessDataWithKFParticle || doprocessMcWithKFParticle)) { + LOGP(fatal, "DCAFitterN and KFParticle can not be enabled at a time."); } auto vbins = (std::vector)binsPt; @@ -435,19 +438,19 @@ struct HfTaskD0 { } } - void processMcWithDCAFitterN(soa::Join const& particlesMC, + void processMcWithDCAFitterN(soa::Join const&, + soa::Join const& mcParticles, aod::TracksWMc const& tracks) { - processMc(selectedD0CandidatesMc, particlesMC, tracks); + processMc(selectedD0CandidatesMc, mcParticles, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); - void processMcWithKFParticle(soa::Join const& particlesMC, + void processMcWithKFParticle(soa::Join const&, + soa::Join const& mcParticles, aod::TracksWMc const& tracks) { - processMc(selectedD0CandidatesMcKF, particlesMC, tracks); + processMc(selectedD0CandidatesMcKF, mcParticles, tracks); } PROCESS_SWITCH(HfTaskD0, processMcWithKFParticle, "Process MC with KFParticle", false); }; From 55ba45d2e9dec2a9f23e4b0e9a1dc5c9969d2f51 Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Sun, 1 Oct 2023 11:08:11 +0200 Subject: [PATCH 22/24] PWGHF: update PR --- PWGHF/D2H/Tasks/taskD0.cxx | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/PWGHF/D2H/Tasks/taskD0.cxx b/PWGHF/D2H/Tasks/taskD0.cxx index dd3b5139201..d02618ea7ce 100644 --- a/PWGHF/D2H/Tasks/taskD0.cxx +++ b/PWGHF/D2H/Tasks/taskD0.cxx @@ -41,10 +41,15 @@ struct HfTaskD0 { Configurable selectionPid{"selectionPid", 1, "Selection Flag for reco PID candidates"}; Configurable> binsPt{"binsPt", std::vector{hf_cuts_d0_to_pi_k::vecBinsPt}, "pT bin limits"}; - Partition> selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> selectedD0CandidatesKF = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; - Partition> selectedD0CandidatesMc = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; - Partition> selectedD0CandidatesMcKF = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; + using D0Candidates = soa::Join; + using D0CandidatesMc = soa::Join; + using D0CandidatesKF = soa::Join; + using D0CandidatesMcKF = soa::Join; + + Partition selectedD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + Partition selectedD0CandidatesKF = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlagD0 || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlagD0bar; + Partition selectedD0CandidatesMc = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; + Partition selectedD0CandidatesMcKF = aod::hf_sel_candidate_d0::isRecoHfFlag >= selectionFlagHf; HistogramRegistry registry{ "registry", @@ -233,13 +238,13 @@ struct HfTaskD0 { registry.fill(HIST("hCPAXYFinerBinning"), candidate.cpaXY(), ptCandidate); } } - void processDataWithDCAFitterN(soa::Join const&) + void processDataWithDCAFitterN(D0Candidates const&) { processData(selectedD0Candidates); } PROCESS_SWITCH(HfTaskD0, processDataWithDCAFitterN, "process taskD0 with DCAFitterN", true); - void processDataWithKFParticle(soa::Join const&) + void processDataWithKFParticle(D0CandidatesKF const&) { processData(selectedD0CandidatesKF); } @@ -438,7 +443,7 @@ struct HfTaskD0 { } } - void processMcWithDCAFitterN(soa::Join const&, + void processMcWithDCAFitterN(D0CandidatesMc const&, soa::Join const& mcParticles, aod::TracksWMc const& tracks) { @@ -446,7 +451,7 @@ struct HfTaskD0 { } PROCESS_SWITCH(HfTaskD0, processMcWithDCAFitterN, "Process MC with DCAFitterN", false); - void processMcWithKFParticle(soa::Join const&, + void processMcWithKFParticle(D0CandidatesMcKF const&, soa::Join const& mcParticles, aod::TracksWMc const& tracks) { From 4ea01a5fb03cf58ef57e45e58b5df698b5938bdc Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Wed, 4 Oct 2023 10:41:53 +0200 Subject: [PATCH 23/24] Remove unused variables in runCreator2ProngWithKFParticle --- PWGHF/TableProducer/candidateCreator2Prong.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/PWGHF/TableProducer/candidateCreator2Prong.cxx b/PWGHF/TableProducer/candidateCreator2Prong.cxx index d576a483a1f..418a950e013 100644 --- a/PWGHF/TableProducer/candidateCreator2Prong.cxx +++ b/PWGHF/TableProducer/candidateCreator2Prong.cxx @@ -244,8 +244,6 @@ struct HfCandidateCreator2Prong { for (const auto& rowTrackIndexProng2 : rowsTrackIndexProng2) { auto track0 = rowTrackIndexProng2.template prong0_as(); auto track1 = rowTrackIndexProng2.template prong1_as(); - auto trackParVarPos1 = getTrackParCov(track0); - auto trackParVarNeg1 = getTrackParCov(track1); auto collision = rowTrackIndexProng2.collision(); /// Set the magnetic field from ccdb. From 3b622f5e7b72b78b61c8cc9522dede309f24551c Mon Sep 17 00:00:00 2001 From: lupengzhong Date: Fri, 17 Nov 2023 00:21:52 +0800 Subject: [PATCH 24/24] PWGDQ: Add "processDecayToEESkimmedWithColl" to the init function of the AnalysisSameEventPairing task in tableReader.cxx --- PWGDQ/Tasks/tableReader.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGDQ/Tasks/tableReader.cxx b/PWGDQ/Tasks/tableReader.cxx index a769ba715f5..0989bb74ffe 100644 --- a/PWGDQ/Tasks/tableReader.cxx +++ b/PWGDQ/Tasks/tableReader.cxx @@ -809,7 +809,7 @@ struct AnalysisSameEventPairing { } } - if (context.mOptions.get("processDecayToEESkimmed") || context.mOptions.get("processDecayToEESkimmedWithCov") || context.mOptions.get("processDecayToEEVertexingSkimmed") || context.mOptions.get("processVnDecayToEESkimmed") || context.mOptions.get("processDecayToEEPrefilterSkimmed") || context.mOptions.get("processDecayToPiPiSkimmed") || context.mOptions.get("processAllSkimmed")) { + if (context.mOptions.get("processDecayToEESkimmed") || context.mOptions.get("processDecayToEESkimmedWithCov") || context.mOptions.get("processDecayToEEVertexingSkimmed") || context.mOptions.get("processVnDecayToEESkimmed") || context.mOptions.get("processDecayToEEPrefilterSkimmed") || context.mOptions.get("processDecayToEESkimmedWithColl") || context.mOptions.get("processDecayToPiPiSkimmed") || context.mOptions.get("processAllSkimmed")) { TString cutNames = fConfigTrackCuts.value; if (!cutNames.IsNull()) { // if track cuts std::unique_ptr objArray(cutNames.Tokenize(","));