Project
Loading...
Searching...
No Matches
Propagator.cxx
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
13#include "GPUCommonLogger.h"
14#include "GPUCommonConstants.h"
15#include "GPUCommonMath.h"
17#include "MathUtils/Utils.h"
20
21using namespace o2::base;
22using namespace o2::gpu;
23
24#if !defined(GPUCA_GPUCODE)
25#include "Field/MagFieldFast.h" // Don't use this on the GPU
26#endif
27
28#if !defined(GPUCA_STANDALONE) && !defined(GPUCA_GPUCODE)
29#include "Field/MagneticField.h"
33#include <FairRunAna.h> // eventually will get rid of it
34#include <TGeoGlobalMagField.h>
35
36template <typename value_T>
38{
39 if (uninitialized) {
40 return;
41 }
43 updateField();
44}
45
46//____________________________________________________________
47template <typename value_T>
49{
50 if (!mField) {
51 mField = static_cast<o2::field::MagneticField*>(TGeoGlobalMagField::Instance()->GetField());
52 if (!mField) {
53 LOG(warning) << "No Magnetic Field in TGeoGlobalMagField, checking legacy FairRunAna";
54 mField = dynamic_cast<o2::field::MagneticField*>(FairRunAna::Instance()->GetField());
55 }
56 if (!mField) {
57 LOG(fatal) << "Magnetic field is not initialized!";
58 }
59 if (!mField->getFastField() && mField->fastFieldExists()) {
60 mField->AllowFastField(true);
61 mFieldFast = mField->getFastField();
62 }
63 }
64 const value_type xyz[3] = {0.};
65 if (mFieldFast) {
66 mFieldFast->GetBz(xyz, mNominalBz);
67 } else {
68 mNominalBz = mField->GetBz(xyz[0], xyz[1], xyz[2]);
69 }
70}
71
72//____________________________________________________________
73template <typename value_T>
74int PropagatorImpl<value_T>::initFieldFromGRP(const std::string grpFileName, bool verbose)
75{
77 if (verbose) {
78 LOG(info) << "Loading field from GRP of " << grpFileName;
79 }
80 const auto grp = o2::parameters::GRPObject::loadFrom(grpFileName);
81 if (!grp) {
82 return -1;
83 }
84 if (verbose) {
85 grp->print();
86 }
87
88 return initFieldFromGRP(grp);
89}
90
91//____________________________________________________________
92template <typename value_T>
94{
96
97 if (TGeoGlobalMagField::Instance()->IsLocked()) {
98 if (TGeoGlobalMagField::Instance()->GetField()->TestBit(o2::field::MagneticField::kOverrideGRP)) {
99 LOG(warning) << "ExpertMode!!! GRP information will be ignored";
100 LOG(warning) << "ExpertMode!!! Running with the externally locked B field";
101 return 0;
102 } else {
103 LOG(info) << "Destroying existing B field instance";
104 delete TGeoGlobalMagField::Instance();
105 }
106 }
108 TGeoGlobalMagField::Instance()->SetField(fld);
109 TGeoGlobalMagField::Instance()->Lock();
110 if (verbose) {
111 LOG(info) << "Running with the B field constructed out of GRP";
112 LOG(info) << "Access field via TGeoGlobalMagField::Instance()->Field(xyz,bxyz) or via";
113 LOG(info) << "auto o2field = static_cast<o2::field::MagneticField*>( TGeoGlobalMagField::Instance()->GetField() )";
114 }
115 return 0;
116}
117
118//____________________________________________________________
119template <typename value_T>
121{
123
124 if (TGeoGlobalMagField::Instance()->IsLocked()) {
125 if (TGeoGlobalMagField::Instance()->GetField()->TestBit(o2::field::MagneticField::kOverrideGRP)) {
126 LOG(warning) << "ExpertMode!!! GRP information will be ignored";
127 LOG(warning) << "ExpertMode!!! Running with the externally locked B field";
128 return 0;
129 } else {
130 LOG(info) << "Destroying existing B field instance";
131 delete TGeoGlobalMagField::Instance();
132 }
133 }
135 TGeoGlobalMagField::Instance()->SetField(fld);
136 TGeoGlobalMagField::Instance()->Lock();
137 if (verbose) {
138 LOG(info) << "Running with the B field constructed out of GRP";
139 LOG(info) << "Access field via TGeoGlobalMagField::Instance()->Field(xyz,bxyz) or via";
140 LOG(info) << "auto o2field = static_cast<o2::field::MagneticField*>( TGeoGlobalMagField::Instance()->GetField() )";
141 }
142 return 0;
143}
144
145#elif !defined(GPUCA_GPUCODE)
146template <typename value_T>
148{
149} // empty dummy constructor for standalone benchmark
150#endif
151
152//_______________________________________________________________________
153template <typename value_T>
154GPUd() bool PropagatorImpl<value_T>::PropagateToXBxByBz(TrackParCov_t& track, value_type xToGo, value_type maxSnp, value_type maxStep,
155 PropagatorImpl<value_T>::MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const
156{
157 //----------------------------------------------------------------
158 //
159 // Propagates the track to the plane X=xk (cm)
160 // taking into account all the three components of the magnetic field
161 // and correcting for the crossed material.
162 //
163 // maxStep - maximal step for propagation
164 // tofInfo - optional container for track length and PID-dependent TOF integration
165 //
166 // matCorr - material correction type, it is up to the user to make sure the pointer is attached (if LUT is requested)
167 //----------------------------------------------------------------
168 auto dx = xToGo - track.getX();
169 int dir = dx > 0.f ? 1 : -1;
170 if (!signCorr) {
171 signCorr = -dir; // sign of eloss correction is not imposed
172 }
173
174 std::array<value_type, 3> b{};
175 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
176 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
177 if (dir < 0) {
178 step = -step;
179 }
180 auto x = track.getX() + step;
181 auto xyz0 = track.getXYZGlo();
182 getFieldXYZ(xyz0, &b[0]);
183
184 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr, this]() {
185 bool res = true;
186 if (matCorr != MatCorrType::USEMatCorrNONE) {
187 auto xyz1 = track.getXYZGlo();
188 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
189 if (!track.correctForMaterial(mb.meanX2X0, mb.getXRho(signCorr))) {
190 res = false;
191 }
192 if (tofInfo) {
193 tofInfo->addStep(mb.length, track.getQ2P2()); // fill L,ToF info using already calculated step length
194 tofInfo->addX2X0(mb.meanX2X0);
195 tofInfo->addXRho(mb.getXRho(signCorr));
196 }
197 } else if (tofInfo) { // if tofInfo filling was requested w/o material correction, we need to calculate the step lenght
198 auto xyz1 = track.getXYZGlo();
199 math_utils::Vector3D<value_type> stepV(xyz1.X() - xyz0.X(), xyz1.Y() - xyz0.Y(), xyz1.Z() - xyz0.Z());
200 tofInfo->addStep(stepV.R(), track.getQ2P2());
201 }
202 return res;
203 };
204
205 if (!track.propagateTo(x, b)) {
206 return false;
207 }
208 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
209 correct();
210 return false;
211 }
212 if (!correct()) {
213 return false;
214 }
215 dx = xToGo - track.getX();
216 }
217 track.setX(xToGo);
218 return true;
219}
220
221//_______________________________________________________________________
222template <typename value_T>
223GPUd() bool PropagatorImpl<value_T>::PropagateToXBxByBz(TrackPar_t& track, value_type xToGo, value_type maxSnp, value_type maxStep,
224 PropagatorImpl<value_T>::MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const
225{
226 //----------------------------------------------------------------
227 //
228 // Propagates the track params to the plane X=xk (cm), NO error evaluation
229 // taking into account all the three components of the magnetic field
230 // and optionally correcting for the e.loss crossed material.
231 //
232 // maxStep - maximal step for propagation
233 // tofInfo - optional container for track length and PID-dependent TOF integration
234 //
235 // matCorr - material correction type, it is up to the user to make sure the pointer is attached (if LUT is requested)
236 //----------------------------------------------------------------
237 auto dx = xToGo - track.getX();
238 int dir = dx > 0.f ? 1 : -1;
239 if (!signCorr) {
240 signCorr = -dir; // sign of eloss correction is not imposed
241 }
242
243 std::array<value_type, 3> b{};
244 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
245 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
246 if (dir < 0) {
247 step = -step;
248 }
249 auto x = track.getX() + step;
250 auto xyz0 = track.getXYZGlo();
251 getFieldXYZ(xyz0, &b[0]);
252
253 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr, this]() {
254 bool res = true;
255 if (matCorr != MatCorrType::USEMatCorrNONE) {
256 auto xyz1 = track.getXYZGlo();
257 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
258 if (!track.correctForELoss(((signCorr < 0) ? -mb.length : mb.length) * mb.meanRho)) {
259 res = false;
260 }
261 if (tofInfo) {
262 tofInfo->addStep(mb.length, track.getQ2P2()); // fill L,ToF info using already calculated step length
263 tofInfo->addX2X0(mb.meanX2X0);
264 tofInfo->addXRho(mb.getXRho(signCorr));
265 }
266 } else if (tofInfo) { // if tofInfo filling was requested w/o material correction, we need to calculate the step lenght
267 auto xyz1 = track.getXYZGlo();
268 math_utils::Vector3D<value_type> stepV(xyz1.X() - xyz0.X(), xyz1.Y() - xyz0.Y(), xyz1.Z() - xyz0.Z());
269 tofInfo->addStep(stepV.R(), track.getQ2P2());
270 }
271 return res;
272 };
273
274 if (!track.propagateParamTo(x, b)) {
275 return false;
276 }
277 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
278 correct();
279 return false;
280 }
281 if (!correct()) {
282 return false;
283 }
284 dx = xToGo - track.getX();
285 }
286 track.setX(xToGo);
287 return true;
288}
289
290//_______________________________________________________________________
291template <typename value_T>
292GPUd() bool PropagatorImpl<value_T>::propagateToX(TrackParCov_t& track, value_type xToGo, value_type bZ, value_type maxSnp, value_type maxStep,
293 PropagatorImpl<value_T>::MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const
294{
295 //----------------------------------------------------------------
296 //
297 // Propagates the track to the plane X=xk (cm)
298 // taking into account all the three components of the magnetic field
299 // and correcting for the crossed material.
300 //
301 // maxStep - maximal step for propagation
302 // tofInfo - optional container for track length and PID-dependent TOF integration
303 //
304 // matCorr - material correction type, it is up to the user to make sure the pointer is attached (if LUT is requested)
305 //----------------------------------------------------------------
306 auto dx = xToGo - track.getX();
307 int dir = dx > 0.f ? 1 : -1;
308 if (!signCorr) {
309 signCorr = -dir; // sign of eloss correction is not imposed
310 }
311
312 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
313 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
314 if (dir < 0) {
315 step = -step;
316 }
317 auto x = track.getX() + step;
318 auto xyz0 = track.getXYZGlo();
319 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr, this]() {
320 bool res = true;
321 if (matCorr != MatCorrType::USEMatCorrNONE) {
322 auto xyz1 = track.getXYZGlo();
323 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
324 if (!track.correctForMaterial(mb.meanX2X0, mb.getXRho(signCorr))) {
325 res = false;
326 }
327 if (tofInfo) {
328 tofInfo->addStep(mb.length, track.getQ2P2()); // fill L,ToF info using already calculated step length
329 tofInfo->addX2X0(mb.meanX2X0);
330 tofInfo->addXRho(mb.getXRho(signCorr));
331 }
332 } else if (tofInfo) { // if tofInfo filling was requested w/o material correction, we need to calculate the step lenght
333 auto xyz1 = track.getXYZGlo();
334 math_utils::Vector3D<value_type> stepV(xyz1.X() - xyz0.X(), xyz1.Y() - xyz0.Y(), xyz1.Z() - xyz0.Z());
335 tofInfo->addStep(stepV.R(), track.getQ2P2());
336 }
337 return res;
338 };
339 if (!track.propagateTo(x, bZ)) {
340 return false;
341 }
342 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
343 correct();
344 return false;
345 }
346 if (!correct()) {
347 return false;
348 }
349 dx = xToGo - track.getX();
350 }
351 track.setX(xToGo);
352 return true;
353}
354
355//_______________________________________________________________________
356template <typename value_T>
357GPUd() bool PropagatorImpl<value_T>::propagateToX(TrackPar_t& track, value_type xToGo, value_type bZ, value_type maxSnp, value_type maxStep,
358 PropagatorImpl<value_T>::MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const
359{
360 //----------------------------------------------------------------
361 //
362 // Propagates the track parameters only to the plane X=xk (cm)
363 // taking into account all the three components of the magnetic field
364 // and correcting for the crossed material.
365 //
366 // maxStep - maximal step for propagation
367 // tofInfo - optional container for track length and PID-dependent TOF integration
368 //
369 // matCorr - material correction type, it is up to the user to make sure the pointer is attached (if LUT is requested)
370 //----------------------------------------------------------------
371 auto dx = xToGo - track.getX();
372 int dir = dx > 0.f ? 1 : -1;
373 if (!signCorr) {
374 signCorr = -dir; // sign of eloss correction is not imposed
375 }
376
377 while (math_utils::detail::abs<value_type>(dx) > Epsilon) {
378 auto step = math_utils::detail::min<value_type>(math_utils::detail::abs<value_type>(dx), maxStep);
379 if (dir < 0) {
380 step = -step;
381 }
382 auto x = track.getX() + step;
383 auto xyz0 = track.getXYZGlo();
384
385 auto correct = [&track, &xyz0, tofInfo, matCorr, signCorr, this]() {
386 bool res = true;
387 if (matCorr != MatCorrType::USEMatCorrNONE) {
388 auto xyz1 = track.getXYZGlo();
389 auto mb = this->getMatBudget(matCorr, xyz0, xyz1);
390 if (!track.correctForELoss(mb.getXRho(signCorr))) {
391 res = false;
392 }
393 if (tofInfo) {
394 tofInfo->addStep(mb.length, track.getQ2P2()); // fill L,ToF info using already calculated step length
395 tofInfo->addX2X0(mb.meanX2X0);
396 tofInfo->addXRho(mb.getXRho(signCorr));
397 }
398 } else if (tofInfo) { // if tofInfo filling was requested w/o material correction, we need to calculate the step lenght
399 auto xyz1 = track.getXYZGlo();
400 math_utils::Vector3D<value_type> stepV(xyz1.X() - xyz0.X(), xyz1.Y() - xyz0.Y(), xyz1.Z() - xyz0.Z());
401 tofInfo->addStep(stepV.R(), track.getQ2P2());
402 }
403 return res;
404 };
405
406 if (!track.propagateParamTo(x, bZ)) {
407 return false;
408 }
409 if (maxSnp > 0 && math_utils::detail::abs<value_type>(track.getSnp()) >= maxSnp) {
410 correct();
411 return false;
412 }
413 if (!correct()) {
414 return false;
415 }
416 dx = xToGo - track.getX();
417 }
418 track.setX(xToGo);
419 return true;
420}
421
422//_______________________________________________________________________
423template <typename value_T>
424template <typename track_T>
425GPUd() bool PropagatorImpl<value_T>::propagateToR(track_T& track, value_type r, bool bzOnly, value_type maxSnp, value_type maxStep,
426 MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const
427{
428 const value_T MaxPhiLoc = math_utils::detail::asin<value_T>(maxSnp), MaxPhiLocSafe = 0.95 * MaxPhiLoc;
429 auto bz = getNominalBz();
430 if (math_utils::detail::abs(bz) > constants::math::Almost0) {
431 o2::track::TrackAuxPar traux(track, bz);
433 value_type r0 = math_utils::detail::sqrt<value_T>(track.getX() * track.getX() + track.getY() * track.getY());
434 value_type dr = (r - r0);
435 value_type rTmp = r - (math_utils::detail::abs<value_T>(dr) > 1. ? (dr > 0 ? 0.5 : -0.5) : 0.5 * dr); // 1st propagate a few mm short of the targer R
436 crad.rC = rTmp;
437 crad.c = crad.cc = 1.f;
438 crad.s = crad.ss = crad.cs = 0.f;
440 cross.circlesCrossInfo(crad, traux, 0.);
441 if (cross.nDCA < 1) {
442 return false;
443 }
444 double phiCross[2] = {}, dphi[2] = {};
445 auto curv = track.getCurvature(bz);
446 bool clockwise = curv < 0; // q+ in B+ or q- in B- goes clockwise
447 auto phiLoc = math_utils::detail::asin<double>(track.getSnp());
448 auto phi0 = phiLoc + track.getAlpha();
450 for (int i = 0; i < cross.nDCA; i++) {
451 // track pT direction angle at crossing points:
452 // == angle of the tangential to track circle at the crossing point X,Y
453 // == normal to the radial vector from the track circle center {X-cX, Y-cY}
454 // i.e. the angle of the vector {Y-cY, -(X-cx)}
455 auto normX = double(cross.yDCA[i]) - double(traux.yC), normY = -(double(cross.xDCA[i]) - double(traux.xC));
456 if (!clockwise) {
457 normX = -normX;
458 normY = -normY;
459 }
460 phiCross[i] = math_utils::detail::atan2<double>(normY, normX);
462 dphi[i] = phiCross[i] - phi0;
463 if (dphi[i] > o2::constants::math::PI) {
465 } else if (dphi[i] < -o2::constants::math::PI) {
467 }
468 }
469 int sel = cross.nDCA == 1 ? 0 : (clockwise ? (dphi[0] < dphi[1] ? 0 : 1) : (dphi[1] < dphi[0] ? 0 : 1));
470 auto deltaPhi = dphi[sel];
471
472 while (1) {
473 auto phiLocFin = phiLoc + deltaPhi;
474 // case1
475 if (math_utils::detail::abs<value_type>(phiLocFin) < MaxPhiLocSafe) { // just 1 step propagation
476 auto deltaX = (math_utils::detail::sin<double>(phiLocFin) - track.getSnp()) / track.getCurvature(bz);
477 if (!track.propagateTo(track.getX() + deltaX, bz)) {
478 return false;
479 }
480 break;
481 }
482 if (math_utils::detail::abs<value_type>(deltaPhi) < (2 * MaxPhiLocSafe)) { // still can go in 1 step with one extra rotation
483 auto rot = phiLoc + 0.5 * deltaPhi;
484 if (!track.rotate(track.getAlpha() + rot)) {
485 return false;
486 }
487 phiLoc -= rot;
488 continue; // should be ok for the case 1 now.
489 }
490
491 auto rot = phiLoc + (deltaPhi > 0 ? MaxPhiLocSafe : -MaxPhiLocSafe);
492 if (!track.rotate(track.getAlpha() + rot)) {
493 return false;
494 }
495 phiLoc -= rot; // = +- MaxPhiLocSafe
496
497 // propagate to phiLoc = +-MaxPhiLocSafe
498 auto tgtPhiLoc = deltaPhi > 0 ? MaxPhiLocSafe : -MaxPhiLocSafe;
499 auto deltaX = (math_utils::detail::sin<double>(tgtPhiLoc) - track.getSnp()) / track.getCurvature(bz);
500 if (!track.propagateTo(track.getX() + deltaX, bz)) {
501 return false;
502 }
503 deltaPhi -= tgtPhiLoc - phiLoc;
504 phiLoc = deltaPhi > 0 ? MaxPhiLocSafe : -MaxPhiLocSafe;
505 continue; // should be of for the case 1 now.
506 }
507 bz = getBz(math_utils::Point3D<value_type>{value_type(cross.xDCA[sel]), value_type(cross.yDCA[sel]), value_type(track.getZ())});
508 }
509 // do final step till target R, also covers Bz = 0;
510 value_type xfin;
511 if (!track.getXatLabR(r, xfin, bz)) {
512 return false;
513 }
514 return propagateToX(track, xfin, bzOnly, maxSnp, maxStep, matCorr, tofInfo, signCorr);
515}
516
517//_______________________________________________________________________
518template <typename value_T>
519template <typename track_T>
520GPUd() bool PropagatorImpl<value_T>::propagateToAlphaX(track_T& track, value_type alpha, value_type x, bool bzOnly, value_type maxSnp, value_type maxStep, int minSteps,
521 MatCorrType matCorr, track::TrackLTIntegral* tofInfo, int signCorr) const
522{
523 // propagate to alpha,X, if needed in a few steps
524 auto snp = track.getSnpAt(alpha, x, getNominalBz());
525 // apply safety factor 0.9 for crude rotation estimate
526 if (math_utils::detail::abs<value_type>(snp) < maxSnp * 0.9 && track.rotate(alpha)) {
527 auto dx = math_utils::detail::abs<value_type>(x - track.getX());
528 if (dx < Epsilon) {
529 return true;
530 }
531 return propagateTo(track, x, bzOnly, maxSnp, math_utils::detail::min<value_type>(dx / minSteps, maxStep), matCorr, tofInfo, signCorr);
532 }
533 return false;
534 /*
535 // try to go in a few steps with intermediate rotations
536
537
538 auto alphaTrg = alpha;
539 math_utils::detail::bringToPMPi<value_type>(alphaTrg);
540 auto alpCurr = track.getAlpha();
541 math_utils::detail::bringToPMPi<value_type>(alpCurr);
542 int nsteps = minSteps > 2 ? minSteps : 2;
543 auto dalp = math_utils::detail::deltaPhiSmall<value_type>(alpCurr, alphaTrg) / nsteps; // effective (alpha - alphaCurr)/nsteps
544 auto xtmp = (track.getX() + x) / nsteps;
545 return track.rotate(alpCurr + dalp) && propagateTo(track, xtmp, bzOnly, maxSnp, maxStep, matCorr, tofInfo, signCorr) &&
546 track.rotate(alpha) && propagateTo(track, x, bzOnly, maxSnp, maxStep, matCorr, tofInfo, signCorr);
547 */
548}
549
550//_______________________________________________________________________
551template <typename value_T>
552GPUd() bool PropagatorImpl<value_T>::propagateToDCA(const o2::dataformats::VertexBase& vtx, TrackParCov_t& track, value_type bZ,
553 value_type maxStep, PropagatorImpl<value_type>::MatCorrType matCorr,
554 o2::dataformats::DCA* dca, track::TrackLTIntegral* tofInfo,
555 int signCorr, value_type maxD) const
556{
557 // propagate track to DCA to the vertex
558 value_type sn, cs, alp = track.getAlpha();
559 math_utils::detail::sincos<value_type>(alp, sn, cs);
560 value_type x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
561 value_type xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
562 x -= xv;
563 y -= yv;
564 // Estimate the impact parameter neglecting the track curvature
565 value_type d = math_utils::detail::abs<value_type>(x * snp - y * csp);
566 if (d > maxD) {
567 return false;
568 }
569 value_type crv = track.getCurvature(bZ);
570 value_type tgfv = -(crv * x - snp) / (crv * y + csp);
571 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
572 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
573 cs = (math_utils::detail::abs<value_type>(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
574
575 x = xv * cs + yv * sn;
576 yv = -xv * sn + yv * cs;
577 xv = x;
578
579 auto tmpT(track); // operate on the copy to recover after the failure
580 alp += math_utils::detail::asin<value_type>(sn);
581 if (!tmpT.rotate(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
582#ifndef GPUCA_ALIGPUCODE
583 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx << " | Track is: " << tmpT.asString();
584#elif !defined(GPUCA_NO_FMT)
585 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx;
586#endif
587 return false;
588 }
589 track = tmpT;
590 if (dca) {
591 math_utils::detail::sincos<value_type>(alp, sn, cs);
592 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
593 dca->set(track.getY() - yv, track.getZ() - zv,
594 track.getSigmaY2() + s2ylocvtx, track.getSigmaZY(), track.getSigmaZ2() + vtx.getSigmaZ2());
595 }
596 return true;
597}
598
599//_______________________________________________________________________
600template <typename value_T>
601GPUd() bool PropagatorImpl<value_T>::propagateToDCABxByBz(const o2::dataformats::VertexBase& vtx, TrackParCov_t& track,
602 value_type maxStep, PropagatorImpl<value_type>::MatCorrType matCorr,
603 o2::dataformats::DCA* dca, track::TrackLTIntegral* tofInfo,
604 int signCorr, value_type maxD) const
605{
606 // propagate track to DCA to the vertex
607 value_type sn, cs, alp = track.getAlpha();
608 math_utils::detail::sincos<value_type>(alp, sn, cs);
609 value_type x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
610 value_type xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
611 x -= xv;
612 y -= yv;
613 // Estimate the impact parameter neglecting the track curvature
614 value_type d = math_utils::detail::abs<value_type>(x * snp - y * csp);
615 if (d > maxD) {
616 return false;
617 }
618 value_type crv = track.getCurvature(mNominalBz);
619 value_type tgfv = -(crv * x - snp) / (crv * y + csp);
620 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
621 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
622 cs = (math_utils::detail::abs<value_type>(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
623
624 x = xv * cs + yv * sn;
625 yv = -xv * sn + yv * cs;
626 xv = x;
627
628 auto tmpT(track); // operate on the copy to recover after the failure
629 alp += math_utils::detail::asin<value_type>(sn);
630 if (!tmpT.rotate(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
631#ifndef GPUCA_ALIGPUCODE
632 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx << " | Track is: " << tmpT.asString();
633#elif !defined(GPUCA_NO_FMT)
634 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx;
635#endif
636 return false;
637 }
638 track = tmpT;
639 if (dca) {
640 math_utils::detail::sincos<value_type>(alp, sn, cs);
641 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
642 dca->set(track.getY() - yv, track.getZ() - zv,
643 track.getSigmaY2() + s2ylocvtx, track.getSigmaZY(), track.getSigmaZ2() + vtx.getSigmaZ2());
644 }
645 return true;
646}
647
648//_______________________________________________________________________
649template <typename value_T>
650GPUd() bool PropagatorImpl<value_T>::propagateToDCA(const math_utils::Point3D<value_type>& vtx, TrackPar_t& track, value_type bZ,
651 value_type maxStep, PropagatorImpl<value_T>::MatCorrType matCorr,
652 std::array<value_type, 2>* dca, track::TrackLTIntegral* tofInfo,
653 int signCorr, value_type maxD) const
654{
655 // propagate track to DCA to the vertex
656 value_type sn, cs, alp = track.getAlpha();
657 math_utils::detail::sincos<value_type>(alp, sn, cs);
658 value_type x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
659 value_type xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
660 x -= xv;
661 y -= yv;
662 // Estimate the impact parameter neglecting the track curvature
663 value_type d = math_utils::detail::abs<value_type>(x * snp - y * csp);
664 if (d > maxD) {
665 return false;
666 }
667 value_type crv = track.getCurvature(bZ);
668 value_type tgfv = -(crv * x - snp) / (crv * y + csp);
669 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
670 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
671 cs = (math_utils::detail::abs<value_type>(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
672
673 x = xv * cs + yv * sn;
674 yv = -xv * sn + yv * cs;
675 xv = x;
676
677 auto tmpT(track); // operate on the copy to recover after the failure
678 alp += math_utils::detail::asin<value_type>(sn);
679 if (!tmpT.rotateParam(alp) || !propagateToX(tmpT, xv, bZ, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
680#ifndef GPUCA_ALIGPUCODE
681 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex "
682 << vtx.X() << ' ' << vtx.Y() << ' ' << vtx.Z() << " | Track is: " << tmpT.asString();
683#else
684 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex " << vtx.X() << ' ' << vtx.Y() << ' ' << vtx.Z();
685#endif
686 return false;
687 }
688 track = tmpT;
689 if (dca) {
690 (*dca)[0] = track.getY() - yv;
691 (*dca)[1] = track.getZ() - zv;
692 }
693 return true;
694}
695
696//_______________________________________________________________________
697template <typename value_T>
698GPUd() bool PropagatorImpl<value_T>::propagateToDCABxByBz(const math_utils::Point3D<value_type>& vtx, TrackPar_t& track,
699 value_type maxStep, PropagatorImpl<value_T>::MatCorrType matCorr,
700 std::array<value_type, 2>* dca, track::TrackLTIntegral* tofInfo,
701 int signCorr, value_type maxD) const
702{
703 // propagate track to DCA to the vertex
704 value_type sn, cs, alp = track.getAlpha();
705 math_utils::detail::sincos<value_type>(alp, sn, cs);
706 value_type x = track.getX(), y = track.getY(), snp = track.getSnp(), csp = math_utils::detail::sqrt<value_type>((1.f - snp) * (1.f + snp));
707 value_type xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
708 x -= xv;
709 y -= yv;
710 // Estimate the impact parameter neglecting the track curvature
711 value_type d = math_utils::detail::abs<value_type>(x * snp - y * csp);
712 if (d > maxD) {
713 return false;
714 }
715 value_type crv = track.getCurvature(mNominalBz);
716 value_type tgfv = -(crv * x - snp) / (crv * y + csp);
717 sn = tgfv / math_utils::detail::sqrt<value_type>(1.f + tgfv * tgfv);
718 cs = math_utils::detail::sqrt<value_type>((1. - sn) * (1. + sn));
719 cs = (math_utils::detail::abs<value_type>(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
720
721 x = xv * cs + yv * sn;
722 yv = -xv * sn + yv * cs;
723 xv = x;
724
725 auto tmpT(track); // operate on the copy to recover after the failure
726 alp += math_utils::detail::asin<value_type>(sn);
727 if (!tmpT.rotateParam(alp) || !PropagateToXBxByBz(tmpT, xv, 0.85, maxStep, matCorr, tofInfo, signCorr)) {
728#ifndef GPUCA_ALIGPUCODE
729 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex "
730 << vtx.X() << ' ' << vtx.Y() << ' ' << vtx.Z() << " | Track is: " << tmpT.asString();
731#else
732 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex " << vtx.X() << ' ' << vtx.Y() << ' ' << vtx.Z();
733#endif
734 return false;
735 }
736 track = tmpT;
737 if (dca) {
738 (*dca)[0] = track.getY() - yv;
739 (*dca)[1] = track.getZ() - zv;
740 }
741 return true;
742}
743
744//____________________________________________________________
745template <typename value_T>
746GPUd() float PropagatorImpl<value_T>::estimateLTIncrement(const o2::track::TrackParametrization<value_type>& trc,
747 const o2::math_utils::Point3D<value_type>& posStart,
748 const o2::math_utils::Point3D<value_type>& posEnd) const
749{
750 // estimate helical step increment between 2 point
751 float dX = posEnd.X() - posStart.X(), dY = posEnd.Y() - posStart.Y(), dZ = posEnd.Z() - posStart.Z(), d2XY = dX * dX + dY * dY;
752 if (getNominalBz() != 0) { // circular arc = 2*R*asin(dXY/2R)
753 float b[3];
754 o2::math_utils::Point3D<float> posAv(0.5 * (posEnd.X() + posStart.X()), 0.5 * (posEnd.Y() + posStart.Y()), 0.5 * (posEnd.Z() + posStart.Z()));
755 getFieldXYZ(posAv, b);
756 float curvH = math_utils::detail::abs<value_type>(0.5f * trc.getCurvature(b[2])), asArg = curvH * math_utils::detail::sqrt<value_type>(d2XY);
757 if (curvH > 0.f) {
758 d2XY = asArg < 1.f ? math_utils::detail::asin<value_type>(asArg) / curvH : o2::constants::math::PIHalf / curvH;
759 d2XY *= d2XY;
760 }
761 }
762 return math_utils::detail::sqrt<value_type>(d2XY + dZ * dZ);
763}
764
765//____________________________________________________________
766template <typename value_T>
767GPUd() value_T PropagatorImpl<value_T>::estimateLTFast(o2::track::TrackLTIntegral& lt, const o2::track::TrackParametrization<value_type>& trc) const
768{
769 value_T xdca = 0., ydca = 0., length = 0.; // , zdca = 0. // zdca might be used in future
771 constexpr float TinyF = 1e-9;
772 auto straigh_line_approx = [&]() {
773 auto csp2 = (1.f - trc.getSnp()) * (1.f + trc.getSnp());
774 if (csp2 > TinyF) {
775 auto csp = math_utils::detail::sqrt<value_type>(csp2);
776 auto tgp = trc.getSnp() / csp, f = trc.getX() * tgp - trc.getY();
777 xdca = tgp * f * csp2;
778 ydca = -f * csp2;
779 auto dx = xdca - trc.getX(), dy = ydca - trc.getY(), dz = dx * trc.getTgl() / csp;
780 return math_utils::detail::sqrt<value_type>(dx * dx + dy * dy + dz * dz);
781 } else { // track is parallel to Y axis
782 xdca = trc.getX(); // ydca = 0
783 return math_utils::detail::abs<value_type>(trc.getY() * math_utils::detail::sqrt<value_type>(1. + trc.getTgl() * trc.getTgl())); // distance from the current point to DCA
784 }
785 };
786 trc.getCircleParamsLoc(mNominalBz, c);
787 if (c.rC != 0.) { // helix
788 auto distC = math_utils::detail::sqrt<value_type>(c.getCenterD2()); // distance from the circle center to origin
789 if (distC > 1.e-3) {
790 auto nrm = (distC - c.rC) / distC;
791 xdca = nrm * c.xC; // coordinates of the DCA to 0,0 in the local frame
792 ydca = nrm * c.yC;
793 auto v0x = trc.getX() - c.xC, v0y = trc.getY() - c.yC, v1x = xdca - c.xC, v1y = ydca - c.yC;
794 auto angcos = (v0x * v1x + v0y * v1y) / (c.rC * c.rC);
795 if (math_utils::detail::abs<value_type>(angcos) < 1.f) {
796 auto ang = math_utils::detail::acos<value_type>(angcos);
797 if ((trc.getSign() > 0.f) == (mNominalBz > 0.f)) {
798 ang = -ang; // we need signeg angle
799 c.rC = -c.rC; // we need signed curvature for zdca
800 }
801 // zdca = trc.getZ() + (trc.getSign() > 0. ? c.rC : -c.rC) * trc.getTgl() * ang;
802 length = math_utils::detail::abs<value_type>(c.rC * ang * math_utils::detail::sqrt<value_type>(1. + trc.getTgl() * trc.getTgl()));
803 } else { // calculation of the arc length between the position and DCA makes no sense
804 length = straigh_line_approx();
805 }
806 } else { // track with circle center at the origin, and LT makes no sense, take direct distance
807 xdca = trc.getX();
808 ydca = trc.getY();
809 }
810 } else { // straight line
811 length = straigh_line_approx();
812 }
813 // since we assume the track or its parent comes from the beam-line or decay, add XY(?) distance to it
814 value_T dcaT = math_utils::detail::sqrt<value_type>(xdca * xdca + ydca * ydca);
815 length += dcaT;
816 lt.addStep(length, trc.getQ2P2());
817 return dcaT;
818}
819
820//____________________________________________________________
821template <typename value_T>
822GPUd() MatBudget PropagatorImpl<value_T>::getMatBudget(PropagatorImpl<value_type>::MatCorrType corrType, const math_utils::Point3D<value_type>& p0, const math_utils::Point3D<value_type>& p1) const
823{
824#if !defined(GPUCA_STANDALONE) && !defined(GPUCA_GPUCODE)
825 if (corrType == MatCorrType::USEMatCorrTGeo) {
827 }
828 if (!mMatLUT) {
829 if (mTGeoFallBackAllowed) {
831 } else {
832 throw std::runtime_error("requested MatLUT is absent and fall-back to TGeo is disabled");
833 }
834 }
835#endif
836 return mMatLUT->getMatBudget(p0.X(), p0.Y(), p0.Z(), p1.X(), p1.Y(), p1.Z());
837}
838
839template <typename value_T>
840template <typename T>
841GPUd() void PropagatorImpl<value_T>::getFieldXYZImpl(const math_utils::Point3D<T> xyz, T* bxyz) const
842{
843 if (mGPUField) {
844#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
845 const auto* f = &GPUCA_CONSMEM.param.polynomialField; // Access directly from constant memory on GPU (copied here to avoid complicated header dependencies)
846#else
847 const auto* f = mGPUField;
848#endif
849 float bxyzF[3] = {};
850 f->GetField(xyz.X(), xyz.Y(), xyz.Z(), bxyzF);
851 // copy and convert
852 constexpr value_type kCLight1 = 1. / o2::gpu::gpu_common_constants::kCLight;
853 for (uint i = 0; i < 3; ++i) {
854 bxyz[i] = static_cast<value_type>(bxyzF[i]) * kCLight1;
855 }
856 } else {
857#ifndef GPUCA_GPUCODE
858 if (mFieldFast) {
859 mFieldFast->Field(xyz, bxyz); // Must not call the host-only function in GPU compilation
860 } else {
861#ifdef GPUCA_STANDALONE
862 LOG(fatal) << "Normal field cannot be used in standalone benchmark";
863#else
864 mField->field(xyz, bxyz);
865#endif
866 }
867#endif
868 }
869}
870
871template <typename value_T>
872template <typename T>
873GPUd() T PropagatorImpl<value_T>::getBzImpl(const math_utils::Point3D<T> xyz) const
874{
875 T bz = 0;
876 if (mGPUField) {
877#if defined(GPUCA_GPUCODE_DEVICE) && defined(GPUCA_HAS_GLOBAL_SYMBOL_CONSTANT_MEM)
878 const auto* f = &GPUCA_CONSMEM.param.polynomialField; // Access directly from constant memory on GPU (copied here to avoid complicated header dependencies)
879#else
880 const auto* f = mGPUField;
881#endif
882 constexpr value_type kCLight1 = 1. / o2::gpu::gpu_common_constants::kCLight;
883 bz = f->GetFieldBz(xyz.X(), xyz.Y(), xyz.Z()) * kCLight1;
884 } else {
885#ifndef GPUCA_GPUCODE
886 if (mFieldFast) {
887 mFieldFast->GetBz(xyz, bz); // Must not call the host-only function in GPU compilation
888 } else {
889#ifdef GPUCA_STANDALONE
890 LOG(fatal) << "Normal field cannot be used in standalone benchmark";
891#else
892 bz = mField->GetBz(xyz.X(), xyz.Y(), xyz.Z());
893#endif
894 }
895#endif
896 }
897 return bz;
898}
899
900template <typename value_T>
901GPUd() void PropagatorImpl<value_T>::getFieldXYZ(const math_utils::Point3D<float> xyz, float* bxyz) const
902{
903 getFieldXYZImpl<float>(xyz, bxyz);
904}
905
906template <typename value_T>
907GPUd() void PropagatorImpl<value_T>::getFieldXYZ(const math_utils::Point3D<double> xyz, double* bxyz) const
908{
909 getFieldXYZImpl<double>(xyz, bxyz);
910}
911
912template <typename value_T>
913GPUd() float PropagatorImpl<value_T>::getBz(const math_utils::Point3D<float> xyz) const
914{
915 return getBzImpl<float>(xyz);
916}
917
918template <typename value_T>
919GPUd() double PropagatorImpl<value_T>::getBz(const math_utils::Point3D<double> xyz) const
920{
921 return getBzImpl<double>(xyz);
922}
923
924namespace o2::base
925{
926#if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code.
927template class PropagatorImpl<float>;
928template bool GPUdni() PropagatorImpl<float>::propagateToAlphaX<PropagatorImpl<float>::TrackPar_t>(PropagatorImpl<float>::TrackPar_t&, float, float, bool, float, float, int, PropagatorImpl<float>::MatCorrType matCorr, track::TrackLTIntegral*, int) const;
929template bool GPUdni() PropagatorImpl<float>::propagateToAlphaX<PropagatorImpl<float>::TrackParCov_t>(PropagatorImpl<float>::TrackParCov_t&, float, float, bool, float, float, int, PropagatorImpl<float>::MatCorrType matCorr, track::TrackLTIntegral*, int) const;
930template bool GPUdni() PropagatorImpl<float>::propagateToR<PropagatorImpl<float>::TrackPar_t>(PropagatorImpl<float>::TrackPar_t&, float, bool, float, float, PropagatorImpl<float>::MatCorrType matCorr, track::TrackLTIntegral*, int) const;
931template bool GPUdni() PropagatorImpl<float>::propagateToR<PropagatorImpl<float>::TrackParCov_t>(PropagatorImpl<float>::TrackParCov_t&, float, bool, float, float, PropagatorImpl<float>::MatCorrType matCorr, track::TrackLTIntegral*, int) const;
932#endif
933#ifndef GPUCA_GPUCODE
934template class PropagatorImpl<double>;
937#endif
938} // namespace o2::base
General auxilliary methods.
Helper classes for helical tracks manipulations.
Definition of the GeometryManager class.
int32_t i
#define GPUdni()
#define GPUCA_CONSMEM
Header of the General Run Parameters object for B field values.
Header of the General Run Parameters object.
constexpr int p1()
constexpr to accelerate the coordinates changing
Definition of the fast magnetic field parametrization MagFieldFast.
Definition of the MagF class.
uint32_t res
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
std::ostringstream debug
static o2::base::MatBudget meanMaterialBudget(float x0, float y0, float z0, float x1, float y1, float z1)
GPUd() bool propagateToDCA(const o2 GPUd() bool propagateToDCABxByBz(const o2 GPUd() bool propagateToDCA(const o2 GPUd() bool propagateToDCABxByBz(const o2 PropagatorImpl(PropagatorImpl const &)=delete
static int initFieldFromGRP(const o2::parameters::GRPMagField *grp, bool verbose=false)
static MagneticField * createFieldMap(float l3Current=-30000., float diCurrent=-6000., Int_t convention=0, Bool_t uniform=kFALSE, float beamenergy=7000, const Char_t *btype="pp", const std::string path=std::string(gSystem->Getenv("VMCWORKDIR"))+std::string("/Common/maps/mfchebKGI_sym.root"))
void AllowFastField(bool v=true)
allow fast field param
o2::units::Current_t getDipoleCurrent() const
Definition GRPMagField.h:43
bool getFieldUniformity() const
Definition GRPMagField.h:49
o2::units::Current_t getL3Current() const
getters/setters for magnets currents
Definition GRPMagField.h:37
o2::units::Current_t getDipoleCurrent() const
Definition GRPObject.h:80
bool getFieldUniformity() const
Definition GRPObject.h:81
static GRPObject * loadFrom(const std::string &grpFileName="")
o2::units::Current_t getL3Current() const
getters/setters for magnets currents
Definition GRPObject.h:79
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
GLenum array
Definition glcorearb.h:4274
GLdouble f
Definition glcorearb.h:310
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLuint GLsizei GLsizei * length
Definition glcorearb.h:790
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLboolean r
Definition glcorearb.h:1233
constexpr float Almost0
constexpr float TwoPI
constexpr float PI
constexpr float PIHalf
constexpr float Almost1
GPUd() const expr uint32_t MultivariatePolynomialHelper< Dim
std::array< int, 24 > p0
value_T step
Definition TrackUtils.h:42
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Defining DataPointCompositeObject explicitly as copiable.
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"