Skip to content

Commit

Permalink
Add collinear DCAFitter and use it in SVertexer. (#12770)
Browse files Browse the repository at this point in the history 
* 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>
  • Loading branch information
@f3sch @noferini
f3sch authored and noferini committed Mar 15, 2024
1 parent 9f9a131 commit 06c8e24
Show file tree
Hide file tree
Showing 5 changed files with 57 additions and 25 deletions.
8 changes: 5 additions & 3 deletions Common/DCAFitter/include/DCAFitter/DCAFitterN.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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; }
Expand Down Expand Up @@ -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
Expand All @@ -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);
};

///_________________________________________________________________________
Expand All @@ -355,8 +357,8 @@ int DCAFitterN<N, Args...>::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;
Expand Down
60 changes: 40 additions & 20 deletions Common/DCAFitter/include/DCAFitter/HelixHelper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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.) {
Expand Down Expand Up @@ -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 <typename T>
Expand Down Expand Up @@ -251,12 +271,12 @@ struct CrossInfo {
}

template <typename T>
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 {
Expand All @@ -269,9 +289,9 @@ struct CrossInfo {
CrossInfo() = default;

template <typename T>
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);
};
Expand Down
2 changes: 2 additions & 0 deletions DataFormats/Reconstruction/include/ReconstructionDataFormats/DecayNBodyIndex.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
};

Expand Down
5 changes: 5 additions & 0 deletions Detectors/Vertexing/src/SVertexer.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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) {
Expand Down
7 changes: 5 additions & 2 deletions Framework/Core/include/Framework/AnalysisDataModel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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::V0Type>,
v0::IsPhotonV0<v0::V0Type>);
v0::IsPhotonV0<v0::V0Type>,
v0::IsCollinearV0<v0::V0Type>);

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

0 comments on commit 06c8e24

Please sign in to comment.