From 81003f467dc039160dde1b79bdecf366f68eef31 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Mon, 20 Nov 2023 09:35:22 +0100 Subject: [PATCH] SVertexer: V0 Photon tuning Signed-off-by: Felix Schlepper --- .../StrangenessTrackingConfigParam.h | 3 +- .../src/StrangenessTracker.cxx | 11 +- .../include/DetectorsVertexing/SVertexer.h | 10 +- .../DetectorsVertexing/SVertexerParams.h | 56 ++++-- Detectors/Vertexing/src/SVertexer.cxx | 179 ++++++++++++++---- 5 files changed, 198 insertions(+), 61 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..0d7453d96086a 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexer.h @@ -100,8 +100,12 @@ class SVertexer VBracket vBracket{}; float minR = 0; // track lowest point r bool hasTPC = false; - uint8_t nITSclu = -1; + int8_t nITSclu = -1; bool compatibleProton = false; // dE/dx compatibility with proton hypothesis (FIXME: use better, uint8_t compat mask?) + bool hasITS() const + { + return nITSclu > 0; + } }; SVertexer(bool enabCascades = true, bool enab3body = false) : mEnableCascades{enabCascades}, mEnable3BodyDecays{enab3body} @@ -112,6 +116,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 +157,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; @@ -193,6 +198,7 @@ class SVertexer int mNThreads = 1; int mNV0s = 0, mNCascades = 0, mN3Bodies = 0, mNStrangeTracks = 0; + float mBz = 0; float mMinR2ToMeanVertex = 0; float mMaxDCAXY2ToMeanVertex = 0; float mMaxDCAXY2ToMeanVertexV0Casc = 0; diff --git a/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h b/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h index 47c56a400a289..32b2b1f7fad5a 100644 --- a/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h +++ b/Detectors/Vertexing/include/DetectorsVertexing/SVertexerParams.h @@ -40,6 +40,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 = -1.; // don't use TPC standalone tracks with X exceeding this; + float mTPCTrack2Beam = 21.f; // straight line for TPC track back to beamline + bool mTPCTrackPhotonTune = true; // use TPC-only photon tuning + int mTPCTrackMinNClusters = -1; // minimum number of clusters + float mTPCTrackXY2Radius = 90.f; // check for realistic conversion point, e.g., in material + float mTPCTrackD2R = 4.f; // check for maximal distance of pairs and their radii + float mTPCTrackDR = 90.f; // check if preliminary conversion point can lie somewhere reasonable + 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 + float maxV0TglAbsDiff = 0.3; ///< max absolute difference in Tgl for V0 for photons only 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 +117,15 @@ struct SVertexerParams : public o2::conf::ConfigurableParamHelper #include "DetectorsBase/Propagator.h" #include "TPCReconstruction/TPCFastTransformHelperO2.h" #include "DataFormatsTPC/WorkflowHelper.h" @@ -51,13 +52,13 @@ void SVertexer::process(const o2::globaltracking::RecoContainer& recoData, o2::f mStrTracker->loadData(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; @@ -72,14 +73,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; @@ -109,8 +114,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}); @@ -237,7 +242,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(); @@ -334,17 +339,18 @@ void SVertexer::setupThreads() mCascadesIdxTmp.resize(mNThreads); m3bodyIdxTmp.resize(mNThreads); mFitterV0.resize(mNThreads); - auto bz = o2::base::Propagator::Instance()->getNominalBz(); + mBz = o2::base::Propagator::Instance()->getNominalBz(); int fitCounter = 0; for (auto& fitter : mFitterV0) { fitter.setFitterID(fitCounter++); - fitter.setBz(bz); + fitter.setBz(mBz); fitter.setUseAbsDCA(mSVParams->useAbsDCA); fitter.setPropagateToPCA(false); fitter.setMaxR(mSVParams->maxRIni); 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); @@ -357,13 +363,14 @@ void SVertexer::setupThreads() fitCounter = 1000; for (auto& fitter : mFitterCasc) { fitter.setFitterID(fitCounter++); - fitter.setBz(bz); + fitter.setBz(mBz); fitter.setUseAbsDCA(mSVParams->useAbsDCA); fitter.setPropagateToPCA(false); fitter.setMaxR(mSVParams->maxRIniCasc); 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); @@ -377,13 +384,14 @@ void SVertexer::setupThreads() fitCounter = 2000; for (auto& fitter : mFitter3body) { fitter.setFitterID(fitCounter++); - fitter.setBz(bz); + fitter.setBz(mBz); fitter.setUseAbsDCA(mSVParams->useAbsDCA); fitter.setPropagateToPCA(false); fitter.setMaxR(mSVParams->maxRIni3body); 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); @@ -395,7 +403,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; @@ -503,7 +511,7 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a } // get Nclusters in the ITS if available - uint8_t nITSclu = -1; + int8_t nITSclu = -1; auto itsGID = recoData.getITSContributorGID(tvid); if (itsGID.getSource() == GIndex::ITS) { if (isITSloaded) { @@ -539,7 +547,6 @@ void SVertexer::buildT2V(const o2::globaltracking::RecoContainer& recoData) // a } } // register 1st track of each charge for each vertex - for (int pn = 0; pn < 2; pn++) { auto& vtxFirstT = mVtxFirstTrack[pn]; const auto& tracksPool = mTracksPool[pn]; @@ -559,15 +566,50 @@ 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 isTPConly = (seedP.gid.getSource() == GIndex::TPC || seedN.gid.getSource() == GIndex::TPC); + if (mSVParams->mTPCTrackPhotonTune && isTPConly) { + // Check if Tgl is close enough + if (std::abs(seedP.getTgl() - seedN.getTgl()) > mSVParams->maxV0TglAbsDiff) { + LOG(debug) << "RejTglPhi"; + return false; + } + // Check in transverse plane + float sna, csa; + o2::math_utils::CircleXYf_t trkPosCircle; + seedP.getCircleParams(mBz, trkPosCircle, sna, csa); + o2::math_utils::CircleXYf_t trkEleCircle; + seedN.getCircleParams(mBz, trkEleCircle, sna, csa); + // Does the radius of both tracks compare to their absolute circle center distance + float c2c = std::hypot(trkPosCircle.xC - trkEleCircle.xC, + trkPosCircle.yC - trkEleCircle.yC); + float r2r = trkPosCircle.rC + trkEleCircle.rC; + float dcr = c2c - r2r; + if (std::abs(dcr) > mSVParams->mTPCTrackD2R) { + LOG(debug) << "RejD2R " << c2c << " " << r2r << " " << dcr; + return false; + } + // Will the conversion point look somewhat reasonable + float r1_r = trkPosCircle.rC / r2r; + float r2_r = trkEleCircle.rC / r2r; + float dR = std::hypot(r2_r * trkPosCircle.xC + r1_r * trkEleCircle.xC, r2_r * trkPosCircle.yC + r1_r * trkEleCircle.yC); + if (dR > mSVParams->mTPCTrackDR) { + LOG(debug) << "RejDR" << dR; + return false; + } + } + + // feed DCAFitter 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; } @@ -577,13 +619,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 @@ -606,15 +649,16 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, bool goodHyp = false; std::array hypCheckStatus{}; - for (int ipid = 0; ipid < NHypV0; ipid++) { + int nPID = (mSVParams->mTPCTrackPhotonTune && isTPConly) ? (Photon + 1) : NHypV0; + for (int ipid = 0; (ipid < nPID) && 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); + bool usesTPCOnly = (seedP.hasTPC && !seedP.hasITS()) || (seedN.hasTPC && !seedN.hasITS()); if (ptV0 > mSVParams->minPtV0FromCascade && (!mSVParams->mSkipTPCOnlyCascade || !usesTPCOnly) && !usesShortITSOnly) { if (mV0Hyps[Lambda].checkTight(p2Pos, p2Neg, p2V0, ptV0) && (!mSVParams->mRequireTPCforCascBaryons || seedP.hasTPC) && seedP.compatibleProton) { goodLamForCascade = true; @@ -624,6 +668,40 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, } } + // apply loose Armenteros-Podolanski cut to photons and K0 + // the AP space and mass space correlated very well, and should in theory already be detected above but every bit helps + 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; + } + // Check if any good hypothesis remains + goodHyp = false; + for (int ipid = 0; (ipid < nPID) && mSVParams->checkV0Hypothesis; ipid++) { + if (hypCheckStatus[ipid]) { + goodHyp = true; + break; + } + } + } + // apply mass selections for 3-body decay bool good3bodyV0Hyp = false; for (int ipid = 2; ipid < 4; ipid++) { @@ -635,9 +713,16 @@ 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 checkFor3BodyDecays = mEnable3BodyDecays && + (!mSVParams->checkV0Hypothesis || good3bodyV0Hyp) && + (pt2V0 > 0.5) && + (!mSVParams->mSkipTPCOnly3Body || !isTPConly); 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) && + (!isTPConly || !hypCheckStatus[HypV0::Photon])); bool rejectIfNotCascade = false; if (!goodHyp && mSVParams->checkV0Hypothesis) { @@ -688,7 +773,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); @@ -697,9 +782,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); @@ -764,7 +849,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]; @@ -783,7 +867,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.gid.getSource() == GIndex::TPC)) { + continue; // reject TPC-only bachelors + } if (!bach.hasTPC && bach.nITSclu < mSVParams->mITSSAminNcluCascades) { continue; // reject short ITS-only } @@ -952,6 +1038,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.gid.getSource() == GIndex::TPC)) { + 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 @@ -1178,11 +1267,32 @@ 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) { mTracksPool[posneg].pop_back(); // discard + return true; } + + if (mSVParams->mTPCTrackPhotonTune) { + // require minimum of tpc clusters + bool dCls = trTPC.getNClusters() < mSVParams->mTPCTrackMinNClusters; + // check track z cuts + bool dDPV = std::abs(trLoc.getX() * trLoc.getTgl() - trLoc.getZ() - vtx.getZ()) > mSVParams->mTPCTrack2Beam; + // check track transveres cuts + float sna{0}, csa{0}; + o2::math_utils::CircleXYf_t trkCircle; + trLoc.getCircleParams(mBz, trkCircle, sna, csa); + float cR = std::hypot(trkCircle.xC, trkCircle.yC); + float drd2 = std::sqrt(cR * cR - trkCircle.rC * trkCircle.rC); + bool dRD2 = drd2 > mSVParams->mTPCTrackXY2Radius; + + if (dCls || dDPV || dRD2) { + mTracksPool[posneg].pop_back(); + return true; + } + } + return true; } @@ -1205,6 +1315,9 @@ float SVertexer::correctTPCTrack(SVertexer::TrackCand& trc, const o2::tpc::Track } float dDrift = (tTB - tTPC.getTime0()) * mTPCBin2Z; float driftErr = tTBErr * mTPCBin2Z; + if (driftErr < 0.) { // early return 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);