Project
Loading...
Searching...
No Matches
TrackExtrap.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
16
18
19#include <TGeoGlobalMagField.h>
20#include <TGeoManager.h>
21#include <TGeoMaterial.h>
22#include <TGeoNode.h>
23#include <TGeoShape.h>
24#include <TMath.h>
25
26#include "Framework/Logger.h"
27
29
30namespace o2
31{
32namespace mch
33{
34
35bool TrackExtrap::sExtrapV2 = false;
36double TrackExtrap::sSimpleBValue = 0.;
37bool TrackExtrap::sFieldON = false;
38std::size_t TrackExtrap::sNCallExtrapToZCov = 0;
39std::size_t TrackExtrap::sNCallField = 0;
40
41//__________________________________________________________________________
43{
46 const double x[3] = {50., 50., SSimpleBPosition};
47 double b[3] = {0., 0., 0.};
48 TGeoGlobalMagField::Instance()->Field(x, b);
49 sSimpleBValue = b[0];
50 sFieldON = (TMath::Abs(sSimpleBValue) > 1.e-10) ? true : false;
51 LOG(info) << "Track extrapolation with magnetic field " << (sFieldON ? "ON" : "OFF");
52}
53
54//__________________________________________________________________________
56{
61 const double correctionFactor = 1.1; // impact parameter is 10% underestimated
62 if (bendingMomentum == 0.) {
63 return 1.e10;
64 }
65 return correctionFactor * (-0.0003 * sSimpleBValue * SSimpleBLength * SSimpleBPosition / bendingMomentum);
66}
67
68//__________________________________________________________________________
70{
75 const double correctionFactor = 1.1; // bending momentum is 10% underestimated
76 if (impactParam == 0.) {
77 return 1.e10;
78 }
79 if (sFieldON) {
80 return correctionFactor * (-0.0003 * sSimpleBValue * SSimpleBLength * SSimpleBPosition / impactParam);
81 } else {
82 return SMostProbBendingMomentum;
83 }
84}
85
86//__________________________________________________________________________
87void TrackExtrap::linearExtrapToZ(TrackParam& trackParam, double zEnd)
88{
91
92 if (trackParam.getZ() == zEnd) {
93 return; // nothing to be done if same z
94 }
95
96 // Compute track parameters
97 double dZ = zEnd - trackParam.getZ();
98 trackParam.setNonBendingCoor(trackParam.getNonBendingCoor() + trackParam.getNonBendingSlope() * dZ);
99 trackParam.setBendingCoor(trackParam.getBendingCoor() + trackParam.getBendingSlope() * dZ);
100 trackParam.setZ(zEnd);
101}
102
103//__________________________________________________________________________
104void TrackExtrap::linearExtrapToZCov(TrackParam& trackParam, double zEnd, bool updatePropagator)
105{
108
109 if (trackParam.getZ() == zEnd) {
110 return; // nothing to be done if same z
111 }
112
113 // No need to propagate the covariance matrix if it does not exist
114 if (!trackParam.hasCovariances()) {
115 LOG(warning) << "Covariance matrix does not exist";
116 // Extrapolate linearly track parameters to "zEnd"
117 linearExtrapToZ(trackParam, zEnd);
118 return;
119 }
120
121 // Compute track parameters
122 double dZ = zEnd - trackParam.getZ();
123 trackParam.setNonBendingCoor(trackParam.getNonBendingCoor() + trackParam.getNonBendingSlope() * dZ);
124 trackParam.setBendingCoor(trackParam.getBendingCoor() + trackParam.getBendingSlope() * dZ);
125 trackParam.setZ(zEnd);
126
127 // Calculate the jacobian related to the track parameters linear extrapolation to "zEnd"
128 TMatrixD jacob(5, 5);
129 jacob.UnitMatrix();
130 jacob(0, 1) = dZ;
131 jacob(2, 3) = dZ;
132
133 // Extrapolate track parameter covariances to "zEnd"
134 TMatrixD tmp(trackParam.getCovariances(), TMatrixD::kMultTranspose, jacob);
135 TMatrixD tmp2(jacob, TMatrixD::kMult, tmp);
136 trackParam.setCovariances(tmp2);
137
138 // Update the propagator if required
139 if (updatePropagator) {
140 trackParam.updatePropagator(jacob);
141 }
142}
143
144//__________________________________________________________________________
145bool TrackExtrap::extrapToZ(TrackParam& trackParam, double zEnd)
146{
149 if (!sFieldON) {
150 linearExtrapToZ(trackParam, zEnd);
151 return true;
152 } else if (sExtrapV2) {
153 return extrapToZRungekuttaV2(trackParam, zEnd);
154 } else {
155 return extrapToZRungekutta(trackParam, zEnd);
156 }
157}
158
159//__________________________________________________________________________
160bool TrackExtrap::extrapToZCov(TrackParam& trackParam, double zEnd, bool updatePropagator)
161{
164
165 ++sNCallExtrapToZCov;
166
167 if (trackParam.getZ() == zEnd) {
168 return true; // nothing to be done if same z
169 }
170
171 if (!sFieldON) { // linear extrapolation if no magnetic field
172 linearExtrapToZCov(trackParam, zEnd, updatePropagator);
173 return true;
174 }
175
176 // No need to propagate the covariance matrix if it does not exist
177 if (!trackParam.hasCovariances()) {
178 LOG(warning) << "Covariance matrix does not exist";
179 // Extrapolate track parameters to "zEnd"
180 return extrapToZ(trackParam, zEnd);
181 }
182
183 // Save the actual track parameters
184 TrackParam trackParamSave(trackParam);
185 TMatrixD paramSave(trackParamSave.getParameters());
186 double zBegin = trackParamSave.getZ();
187
188 // Get reference to the parameter covariance matrix
189 const TMatrixD& kParamCov = trackParam.getCovariances();
190
191 // Extrapolate track parameters to "zEnd"
192 // Do not update the covariance matrix if the extrapolation failed
193 if (!extrapToZ(trackParam, zEnd)) {
194 return false;
195 }
196
197 // Get reference to the extrapolated parameters
198 const TMatrixD& extrapParam = trackParam.getParameters();
199
200 // Calculate the jacobian related to the track parameters extrapolation to "zEnd"
201 TMatrixD jacob(5, 5);
202 jacob.Zero();
203 TMatrixD dParam(5, 1);
204 double direction[5] = {-1., -1., 1., 1., -1.};
205 for (int i = 0; i < 5; i++) {
206 // Skip jacobian calculation for parameters with no associated error
207 if (kParamCov(i, i) <= 0.) {
208 continue;
209 }
210
211 // Small variation of parameter i only
212 for (int j = 0; j < 5; j++) {
213 if (j == i) {
214 dParam(j, 0) = TMath::Sqrt(kParamCov(i, i));
215 dParam(j, 0) *= TMath::Sign(1., direction[j] * paramSave(j, 0)); // variation always in the same direction
216 } else {
217 dParam(j, 0) = 0.;
218 }
219 }
220
221 // Set new parameters
222 trackParamSave.setParameters(paramSave);
223 trackParamSave.addParameters(dParam);
224 trackParamSave.setZ(zBegin);
225
226 // Extrapolate new track parameters to "zEnd"
227 if (!extrapToZ(trackParamSave, zEnd)) {
228 LOG(warning) << "Bad covariance matrix";
229 return false;
230 }
231
232 // Calculate the jacobian
233 TMatrixD jacobji(trackParamSave.getParameters(), TMatrixD::kMinus, extrapParam);
234 jacobji *= 1. / dParam(i, 0);
235 jacob.SetSub(0, i, jacobji);
236 }
237
238 // Extrapolate track parameter covariances to "zEnd"
239 TMatrixD tmp(kParamCov, TMatrixD::kMultTranspose, jacob);
240 TMatrixD tmp2(jacob, TMatrixD::kMult, tmp);
241 trackParam.setCovariances(tmp2);
242
243 // Update the propagator if required
244 if (updatePropagator) {
245 trackParam.updatePropagator(jacob);
246 }
247
248 return true;
249}
250
251//__________________________________________________________________________
253{
256
257 if (trackParam.getZ() == SMIDZ) {
258 return true; // nothing to be done
259 }
260
261 // check the track position with respect to the muon filter (spectro z<0)
262 if (trackParam.getZ() < SMuonFilterZBeg) {
263 LOG(warning) << "The track already passed the beginning of the muon filter";
264 return false;
265 }
266
267 // propagate through the muon filter and add MCS effets
268 if (!extrapToZCov(trackParam, SMuonFilterZEnd)) {
269 return false;
270 }
271 addMCSEffect(trackParam, -SMuonFilterThickness, SMuonFilterX0);
272
273 // propagate to the first MID chamber
274 return extrapToZCov(trackParam, SMIDZ);
275}
276
277//__________________________________________________________________________
278bool TrackExtrap::extrapToVertex(TrackParam& trackParam, double xVtx, double yVtx, double zVtx,
279 double errXVtx, double errYVtx, bool correctForMCS, bool correctForEnergyLoss,
280 double xUpstream, double yUpstream, std::optional<double> zUpstream)
281{
291
292 if (trackParam.getZ() == zVtx) {
293 return true; // nothing to be done if already at vertex
294 }
295
296 // Check the vertex position with respect to the absorber (spectro z<0)
297 if (zVtx < SAbsZBeg) {
298 if (zVtx < SAbsZEnd) {
299 LOG(warning) << "Ending Z (" << zVtx << ") downstream the front absorber (zAbsorberEnd = " << SAbsZEnd << ")";
300 return false;
301 } else {
302 LOG(warning) << "Ending Z (" << zVtx << ") inside the front absorber (" << SAbsZBeg << ", " << SAbsZEnd << ")";
303 return false;
304 }
305 }
306
307 // check the upstream track position with respect to the absorber if provided and used (spectro z<0)
308 // zUpstream must be >= SAbsZBeg with 100 µm tolerance to account for numerical precision
309 if (!correctForMCS && zUpstream && *zUpstream < SAbsZBeg - 0.01) {
310 if (*zUpstream < SAbsZEnd) {
311 LOG(warning) << "Upstream Z (" << *zUpstream << ") downstream the front absorber (zAbsorberEnd = " << SAbsZEnd << ")";
312 return false;
313 } else {
314 LOG(warning) << "Upstream Z (" << *zUpstream << ") inside the front absorber (" << SAbsZBeg << ", " << SAbsZEnd << ")";
315 return false;
316 }
317 }
318
319 // Check the track position with respect to the vertex and the absorber (spectro z<0)
320 if (trackParam.getZ() > SAbsZEnd) {
321 if (trackParam.getZ() > zVtx) {
322 LOG(warning) << "Starting Z (" << trackParam.getZ() << ") upstream the vertex (zVtx = " << zVtx << ")";
323 return false;
324 } else if (trackParam.getZ() > SAbsZBeg) {
325 LOG(warning) << "Starting Z (" << trackParam.getZ() << ") upstream the front absorber (zAbsorberBegin = " << SAbsZBeg << ")";
326 return false;
327 } else {
328 LOG(warning) << "Starting Z (" << trackParam.getZ() << ") inside the front absorber (" << SAbsZBeg << ", " << SAbsZEnd << ")";
329 return false;
330 }
331 }
332
333 // Extrapolate track parameters (and covariances if any) to the end of the absorber
334 if ((trackParam.hasCovariances() && !extrapToZCov(trackParam, SAbsZEnd)) ||
335 (!trackParam.hasCovariances() && !extrapToZ(trackParam, SAbsZEnd))) {
336 return false;
337 }
338
339 // Get absorber correction parameters assuming linear propagation in absorber
340 double trackXYZOut[3] = {trackParam.getNonBendingCoor(), trackParam.getBendingCoor(), trackParam.getZ()};
341 double trackXYZIn[3] = {0., 0., 0.};
342 if (correctForMCS) { // assume linear propagation to the vertex
343 trackXYZIn[2] = SAbsZBeg;
344 trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
345 trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
346 } else if (zUpstream) { // or linear propagation to the upstream track position
347 trackXYZIn[2] = SAbsZBeg;
348 trackXYZIn[0] = trackXYZOut[0] + (xUpstream - trackXYZOut[0]) / (*zUpstream - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
349 trackXYZIn[1] = trackXYZOut[1] + (yUpstream - trackXYZOut[1]) / (*zUpstream - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
350 } else { // or standard propagation without vertex constraint
351 TrackParam trackParamIn(trackParam);
352 if (!extrapToZ(trackParamIn, SAbsZBeg)) {
353 return false;
354 }
355 trackXYZIn[0] = trackParamIn.getNonBendingCoor();
356 trackXYZIn[1] = trackParamIn.getBendingCoor();
357 trackXYZIn[2] = trackParamIn.getZ();
358 }
359 double pTot = trackParam.p();
360 double pathLength(0.), f0(0.), f1(0.), f2(0.), meanRho(0.), totalELoss(0.), sigmaELoss2(0.);
361 if (!getAbsorberCorrectionParam(trackXYZIn, trackXYZOut, pTot, pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2)) {
362 return false;
363 }
364
365 // Compute track parameters and covariances at vertex according to correctForMCS and correctForEnergyLoss flags
366 if (correctForMCS) {
367 if (correctForEnergyLoss) {
368 // Correct for multiple scattering and energy loss
369 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
370 if (!correctMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, trackXYZIn[2], pathLength, f0, f1, f2)) {
371 return false;
372 }
373 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
374 } else {
375 // Correct for multiple scattering
376 if (!correctMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, trackXYZIn[2], pathLength, f0, f1, f2)) {
377 return false;
378 }
379 }
380 } else {
381 if (correctForEnergyLoss) {
382 // Correct for energy loss add multiple scattering dispersion in covariance matrix
383 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
384 addMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
385 if (!extrapToZCov(trackParam, trackXYZIn[2])) {
386 return false;
387 }
388 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
389 if (!extrapToZCov(trackParam, zVtx)) {
390 return false;
391 }
392 } else {
393 // add multiple scattering dispersion in covariance matrix
394 addMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
395 if (!extrapToZCov(trackParam, zVtx)) {
396 return false;
397 }
398 }
399 }
400
401 return true;
402}
403
404//__________________________________________________________________________
405bool TrackExtrap::getAbsorberCorrectionParam(double trackXYZIn[3], double trackXYZOut[3], double pTotal,
406 double& pathLength, double& f0, double& f1, double& f2,
407 double& meanRho, double& totalELoss, double& sigmaELoss2)
408{
418
419 // Check whether the geometry is available
420 if (!gGeoManager) {
421 LOG(warning) << "geometry is missing";
422 return false;
423 }
424
425 // Initialize starting point and direction
426 pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0]) * (trackXYZOut[0] - trackXYZIn[0]) +
427 (trackXYZOut[1] - trackXYZIn[1]) * (trackXYZOut[1] - trackXYZIn[1]) +
428 (trackXYZOut[2] - trackXYZIn[2]) * (trackXYZOut[2] - trackXYZIn[2]));
429 if (pathLength < TGeoShape::Tolerance()) {
430 LOG(warning) << "path length is too small";
431 return false;
432 }
433 double b[3] = {(trackXYZOut[0] - trackXYZIn[0]) / pathLength, (trackXYZOut[1] - trackXYZIn[1]) / pathLength, (trackXYZOut[2] - trackXYZIn[2]) / pathLength};
434 TGeoNode* currentnode = gGeoManager->InitTrack(trackXYZIn, b);
435 if (!currentnode) {
436 LOG(warning) << "starting point out of geometry";
437 return false;
438 }
439
440 // loop over absorber slices and calculate absorber's parameters
441 f0 = f1 = f2 = meanRho = totalELoss = 0.;
442 double sigmaELoss(0.);
443 double zB = trackXYZIn[2];
444 double remainingPathLength = pathLength;
445 do {
446
447 // Get material properties
448 TGeoMaterial* material = currentnode->GetVolume()->GetMedium()->GetMaterial();
449 double rho = material->GetDensity(); // material density (g/cm3)
450 double x0 = material->GetRadLen(); // radiation-length (cm-1)
451 double atomicZ = material->GetZ(); // Z of material
452 double atomicZoverA(0.); // Z/A of material
453 if (material->IsMixture()) {
454 TGeoMixture* mixture = static_cast<TGeoMixture*>(material);
455 double sum(0.);
456 for (int iel = 0; iel < mixture->GetNelements(); ++iel) {
457 sum += mixture->GetWmixt()[iel];
458 atomicZoverA += mixture->GetWmixt()[iel] * mixture->GetZmixt()[iel] / mixture->GetAmixt()[iel];
459 }
460 atomicZoverA /= sum;
461 } else {
462 atomicZoverA = atomicZ / material->GetA();
463 }
464
465 // Get path length within this material
466 gGeoManager->FindNextBoundary(remainingPathLength);
467 double localPathLength = gGeoManager->GetStep() + 1.e-6;
468 // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step
469 if (localPathLength >= remainingPathLength) {
470 localPathLength = remainingPathLength;
471 } else {
472 currentnode = gGeoManager->Step();
473 if (!currentnode) {
474 LOG(warning) << "navigation failed";
475 return false;
476 }
477 if (!gGeoManager->IsEntering()) {
478 // make another small step to try to enter in new absorber slice
479 gGeoManager->SetStep(0.001);
480 currentnode = gGeoManager->Step();
481 if (!gGeoManager->IsEntering() || !currentnode) {
482 LOG(warning) << "navigation failed";
483 return false;
484 }
485 localPathLength += 0.001;
486 }
487 }
488
489 // calculate absorber's parameters
490 double zE = b[2] * localPathLength + zB;
491 double dzB = zB - trackXYZIn[2];
492 double dzE = zE - trackXYZIn[2];
493 f0 += localPathLength / x0;
494 f1 += (dzE * dzE - dzB * dzB) / b[2] / b[2] / x0 / 2.;
495 f2 += (dzE * dzE * dzE - dzB * dzB * dzB) / b[2] / b[2] / b[2] / x0 / 3.;
496 meanRho += localPathLength * rho;
497 totalELoss += betheBloch(pTotal, localPathLength, rho, atomicZ, atomicZoverA);
498 sigmaELoss += energyLossFluctuation(pTotal, localPathLength, rho, atomicZoverA);
499
500 // prepare next step
501 zB = zE;
502 remainingPathLength -= localPathLength;
503 } while (remainingPathLength > TGeoShape::Tolerance());
504
505 meanRho /= pathLength;
506 sigmaELoss2 = sigmaELoss * sigmaELoss;
507
508 return true;
509}
510
511//__________________________________________________________________________
512double TrackExtrap::getMCSAngle2(const TrackParam& param, double dZ, double x0)
513{
517 double bendingSlope = param.getBendingSlope();
518 double nonBendingSlope = param.getNonBendingSlope();
519 double inverseTotalMomentum2 = param.getInverseBendingMomentum() * param.getInverseBendingMomentum() *
520 (1.0 + bendingSlope * bendingSlope) /
521 (1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
522 // Path length in the material
523 double pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
524 // relativistic velocity
525 double velo = 1.;
526 // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
527 double theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength / x0));
528 return theta02 * theta02 * inverseTotalMomentum2 * pathLength / x0;
529}
530
531//__________________________________________________________________________
532void TrackExtrap::addMCSEffect(TrackParam& trackParam, double dZ, double x0)
533{
539
540 double bendingSlope = trackParam.getBendingSlope();
541 double nonBendingSlope = trackParam.getNonBendingSlope();
542 double inverseBendingMomentum = trackParam.getInverseBendingMomentum();
543 double inverseTotalMomentum2 = inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
544 (1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
545 // Path length in the material
546 double signedPathLength = dZ * TMath::Sqrt(1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
547 double pathLengthOverX0 = (x0 > 0.) ? TMath::Abs(signedPathLength) / x0 : TMath::Abs(signedPathLength);
548 // relativistic velocity
549 double velo = 1.;
550 // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
551 double theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLengthOverX0));
552 theta02 *= theta02 * inverseTotalMomentum2 * pathLengthOverX0;
553
554 double varCoor = (x0 > 0.) ? signedPathLength * signedPathLength * theta02 / 3. : 0.;
555 double varSlop = theta02;
556 double covCorrSlope = (x0 > 0.) ? signedPathLength * theta02 / 2. : 0.;
557
558 // Set MCS covariance matrix
559 TMatrixD newParamCov(trackParam.getCovariances());
560 // Non bending plane
561 newParamCov(0, 0) += varCoor;
562 newParamCov(0, 1) += covCorrSlope;
563 newParamCov(1, 0) += covCorrSlope;
564 newParamCov(1, 1) += varSlop;
565 // Bending plane
566 newParamCov(2, 2) += varCoor;
567 newParamCov(2, 3) += covCorrSlope;
568 newParamCov(3, 2) += covCorrSlope;
569 newParamCov(3, 3) += varSlop;
570
571 // Set momentum related covariances if B!=0
572 if (sFieldON) {
573 // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeY
574 double dqPxydSlopeX =
575 inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
576 double dqPxydSlopeY = -inverseBendingMomentum * nonBendingSlope * nonBendingSlope * bendingSlope /
577 (1. + bendingSlope * bendingSlope) /
578 (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
579 // Inverse bending momentum (due to dependences with bending and non bending slopes)
580 newParamCov(4, 0) += dqPxydSlopeX * covCorrSlope;
581 newParamCov(0, 4) += dqPxydSlopeX * covCorrSlope;
582 newParamCov(4, 1) += dqPxydSlopeX * varSlop;
583 newParamCov(1, 4) += dqPxydSlopeX * varSlop;
584 newParamCov(4, 2) += dqPxydSlopeY * covCorrSlope;
585 newParamCov(2, 4) += dqPxydSlopeY * covCorrSlope;
586 newParamCov(4, 3) += dqPxydSlopeY * varSlop;
587 newParamCov(3, 4) += dqPxydSlopeY * varSlop;
588 newParamCov(4, 4) += (dqPxydSlopeX * dqPxydSlopeX + dqPxydSlopeY * dqPxydSlopeY) * varSlop;
589 }
590
591 // Set new covariances
592 trackParam.setCovariances(newParamCov);
593}
594
595//__________________________________________________________________________
596void TrackExtrap::addMCSEffectInAbsorber(TrackParam& param, double signedPathLength, double f0, double f1, double f2)
597{
600
601 // absorber related covariance parameters
602 double bendingSlope = param.getBendingSlope();
603 double nonBendingSlope = param.getNonBendingSlope();
604 double inverseBendingMomentum = param.getInverseBendingMomentum();
605 double alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
606 (1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope); // velocity = 1
607 double pathLength = TMath::Abs(signedPathLength);
608 double varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2);
609 double covCorrSlope = TMath::Sign(1., signedPathLength) * alpha2 * (pathLength * f0 - f1);
610 double varSlop = alpha2 * f0;
611
612 // Set MCS covariance matrix
613 TMatrixD newParamCov(param.getCovariances());
614 // Non bending plane
615 newParamCov(0, 0) += varCoor;
616 newParamCov(0, 1) += covCorrSlope;
617 newParamCov(1, 0) += covCorrSlope;
618 newParamCov(1, 1) += varSlop;
619 // Bending plane
620 newParamCov(2, 2) += varCoor;
621 newParamCov(2, 3) += covCorrSlope;
622 newParamCov(3, 2) += covCorrSlope;
623 newParamCov(3, 3) += varSlop;
624
625 // Set momentum related covariances if B!=0
626 if (sFieldON) {
627 // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeY
628 double dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
629 double dqPxydSlopeY = -inverseBendingMomentum * nonBendingSlope * nonBendingSlope * bendingSlope /
630 (1. + bendingSlope * bendingSlope) / (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
631 // Inverse bending momentum (due to dependences with bending and non bending slopes)
632 newParamCov(4, 0) += dqPxydSlopeX * covCorrSlope;
633 newParamCov(0, 4) += dqPxydSlopeX * covCorrSlope;
634 newParamCov(4, 1) += dqPxydSlopeX * varSlop;
635 newParamCov(1, 4) += dqPxydSlopeX * varSlop;
636 newParamCov(4, 2) += dqPxydSlopeY * covCorrSlope;
637 newParamCov(2, 4) += dqPxydSlopeY * covCorrSlope;
638 newParamCov(4, 3) += dqPxydSlopeY * varSlop;
639 newParamCov(3, 4) += dqPxydSlopeY * varSlop;
640 newParamCov(4, 4) += (dqPxydSlopeX * dqPxydSlopeX + dqPxydSlopeY * dqPxydSlopeY) * varSlop;
641 }
642
643 // Set new covariances
644 param.setCovariances(newParamCov);
645}
646
647//__________________________________________________________________________
648double TrackExtrap::betheBloch(double pTotal, double pathLength, double rho, double atomicZ, double atomicZoverA)
649{
653
654 static constexpr double mK = 0.307075e-3; // [GeV*cm^2/g]
655 static constexpr double me = 0.511e-3; // [GeV/c^2]
656 static constexpr double x0 = 0.2 * 2.303; // density effect first junction point * 2.303
657 static constexpr double x1 = 3. * 2.303; // density effect second junction point * 2.303
658
659 double bg = pTotal / SMuMass; // beta*gamma
660 double bg2 = bg * bg;
661 double maxT = 2 * me * bg2; // neglecting the electron mass
662
663 // mean exitation energy (GeV)
664 double mI = (atomicZ < 13) ? (12. * atomicZ + 7.) * 1.e-9 : (9.76 * atomicZ + 58.8 * TMath::Power(atomicZ, -0.19)) * 1.e-9;
665
666 // density effect
667 double x = TMath::Log(bg);
668 double lhwI = TMath::Log(28.816 * 1e-9 * TMath::Sqrt(rho * atomicZoverA) / mI);
669 double d2(0.);
670 if (x > x1) {
671 d2 = lhwI + x - 0.5;
672 } else if (x > x0) {
673 double r = (x1 - x) / (x1 - x0);
674 d2 = lhwI + x - 0.5 + (0.5 - lhwI - x0) * r * r * r;
675 }
676
677 return pathLength * rho * (mK * atomicZoverA * (1 + bg2) / bg2 * (0.5 * TMath::Log(2 * me * bg2 * maxT / (mI * mI)) - bg2 / (1 + bg2) - d2));
678}
679
680//__________________________________________________________________________
681double TrackExtrap::energyLossFluctuation(double pTotal, double pathLength, double rho, double atomicZoverA)
682{
685
686 static constexpr double mK = 0.307075e-3; // GeV.g^-1.cm^2
687
688 double p2 = pTotal * pTotal;
689 double beta2 = p2 / (p2 + SMuMass * SMuMass);
690
691 double fwhm = 2. * mK * rho * pathLength * atomicZoverA / beta2; // FWHM of the energy loss Landau distribution
692
693 return fwhm / TMath::Sqrt(8. * log(2.)); // gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma)
694}
695
696//__________________________________________________________________________
697bool TrackExtrap::correctMCSEffectInAbsorber(TrackParam& param, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx,
698 double absZBeg, double pathLength, double f0, double f1, double f2)
699{
704
705 // Position of the Branson plane (spectro. (z<0))
706 double zB = (f1 > 0.) ? absZBeg - f2 / f1 : 0.;
707
708 // Add MCS effects to current parameter covariances (spectro. (z<0))
709 addMCSEffectInAbsorber(param, -pathLength, f0, f1, f2);
710
711 // Get track parameters and covariances in the Branson plane corrected for magnetic field effect
712 if (!extrapToZCov(param, zVtx)) {
713 return false;
714 }
716
717 // compute track parameters at vertex
718 TMatrixD newParam(5, 1);
719 newParam(0, 0) = xVtx;
720 newParam(1, 0) = (param.getNonBendingCoor() - xVtx) / (zB - zVtx);
721 newParam(2, 0) = yVtx;
722 newParam(3, 0) = (param.getBendingCoor() - yVtx) / (zB - zVtx);
723 newParam(4, 0) = param.getCharge() / param.p() *
724 TMath::Sqrt(1.0 + newParam(1, 0) * newParam(1, 0) + newParam(3, 0) * newParam(3, 0)) /
725 TMath::Sqrt(1.0 + newParam(3, 0) * newParam(3, 0));
726
727 // Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
728 TMatrixD paramCovP(param.getCovariances());
729 cov2CovP(param.getParameters(), paramCovP);
730
731 // Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system
732 TMatrixD paramCovVtx(5, 5);
733 paramCovVtx.Zero();
734 paramCovVtx(0, 0) = errXVtx * errXVtx;
735 paramCovVtx(1, 1) = paramCovP(0, 0);
736 paramCovVtx(2, 2) = errYVtx * errYVtx;
737 paramCovVtx(3, 3) = paramCovP(2, 2);
738 paramCovVtx(4, 4) = paramCovP(4, 4);
739 paramCovVtx(1, 3) = paramCovP(0, 2);
740 paramCovVtx(3, 1) = paramCovP(2, 0);
741 paramCovVtx(1, 4) = paramCovP(0, 4);
742 paramCovVtx(4, 1) = paramCovP(4, 0);
743 paramCovVtx(3, 4) = paramCovP(2, 4);
744 paramCovVtx(4, 3) = paramCovP(4, 2);
745
746 // Jacobian of the transformation (XVtx, X, YVtx, Y, q*PTot) -> (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx)
747 TMatrixD jacob(5, 5);
748 jacob.UnitMatrix();
749 jacob(1, 0) = -1. / (zB - zVtx);
750 jacob(1, 1) = 1. / (zB - zVtx);
751 jacob(3, 2) = -1. / (zB - zVtx);
752 jacob(3, 3) = 1. / (zB - zVtx);
753
754 // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx) coordinate system
755 TMatrixD tmp(paramCovVtx, TMatrixD::kMultTranspose, jacob);
756 TMatrixD newParamCov(jacob, TMatrixD::kMult, tmp);
757
758 // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q/PyzVtx) coordinate system
759 covP2Cov(newParam, newParamCov);
760
761 // Set parameters and covariances at vertex
762 param.setParameters(newParam);
763 param.setZ(zVtx);
764 param.setCovariances(newParamCov);
765
766 return true;
767}
768
769//__________________________________________________________________________
770void TrackExtrap::correctELossEffectInAbsorber(TrackParam& param, double eLoss, double sigmaELoss2)
771{
773
774 // Get parameter covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
775 TMatrixD newParamCov(param.getCovariances());
776 cov2CovP(param.getParameters(), newParamCov);
777
778 // Compute new parameters corrected for energy loss
779 double p = param.p();
780 double e = TMath::Sqrt(p * p + SMuMass * SMuMass);
781 double eCorr = e + eLoss;
782 double pCorr = TMath::Sqrt(eCorr * eCorr - SMuMass * SMuMass);
783 double nonBendingSlope = param.getNonBendingSlope();
784 double bendingSlope = param.getBendingSlope();
785 param.setInverseBendingMomentum(param.getCharge() / pCorr *
786 TMath::Sqrt(1.0 + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope) /
787 TMath::Sqrt(1.0 + bendingSlope * bendingSlope));
788
789 // Add effects of energy loss fluctuation to covariances
790 newParamCov(4, 4) += eCorr * eCorr / pCorr / pCorr * sigmaELoss2;
791
792 // Get new parameter covariances in (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
793 covP2Cov(param.getParameters(), newParamCov);
794
795 // Set new parameter covariances
796 param.setCovariances(newParamCov);
797}
798
799//__________________________________________________________________________
800void TrackExtrap::cov2CovP(const TMatrixD& param, TMatrixD& cov)
801{
804
805 // charge * total momentum
806 double qPTot = TMath::Sqrt(1. + param(1, 0) * param(1, 0) + param(3, 0) * param(3, 0)) /
807 TMath::Sqrt(1. + param(3, 0) * param(3, 0)) / param(4, 0);
808
809 // Jacobian of the opposite transformation
810 TMatrixD jacob(5, 5);
811 jacob.UnitMatrix();
812 jacob(4, 1) = qPTot * param(1, 0) / (1. + param(1, 0) * param(1, 0) + param(3, 0) * param(3, 0));
813 jacob(4, 3) = -qPTot * param(1, 0) * param(1, 0) * param(3, 0) /
814 (1. + param(3, 0) * param(3, 0)) / (1. + param(1, 0) * param(1, 0) + param(3, 0) * param(3, 0));
815 jacob(4, 4) = -qPTot / param(4, 0);
816
817 // compute covariances in new coordinate system
818 TMatrixD tmp(cov, TMatrixD::kMultTranspose, jacob);
819 cov.Mult(jacob, tmp);
820}
821
822//__________________________________________________________________________
823void TrackExtrap::covP2Cov(const TMatrixD& param, TMatrixD& covP)
824{
827
828 // charge * total momentum
829 double qPTot = TMath::Sqrt(1. + param(1, 0) * param(1, 0) + param(3, 0) * param(3, 0)) /
830 TMath::Sqrt(1. + param(3, 0) * param(3, 0)) / param(4, 0);
831
832 // Jacobian of the transformation
833 TMatrixD jacob(5, 5);
834 jacob.UnitMatrix();
835 jacob(4, 1) = param(4, 0) * param(1, 0) / (1. + param(1, 0) * param(1, 0) + param(3, 0) * param(3, 0));
836 jacob(4, 3) = -param(4, 0) * param(1, 0) * param(1, 0) * param(3, 0) /
837 (1. + param(3, 0) * param(3, 0)) / (1. + param(1, 0) * param(1, 0) + param(3, 0) * param(3, 0));
838 jacob(4, 4) = -param(4, 0) / qPTot;
839
840 // compute covariances in new coordinate system
841 TMatrixD tmp(covP, TMatrixD::kMultTranspose, jacob);
842 covP.Mult(jacob, tmp);
843}
844
845//__________________________________________________________________________
846void TrackExtrap::convertTrackParamForExtrap(TrackParam& trackParam, double forwardBackward, double* v3)
847{
851 v3[0] = trackParam.getNonBendingCoor(); // X
852 v3[1] = trackParam.getBendingCoor(); // Y
853 v3[2] = trackParam.getZ(); // Z
854 double pYZ = TMath::Abs(1.0 / trackParam.getInverseBendingMomentum());
855 double pZ = pYZ / TMath::Sqrt(1.0 + trackParam.getBendingSlope() * trackParam.getBendingSlope());
856 v3[6] = TMath::Sqrt(pYZ * pYZ + pZ * pZ * trackParam.getNonBendingSlope() * trackParam.getNonBendingSlope()); // P
857 v3[5] = -forwardBackward * pZ / v3[6]; // PZ/P spectro. z<0
858 v3[3] = trackParam.getNonBendingSlope() * v3[5]; // PX/P
859 v3[4] = trackParam.getBendingSlope() * v3[5]; // PY/P
860}
861
862//__________________________________________________________________________
863void TrackExtrap::recoverTrackParam(double* v3, double charge, TrackParam& trackParam)
864{
868 trackParam.setNonBendingCoor(v3[0]); // X
869 trackParam.setBendingCoor(v3[1]); // Y
870 trackParam.setZ(v3[2]); // Z
871 double pYZ = v3[6] * TMath::Sqrt((1. - v3[3]) * (1. + v3[3]));
872 trackParam.setInverseBendingMomentum(charge / pYZ);
873 trackParam.setBendingSlope(v3[4] / v3[5]);
874 trackParam.setNonBendingSlope(v3[3] / v3[5]);
875}
876
877//__________________________________________________________________________
878bool TrackExtrap::extrapToZRungekutta(TrackParam& trackParam, double zEnd)
879{
883 if (trackParam.getZ() == zEnd) {
884 return true; // nothing to be done if same Z
885 }
886 double forwardBackward; // +1 if forward, -1 if backward
887 if (zEnd < trackParam.getZ()) {
888 forwardBackward = 1.0; // spectro. z<0
889 } else {
890 forwardBackward = -1.0;
891 }
892 // sign of charge (sign of fInverseBendingMomentum if forward motion)
893 // must be changed if backward extrapolation
894 double chargeExtrap = forwardBackward * TMath::Sign(double(1.0), trackParam.getInverseBendingMomentum());
895 double v3[7] = {0.}, v3New[7] = {0.};
896 double dZ(0.), step(0.);
897 int stepNumber = 0;
898
899 // Extrapolation loop (until within tolerance or the track turn around)
900 double residue = zEnd - trackParam.getZ();
901 while (TMath::Abs(residue) > SRungeKuttaMaxResidue) {
902
903 dZ = zEnd - trackParam.getZ();
904 // step length assuming linear trajectory
905 step = dZ * TMath::Sqrt(1.0 + trackParam.getBendingSlope() * trackParam.getBendingSlope() +
906 trackParam.getNonBendingSlope() * trackParam.getNonBendingSlope());
907 convertTrackParamForExtrap(trackParam, forwardBackward, v3);
908
909 do { // reduce step length while zEnd overstepped
910 if (++stepNumber > SMaxStepNumber) {
911 LOG(warning) << "Too many trials";
912 return false;
913 }
914 step = TMath::Abs(step);
915 if (!extrapOneStepRungekutta(chargeExtrap, step, v3, v3New)) {
916 return false;
917 }
918 residue = zEnd - v3New[2];
919 step *= dZ / (v3New[2] - trackParam.getZ());
920 } while (residue * dZ < 0 && TMath::Abs(residue) > SRungeKuttaMaxResidue);
921
922 if (v3New[5] * v3[5] < 0) { // the track turned around
923 LOG(warning) << "The track turned around";
924 return false;
925 }
926
927 recoverTrackParam(v3New, chargeExtrap * forwardBackward, trackParam);
928 }
929
930 // terminate the extropolation with a straight line up to the exact "zEnd" value
931 trackParam.setNonBendingCoor(trackParam.getNonBendingCoor() + residue * trackParam.getNonBendingSlope());
932 trackParam.setBendingCoor(trackParam.getBendingCoor() + residue * trackParam.getBendingSlope());
933 trackParam.setZ(zEnd);
934
935 return true;
936}
937
938//__________________________________________________________________________
939bool TrackExtrap::extrapToZRungekuttaV2(TrackParam& trackParam, double zEnd)
940{
944
945 if (trackParam.getZ() == zEnd) {
946 return true; // nothing to be done if same Z
947 }
948
949 double residue = zEnd - trackParam.getZ();
950 double forwardBackward = (residue < 0) ? 1. : -1.; // +1 if forward, -1 if backward
951 double v3[7] = {0.};
952 convertTrackParamForExtrap(trackParam, forwardBackward, v3);
953 double charge = TMath::Sign(double(1.), trackParam.getInverseBendingMomentum());
954
955 // Extrapolation loop (until within tolerance or the track turn around)
956 double v3New[7] = {0.};
957 int stepNumber = 0;
958 while (true) {
959
960 if (++stepNumber > SMaxStepNumber) {
961 LOG(warning) << "Too many trials";
962 return false;
963 }
964
965 // step length assuming linear trajectory
966 double slopeX = v3[3] / v3[5];
967 double slopeY = v3[4] / v3[5];
968 double step = TMath::Abs(residue) * TMath::Sqrt(1.0 + slopeX * slopeX + slopeY * slopeY);
969
970 if (!extrapOneStepRungekutta(forwardBackward * charge, step, v3, v3New)) {
971 return false;
972 }
973
974 if (v3New[5] * v3[5] < 0) {
975 LOG(warning) << "The track turned around";
976 return false;
977 }
978
979 residue = zEnd - v3New[2];
980 if (TMath::Abs(residue) < SRungeKuttaMaxResidueV2) {
981 break;
982 }
983
984 for (int i = 0; i < 7; ++i) {
985 v3[i] = v3New[i];
986 }
987
988 // invert the sens of propagation if the track went too far
989 if (forwardBackward * residue > 0) {
990 forwardBackward = -forwardBackward;
991 v3[3] = -v3[3];
992 v3[4] = -v3[4];
993 v3[5] = -v3[5];
994 }
995 }
996
997 recoverTrackParam(v3New, charge, trackParam);
998
999 // terminate the extrapolation with a straight line up to the exact "zEnd" value
1000 trackParam.setNonBendingCoor(trackParam.getNonBendingCoor() + residue * trackParam.getNonBendingSlope());
1001 trackParam.setBendingCoor(trackParam.getBendingCoor() + residue * trackParam.getBendingSlope());
1002 trackParam.setZ(zEnd);
1003
1004 return true;
1005}
1006
1007//__________________________________________________________________________
1008bool TrackExtrap::extrapOneStepRungekutta(double charge, double step, const double* vect, double* vout)
1009{
1032
1033 double h2(0.), h4(0.), f[4] = {0.};
1034 double xyzt[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
1035 double a(0.), b(0.), c(0.), ph(0.), ph2(0.);
1036 double secxs[4] = {0.}, secys[4] = {0.}, seczs[4] = {0.}, hxp[3] = {0.};
1037 double g1(0.), g2(0.), g3(0.), g4(0.), g5(0.), g6(0.), ang2(0.), dxt(0.), dyt(0.), dzt(0.);
1038 double est(0.), at(0.), bt(0.), ct(0.), cba(0.);
1039 double f1(0.), f2(0.), f3(0.), f4(0.), rho(0.), tet(0.), hnorm(0.), hp(0.), rho1(0.), sint(0.), cost(0.);
1040
1041 double x(0.);
1042 double y(0.);
1043 double z(0.);
1044
1045 double xt(0.);
1046 double yt(0.);
1047 double zt(0.);
1048
1049 double maxit = 1992;
1050 double maxcut = 11;
1051
1052 const double kdlt = 1e-4;
1053 const double kdlt32 = kdlt / 32.;
1054 const double kthird = 1. / 3.;
1055 const double khalf = 0.5;
1056 const double kec = 2.9979251e-4;
1057
1058 const double kpisqua = 9.86960440109;
1059 const int kix = 0;
1060 const int kiy = 1;
1061 const int kiz = 2;
1062 const int kipx = 3;
1063 const int kipy = 4;
1064 const int kipz = 5;
1065
1066 // *.
1067 // *. ------------------------------------------------------------------
1068 // *.
1069 // * this constant is for units cm,gev/c and kgauss
1070 // *
1071 int iter = 0;
1072 int ncut = 0;
1073 for (int j = 0; j < 7; j++) {
1074 vout[j] = vect[j];
1075 }
1076
1077 double pinv = kec * charge / vect[6];
1078 double tl = 0.;
1079 double h = step;
1080 double rest(0.);
1081
1082 do {
1083 rest = step - tl;
1084 if (TMath::Abs(h) > TMath::Abs(rest)) {
1085 h = rest;
1086 }
1087 // cmodif: call gufld(vout,f) changed into:
1088 TGeoGlobalMagField::Instance()->Field(vout, f);
1089 ++sNCallField;
1090
1091 // *
1092 // * start of integration
1093 // *
1094 x = vout[0];
1095 y = vout[1];
1096 z = vout[2];
1097 a = vout[3];
1098 b = vout[4];
1099 c = vout[5];
1100
1101 h2 = khalf * h;
1102 h4 = khalf * h2;
1103 ph = pinv * h;
1104 ph2 = khalf * ph;
1105 secxs[0] = (b * f[2] - c * f[1]) * ph2;
1106 secys[0] = (c * f[0] - a * f[2]) * ph2;
1107 seczs[0] = (a * f[1] - b * f[0]) * ph2;
1108 ang2 = (secxs[0] * secxs[0] + secys[0] * secys[0] + seczs[0] * seczs[0]);
1109 if (ang2 > kpisqua) {
1110 break;
1111 }
1112
1113 dxt = h2 * a + h4 * secxs[0];
1114 dyt = h2 * b + h4 * secys[0];
1115 dzt = h2 * c + h4 * seczs[0];
1116 xt = x + dxt;
1117 yt = y + dyt;
1118 zt = z + dzt;
1119 // *
1120 // * second intermediate point
1121 // *
1122
1123 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1124 if (est > h) {
1125 if (ncut++ > maxcut) {
1126 break;
1127 }
1128 h *= khalf;
1129 continue;
1130 }
1131
1132 xyzt[0] = xt;
1133 xyzt[1] = yt;
1134 xyzt[2] = zt;
1135
1136 // cmodif: call gufld(xyzt,f) changed into:
1137 TGeoGlobalMagField::Instance()->Field(xyzt, f);
1138 ++sNCallField;
1139
1140 at = a + secxs[0];
1141 bt = b + secys[0];
1142 ct = c + seczs[0];
1143
1144 secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
1145 secys[1] = (ct * f[0] - at * f[2]) * ph2;
1146 seczs[1] = (at * f[1] - bt * f[0]) * ph2;
1147 at = a + secxs[1];
1148 bt = b + secys[1];
1149 ct = c + seczs[1];
1150 secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
1151 secys[2] = (ct * f[0] - at * f[2]) * ph2;
1152 seczs[2] = (at * f[1] - bt * f[0]) * ph2;
1153 dxt = h * (a + secxs[2]);
1154 dyt = h * (b + secys[2]);
1155 dzt = h * (c + seczs[2]);
1156 xt = x + dxt;
1157 yt = y + dyt;
1158 zt = z + dzt;
1159 at = a + 2. * secxs[2];
1160 bt = b + 2. * secys[2];
1161 ct = c + 2. * seczs[2];
1162
1163 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1164 if (est > 2. * TMath::Abs(h)) {
1165 if (ncut++ > maxcut) {
1166 break;
1167 }
1168 h *= khalf;
1169 continue;
1170 }
1171
1172 xyzt[0] = xt;
1173 xyzt[1] = yt;
1174 xyzt[2] = zt;
1175
1176 // cmodif: call gufld(xyzt,f) changed into:
1177 TGeoGlobalMagField::Instance()->Field(xyzt, f);
1178 ++sNCallField;
1179
1180 z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
1181 y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
1182 x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
1183
1184 secxs[3] = (bt * f[2] - ct * f[1]) * ph2;
1185 secys[3] = (ct * f[0] - at * f[2]) * ph2;
1186 seczs[3] = (at * f[1] - bt * f[0]) * ph2;
1187 a = a + (secxs[0] + secxs[3] + 2. * (secxs[1] + secxs[2])) * kthird;
1188 b = b + (secys[0] + secys[3] + 2. * (secys[1] + secys[2])) * kthird;
1189 c = c + (seczs[0] + seczs[3] + 2. * (seczs[1] + seczs[2])) * kthird;
1190
1191 est = TMath::Abs(secxs[0] + secxs[3] - (secxs[1] + secxs[2])) +
1192 TMath::Abs(secys[0] + secys[3] - (secys[1] + secys[2])) +
1193 TMath::Abs(seczs[0] + seczs[3] - (seczs[1] + seczs[2]));
1194
1195 if (est > kdlt && TMath::Abs(h) > 1.e-4) {
1196 if (ncut++ > maxcut) {
1197 break;
1198 }
1199 h *= khalf;
1200 continue;
1201 }
1202
1203 ncut = 0;
1204 // * if too many iterations, go to helix
1205 if (iter++ > maxit) {
1206 break;
1207 }
1208
1209 tl += h;
1210 if (est < kdlt32) {
1211 h *= 2.;
1212 }
1213 cba = 1. / TMath::Sqrt(a * a + b * b + c * c);
1214 vout[0] = x;
1215 vout[1] = y;
1216 vout[2] = z;
1217 vout[3] = cba * a;
1218 vout[4] = cba * b;
1219 vout[5] = cba * c;
1220 rest = step - tl;
1221 if (step < 0.) {
1222 rest = -rest;
1223 }
1224 if (rest < 1.e-5 * TMath::Abs(step)) {
1225 return true;
1226 }
1227
1228 } while (1);
1229
1230 // angle too big, use helix
1231 LOG(warning) << "Ruge-Kutta failed: switch to helix";
1232
1233 f1 = f[0];
1234 f2 = f[1];
1235 f3 = f[2];
1236 f4 = TMath::Sqrt(f1 * f1 + f2 * f2 + f3 * f3);
1237 if (f4 < 1.e-10) {
1238 LOG(warning) << "magnetic field at (" << xyzt[0] << ", " << xyzt[1] << ", " << xyzt[2] << ") = " << f4 << ": giving up";
1239 return false;
1240 }
1241 rho = -f4 * pinv;
1242 tet = rho * step;
1243
1244 hnorm = 1. / f4;
1245 f1 = f1 * hnorm;
1246 f2 = f2 * hnorm;
1247 f3 = f3 * hnorm;
1248
1249 hxp[0] = f2 * vect[kipz] - f3 * vect[kipy];
1250 hxp[1] = f3 * vect[kipx] - f1 * vect[kipz];
1251 hxp[2] = f1 * vect[kipy] - f2 * vect[kipx];
1252
1253 hp = f1 * vect[kipx] + f2 * vect[kipy] + f3 * vect[kipz];
1254
1255 rho1 = 1. / rho;
1256 sint = TMath::Sin(tet);
1257 cost = 2. * TMath::Sin(khalf * tet) * TMath::Sin(khalf * tet);
1258
1259 g1 = sint * rho1;
1260 g2 = cost * rho1;
1261 g3 = (tet - sint) * hp * rho1;
1262 g4 = -cost;
1263 g5 = sint;
1264 g6 = cost * hp;
1265
1266 vout[kix] = vect[kix] + g1 * vect[kipx] + g2 * hxp[0] + g3 * f1;
1267 vout[kiy] = vect[kiy] + g1 * vect[kipy] + g2 * hxp[1] + g3 * f2;
1268 vout[kiz] = vect[kiz] + g1 * vect[kipz] + g2 * hxp[2] + g3 * f3;
1269
1270 vout[kipx] = vect[kipx] + g4 * vect[kipx] + g5 * hxp[0] + g6 * f1;
1271 vout[kipy] = vect[kipy] + g4 * vect[kipy] + g5 * hxp[1] + g6 * f2;
1272 vout[kipz] = vect[kipz] + g4 * vect[kipz] + g5 * hxp[2] + g6 * f3;
1273
1274 return true;
1275}
1276
1277//__________________________________________________________________________
1279{
1281 LOG(info) << "number of times extrapToZCov() is called = " << sNCallExtrapToZCov;
1282 LOG(info) << "number of times Field() is called = " << sNCallField;
1283}
1284
1285} // namespace mch
1286} // namespace o2
int16_t charge
Definition RawEventData.h:5
int32_t i
constexpr int p2()
uint32_t j
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
Definition of tools for track extrapolation.
Definition of the MCH track parameters for internal use.
Class for time synchronization of RawReader instances.
static bool extrapToZ(TrackParam &trackParam, double zEnd)
static void addMCSEffect(TrackParam &trackParam, double dZ, double x0)
static void setField()
static bool extrapToVertex(TrackParam &trackParam, double xVtx, double yVtx, double zVtx, double errXVtx, double errYVtx)
Definition TrackExtrap.h:62
static void linearExtrapToZ(TrackParam &trackParam, double zEnd)
static double getMCSAngle2(const TrackParam &param, double dZ, double x0)
static bool extrapToZCov(TrackParam &trackParam, double zEnd, bool updatePropagator=false)
static double getBendingMomentumFromImpactParam(double impactParam)
static bool extrapToMID(TrackParam &trackParam)
static double getImpactParamFromBendingMomentum(double bendingMomentum)
static void printNCalls()
static void linearExtrapToZCov(TrackParam &trackParam, double zEnd, bool updatePropagator=false)
track parameters for internal use
Definition TrackParam.h:34
Double_t getInverseBendingMomentum() const
return inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
Definition TrackParam.h:67
void setZ(Double_t z)
set Z coordinate (cm)
Definition TrackParam.h:49
Bool_t hasCovariances() const
return kTRUE if the covariance matrix exist, kFALSE if not
Definition TrackParam.h:95
Double_t getNonBendingCoor() const
return non bending coordinate (cm)
Definition TrackParam.h:51
Double_t getZ() const
return Z coordinate (cm)
Definition TrackParam.h:47
void addParameters(const TMatrixD &parameters)
add track parameters
Definition TrackParam.h:87
const TMatrixD & getParameters() const
return track parameters
Definition TrackParam.h:81
Double_t getNonBendingSlope() const
return non bending slope (cm ** -1)
Definition TrackParam.h:55
void setCovariances(const TMatrixD &covariances)
const TMatrixD & getCovariances() const
void setBendingCoor(Double_t bendingCoor)
set bending coordinate (cm)
Definition TrackParam.h:61
Double_t getBendingCoor() const
return bending coordinate (cm)
Definition TrackParam.h:59
void setParameters(const TMatrixD &parameters)
set track parameters
Definition TrackParam.h:83
void updatePropagator(const TMatrixD &propagator)
Double_t p() const
Double_t getBendingSlope() const
return bending slope (cm ** -1)
Definition TrackParam.h:63
void setNonBendingCoor(Double_t nonBendingCoor)
set non bending coordinate (cm)
Definition TrackParam.h:53
float sum(float s, o2::dcs::DataPointValue v)
Definition dcs-ccdb.cxx:39
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint GLfloat GLfloat GLfloat x1
Definition glcorearb.h:5034
GLdouble f
Definition glcorearb.h:310
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLint y
Definition glcorearb.h:270
GLuint GLfloat x0
Definition glcorearb.h:5034
GLboolean r
Definition glcorearb.h:1233
GLfloat GLfloat GLfloat GLfloat v3
Definition glcorearb.h:814
GLenum GLfloat param
Definition glcorearb.h:271
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
value_T f3
Definition TrackUtils.h:93
value_T sint
Definition TrackUtils.h:77
kp1 *kp2 *value_T bg2
Definition TrackUtils.h:131
value_T f4
Definition TrackUtils.h:94
value_T step
Definition TrackUtils.h:42
value_T f1
Definition TrackUtils.h:91
value_T d2
Definition TrackUtils.h:135
value_T gpu::gpustd::array< value_T, 7 > & vect
Definition TrackUtils.h:42
kp1 *kp2 *value_T beta2
Definition TrackUtils.h:131
constexpr value_T mK
Definition TrackUtils.h:127
value_T rho
Definition TrackUtils.h:36
value_T f2
Definition TrackUtils.h:92
value_T tet
Definition TrackUtils.h:75
const value_T lhwI
Definition TrackUtils.h:137
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"