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