Skip to content

Commit a42b89c

Browse files
f3schnoferini
f3sch
authored and
noferini
committed
Add collinear DCAFitter and use it in SVertexer. (#12770)
* DCAFitter: Add option to assume collinearity of prongs Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch> * SVertexer: Fit TPConly collinearly Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch> * AnalysisDataModel: Add V0 collinear flag Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch> --------- Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>
1 parent 23e89f2 commit a42b89c

File tree

5 files changed

+57
-25
lines changed

5 files changed

+57
-25
lines changed

‎Common/DCAFitter/include/DCAFitter/DCAFitterN.h‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -188,6 +188,7 @@ class DCAFitterN
188188
void setMaxSnp(float s) { mMaxSnp = s; }
189189
void setMaxStep(float s) { mMaxStep = s; }
190190
void setMinXSeed(float x) { mMinXSeed = x; }
191+
void setCollinear(bool isCollinear) { mIsCollinear = isCollinear; }
191192

192193
int getNCandidates() const { return mCurHyp; }
193194
int getMaxIter() const { return mMaxIter; }
@@ -324,6 +325,7 @@ class DCAFitterN
324325
bool mPropagateToPCA = true; // create tracks version propagated to PCA
325326
bool mUsePropagator = false; // use propagator with 3D B-field, set automatically if material correction is requested
326327
bool mRefitWithMatCorr = false; // when doing propagateTracksToVertex, propagate tracks to V0 with material corrections and rerun minimization again
328+
bool mIsCollinear = false; // use collinear fits when there 2 crossing points
327329
o2::base::Propagator::MatCorrType mMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; // material corrections type
328330
int mMaxIter = 20; // max number of iterations
329331
float mBz = 0; // bz field, to be set by user
@@ -339,7 +341,7 @@ class DCAFitterN
339341
float mMaxStep = 2.0; // Max step for propagation with Propagator
340342
int mFitterID = 0; // locat fitter ID (mostly for debugging)
341343
size_t mCallID = 0;
342-
ClassDefNV(DCAFitterN, 1);
344+
ClassDefNV(DCAFitterN, 2);
343345
};
344346

345347
///_________________________________________________________________________
@@ -355,8 +357,8 @@ int DCAFitterN<N, Args...>::process(const Tr&... args)
355357
for (int i = 0; i < N; i++) {
356358
mTrAux[i].set(*mOrigTrPtr[i], mBz);
357359
}
358-
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
359-
return 0; // no crossing
360+
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
361+
return 0; // no crossing
360362
}
361363
for (int ih = 0; ih < MAXHYP; ih++) {
362364
mPropFailed[ih] = false;

‎Common/DCAFitter/include/DCAFitter/HelixHelper.h‎

Lines changed: 40 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ struct CrossInfo {
5959
float yDCA[2] = {};
6060
int nDCA = 0;
6161

62-
int circlesCrossInfo(const TrackAuxPar& trax0, const TrackAuxPar& trax1, float maxDistXY = MaxDistXYDef)
62+
int circlesCrossInfo(const TrackAuxPar& trax0, const TrackAuxPar& trax1, float maxDistXY = MaxDistXYDef, bool isCollinear = false)
6363
{
6464
const auto& trcA = trax0.rC > trax1.rC ? trax0 : trax1; // designate the largest circle as A
6565
const auto& trcB = trax0.rC > trax1.rC ? trax1 : trax0;
@@ -74,14 +74,24 @@ struct CrossInfo {
7474
if (dist - rsum > maxDistXY) { // too large distance
7575
return nDCA;
7676
}
77-
notTouchingXY(dist, xDist, yDist, trcA, trcB.rC);
77+
notTouchingXY(dist, xDist, yDist, trcA, trcB.rC, isCollinear);
7878
} else if (dist + trcB.rC < trcA.rC) { // the small circle is nestled into large one w/o touching
7979
// select the point of closest approach of 2 circles
80-
notTouchingXY(dist, xDist, yDist, trcA, -trcB.rC);
80+
notTouchingXY(dist, xDist, yDist, trcA, -trcB.rC, isCollinear);
8181
} else { // 2 intersection points
82-
// to simplify calculations, we move to new frame x->x+Xc0, y->y+Yc0, so that
83-
// the 1st one is centered in origin
84-
if (std::abs(xDist) < std::abs(yDist)) {
82+
if (isCollinear) {
83+
/// collinear tracks, e.g. electrons from photon conversion
84+
/// if there are 2 crossings of the circle it is better to take
85+
/// a weighted average of the crossing points as a radius
86+
float r2r = trcA.rC + trcB.rC;
87+
float r1_r = trcA.rC / r2r;
88+
float r2_r = trcB.rC / r2r;
89+
xDCA[0] = r2_r * trcA.xC + r1_r * trcB.xC;
90+
yDCA[0] = r2_r * trcA.yC + r1_r * trcB.yC;
91+
nDCA = 1;
92+
} else if (std::abs(xDist) < std::abs(yDist)) {
93+
// to simplify calculations, we move to new frame x->x+Xc0, y->y+Yc0, so that
94+
// the 1st one is centered in origin
8595
float a = (trcA.rC * trcA.rC - trcB.rC * trcB.rC + dist2) / (2. * yDist), b = -xDist / yDist, ab = a * b, bb = b * b;
8696
float det = ab * ab - (1. + bb) * (a * a - trcA.rC * trcA.rC);
8797
if (det > 0.) {
@@ -116,18 +126,28 @@ struct CrossInfo {
116126
return nDCA;
117127
}
118128

119-
void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar& trcA, float rBSign)
129+
void notTouchingXY(float dist, float xDist, float yDist, const TrackAuxPar& trcA, float rBSign, bool isCollinear = false)
120130
{
121-
// fast method to calculate DCA between 2 circles, assuming that they don't touch each outer:
122-
// the parametric equation of lines connecting the centers is x = xA + t/dist * xDist, y = yA + t/dist * yDist
123-
// with xA,yY being the center of the circle A ( = trcA.xC, trcA.yC ), xDist = trcB.xC = trcA.xC ...
124-
// There are 2 special cases:
125-
// (a) small circle is inside the large one: provide rBSign as -trcB.rC
126-
// (b) circle are side by side: provide rBSign as trcB.rC
131+
if (isCollinear) {
132+
/// for collinear tracks it is better to take
133+
/// a weighted average of the crossing points as a radius
134+
float r2r = trcA.rC + std::abs(rBSign);
135+
float r1_r = trcA.rC / r2r;
136+
float r2_r = std::abs(rBSign) / r2r;
137+
xDCA[0] = r2_r * trcA.xC + r1_r * (xDist + trcA.xC);
138+
yDCA[0] = r2_r * trcA.yC + r1_r * (yDist + trcA.yC);
139+
} else {
140+
// fast method to calculate DCA between 2 circles, assuming that they don't touch each outer:
141+
// the parametric equation of lines connecting the centers is x = xA + t/dist * xDist, y = yA + t/dist * yDist
142+
// with xA,yY being the center of the circle A ( = trcA.xC, trcA.yC ), xDist = trcB.xC = trcA.xC ...
143+
// There are 2 special cases:
144+
// (a) small circle is inside the large one: provide rBSign as -trcB.rC
145+
// (b) circle are side by side: provide rBSign as trcB.rC
146+
auto t2d = (dist + trcA.rC - rBSign) / dist;
147+
xDCA[0] = trcA.xC + 0.5 * (xDist * t2d);
148+
yDCA[0] = trcA.yC + 0.5 * (yDist * t2d);
149+
}
127150
nDCA = 1;
128-
auto t2d = (dist + trcA.rC - rBSign) / dist;
129-
xDCA[0] = trcA.xC + 0.5 * (xDist * t2d);
130-
yDCA[0] = trcA.yC + 0.5 * (yDist * t2d);
131151
}
132152

133153
template <typename T>
@@ -251,12 +271,12 @@ struct CrossInfo {
251271
}
252272

253273
template <typename T>
254-
int set(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef)
274+
int set(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef, bool isCollinear = false)
255275
{
256276
// calculate up to 2 crossings between 2 circles
257277
nDCA = 0;
258278
if (trax0.rC > o2::constants::math::Almost0 && trax1.rC > o2::constants::math::Almost0) { // both are not straight lines
259-
nDCA = circlesCrossInfo(trax0, trax1, maxDistXY);
279+
nDCA = circlesCrossInfo(trax0, trax1, maxDistXY, isCollinear);
260280
} else if (trax0.rC < o2::constants::math::Almost0 && trax1.rC < o2::constants::math::Almost0) { // both are straigt lines
261281
nDCA = linesCrossInfo(trax0, tr0, trax1, tr1, maxDistXY);
262282
} else {
@@ -269,9 +289,9 @@ struct CrossInfo {
269289
CrossInfo() = default;
270290

271291
template <typename T>
272-
CrossInfo(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef)
292+
CrossInfo(const TrackAuxPar& trax0, const T& tr0, const TrackAuxPar& trax1, const T& tr1, float maxDistXY = MaxDistXYDef, bool isCollinear = false)
273293
{
274-
set(trax0, tr0, trax1, tr1, maxDistXY);
294+
set(trax0, tr0, trax1, tr1, maxDistXY, isCollinear);
275295
}
276296
ClassDefNV(CrossInfo, 1);
277297
};

‎DataFormats/Reconstruction/include/ReconstructionDataFormats/DecayNBodyIndex.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,8 +59,10 @@ class V0Index : public DecayNBodyIndex<2>
5959
V0Index(int v, GIndex p, GIndex n) : DecayNBodyIndex<2>(v, {p, n}) {}
6060
bool isStandaloneV0() const { return testBit(0); }
6161
bool isPhotonOnly() const { return testBit(1); }
62+
bool isCollinear() const { return testBit(2); }
6263
void setStandaloneV0() { setBit(0); }
6364
void setPhotonOnly() { setBit(1); }
65+
void setCollinear() { setBit(2); }
6466
ClassDefNV(V0Index, 1);
6567
};
6668

‎Detectors/Vertexing/src/SVertexer.cxx‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -608,6 +608,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
608608
fitterV0.setMaxDZIni(mSVParams->mTPCTrackMaxDZIni);
609609
fitterV0.setMaxDXYIni(mSVParams->mTPCTrackMaxDXYIni);
610610
fitterV0.setMaxChi2(mSVParams->mTPCTrackMaxChi2);
611+
fitterV0.setCollinear(true);
611612
}
612613

613614
// feed DCAFitter
@@ -617,7 +618,9 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
617618
fitterV0.setMaxDZIni(mSVParams->maxDZIni);
618619
fitterV0.setMaxDXYIni(mSVParams->maxDXYIni);
619620
fitterV0.setMaxChi2(mSVParams->maxChi2);
621+
fitterV0.setCollinear(false);
620622
}
623+
621624
if (nCand == 0) { // discard this pair
622625
LOG(debug) << "RejDCAFitter";
623626
return false;
@@ -627,6 +630,7 @@ bool SVertexer::checkV0(const TrackCand& seedP, const TrackCand& seedN, int iP,
627630
// check closeness to the beam-line
628631
float dxv0 = v0XYZ[0] - mMeanVertex.getX(), dyv0 = v0XYZ[1] - mMeanVertex.getY(), r2v0 = dxv0 * dxv0 + dyv0 * dyv0;
629632
if (r2v0 < mMinR2ToMeanVertex) {
633+
LOG(debug) << "RejMinR2ToMeanVertex";
630634
return false;
631635
}
632636
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,
860864
}
861865
if (photonOnly) {
862866
mV0sIdxTmp[ithread].back().setPhotonOnly();
867+
mV0sIdxTmp[ithread].back().setCollinear();
863868
}
864869

865870
if (mSVParams->createFullV0s) {

‎Framework/Core/include/Framework/AnalysisDataModel.h‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1305,8 +1305,10 @@ DECLARE_SOA_COLUMN(V0Type, v0Type, uint8_t); //! cust
13051305

13061306
DECLARE_SOA_DYNAMIC_COLUMN(IsStandardV0, isStandardV0, //! is standard V0
13071307
[](uint8_t V0Type) -> bool { return V0Type & (1 << 0); });
1308-
DECLARE_SOA_DYNAMIC_COLUMN(IsPhotonV0, isPhotonV0, //! is standard V0
1308+
DECLARE_SOA_DYNAMIC_COLUMN(IsPhotonV0, isPhotonV0, //! is TPC-only V0 for which the photon-mass-hypothesis was good
13091309
[](uint8_t V0Type) -> bool { return V0Type & (1 << 1); });
1310+
DECLARE_SOA_DYNAMIC_COLUMN(IsCollinearV0, isCollinearV0, //! is V0 for which the photon-mass-hypothesis was good and was fitted collinearly
1311+
[](uint8_t V0Type) -> bool { return V0Type & (1 << 2); });
13101312

13111313
} // namespace v0
13121314

@@ -1321,7 +1323,8 @@ DECLARE_SOA_TABLE_VERSIONED(V0s_002, "AOD", "V0", 2, //! Run 3 V0 table (version
13211323
v0::PosTrackId, v0::NegTrackId,
13221324
v0::V0Type,
13231325
v0::IsStandardV0<v0::V0Type>,
1324-
v0::IsPhotonV0<v0::V0Type>);
1326+
v0::IsPhotonV0<v0::V0Type>,
1327+
v0::IsCollinearV0<v0::V0Type>);
13251328

13261329
using V0s = V0s_002; //! this defines the current default version
13271330
using V0 = V0s::iterator;

0 commit comments

Comments
 (0)