From c6a933b048f12d45c072d18d8ab713f36ab2a33f Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Mon, 20 Nov 2023 09:35:22 +0100 Subject: [PATCH] SVERTEXER: V0 tuning Signed-off-by: Felix Schlepper --- .../StrangenessTrackingConfigParam.h | 3 +- .../src/StrangenessTracker.cxx | 11 +- .../include/DetectorsVertexing/SVertexer.h | 3 +- .../DetectorsVertexing/SVertexerParams.h | 51 +++++---- Detectors/Vertexing/src/SVertexer.cxx | 103 +++++++++++++----- 5 files changed, 118 insertions(+), 53 deletions(-) diff --git a/Detectors/Vertexing/StrangenessTracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h b/Detectors/Vertexing/StrangenessTracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h index 3b0226dbcdba7..80c7a64c8f7ff 100644 --- a/Detectors/Vertexing/StrangenessTracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h +++ b/Detectors/Vertexing/StrangenessTracking/include/StrangenessTracking/StrangenessTrackingConfigParam.h @@ -32,10 +32,11 @@ struct StrangenessTrackingParamConfig : public o2::conf::ConfigurableParamHelper float mMinMotherClus = 3.; // minimum number of cluster to be attached to the mother float mMaxChi2 = 50; // Maximum matching chi2 bool mVertexMatching = true; // Flag to enable/disable vertex matching + bool mSkipTPC = true; // Flag to enable/disable TPC only tracks O2ParamDef(StrangenessTrackingParamConfig, "strtracker"); }; } // namespace strangeness_tracking } // namespace o2 -#endif \ No newline at end of file +#endif diff --git a/Detectors/Vertexing/StrangenessTracking/src/StrangenessTracker.cxx b/Detectors/Vertexing/StrangenessTracking/src/StrangenessTracker.cxx index e1ea987e6dbc8..56daede8171ec 100644 --- a/Detectors/Vertexing/StrangenessTracking/src/StrangenessTracker.cxx +++ b/Detectors/Vertexing/StrangenessTracking/src/StrangenessTracker.cxx @@ -112,6 +112,9 @@ void StrangenessTracker::prepareITStracks() // sort tracks by eta and phi and se void StrangenessTracker::processV0(int iv0, const V0& v0, const V0Index& v0Idx, int iThread) { + if (mStrParams->mSkipTPC && (v0Idx.getProngID(kV0DauNeg).getSource() == GIndex::TPC || v0Idx.getProngID(kV0DauPos).getSource() == GIndex::TPC)) { + return; + } ClusAttachments structClus; auto& daughterTracks = mDaughterTracks[iThread]; daughterTracks.resize(2); // resize to 2 prongs: first positive second negative @@ -205,10 +208,10 @@ void StrangenessTracker::processCascade(int iCasc, const Cascade& casc, const Ca } decayVtxTrackClone.getPxPyPzGlo(strangeTrack.mDecayMom); std::array momV0, mombach; - mFitter3Body[iThread].getTrack(0).getPxPyPzGlo(momV0); // V0 momentum at decay vertex - mFitter3Body[iThread].getTrack(1).getPxPyPzGlo(mombach); // bachelor momentum at decay vertex - strangeTrack.mMasses[0] = calcMotherMass(momV0, mombach, PID::Lambda, PID::Pion); // Xi invariant mass at decay vertex - strangeTrack.mMasses[1] = calcMotherMass(momV0, mombach, PID::Lambda, PID::Kaon); // Omega invariant mass at decay vertex + mFitter3Body[iThread].getTrack(0).getPxPyPzGlo(momV0); // V0 momentum at decay vertex + mFitter3Body[iThread].getTrack(1).getPxPyPzGlo(mombach); // bachelor momentum at decay vertex + strangeTrack.mMasses[0] = calcMotherMass(momV0, mombach, PID::Lambda, PID::Pion); // Xi invariant mass at decay vertex + strangeTrack.mMasses[1] = calcMotherMass(momV0, mombach, PID::Lambda, PID::Kaon); // Omega invariant mass at decay vertex LOG(debug) << "ITS Track matched with a Cascade decay topology ...."; LOG(debug) << "Number of ITS track clusters attached: " << itsTrack.getNumberOfClusters(); diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h index 42f67378e2e04..bde5079b20b0b 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h @@ -112,6 +112,7 @@ class SVertexer void setEnable3BodyDecays(bool v) { mEnable3BodyDecays = v; } void init(); void process(const o2::globaltracking::RecoContainer& recoTracks, o2::framework::ProcessingContext& pc); + void produceOutput(o2::framework::ProcessingContext& pc); int getNV0s() const { return mNV0s; } int getNCascades() const { return mNCascades; } int getN3Bodies() const { return mN3Bodies; } @@ -152,7 +153,7 @@ class SVertexer void setupThreads(); void buildT2V(const o2::globaltracking::RecoContainer& recoTracks); void updateTimeDependentParams(); - bool acceptTrack(GIndex gid, const o2::track::TrackParCov& trc) const; + bool acceptTrack(const GIndex& gid, const o2::track::TrackParCov& trc) const; bool processTPCTrack(const o2::tpc::TrackTPC& trTPC, GIndex gid, int vtxid); float correctTPCTrack(TrackCand& trc, const o2::tpc::TrackTPC& tTPC, float tmus, float tmusErr) const; diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h index 47c56a400a289..433e370f831fe 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h @@ -19,6 +19,7 @@ #include "CommonUtils/ConfigurableParam.h" #include "CommonUtils/ConfigurableParamHelper.h" #include "DetectorsVertexing/SVertexHypothesis.h" +#include "CommonConstants/GeomConstants.h" #include "DetectorsBase/Propagator.h" namespace o2 @@ -40,6 +41,7 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper minMomTPCdEdx are always accepted - float minMomTPCdEdx = 0.8; // minimum p for tracks with dEdx > mMinTPCdEdx to be accepted + // Cuts for TPC only tracks + bool mExcludeTPCtracks = false; // don't loop over TPC tracks if true (if loaded, dEdx info is used instead) + float mTPCTrackMaxX = constants::geom::XTPCInnerRef; // don't use TPC standalone tracks with X exceeding this; + float mTPCTrack2Beam = 21.f; // Straight line for TPC track back to beamline + float minTPCdEdx = 250; // starting from this dEdx value, tracks with p > minMomTPCdEdx are always accepted + float minMomTPCdEdx = 0.8; // minimum p for tracks with dEdx > mMinTPCdEdx to be accepted uint8_t mITSSAminNclu = 6; // global requirement of at least this many ITS clusters if no TPC info present (N.B.: affects all secondary vertexing) uint8_t mITSSAminNcluCascades = 6; // require at least this many ITS clusters if no TPC info present for cascade finding. @@ -108,10 +114,15 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelperloadData(recoData); mStrTracker->prepareITStracks(); } -#ifdef WITH_OPENMP - int dynGrp = std::min(4, std::max(1, mNThreads / 2)); -#pragma omp parallel for schedule(dynamic, dynGrp) num_threads(mNThreads) -#endif + + int iThread = 0; +#pragma omp parallel for schedule(dynamic, std::min(4, std::max(1, mNThreads / 2))) \ + num_threads(mNThreads) if (mNThreads > 1) private(iThread) for (int itp = 0; itp < ntrP; itp++) { auto& seedP = mTracksPool[POS][itp]; - int firstN = mVtxFirstTrack[NEG][seedP.vBracket.getMin()]; + const int firstN = mVtxFirstTrack[NEG][seedP.vBracket.getMin()]; if (firstN < 0) { LOG(debug) << "No partner is found for pos.track " << itp << " out of " << ntrP; continue; @@ -71,14 +71,18 @@ void SVertexer::process(const o2::globaltracking::RecoContainer& recoData, o2::f continue; } #ifdef WITH_OPENMP - int iThread = omp_get_thread_num(); -#else - int iThread = 0; + iThread = omp_get_thread_num(); #endif checkV0(seedP, seedN, itp, itn, iThread); } } + produceOutput(pc); +} + +//__________________________________________________________________ +void SVertexer::produceOutput(o2::framework::ProcessingContext& pc) +{ // sort V0s and Cascades in vertex id struct vid { int thrID; @@ -108,8 +112,8 @@ void SVertexer::process(const o2::globaltracking::RecoContainer& recoData, o2::f std::sort(v0SortID.begin(), v0SortID.end(), [](const vid& a, const vid& b) { return a.vtxID < b.vtxID; }); std::sort(cascSortID.begin(), cascSortID.end(), [](const vid& a, const vid& b) { return a.vtxID < b.vtxID; }); std::sort(nbodySortID.begin(), nbodySortID.end(), [](const vid& a, const vid& b) { return a.vtxID < b.vtxID; }); - // sorted V0s + // dpl outpu auto& v0sIdx = pc.outputs().make>(o2f::Output{"GLO", "V0S_IDX", 0, o2f::Lifetime::Timeframe}); auto& cascsIdx = pc.outputs().make>(o2f::Output{"GLO", "CASCS_IDX", 0, o2f::Lifetime::Timeframe}); auto& body3Idx = pc.outputs().make>(o2f::Output{"GLO", "DECAYS3BODY_IDX", 0, o2f::Lifetime::Timeframe}); @@ -236,7 +240,7 @@ void SVertexer::process(const o2::globaltracking::RecoContainer& recoData, o2::f strTrMCLableOut.swap(mcLabsOut); } } - // + for (int ith = 0; ith < mNThreads; ith++) { // clean unneeded s.vertices mV0sTmp[ith].clear(); mCascadesTmp[ith].clear(); @@ -344,6 +348,7 @@ void SVertexer::setupThreads() fitter.setMinParamChange(mSVParams->minParamChange); fitter.setMinRelChi2Change(mSVParams->minRelChi2Change); fitter.setMaxDZIni(mSVParams->maxDZIni); + fitter.setMaxDXYIni(mSVParams->maxDXYIni); fitter.setMaxChi2(mSVParams->maxChi2); fitter.setMatCorrType(o2::base::Propagator::MatCorrType(mSVParams->matCorr)); fitter.setUsePropagator(mSVParams->usePropagator); @@ -363,6 +368,7 @@ void SVertexer::setupThreads() fitter.setMinParamChange(mSVParams->minParamChange); fitter.setMinRelChi2Change(mSVParams->minRelChi2Change); fitter.setMaxDZIni(mSVParams->maxDZIni); + fitter.setMaxDXYIni(mSVParams->maxDXYIni); fitter.setMaxChi2(mSVParams->maxChi2); fitter.setMatCorrType(o2::base::Propagator::MatCorrType(mSVParams->matCorr)); fitter.setUsePropagator(mSVParams->usePropagator); @@ -383,6 +389,7 @@ void SVertexer::setupThreads() fitter.setMinParamChange(mSVParams->minParamChange); fitter.setMinRelChi2Change(mSVParams->minRelChi2Change); fitter.setMaxDZIni(mSVParams->maxDZIni); + fitter.setMaxDXYIni(mSVParams->maxDXYIni); fitter.setMaxChi2(mSVParams->maxChi2); fitter.setMatCorrType(o2::base::Propagator::MatCorrType(mSVParams->matCorr)); fitter.setUsePropagator(mSVParams->usePropagator); @@ -394,7 +401,7 @@ void SVertexer::setupThreads() } //__________________________________________________________________ -bool SVertexer::acceptTrack(GIndex gid, const o2::track::TrackParCov& trc) const +bool SVertexer::acceptTrack(const GIndex& gid, const o2::track::TrackParCov& trc) const { if (gid.isPVContributor() && mSVParams->maxPVContributors < 1) { return false; @@ -551,15 +558,24 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a //__________________________________________________________________ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, int iN, int ithread) { + // Fast rough cuts on pairs before feeding to DCAFitter, tracks are not in the same Frame or at same X + bool usesTPCOnly = (seedP.hasTPC && seedP.nITSclu < 1) || (seedN.hasTPC && seedN.nITSclu < 1); + if (auto deta = std::abs(seedP.getEta() - seedN.getEta()), dphi = std::abs(seedP.getPhi() - seedN.getPhi()); + usesTPCOnly && (deta > mSVParams->maxV0EtaAbsDiff || dphi > mSVParams->maxV0PhiAbsDiff)) { + LOG(debug) << "RejEtaPhi " << deta << " " << dphi; + return false; + } + // auto& fitterV0 = mFitterV0[ithread]; int nCand = fitterV0.process(seedP, seedN); if (nCand == 0) { // discard this pair + LOG(debug) << "RejDCAFitter"; return false; } const auto& v0XYZ = fitterV0.getPCACandidate(); // validate V0 radial position // check closeness to the beam-line - float dxv0 = v0XYZ[0] - mMeanVertex.getX(), dyv0 = v0XYZ[1] - mMeanVertex.getY(), r2v0 = dxv0 * dxv0 + dyv0 * dyv0; + float dxv0 = static_cast(v0XYZ[0]) - mMeanVertex.getX(), dyv0 = v0XYZ[1] - mMeanVertex.getY(), r2v0 = dxv0 * dxv0 + dyv0 * dyv0; if (r2v0 < mMinR2ToMeanVertex) { return false; } @@ -569,13 +585,14 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, LOG(debug) << "RejCausality " << drv0P << " " << drv0N; return false; } - const int cand = 0; - if (!fitterV0.isPropagateTracksToVertexDone(cand) && !fitterV0.propagateTracksToVertex(cand)) { + + if (!fitterV0.isPropagateTracksToVertexDone() && !fitterV0.propagateTracksToVertex()) { + LOG(debug) << "RejProp failed"; return false; } - auto& trPProp = fitterV0.getTrack(0, cand); - auto& trNProp = fitterV0.getTrack(1, cand); - std::array pP, pN; + const auto& trPProp = fitterV0.getTrack(0); + const auto& trNProp = fitterV0.getTrack(1); + std::array pP{}, pN{}; trPProp.getPxPyPzGlo(pP); trNProp.getPxPyPzGlo(pN); // estimate DCA of neutral V0 track to beamline: straight line with parametric equation @@ -598,14 +615,13 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, bool goodHyp = false; std::array hypCheckStatus{}; - for (int ipid = 0; ipid < NHypV0; ipid++) { + for (int ipid = 0; ipid < NHypV0 && mSVParams->checkV0Hypothesis; ipid++) { if (mV0Hyps[ipid].check(p2Pos, p2Neg, p2V0, ptV0)) { goodHyp = hypCheckStatus[ipid] = true; } } // check tight lambda mass only bool goodLamForCascade = false, goodALamForCascade = false; - bool usesTPCOnly = (seedP.hasTPC && seedP.nITSclu < 1) || (seedN.hasTPC && seedN.nITSclu < 1); bool usesShortITSOnly = (!seedP.hasTPC && seedP.nITSclu < mSVParams->mITSSAminNcluCascades) || (!seedN.hasTPC && seedN.nITSclu < mSVParams->mITSSAminNcluCascades); if (ptV0 > mSVParams->minPtV0FromCascade && (!mSVParams->mSkipTPCOnlyCascade || !usesTPCOnly) && !usesShortITSOnly) { if (mV0Hyps[Lambda].checkTight(p2Pos, p2Neg, p2V0, ptV0) && (!mSVParams->mRequireTPCforCascBaryons || seedP.hasTPC) && seedP.compatibleProton) { @@ -616,6 +632,31 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, } } + // apply loose Armenteros-Podolanski cut to photons and K0 + if (mSVParams->checkV0Hypothesis && mSVParams->checkV0AP && + (hypCheckStatus[HypV0::Photon] || hypCheckStatus[HypV0::K0] || goodLamForCascade || goodALamForCascade)) { + float pV0Abs = std::sqrt(p2V0); + float pPos = std::sqrt(p2Pos); + float p1 = (pV0[0] * pP[0] + pV0[1] * pP[1] + pV0[2] * pP[2]) / + (pV0Abs * pPos); + float p2 = (pV0[0] * pN[0] + pV0[1] * pN[1] + pV0[2] * pN[2]) / + (pV0Abs * seedN.getP()); + float pL1 = p1 * seedP.getP(); + float pL2 = p2 * seedN.getP(); + float alpha = (pL1 - pL2) / (pL1 + pL2); + float pT1 = pPos * std::sqrt(1 - p1); + if (hypCheckStatus[HypV0::Photon] && pT1 > mSVParams->pidCutsPhotonAP_qT) { + hypCheckStatus[HypV0::Photon] = false; + } + if (hypCheckStatus[HypV0::K0] && pT1 < mSVParams->pidCutsK0AP_qT) { + hypCheckStatus[HypV0::K0] = false; + } + if ((goodLamForCascade || goodALamForCascade) && std::abs(alpha) > mSVParams->pidCutsLambdaAP_a && pT1 > mSVParams->pidCutsLambdaAP_qT) { + goodLamForCascade = false; + goodALamForCascade = false; + } + } + // apply mass selections for 3-body decay bool good3bodyV0Hyp = false; for (int ipid = 2; ipid < 4; ipid++) { @@ -629,7 +670,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, // we want to reconstruct the 3 body decay of hypernuclei starting from the V0 of a proton and a pion (e.g. H3L->d + (p + pi-), or He4L->He3 + (p + pi-))) bool checkFor3BodyDecays = mEnable3BodyDecays && (!mSVParams->checkV0Hypothesis || good3bodyV0Hyp) && (pt2V0 > 0.5) && (!mSVParams->mSkipTPCOnly3Body || !usesTPCOnly); bool rejectAfter3BodyCheck = false; // To reject v0s which can be 3-body decay candidates but not cascade or v0 - bool checkForCascade = mEnableCascades && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (goodLamForCascade || goodALamForCascade)); + bool checkForCascade = mEnableCascades && (!mSVParams->mSkipTPCOnlyCascade || !usesTPCOnly) && r2v0 < mMaxR2ToMeanVertexCascV0 && (!mSVParams->checkV0Hypothesis || (goodLamForCascade || goodALamForCascade) && !hypCheckStatus[HypV0::Photon]); bool rejectIfNotCascade = false; if (!goodHyp && mSVParams->checkV0Hypothesis) { @@ -680,7 +721,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, for (int iv = vlist.getMin(); iv <= vlist.getMax(); iv++) { const auto& pv = mPVertices[iv]; - const auto v0XYZ = fitterV0.getPCACandidatePos(cand); + const auto v0XYZ = fitterV0.getPCACandidatePos(); // check cos of pointing angle float dx = v0XYZ[0] - pv.getX(), dy = v0XYZ[1] - pv.getY(), dz = v0XYZ[2] - pv.getZ(), prodXYZv0 = dx * pV0[0] + dy * pV0[1] + dz * pV0[2]; float cosPA = prodXYZv0 / std::sqrt((dx * dx + dy * dy + dz * dz) * p2V0); @@ -689,9 +730,9 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, continue; } if (!candFound) { - new (&v0new) V0(v0XYZ, pV0, fitterV0.calcPCACovMatrixFlat(cand), trPProp, trNProp); + new (&v0new) V0(v0XYZ, pV0, fitterV0.calcPCACovMatrixFlat(), trPProp, trNProp); new (&v0Idxnew) V0Index(-1, seedP.gid, seedN.gid); - v0new.setDCA(fitterV0.getChi2AtPCACandidate(cand)); + v0new.setDCA(fitterV0.getChi2AtPCACandidate()); candFound = true; } v0new.setCosPA(cosPA); @@ -756,7 +797,6 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, //__________________________________________________________________ int SVertexer::checkCascades(const V0Index& v0Idx, const V0& v0, float rv0, std::array pV0, float p2V0, int avoidTrackID, int posneg, VBracket v0vlist, int ithread) { - // check last added V0 for belonging to cascade auto& fitterCasc = mFitterCasc[ithread]; auto& tracks = mTracksPool[posneg]; @@ -775,7 +815,9 @@ int SVertexer::checkCascades(const V0Index& v0Idx, const V0& v0, float rv0, std: continue; // skip the track used by V0 } auto& bach = tracks[it]; - + if (mSVParams->mSkipTPCOnlyCascade && bach.hasTPC) { + continue; // reject TPC-only bachelors + } if (!bach.hasTPC && bach.nITSclu < mSVParams->mITSSAminNcluCascades) { continue; // reject short ITS-only } @@ -944,6 +986,9 @@ int SVertexer::check3bodyDecays(const V0Index& v0Idx, const V0& v0, float rv0, s continue; // skip the track used by V0 } auto& bach = tracks[it]; + if (mSVParams->mSkipTPCOnly3Body && bach.hasTPC) { + continue; // reject TPC-only bachelors + } if (bach.vBracket > v0vlist.getMax()) { LOG(debug) << "Skipping"; break; // all other bachelor candidates will be also not compatible with this PV @@ -1170,10 +1215,11 @@ bool SVertexer::processTPCTrack(const o2::tpc::TrackTPC& trTPC, GIndex gid, int const auto& vtx = mPVertices[vtxid]; auto twe = vtx.getTimeStamp(); int posneg = trTPC.getSign() < 0 ? 1 : 0; - auto& trLoc = mTracksPool[posneg].emplace_back(TrackCand{trTPC, gid, {vtxid, vtxid}, 0.}); + auto& trLoc = mTracksPool[posneg].emplace_back(TrackCand{trTPC, gid, {vtxid, vtxid}, 0., true}); auto err = correctTPCTrack(trLoc, trTPC, twe.getTimeStamp(), twe.getTimeStampError()); - if (err < 0) { + if (err < 0 || std::abs(trLoc.getX() * trLoc.getTgl() - trLoc.getZ() - vtx.getZ()) > mSVParams->mTPCTrack2Beam) { mTracksPool[posneg].pop_back(); // discard + return true; // skip minR calculation to 'invalid' memory (trLoc will probably be still there and not be invalidated but we discarded it) } return true; } @@ -1197,6 +1243,9 @@ float SVertexer::correctTPCTrack(SVertexer::TrackCand& trc, const o2::tpc::Track } float dDrift = (tTB - tTPC.getTime0()) * mTPCBin2Z; float driftErr = tTBErr * mTPCBin2Z; + if (driftErr < 0.) { // will be discarded anyway + return driftErr; + } // eventually should be refitted, at the moment we simply shift... trc.setZ(tTPC.getZ() + (tTPC.hasASideClustersOnly() ? dDrift : -dDrift)); trc.setCov(trc.getSigmaZ2() + driftErr * driftErr, o2::track::kSigZ2);