From 06c8e24b9232945fbdb03b9d183f965dd361a565 Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Fri, 15 Mar 2024 16:42:29 +0100 Subject: [PATCH] Add collinear DCAFitter and use it in SVertexer. (#12770) * DCAFitter: Add option to assume collinearity of prongs Signed-off-by: Felix Schlepper * SVertexer: Fit TPConly collinearly Signed-off-by: Felix Schlepper * AnalysisDataModel: Add V0 collinear flag Signed-off-by: Felix Schlepper --------- Signed-off-by: Felix Schlepper --- .../DCAFitter/include/DCAFitter/DCAFitterN.h | 8 ++- .../DCAFitter/include/DCAFitter/HelixHelper.h | 60 ++++++++++++------- .../DecayNBodyIndex.h | 2 + Detectors/Vertexing/src/SVertexer.cxx | 5 ++ .../include/Framework/AnalysisDataModel.h | 7 ++- 5 files changed, 57 insertions(+), 25 deletions(-) diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h index c7f6631de5abe..548cb7321e104 100644 --- a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h +++ b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h @@ -188,6 +188,7 @@ class DCAFitterN void setMaxSnp(float s) { mMaxSnp = s; } void setMaxStep(float s) { mMaxStep = s; } void setMinXSeed(float x) { mMinXSeed = x; } + void setCollinear(bool isCollinear) { mIsCollinear = isCollinear; } int getNCandidates() const { return mCurHyp; } int getMaxIter() const { return mMaxIter; } @@ -324,6 +325,7 @@ class DCAFitterN bool mPropagateToPCA = true; // create tracks version propagated to PCA bool mUsePropagator = false; // use propagator with 3D B-field, set automatically if material correction is requested bool mRefitWithMatCorr = false; // when doing propagateTracksToVertex, propagate tracks to V0 with material corrections and rerun minimization again + bool mIsCollinear = false; // use collinear fits when there 2 crossing points o2::base::Propagator::MatCorrType mMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; // material corrections type int mMaxIter = 20; // max number of iterations float mBz = 0; // bz field, to be set by user @@ -339,7 +341,7 @@ class DCAFitterN float mMaxStep = 2.0; // Max step for propagation with Propagator int mFitterID = 0; // locat fitter ID (mostly for debugging) size_t mCallID = 0; - ClassDefNV(DCAFitterN, 1); + ClassDefNV(DCAFitterN, 2); }; ///_________________________________________________________________________ @@ -355,8 +357,8 @@ int DCAFitterN::process(const Tr&... args) for (int i = 0; i < N; i++) { mTrAux[i].set(*mOrigTrPtr[i], mBz); } - if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1], mMaxDXYIni)) { // even for N>2 it should be enough to test just 1 loop - return 0; // no crossing + if (!mCrossings.set(mTrAux[0], *mOrigTrPtr[0], mTrAux[1], *mOrigTrPtr[1], mMaxDXYIni, mIsCollinear)) { // even for N>2 it should be enough to test just 1 loop + return 0; // no crossing } for (int ih = 0; ih < MAXHYP; ih++) { mPropFailed[ih] = false; diff --git a/Common/DCAFitter/include/DCAFitter/HelixHelper.h b/Common/DCAFitter/include/DCAFitter/HelixHelper.h index 87575367e1af5..39dee1d399848 100644 --- a/Common/DCAFitter/include/DCAFitter/HelixHelper.h +++ b/Common/DCAFitter/include/DCAFitter/HelixHelper.h @@ -59,7 +59,7 @@ struct CrossInfo { float yDCA[2] = {}; int nDCA = 0; - int circlesCrossInfo(const TrackAuxPar& trax0, const TrackAuxPar& trax1, float maxDistXY = MaxDistXYDef) + int circlesCrossInfo(const TrackAuxPar& trax0, const TrackAuxPar& trax1, float maxDistXY = MaxDistXYDef, bool isCollinear = false) { const auto& trcA = trax0.rC > trax1.rC ? trax0 : trax1; // designate the largest circle as A const auto& trcB = trax0.rC > trax1.rC ? trax1 : trax0; @@ -74,14 +74,24 @@ struct CrossInfo { if (dist - rsum > maxDistXY) { // too large distance return nDCA; } - notTouchingXY(dist, xDist, yDist, trcA, trcB.rC); + notTouchingXY(dist, xDist, yDist, trcA, trcB.rC, isCollinear); } else if (dist + trcB.rC < trcA.rC) { // the small circle is nestled into large one w/o touching // select the point of closest approach of 2 circles - notTouchingXY(dist, xDist, yDist, trcA, -trcB.rC); + notTouchingXY(dist, xDist, yDist, trcA, -trcB.rC, isCollinear); } else { // 2 intersection points - // to simplify calculations, we move to new frame x->x+Xc0, y->y+Yc0, so that - // the 1st one is centered in origin - if (std::abs(xDist) < std::abs(yDist)) { + if (isCollinear) { + /// collinear tracks, e.g. electrons from photon conversion + /// if there are 2 crossings of the circle it is better to take + /// a weighted average of the crossing points as a radius + float r2r = trcA.rC + trcB.rC; + float r1_r = trcA.rC / r2r; + float r2_r = trcB.rC / r2r; + xDCA[0] = r2_r * trcA.xC + r1_r * trcB.xC; + yDCA[0] = r2_r * trcA.yC + r1_r * trcB.yC; + nDCA = 1; + } else if (std::abs(xDist) < std::abs(yDist)) { + // to simplify calculations, we move to new frame x->x+Xc0, y->y+Yc0, so that + // the 1st one is centered in origin float a = (trcA.rC * trcA.rC - trcB.rC * trcB.rC + dist2) / (2. * yDist), b = -xDist / yDist, ab = a * b, bb = b * b; float det = ab * ab - (1. + bb) * (a * a - trcA.rC * trcA.rC); if (det > 0.) { @@ -116,18 +126,28 @@ struct CrossInfo { return nDCA; } - void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar& trcA, float rBSign) + void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar& trcA, float rBSign, bool isCollinear = false) { - // fast method to calculate DCA between 2 circles, assuming that they don't touch each outer: - // the parametric equation of lines connecting the centers is x = xA + t/dist * xDist, y = yA + t/dist * yDist - // with xA,yY being the center of the circle A ( = trcA.xC, trcA.yC ), xDist = trcB.xC = trcA.xC ... - // There are 2 special cases: - // (a) small circle is inside the large one: provide rBSign as -trcB.rC - // (b) circle are side by side: provide rBSign as trcB.rC + if (isCollinear) { + /// for collinear tracks it is better to take + /// a weighted average of the crossing points as a radius + float r2r = trcA.rC + std::abs(rBSign); + float r1_r = trcA.rC / r2r; + float r2_r = std::abs(rBSign) / r2r; + xDCA[0] = r2_r * trcA.xC + r1_r * (xDist + trcA.xC); + yDCA[0] = r2_r * trcA.yC + r1_r * (yDist + trcA.yC); + } else { + // fast method to calculate DCA between 2 circles, assuming that they don't touch each outer: + // the parametric equation of lines connecting the centers is x = xA + t/dist * xDist, y = yA + t/dist * yDist + // with xA,yY being the center of the circle A ( = trcA.xC, trcA.yC ), xDist = trcB.xC = trcA.xC ... + // There are 2 special cases: + // (a) small circle is inside the large one: provide rBSign as -trcB.rC + // (b) circle are side by side: provide rBSign as trcB.rC + auto t2d = (dist + trcA.rC - rBSign) / dist; + xDCA[0] = trcA.xC + 0.5 * (xDist * t2d); + yDCA[0] = trcA.yC + 0.5 * (yDist * t2d); + } nDCA = 1; - auto t2d = (dist + trcA.rC - rBSign) / dist; - xDCA[0] = trcA.xC + 0.5 * (xDist * t2d); - yDCA[0] = trcA.yC + 0.5 * (yDist * t2d); } template @@ -251,12 +271,12 @@ struct CrossInfo { } template - int set(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef) + int set(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef, bool isCollinear = false) { // calculate up to 2 crossings between 2 circles nDCA = 0; if (trax0.rC > o2::constants::math::Almost0 && trax1.rC > o2::constants::math::Almost0) { // both are not straight lines - nDCA = circlesCrossInfo(trax0, trax1, maxDistXY); + nDCA = circlesCrossInfo(trax0, trax1, maxDistXY, isCollinear); } else if (trax0.rC < o2::constants::math::Almost0 && trax1.rC < o2::constants::math::Almost0) { // both are straigt lines nDCA = linesCrossInfo(trax0, tr0, trax1, tr1, maxDistXY); } else { @@ -269,9 +289,9 @@ struct CrossInfo { CrossInfo() = default; template - CrossInfo(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef) + CrossInfo(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef, bool isCollinear = false) { - set(trax0, tr0, trax1, tr1, maxDistXY); + set(trax0, tr0, trax1, tr1, maxDistXY, isCollinear); } ClassDefNV(CrossInfo, 1); }; diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/DecayNBodyIndex.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/DecayNBodyIndex.h index 590c1b9e3ffff..31a4b8ebc44b3 100644 --- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/DecayNBodyIndex.h +++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/DecayNBodyIndex.h @@ -59,8 +59,10 @@ class V0Index : public DecayNBodyIndex<2> V0Index(int v, GIndex p, GIndex n) : DecayNBodyIndex<2>(v, {p, n}) {} bool isStandaloneV0() const { return testBit(0); } bool isPhotonOnly() const { return testBit(1); } + bool isCollinear() const { return testBit(2); } void setStandaloneV0() { setBit(0); } void setPhotonOnly() { setBit(1); } + void setCollinear() { setBit(2); } ClassDefNV(V0Index, 1); }; diff --git a/Detectors/Vertexing/src/SVertexer.cxx b/Detectors/Vertexing/src/SVertexer.cxx index 349697f5e2563..1f9eafc668cfa 100644 --- a/Detectors/Vertexing/src/SVertexer.cxx +++ b/Detectors/Vertexing/src/SVertexer.cxx @@ -608,6 +608,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, fitterV0.setMaxDZIni(mSVParams->mTPCTrackMaxDZIni); fitterV0.setMaxDXYIni(mSVParams->mTPCTrackMaxDXYIni); fitterV0.setMaxChi2(mSVParams->mTPCTrackMaxChi2); + fitterV0.setCollinear(true); } // feed DCAFitter @@ -617,7 +618,9 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, fitterV0.setMaxDZIni(mSVParams->maxDZIni); fitterV0.setMaxDXYIni(mSVParams->maxDXYIni); fitterV0.setMaxChi2(mSVParams->maxChi2); + fitterV0.setCollinear(false); } + if (nCand == 0) { // discard this pair LOG(debug) << "RejDCAFitter"; return false; @@ -627,6 +630,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, // check closeness to the beam-line float dxv0 = v0XYZ[0] - mMeanVertex.getX(), dyv0 = v0XYZ[1] - mMeanVertex.getY(), r2v0 = dxv0 * dxv0 + dyv0 * dyv0; if (r2v0 < mMinR2ToMeanVertex) { + LOG(debug) << "RejMinR2ToMeanVertex"; return false; } float rv0 = std::sqrt(r2v0), drv0P = rv0 - seedP.minR, drv0N = rv0 - seedN.minR; @@ -860,6 +864,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP, } if (photonOnly) { mV0sIdxTmp[ithread].back().setPhotonOnly(); + mV0sIdxTmp[ithread].back().setCollinear(); } if (mSVParams->createFullV0s) { diff --git a/Framework/Core/include/Framework/AnalysisDataModel.h b/Framework/Core/include/Framework/AnalysisDataModel.h index 784ff55f1402c..03221a52e114c 100644 --- a/Framework/Core/include/Framework/AnalysisDataModel.h +++ b/Framework/Core/include/Framework/AnalysisDataModel.h @@ -1305,8 +1305,10 @@ DECLARE_SOA_COLUMN(V0Type, v0Type, uint8_t); //! cust DECLARE_SOA_DYNAMIC_COLUMN(IsStandardV0, isStandardV0, //! is standard V0 [](uint8_t V0Type) -> bool { return V0Type & (1 << 0); }); -DECLARE_SOA_DYNAMIC_COLUMN(IsPhotonV0, isPhotonV0, //! is standard V0 +DECLARE_SOA_DYNAMIC_COLUMN(IsPhotonV0, isPhotonV0, //! is TPC-only V0 for which the photon-mass-hypothesis was good [](uint8_t V0Type) -> bool { return V0Type & (1 << 1); }); +DECLARE_SOA_DYNAMIC_COLUMN(IsCollinearV0, isCollinearV0, //! is V0 for which the photon-mass-hypothesis was good and was fitted collinearly + [](uint8_t V0Type) -> bool { return V0Type & (1 << 2); }); } // namespace v0 @@ -1321,7 +1323,8 @@ DECLARE_SOA_TABLE_VERSIONED(V0s_002, "AOD", "V0", 2, //! Run 3 V0 table (version v0::PosTrackId, v0::NegTrackId, v0::V0Type, v0::IsStandardV0, - v0::IsPhotonV0); + v0::IsPhotonV0, + v0::IsCollinearV0); using V0s = V0s_002; //! this defines the current default version using V0 = V0s::iterator;