From 8b19ab7e36530d2d8b8826fb5c7ce9dce46372a0 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Tue, 31 Dec 2024 16:09:31 +0100 Subject: [PATCH 1/4] Remove dependence on lambdakzerobuilder --- EventFiltering/PWGHF/HFFilter.cxx | 87 +++----- EventFiltering/PWGHF/HFFilterHelpers.h | 282 +++++++++++++++++++++---- 2 files changed, 268 insertions(+), 101 deletions(-) diff --git a/EventFiltering/PWGHF/HFFilter.cxx b/EventFiltering/PWGHF/HFFilter.cxx index 4b414ed814b..a23cdd31f73 100644 --- a/EventFiltering/PWGHF/HFFilter.cxx +++ b/EventFiltering/PWGHF/HFFilter.cxx @@ -155,6 +155,7 @@ struct HfFilter { // Main struct for HF triggers o2::vertexing::DCAFitterN<3> df3; // fitter for Charm Hadron vertex (3-prong vertex fitter) o2::vertexing::DCAFitterN<2> dfB; // fitter for Beauty Hadron vertex (2-prong vertex fitter) o2::vertexing::DCAFitterN<3> dfBtoDstar; // fitter for Beauty Hadron to D* vertex (3-prong vertex fitter) + o2::vertexing::DCAFitterN<2> dfV0; // fitter for V0s (3-prong vertex fitter) HistogramRegistry registry{"registry"}; std::shared_ptr hProcessedEvents; @@ -215,6 +216,8 @@ struct HfFilter { // Main struct for HF triggers helper.setPtRangeSoftPiSigmaC(ptCuts->get(0u, 4u), ptCuts->get(1u, 4u)); helper.setPtDeltaMassRangeSigmaC(cutsPtDeltaMassCharmReso->get(0u, 6u), cutsPtDeltaMassCharmReso->get(1u, 6u), cutsPtDeltaMassCharmReso->get(0u, 7u), cutsPtDeltaMassCharmReso->get(1u, 7u), cutsPtDeltaMassCharmReso->get(0u, 8u), cutsPtDeltaMassCharmReso->get(1u, 8u), cutsPtDeltaMassCharmReso->get(0u, 9u), cutsPtDeltaMassCharmReso->get(1u, 9u), cutsPtDeltaMassCharmReso->get(2u, 6u), cutsPtDeltaMassCharmReso->get(2u, 7u), cutsPtDeltaMassCharmReso->get(2u, 8u), cutsPtDeltaMassCharmReso->get(2u, 9u)); helper.setPtRangeSoftKaonXicResoToSigmaC(ptCuts->get(0u, 5u), ptCuts->get(1u, 5u)); + helper.setVtxConfiguration(dfV0, true); // (DCAFitterN, useAbsDCA) + dfV0.setMatCorrType(matCorr); if (activateSecVtxForB) { helper.setVtxConfiguration(df2, false); // (DCAFitterN, useAbsDCA) helper.setVtxConfiguration(df3, false); @@ -327,13 +330,14 @@ struct HfFilter { // Main struct for HF triggers using BigTracksMCPID = soa::Join; using BigTracksPID = soa::Join; + using TracksIUPID = soa::Join; using CollsWithEvSel = soa::Join; using Hf2ProngsWithMl = soa::Join; using Hf3ProngsWithMl = soa::Join; Preslice trackIndicesPerCollision = aod::track_association::collisionId; - Preslice v0sPerCollision = aod::v0data::collisionId; + Preslice v0sPerCollision = aod::v0data::collisionId; Preslice hf2ProngPerCollision = aod::track_association::collisionId; Preslice hf3ProngPerCollision = aod::track_association::collisionId; Preslice cascPerCollision = aod::cascdata::collisionId; @@ -341,12 +345,13 @@ struct HfFilter { // Main struct for HF triggers void process(CollsWithEvSel const& collisions, aod::BCsWithTimestamps const&, - aod::V0Datas const& v0s, + aod::V0s const& v0s, aod::CascDatas const& cascades, Hf2ProngsWithMl const& cand2Prongs, Hf3ProngsWithMl const& cand3Prongs, aod::TrackAssoc const& trackIndices, BigTracksPID const&, + TracksIUPID const& tracksIU, aod::V0PhotonsKF const& photons, aod::V0Legs const&) { @@ -702,30 +707,21 @@ struct HfFilter { // Main struct for HF triggers if (!keepEvent[kV0Charm2P] && isSignalTagged && (TESTBIT(selD0, 0) || TESTBIT(selD0, 1))) { auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, thisCollId); for (const auto& v0 : v0sThisCollision) { - auto posTrack = v0.posTrack_as(); - auto negTrack = v0.negTrack_as(); - auto selV0 = helper.isSelectedV0(v0, std::array{posTrack, negTrack}, collision, activateQA, hV0Selected, hArmPod); + V0Cand v0Cand; + if (!buildV0(v0, tracksIU, collision, dfV0, std::vector{cand2Prong.prong0Id(), cand2Prong.prong1Id()}, v0Cand)) { + continue; + } + auto selV0 = helper.isSelectedV0(v0Cand, activateQA, hV0Selected, hArmPod); if (!selV0) { continue; } - // propagate to PV - gpu::gpustd::array dcaInfo; - std::array pVecV0 = {v0.px(), v0.py(), v0.pz()}; - std::array pVecV0Orig = {v0.px(), v0.py(), v0.pz()}; - std::array posVecV0 = {v0.x(), v0.y(), v0.z()}; if (!keepEvent[kV0Charm2P] && TESTBIT(selV0, kK0S)) { - auto trackParK0 = o2::track::TrackPar(posVecV0, pVecV0Orig, 0, true); - trackParK0.setPID(o2::track::PID::K0); - trackParK0.setAbsCharge(0); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParK0, 2.f, matCorr, &dcaInfo); - getPxPyPz(trackParK0, pVecV0); - // we first look for a D*+ for (const auto& trackBachelorId : trackIdsThisCollision) { // start loop over tracks auto trackBachelor = trackBachelorId.track_as(); - if (trackBachelor.globalIndex() == trackPos.globalIndex() || trackBachelor.globalIndex() == trackNeg.globalIndex()) { + if (trackBachelor.globalIndex() == trackPos.globalIndex() || trackBachelor.globalIndex() == trackNeg.globalIndex() || trackBachelor.globalIndex() == v0.posTrackId() || trackBachelor.globalIndex() == v0.negTrackId()) { continue; } @@ -756,10 +752,10 @@ struct HfFilter { // Main struct for HF triggers if (activateQA) { hMassVsPtC[kNCharmParticles]->Fill(ptDStarCand, massDiffDstar); } - auto pVecReso2Prong = RecoDecay::pVec(pVecDStarCand, pVecV0); + auto pVecReso2Prong = RecoDecay::pVec(pVecDStarCand, v0Cand.mom); auto ptCand = RecoDecay::pt(pVecReso2Prong); if (ptCand > cutsPtDeltaMassCharmReso->get(2u, 3u)) { - auto massDStarK0S = RecoDecay::m(std::array{pVecPos, pVecNeg, pVecBachelor, pVecV0}, std::array{massDausD0[0], massDausD0[1], massPi, massK0S}); + auto massDStarK0S = RecoDecay::m(std::array{pVecPos, pVecNeg, pVecBachelor, v0Cand.mom}, std::array{massDausD0[0], massDausD0[1], massPi, massK0S}); auto massDiffDsReso = massDStarK0S - massDStarCand; if (cutsPtDeltaMassCharmReso->get(0u, 3u) < massDiffDsReso && massDiffDsReso < cutsPtDeltaMassCharmReso->get(1u, 3u)) { if (activateQA) { @@ -775,24 +771,19 @@ struct HfFilter { // Main struct for HF triggers } } if (!keepEvent[kV0Charm2P] && (TESTBIT(selV0, kLambda) || TESTBIT(selV0, kAntiLambda))) { // Xic(3055) and Xic(3080) --> since it occupies only a small bandwidth, we might want to keep also wrong sign pairs - auto trackParLambda = o2::track::TrackPar(posVecV0, pVecV0Orig, 0, true); - trackParLambda.setAbsCharge(0); - trackParLambda.setPID(o2::track::PID::Lambda); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParLambda, 2.f, matCorr, &dcaInfo); - getPxPyPz(trackParLambda, pVecV0); float massXicStarCand{-999.}, massXicStarBarCand{-999.}; float massDiffXicStarCand{-999.}, massDiffXicStarBarCand{-999.}; bool isRightSignXicStar{false}, isRightSignXicStarBar{false}; - auto pVecReso2Prong = RecoDecay::pVec(pVec2Prong, pVecV0); + auto pVecReso2Prong = RecoDecay::pVec(pVec2Prong, v0Cand.mom); auto ptCand = RecoDecay::pt(pVecReso2Prong); if (ptCand > cutsPtDeltaMassCharmReso->get(2u, 5u)) { if (TESTBIT(selD0, 0)) { - massXicStarCand = RecoDecay::m(std::array{pVecPos, pVecNeg, pVecV0}, std::array{massPi, massKa, massLambda}); + massXicStarCand = RecoDecay::m(std::array{pVecPos, pVecNeg, v0Cand.mom}, std::array{massPi, massKa, massLambda}); massDiffXicStarCand = massXicStarCand - massD0Cand; isRightSignXicStar = TESTBIT(selV0, kLambda); // right sign if Lambda } if (TESTBIT(selD0, 1)) { - massXicStarBarCand = RecoDecay::m(std::array{pVecPos, pVecNeg, pVecV0}, std::array{massKa, massPi, massLambda}); + massXicStarBarCand = RecoDecay::m(std::array{pVecPos, pVecNeg, v0Cand.mom}, std::array{massKa, massPi, massLambda}); massDiffXicStarBarCand = massXicStarBarCand - massD0BarCand; isRightSignXicStarBar = TESTBIT(selV0, kAntiLambda); // right sign if AntiLambda } @@ -1249,28 +1240,21 @@ struct HfFilter { // Main struct for HF triggers if ((!keepEvent[kV0Charm3P] && isGoodDPlus) || (!keepEvent[kSigmaC0K0] && (isGoodLcToPKPi || isGoodLcToPiKP))) { for (const auto& v0 : v0sThisCollision) { - auto posTrack = v0.posTrack_as(); - auto negTrack = v0.negTrack_as(); - auto selV0 = helper.isSelectedV0(v0, std::array{posTrack, negTrack}, collision, activateQA, hV0Selected, hArmPod); + V0Cand v0Cand; + if (!buildV0(v0, tracksIU, collision, dfV0, std::vector{cand3Prong.prong0Id(), cand3Prong.prong1Id(), cand3Prong.prong2Id()}, v0Cand)) { + continue; + } + auto selV0 = helper.isSelectedV0(v0Cand, activateQA, hV0Selected, hArmPod); if (!selV0) { continue; } - gpu::gpustd::array dcaInfo; - std::array pVecV0Orig = {v0.px(), v0.py(), v0.pz()}; - std::array pVecV0 = {v0.px(), v0.py(), v0.pz()}; - std::array posVecV0 = {v0.x(), v0.y(), v0.z()}; // we pair D+ with V0 if (!keepEvent[kV0Charm3P] && isGoodDPlus) { if (!keepEvent[kV0Charm3P] && TESTBIT(selV0, kK0S)) { // Ds2* - auto trackParK0S = o2::track::TrackPar(posVecV0, pVecV0Orig, 0, true); - trackParK0S.setAbsCharge(0); - trackParK0S.setPID(o2::track::PID::K0); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParK0S, 2.f, matCorr, &dcaInfo); - getPxPyPz(trackParK0S, pVecV0); - auto massDsStarCand = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecV0}, std::array{massPi, massKa, massPi, massK0S}); + auto massDsStarCand = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, v0Cand.mom}, std::array{massPi, massKa, massPi, massK0S}); auto massDiffDsStar = massDsStarCand - massDPlusCand; - auto pVecReso3Prong = RecoDecay::pVec(pVec3Prong, pVecV0); + auto pVecReso3Prong = RecoDecay::pVec(pVec3Prong, v0Cand.mom); auto ptCand = RecoDecay::pt(pVecReso3Prong); if (ptCand > cutsPtDeltaMassCharmReso->get(2u, 4u)) { if (cutsPtDeltaMassCharmReso->get(0u, 4u) < massDiffDsStar && massDiffDsStar < cutsPtDeltaMassCharmReso->get(1u, 4u)) { @@ -1282,15 +1266,10 @@ struct HfFilter { // Main struct for HF triggers } } if (!keepEvent[kV0Charm3P] && (TESTBIT(selV0, kLambda) || TESTBIT(selV0, kAntiLambda))) { // Xic(3055) and Xic(3080) --> since it occupies only a small bandwidth, we might want to keep also wrong sign pairs - auto trackParLambda = o2::track::TrackPar(posVecV0, pVecV0Orig, 0, true); - trackParLambda.setAbsCharge(0); - trackParLambda.setPID(o2::track::PID::Lambda); - o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParLambda, 2.f, matCorr, &dcaInfo); - getPxPyPz(trackParLambda, pVecV0); - auto pVecReso3Prong = RecoDecay::pVec(pVec3Prong, pVecV0); + auto pVecReso3Prong = RecoDecay::pVec(pVec3Prong, v0Cand.mom); auto ptCand = RecoDecay::pt(pVecReso3Prong); if (ptCand > cutsPtDeltaMassCharmReso->get(2u, 5u)) { - auto massXicStarCand = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecV0}, std::array{massPi, massKa, massPi, massLambda}); + auto massXicStarCand = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, v0Cand.mom}, std::array{massPi, massKa, massPi, massLambda}); auto massDiffXicStar = massXicStarCand - massDPlusCand; bool isRightSign = ((TESTBIT(selV0, kLambda) && sign3Prong > 0) || (TESTBIT(selV0, kAntiLambda) && sign3Prong < 0)); if (cutsPtDeltaMassCharmReso->get(0u, 5u) < massDiffXicStar && massDiffXicStar < cutsPtDeltaMassCharmReso->get(1u, 5u)) { @@ -1319,8 +1298,8 @@ struct HfFilter { // Main struct for HF triggers auto globalIndexSoftPi = trackSoftPi.globalIndex(); // exclude tracks already used to build the 3-prong candidate - if (globalIndexSoftPi == trackFirst.globalIndex() || globalIndexSoftPi == trackSecond.globalIndex() || globalIndexSoftPi == trackThird.globalIndex()) { - // do not consider as candidate soft pion a track already used to build the current 3-prong candidate + if (globalIndexSoftPi == trackFirst.globalIndex() || globalIndexSoftPi == trackSecond.globalIndex() || globalIndexSoftPi == trackThird.globalIndex() || globalIndexSoftPi == v0.posTrackId() || globalIndexSoftPi == v0.negTrackId()) { + // do not consider as candidate soft pion a track already used to build the current 3-prong candidate / V0 candidate continue; } @@ -1353,17 +1332,15 @@ struct HfFilter { // Main struct for HF triggers /// and keep it only if it is in the correct mass range float massSigmaCPKPi{-999.}, massSigmaCPiKP{-999.}, deltaMassXicResoPKPi{-999.}, deltaMassXicResoPiKP{-999.}; - std::array pVecPiPosK0s = posTrack.pVector(); - std::array pVecPiNegK0s = negTrack.pVector(); - float ptSigmaCKaon = RecoDecay::pt(pVecSigmaC, pVecPiPosK0s, pVecPiNegK0s); + float ptSigmaCKaon = RecoDecay::pt(pVecSigmaC, v0Cand.mom); if (ptSigmaCKaon > cutsPtDeltaMassCharmReso->get(2u, 10u)) { if (TESTBIT(whichSigmaC, 0)) { massSigmaCPKPi = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecSoftPi}, std::array{massProton, massKa, massPi, massPi}); - deltaMassXicResoPKPi = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecSoftPi, pVecPiPosK0s, pVecPiNegK0s}, std::array{massProton, massKa, massPi, massPi, massPi, massPi}) - massSigmaCPKPi; + deltaMassXicResoPKPi = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecSoftPi, v0Cand.mom}, std::array{massProton, massKa, massPi, massPi, massK0S}) - massSigmaCPKPi; } if (TESTBIT(whichSigmaC, 1)) { massSigmaCPiKP = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecSoftPi}, std::array{massPi, massKa, massProton, massPi}); - deltaMassXicResoPiKP = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecSoftPi, pVecPiPosK0s, pVecPiNegK0s}, std::array{massPi, massKa, massProton, massPi, massPi, massPi}) - massSigmaCPiKP; + deltaMassXicResoPiKP = RecoDecay::m(std::array{pVecFirst, pVecSecond, pVecThird, pVecSoftPi, v0Cand.mom}, std::array{massPi, massKa, massProton, massPi, massK0S}) - massSigmaCPiKP; } bool isPKPiOk = (cutsPtDeltaMassCharmReso->get(0u, 10u) < deltaMassXicResoPKPi && deltaMassXicResoPKPi < cutsPtDeltaMassCharmReso->get(1u, 10u)); diff --git a/EventFiltering/PWGHF/HFFilterHelpers.h b/EventFiltering/PWGHF/HFFilterHelpers.h index fe38a740903..5d1a2b1f49b 100644 --- a/EventFiltering/PWGHF/HFFilterHelpers.h +++ b/EventFiltering/PWGHF/HFFilterHelpers.h @@ -39,6 +39,8 @@ #include "CommonConstants/MathConstants.h" #include "CommonConstants/PhysicsConstants.h" #include "DataFormatsTPC/BetheBlochAleph.h" +#include "DCAFitter/DCAFitterN.h" +#include "DetectorsBase/Propagator.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/DataTypes.h" @@ -143,6 +145,39 @@ enum HfVtxStage : uint8_t { kNHfVtxStage }; +// Helper struct to pass V0 informations +struct V0Cand { + std::array mom; + float etaPos; + float etaNeg; + float ptPos; + float ptNeg; + float pinTpcPos; + float pinTpcNeg; + float nClsFoundTpcPos; + float nClsFoundTpcNeg; + float signalTpcPos; + float signalTpcNeg; + float nClsTpcNeg; + float v0cosPA; + float dcav0topv; + float dcaV0daughters; + float dcapostopv; + float dcanegtopv; + float alpha; + float qtarm; + float v0radius; + float mK0Short; + float mLambda; + float mAntiLambda; + float nSigmaPrTpcPos; + float nSigmaPrTofPos; + float nSigmaPrTpcNeg; + float nSigmaPrTofNeg; + bool hasTofPos; + bool hasTofNeg; +}; + static const std::array charmParticleNames{"D0", "Dplus", "Ds", "Lc", "Xic"}; static const std::array beautyParticleNames{"Bplus", "B0toDStar", "B0", "Bs", "Lb", "Xib"}; static const std::array pdgCodesCharm{421, 411, 431, 4122, 4232}; @@ -460,10 +495,10 @@ class HfFilterHelper int8_t isSelectedSigmaCInDeltaMassRange(const T& pTrackSameChargeFirst, const T& pTrackSameChargeSecond, const T& pTrackOppositeCharge, const T& pTrackSoftPi, const float ptSigmaC, const int8_t isSelectedLc, H2 hMassVsPt, const int& activateQA); template int8_t isSelectedXicInMassRange(const T& pTrackSameChargeFirst, const T& pTrackSameChargeSecond, const T& pTrackOppositeCharge, const float& ptXic, const int8_t isSelected, const int& activateQA, H2 hMassVsPt); - template - int8_t isSelectedV0(const V0& v0, const std::array& dauTracks, const Coll& collision, const int& activateQA, H2 hV0Selected, std::array& hArmPod); + template + int8_t isSelectedV0(const V0& v0, const int& activateQA, H2 hV0Selected, std::array& hArmPod); template - inline bool isSelectedPhoton(const Photon& photon, const std::array& dauTracks, const int& activateQA, H2 hV0Selected, std::array& hArmPod); + bool isSelectedPhoton(const Photon& photon, const std::array& dauTracks, const int& activateQA, H2 hV0Selected, std::array& hArmPod); template bool isSelectedCascade(const Casc& casc, const std::array& dauTracks, const Coll& collision); template @@ -473,13 +508,13 @@ class HfFilterHelper template bool isSelectedKaonFromXicResoToSigmaC(const T& track); template - inline bool isSelectedBhadron(T1 const& pVecTrack0, T1 const& pVecTrack1, T2 const& dcaTrack0, T2 const& dcaTrack1, const T3& primVtx, const T4& secVtx, const int whichB); + bool isSelectedBhadron(T1 const& pVecTrack0, T1 const& pVecTrack1, T2 const& dcaTrack0, T2 const& dcaTrack1, const T3& primVtx, const T4& secVtx, const int whichB); template - inline bool isSelectedBhadronInMassRange(T1 const& ptCand, T2 const& massCand, const int whichB); + bool isSelectedBhadronInMassRange(T1 const& ptCand, T2 const& massCand, const int whichB); template - inline bool isSelectedBzeroToDstar(T1 const& pVecTrack0, T1 const& pVecTrack1, T1 const& pVecTrack2, const T2& primVtx, const T3& secVtx); + bool isSelectedBzeroToDstar(T1 const& pVecTrack0, T1 const& pVecTrack1, T1 const& pVecTrack2, const T2& primVtx, const T3& secVtx); template - inline bool isCharmHadronMassInSbRegions(T1 const& massHypo1, T1 const& massHypo2, const float& lowLimitSB, const float& upLimitSB); + bool isCharmHadronMassInSbRegions(T1 const& massHypo1, T1 const& massHypo2, const float& lowLimitSB, const float& upLimitSB); // helpers template @@ -487,7 +522,7 @@ class HfFilterHelper template int computeNumberOfCandidates(std::vector> indices); template - inline int setVtxConfiguration(T1 vertexer, bool useAbsDCA); + int setVtxConfiguration(T1 vertexer, bool useAbsDCA); // PID void setValuesBB(o2::ccdb::CcdbApi& ccdbApi, aod::BCsWithTimestamps::iterator const& bunchCrossing, const std::array& ccdbPaths); @@ -501,14 +536,20 @@ class HfFilterHelper bool isSelectedProton4CharmBaryons(const T& track); // PID + float getTPCSplineCalib(const float tpcPin, const float dEdx, const int& pidSpecies); template - double getTPCSplineCalib(const T& track, const int& pidSpecies); + float getTPCSplineCalib(const T& track, const int& pidSpecies); + float getTPCPostCalib(const float tpcPin, const float tpcNCls, const float eta, const float tpcNSigma, const int& pidSpecies); template float getTPCPostCalib(const T& track, const int& pidSpecies); // helpers template int findBin(T1 const& binsPt, T2 value); + template + std::array alphaAndQtAP(std::array const& momPos, std::array const& momNeg); + template + bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, V0Cand& v0Cand); // selections std::vector mPtBinsTracks{}; // vector of pT bins for single track cuts @@ -1154,14 +1195,12 @@ inline int8_t HfFilterHelper::isSelectedXicInMassRange(const T& pTrackSameCharge /// Basic selection of V0 candidates /// \param v0 is the v0 candidate -/// \param dauTracks is a 2-element array with positive and negative V0 daughter tracks -/// \param collision is the current collision /// \param activateQA flag to fill QA histos /// \param hV0Selected is the pointer to the QA histo for selected V0S /// \param hArmPod is the pointer to an array of QA histo AP plot after selection /// \return an integer passes all cuts -template -inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& dauTracks, const Coll& /*collision*/, const int& activateQA, H2 hV0Selected, std::array& hArmPod) +template +inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const int& activateQA, H2 hV0Selected, std::array& hArmPod) { int8_t isSelected{BIT(kK0S) | BIT(kLambda) | BIT(kAntiLambda)}; @@ -1172,7 +1211,7 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& } // eta of daughters - if (std::fabs(dauTracks[0].eta()) > 1. || std::fabs(dauTracks[1].eta()) > 1.) { // cut all V0 daughters with |eta| > 1. + if (std::fabs(v0.etaPos) > 1. || std::fabs(v0.etaNeg) > 1.) { // cut all V0 daughters with |eta| > 1. if (activateQA > 1) { for (int iV0{kK0S}; iV0 < kNV0; ++iV0) { hV0Selected->Fill(1., iV0); @@ -1182,7 +1221,7 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& } // V0 radius - if (v0.v0radius() < mMinK0sLambdaRadius) { + if (v0.v0radius < mMinK0sLambdaRadius) { for (int iV0{kK0S}; iV0 < kNV0; ++iV0) { CLRBIT(isSelected, iV0); if (activateQA > 1) { @@ -1191,9 +1230,8 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& } } - auto v0CosinePa = v0.v0cosPA(); for (int iV0{kK0S}; iV0 < kNV0; ++iV0) { - if (TESTBIT(isSelected, iV0) && v0CosinePa < mMinK0sLambdaCosinePa) { + if (TESTBIT(isSelected, iV0) && v0.v0cosPA < mMinK0sLambdaCosinePa) { CLRBIT(isSelected, iV0); if (activateQA > 1) { hV0Selected->Fill(3., iV0); @@ -1202,19 +1240,19 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& } // armenteros-podolanski / mass - if (TESTBIT(isSelected, kK0S) && std::fabs(v0.mK0Short() - massK0S) > mDeltaMassK0s) { + if (TESTBIT(isSelected, kK0S) && std::fabs(v0.mK0Short - massK0S) > mDeltaMassK0s) { CLRBIT(isSelected, kK0S); if (activateQA > 1) { hV0Selected->Fill(4., kK0S); } } - if (TESTBIT(isSelected, kLambda) && std::fabs(v0.mLambda() - massLambda) > mDeltaMassLambda) { + if (TESTBIT(isSelected, kLambda) && std::fabs(v0.mLambda - massLambda) > mDeltaMassLambda) { CLRBIT(isSelected, kLambda); if (activateQA > 1) { hV0Selected->Fill(4., kLambda); } } - if (TESTBIT(isSelected, kAntiLambda) && std::fabs(v0.mAntiLambda() - massLambda) > mDeltaMassLambda) { + if (TESTBIT(isSelected, kAntiLambda) && std::fabs(v0.mAntiLambda - massLambda) > mDeltaMassLambda) { CLRBIT(isSelected, kAntiLambda); if (activateQA > 1) { hV0Selected->Fill(4., kAntiLambda); @@ -1223,13 +1261,13 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& // DCA V0 and V0 daughters for (int iV0{kK0S}; iV0 < kNV0; ++iV0) { - if (TESTBIT(isSelected, iV0) && v0.dcav0topv() > 0.1f) { // we want only primary V0s + if (TESTBIT(isSelected, iV0) && v0.dcav0topv > 0.1f) { // we want only primary V0s CLRBIT(isSelected, iV0); if (activateQA > 1) { hV0Selected->Fill(5., iV0); } } - if (TESTBIT(isSelected, iV0) && (v0.dcaV0daughters() > 1.f || std::fabs(v0.dcapostopv()) < 0.05f || std::fabs(v0.dcanegtopv()) < 0.05f)) { + if (TESTBIT(isSelected, iV0) && (v0.dcaV0daughters > 1.f || std::fabs(v0.dcapostopv) < 0.05f || std::fabs(v0.dcanegtopv) < 0.05f)) { CLRBIT(isSelected, iV0); if (activateQA > 1) { hV0Selected->Fill(6., iV0); @@ -1238,25 +1276,29 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& } // PID (Lambda/AntiLambda only) - float nSigmaPrTpc[2] = {dauTracks[0].tpcNSigmaPr(), dauTracks[1].tpcNSigmaPr()}; - float nSigmaPrTof[2] = {dauTracks[0].tofNSigmaPr(), dauTracks[1].tofNSigmaPr()}; + float nSigmaPrTpc[2] = {v0.nSigmaPrTpcPos, v0.nSigmaPrTpcNeg}; + float nSigmaPrTof[2] = {v0.nSigmaPrTofPos, v0.nSigmaPrTofNeg}; + float pInTpc[2] = {v0.pinTpcPos, v0.pinTpcNeg}; + float nClsTpc[2] = {v0.nClsFoundTpcPos, v0.nClsFoundTpcNeg}; + float etaDaus[2] = {v0.etaPos, v0.etaNeg}; + float signalTpc[2] = {v0.signalTpcPos, v0.signalTpcNeg}; if (mTpcPidCalibrationOption == 1) { for (int iDau{0}; iDau < 2; ++iDau) { - nSigmaPrTpc[iDau] = getTPCPostCalib(dauTracks[iDau], kPr); + nSigmaPrTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[2], kPr); } } else if (mTpcPidCalibrationOption == 2) { for (int iDau{0}; iDau < 2; ++iDau) { - nSigmaPrTpc[iDau] = getTPCSplineCalib(dauTracks[iDau], (iDau == 0) ? kPr : kAntiPr); + nSigmaPrTpc[iDau] = getTPCSplineCalib(pInTpc[iDau], signalTpc[iDau], (iDau == 0) ? kPr : kAntiPr); } } - if (TESTBIT(isSelected, kLambda) && ((dauTracks[0].hasTPC() && std::fabs(nSigmaPrTpc[0]) > mMaxNsigmaPrForLambda) || (dauTracks[0].hasTOF() && std::fabs(nSigmaPrTof[0]) > mMaxNsigmaPrForLambda))) { + if (TESTBIT(isSelected, kLambda) && (std::fabs(nSigmaPrTpc[0]) > mMaxNsigmaPrForLambda || (v0.hasTofPos && std::fabs(nSigmaPrTof[0]) > mMaxNsigmaPrForLambda))) { CLRBIT(isSelected, kLambda); if (activateQA > 1) { hV0Selected->Fill(7., kLambda); } } - if (TESTBIT(isSelected, kAntiLambda) && ((dauTracks[1].hasTPC() && std::fabs(nSigmaPrTpc[1]) > mMaxNsigmaPrForLambda) || (dauTracks[1].hasTOF() && std::fabs(nSigmaPrTof[1]) > mMaxNsigmaPrForLambda))) { + if (TESTBIT(isSelected, kAntiLambda) && (std::fabs(nSigmaPrTpc[1]) > mMaxNsigmaPrForLambda || (v0.hasTofNeg && std::fabs(nSigmaPrTof[1]) > mMaxNsigmaPrForLambda))) { CLRBIT(isSelected, kAntiLambda); if (activateQA > 1) { hV0Selected->Fill(7., kAntiLambda); @@ -1266,7 +1308,7 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const std::array& if (activateQA) { for (int iV0{kK0S}; iV0 < kNV0; ++iV0) { if (TESTBIT(isSelected, iV0)) { - hArmPod[iV0]->Fill(v0.alpha(), v0.qtarm()); + hArmPod[iV0]->Fill(v0.alpha, v0.qtarm); if (activateQA > 1) { hV0Selected->Fill(8., iV0); } @@ -1877,7 +1919,20 @@ inline bool isCharmHadronMassInSbRegions(T1 const& massHypo1, T1 const& massHypo /// \param pidSpecies is the particle species to be considered /// \return updated nsigma value for TPC PID template -inline double HfFilterHelper::getTPCSplineCalib(const T& track, const int& pidSpecies) +inline float HfFilterHelper::getTPCSplineCalib(const T& track, const int& pidSpecies) +{ + float tpcPin = track.tpcInnerParam(); + float dEdx = track.tpcSignal(); + + return getTPCSplineCalib(tpcPin, dEdx, pidSpecies); +} + +/// Update the TPC PID baesd on the spline of particles +/// \param tpcPin is the TPC momentum at innermost update +/// \param dEdx is the TPC dEdx +/// \param pidSpecies is the particle species to be considered +/// \return updated nsigma value for TPC PID +inline float HfFilterHelper::getTPCSplineCalib(const float tpcPin, const float dEdx, const int& pidSpecies) { float mMassPar{0.}; if (pidSpecies == kPi || pidSpecies == kAntiPi) { @@ -1894,38 +1949,31 @@ inline double HfFilterHelper::getTPCSplineCalib(const T& track, const int& pidSp } auto bgScaling = 1 / mMassPar; - double expBethe = tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * bgScaling), mBetheBlochPiKaPrDe[pidSpecies][0], mBetheBlochPiKaPrDe[pidSpecies][1], mBetheBlochPiKaPrDe[pidSpecies][2], mBetheBlochPiKaPrDe[pidSpecies][3], mBetheBlochPiKaPrDe[pidSpecies][4]); + double expBethe = tpc::BetheBlochAleph(static_cast(tpcPin * bgScaling), mBetheBlochPiKaPrDe[pidSpecies][0], mBetheBlochPiKaPrDe[pidSpecies][1], mBetheBlochPiKaPrDe[pidSpecies][2], mBetheBlochPiKaPrDe[pidSpecies][3], mBetheBlochPiKaPrDe[pidSpecies][4]); double expSigma = expBethe * mBetheBlochPiKaPrDe[pidSpecies][5]; - return static_cast((track.tpcSignal() - expBethe) / expSigma); + return static_cast((dEdx - expBethe) / expSigma); } /// compute TPC postcalibrated nsigma based on calibration histograms from CCDB -/// \param hCalibMean calibration histograms of mean from CCDB -/// \param hCalibSigma calibration histograms of sigma from CCDB -/// \param track is the track +/// \param tpcPin is the TPC momentum at innermost update +/// \param tpcNCls is the number of found TPC clusters +/// \param eta is the pseudorapidity +/// \param tpcNSigma is the original Nsigma /// \param pidSpecies is the PID species /// \return the corrected Nsigma value for the PID species -template -inline float HfFilterHelper::getTPCPostCalib(const T& track, const int& pidSpecies) +inline float HfFilterHelper::getTPCPostCalib(const float tpcPin, const float tpcNCls, const float eta, const float tpcNSigma, const int& pidSpecies) { - float tpcNCls = track.tpcNClsFound(); - float tpcPin = track.tpcInnerParam(); - float eta = track.eta(); - float tpcNSigma{0.}; int iHist{0}; - if (pidSpecies == kPi) { - tpcNSigma = track.tpcNSigmaPi(); iHist = 0; } else if (pidSpecies == kKa) { - tpcNSigma = track.tpcNSigmaKa(); iHist = 2; } else if (pidSpecies == kPr) { - tpcNSigma = track.tpcNSigmaPr(); iHist = 4; } else { LOG(fatal) << "Wrong PID Species be selected, please check!"; } + if (!mHistMapPiPrKaDe[iHist] || !mHistMapPiPrKaDe[iHist + 1]) { LOGP(warn, "Postcalibration TPC PID histograms not set. Use default Nsigma values."); } @@ -1946,6 +1994,31 @@ inline float HfFilterHelper::getTPCPostCalib(const T& track, const int& pidSpeci return (tpcNSigma - mean) / width; } +/// compute TPC postcalibrated nsigma based on calibration histograms from CCDB +/// \param track is the track +/// \param pidSpecies is the PID species +/// \return the corrected Nsigma value for the PID species +template +inline float HfFilterHelper::getTPCPostCalib(const T& track, const int& pidSpecies) +{ + float tpcNCls = track.tpcNClsFound(); + float tpcPin = track.tpcInnerParam(); + float eta = track.eta(); + float tpcNSigma{-999.}; + + if (pidSpecies == kPi) { + tpcNSigma = track.tpcNSigmaPi(); + } else if (pidSpecies == kKa) { + tpcNSigma = track.tpcNSigmaKa(); + } else if (pidSpecies == kPr) { + tpcNSigma = track.tpcNSigmaPr(); + } else { + LOG(fatal) << "Wrong PID Species be selected, please check!"; + } + + return getTPCPostCalib(tpcPin, tpcNCls, eta, tpcNSigma, pidSpecies); +} + /// Finds pT bin in an array. /// \param bins array of pT bins /// \param value pT @@ -1971,14 +2044,131 @@ inline int HfFilterHelper::setVtxConfiguration(T1 vertexer, bool useAbsDCA) // Fitter initialisation vertexer.setPropagateToPCA(true); vertexer.setMaxR(200.); - vertexer.setMaxDZIni(4.); + vertexer.setMaxDZIni(1.e9); + vertexer.setMaxDXYIni(4.); vertexer.setMinParamChange(1.e-3); vertexer.setMinRelChi2Change(0.9); + vertexer.setMaxChi2(0.9); vertexer.setUseAbsDCA(useAbsDCA); vertexer.setWeightedFinalPCA(false); return 1; } +/// Utility to compute AP alpha and qt +/// \param momPos momentum array of positive daughter +/// \param momNeg momentum array of negative daughter +template +inline std::array alphaAndQtAP(std::array const& momPos, std::array const& momNeg) +{ + float momTotSq = RecoDecay::p2(momPos[0] + momNeg[0], momPos[1] + momNeg[1], momPos[2] + momNeg[2]); + float momTot = std::sqrt(momTotSq); + float lQlNeg = RecoDecay::dotProd(momNeg, std::array{momPos[0] + momNeg[0], momPos[1] + momNeg[1], momPos[2] + momNeg[2]}) / momTot; + float lQlPos = RecoDecay::dotProd(momPos, std::array{momPos[0] + momNeg[0], momPos[1] + momNeg[1], momPos[2] + momNeg[2]}) / momTot; + float alpha = (lQlPos - lQlNeg) / (lQlPos + lQlNeg); + float qtarm = std::sqrt(RecoDecay::p2(momNeg) - lQlNeg * lQlNeg / momTotSq); + + std::array alphaAndQt = {alpha, qtarm}; + return alphaAndQt; +} + +/// build V0 candidate from table with track indices +/// \param v0Indices V0 candidate from AO2D table (track indices) +/// \param tracks track table +/// \param collision collision +/// \param dcaFitter DCA fitter to be used +/// \param vetoedTrackIds vector with forbidden track indices, if any +template +inline bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, V0Cand& v0Cand) +{ + auto trackPos = tracks.rawIteratorAt(v0Indices.posTrackId()); + auto trackNeg = tracks.rawIteratorAt(v0Indices.negTrackId()); + + // minimal track cuts + if (!trackPos.hasTPC() || !trackNeg.hasTPC()) { + return false; + } + if (!(trackPos.trackType() & o2::aod::track::TPCrefit) || !(trackNeg.trackType() & o2::aod::track::TPCrefit)) { + return false; + } + if (trackPos.tpcNClsCrossedRows() < 50 || trackNeg.tpcNClsCrossedRows() < 50) { + return false; + } + + if (std::find(vetoedTrackIds.begin(), vetoedTrackIds.end(), trackPos.globalIndex()) != vetoedTrackIds.end() || std::find(vetoedTrackIds.begin(), vetoedTrackIds.end(), trackNeg.globalIndex()) != vetoedTrackIds.end()) { + return false; + } + + auto trackParCovPos = getTrackParCov(trackPos); + auto trackParCovNeg = getTrackParCov(trackNeg); + std::array primVtx = {collision.posX(), collision.posY(), collision.posZ()}; + gpu::gpustd::array dcaInfoPos, dcaInfoNeg; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovPos, 2.f, dcaFitter.getMatCorrType(), &dcaInfoPos); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCovNeg, 2.f, dcaFitter.getMatCorrType(), &dcaInfoNeg); + + // reconstruct vertex + int nCand = 0; + try { + nCand = dcaFitter.process(trackParCovPos, trackParCovNeg); + } catch (...) { + LOG(error) << "Exception caught in DCA fitter process call!"; + return false; + } + if (nCand == 0) { + return false; + } + + // compute candidate momentum from tracks propagated to decay vertex + auto& trackPosProp = dcaFitter.getTrack(0); + auto& trackNegProp = dcaFitter.getTrack(1); + std::array momPos{}, momNeg{}; + trackPosProp.getPxPyPzGlo(momPos); + trackNegProp.getPxPyPzGlo(momNeg); + v0Cand.mom = RecoDecay::pVec(momPos, momNeg); + + // fill V0 quantities + v0Cand.dcapostopv = dcaInfoPos[0]; + v0Cand.dcanegtopv = dcaInfoNeg[0]; + v0Cand.ptPos = RecoDecay::pt(momPos); + v0Cand.ptNeg = RecoDecay::pt(momNeg); + v0Cand.pinTpcPos = trackPos.tpcInnerParam(); + v0Cand.pinTpcNeg = trackNeg.tpcInnerParam(); + v0Cand.nClsFoundTpcPos = trackPos.tpcNClsFound(); + v0Cand.nClsFoundTpcNeg = trackNeg.tpcNClsFound(); + v0Cand.signalTpcPos = trackPos.tpcSignal(); + v0Cand.signalTpcNeg = trackNeg.tpcSignal(); + v0Cand.etaPos = RecoDecay::eta(momPos); + v0Cand.etaNeg = RecoDecay::eta(momNeg); + v0Cand.dcaV0daughters = std::sqrt(dcaFitter.getChi2AtPCACandidate()); + + const auto& vtx = dcaFitter.getPCACandidate(); + v0Cand.v0radius = std::hypot(vtx[0], vtx[1]); + v0Cand.v0cosPA = RecoDecay::cpa(primVtx, vtx, v0Cand.mom); + + auto trackParV0 = dcaFitter.createParentTrackParCov(); + trackParV0.setAbsCharge(0); + trackParV0.setPID(o2::track::PID::K0); + gpu::gpustd::array dcaInfoV0; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParV0, 2.f, dcaFitter.getMatCorrType(), &dcaInfoV0); + v0Cand.dcav0topv = dcaInfoV0[0]; + + v0Cand.mK0Short = RecoDecay::m(std::array{momPos, momNeg}, std::array{massPi, massPi}); + v0Cand.mLambda = RecoDecay::m(std::array{momPos, momNeg}, std::array{massProton, massPi}); + v0Cand.mAntiLambda = RecoDecay::m(std::array{momPos, momNeg}, std::array{massPi, massProton}); + + auto alphaAndQt = alphaAndQtAP(momPos, momNeg); + v0Cand.alpha = alphaAndQt[0]; + v0Cand.qtarm = alphaAndQt[1]; + + v0Cand.hasTofPos = trackPos.hasTOF(); + v0Cand.hasTofNeg = trackNeg.hasTOF(); + v0Cand.nSigmaPrTpcPos = trackPos.tpcNSigmaPr(); + v0Cand.nSigmaPrTofPos = trackPos.tofNSigmaPr(); + v0Cand.nSigmaPrTpcNeg = trackNeg.tpcNSigmaPr(); + v0Cand.nSigmaPrTofNeg = trackNeg.tofNSigmaPr(); + + return true; +} + } // namespace hffilters /// definition of tables From f0473dfc1d6ce39f37b2987054b5418cc1489fba Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Wed, 1 Jan 2025 18:54:20 +0100 Subject: [PATCH 2/4] Remove dependence on cascadebuilder --- EventFiltering/PWGHF/HFFilter.cxx | 51 +++-- EventFiltering/PWGHF/HFFilterHelpers.h | 247 +++++++++++++++++++++---- 2 files changed, 240 insertions(+), 58 deletions(-) diff --git a/EventFiltering/PWGHF/HFFilter.cxx b/EventFiltering/PWGHF/HFFilter.cxx index a23cdd31f73..7436e21da99 100644 --- a/EventFiltering/PWGHF/HFFilter.cxx +++ b/EventFiltering/PWGHF/HFFilter.cxx @@ -155,7 +155,7 @@ struct HfFilter { // Main struct for HF triggers o2::vertexing::DCAFitterN<3> df3; // fitter for Charm Hadron vertex (3-prong vertex fitter) o2::vertexing::DCAFitterN<2> dfB; // fitter for Beauty Hadron vertex (2-prong vertex fitter) o2::vertexing::DCAFitterN<3> dfBtoDstar; // fitter for Beauty Hadron to D* vertex (3-prong vertex fitter) - o2::vertexing::DCAFitterN<2> dfV0; // fitter for V0s (3-prong vertex fitter) + o2::vertexing::DCAFitterN<2> dfStrangeness; // fitter for V0s and cascades (2-prong vertex fitter) HistogramRegistry registry{"registry"}; std::shared_ptr hProcessedEvents; @@ -216,8 +216,8 @@ struct HfFilter { // Main struct for HF triggers helper.setPtRangeSoftPiSigmaC(ptCuts->get(0u, 4u), ptCuts->get(1u, 4u)); helper.setPtDeltaMassRangeSigmaC(cutsPtDeltaMassCharmReso->get(0u, 6u), cutsPtDeltaMassCharmReso->get(1u, 6u), cutsPtDeltaMassCharmReso->get(0u, 7u), cutsPtDeltaMassCharmReso->get(1u, 7u), cutsPtDeltaMassCharmReso->get(0u, 8u), cutsPtDeltaMassCharmReso->get(1u, 8u), cutsPtDeltaMassCharmReso->get(0u, 9u), cutsPtDeltaMassCharmReso->get(1u, 9u), cutsPtDeltaMassCharmReso->get(2u, 6u), cutsPtDeltaMassCharmReso->get(2u, 7u), cutsPtDeltaMassCharmReso->get(2u, 8u), cutsPtDeltaMassCharmReso->get(2u, 9u)); helper.setPtRangeSoftKaonXicResoToSigmaC(ptCuts->get(0u, 5u), ptCuts->get(1u, 5u)); - helper.setVtxConfiguration(dfV0, true); // (DCAFitterN, useAbsDCA) - dfV0.setMatCorrType(matCorr); + helper.setVtxConfiguration(dfStrangeness, true); // (DCAFitterN, useAbsDCA) + dfStrangeness.setMatCorrType(matCorr); if (activateSecVtxForB) { helper.setVtxConfiguration(df2, false); // (DCAFitterN, useAbsDCA) helper.setVtxConfiguration(df3, false); @@ -330,7 +330,7 @@ struct HfFilter { // Main struct for HF triggers using BigTracksMCPID = soa::Join; using BigTracksPID = soa::Join; - using TracksIUPID = soa::Join; + using TracksIUPID = soa::Join; using CollsWithEvSel = soa::Join; using Hf2ProngsWithMl = soa::Join; @@ -340,13 +340,13 @@ struct HfFilter { // Main struct for HF triggers Preslice v0sPerCollision = aod::v0data::collisionId; Preslice hf2ProngPerCollision = aod::track_association::collisionId; Preslice hf3ProngPerCollision = aod::track_association::collisionId; - Preslice cascPerCollision = aod::cascdata::collisionId; + Preslice cascPerCollision = aod::cascdata::collisionId; Preslice photonsPerCollision = aod::v0photonkf::collisionId; void process(CollsWithEvSel const& collisions, aod::BCsWithTimestamps const&, aod::V0s const& v0s, - aod::CascDatas const& cascades, + aod::Cascades const& cascades, Hf2ProngsWithMl const& cand2Prongs, Hf3ProngsWithMl const& cand3Prongs, aod::TrackAssoc const& trackIndices, @@ -385,6 +385,15 @@ struct HfFilter { // Main struct for HF triggers helper.setValuesBB(ccdbApi, bc, std::array{ccdbBBPion.value, ccdbBBAntiPion.value, ccdbBBKaon.value, ccdbBBAntiKaon.value, ccdbBBProton.value, ccdbBBAntiProton.value, ccdbBBProton.value, ccdbBBAntiProton.value}); // dummy for deuteron } + auto bz = o2::base::Propagator::Instance()->getNominalBz(); + dfStrangeness.setBz(bz); + if (activateSecVtxForB) { + df2.setBz(bz); + df3.setBz(bz); + dfB.setBz(bz); + dfBtoDstar.setBz(bz); + } + currentRun = bc.runNumber(); } @@ -708,7 +717,7 @@ struct HfFilter { // Main struct for HF triggers auto v0sThisCollision = v0s.sliceBy(v0sPerCollision, thisCollId); for (const auto& v0 : v0sThisCollision) { V0Cand v0Cand; - if (!buildV0(v0, tracksIU, collision, dfV0, std::vector{cand2Prong.prong0Id(), cand2Prong.prong1Id()}, v0Cand)) { + if (!helper.buildV0(v0, tracksIU, collision, dfStrangeness, std::vector{cand2Prong.prong0Id(), cand2Prong.prong1Id()}, v0Cand)) { continue; } auto selV0 = helper.isSelectedV0(v0Cand, activateQA, hV0Selected, hArmPod); @@ -1241,7 +1250,7 @@ struct HfFilter { // Main struct for HF triggers if ((!keepEvent[kV0Charm3P] && isGoodDPlus) || (!keepEvent[kSigmaC0K0] && (isGoodLcToPKPi || isGoodLcToPiKP))) { for (const auto& v0 : v0sThisCollision) { V0Cand v0Cand; - if (!buildV0(v0, tracksIU, collision, dfV0, std::vector{cand3Prong.prong0Id(), cand3Prong.prong1Id(), cand3Prong.prong2Id()}, v0Cand)) { + if (!helper.buildV0(v0, tracksIU, collision, dfStrangeness, std::vector{cand3Prong.prong0Id(), cand3Prong.prong1Id(), cand3Prong.prong2Id()}, v0Cand)) { continue; } auto selV0 = helper.isSelectedV0(v0Cand, activateQA, hV0Selected, hArmPod); @@ -1381,35 +1390,43 @@ struct HfFilter { // Main struct for HF triggers if (!keepEvent[kCharmBarToXiBach]) { auto cascThisColl = cascades.sliceBy(cascPerCollision, thisCollId); for (const auto& casc : cascThisColl) { - auto bachelorCasc = casc.bachelor_as(); - auto v0DauPos = casc.posTrack_as(); - auto v0DauNeg = casc.negTrack_as(); - if (!helper.isSelectedCascade(casc, std::array{bachelorCasc, v0DauPos, v0DauNeg}, collision)) { + CascCand cascCand; + if(!helper.buildCascade(casc, v0s, tracksIU, collision, dfStrangeness, {}, cascCand)) { continue; } + + if (!helper.isSelectedCascade(cascCand)) { + continue; + } + if (activateQA) { - hMassXi->Fill(casc.mXi()); + hMassXi->Fill(cascCand.mXi); } + auto bachelorCascId = casc.bachelorId(); + auto v0 = v0s.rawIteratorAt(casc.v0Id()); + auto v0DauPosId = v0.posTrackId(); + auto v0DauNegId = v0.negTrackId(); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, thisCollId); for (const auto& trackId : trackIdsThisCollision) { // start loop over tracks auto track = trackId.track_as(); // ask for opposite sign daughters (omegac daughters) - if (track.sign() * bachelorCasc.sign() >= 0) { + if (track.sign() * cascCand.sign > 0) { continue; } // check if track is one of the Xi daughters - if (track.globalIndex() == bachelorCasc.globalIndex() || track.globalIndex() == v0DauPos.globalIndex() || track.globalIndex() == v0DauNeg.globalIndex()) { + if (track.globalIndex() == bachelorCascId || track.globalIndex() == v0DauPosId || track.globalIndex() == v0DauNegId) { continue; } // propagate to PV gpu::gpustd::array dcaInfo; - std::array pVecCascade = {casc.px(), casc.py(), casc.pz()}; - auto trackParCasc = o2::track::TrackPar(std::array{casc.x(), casc.y(), casc.z()}, pVecCascade, bachelorCasc.sign(), true); + std::array pVecCascade = cascCand.mom; + auto trackParCasc = o2::track::TrackPar(cascCand.vtx, cascCand.mom, cascCand.sign, true); trackParCasc.setPID(o2::track::PID::XiMinus); o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCasc, 2.f, matCorr, &dcaInfo); getPxPyPz(trackParCasc, pVecCascade); diff --git a/EventFiltering/PWGHF/HFFilterHelpers.h b/EventFiltering/PWGHF/HFFilterHelpers.h index 5d1a2b1f49b..7999e66dcaf 100644 --- a/EventFiltering/PWGHF/HFFilterHelpers.h +++ b/EventFiltering/PWGHF/HFFilterHelpers.h @@ -148,6 +148,8 @@ enum HfVtxStage : uint8_t { // Helper struct to pass V0 informations struct V0Cand { std::array mom; + std::array vtx; + std::array cov; float etaPos; float etaNeg; float ptPos; @@ -156,9 +158,12 @@ struct V0Cand { float pinTpcNeg; float nClsFoundTpcPos; float nClsFoundTpcNeg; + float nClsCrossedRowsTpcPos; + float nClsCrossedRowsTpcNeg; + float crossedRowsOverFindableClsTpcPos; + float crossedRowsOverFindableClsTpcNeg; float signalTpcPos; float signalTpcNeg; - float nClsTpcNeg; float v0cosPA; float dcav0topv; float dcaV0daughters; @@ -174,10 +179,39 @@ struct V0Cand { float nSigmaPrTofPos; float nSigmaPrTpcNeg; float nSigmaPrTofNeg; + float nSigmaPiTpcPos; + float nSigmaPiTofPos; + float nSigmaPiTpcNeg; + float nSigmaPiTofNeg; bool hasTofPos; bool hasTofNeg; }; +// Helper struct to pass Cascade informations +struct CascCand { + std::array mom; + std::array vtx; + V0Cand v0; + float ptBach; + float etaBach; + float pinTpcBach; + float nClsFoundTpcBach; + float nClsCrossedRowsTpcBach; + float crossedRowsOverFindableClsTpcBach; + float signalTpcBach; + float pt; + float casccosPA; + float cascradius; + float dcaXYCascToPV; + float dcacascdaughters; + float mXi; + float mOmega; + float nSigmaPiTpcBach; + float nSigmaPiTofBach; + bool hasTofBach; + int sign; +}; + static const std::array charmParticleNames{"D0", "Dplus", "Ds", "Lc", "Xic"}; static const std::array beautyParticleNames{"Bplus", "B0toDStar", "B0", "Bs", "Lb", "Xib"}; static const std::array pdgCodesCharm{421, 411, 431, 4122, 4232}; @@ -499,8 +533,8 @@ class HfFilterHelper int8_t isSelectedV0(const V0& v0, const int& activateQA, H2 hV0Selected, std::array& hArmPod); template bool isSelectedPhoton(const Photon& photon, const std::array& dauTracks, const int& activateQA, H2 hV0Selected, std::array& hArmPod); - template - bool isSelectedCascade(const Casc& casc, const std::array& dauTracks, const Coll& collision); + template + bool isSelectedCascade(const Casc& casc); template int8_t isSelectedBachelorForCharmBaryon(const T& track, const T2& dca); template @@ -523,6 +557,10 @@ class HfFilterHelper int computeNumberOfCandidates(std::vector> indices); template int setVtxConfiguration(T1 vertexer, bool useAbsDCA); + template + bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, V0Cand& v0Cand); + template + bool buildCascade(Casc const& cascIndices, V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, CascCand& cascCand); // PID void setValuesBB(o2::ccdb::CcdbApi& ccdbApi, aod::BCsWithTimestamps::iterator const& bunchCrossing, const std::array& ccdbPaths); @@ -548,8 +586,6 @@ class HfFilterHelper int findBin(T1 const& binsPt, T2 value); template std::array alphaAndQtAP(std::array const& momPos, std::array const& momNeg); - template - bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, V0Cand& v0Cand); // selections std::vector mPtBinsTracks{}; // vector of pT bins for single track cuts @@ -1370,123 +1406,127 @@ inline bool HfFilterHelper::isSelectedPhoton(const Photon& photon, const std::ar /// Basic selection of cascade candidates /// \param casc is the cascade candidate -/// \param dauTracks is a 3-element array with bachelor, positive and negative V0 daughter tracks -/// \param collision is the collision /// \return true if cascade passes all cuts -template -inline bool HfFilterHelper::isSelectedCascade(const Casc& casc, const std::array& dauTracks, const Coll& collision) +template +inline bool HfFilterHelper::isSelectedCascade(const Casc& casc) { // Xi min pT - if (casc.pt() < mMinPtXi) { + if (casc.pt < mMinPtXi) { return false; } // eta of daughters - if (std::fabs(dauTracks[0].eta()) > 1. || std::fabs(dauTracks[1].eta()) > 1. || std::fabs(dauTracks[2].eta()) > 1.) { // cut all V0 daughters with |eta| > 1. + if (std::fabs(casc.v0.etaPos) > 1. || std::fabs(casc.v0.etaNeg) > 1. || std::fabs(casc.etaBach) > 1.) { // cut all V0 daughters with |eta| > 1. return false; } // V0 radius - if (casc.v0radius() < 1.2) { + if (casc.v0.v0radius < 1.2) { return false; } // cascade radius - if (casc.cascradius() < 0.6) { + if (casc.cascradius < 0.6) { return false; } // V0 cosp - if (casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()) < mCosPaLambdaFromXi) { + if (casc.v0.v0cosPA < mCosPaLambdaFromXi) { return false; } // cascade cosp - if (casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()) < mCosPaXi) { + if (casc.casccosPA < mCosPaXi) { return false; } // cascade DCAxy to PV - if (std::fabs(casc.dcaXYCascToPV()) > mMaxDcaXyXi) { + if (std::fabs(casc.dcaXYCascToPV) > mMaxDcaXyXi) { return false; } // Xi bachelor min pT - if (dauTracks[0].pt() < mMinPtXiBachelor) { + if (casc.ptBach < mMinPtXiBachelor) { return false; } // dau dca - if (std::fabs(casc.dcaV0daughters()) > 1.f || std::fabs(casc.dcacascdaughters()) > 1.f) { + if (std::fabs(casc.v0.dcaV0daughters) > 1.f || std::fabs(casc.dcacascdaughters) > 1.f) { return false; } // cascade mass - if (std::fabs(casc.mXi() - massXi) > mDeltaMassXi) { + if (std::fabs(casc.mXi - massXi) > mDeltaMassXi) { return false; } // V0 mass - if (std::fabs(casc.mLambda() - massLambda) > mDeltaMassLambdaFromXi) { + if (std::fabs(casc.v0.mLambda - massLambda) > mDeltaMassLambdaFromXi) { return false; } // PID - float nSigmaPrTpc[3] = {-999., dauTracks[1].tpcNSigmaPr(), dauTracks[2].tpcNSigmaPr()}; - float nSigmaPrTof[3] = {-999., dauTracks[1].tofNSigmaPr(), dauTracks[2].tofNSigmaPr()}; - float nSigmaPiTpc[3] = {dauTracks[0].tpcNSigmaPi(), dauTracks[1].tpcNSigmaPi(), dauTracks[2].tpcNSigmaPi()}; - float nSigmaPiTof[3] = {dauTracks[0].tofNSigmaPi(), dauTracks[1].tofNSigmaPi(), dauTracks[2].tofNSigmaPi()}; + float nSigmaPrTpc[3] = {-999, casc.v0.nSigmaPrTpcPos, casc.v0.nSigmaPrTpcNeg}; + float nSigmaPrTof[3] = {-999., casc.v0.nSigmaPrTofPos, casc.v0.nSigmaPrTofNeg}; + float nSigmaPiTpc[3] = {casc.nSigmaPiTpcBach, casc.v0.nSigmaPiTpcPos, casc.v0.nSigmaPiTpcNeg}; + float nSigmaPiTof[3] = {casc.nSigmaPiTofBach, casc.v0.nSigmaPiTofPos, casc.v0.nSigmaPiTofNeg}; + float pInTpc[3] = {casc.pinTpcBach, casc.v0.pinTpcPos, casc.v0.pinTpcNeg}; + float nClsTpc[3] = {casc.nClsFoundTpcBach, casc.v0.nClsFoundTpcPos, casc.v0.nClsFoundTpcNeg}; + float nCrossedRowsTpc[3] = {casc.nClsCrossedRowsTpcBach, casc.v0.nClsCrossedRowsTpcPos, casc.v0.nClsCrossedRowsTpcNeg}; + float crossedRowsOverFindableClsTpc[3] = {casc.crossedRowsOverFindableClsTpcBach, casc.v0.crossedRowsOverFindableClsTpcPos, casc.v0.crossedRowsOverFindableClsTpcNeg}; + float etaDaus[3] = {casc.etaBach, casc.v0.etaPos, casc.v0.etaNeg}; + float signalTpc[3] = {casc.signalTpcBach, casc.v0.signalTpcPos, casc.v0.signalTpcNeg}; if (mTpcPidCalibrationOption == 1) { for (int iDau{0}; iDau < 3; ++iDau) { - nSigmaPiTpc[iDau] = getTPCPostCalib(dauTracks[iDau], kPi); + nSigmaPiTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[2], kPi); if (iDau == 0) { continue; } - nSigmaPrTpc[iDau] = getTPCPostCalib(dauTracks[iDau], kPr); + nSigmaPrTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[2], kPr); } } else if (mTpcPidCalibrationOption == 2) { for (int iDau{0}; iDau < 3; ++iDau) { - nSigmaPiTpc[iDau] = getTPCSplineCalib(dauTracks[iDau], (dauTracks[iDau].sign() > 0) ? kPi : kAntiPi); + nSigmaPiTpc[iDau] = getTPCSplineCalib(pInTpc[iDau], signalTpc[iDau], (iDau == 0) ? kPi : kAntiPi); if (iDau == 0) { continue; } - nSigmaPrTpc[iDau] = getTPCSplineCalib(dauTracks[iDau], (dauTracks[iDau].sign() > 0) ? kPr : kAntiPr); + nSigmaPrTpc[iDau] = getTPCSplineCalib(pInTpc[iDau], signalTpc[iDau], (iDau == 0) ? kPr : kAntiPr); } } // PID to V0 tracks - if (dauTracks[0].sign() < 0) { // Xi- - if ((dauTracks[1].hasTPC() && std::fabs(nSigmaPrTpc[1]) > mMaxNsigmaXiDau) && (dauTracks[1].hasTOF() && std::fabs(nSigmaPrTof[1]) > mMaxNsigmaXiDau)) { + if (casc.sign < 0) { // Xi- + if (std::fabs(nSigmaPrTpc[1]) > mMaxNsigmaXiDau && (casc.v0.hasTofPos && std::fabs(nSigmaPrTof[1]) > mMaxNsigmaXiDau)) { return false; } - if ((dauTracks[2].hasTPC() && std::fabs(nSigmaPiTpc[2]) > mMaxNsigmaXiDau) && (dauTracks[2].hasTOF() && std::fabs(nSigmaPiTof[2]) > mMaxNsigmaXiDau)) { + if (std::fabs(nSigmaPiTpc[2]) > mMaxNsigmaXiDau && (casc.v0.hasTofNeg && std::fabs(nSigmaPiTof[2]) > mMaxNsigmaXiDau)) { return false; } - } else if (dauTracks[0].sign() > 0) { // Xi+ - if ((dauTracks[2].hasTPC() && std::fabs(nSigmaPrTpc[2]) > mMaxNsigmaXiDau) && (dauTracks[2].hasTOF() && std::fabs(nSigmaPrTof[2]) > mMaxNsigmaXiDau)) { + } else if (casc.sign > 0) { // Xi+ + if (std::fabs(nSigmaPrTpc[2]) > mMaxNsigmaXiDau && (casc.v0.hasTofNeg && std::fabs(nSigmaPrTof[2]) > mMaxNsigmaXiDau)) { return false; } - if ((dauTracks[1].hasTPC() && std::fabs(nSigmaPiTpc[1]) > mMaxNsigmaXiDau) && (dauTracks[1].hasTOF() && std::fabs(nSigmaPiTof[1]) > mMaxNsigmaXiDau)) { + if (std::fabs(nSigmaPiTpc[1]) > mMaxNsigmaXiDau && (casc.v0.hasTofPos && std::fabs(nSigmaPiTof[1]) > mMaxNsigmaXiDau)) { return false; } } // bachelor PID - if ((dauTracks[0].hasTPC() && std::fabs(nSigmaPiTpc[0]) > mMaxNsigmaXiDau) && (dauTracks[0].hasTOF() && std::fabs(nSigmaPiTof[0]) > mMaxNsigmaXiDau)) { + if (std::fabs(nSigmaPiTpc[0]) > mMaxNsigmaXiDau && (casc.hasTofBach && std::fabs(nSigmaPiTof[0]) > mMaxNsigmaXiDau)) { return false; } // additional track cuts - for (const auto& dauTrack : dauTracks) { + for (int iTrack{0}; iTrack < 3; ++iTrack) { // TPC clusters selections - if (dauTrack.tpcNClsFound() < 70) { // TODO: put me as a configurable please + if (nClsTpc[iTrack] < 70) { // TODO: put me as a configurable please return false; } - if (dauTrack.tpcNClsCrossedRows() < 70) { + if (nCrossedRowsTpc[iTrack] < 70) { return false; } - if (dauTrack.tpcCrossedRowsOverFindableCls() < 0.8) { + if (crossedRowsOverFindableClsTpc[iTrack] < 0.8) { return false; } } @@ -2058,7 +2098,7 @@ inline int HfFilterHelper::setVtxConfiguration(T1 vertexer, bool useAbsDCA) /// \param momPos momentum array of positive daughter /// \param momNeg momentum array of negative daughter template -inline std::array alphaAndQtAP(std::array const& momPos, std::array const& momNeg) +inline std::array HfFilterHelper::alphaAndQtAP(std::array const& momPos, std::array const& momNeg) { float momTotSq = RecoDecay::p2(momPos[0] + momNeg[0], momPos[1] + momNeg[1], momPos[2] + momNeg[2]); float momTot = std::sqrt(momTotSq); @@ -2078,7 +2118,7 @@ inline std::array alphaAndQtAP(std::array const& momPos, std::array< /// \param dcaFitter DCA fitter to be used /// \param vetoedTrackIds vector with forbidden track indices, if any template -inline bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, V0Cand& v0Cand) +inline bool HfFilterHelper::buildV0(V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, V0Cand& v0Cand) { auto trackPos = tracks.rawIteratorAt(v0Indices.posTrackId()); auto trackNeg = tracks.rawIteratorAt(v0Indices.negTrackId()); @@ -2134,6 +2174,10 @@ inline bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2: v0Cand.pinTpcNeg = trackNeg.tpcInnerParam(); v0Cand.nClsFoundTpcPos = trackPos.tpcNClsFound(); v0Cand.nClsFoundTpcNeg = trackNeg.tpcNClsFound(); + v0Cand.nClsCrossedRowsTpcPos = trackPos.tpcNClsCrossedRows(); + v0Cand.nClsCrossedRowsTpcNeg = trackNeg.tpcNClsCrossedRows(); + v0Cand.crossedRowsOverFindableClsTpcPos = trackPos.tpcCrossedRowsOverFindableCls(); + v0Cand.crossedRowsOverFindableClsTpcNeg = trackNeg.tpcCrossedRowsOverFindableCls(); v0Cand.signalTpcPos = trackPos.tpcSignal(); v0Cand.signalTpcNeg = trackNeg.tpcSignal(); v0Cand.etaPos = RecoDecay::eta(momPos); @@ -2141,6 +2185,25 @@ inline bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2: v0Cand.dcaV0daughters = std::sqrt(dcaFitter.getChi2AtPCACandidate()); const auto& vtx = dcaFitter.getPCACandidate(); + for (int iCoord{0}; iCoord < 3; ++iCoord) { + v0Cand.vtx[iCoord] = vtx[iCoord]; + } + auto covVtxV = dcaFitter.calcPCACovMatrix(0); + v0Cand.cov = {}; + v0Cand.cov[0] = covVtxV(0, 0); + v0Cand.cov[1] = covVtxV(1, 0); + v0Cand.cov[2] = covVtxV(1, 1); + v0Cand.cov[3] = covVtxV(2, 0); + v0Cand.cov[4] = covVtxV(2, 1); + v0Cand.cov[5] = covVtxV(2, 2); + std::array covTpositive = {0.}; + std::array covTnegative = {0.}; + trackPosProp.getCovXYZPxPyPzGlo(covTpositive); + trackNegProp.getCovXYZPxPyPzGlo(covTnegative); + constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component + for (int iCoord{0}; iCoord < 6; ++iCoord) { + v0Cand.cov[MomInd[iCoord]] = covTpositive[MomInd[iCoord]] + covTnegative[MomInd[iCoord]]; + } v0Cand.v0radius = std::hypot(vtx[0], vtx[1]); v0Cand.v0cosPA = RecoDecay::cpa(primVtx, vtx, v0Cand.mom); @@ -2165,6 +2228,108 @@ inline bool buildV0(V const& v0Indices, T const& tracks, C const& collision, o2: v0Cand.nSigmaPrTofPos = trackPos.tofNSigmaPr(); v0Cand.nSigmaPrTpcNeg = trackNeg.tpcNSigmaPr(); v0Cand.nSigmaPrTofNeg = trackNeg.tofNSigmaPr(); + v0Cand.nSigmaPiTpcPos = trackPos.tpcNSigmaPi(); + v0Cand.nSigmaPiTofPos = trackPos.tofNSigmaPi(); + v0Cand.nSigmaPiTpcNeg = trackNeg.tpcNSigmaPi(); + v0Cand.nSigmaPiTofNeg = trackNeg.tofNSigmaPi(); + + return true; +} + +/// build cascade candidate from table with track indices +/// \param cascIndices cascade candidate from AO2D table (track indices) +/// \param v0Indices V0 candidate from AO2D table (track indices) +/// \param tracks track table +/// \param collision collision +/// \param dcaFitter DCA fitter to be used +/// \param vetoedTrackIds vector with forbidden track indices, if any +template +inline bool HfFilterHelper::buildCascade(Casc const& cascIndices, V const& v0Indices, T const& tracks, C const& collision, o2::vertexing::DCAFitterN<2>& dcaFitter, const std::vector& vetoedTrackIds, CascCand& cascCand) +{ + auto v0 = v0Indices.rawIteratorAt(cascIndices.v0Id()); + auto trackBachelor = tracks.rawIteratorAt(cascIndices.bachelorId()); + + // minimal track cuts + if (!trackBachelor.hasTPC()) { + return false; + } + if (!(trackBachelor.trackType() & o2::aod::track::TPCrefit)) { + return false; + } + if (trackBachelor.tpcNClsCrossedRows() < 50) { + return false; + } + + if (std::find(vetoedTrackIds.begin(), vetoedTrackIds.end(), trackBachelor.globalIndex()) != vetoedTrackIds.end()) { + return false; + } + + gpu::gpustd::array dcaInfoBach; + auto bachTrackParCov = getTrackParCov(trackBachelor); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, bachTrackParCov, 2.f, dcaFitter.getMatCorrType(), &dcaInfoBach); + + // first we build V0 candidate + V0Cand v0Cand; + if (!buildV0(v0, tracks, collision, dcaFitter, vetoedTrackIds, v0Cand)) { + return false; + } + + // Set up covariance matrices (should in fact be optional) + auto v0TrackParCov = o2::track::TrackParCov({v0Cand.vtx[0], v0Cand.vtx[1], v0Cand.vtx[2]}, {v0Cand.mom[0], v0Cand.mom[1], v0Cand.mom[2]}, v0Cand.cov, 0, true); + v0TrackParCov.setAbsCharge(0); + v0TrackParCov.setPID(o2::track::PID::Lambda); + + int nCand = 0; + try { + nCand = dcaFitter.process(v0TrackParCov, bachTrackParCov); + } catch (...) { + LOG(error) << "Exception caught in DCA fitter process call!"; + return false; + } + if (nCand == 0) { + return false; + } + + // compute candidate momentum from tracks propagated to decay vertex + auto& trackV0Prop = dcaFitter.getTrack(0); + auto& trackBachProp = dcaFitter.getTrack(1); + std::array momV0{}, momBach{}; + trackV0Prop.getPxPyPzGlo(momV0); + trackBachProp.getPxPyPzGlo(momBach); + cascCand.mom = RecoDecay::pVec(momV0, momBach); + cascCand.sign = trackBachelor.sign(); + + cascCand.v0 = v0Cand; + cascCand.ptBach = RecoDecay::pt(momBach); + cascCand.etaBach = RecoDecay::eta(momBach); + cascCand.pinTpcBach = trackBachelor.tpcInnerParam(); + cascCand.nClsFoundTpcBach = trackBachelor.tpcNClsFound(); + cascCand.nClsCrossedRowsTpcBach = trackBachelor.tpcNClsCrossedRows(); + cascCand.crossedRowsOverFindableClsTpcBach = trackBachelor.tpcCrossedRowsOverFindableCls(); + cascCand.signalTpcBach = trackBachelor.tpcSignal(); + cascCand.pt = RecoDecay::pt(cascCand.mom); + + std::array primVtx = {collision.posX(), collision.posY(), collision.posZ()}; + const auto& vtx = dcaFitter.getPCACandidate(); + for (int iCoord{0}; iCoord < 3; ++iCoord) { + cascCand.vtx[iCoord] = vtx[iCoord]; + } + cascCand.cascradius = std::hypot(vtx[0], vtx[1]);; + cascCand.casccosPA = RecoDecay::cpa(primVtx, vtx, cascCand.mom); + + auto trackParCasc = dcaFitter.createParentTrackParCov(); + trackParCasc.setAbsCharge(1); + trackParCasc.setPID(o2::track::PID::XiMinus); + gpu::gpustd::array dcaInfoCasc; + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, trackParCasc, 2.f, dcaFitter.getMatCorrType(), &dcaInfoCasc); + cascCand.dcaXYCascToPV = dcaInfoCasc[0]; + cascCand.dcacascdaughters = std::sqrt(dcaFitter.getChi2AtPCACandidate()); + cascCand.mXi = RecoDecay::m(std::array{momBach, momV0}, std::array{massPi, massLambda}); + cascCand.mOmega = RecoDecay::m(std::array{momBach, momV0}, std::array{massKa, massLambda});; + + cascCand.hasTofBach = trackBachelor.hasTOF(); + cascCand.nSigmaPiTpcBach = trackBachelor.tpcNSigmaPi(); + cascCand.nSigmaPiTofBach = trackBachelor.tofNSigmaPi(); return true; } From 20f94165cdb43d86fc1ae81041845fa21bfde429 Mon Sep 17 00:00:00 2001 From: Fabrizio Grosa Date: Wed, 1 Jan 2025 18:55:31 +0100 Subject: [PATCH 3/4] Fix typo --- EventFiltering/PWGHF/HFFilterHelpers.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/EventFiltering/PWGHF/HFFilterHelpers.h b/EventFiltering/PWGHF/HFFilterHelpers.h index 7999e66dcaf..34099d5d0e5 100644 --- a/EventFiltering/PWGHF/HFFilterHelpers.h +++ b/EventFiltering/PWGHF/HFFilterHelpers.h @@ -1320,7 +1320,7 @@ inline int8_t HfFilterHelper::isSelectedV0(const V0& v0, const int& activateQA, float signalTpc[2] = {v0.signalTpcPos, v0.signalTpcNeg}; if (mTpcPidCalibrationOption == 1) { for (int iDau{0}; iDau < 2; ++iDau) { - nSigmaPrTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[2], kPr); + nSigmaPrTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[iDau], kPr); } } else if (mTpcPidCalibrationOption == 2) { for (int iDau{0}; iDau < 2; ++iDau) { @@ -1479,11 +1479,11 @@ inline bool HfFilterHelper::isSelectedCascade(const Casc& casc) float signalTpc[3] = {casc.signalTpcBach, casc.v0.signalTpcPos, casc.v0.signalTpcNeg}; if (mTpcPidCalibrationOption == 1) { for (int iDau{0}; iDau < 3; ++iDau) { - nSigmaPiTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[2], kPi); + nSigmaPiTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[iDau], kPi); if (iDau == 0) { continue; } - nSigmaPrTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[2], kPr); + nSigmaPrTpc[iDau] = getTPCPostCalib(pInTpc[iDau], nClsTpc[iDau], etaDaus[iDau], nSigmaPrTpc[iDau], kPr); } } else if (mTpcPidCalibrationOption == 2) { for (int iDau{0}; iDau < 3; ++iDau) { From fb04b0f20c774e0eb93a10362f4c874bfe32a336 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 1 Jan 2025 20:19:05 +0000 Subject: [PATCH 4/4] Please consider the following formatting changes --- EventFiltering/PWGHF/HFFilter.cxx | 2 +- EventFiltering/PWGHF/HFFilterHelpers.h | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/EventFiltering/PWGHF/HFFilter.cxx b/EventFiltering/PWGHF/HFFilter.cxx index 7436e21da99..1062635ce39 100644 --- a/EventFiltering/PWGHF/HFFilter.cxx +++ b/EventFiltering/PWGHF/HFFilter.cxx @@ -1392,7 +1392,7 @@ struct HfFilter { // Main struct for HF triggers for (const auto& casc : cascThisColl) { CascCand cascCand; - if(!helper.buildCascade(casc, v0s, tracksIU, collision, dfStrangeness, {}, cascCand)) { + if (!helper.buildCascade(casc, v0s, tracksIU, collision, dfStrangeness, {}, cascCand)) { continue; } diff --git a/EventFiltering/PWGHF/HFFilterHelpers.h b/EventFiltering/PWGHF/HFFilterHelpers.h index 34099d5d0e5..2a41c163542 100644 --- a/EventFiltering/PWGHF/HFFilterHelpers.h +++ b/EventFiltering/PWGHF/HFFilterHelpers.h @@ -2314,7 +2314,8 @@ inline bool HfFilterHelper::buildCascade(Casc const& cascIndices, V const& v0Ind for (int iCoord{0}; iCoord < 3; ++iCoord) { cascCand.vtx[iCoord] = vtx[iCoord]; } - cascCand.cascradius = std::hypot(vtx[0], vtx[1]);; + cascCand.cascradius = std::hypot(vtx[0], vtx[1]); + ; cascCand.casccosPA = RecoDecay::cpa(primVtx, vtx, cascCand.mom); auto trackParCasc = dcaFitter.createParentTrackParCov(); @@ -2325,7 +2326,8 @@ inline bool HfFilterHelper::buildCascade(Casc const& cascIndices, V const& v0Ind cascCand.dcaXYCascToPV = dcaInfoCasc[0]; cascCand.dcacascdaughters = std::sqrt(dcaFitter.getChi2AtPCACandidate()); cascCand.mXi = RecoDecay::m(std::array{momBach, momV0}, std::array{massPi, massLambda}); - cascCand.mOmega = RecoDecay::m(std::array{momBach, momV0}, std::array{massKa, massLambda});; + cascCand.mOmega = RecoDecay::m(std::array{momBach, momV0}, std::array{massKa, massLambda}); + ; cascCand.hasTofBach = trackBachelor.hasTOF(); cascCand.nSigmaPiTpcBach = trackBachelor.tpcNSigmaPi();