diff --git a/DPG/Tasks/AOTTrack/D0CalibTables.h b/DPG/Tasks/AOTTrack/D0CalibTables.h index 81e16201b63..f868e08d3dc 100644 --- a/DPG/Tasks/AOTTrack/D0CalibTables.h +++ b/DPG/Tasks/AOTTrack/D0CalibTables.h @@ -242,6 +242,20 @@ static const std::vector labelsPtCand = { // column labels static const std::vector labelsCutVarCand = {"delta inv. mass", "max d0d0", "max pointing angle", "max pointing angle XY", "min cos pointing angle", "min cos pointing angle XY", "min norm decay length", "min norm decay length XY", "min decay length", "min decay length XY"}; +// var indices for the cuts +enum CutsIds : uint8_t { + DeltaMass = 0, + MaxD0D0 = 1, + MaxPointingAngle = 2, + MaxPointingAngleXY = 3, + MinCosPointingAngle = 4, + MinCosPointingAngleXY = 5, + MinNormDecayLength = 6, + MinNormDecayLengthXY = 7, + MinDecayLength = 8, + MinDecayLengthXY = 9 +}; + static constexpr int NBinsPtMl = 10; // default values for the pT bin edges (can be used to configure histogram axis) // offset by 1 from the bin numbers in cuts array @@ -324,10 +338,10 @@ DECLARE_SOA_COLUMN(TpcNumSigmaPi, tpcNumSigmaPi, int8_t); //! DECLARE_SOA_COLUMN(TpcNumSigmaKa, tpcNumSigmaKa, int8_t); //! compressed NsigmaTPC for kaons DECLARE_SOA_COLUMN(TofNumSigmaPi, tofNumSigmaPi, int8_t); //! compressed NsigmaTOF for pions DECLARE_SOA_COLUMN(TofNumSigmaKa, tofNumSigmaKa, int8_t); //! compressed NsigmaTOF for kaons -DECLARE_SOA_COLUMN(ITSChi2NCl, itsChi2NCl, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column -DECLARE_SOA_COLUMN(TPCChi2NCl, tpcChi2NCl, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column -DECLARE_SOA_COLUMN(TRDChi2, trdChi2, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column -DECLARE_SOA_COLUMN(TOFChi2, tofChi2, uint8_t); //! compressed NsigmaTOF for kaons // o2-linter: disable=name/o2-column +DECLARE_SOA_COLUMN(ITSChi2NCl, itsChi2NCl, uint8_t); // o2-linter: disable=name/o2-column +DECLARE_SOA_COLUMN(TPCChi2NCl, tpcChi2NCl, uint8_t); // o2-linter: disable=name/o2-column +DECLARE_SOA_COLUMN(TRDChi2, trdChi2, uint8_t); // o2-linter: disable=name/o2-column +DECLARE_SOA_COLUMN(TOFChi2, tofChi2, uint8_t); // o2-linter: disable=name/o2-column DECLARE_SOA_COLUMN(CmoPrimUnfm80, cmoPrimUnfm80, uint8_t); DECLARE_SOA_COLUMN(CmoFV0AUnfm80, cmoFV0AUnfm80, uint8_t); DECLARE_SOA_COLUMN(CmoFT0AUnfm80, cmoFT0AUnfm80, uint8_t); diff --git a/DPG/Tasks/AOTTrack/derivedDataCreatorD0Calibration.cxx b/DPG/Tasks/AOTTrack/derivedDataCreatorD0Calibration.cxx index 4e41c0e44c6..7a45020ba65 100644 --- a/DPG/Tasks/AOTTrack/derivedDataCreatorD0Calibration.cxx +++ b/DPG/Tasks/AOTTrack/derivedDataCreatorD0Calibration.cxx @@ -60,8 +60,8 @@ #include #include #include -#include #include +#include #include using namespace o2; @@ -192,8 +192,18 @@ struct DerivedDataCreatorD0Calibration { TTrackMeanOccs const&, TTrackQa const&) { - std::map selectedCollisions; // map with indices of selected collisions (key: original AOD Collision table index, value: D0 collision index) - std::map selectedTracks; // map with indices of selected tracks (key: original AOD Track table index, value: D0 daughter track index) + std::unordered_map selectedCollisions; // map with indices of selected collisions (key: original AOD Collision table index, value: D0 collision index) + std::unordered_map selectedTracks; // map with indices of selected tracks (key: original AOD Track table index, value: D0 daughter track index) + + std::vector bdtScoresD0; + std::vector bdtScoresD0bar; + std::vector featuresCandD0; + std::vector featuresCandD0bar; + bdtScoresD0.resize(3); + bdtScoresD0bar.resize(3); + featuresCandD0.resize(16); + featuresCandD0bar.resize(16); + auto const& cuts = cfgCandCuts.topologicalCuts; for (auto const& collision : collisions) { @@ -217,6 +227,13 @@ struct DerivedDataCreatorD0Calibration { eventIR.setFromLong(bc.globalBC()); auto groupedTrackIndices = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); + std::vector dcaNeg; // we cache the negative track, to avoid the same computation for each positive track + dcaNeg.resize(groupedTrackIndices.size()); + std::vector isSelectedNeg(groupedTrackIndices.size(), true); + std::vector isPropagatedNeg(groupedTrackIndices.size(), false); + std::vector pidTrackNegKaon(groupedTrackIndices.size(), -1); + std::vector pidTrackNegPion(groupedTrackIndices.size(), -1); + for (auto const& trackIndexPos : groupedTrackIndices) { auto trackPos = trackIndexPos.template track_as(); // track selections @@ -255,27 +272,40 @@ struct DerivedDataCreatorD0Calibration { pidTrackPosPion = selectorPion.statusTpcAndTof(trackPos); } + int iTrack{-1}; for (auto const& trackIndexNeg : groupedTrackIndices) { + iTrack++; + if (!isSelectedNeg[iTrack]) { + continue; + } auto trackNeg = trackIndexNeg.template track_as(); // track selections if (trackNeg.sign() > 0) { // second negative track + isSelectedNeg[iTrack] = false; continue; } if (!trackNeg.isGlobalTrackWoDCA()) { + isSelectedNeg[iTrack] = false; continue; } if (trackNeg.pt() < cfgTrackCuts.ptMin) { + isSelectedNeg[iTrack] = false; continue; } if (std::abs(trackNeg.eta()) > cfgTrackCuts.absEtaMax) { + isSelectedNeg[iTrack] = false; continue; } auto trackParCovNeg = getTrackParCov(trackNeg); - o2::dataformats::DCA dcaNeg(trackNeg.dcaXY(), trackNeg.dcaZ(), trackNeg.cYY(), trackNeg.cZY(), trackNeg.cZZ()); - if (trackNeg.collisionId() != collision.globalIndex()) { - trackParCovNeg.propagateToDCA(primaryVertex, bz, &dcaNeg); + if (!isPropagatedNeg[iTrack]) { + dcaNeg[iTrack] = o2::dataformats::DCA(trackNeg.dcaXY(), trackNeg.dcaZ(), trackNeg.cYY(), trackNeg.cZY(), trackNeg.cZZ()); + if (trackNeg.collisionId() != collision.globalIndex()) { + trackParCovNeg.propagateToDCA(primaryVertex, bz, &dcaNeg[iTrack]); + } + isPropagatedNeg[iTrack] = true; } - if (!isSelectedTrackDca(cfgTrackCuts.binsPt, cfgTrackCuts.limitsDca, trackParCovNeg.getPt(), dcaNeg.getY(), dcaNeg.getZ())) { + if (!isSelectedTrackDca(cfgTrackCuts.binsPt, cfgTrackCuts.limitsDca, trackParCovNeg.getPt(), dcaNeg[iTrack].getY(), dcaNeg[iTrack].getZ())) { + isSelectedNeg[iTrack] = false; continue; } @@ -307,26 +337,24 @@ struct DerivedDataCreatorD0Calibration { } } - int pidTrackNegKaon{-1}; - int pidTrackNegPion{-1}; if (cfgTrackCuts.usePidTpcOnly) { /// kaon TPC PID negative daughter - pidTrackNegKaon = selectorKaon.statusTpc(trackNeg); + pidTrackNegKaon[iTrack] = selectorKaon.statusTpc(trackNeg); /// pion TPC PID negative daughter - pidTrackNegPion = selectorPion.statusTpc(trackNeg); + pidTrackNegPion[iTrack] = selectorPion.statusTpc(trackNeg); } else { /// kaon TPC, TOF PID negative daughter - pidTrackNegKaon = selectorKaon.statusTpcAndTof(trackNeg); + pidTrackNegKaon[iTrack] = selectorKaon.statusTpcAndTof(trackNeg); /// pion TPC, TOF PID negative daughter - pidTrackNegPion = selectorPion.statusTpcAndTof(trackNeg); + pidTrackNegPion[iTrack] = selectorPion.statusTpcAndTof(trackNeg); } // PID uint8_t massHypo{D0MassHypo::D0AndD0Bar}; // both mass hypotheses a priori - if (pidTrackPosPion == TrackSelectorPID::Rejected || pidTrackNegKaon == TrackSelectorPID::Rejected) { + if (pidTrackPosPion == TrackSelectorPID::Rejected || pidTrackNegKaon[iTrack] == TrackSelectorPID::Rejected) { massHypo -= D0MassHypo::D0; // exclude D0 } - if (pidTrackNegPion == TrackSelectorPID::Rejected || pidTrackPosKaon == TrackSelectorPID::Rejected) { + if (pidTrackNegPion[iTrack] == TrackSelectorPID::Rejected || pidTrackPosKaon == TrackSelectorPID::Rejected) { massHypo -= D0MassHypo::D0Bar; // exclude D0Bar } if (massHypo == 0) { @@ -334,20 +362,20 @@ struct DerivedDataCreatorD0Calibration { } // d0xd0 - if (dcaPos.getY() * dcaNeg.getY() > cfgCandCuts.topologicalCuts->get(ptBinNoVtxD0, "max d0d0")) { + if (dcaPos.getY() * dcaNeg[iTrack].getY() > cuts->get(ptBinNoVtxD0, CutsIds::MaxD0D0)) { continue; } // invariant mass if (massHypo == D0MassHypo::D0 || massHypo == D0MassHypo::D0AndD0Bar) { float invMassNoVtxD0 = RecoDecay::m(std::array{trackPos.pVector(), trackNeg.pVector()}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus}); - if (std::abs(invMassNoVtxD0 - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinNoVtxD0, "delta inv. mass") + invMassTolerance) { + if (std::abs(invMassNoVtxD0 - o2::constants::physics::MassD0) > cuts->get(ptBinNoVtxD0, CutsIds::DeltaMass) + invMassTolerance) { massHypo -= D0MassHypo::D0; } } if (massHypo >= D0MassHypo::D0Bar) { float invMassNoVtxD0bar = RecoDecay::m(std::array{trackNeg.pVector(), trackPos.pVector()}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus}); - if (std::abs(invMassNoVtxD0bar - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinNoVtxD0, "delta inv. mass") + invMassTolerance) { + if (std::abs(invMassNoVtxD0bar - o2::constants::physics::MassD0) > cuts->get(ptBinNoVtxD0, CutsIds::DeltaMass) + invMassTolerance) { massHypo -= D0MassHypo::D0Bar; } } @@ -383,64 +411,72 @@ struct DerivedDataCreatorD0Calibration { } // d0xd0 - if (dcaPos.getY() * dcaNeg.getY() > cfgCandCuts.topologicalCuts->get(ptBinD0, "max d0d0")) { + if (dcaPos.getY() * dcaNeg[iTrack].getY() > cuts->get(ptBinD0, CutsIds::MaxD0D0)) { continue; } // cospa float cosPaD0 = RecoDecay::cpa(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, pVecD0); - if (cosPaD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min cos pointing angle")) { + if (cosPaD0 < cuts->get(ptBinD0, CutsIds::MinCosPointingAngle)) { continue; } // pointing angle float paD0 = std::acos(cosPaD0); - if (paD0 > cfgCandCuts.topologicalCuts->get(ptBinD0, "max pointing angle")) { + if (paD0 > cuts->get(ptBinD0, CutsIds::MaxPointingAngle)) { continue; } // cospa XY float cosPaXYD0 = RecoDecay::cpaXY(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, pVecD0); - if (cosPaXYD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min cos pointing angle XY")) { + if (cosPaXYD0 < cuts->get(ptBinD0, CutsIds::MinCosPointingAngleXY)) { continue; } // pointing angle XY float paXYD0 = std::acos(cosPaXYD0); - if (paXYD0 > cfgCandCuts.topologicalCuts->get(ptBinD0, "max pointing angle XY")) { + if (paXYD0 > cuts->get(ptBinD0, CutsIds::MaxPointingAngleXY)) { continue; } // decay length float decLenD0 = RecoDecay::distance(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex); - if (decLenD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min decay length")) { + if (decLenD0 < cuts->get(ptBinD0, CutsIds::MinDecayLength)) { continue; } // decay length XY float decLenXYD0 = RecoDecay::distanceXY(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex); - if (decLenXYD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min decay length XY")) { + if (decLenXYD0 < cuts->get(ptBinD0, CutsIds::MinDecayLengthXY)) { continue; } // normalised decay length float phi{0.f}, theta{0.f}; getPointDirection(std::array{primaryVertex.getX(), primaryVertex.getY(), primaryVertex.getZ()}, secondaryVertex, phi, theta); float errorDecayLengthD0 = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); - if (decLenD0 / errorDecayLengthD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min norm decay length")) { + if (decLenD0 / errorDecayLengthD0 < cuts->get(ptBinD0, CutsIds::MinNormDecayLength)) { continue; } // normalised decay length XY float errorDecayLengthXYD0 = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.f) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.f)); - if (decLenXYD0 / errorDecayLengthXYD0 < cfgCandCuts.topologicalCuts->get(ptBinD0, "min norm decay length XY")) { + if (decLenXYD0 / errorDecayLengthXYD0 < cuts->get(ptBinD0, CutsIds::MinNormDecayLengthXY)) { continue; } float invMassD0{0.f}, invMassD0bar{0.f}; float cosThetaStarD0{0.f}, cosThetaStarD0bar{0.f}; - std::vector bdtScoresD0{0.f, 1.f, 1.f}, bdtScoresD0bar{0.f, 1.f, 1.f}; // always selected a priori + // always selected a priori + bdtScoresD0[0] = 0.f; + bdtScoresD0[1] = 1.f; + bdtScoresD0[2] = 1.f; + bdtScoresD0bar[0] = 0.f; + bdtScoresD0bar[1] = 1.f; + bdtScoresD0bar[2] = 1.f; if (massHypo == D0MassHypo::D0 || massHypo == D0MassHypo::D0AndD0Bar) { invMassD0 = RecoDecay::m(std::array{pVecPos, pVecNeg}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus}); - if (std::abs(invMassD0 - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinD0, "delta inv. mass")) { + if (std::abs(invMassD0 - o2::constants::physics::MassD0) > cuts->get(ptBinD0, CutsIds::DeltaMass)) { massHypo -= D0MassHypo::D0; - bdtScoresD0 = std::vector{1.f, 0.f, 0.f}; + bdtScoresD0[0] = 1.f; + bdtScoresD0[1] = 0.f; + bdtScoresD0[2] = 0.f; } else { // apply BDT models if (cfgMl.apply) { - std::vector featuresCandD0 = {dcaPos.getY(), dcaNeg.getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg.getY(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi())}; + featuresCandD0 = {dcaPos.getY(), dcaNeg[iTrack].getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg[iTrack].getY(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi())}; if (!mlResponse.isSelectedMl(featuresCandD0, ptD0, bdtScoresD0)) { massHypo -= D0MassHypo::D0; } else { // selected, we compute cost* @@ -457,13 +493,15 @@ struct DerivedDataCreatorD0Calibration { } if (massHypo >= D0MassHypo::D0Bar) { invMassD0bar = RecoDecay::m(std::array{pVecNeg, pVecPos}, std::array{o2::constants::physics::MassPiPlus, o2::constants::physics::MassKPlus}); - if (std::abs(invMassD0bar - o2::constants::physics::MassD0) > cfgCandCuts.topologicalCuts->get(ptBinD0, "delta inv. mass")) { + if (std::abs(invMassD0bar - o2::constants::physics::MassD0) > cuts->get(ptBinD0, CutsIds::DeltaMass)) { massHypo -= D0MassHypo::D0Bar; - bdtScoresD0bar = std::vector{1.f, 0.f, 0.f}; + bdtScoresD0bar[0] = 1.f; + bdtScoresD0bar[1] = 0.f; + bdtScoresD0bar[2] = 0.f; } else { // apply BDT models if (cfgMl.apply) { - std::vector featuresCandD0bar = {dcaPos.getY(), dcaNeg.getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg.getY(), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi())}; + featuresCandD0bar = {dcaPos.getY(), dcaNeg[iTrack].getY(), chi2PCA, cosPaD0, cosPaXYD0, decLenXYD0, decLenD0, dcaPos.getY() * dcaNeg[iTrack].getY(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaPi(), trackPos.tofNSigmaPi()), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaKa(), trackNeg.tofNSigmaKa()), trackPos.tpcNSigmaPi(), trackPos.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackPos.tpcNSigmaKa(), trackPos.tofNSigmaKa()), trackNeg.tpcNSigmaPi(), trackNeg.tpcNSigmaKa(), aod::pid_tpc_tof_utils::combineNSigma(trackNeg.tpcNSigmaPi(), trackNeg.tofNSigmaPi())}; if (!mlResponse.isSelectedMl(featuresCandD0bar, ptD0, bdtScoresD0bar)) { massHypo -= D0MassHypo::D0Bar; } else { // selected, we compute cost* @@ -823,8 +861,8 @@ struct DerivedDataCreatorD0Calibration { deltaRefGloParamQ2Pt, deltaTOFdX, deltaTOFdZ, - o2::math_utils::detail::truncateFloatFraction(dcaNeg.getY(), precisionDcaMask), - o2::math_utils::detail::truncateFloatFraction(dcaNeg.getZ(), precisionDcaMask), + o2::math_utils::detail::truncateFloatFraction(dcaNeg[iTrack].getY(), precisionDcaMask), + o2::math_utils::detail::truncateFloatFraction(dcaNeg[iTrack].getZ(), precisionDcaMask), getCompressedNumSigmaPid(trackNeg.tpcNSigmaPi()), getCompressedNumSigmaPid(trackNeg.tpcNSigmaKa()), getCompressedNumSigmaPid(trackNeg.tofNSigmaPi()),