Project
Loading...
Searching...
No Matches
DCAFitterN.h
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
16
17#ifndef _ALICEO2_DCA_FITTERN_
18#define _ALICEO2_DCA_FITTERN_
19
22#include "MathUtils/Cartesian.h"
24
25namespace o2
26{
27namespace vertexing
28{
29
32struct TrackCovI {
33 float sxx, syy, syz, szz;
34
35 GPUdDefault() TrackCovI() = default;
36
37 GPUd() bool set(const o2::track::TrackParCov& trc, float xerrFactor = 1.f)
38 {
39 // we assign Y error to X for DCA calculation
40 // (otherwise for quazi-collinear tracks the X will not be constrained)
41 float cyy = trc.getSigmaY2(), czz = trc.getSigmaZ2(), cyz = trc.getSigmaZY(), cxx = cyy * xerrFactor;
42 float detYZ = cyy * czz - cyz * cyz;
43 bool res = true;
44 if (detYZ <= 0.) {
45 cyz = o2::gpu::GPUCommonMath::Sqrt(cyy * czz) * (cyz > 0 ? 0.98f : -0.98f);
46 detYZ = cyy * czz - cyz * cyz;
47 res = false;
48 }
49 auto detYZI = 1. / detYZ;
50 sxx = 1. / cxx;
51 syy = czz * detYZI;
52 syz = -cyz * detYZI;
53 szz = cyy * detYZI;
54 return res;
55 }
56};
57
60struct TrackDeriv {
62 GPUdDefault() TrackDeriv() = default;
63 GPUd() TrackDeriv(const o2::track::TrackPar& trc, float bz) { set(trc, bz); }
64 GPUd() void set(const o2::track::TrackPar& trc, float bz)
65 {
66 float snp = trc.getSnp(), csp = o2::gpu::GPUCommonMath::Sqrt((1. - snp) * (1. + snp)), cspI = 1. / csp, crv2c = trc.getCurvature(bz) * cspI;
67 dydx = snp * cspI; // = snp/csp
68 dzdx = trc.getTgl() * cspI; // = tgl/csp
69 d2ydx2 = crv2c * cspI * cspI; // = crv/csp^3
70 d2zdx2 = crv2c * dzdx * dydx; // = crv*tgl*snp/csp^3
71 }
72};
73
77 size_t evCount{0};
78 size_t nextLog{1};
79 GPUdi() bool needToLog()
80 {
81 if (++evCount > nextLog) {
82 nextLog *= 2;
83 return true;
84 }
85 return false;
86 }
88 {
89 evCount = 0;
90 nextLog = 1;
91 }
92};
93
94template <int N, typename... Args>
96{
97 static constexpr double NMin = 2;
98 static constexpr double NMax = 4;
99 static constexpr double NInv = 1. / N;
100 static constexpr int MAXHYP = 2;
101 static constexpr float XerrFactor = 5.; // factor for conversion of track covYY to dummy covXX
102 using Track = o2::track::TrackParCov;
105
112 using TrackCoefVtx = MatStd3D;
113 using ArrTrack = std::array<Track, N>; // container for prongs (tracks) at single vertex cand.
114 using ArrTrackCovI = std::array<TrackCovI, N>; // container for inv.cov.matrices at single vertex cand.
115 using ArrTrCoef = std::array<TrackCoefVtx, N>; // container of TrackCoefVtx coefficients at single vertex cand.
116 using ArrTrDer = std::array<TrackDeriv, N>; // container of Track 1st and 2nd derivative over their X param
117 using ArrTrPos = std::array<Vec3D, N>; // container of Track positions
118
119 public:
120 enum BadCovPolicy : uint8_t { // if encountering non-positive defined cov. matrix, the choice is:
121 Discard = 0, // stop evaluation
122 Override = 1, // override correlation coef. to have cov.matrix pos.def and continue
123 OverrideAndFlag = 2 // override correlation coef. to have cov.matrix pos.def, set mPropFailed flag of corresponding candidate to true and continue (up to the user to check the flag)
124 };
125
126 enum FitStatus : uint8_t { // fit status of crossing hypothesis
127 None, // no status set (should not be possible!)
128
129 /* Good Conditions */
130 Converged, // fit converged
131 MaxIter, // max iterations reached before fit convergence
132
133 /* Error Conditions */
134 NoCrossing, // no reasaonable crossing was found
135 RejRadius, // radius of crossing was not acceptable
136 RejTrackX, // one candidate track x was below the mimimum required radius
137 RejTrackRoughZ, // rejected by rough cut on tracks Z difference
138 RejChi2Max, // rejected by maximum chi2 cut
139 FailProp, // propagation of at least prong to PCA failed
140 FailInvCov, // inversion of cov.-matrix failed
141 FailInvWeight, // inversion of Ti weight matrix failed
142 FailInv2ndDeriv, // inversion of 2nd derivatives failed
143 FailCorrTracks, // correction of tracks to updated x failed
144 FailCloserAlt, // alternative PCA is closer
145 //
147 };
148
149 static constexpr int getNProngs() { return N; }
150
151 DCAFitterN() = default;
152 DCAFitterN(float bz, bool useAbsDCA, bool prop2DCA) : mBz(bz), mUseAbsDCA(useAbsDCA), mPropagateToPCA(prop2DCA)
153 {
154 static_assert(N >= NMin && N <= NMax, "N prongs outside of allowed range");
155 }
156
157 //=========================================================================
159 GPUd() const Vec3D& getPCACandidate(int cand = 0) const { return mPCA[mOrder[cand]]; }
160 GPUd() const auto getPCACandidatePos(int cand = 0) const
161 {
162 const auto& vd = mPCA[mOrder[cand]];
163 return std::array<float, 3>{static_cast<float>(vd[0]), static_cast<float>(vd[1]), static_cast<float>(vd[2])};
164 }
165
167 int getCandidatePosition(int cand = 0) const { return mOrder[cand]; }
168
170 float getChi2AtPCACandidate(int cand = 0) const { return mChi2[mOrder[cand]]; }
171
174 GPUd() bool propagateTracksToVertex(int cand = 0);
175
177 GPUd() bool isPropagateTracksToVertexDone(int cand = 0) const { return mTrPropDone[mOrder[cand]]; }
178
180 bool isPropagationFailure(int cand = 0) const { return mPropFailed[mOrder[cand]]; }
181
184 Track& getTrack(int i, int cand = 0)
185 {
186 if (!mTrPropDone[mOrder[cand]]) {
187#ifndef GPUCA_GPUCODE_DEVICE
188 throw std::runtime_error("propagateTracksToVertex was not called yet");
189#endif
190 }
191 return mCandTr[mOrder[cand]][i];
192 }
193
194 const Track& getTrack(int i, int cand = 0) const
195 {
196 if (!mTrPropDone[mOrder[cand]]) {
197#ifndef GPUCA_GPUCODE_DEVICE
198 throw std::runtime_error("propagateTracksToVertex was not called yet");
199#endif
200 }
201 return mCandTr[mOrder[cand]][i];
202 }
203
205 GPUd() o2::track::TrackParCov createParentTrackParCov(int cand = 0, bool sectorAlpha = true) const;
206
208 GPUd() o2::track::TrackPar createParentTrackPar(int cand = 0, bool sectorAlpha = true) const;
209
211 GPUd() o2::track::TrackPar getTrackParamAtPCA(int i, int cand = 0);
212
214 GPUd() bool recalculatePCAWithErrors(int cand = 0);
215
216 GPUd() MatSym3D calcPCACovMatrix(int cand = 0) const;
217
218 std::array<float, 6> calcPCACovMatrixFlat(int cand = 0) const
219 {
220 auto m = calcPCACovMatrix(cand);
221 return {static_cast<float>(m(0, 0)), static_cast<float>(m(1, 0)), static_cast<float>(m(1, 1)), static_cast<float>(m(2, 0)), static_cast<float>(m(2, 1)), static_cast<float>(m(2, 2))};
222 }
223
224 const Track* getOrigTrackPtr(int i) const { return mOrigTrPtr[i]; }
225
226 GPUdi() FitStatus getFitStatus(int cand = 0) const noexcept { return mFitStatus[mOrder[cand]]; }
227
229 GPUdi() int getNIterations(int cand = 0) const { return mNIters[mOrder[cand]]; }
230 GPUdi() void setPropagateToPCA(bool v = true) { mPropagateToPCA = v; }
231 GPUdi() void setMaxIter(int n = 20) { mMaxIter = n > 2 ? n : 2; }
232 GPUdi() void setMaxR(float r = 200.) { mMaxR2 = r * r; }
233 GPUdi() void setMaxDZIni(float d = 4.) { mMaxDZIni = d; }
234 GPUdi() void setMaxDXYIni(float d = 4.) { mMaxDXYIni = d > 0 ? d : 1e9; }
235 GPUdi() void setMaxChi2(float chi2 = 999.) { mMaxChi2 = chi2; }
236 GPUdi() void setBz(float bz) { mBz = o2::gpu::GPUCommonMath::Abs(bz) > o2::constants::math::Almost0 ? bz : 0.f; }
237 GPUdi() void setMinParamChange(float x = 1e-3) { mMinParamChange = x > 1e-4 ? x : 1.e-4; }
238 GPUdi() void setMinRelChi2Change(float r = 0.9) { mMinRelChi2Change = r > 0.1 ? r : 999.; }
239 GPUdi() void setUseAbsDCA(bool v) { mUseAbsDCA = v; }
240 GPUdi() void setWeightedFinalPCA(bool v) { mWeightedFinalPCA = v; }
241 GPUdi() void setMaxDistance2ToMerge(float v) { mMaxDist2ToMergeSeeds = v; }
242 GPUdi() void setMatCorrType(o2::base::Propagator::MatCorrType m = o2::base::Propagator::MatCorrType::USEMatCorrLUT) { mMatCorr = m; }
243 GPUdi() void setUsePropagator(bool v) { mUsePropagator = v; }
244 GPUdi() void setRefitWithMatCorr(bool v) { mRefitWithMatCorr = v; }
245 GPUdi() void setMaxSnp(float s) { mMaxSnp = s; }
246 GPUdi() void setMaxStep(float s) { mMaxStep = s; }
247 GPUdi() void setMinXSeed(float x) { mMinXSeed = x; }
248 GPUdi() void setCollinear(bool isCollinear) { mIsCollinear = isCollinear; }
249
250 GPUdi() int getNCandidates() const { return mCurHyp; }
251 GPUdi() int getMaxIter() const { return mMaxIter; }
252 GPUdi() float getMaxR() const { return o2::gpu::GPUCommonMath::Sqrt(mMaxR2); }
253 GPUdi() float getMaxDZIni() const { return mMaxDZIni; }
254 GPUdi() float getMaxDXYIni() const { return mMaxDXYIni; }
255 GPUdi() float getMaxChi2() const { return mMaxChi2; }
256 GPUdi() float getMinParamChange() const { return mMinParamChange; }
257 GPUdi() float getBz() const { return mBz; }
258 GPUdi() float getMaxDistance2ToMerge() const { return mMaxDist2ToMergeSeeds; }
259 GPUdi() bool getUseAbsDCA() const { return mUseAbsDCA; }
260 GPUdi() bool getWeightedFinalPCA() const { return mWeightedFinalPCA; }
261 GPUdi() bool getPropagateToPCA() const { return mPropagateToPCA; }
262 GPUdi() o2::base::Propagator::MatCorrType getMatCorrType() const { return mMatCorr; }
263 GPUdi() bool getUsePropagator() const { return mUsePropagator; }
264 GPUdi() bool getRefitWithMatCorr() const { return mRefitWithMatCorr; }
265 GPUdi() float getMaxSnp() const { return mMaxSnp; }
266 GPUdi() float getMasStep() const { return mMaxStep; }
267 GPUdi() float getMinXSeed() const { return mMinXSeed; }
268
269 template <class... Tr>
270 GPUd() int process(const Tr&... args);
271 GPUd() void print() const;
272
273 GPUdi() int getFitterID() const { return mFitterID; }
274 GPUdi() void setFitterID(int i) { mFitterID = i; }
275 GPUdi() size_t getCallID() const { return mCallID; }
276
277 protected:
278 GPUd() bool calcPCACoefs();
279 GPUd() bool calcInverseWeight();
280 GPUd() void calcResidDerivatives();
281 GPUd() void calcResidDerivativesNoErr();
282 GPUd() void calcRMatrices();
283 GPUd() void calcChi2Derivatives();
284 GPUd() void calcChi2DerivativesNoErr();
285 GPUd() void calcPCA();
286 GPUd() void calcPCANoErr();
287 GPUd() void calcTrackResiduals();
288 GPUd() void calcTrackDerivatives();
289 GPUd() double calcChi2() const;
290 GPUd() double calcChi2NoErr() const;
291 GPUd() bool correctTracks(const VecND& corrX);
292 GPUd() bool minimizeChi2();
293 GPUd() bool minimizeChi2NoErr();
294 GPUd() bool roughDZCut() const;
295 GPUd() bool closerToAlternative() const;
296 GPUd() bool propagateToX(o2::track::TrackParCov& t, float x);
297 GPUd() bool propagateParamToX(o2::track::TrackPar& t, float x);
298
299 GPUd() static double getAbsMax(const VecND& v);
301 GPUdi() const Vec3D& getTrackPos(int i, int cand = 0) const { return mTrPos[mOrder[cand]][i]; }
302
304 GPUd() float getTrackX(int i, int cand = 0) const { return getTrackPos(i, cand)[0]; }
305
306 GPUd() MatStd3D getTrackRotMatrix(int i) const // generate 3D matrix for track rotation to global frame
307 {
308 MatStd3D mat;
309 mat(2, 2) = 1;
310 mat(0, 0) = mat(1, 1) = mTrAux[i].c;
311 mat(0, 1) = -mTrAux[i].s;
312 mat(1, 0) = mTrAux[i].s;
313 return mat;
314 }
315
316 GPUd() MatSym3D getTrackCovMatrix(int i, int cand = 0) const // generate covariance matrix of track position, adding fake X error
317 {
318 const auto& trc = mCandTr[mOrder[cand]][i];
319 MatSym3D mat;
320 mat(0, 0) = trc.getSigmaY2() * XerrFactor;
321 mat(1, 1) = trc.getSigmaY2();
322 mat(2, 2) = trc.getSigmaZ2();
323 mat(2, 1) = trc.getSigmaZY();
324 return mat;
325 }
326
327 GPUd() void assign(int) {}
328 template <class T, class... Tr>
329 GPUd() void assign(int i, const T& t, const Tr&... args)
330 {
331#ifndef GPUCA_GPUCODE_DEVICE
332 static_assert(std::is_convertible<T, Track>(), "Wrong track type");
333#endif
335 assign(i + 1, args...);
336 }
337
339 {
340 mCurHyp = 0;
341 mAllowAltPreference = true;
342 mOrder.fill(0);
343 mPropFailed.fill(false);
344 mTrPropDone.fill(false);
345 mNIters.fill(0);
346 mChi2.fill(-1);
347 mFitStatus.fill(FitStatus::None);
348 }
349
350 GPUdi() static void setTrackPos(Vec3D& pnt, const Track& tr)
351 {
352 pnt[0] = tr.getX();
353 pnt[1] = tr.getY();
354 pnt[2] = tr.getZ();
355 }
356
357 GPUdi() void clearLogThrottlers()
358 {
359 mLoggerBadCov.clear();
360 mLoggerBadInv.clear();
361 mLoggerBadProp.clear();
362 }
363
364 void setBadCovPolicy(BadCovPolicy v) { mBadCovPolicy = v; }
365 BadCovPolicy getBadCovPolicy() const { return mBadCovPolicy; }
366
367 private:
368 // vectors of 1st derivatives of track local residuals over X parameters
369 std::array<std::array<Vec3D, N>, N> mDResidDx;
370 // vectors of 1nd derivatives of track local residuals over X parameters
371 // (cross-derivatives DR/(dx_j*dx_k) = 0 for j!=k, therefore the hessian is diagonal)
372 std::array<std::array<Vec3D, N>, N> mD2ResidDx2;
373 VecND mDChi2Dx; // 1st derivatives of chi2 over tracks X params
374 MatSymND mD2Chi2Dx2; // 2nd derivatives of chi2 over tracks X params (symmetric matrix)
375 MatSymND mCosDif; // matrix with cos(alp_j-alp_i) for j<i
376 MatSymND mSinDif; // matrix with sin(alp_j-alp_i) for j<i
377 std::array<const Track*, N> mOrigTrPtr;
378 std::array<TrackAuxPar, N> mTrAux; // Aux track info for each track at each cand. vertex
379 CrossInfo mCrossings; // info on track crossing
380
381 std::array<ArrTrackCovI, MAXHYP> mTrcEInv; // errors for each track at each cand. vertex
382 std::array<ArrTrack, MAXHYP> mCandTr; // tracks at each cond. vertex (Note: Errors are at seed XY point)
383 std::array<ArrTrCoef, MAXHYP> mTrCFVT; // TrackCoefVtx for each track at each cand. vertex
384 std::array<ArrTrDer, MAXHYP> mTrDer; // Track derivativse
385 std::array<ArrTrPos, MAXHYP> mTrPos; // Track positions
386 std::array<ArrTrPos, MAXHYP> mTrRes; // Track residuals
387 std::array<Vec3D, MAXHYP> mPCA; // PCA for each vertex candidate
388 std::array<float, MAXHYP> mChi2 = {0}; // Chi2 at PCA candidate
389 std::array<int, MAXHYP> mNIters; // number of iterations for each seed
390 std::array<bool, MAXHYP> mTrPropDone{}; // Flag that the tracks are fully propagated to PCA
391 std::array<bool, MAXHYP> mPropFailed{}; // Flag that some propagation failed for this PCA candidate
392 LogLogThrottler mLoggerBadCov{};
393 LogLogThrottler mLoggerBadInv{};
394 LogLogThrottler mLoggerBadProp{};
395 MatSym3D mWeightInv; // inverse weight of single track, [sum{M^T E M}]^-1 in EQ.T
396 std::array<int, MAXHYP> mOrder{0};
397 int mCurHyp = 0;
398 int mCrossIDCur = 0;
399 int mCrossIDAlt = -1;
400 BadCovPolicy mBadCovPolicy{BadCovPolicy::Discard}; // what to do in case of non-pos-def. cov. matrix, see BadCovPolicy enum
401 std::array<FitStatus, MAXHYP> mFitStatus{}; // fit status of each hypothesis fit
402 bool mAllowAltPreference = true; // if the fit converges to alternative PCA seed, abandon the current one
403 bool mUseAbsDCA = false; // use abs. distance minimization rather than chi2
404 bool mWeightedFinalPCA = false; // recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used
405 bool mPropagateToPCA = true; // create tracks version propagated to PCA
406 bool mUsePropagator = false; // use propagator with 3D B-field, set automatically if material correction is requested
407 bool mRefitWithMatCorr = false; // when doing propagateTracksToVertex, propagate tracks to V0 with material corrections and rerun minimization again
408 bool mIsCollinear = false; // use collinear fits when there 2 crossing points
409 o2::base::Propagator::MatCorrType mMatCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; // material corrections type
410 int mMaxIter = 20; // max number of iterations
411 float mBz = 0; // bz field, to be set by user
412 float mMaxR2 = 200. * 200.; // reject PCA's above this radius
413 float mMinXSeed = -50.; // reject seed if it corresponds to X-param < mMinXSeed for one of candidates (e.g. X becomes strongly negative)
414 float mMaxDZIni = 4.; // reject (if>0) PCA candidate if tracks DZ exceeds threshold
415 float mMaxDXYIni = 4.; // reject (if>0) PCA candidate if tracks dXY exceeds threshold
416 float mMinParamChange = 1e-3; // stop iterations if largest change of any X is smaller than this
417 float mMinRelChi2Change = 0.9; // stop iterations is chi2/chi2old > this
418 float mMaxChi2 = 100; // abs cut on chi2 or abs distance
419 float mMaxDist2ToMergeSeeds = 1.; // merge 2 seeds to their average if their distance^2 is below the threshold
420 float mMaxSnp = 0.95; // Max snp for propagation with Propagator
421 float mMaxStep = 2.0; // Max step for propagation with Propagator
422 int mFitterID = 0; // locat fitter ID (mostly for debugging)
423 size_t mCallID = 0;
424 ClassDefNV(DCAFitterN, 3);
425};
426
428template <int N, typename... Args>
429template <class... Tr>
430GPUd() int DCAFitterN<N, Args...>::process(const Tr&... args)
431{
432 // This is a main entry point: fit PCA of N tracks
433 mCallID++;
434 static_assert(sizeof...(args) == N, "incorrect number of input tracks");
435 assign(0, args...);
436 clear();
437 for (int i = 0; i < N; i++) {
438 mTrAux[i].set(*mOrigTrPtr[i], mBz);
439 }
440 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
441 mFitStatus[mCurHyp] = FitStatus::NoCrossing;
442 return 0;
443 }
444 if (mUseAbsDCA) {
445 calcRMatrices(); // needed for fast residuals derivatives calculation in case of abs. distance minimization
446 }
447 if (mCrossings.nDCA == MAXHYP) { // if there are 2 candidates and they are too close, chose their mean as a starting point
448 auto dst2 = (mCrossings.xDCA[0] - mCrossings.xDCA[1]) * (mCrossings.xDCA[0] - mCrossings.xDCA[1]) +
449 (mCrossings.yDCA[0] - mCrossings.yDCA[1]) * (mCrossings.yDCA[0] - mCrossings.yDCA[1]);
450 if (dst2 < mMaxDist2ToMergeSeeds) {
451 mCrossings.nDCA = 1;
452 mCrossings.xDCA[0] = 0.5 * (mCrossings.xDCA[0] + mCrossings.xDCA[1]);
453 mCrossings.yDCA[0] = 0.5 * (mCrossings.yDCA[0] + mCrossings.yDCA[1]);
454 }
455 }
456 // check all crossings
457 for (int ic = 0; ic < mCrossings.nDCA; ic++) {
458 // check if radius is acceptable
459 if (mCrossings.xDCA[ic] * mCrossings.xDCA[ic] + mCrossings.yDCA[ic] * mCrossings.yDCA[ic] > mMaxR2) {
460 mFitStatus[mCurHyp] = FitStatus::RejRadius;
461 continue;
462 }
463 mCrossIDCur = ic;
464 mCrossIDAlt = (mCrossings.nDCA == 2 && mAllowAltPreference) ? 1 - ic : -1; // works for max 2 crossings
465 mPCA[mCurHyp][0] = mCrossings.xDCA[ic];
466 mPCA[mCurHyp][1] = mCrossings.yDCA[ic];
467
468 if (mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2()) {
469 mOrder[mCurHyp] = mCurHyp;
470 if (mPropagateToPCA && !propagateTracksToVertex(mCurHyp)) {
471 continue; // discard candidate if failed to propagate to it
472 }
473 mCurHyp++;
474 }
475 }
476
477 for (int i = mCurHyp; i--;) { // order in quality
478 for (int j = i; j--;) {
479 if (mChi2[mOrder[i]] < mChi2[mOrder[j]]) {
480 o2::gpu::GPUCommonMath::Swap(mOrder[i], mOrder[j]);
481 }
482 }
483 }
484 if (mUseAbsDCA && mWeightedFinalPCA) {
485 for (int i = mCurHyp; i--;) {
486 recalculatePCAWithErrors(i);
487 }
488 }
489 return mCurHyp;
490}
491
492//__________________________________________________________________________
493template <int N, typename... Args>
494GPUd() bool DCAFitterN<N, Args...>::calcPCACoefs()
495{
496 //< calculate Ti matrices for global vertex decomposition to V = sum_{0<i<N} Ti pi, see EQ.T in the ref
497 if (!calcInverseWeight()) {
498 mFitStatus[mCurHyp] = FitStatus::FailInvWeight;
499 return false;
500 }
501 for (int i = N; i--;) { // build Mi*Ei matrix
502 const auto& taux = mTrAux[i];
503 const auto& tcov = mTrcEInv[mCurHyp][i];
504 MatStd3D miei;
505 miei[0][0] = taux.c * tcov.sxx;
506 miei[0][1] = -taux.s * tcov.syy;
507 miei[0][2] = -taux.s * tcov.syz;
508 miei[1][0] = taux.s * tcov.sxx;
509 miei[1][1] = taux.c * tcov.syy;
510 miei[1][2] = taux.c * tcov.syz;
511 miei[2][0] = 0;
512 miei[2][1] = tcov.syz;
513 miei[2][2] = tcov.szz;
514 mTrCFVT[mCurHyp][i] = mWeightInv * miei;
515 }
516 return true;
517}
518
519//__________________________________________________________________________
520template <int N, typename... Args>
521GPUd() bool DCAFitterN<N, Args...>::calcInverseWeight()
522{
523 //< calculate [sum_{0<j<N} M_j*E_j*M_j^T]^-1 used for Ti matrices, see EQ.T
524 auto* arrmat = mWeightInv.Array();
525 memset(arrmat, 0, sizeof(mWeightInv));
526 enum { XX,
527 XY,
528 YY,
529 XZ,
530 YZ,
531 ZZ };
532 for (int i = N; i--;) {
533 const auto& taux = mTrAux[i];
534 const auto& tcov = mTrcEInv[mCurHyp][i];
535 arrmat[XX] += taux.cc * tcov.sxx + taux.ss * tcov.syy;
536 arrmat[XY] += taux.cs * (tcov.sxx - tcov.syy);
537 arrmat[XZ] += -taux.s * tcov.syz;
538 arrmat[YY] += taux.cc * tcov.syy + taux.ss * tcov.sxx;
539 arrmat[YZ] += taux.c * tcov.syz;
540 arrmat[ZZ] += tcov.szz;
541 }
542 // invert 3x3 symmetrix matrix
543 return mWeightInv.Invert();
544}
545
546//__________________________________________________________________________
547template <int N, typename... Args>
548GPUd() void DCAFitterN<N, Args...>::calcResidDerivatives()
549{
550 //< calculate matrix of derivatives for weighted chi2: residual i vs parameter X of track j
551 MatStd3D matMT;
552 for (int i = N; i--;) { // residual being differentiated
553 const auto& taux = mTrAux[i];
554 for (int j = N; j--;) { // track over which we differentiate
555 const auto& matT = mTrCFVT[mCurHyp][j]; // coefficient matrix for track J
556 const auto& trDx = mTrDer[mCurHyp][j]; // track point derivs over track X param
557 auto& dr1 = mDResidDx[i][j];
558 auto& dr2 = mD2ResidDx2[i][j];
559 // calculate M_i^tr * T_j
560 matMT[0][0] = taux.c * matT[0][0] + taux.s * matT[1][0];
561 matMT[0][1] = taux.c * matT[0][1] + taux.s * matT[1][1];
562 matMT[0][2] = taux.c * matT[0][2] + taux.s * matT[1][2];
563 matMT[1][0] = -taux.s * matT[0][0] + taux.c * matT[1][0];
564 matMT[1][1] = -taux.s * matT[0][1] + taux.c * matT[1][1];
565 matMT[1][2] = -taux.s * matT[0][2] + taux.c * matT[1][2];
566 matMT[2][0] = matT[2][0];
567 matMT[2][1] = matT[2][1];
568 matMT[2][2] = matT[2][2];
569
570 // calculate DResid_i/Dx_j = (delta_ij - M_i^tr * T_j) * DTrack_k/Dx_k
571 dr1[0] = -(matMT[0][0] + matMT[0][1] * trDx.dydx + matMT[0][2] * trDx.dzdx);
572 dr1[1] = -(matMT[1][0] + matMT[1][1] * trDx.dydx + matMT[1][2] * trDx.dzdx);
573 dr1[2] = -(matMT[2][0] + matMT[2][1] * trDx.dydx + matMT[2][2] * trDx.dzdx);
574
575 // calculate D2Resid_I/(Dx_J Dx_K) = (delta_ijk - M_i^tr * T_j * delta_jk) * D2Track_k/dx_k^2
576 dr2[0] = -(matMT[0][1] * trDx.d2ydx2 + matMT[0][2] * trDx.d2zdx2);
577 dr2[1] = -(matMT[1][1] * trDx.d2ydx2 + matMT[1][2] * trDx.d2zdx2);
578 dr2[2] = -(matMT[2][1] * trDx.d2ydx2 + matMT[2][2] * trDx.d2zdx2);
579
580 if (i == j) {
581 dr1[0] += 1.;
582 dr1[1] += trDx.dydx;
583 dr1[2] += trDx.dzdx;
584
585 dr2[1] += trDx.d2ydx2;
586 dr2[2] += trDx.d2zdx2;
587 }
588 } // track over which we differentiate
589 } // residual being differentiated
590}
591
592//__________________________________________________________________________
593template <int N, typename... Args>
594GPUd() void DCAFitterN<N, Args...>::calcResidDerivativesNoErr()
595{
596 //< calculate matrix of derivatives for absolute distance chi2: residual i vs parameter X of track j
597 constexpr double NInv1 = 1. - NInv; // profit from Rii = I/Ninv
598 for (int i = N; i--;) { // residual being differentiated
599 const auto& trDxi = mTrDer[mCurHyp][i]; // track point derivs over track X param
600 auto& dr1ii = mDResidDx[i][i];
601 auto& dr2ii = mD2ResidDx2[i][i];
602 dr1ii[0] = NInv1;
603 dr1ii[1] = NInv1 * trDxi.dydx;
604 dr1ii[2] = NInv1 * trDxi.dzdx;
605
606 dr2ii[0] = 0;
607 dr2ii[1] = NInv1 * trDxi.d2ydx2;
608 dr2ii[2] = NInv1 * trDxi.d2zdx2;
609
610 for (int j = i; j--;) { // track over which we differentiate
611 auto& dr1ij = mDResidDx[i][j];
612 auto& dr1ji = mDResidDx[j][i];
613 const auto& trDxj = mTrDer[mCurHyp][j]; // track point derivs over track X param
614 auto cij = mCosDif[i][j], sij = mSinDif[i][j]; // M_i^T*M_j / N matrices non-trivial elements = {ci*cj+si*sj , si*cj-ci*sj }, see 5 in ref.
615
616 // calculate DResid_i/Dx_j = (delta_ij - R_ij) * DTrack_j/Dx_j for j<i
617 dr1ij[0] = -(cij + sij * trDxj.dydx);
618 dr1ij[1] = -(-sij + cij * trDxj.dydx);
619 dr1ij[2] = -trDxj.dzdx * NInv;
620
621 // calculate DResid_j/Dx_i = (delta_ij - R_ji) * DTrack_i/Dx_i for j<i
622 dr1ji[0] = -(cij - sij * trDxi.dydx);
623 dr1ji[1] = -(sij + cij * trDxi.dydx);
624 dr1ji[2] = -trDxi.dzdx * NInv;
625
626 auto& dr2ij = mD2ResidDx2[i][j];
627 auto& dr2ji = mD2ResidDx2[j][i];
628 // calculate D2Resid_I/(Dx_J Dx_K) = (delta_ij - Rij) * D2Track_j/dx_j^2 * delta_jk for j<i
629 dr2ij[0] = -sij * trDxj.d2ydx2;
630 dr2ij[1] = -cij * trDxj.d2ydx2;
631 dr2ij[2] = -trDxj.d2zdx2 * NInv;
632
633 // calculate D2Resid_j/(Dx_i Dx_k) = (delta_ij - Rji) * D2Track_i/dx_i^2 * delta_ik for j<i
634 dr2ji[0] = sij * trDxi.d2ydx2;
635 dr2ji[1] = -cij * trDxi.d2ydx2;
636 dr2ji[2] = -trDxi.d2zdx2 * NInv;
637
638 } // track over which we differentiate
639 } // residual being differentiated
640}
641
642//__________________________________________________________________________
643template <int N, typename... Args>
644GPUd() void DCAFitterN<N, Args...>::calcRMatrices()
645{
646 //< calculate Rij = 1/N M_i^T * M_j matrices (rotation from j-th track to i-th track frame)
647 for (int i = N; i--;) {
648 const auto& mi = mTrAux[i];
649 for (int j = i; j--;) {
650 const auto& mj = mTrAux[j];
651 mCosDif[i][j] = (mi.c * mj.c + mi.s * mj.s) * NInv; // cos(alp_i-alp_j) / N
652 mSinDif[i][j] = (mi.s * mj.c - mi.c * mj.s) * NInv; // sin(alp_i-alp_j) / N
653 }
654 }
655}
656
657//__________________________________________________________________________
658template <int N, typename... Args>
659GPUd() void DCAFitterN<N, Args...>::calcChi2Derivatives()
660{
661 //< calculate 1st and 2nd derivatives of wighted DCA (chi2) over track parameters X, see EQ.Chi2 in the ref
662 std::array<std::array<Vec3D, N>, N> covIDrDx; // tempory vectors of covI_j * dres_j/dx_i
663
664 // chi2 1st derivative
665 for (int i = N; i--;) {
666 auto& dchi1 = mDChi2Dx[i]; // DChi2/Dx_i = sum_j { res_j * covI_j * Dres_j/Dx_i }
667 dchi1 = 0;
668 for (int j = N; j--;) {
669 const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j
670 const auto& covI = mTrcEInv[mCurHyp][j]; // inverse cov matrix of track j
671 const auto& dr1 = mDResidDx[j][i]; // vector of j-th residuals 1st derivative over X param of track i
672 auto& cidr = covIDrDx[i][j]; // vector covI_j * dres_j/dx_i, save for 2nd derivative calculation
673 cidr[0] = covI.sxx * dr1[0];
674 cidr[1] = covI.syy * dr1[1] + covI.syz * dr1[2];
675 cidr[2] = covI.syz * dr1[1] + covI.szz * dr1[2];
676 // calculate res_i * covI_j * dres_j/dx_i
677 dchi1 += o2::math_utils::Dot(res, cidr);
678 }
679 }
680 // chi2 2nd derivative
681 for (int i = N; i--;) {
682 for (int j = i + 1; j--;) { // symmetric matrix
683 auto& dchi2 = mD2Chi2Dx2[i][j]; // D2Chi2/Dx_i/Dx_j = sum_k { Dres_k/Dx_j * covI_k * Dres_k/Dx_i + res_k * covI_k * D2res_k/Dx_i/Dx_j }
684 dchi2 = 0;
685 for (int k = N; k--;) {
686 const auto& dr1j = mDResidDx[k][j]; // vector of k-th residuals 1st derivative over X param of track j
687 const auto& cidrkj = covIDrDx[i][k]; // vector covI_k * dres_k/dx_i
688 dchi2 += o2::math_utils::Dot(dr1j, cidrkj);
689 if (k == j) {
690 const auto& res = mTrRes[mCurHyp][k]; // vector of residuals of track k
691 const auto& covI = mTrcEInv[mCurHyp][k]; // inverse cov matrix of track k
692 const auto& dr2ij = mD2ResidDx2[k][j]; // vector of k-th residuals 2nd derivative over X params of track j
693 dchi2 += res[0] * covI.sxx * dr2ij[0] + res[1] * (covI.syy * dr2ij[1] + covI.syz * dr2ij[2]) + res[2] * (covI.syz * dr2ij[1] + covI.szz * dr2ij[2]);
694 }
695 }
696 }
697 }
698}
699
700//__________________________________________________________________________
701template <int N, typename... Args>
702GPUd() void DCAFitterN<N, Args...>::calcChi2DerivativesNoErr()
703{
704 //< calculate 1st and 2nd derivatives of abs DCA (chi2) over track parameters X, see (6) in the ref
705 for (int i = N; i--;) {
706 auto& dchi1 = mDChi2Dx[i]; // DChi2/Dx_i = sum_j { res_j * Dres_j/Dx_i }
707 dchi1 = 0; // chi2 1st derivative
708 for (int j = N; j--;) {
709 const auto& res = mTrRes[mCurHyp][j]; // vector of residuals of track j
710 const auto& dr1 = mDResidDx[j][i]; // vector of j-th residuals 1st derivative over X param of track i
711 dchi1 += o2::math_utils::Dot(res, dr1);
712 if (i >= j) { // symmetrix matrix
713 // chi2 2nd derivative
714 auto& dchi2 = mD2Chi2Dx2[i][j]; // D2Chi2/Dx_i/Dx_j = sum_k { Dres_k/Dx_j * covI_k * Dres_k/Dx_i + res_k * covI_k * D2res_k/Dx_i/Dx_j }
715 dchi2 = o2::math_utils::Dot(mTrRes[mCurHyp][i], mD2ResidDx2[i][j]);
716 for (int k = N; k--;) {
717 dchi2 += o2::math_utils::Dot(mDResidDx[k][i], mDResidDx[k][j]);
718 }
719 }
720 }
721 }
722}
723
724//___________________________________________________________________
725template <int N, typename... Args>
726GPUd() void DCAFitterN<N, Args...>::calcPCA()
727{
728 // calculate point of closest approach for N prongs
729 mPCA[mCurHyp] = mTrCFVT[mCurHyp][N - 1] * mTrPos[mCurHyp][N - 1];
730 for (int i = N - 1; i--;) {
731 mPCA[mCurHyp] += mTrCFVT[mCurHyp][i] * mTrPos[mCurHyp][i];
732 }
733}
734
735//___________________________________________________________________
736template <int N, typename... Args>
737GPUd() bool DCAFitterN<N, Args...>::recalculatePCAWithErrors(int cand)
738{
739 // recalculate PCA as a cov-matrix weighted mean, even if absDCA method was used
740 if (isPropagateTracksToVertexDone(cand) && !propagateTracksToVertex(cand)) {
741 return false;
742 }
743 int saveCurHyp = mCurHyp;
744 mCurHyp = mOrder[cand];
745 if (mUseAbsDCA) {
746 for (int i = N; i--;) {
747 if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point
748 if (mLoggerBadCov.needToLog()) {
749#ifndef GPUCA_GPUCODE
750 printf("fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
751 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][i].asString().c_str());
752#else
753 printf("fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
754 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][i].getSigmaY2(), mCandTr[mCurHyp][i].getSigmaZ2(), mCandTr[mCurHyp][i].getSigmaZY());
755#endif
756 }
757 mFitStatus[mCurHyp] = FitStatus::FailInvCov;
758 if (mBadCovPolicy == Discard) {
759 return false;
760 } else if (mBadCovPolicy == OverrideAndFlag) {
761 mPropFailed[mCurHyp] = true;
762 } // otherwise, just use overridden errors w/o flagging
763 }
764 }
765 if (!calcPCACoefs()) {
766 mCurHyp = saveCurHyp;
767 return false;
768 }
769 }
770 auto oldPCA = mPCA[mOrder[cand]];
771 calcPCA();
772 mCurHyp = saveCurHyp;
773 return true;
774}
775
776//___________________________________________________________________
777template <int N, typename... Args>
778GPUd() void DCAFitterN<N, Args...>::calcPCANoErr()
779{
780 // calculate point of closest approach for N prongs w/o errors
781 auto& pca = mPCA[mCurHyp];
782 o2::math_utils::rotateZd(mTrPos[mCurHyp][N - 1][0], mTrPos[mCurHyp][N - 1][1], pca[0], pca[1], mTrAux[N - 1].s, mTrAux[N - 1].c);
783 // RRRR mTrAux[N-1].loc2glo(mTrPos[mCurHyp][N-1][0], mTrPos[mCurHyp][N-1][1], pca[0], pca[1] );
784 pca[2] = mTrPos[mCurHyp][N - 1][2];
785 for (int i = N - 1; i--;) {
786 double x, y;
787 o2::math_utils::rotateZd(mTrPos[mCurHyp][i][0], mTrPos[mCurHyp][i][1], x, y, mTrAux[i].s, mTrAux[i].c);
788 // RRRR mTrAux[i].loc2glo(mTrPos[mCurHyp][i][0], mTrPos[mCurHyp][i][1], x, y );
789 pca[0] += x;
790 pca[1] += y;
791 pca[2] += mTrPos[mCurHyp][i][2];
792 }
793 pca[0] *= NInv;
794 pca[1] *= NInv;
795 pca[2] *= NInv;
796}
797
798//___________________________________________________________________
799template <int N, typename... Args>
800GPUd() o2::math_utils::SMatrix<double, 3, 3, o2::math_utils::MatRepSym<double, 3>> DCAFitterN<N, Args...>::calcPCACovMatrix(int cand) const
801{
802 // calculate covariance matrix for the point of closest approach
803 MatSym3D covm;
804 int nAdded = 0;
805 for (int i = N; i--;) { // calculate sum of inverses
806 // MatSym3D covTr = o2::math_utils::Similarity(mUseAbsDCA ? getTrackRotMatrix(i) : mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand));
807 // RS by using Similarity(mTrCFVT[mOrder[cand]][i], getTrackCovMatrix(i, cand)) we underestimate the error, use simple rotation
808 MatSym3D covTr = o2::math_utils::Similarity(getTrackRotMatrix(i), getTrackCovMatrix(i, cand));
809 if (covTr.Invert()) {
810 covm += covTr;
811 nAdded++;
812 }
813 }
814 if (nAdded && covm.Invert()) {
815 return covm;
816 }
817 // correct way has failed, use simple sum
818 MatSym3D covmSum;
819 for (int i = N; i--;) {
820 MatSym3D covTr = o2::math_utils::Similarity(getTrackRotMatrix(i), getTrackCovMatrix(i, cand));
821 }
822 return covmSum;
823}
824
825//___________________________________________________________________
826template <int N, typename... Args>
827GPUd() void DCAFitterN<N, Args...>::calcTrackResiduals()
828{
829 // calculate residuals
830 Vec3D vtxLoc;
831 for (int i = N; i--;) {
832 mTrRes[mCurHyp][i] = mTrPos[mCurHyp][i];
833 vtxLoc = mPCA[mCurHyp];
834 o2::math_utils::rotateZInvd(vtxLoc[0], vtxLoc[1], vtxLoc[0], vtxLoc[1], mTrAux[i].s, mTrAux[i].c); // glo->loc
835 mTrRes[mCurHyp][i] -= vtxLoc;
836 }
837}
838
839//___________________________________________________________________
840template <int N, typename... Args>
841GPUdi() void DCAFitterN<N, Args...>::calcTrackDerivatives()
842{
843 // calculate track derivatives over X param
844 for (int i = N; i--;) {
845 mTrDer[mCurHyp][i].set(mCandTr[mCurHyp][i], mBz);
846 }
847}
848
849//___________________________________________________________________
850template <int N, typename... Args>
851GPUdi() double DCAFitterN<N, Args...>::calcChi2() const
852{
853 // calculate current chi2
854 double chi2 = 0;
855 for (int i = N; i--;) {
856 const auto& res = mTrRes[mCurHyp][i];
857 const auto& covI = mTrcEInv[mCurHyp][i];
858 chi2 += res[0] * res[0] * covI.sxx + res[1] * res[1] * covI.syy + res[2] * res[2] * covI.szz + 2. * res[1] * res[2] * covI.syz;
859 }
860 return chi2;
861}
862
863//___________________________________________________________________
864template <int N, typename... Args>
865GPUdi() double DCAFitterN<N, Args...>::calcChi2NoErr() const
866{
867 // calculate current chi2 of abs. distance minimization
868 double chi2 = 0;
869 for (int i = N; i--;) {
870 const auto& res = mTrRes[mCurHyp][i];
871 chi2 += res[0] * res[0] + res[1] * res[1] + res[2] * res[2];
872 }
873 return chi2;
874}
875
876//___________________________________________________________________
877template <int N, typename... Args>
878GPUd() bool DCAFitterN<N, Args...>::correctTracks(const VecND& corrX)
879{
880 // propagate tracks to updated X
881 for (int i = N; i--;) {
882 const auto& trDer = mTrDer[mCurHyp][i];
883 auto dx2h = 0.5 * corrX[i] * corrX[i];
884 mTrPos[mCurHyp][i][0] -= corrX[i];
885 mTrPos[mCurHyp][i][1] -= trDer.dydx * corrX[i] - dx2h * trDer.d2ydx2;
886 mTrPos[mCurHyp][i][2] -= trDer.dzdx * corrX[i] - dx2h * trDer.d2zdx2;
887 }
888 return true;
889}
890
891//___________________________________________________________________
892template <int N, typename... Args>
893GPUd() bool DCAFitterN<N, Args...>::propagateTracksToVertex(int icand)
894{
895 // propagate tracks to current vertex
896 int ord = mOrder[icand];
897 if (mTrPropDone[ord]) {
898 return true;
899 }
900
901 // need to refit taking as a seed already found vertex
902 if (mRefitWithMatCorr) {
903 int curHypSav = mCurHyp, curCrosIDAlt = mCrossIDAlt; // save
904 mCurHyp = ord;
905 mCrossIDAlt = -1; // disable alternative check
906 auto restore = [this, curHypSav, curCrosIDAlt]() { this->mCurHyp = curHypSav; this->mCrossIDAlt = curCrosIDAlt; };
907 if (!(mUseAbsDCA ? minimizeChi2NoErr() : minimizeChi2())) { // do final propagation
908 restore();
909 return false;
910 }
911 restore();
912 }
913
914 for (int i = N; i--;) {
915 if (mUseAbsDCA || mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
916 mCandTr[ord][i] = *mOrigTrPtr[i]; // fetch the track again, as mCandTr might have been propagated w/o errors or material corrections might be wrong
917 }
918 auto x = mTrAux[i].c * mPCA[ord][0] + mTrAux[i].s * mPCA[ord][1]; // X of PCA in the track frame
919 if (!propagateToX(mCandTr[ord][i], x)) {
920 return false;
921 }
922 }
923
924 mTrPropDone[ord] = true;
925 return true;
926}
927
928//___________________________________________________________________
929template <int N, typename... Args>
930GPUdi() o2::track::TrackPar DCAFitterN<N, Args...>::getTrackParamAtPCA(int i, int icand)
931{
932 // propagate tracks param only to current vertex (if not already done)
933 int ord = mOrder[icand];
934 o2::track::TrackPar trc(mCandTr[ord][i]);
935 if (!mTrPropDone[ord]) {
936 auto x = mTrAux[i].c * mPCA[ord][0] + mTrAux[i].s * mPCA[ord][1]; // X of PCA in the track frame
937 if (!propagateParamToX(trc, x)) {
938 trc.invalidate();
939 }
940 }
941 return trc;
942}
943
944//___________________________________________________________________
945template <int N, typename... Args>
946GPUdi() double DCAFitterN<N, Args...>::getAbsMax(const VecND& v)
947{
948 double mx = -1;
949 for (int i = N; i--;) {
950 auto vai = o2::gpu::GPUCommonMath::Abs(v[i]);
951 if (mx < vai) {
952 mx = vai;
953 }
954 }
955 return mx;
956}
957
958//___________________________________________________________________
959template <int N, typename... Args>
960GPUd() bool DCAFitterN<N, Args...>::minimizeChi2()
961{
962 // find best chi2 (weighted DCA) of N tracks in the vicinity of the seed PCA
963 for (int i = N; i--;) {
964 mCandTr[mCurHyp][i] = *mOrigTrPtr[i];
965 auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame
966 if (x < mMinXSeed) {
967 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
968 return false;
969 }
970 if (!propagateToX(mCandTr[mCurHyp][i], x)) {
971 return false;
972 }
973 setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
974 if (!mTrcEInv[mCurHyp][i].set(mCandTr[mCurHyp][i], XerrFactor)) { // prepare inverse cov.matrices at starting point
975 if (mLoggerBadCov.needToLog()) {
976#ifndef GPUCA_GPUCODE
977 printf("fitter %d: error (%ld muted): overrode invalid track covariance from %s\n",
978 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][i].asString().c_str());
979#else
980 printf("fitter %d: error (%ld muted): overrode invalid track covariance cyy:%e czz:%e cyz:%e\n",
981 mFitterID, mLoggerBadCov.evCount, mCandTr[mCurHyp][i].getSigmaY2(), mCandTr[mCurHyp][i].getSigmaZ2(), mCandTr[mCurHyp][i].getSigmaZY());
982#endif
983 }
984 mFitStatus[mCurHyp] = FitStatus::FailInvCov;
985 if (mBadCovPolicy == Discard) {
986 return false;
987 } else if (mBadCovPolicy == OverrideAndFlag) {
988 mPropFailed[mCurHyp] = true;
989 } // otherwise, just use overridden errors w/o flagging
990 }
991 }
992
993 if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference
994 mFitStatus[mCurHyp] = FitStatus::RejTrackRoughZ;
995 return false;
996 }
997
998 if (!calcPCACoefs()) { // prepare tracks contribution matrices to the global PCA
999 return false;
1000 }
1001 calcPCA(); // current PCA
1002 calcTrackResiduals(); // current track residuals
1003 float chi2Upd, chi2 = calcChi2();
1004 do {
1005 calcTrackDerivatives(); // current track derivatives (1st and 2nd)
1006 calcResidDerivatives(); // current residals derivatives (1st and 2nd)
1007 calcChi2Derivatives(); // current chi2 derivatives (1st and 2nd)
1008
1009 // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1
1010 if (!mD2Chi2Dx2.Invert()) {
1011 if (mLoggerBadInv.needToLog()) {
1012 printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount);
1013 }
1014 mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
1015 return false;
1016 }
1017 VecND dx = mD2Chi2Dx2 * mDChi2Dx;
1018 if (!correctTracks(dx)) {
1019 mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
1020 return false;
1021 }
1022 calcPCA(); // updated PCA
1023 if (mCrossIDAlt >= 0 && closerToAlternative()) {
1024 mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
1025 mAllowAltPreference = false;
1026 return false;
1027 }
1028 calcTrackResiduals(); // updated residuals
1029 chi2Upd = calcChi2(); // updated chi2
1030 if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
1031 chi2 = chi2Upd;
1032 mFitStatus[mCurHyp] = FitStatus::Converged;
1033 break; // converged
1034 }
1035 chi2 = chi2Upd;
1036 } while (++mNIters[mCurHyp] < mMaxIter);
1037 if (mNIters[mCurHyp] == mMaxIter) {
1038 mFitStatus[mCurHyp] = FitStatus::MaxIter;
1039 }
1040 //
1041 mChi2[mCurHyp] = chi2 * NInv;
1042 if (mChi2[mCurHyp] >= mMaxChi2) {
1043 mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1044 return false;
1045 }
1046 return true;
1047}
1048
1049//___________________________________________________________________
1050template <int N, typename... Args>
1051GPUd() bool DCAFitterN<N, Args...>::minimizeChi2NoErr()
1052{
1053 // find best chi2 (absolute DCA) of N tracks in the vicinity of the PCA seed
1054
1055 for (int i = N; i--;) {
1056 mCandTr[mCurHyp][i] = *mOrigTrPtr[i];
1057 auto x = mTrAux[i].c * mPCA[mCurHyp][0] + mTrAux[i].s * mPCA[mCurHyp][1]; // X of PCA in the track frame
1058 if (x < mMinXSeed) {
1059 mFitStatus[mCurHyp] = FitStatus::RejTrackX;
1060 return false;
1061 }
1062 if (!propagateParamToX(mCandTr[mCurHyp][i], x)) {
1063 return false;
1064 }
1065 setTrackPos(mTrPos[mCurHyp][i], mCandTr[mCurHyp][i]); // prepare positions
1066 }
1067 if (mMaxDZIni > 0 && !roughDZCut()) { // apply rough cut on tracks Z difference
1068 mFitStatus[mCurHyp] = FitStatus::RejTrackRoughZ;
1069 return false;
1070 }
1071
1072 calcPCANoErr(); // current PCA
1073 calcTrackResiduals(); // current track residuals
1074 float chi2Upd, chi2 = calcChi2NoErr();
1075 do {
1076 calcTrackDerivatives(); // current track derivatives (1st and 2nd)
1077 calcResidDerivativesNoErr(); // current residals derivatives (1st and 2nd)
1078 calcChi2DerivativesNoErr(); // current chi2 derivatives (1st and 2nd)
1079
1080 // do Newton-Rapson iteration with corrections = - dchi2/d{x0..xN} * [ d^2chi2/d{x0..xN}^2 ]^-1
1081 if (!mD2Chi2Dx2.Invert()) {
1082 if (mLoggerBadInv.needToLog()) {
1083 printf("fitter %d: error (%ld muted): Inversion failed\n", mFitterID, mLoggerBadCov.evCount);
1084 }
1085 mFitStatus[mCurHyp] = FitStatus::FailInv2ndDeriv;
1086 return false;
1087 }
1088 VecND dx = mD2Chi2Dx2 * mDChi2Dx;
1089 if (!correctTracks(dx)) {
1090 mFitStatus[mCurHyp] = FitStatus::FailCorrTracks;
1091 return false;
1092 }
1093 calcPCANoErr(); // updated PCA
1094 if (mCrossIDAlt >= 0 && closerToAlternative()) {
1095 mFitStatus[mCurHyp] = FitStatus::FailCloserAlt;
1096 mAllowAltPreference = false;
1097 return false;
1098 }
1099 calcTrackResiduals(); // updated residuals
1100 chi2Upd = calcChi2NoErr(); // updated chi2
1101 if (getAbsMax(dx) < mMinParamChange || chi2Upd > chi2 * mMinRelChi2Change) {
1102 chi2 = chi2Upd;
1103 mFitStatus[mCurHyp] = FitStatus::Converged;
1104 break; // converged
1105 }
1106 chi2 = chi2Upd;
1107 } while (++mNIters[mCurHyp] < mMaxIter);
1108 if (mNIters[mCurHyp] == mMaxIter) {
1109 mFitStatus[mCurHyp] = FitStatus::MaxIter;
1110 }
1111 //
1112 mChi2[mCurHyp] = chi2 * NInv;
1113 if (mChi2[mCurHyp] >= mMaxChi2) {
1114 mFitStatus[mCurHyp] = FitStatus::RejChi2Max;
1115 return false;
1116 }
1117 return true;
1118}
1119
1120//___________________________________________________________________
1121template <int N, typename... Args>
1122GPUd() bool DCAFitterN<N, Args...>::roughDZCut() const
1123{
1124 // apply rough cut on DZ between the tracks in the seed point
1125 bool accept = true;
1126 for (int i = N; accept && i--;) {
1127 for (int j = i; j--;) {
1128 if (o2::gpu::GPUCommonMath::Abs(mCandTr[mCurHyp][i].getZ() - mCandTr[mCurHyp][j].getZ()) > mMaxDZIni) {
1129 accept = false;
1130 break;
1131 }
1132 }
1133 }
1134 return accept;
1135}
1136
1137//___________________________________________________________________
1138template <int N, typename... Args>
1139GPUd() bool DCAFitterN<N, Args...>::closerToAlternative() const
1140{
1141 // check if the point current PCA point is closer to the seeding XY point being tested or to alternative see (if any)
1142 auto dxCur = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDCur], dyCur = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDCur];
1143 auto dxAlt = mPCA[mCurHyp][0] - mCrossings.xDCA[mCrossIDAlt], dyAlt = mPCA[mCurHyp][1] - mCrossings.yDCA[mCrossIDAlt];
1144 return dxCur * dxCur + dyCur * dyCur > dxAlt * dxAlt + dyAlt * dyAlt;
1145}
1146
1147//___________________________________________________________________
1148template <int N, typename... Args>
1149GPUd() void DCAFitterN<N, Args...>::print() const
1150{
1151#ifndef GPUCA_GPUCODE_DEVICE
1152 LOG(info) << N << "-prong vertex fitter in " << (mUseAbsDCA ? "abs." : "weighted") << " distance minimization mode, collinear tracks mode: " << (mIsCollinear ? "ON" : "OFF");
1153 LOG(info) << "Bz: " << mBz << " MaxIter: " << mMaxIter << " MaxChi2: " << mMaxChi2 << " MatCorrType: " << int(mMatCorr);
1154 LOG(info) << "Stopping condition: Max.param change < " << mMinParamChange << " Rel.Chi2 change > " << mMinRelChi2Change;
1155 LOG(info) << "Discard candidates for : Rvtx > " << getMaxR() << " DZ between tracks > " << mMaxDZIni;
1156 LOG(info) << "PropagateToPCA:" << mPropagateToPCA << " WeightedFinalPCA:" << mWeightedFinalPCA << " UsePropagator:" << mUsePropagator << " RefitWithMatCorr:" << mRefitWithMatCorr;
1157 std::string rep{};
1158 for (int i = 0; i < mCrossings.nDCA; i++) {
1159 rep += fmt::format("seed{}:{}/{} ", i, mTrPropDone[i], mPropFailed[i]);
1160 }
1161 LOG(info) << "Last call: NCand:" << mCurHyp << " from " << mCrossings.nDCA << " seeds, prop.done/failed: " << rep;
1162#else
1163 if (mUseAbsDCA) {
1164 printf("%d-prong vertex fitter in abs. distance minimization mode\n", N);
1165 } else {
1166 printf("%d-prong vertex fitter in weighted distance minimization mode\n", N);
1167 }
1168 printf("Bz: %1.f MaxIter: %3.d MaxChi2: %2.3f\n", mBz, mMaxIter, mMaxChi2);
1169 printf("Stopping condition: Max.param change < %2.3f Rel.Chi2 change > %2.3f\n", mMinParamChange, mMinRelChi2Change);
1170 printf("Discard candidates for : Rvtx > %2.3f DZ between tracks > %2.3f\n", getMaxR(), mMaxDZIni);
1171#endif
1172}
1173
1174//___________________________________________________________________
1175template <int N, typename... Args>
1176GPUd() o2::track::TrackParCov DCAFitterN<N, Args...>::createParentTrackParCov(int cand, bool sectorAlpha) const
1177{
1178 const auto& trP = getTrack(0, cand);
1179 const auto& trN = getTrack(1, cand);
1180 std::array<float, 21> covV = {0.};
1181 std::array<float, 3> pvecV = {0.};
1182 int q = 0;
1183 for (int it = 0; it < N; it++) {
1184 const auto& trc = getTrack(it, cand);
1185 std::array<float, 3> pvecT = {0.};
1186 std::array<float, 21> covT = {0.};
1187 trc.getPxPyPzGlo(pvecT);
1188 trc.getCovXYZPxPyPzGlo(covT);
1189 constexpr int MomInd[6] = {9, 13, 14, 18, 19, 20}; // cov matrix elements for momentum component
1190 for (int i = 0; i < 6; i++) {
1191 covV[MomInd[i]] += covT[MomInd[i]];
1192 }
1193 for (int i = 0; i < 3; i++) {
1194 pvecV[i] += pvecT[i];
1195 }
1196 q += trc.getCharge();
1197 }
1198 auto covVtxV = calcPCACovMatrix(cand);
1199 covV[0] = covVtxV(0, 0);
1200 covV[1] = covVtxV(1, 0);
1201 covV[2] = covVtxV(1, 1);
1202 covV[3] = covVtxV(2, 0);
1203 covV[4] = covVtxV(2, 1);
1204 covV[5] = covVtxV(2, 2);
1205 return o2::track::TrackParCov(getPCACandidatePos(cand), pvecV, covV, q, sectorAlpha);
1206}
1207
1208//___________________________________________________________________
1209template <int N, typename... Args>
1210GPUd() o2::track::TrackPar DCAFitterN<N, Args...>::createParentTrackPar(int cand, bool sectorAlpha) const
1211{
1212 const auto& trP = getTrack(0, cand);
1213 const auto& trN = getTrack(1, cand);
1214 const auto& wvtx = getPCACandidate(cand);
1215 std::array<float, 3> pvecV = {0.};
1216 int q = 0;
1217 for (int it = 0; it < N; it++) {
1218 const auto& trc = getTrack(it, cand);
1219 std::array<float, 3> pvecT = {0.};
1220 trc.getPxPyPzGlo(pvecT);
1221 for (int i = 0; i < 3; i++) {
1222 pvecV[i] += pvecT[i];
1223 }
1224 q += trc.getCharge();
1225 }
1226 const std::array<float, 3> vertex = {(float)wvtx[0], (float)wvtx[1], (float)wvtx[2]};
1227 return o2::track::TrackPar(vertex, pvecV, q, sectorAlpha);
1228}
1229
1230//___________________________________________________________________
1231template <int N, typename... Args>
1232GPUdi() bool DCAFitterN<N, Args...>::propagateParamToX(o2::track::TrackPar& t, float x)
1233{
1234 bool res = true;
1235 if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
1236#ifndef GPUCA_GPUCODE
1237 res = o2::base::Propagator::Instance()->PropagateToXBxByBz(t, x, mMaxSnp, mMaxStep, mMatCorr);
1238#endif
1239 } else {
1240 res = t.propagateParamTo(x, mBz);
1241 }
1242 if (!res) {
1243 mFitStatus[mCurHyp] = FitStatus::FailProp;
1244 mPropFailed[mCurHyp] = true;
1245 if (mLoggerBadProp.needToLog()) {
1246#ifndef GPUCA_GPUCODE
1247 printf("fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str());
1248#else
1249 printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount);
1250#endif
1251 }
1252 }
1253 return res;
1254}
1255
1256//___________________________________________________________________
1257template <int N, typename... Args>
1258GPUdi() bool DCAFitterN<N, Args...>::propagateToX(o2::track::TrackParCov& t, float x)
1259{
1260 bool res = true;
1261 if (mUsePropagator || mMatCorr != o2::base::Propagator::MatCorrType::USEMatCorrNONE) {
1262#ifndef GPUCA_GPUCODE
1263 res = o2::base::Propagator::Instance()->PropagateToXBxByBz(t, x, mMaxSnp, mMaxStep, mMatCorr);
1264#endif
1265 } else {
1266 res = t.propagateTo(x, mBz);
1267 }
1268 if (!res) {
1269 mFitStatus[mCurHyp] = FitStatus::FailProp;
1270 mPropFailed[mCurHyp] = true;
1271 if (mLoggerBadProp.needToLog()) {
1272#ifndef GPUCA_GPUCODE
1273 printf("fitter %d: error (%ld muted): propagation failed for %s\n", mFitterID, mLoggerBadProp.evCount, t.asString().c_str());
1274#else
1275 printf("fitter %d: error (%ld muted): propagation failed\n", mFitterID, mLoggerBadProp.evCount);
1276#endif
1277 }
1278 }
1279 return res;
1280}
1281
1284
1285namespace device
1286{
1287template <typename Fitter>
1288void print(const int nBlocks, const int nThreads, Fitter& ft);
1289
1290template <typename Fitter, class... Tr>
1291int process(const int nBlocks, const int nThreads, Fitter&, Tr&... args);
1292
1293template <class Fitter, class... Tr>
1294void processBulk(const int nBlocks, const int nThreads, const int nBatches, std::vector<Fitter>& fitters, std::vector<int>& results, std::vector<Tr>&... args);
1295} // namespace device
1296
1297} // namespace vertexing
1298} // namespace o2
1299#endif // _ALICEO2_DCA_FITTERN_
Helper classes for helical tracks manipulations.
Base track model for the Barrel, params only, w/o covariance.
uint64_t vertex
Definition RawEventData.h:9
void print() const
int32_t i
uint32_t j
Definition RawData.h:0
uint32_t res
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
o2::track::TrackParCov TrackParCov
Definition Recon.h:39
GPUd() value_type estimateLTFast(o2 static GPUd() float estimateLTIncrement(const o2 PropagatorImpl * Instance(bool uninitialized=false)
Definition Propagator.h:143
const Track & getTrack(int i, int cand=0) const
create parent track param with errors for decay vertex
Definition DCAFitterN.h:194
GPUd() const auto getPCACandidatePos(int cand=0) const
return position of quality-ordered candidate in the internal structures
Definition DCAFitterN.h:160
const Track * getOrigTrackPtr(int i) const
Definition DCAFitterN.h:224
GPUdi() FitStatus getFitStatus(int cand=0) const noexcept
return number of iterations during minimization (no check for its validity)
Definition DCAFitterN.h:226
int class Tr const T & t
Definition DCAFitterN.h:329
void setBadCovPolicy(BadCovPolicy v)
Definition DCAFitterN.h:364
GPUd() bool calcPCACoefs()
GPUdi() size_t getCallID() const
Definition DCAFitterN.h:275
GPUdi() void setFitterID(int i)
Definition DCAFitterN.h:274
int cand
track X-param at V0 candidate (no check for the candidate validity)
Definition DCAFitterN.h:301
GPUdi() int getNIterations(int cand=0) const
Definition DCAFitterN.h:229
GPUd() const Vec3D &getPCACandidate(int cand=0) const
< return PCA candidate, by default best on is provided (no check for the index validity)
Definition DCAFitterN.h:159
GPUdi() void clearLogThrottlers()
Definition DCAFitterN.h:357
int getCandidatePosition(int cand=0) const
return Chi2 at PCA candidate (no check for its validity)
Definition DCAFitterN.h:167
static constexpr int getNProngs()
Definition DCAFitterN.h:149
bool isPropagationFailure(int cand=0) const
Definition DCAFitterN.h:180
DCAFitterN(float bz, bool useAbsDCA, bool prop2DCA)
Definition DCAFitterN.h:152
GPUdi() void setPropagateToPCA(bool v
Track & getTrack(int i, int cand=0)
Definition DCAFitterN.h:184
std::array< float, 6 > calcPCACovMatrixFlat(int cand=0) const
Definition DCAFitterN.h:218
GPUdi() static void setTrackPos(Vec3D &pnt
float getChi2AtPCACandidate(int cand=0) const
Definition DCAFitterN.h:170
GPUd() bool propagateTracksToVertex(int cand=0)
check if propagation of tracks to candidate vertex was done
BadCovPolicy getBadCovPolicy() const
Definition DCAFitterN.h:365
int class Tr const T const Tr & args
Definition DCAFitterN.h:330
GLdouble n
Definition glcorearb.h:1982
GLint GLenum GLint x
Definition glcorearb.h:403
const GLfloat * m
Definition glcorearb.h:4066
const GLdouble * v
Definition glcorearb.h:832
GLenum array
Definition glcorearb.h:4274
GLdouble f
Definition glcorearb.h:310
GLint y
Definition glcorearb.h:270
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLboolean r
Definition glcorearb.h:1233
void set(T *h, size_t v)
Definition PageParser.h:107
constexpr float Almost0
std::tuple< double, double > rotateZd(double xL, double yL, double snAlp, double csAlp)
Definition Utils.h:167
std::tuple< double, double > rotateZInvd(double xG, double yG, double snAlp, double csAlp)
Definition Utils.h:147
SMatrix< T, D1, D1, MatRepSym< T, D1 > > Similarity(const SMatrix< T, D1, D2, R > &lhs, const SMatrix< T, D2, D2, MatRepSym< T, D2 > > &rhs)
Definition Cartesian.h:263
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Definition Cartesian.h:257
TrackParCovF TrackParCov
Definition Track.h:33
TrackParF TrackPar
Definition Track.h:29
void processBulk(const int nBlocks, const int nThreads, const int nBatches, std::vector< Fitter > &fitters, std::vector< int > &results, std::vector< Tr > &... args)
GPUd() int DCAFitterN< N
ROOT::Math::SVector< double, 3 > Vec3D
GPUdi() void DCAFitterN< N
Definition DCAFitterN.h:930
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Defining DataPointCompositeObject explicitly as copiable.
GPUdi() bool needToLog()
Definition DCAFitterN.h:79
GPUd() bool set(const o2
Definition DCAFitterN.h:37
GPUdDefault() TrackCovI()=default
GPUd() TrackDeriv(const o2
Definition DCAFitterN.h:63
GPUd() void set(const o2
Definition DCAFitterN.h:64
GPUdDefault() TrackDeriv()=default
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"
vec clear()