ATLAS Offline Software
STEP_Propagator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // Implementation file for class Trk::STEP_Propagator
8 // (c) ATLAS Detector software
10 // Based on RungeKuttaPropagator Version 1.0 21/05/2004 I.Gavrilenko
11 // Version 0.1 16/2/2005 Esben Lund
12 // Version 1.0 27/7/2006 Esben Lund
14 
16 //
20 
31 #include "TrkSurfaces/Surface.h"
34 // CLHEP
35 #include "CLHEP/Random/RandFlat.h"
36 #include "CLHEP/Random/RandGauss.h"
37 #include "GaudiKernel/PhysicalConstants.h"
38 // Gaudi
43 //
44 #include <cmath>
45 //
46 #include "CxxUtils/inline_hints.h"
47 
48 namespace{
49 
50 using Cache = Trk::STEP_Propagator::Cache;
52 // methods for magnetic field information
54 void
55 getField(Cache& cache, const double* ATH_RESTRICT R, double* ATH_RESTRICT H)
56 {
57  if (cache.m_solenoid) {
58  cache.m_fieldCache.getFieldZR(R, H);
59  } else {
60  cache.m_fieldCache.getField(R, H);
61  }
62 }
63 
64 void
65 getFieldGradient(Cache& cache, const double* ATH_RESTRICT R, double* ATH_RESTRICT H, double* ATH_RESTRICT dH)
66 {
67  if (cache.m_solenoid) {
68  cache.m_fieldCache.getFieldZR(R, H, dH);
69  } else {
70  cache.m_fieldCache.getField(R, H, dH);
71  }
72 }
73 
75 // Get the magnetic field and gradients
76 // Input: Globalposition
77 // Output: BG, which contains Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz, dBy/dx,
78 // dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz
80 
81 void
82 getMagneticField(Cache& cache,
83  const Amg::Vector3D& position,
84  bool getGradients,
85  double* ATH_RESTRICT BG)
86 {
87  constexpr double magScale = 10000.;
88  const double* R = position.data();
89  double H[3];
90  double dH[9];
91 
92  if (getGradients &&
93  cache.m_includeBgradients) { // field gradients needed and available
94 
95  getFieldGradient(cache, R, H, dH);
96  BG[0] = H[0] * magScale;
97  BG[1] = H[1] * magScale;
98  BG[2] = H[2] * magScale;
99  BG[3] = dH[0] * magScale;
100  BG[4] = dH[1] * magScale;
101  BG[5] = dH[2] * magScale;
102  BG[6] = dH[3] * magScale;
103  BG[7] = dH[4] * magScale;
104  BG[8] = dH[5] * magScale;
105  BG[9] = dH[6] * magScale;
106  BG[10] = dH[7] * magScale;
107  BG[11] = dH[8] * magScale;
108 
109  } else { // Homogenous field or no gradients needed, only retrieve the field
110  // strength.
111  getField(cache, R, H);
112  BG[0] = H[0] * magScale;
113  BG[1] = H[1] * magScale;
114  BG[2] = H[2] * magScale;
115 
116  for (int i = 3; i < 12; ++i) { // Set gradients to zero
117  BG[i] = 0.;
118  }
119  }
120 }
121 
123 // Create straight line in case of q/p = 0
125 
126 std::unique_ptr<Trk::TrackParameters>
127 createStraightLine(const Trk::TrackParameters* inputTrackParameters)
128 {
129  AmgVector(5) lp = inputTrackParameters->parameters();
130  lp[Trk::qOverP] = 1. / 1e10;
131 
132  if (inputTrackParameters->type() == Trk::Curvilinear) {
133  return std::make_unique<Trk::CurvilinearParameters>(
134  inputTrackParameters->position(),
135  lp[2],
136  lp[3],
137  lp[4],
138  (inputTrackParameters->covariance() ? std::optional<AmgSymMatrix(5)>(*inputTrackParameters->covariance())
139  : std::nullopt));
140  }
141  return inputTrackParameters->associatedSurface().createUniqueTrackParameters(
142  lp[0],
143  lp[1],
144  lp[2],
145  lp[3],
146  lp[4],
147  (inputTrackParameters->covariance() ? std::optional<AmgSymMatrix(5)>(*inputTrackParameters->covariance())
148  : std::nullopt));
149 }
150 
151 std::unique_ptr<Trk::TrackParameters>
152 propagateNeutral(const Trk::TrackParameters& parm,
153  std::vector<Trk::DestSurf>& targetSurfaces,
154  Trk::PropDirection propDir,
155  std::vector<unsigned int>& solutions,
156  double& path,
157  bool usePathLimit,
158  bool returnCurv)
159 {
160 
161  // find nearest valid intersection
162  double tol = 0.001;
163  solutions.clear();
164  std::vector<std::pair<unsigned int, double>> currentDist;
165  currentDist.clear();
166  std::vector<std::pair<const Trk::Surface*, Trk::BoundaryCheck>>::iterator sIter = targetSurfaces.begin();
167  std::vector<std::pair<unsigned int, double>>::iterator oIter = currentDist.begin();
168  std::vector<Trk::DestSurf>::iterator sBeg = targetSurfaces.begin();
169 
170  const Amg::Vector3D& position(parm.position());
171  Amg::Vector3D direction(parm.momentum().normalized());
172 
173  for (; sIter != targetSurfaces.end(); ++sIter) {
174  Trk::DistanceSolution distSol = (*sIter).first->straightLineDistanceEstimate(position, direction);
175  if (distSol.numberOfSolutions() > 0) {
176  unsigned int numSf = sIter - sBeg;
177  if (distSol.first() > tol) {
178  if (!currentDist.empty()) {
179  oIter = currentDist.end();
180  while (oIter != currentDist.begin() && distSol.first() < (*(oIter - 1)).second)
181  --oIter;
182  oIter = currentDist.insert(oIter, std::pair<unsigned int, double>(numSf, distSol.first()));
183  } else {
184  currentDist.emplace_back(numSf, distSol.first());
185  }
186  }
187  if (distSol.numberOfSolutions() > 1 && distSol.second() > tol) {
188  if (!currentDist.empty()) {
189  oIter = currentDist.end();
190  while (oIter != currentDist.begin() && distSol.second() < (*(oIter - 1)).second)
191  --oIter;
192  oIter = currentDist.insert(oIter, std::pair<unsigned int, double>(numSf, distSol.second()));
193  } else {
194  currentDist.emplace_back(numSf, distSol.second());
195  }
196  }
197  }
198  }
199 
200  // loop over surfaces, find valid intersections
201  for (oIter = currentDist.begin(); oIter != currentDist.end(); ++oIter) {
202  Amg::Vector3D xsct = position + propDir * direction * ((*oIter).second);
203  if (targetSurfaces[(*oIter).first].first->isOnSurface(xsct, (*oIter).second, 0.001, 0.001)) {
204  if (usePathLimit && path > 0. && path < ((*oIter).second)) {
205  Amg::Vector3D endpoint = position + propDir * direction * path;
206  return std::make_unique<Trk::CurvilinearParameters>(endpoint, parm.momentum(), parm.charge());
207  }
208  path = (*oIter).second;
209  solutions.push_back((*oIter).first);
210  const Trk::Surface* sf = targetSurfaces[(*oIter).first].first;
211 
212  if (returnCurv || sf->type() == Trk::SurfaceType::Cone) {
213  return std::make_unique<Trk::CurvilinearParameters>(xsct, parm.momentum(), parm.charge());
214  }
215  return sf->createUniqueTrackParameters(xsct, parm.momentum(), parm.charge(), std::nullopt);
216  }
217  }
218 
219  return nullptr;
220 }
221 
222 void
223 clearMaterialEffects(Trk::STEP_Propagator::Cache& cache)
224 {
225  cache.m_delIoni = 0.;
226  cache.m_delRad = 0.;
227  cache.m_sigmaIoni = 0.;
228  cache.m_sigmaRad = 0.;
229 }
230 
232 // Calculate energy loss in MeV/mm.
234 double
235 dEds(Cache& cache, double p)
236 {
237  cache.m_delIoni = 0.;
238  cache.m_delRad = 0.;
239  cache.m_kazL = 0.;
240 
241  if (cache.m_particle == Trk::geantino || cache.m_particle == Trk::nonInteractingMuon)
242  return 0.;
243  if (cache.m_material->x0() == 0 || cache.m_material->averageZ() == 0)
244  return 0.;
245 
246  cache.m_delIoni =
247  Trk::MaterialInteraction::dEdl_ionization(p, *(cache.m_material), cache.m_particle, cache.m_sigmaIoni, cache.m_kazL);
248 
249  cache.m_delRad = Trk::MaterialInteraction::dEdl_radiation(p, *(cache.m_material), cache.m_particle, cache.m_sigmaRad);
250 
251  double eLoss = cache.m_MPV ? 0.9 * cache.m_delIoni + 0.15 * cache.m_delRad : cache.m_delIoni + cache.m_delRad;
252 
253  return eLoss;
254 }
255 
257 // dg/dlambda for non-electrons (g=dEdX and lambda=q/p).
259 
260 double
261 dgdlambda(Cache& cache, double l)
262 {
263  if (cache.m_particle == Trk::geantino || cache.m_particle == Trk::nonInteractingMuon)
264  return 0.;
265  if (cache.m_material->x0() == 0 || cache.m_material->averageZ() == 0)
266  return 0.;
267  if (cache.m_material->zOverAtimesRho() == 0)
268  return 0.;
269 
270  double p = std::abs(1. / l);
271  double m = Trk::ParticleMasses::mass[cache.m_particle];
273  double E = std::sqrt(p * p + m * m);
274  double beta = p / E;
275  double gamma = E / m;
276  double I = 16.e-6 * std::pow(cache.m_material->averageZ(), 0.9);
277  double kaz = 0.5 * 30.7075 * cache.m_material->zOverAtimesRho();
278 
279  // Bethe-Bloch
280  double lnCore = 4. * me * me / (m * m * m * m * I * I * l * l * l * l) / (1. + 2. * gamma * me / m + me * me / m * m);
281  double lnCore_deriv =
282  -4. * me * me / (m * m * m * m * I * I) *
283  std::pow(l * l * l * l + 2. * gamma * l * l * l * l * me / m + l * l * l * l * me * me / (m * m), -2) *
284  (4. * l * l * l + 8. * me * l * l * l * gamma / m - 2. * me * l / (m * m * m * gamma) +
285  4. * l * l * l * me * me / (m * m));
286  double ln_deriv = 2. * l * m * m * std::log(lnCore) + lnCore_deriv / (lnCore * beta * beta);
287  double Bethe_Bloch_deriv = -kaz * ln_deriv;
288 
289  // density effect, only valid for high energies (gamma > 10 -> p > 1GeV for muons)
290  if (gamma > 10.) {
291  double delta = 2. * std::log(28.816e-6 * std::sqrt(1000. * cache.m_material->zOverAtimesRho()) / I) +
292  2. * std::log(beta * gamma) - 1.;
293  double delta_deriv = -2. / (l * beta * beta) + 2. * delta * l * m * m;
294  Bethe_Bloch_deriv += kaz * delta_deriv;
295  }
296 
297  // Bethe-Heitler
298  double Bethe_Heitler_deriv = me * me / (m * m * cache.m_material->x0() * l * l * l * E);
299 
300  // Radiative corrections (e+e- pair production + photonuclear) for muons at energies above 8 GeV and below 1 TeV
301  double radiative_deriv = 0.;
302  if ((cache.m_particle == Trk::muon) && (E > 8000.)) {
303  if (E < 1.e6) {
304  radiative_deriv = 6.803e-5 / (cache.m_material->x0() * l * l * l * E) +
305  2. * 2.278e-11 / (cache.m_material->x0() * l * l * l) -
306  3. * 9.899e-18 * E / (cache.m_material->x0() * l * l * l);
307  } else {
308  radiative_deriv = 9.253e-5 / (cache.m_material->x0() * l * l * l * E);
309  }
310  }
311 
312  // return the total derivative
313  if (cache.m_MPV) {
314  return 0.9 * Bethe_Bloch_deriv + 0.15 * Bethe_Heitler_deriv + 0.15 * radiative_deriv; // Most probable value
315  } else {
316  return Bethe_Bloch_deriv + Bethe_Heitler_deriv + radiative_deriv; // Mean value
317  }
318 }
319 
321 // update material effects // to follow change of material
323 void
324 updateMaterialEffects(Cache& cache, double mom, double sinTheta, double path)
325 {
326  // collects material effects in between material dumps ( m_covariance )
327 
328  // collect straggling
329  if (cache.m_straggling) {
330  cache.m_covariance(4, 4) += cache.m_stragglingVariance;
331  cache.m_stragglingVariance = 0.;
332  }
333 
334  if (!cache.m_multipleScattering)
335  return;
336 
337  double totalMomentumLoss = mom - cache.m_matupd_lastmom;
338  double pathSinceLastUpdate = path - cache.m_matupd_lastpath;
339 
340  double pathAbs = std::abs(pathSinceLastUpdate);
341 
342  if (pathAbs < 1.e-03)
343  return;
344 
345  double matX0 = cache.m_material->x0();
346 
347  // calculate number of layers : TODO : thickness to be adjusted
348  int msLayers = int(pathAbs / cache.m_layXmax / matX0) + 1;
349 
350  double massSquared = Trk::ParticleMasses::mass[cache.m_particle] * Trk::ParticleMasses::mass[cache.m_particle];
351  double layerThickness = pathAbs / msLayers;
352  double radiationLengths = layerThickness / matX0;
353  sinTheta = std::max(sinTheta, 1e-20); // avoid division by zero
354  double remainingPath = pathAbs;
355  double momentum = 0.;
356  double mom2 = 0.;
357  double E = 0.;
358  double beta = 0.;
359  double dLambdads = 0.;
360  double thetaVariance = 0.;
361 
362  double average_dEds = totalMomentumLoss / pathAbs;
363 
364  double cumulatedVariance = cache.m_inputThetaVariance + cache.m_combinedCovariance(3, 3) + cache.m_covariance(3, 3);
365  if (cumulatedVariance < 0) {
366  cumulatedVariance = 0;
367  }
368  double cumulatedX0 = 0.;
369 
370  bool useCache = cache.m_extrapolationCache != nullptr;
371  if (useCache) {
372  double dX0 = std::abs(cache.m_combinedThickness) - pathAbs / matX0;
373  if (dX0 < 0)
374  dX0 = 0.;
375  if (cache.m_extrapolationCache->x0tot() > 0)
376  cumulatedX0 = cache.m_extrapolationCache->x0tot() + dX0;
377  }
378 
379  // calculate multiple scattering by summing the contributions from the layers
380  for (int layer = 1; layer <= msLayers; ++layer) {
381 
382  // calculate momentum in the middle of the layer by assuming a linear momentum loss
383  momentum = cache.m_matupd_lastmom + totalMomentumLoss * (layer - 0.5) / msLayers;
384 
385  mom2 = momentum * momentum;
386  E = std::sqrt(mom2 + massSquared);
387  beta = momentum / E;
388 
389  double C0 = 13.6 * 13.6 / mom2 / beta / beta;
390  double C1 = 2 * 0.038 * C0;
391  double C2 = 0.038 * 0.038 * C0;
392 
393  double MS2 = radiationLengths;
394 
395  dLambdads = -cache.m_charge * average_dEds * E / (momentum * momentum * momentum);
396  remainingPath -= layerThickness;
397 
398  // simple - step dependent formula
399  // thetaVariance = C0*MS2 + C1*MS2*log(MS2) + C2*MS2*log(MS2)*log(MS2);
400 
401  // replaced by the step-size independent code below
402  if (!useCache) {
403  cumulatedX0 = cumulatedVariance / C0;
404  if (cumulatedX0 > 0.001) {
405  double lX0 = log(cumulatedX0);
406  cumulatedX0 = cumulatedX0 / (1 + 2 * 0.038 * lX0 + 0.038 * 0.038 * lX0 * lX0);
407  }
408  }
409 
410  double MS2s = MS2 + cumulatedX0;
411  thetaVariance = C0 * MS2 + C1 * MS2s * log(MS2s) + C2 * MS2s * log(MS2s) * log(MS2s);
412 
413  if (cumulatedX0 > 0.001) {
414  double lX0 = log(cumulatedX0);
415  thetaVariance += -C1 * cumulatedX0 * lX0 - C2 * cumulatedX0 * lX0 * lX0;
416  }
417 
418  // Calculate ms covariance contributions and scale
419  double varScale = cache.m_scatteringScale * cache.m_scatteringScale;
420  double positionVariance = thetaVariance * (layerThickness * layerThickness / 3. + remainingPath * layerThickness +
421  remainingPath * remainingPath);
422  cache.m_covariance(0, 0) += varScale * positionVariance;
423  cache.m_covariance(1, 1) += varScale * positionVariance;
424  cache.m_covariance.fillSymmetric(
425  2, 0, cache.m_covariance(2, 0) + varScale * thetaVariance / (sinTheta) * (layerThickness / 2. + remainingPath));
426  cache.m_covariance(2, 2) += varScale * thetaVariance / (sinTheta * sinTheta);
427  cache.m_covariance.fillSymmetric(
428  3, 1, cache.m_covariance(3, 1) + varScale * thetaVariance * (-layerThickness / 2. - remainingPath));
429  cache.m_covariance(3, 3) += varScale * thetaVariance;
430  cache.m_covariance(4, 4) +=
431  varScale * 3. * thetaVariance * thetaVariance * layerThickness * layerThickness * dLambdads * dLambdads;
432  cumulatedVariance += varScale * thetaVariance;
433  cumulatedX0 += MS2;
434  }
435 
436  cache.m_matupd_lastmom = mom;
437  cache.m_matupd_lastpath = path;
438 }
439 
441 // Multiple scattering and straggling contribution to the covariance matrix
443 
444 void
445 covarianceContribution(Cache& cache,
446  const Trk::TrackParameters* trackParameters,
447  double path,
448  const Trk::TrackParameters* targetParms,
449  AmgSymMatrix(5)& measurementCovariance)
450 {
451  // kinematics
452  double finalMomentum = targetParms->momentum().mag();
453 
454  // first update to make sure all material counted
455  updateMaterialEffects(cache, finalMomentum, sin(trackParameters->momentum().theta()), path);
456 
457  double Jac[21];
459 
460  // Transform multiple scattering and straggling covariance from curvilinear to local system
461  AmgSymMatrix(5) localMSCov =
462  Trk::RungeKuttaUtils::newCovarianceMatrix(Jac, cache.m_combinedCovariance + cache.m_covariance);
463 
464  measurementCovariance += localMSCov;
465 }
466 
468 // Multiple scattering and straggling contribution to the covariance matrix in curvilinear representation
470 
471 void
472 covarianceContribution(Cache& cache,
473  const Trk::TrackParameters* trackParameters,
474  double path,
475  double finalMomentum,
476  AmgSymMatrix(5)& measurementCovariance)
477 {
478  // first update to make sure all material counted
479  updateMaterialEffects(cache, finalMomentum, sin(trackParameters->momentum().theta()), path);
480 
481  // Add measurement errors and multiple scattering + straggling covariance
482  measurementCovariance += cache.m_combinedCovariance + cache.m_covariance;
483 }
484 
486 // Runge Kutta trajectory model (units->mm,MeV,kGauss)
487 //
488 // Runge Kutta method is adaptive Runge-Kutta-Nystrom.
489 // This is the default STEP method
491 
492 // We compile this package with optimization, even in debug builds; otherwise,
493 // the heavy use of Eigen makes it too slow. However, from here we may call
494 // to out-of-line Eigen code that is linked from other DSOs; in that case,
495 // it would not be optimized. Avoid this by forcing all Eigen code
496 // to be inlined here if possible.
498 bool
499 rungeKuttaStep(Cache& cache,
500  bool errorPropagation,
501  double& h,
502  double* ATH_RESTRICT P,
503  double* ATH_RESTRICT dDir,
504  double* ATH_RESTRICT BG1,
505  bool& firstStep,
506  double& distanceStepped)
507 {
508  constexpr double sol = 0.0299792458; //speed of light * kilogauss in method units
509  double charge = 0;
510  P[6] >= 0. ? charge = 1. : charge = -1.; // Set charge
511  double lambda1 = P[6];
512  double lambda2 = P[6]; // Store inverse momentum for Jacobian transport
513  double lambda3 = P[6];
514  double lambda4 = P[6];
515  double dP1 = 0.;
516  double dP2 = 0.;
517  double dP3 = 0.;
518  double dP4 = 0.; // dp/ds = -g*E/p for positive g=dE/ds
519  double dL1 = 0.;
520  double dL2 = 0.;
521  double dL3 = 0.;
522  double dL4 = 0.; // factor used for calculating dCM/dCM, P[41], in the Jacobian.
523  double initialMomentum = std::abs(1. / P[6]); // Set initial momentum
524  const Amg::Vector3D initialPos(P[0], P[1], P[2]); // Set initial values for position
525  const Amg::Vector3D initialDir(P[3], P[4], P[5]); // Set initial values for direction.
526  if (!Amg::saneVector(initialPos) || !Amg::saneVector(initialDir)) return false;
527  // Directions at the different points. Used by the error propagation
528  Amg::Vector3D dir1;
529  Amg::Vector3D dir2;
530  Amg::Vector3D dir3;
531  Amg::Vector3D dir4;
532  // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz, dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy,dBz/dz
533  double BG23[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
534  // The gradients are used by the error propagation
535  double BG4[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. };
536  double g = 0.; // Energyloss in Mev/mm.
537  double dgdl = 0.; // dg/dlambda in Mev/mm.
538  double particleMass = Trk::ParticleMasses::mass[cache.m_particle]; // Get particle mass from ParticleHypothesis
539  int steps = 0;
540 
541  // POINT 1. Only calculate this once per step, even if step is rejected. This
542  // point is independant of the step length, h
543  double momentum = initialMomentum; // Current momentum
544  if (cache.m_energyLoss && cache.m_matPropOK) {
545  g = dEds(cache, momentum); // Use the same energy loss throughout the step.
546  double E = std::sqrt(momentum * momentum + particleMass * particleMass);
547  dP1 = g * E / momentum;
548  if (errorPropagation) {
549  if (cache.m_includeGgradient) {
550  dgdl = dgdlambda(cache, lambda1); // Use this value throughout the step.
551  }
552  dL1 =
553  -lambda1 * lambda1 * g * E * (3. - (momentum * momentum) / (E * E)) - lambda1 * lambda1 * lambda1 * E * dgdl;
554  }
555  }
556 
557  Amg::Vector3D pos = initialPos; // Coordinate solution
558  Amg::Vector3D dir = initialDir; // Direction solution
559  dir1 = dir;
560  if (firstStep) { // Poll BG1 if this is the first step, else use recycled BG4
561  firstStep = false;
562  // Get the gradients needed for the error propagation if errorPropagation=true
563  getMagneticField(cache,
564  pos,
565  errorPropagation,
566  BG1); // Get the gradients needed for the error propagation if errorPropagation=true
567  }
568  // Lorentz force, d2r/ds2 = lambda * (dir x B)
569  Amg::Vector3D k1(
570  dir.y() * BG1[2] - dir.z() * BG1[1], dir.z() * BG1[0] - dir.x() * BG1[2], dir.x() * BG1[1] - dir.y() * BG1[0]);
571  k1 = sol * lambda1 * k1;
572 
573  while (true) { // Repeat step until error estimate is within the requested tolerance
574  if (steps++ > cache.m_maxSteps)
575  return false; // Abort propagation
576  // POINT 2
577  if (cache.m_energyLoss && cache.m_matPropOK) {
578  momentum = initialMomentum + (h / 2.) * dP1;
579  if (momentum <= cache.m_momentumCutOff)
580  return false; // Abort propagation
581  double E = std::sqrt(momentum * momentum + particleMass * particleMass);
582  dP2 = g * E / momentum;
583  lambda2 = charge / momentum;
584  if (errorPropagation) {
585  dL2 =
586  -lambda2 * lambda2 * g * E * (3. - (momentum * momentum) / (E * E)) - lambda2 * lambda2 * lambda2 * E * dgdl;
587  }
588  }
589  pos = initialPos + (h / 2.) * initialDir + (h * h / 8.) * k1;
590  dir = initialDir + (h / 2.) * k1;
591  dir2 = dir;
592 
593  getMagneticField(cache, pos, errorPropagation, BG23);
594  Amg::Vector3D k2(dir.y() * BG23[2] - dir.z() * BG23[1],
595  dir.z() * BG23[0] - dir.x() * BG23[2],
596  dir.x() * BG23[1] - dir.y() * BG23[0]);
597  k2 = sol * lambda2 * k2;
598 
599  // POINT 3. Same position and B-field as point 2.
600  if (cache.m_energyLoss && cache.m_matPropOK) {
601  momentum = initialMomentum + (h / 2.) * dP2;
602  if (momentum <= cache.m_momentumCutOff)
603  return false; // Abort propagation
604  double E = std::sqrt(momentum * momentum + particleMass * particleMass);
605  dP3 = g * E / momentum;
606  lambda3 = charge / momentum;
607  if (errorPropagation) {
608  dL3 =
609  -lambda3 * lambda3 * g * E * (3. - (momentum * momentum) / (E * E)) - lambda3 * lambda3 * lambda3 * E * dgdl;
610  }
611  }
612  dir = initialDir + (h / 2.) * k2;
613  dir3 = dir;
614  Amg::Vector3D k3(dir.y() * BG23[2] - dir.z() * BG23[1],
615  dir.z() * BG23[0] - dir.x() * BG23[2],
616  dir.x() * BG23[1] - dir.y() * BG23[0]);
617  k3 = sol * lambda3 * k3;
618 
619  // POINT 4
620  if (cache.m_energyLoss && cache.m_matPropOK) {
621  momentum = initialMomentum + h * dP3;
622  if (momentum <= cache.m_momentumCutOff)
623  return false; // Abort propagation
624  double E = std::sqrt(momentum * momentum + particleMass * particleMass);
625  dP4 = g * E / momentum;
626  lambda4 = charge / momentum;
627  if (errorPropagation) {
628  dL4 =
629  -lambda4 * lambda4 * g * E * (3. - (momentum * momentum) / (E * E)) - lambda4 * lambda4 * lambda4 * E * dgdl;
630  }
631  }
632  pos = initialPos + h * initialDir + (h * h / 2.) * k3;
633  dir = initialDir + h * k3;
634  dir4 = dir;
635 
636  // MT Field cache is stored in cache
637  getMagneticField(cache, pos, errorPropagation, BG4);
638  Amg::Vector3D k4(
639  dir.y() * BG4[2] - dir.z() * BG4[1], dir.z() * BG4[0] - dir.x() * BG4[2], dir.x() * BG4[1] - dir.y() * BG4[0]);
640  k4 = sol * lambda4 * k4;
641 
642  // Estimate local error and avoid division by zero
643  Amg::Vector3D errorPos((h * h) * (k1 - k2 - k3 + k4));
644  double errorEstimate = std::max(errorPos.mag(), 1e-20);
645 
646  // Use the error estimate to calculate new step length. h is returned by reference.
647  distanceStepped = h; // Store old step length.
648  h = h * std::min(std::max(0.25, std::sqrt(std::sqrt(cache.m_tolerance / errorEstimate))), 4.);
649 
650  // Repeat step with the new step size if error is too big.
651  if (errorEstimate > 4. * cache.m_tolerance)
652  continue;
653 
654  // Step was ok. Store solutions.
655  // Update positions.
656  pos = initialPos + distanceStepped * initialDir + (distanceStepped * distanceStepped / 6.) * (k1 + k2 + k3);
657  P[0] = pos.x();
658  P[1] = pos.y();
659  P[2] = pos.z();
660 
661  // update directions
662  dir = initialDir + (distanceStepped / 6.) * (k1 + 2. * k2 + 2. * k3 + k4);
663  P[3] = dir.x();
664  P[4] = dir.y();
665  P[5] = dir.z();
666 
667  // Normalize direction
668  double norm = 1. / std::sqrt(P[3] * P[3] + P[4] * P[4] + P[5] * P[5]);
669  P[3] = norm * P[3];
670  P[4] = norm * P[4];
671  P[5] = norm * P[5];
672 
673  // Update inverse momentum if energyloss is switched on
674  if (cache.m_energyLoss && cache.m_matPropOK) {
675  momentum = initialMomentum + (distanceStepped / 6.) * (dP1 + 2. * dP2 + 2. * dP3 + dP4);
676  if (momentum <= cache.m_momentumCutOff)
677  return false; // Abort propagation
678  P[6] = charge / momentum;
679  }
680 
681  // dDir provides a small correction to the final tiny step in PropagateWithJacobian
682  dDir[0] = k4.x();
683  dDir[1] = k4.y();
684  dDir[2] = k4.z();
685 
686  // Transport Jacobian using the same step length, points and magnetic fields as for the track parameters
687  // BG contains Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz, dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz
688  // 0 1 2 3 4 5 6 7 8 9 10 11
689  if (errorPropagation) {
690  double initialjLambda = P[41]; // dLambda/dLambda
691  double jLambda = initialjLambda;
692  double jdL1 = 0.;
693  double jdL2 = 0.;
694  double jdL3 = 0.;
695  double jdL4 = 0.;
696 
697  for (int i = 7; i < 42; i += 7) {
698 
699  // POINT 1
700  const Amg::Vector3D initialjPos(P[i], P[i + 1], P[i + 2]);
701  const Amg::Vector3D initialjDir(P[i + 3], P[i + 4], P[i + 5]);
702 
703  if (!Amg::saneVector(initialjPos) || !Amg::saneVector(initialjDir)) return false;
704 
705  Amg::Vector3D jPos = initialjPos;
706  Amg::Vector3D jDir = initialjDir;
707 
708  // B-field terms
709  Amg::Vector3D jk1(jDir.y() * BG1[2] - jDir.z() * BG1[1],
710  jDir.z() * BG1[0] - jDir.x() * BG1[2],
711  jDir.x() * BG1[1] - jDir.y() * BG1[0]);
712 
713  // B-field gradient terms
714  if (cache.m_includeBgradients) {
716  (dir1.y() * BG1[9] - dir1.z() * BG1[6]) * jPos.x() + (dir1.y() * BG1[10] - dir1.z() * BG1[7]) * jPos.y() +
717  (dir1.y() * BG1[11] - dir1.z() * BG1[8]) * jPos.z(),
718  (dir1.z() * BG1[3] - dir1.x() * BG1[9]) * jPos.x() + (dir1.z() * BG1[4] - dir1.x() * BG1[10]) * jPos.y() +
719  (dir1.z() * BG1[5] - dir1.x() * BG1[11]) * jPos.z(),
720  (dir1.x() * BG1[6] - dir1.y() * BG1[3]) * jPos.x() + (dir1.x() * BG1[7] - dir1.y() * BG1[4]) * jPos.y() +
721  (dir1.x() * BG1[8] - dir1.y() * BG1[5]) * jPos.z());
722  jk1 = jk1 + C;
723  }
724  jk1 = sol * lambda1 * jk1;
725 
726  // Last column of the A-matrix
727  if (i == 35) { // Only dLambda/dLambda, in the last row of the jacobian, is different from zero
728  jdL1 = dL1 * jLambda; // Energy loss term. dL1 = 0 if no material effects -> jLambda = P[41] is unchanged
729  jk1 = jk1 + k1 * jLambda / lambda1; // B-field terms
730  }
731 
732  // POINT 2
733  jPos = initialjPos + (distanceStepped / 2.) * initialjDir + (distanceStepped * distanceStepped / 8.) * jk1;
734  jDir = initialjDir + (distanceStepped / 2.) * jk1;
735 
736  Amg::Vector3D jk2(jDir.y() * BG23[2] - jDir.z() * BG23[1],
737  jDir.z() * BG23[0] - jDir.x() * BG23[2],
738  jDir.x() * BG23[1] - jDir.y() * BG23[0]);
739 
740  if (cache.m_includeBgradients) {
741  Amg::Vector3D C((dir2.y() * BG23[9] - dir2.z() * BG23[6]) * jPos.x() +
742  (dir2.y() * BG23[10] - dir2.z() * BG23[7]) * jPos.y() +
743  (dir2.y() * BG23[11] - dir2.z() * BG23[8]) * jPos.z(),
744  (dir2.z() * BG23[3] - dir2.x() * BG23[9]) * jPos.x() +
745  (dir2.z() * BG23[4] - dir2.x() * BG23[10]) * jPos.y() +
746  (dir2.z() * BG23[5] - dir2.x() * BG23[11]) * jPos.z(),
747  (dir2.x() * BG23[6] - dir2.y() * BG23[3]) * jPos.x() +
748  (dir2.x() * BG23[7] - dir2.y() * BG23[4]) * jPos.y() +
749  (dir2.x() * BG23[8] - dir2.y() * BG23[5]) * jPos.z());
750  jk2 = jk2 + C;
751  }
752  jk2 = sol * lambda2 * jk2;
753 
754  if (i == 35) {
755  jLambda = initialjLambda + (distanceStepped / 2.) * jdL1;
756  jdL2 = dL2 * jLambda;
757  jk2 = jk2 + k2 * jLambda / lambda2;
758  }
759 
760  // POINT 3
761  jDir = initialjDir + (distanceStepped / 2.) * jk2;
762 
763  Amg::Vector3D jk3(jDir.y() * BG23[2] - jDir.z() * BG23[1],
764  jDir.z() * BG23[0] - jDir.x() * BG23[2],
765  jDir.x() * BG23[1] - jDir.y() * BG23[0]);
766 
767  if (cache.m_includeBgradients) {
768  Amg::Vector3D C((dir3.y() * BG23[9] - dir3.z() * BG23[6]) * jPos.x() +
769  (dir3.y() * BG23[10] - dir3.z() * BG23[7]) * jPos.y() +
770  (dir3.y() * BG23[11] - dir3.z() * BG23[8]) * jPos.z(),
771  (dir3.z() * BG23[3] - dir3.x() * BG23[9]) * jPos.x() +
772  (dir3.z() * BG23[4] - dir3.x() * BG23[10]) * jPos.y() +
773  (dir3.z() * BG23[5] - dir3.x() * BG23[11]) * jPos.z(),
774  (dir3.x() * BG23[6] - dir3.y() * BG23[3]) * jPos.x() +
775  (dir3.x() * BG23[7] - dir3.y() * BG23[4]) * jPos.y() +
776  (dir3.x() * BG23[8] - dir3.y() * BG23[5]) * jPos.z());
777  jk3 = jk3 + C;
778  }
779  jk3 = sol * lambda3 * jk3;
780 
781  if (i == 35) {
782  jLambda = initialjLambda + (distanceStepped / 2.) * jdL2;
783  jdL3 = dL3 * jLambda;
784  jk3 = jk3 + k3 * jLambda / lambda3;
785  }
786 
787  // POINT 4
788  jPos = initialjPos + distanceStepped * initialjDir + (distanceStepped * distanceStepped / 2.) * jk3;
789  jDir = initialjDir + distanceStepped * jk3;
790 
791  Amg::Vector3D jk4(jDir.y() * BG4[2] - jDir.z() * BG4[1],
792  jDir.z() * BG4[0] - jDir.x() * BG4[2],
793  jDir.x() * BG4[1] - jDir.y() * BG4[0]);
794 
795  if (cache.m_includeBgradients) {
797  (dir4.y() * BG4[9] - dir4.z() * BG4[6]) * jPos.x() + (dir4.y() * BG4[10] - dir4.z() * BG4[7]) * jPos.y() +
798  (dir4.y() * BG4[11] - dir4.z() * BG4[8]) * jPos.z(),
799  (dir4.z() * BG4[3] - dir4.x() * BG4[9]) * jPos.x() + (dir4.z() * BG4[4] - dir4.x() * BG4[10]) * jPos.y() +
800  (dir4.z() * BG4[5] - dir4.x() * BG4[11]) * jPos.z(),
801  (dir4.x() * BG4[6] - dir4.y() * BG4[3]) * jPos.x() + (dir4.x() * BG4[7] - dir4.y() * BG4[4]) * jPos.y() +
802  (dir4.x() * BG4[8] - dir4.y() * BG4[5]) * jPos.z());
803  jk4 = jk4 + C;
804  }
805  jk4 = sol * lambda4 * jk4;
806 
807  if (i == 35) {
808  jLambda = initialjLambda + distanceStepped * jdL3;
809  jdL4 = dL4 * jLambda;
810  jk4 = jk4 + k4 * jLambda / lambda4;
811  }
812 
813  // solution
814  jPos =
815  initialjPos + distanceStepped * initialjDir + (distanceStepped * distanceStepped / 6.) * (jk1 + jk2 + jk3);
816  jDir = initialjDir + (distanceStepped / 6.) * (jk1 + 2. * jk2 + 2. * jk3 + jk4);
817  if (i == 35) {
818  jLambda = initialjLambda + (distanceStepped / 6.) * (jdL1 + 2. * jdL2 + 2. * jdL3 + jdL4);
819  }
820 
821  // Update positions
822  P[i] = jPos.x();
823  P[i + 1] = jPos.y();
824  P[i + 2] = jPos.z();
825 
826  // update directions
827  P[i + 3] = jDir.x();
828  P[i + 4] = jDir.y();
829  P[i + 5] = jDir.z();
830  }
831  P[41] = jLambda; // update dLambda/dLambda
832  }
833 
834  // Store BG4 for use as BG1 in the next step
835  for (int i = 0; i < 12; ++i) {
836  BG1[i] = BG4[i];
837  }
838  return true;
839  }
840 }
842 // Runge Kutta main program for propagation with or without Jacobian
844 bool
845 propagateWithJacobianImpl(Cache& cache,
846  bool errorPropagation,
847  Trk::SurfaceType surfaceType,
848  double* ATH_RESTRICT targetSurface,
849  double* ATH_RESTRICT P,
850  double& path)
851 {
852  double maxPath = cache.m_maxPath; // Max path allowed
853  double* pos = &P[0]; // Start coordinates
854  double* dir = &P[3]; // Start directions
855  double dDir[3] = { 0., 0., 0. }; // Start directions derivs. Zero in case of no RK steps
856  int targetPassed = 0; // how many times have we passed the target?
857  double previousDistance = 0.;
858  double distanceStepped = 0.;
859  bool distanceEstimationSuccessful = false; // Was the distance estimation successful?
860  double BG1[12] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. }; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
861  // dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz at first point
862  bool firstStep = true; // Poll BG1, else recycle BG4
863  double absolutePath = 0.; // Absolute path to register oscillating behaviour
864  int steps = 0;
865  path = 0.; // signed path of the trajectory
866  double mom = 0.;
867  double beta = 1.;
868 
869  double distanceToTarget = Trk::RungeKuttaUtils::stepEstimator(
870  surfaceType, targetSurface, P, distanceEstimationSuccessful);
871  if (!distanceEstimationSuccessful) {
872  return false;
873  }
874 
875  // Set initial step length to 100mm or the distance to the target surface.
876  double h = 0;
877  distanceToTarget > 100. ? h = 100. : distanceToTarget < -100. ? h = -100. : h = distanceToTarget;
878 
879  // Step to within distanceTolerance, then do the rest as a Taylor expansion.
880  // Keep distanceTolerance within [1 nanometer, 10 microns].
881  // This means that no final Taylor expansions beyond 10 microns and no
882  // Runge-Kutta steps less than 1 nanometer are allowed.
883  double distanceTolerance = std::min(std::max(std::abs(distanceToTarget) * cache.m_tolerance, 1e-6), 1e-2);
884 
885  while (std::abs(distanceToTarget) > distanceTolerance) { // Step until within tolerance
886  // Do the step. Stop the propagation if the energy goes below m_momentumCutOff
887  if (!rungeKuttaStep(cache, errorPropagation, h, P, dDir, BG1, firstStep, distanceStepped))
888  return false;
889  path += distanceStepped;
890  absolutePath += std::abs(distanceStepped);
891 
892  if (std::abs(distanceStepped) > 0.001) {
893  cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL * log(std::abs(distanceStepped)); // the non-linear term
894  }
895  // update straggling covariance
896  if (errorPropagation && cache.m_straggling) {
897  double sigTot2 = cache.m_sigmaIoni * cache.m_sigmaIoni + cache.m_sigmaRad * cache.m_sigmaRad;
898  // /(beta*beta*p*p*p*p) transforms Var(E) to Var(q/p)
899  mom = std::abs(1. / P[6]);
900  beta = mom / std::sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass);
901  double bp2 = beta * mom * mom;
902  cache.m_stragglingVariance += sigTot2 / (bp2 * bp2) * distanceStepped * distanceStepped;
903  }
904  if (cache.m_matstates || errorPropagation)
905  cache.m_combinedEloss.update(cache.m_delIoni * distanceStepped,
906  cache.m_sigmaIoni * std::abs(distanceStepped),
907  cache.m_delRad * distanceStepped,
908  cache.m_sigmaRad * std::abs(distanceStepped),
909  cache.m_MPV);
910  // Calculate new distance to target
911  previousDistance = std::abs(distanceToTarget);
912  distanceToTarget = Trk::RungeKuttaUtils::stepEstimator(
913  surfaceType, targetSurface, P, distanceEstimationSuccessful);
914  if (!distanceEstimationSuccessful)
915  return false;
916 
917  // If h and distance are in opposite directions, target is passed. Flip propagation direction
918  if (h * distanceToTarget < 0.) {
919  h = -h; // Flip direction
920  ++targetPassed;
921  }
922  // don't step beyond surface
923  if (std::abs(h) > std::abs(distanceToTarget))
924  h = distanceToTarget;
925 
926  // Abort if maxPath is reached or solution is diverging
927  if ((targetPassed > 3 && std::abs(distanceToTarget) >= previousDistance) || (absolutePath > maxPath))
928  return false;
929 
930  if (steps++ > cache.m_maxSteps)
931  return false; // Too many steps, something is wrong
932  }
933 
934  if (cache.m_material && cache.m_material->x0() != 0.)
935  cache.m_combinedThickness += std::abs(path) / cache.m_material->x0();
936 
937  // Use Taylor expansions to step the remaining distance (typically microns).
938  path = path + distanceToTarget;
939 
940  // pos = pos + h*dir + 1/2*h*h*dDir. Second order Taylor expansion.
941  pos[0] = pos[0] + distanceToTarget * (dir[0] + 0.5 * distanceToTarget * dDir[0]);
942  pos[1] = pos[1] + distanceToTarget * (dir[1] + 0.5 * distanceToTarget * dDir[1]);
943  pos[2] = pos[2] + distanceToTarget * (dir[2] + 0.5 * distanceToTarget * dDir[2]);
944 
945  // dir = dir + h*dDir. First order Taylor expansion (Euler).
946  dir[0] = dir[0] + distanceToTarget * dDir[0];
947  dir[1] = dir[1] + distanceToTarget * dDir[1];
948  dir[2] = dir[2] + distanceToTarget * dDir[2];
949 
950  // Normalize dir
951  double norm = 1. / std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
952  dir[0] = norm * dir[0];
953  dir[1] = norm * dir[1];
954  dir[2] = norm * dir[2];
955  P[42] = dDir[0];
956  P[43] = dDir[1];
957  P[44] = dDir[2];
958  return true;
959 }
960 
962 // Main function for track propagation with or without jacobian
964 std::unique_ptr<Trk::TrackParameters>
965 propagateRungeKuttaImpl(Cache& cache,
966  bool errorPropagation,
967  const Trk::TrackParameters& inputTrackParameters,
968  const Trk::Surface& targetSurface,
969  Trk::PropDirection propagationDirection,
970  const Trk::MagneticFieldProperties& mft,
972  const Trk::BoundaryCheck& boundaryCheck,
973  double* Jacobian,
974  bool returnCurv)
975 {
976  // Store for later use
977  cache.m_particle = particle;
978  cache.m_charge = inputTrackParameters.charge();
979 
980  std::unique_ptr<Trk::TrackParameters> trackParameters{};
981 
982  // Bfield mode
983  mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
984 
985  // Check inputvalues
986  if (cache.m_tolerance <= 0.)
987  return nullptr;
988  if (cache.m_momentumCutOff < 0.)
989  return nullptr;
990 
991  // Set momentum to 1e10 (straight line) and charge to + if q/p is zero
992  if (inputTrackParameters.parameters()[Trk::qOverP] == 0) {
993  trackParameters = createStraightLine(&inputTrackParameters);
994  if (!trackParameters) {
995  return nullptr;
996  }
997  } else {
998  trackParameters.reset(inputTrackParameters.clone());
999  }
1000 
1001  if (std::abs(1. / trackParameters->parameters()[Trk::qOverP]) <= cache.m_momentumCutOff) {
1002  return nullptr;
1003  }
1004 
1005  // Check for empty volumes. If x != x then x is not a number.
1006  if (cache.m_matPropOK && ((cache.m_material->zOverAtimesRho() == 0.) || (cache.m_material->x0() == 0.) ||
1007  (cache.m_material->zOverAtimesRho() != cache.m_material->zOverAtimesRho()))) {
1008  cache.m_matPropOK = false;
1009  }
1010  if (cache.m_matPropOK && cache.m_straggling)
1011  cache.m_stragglingVariance = 0.;
1012  if (errorPropagation || cache.m_matstates) {
1013  cache.m_combinedCovariance.setZero();
1014  cache.m_covariance.setZero();
1015  // this needs debugging
1016  cache.m_inputThetaVariance = trackParameters->covariance() ? (*trackParameters->covariance())(3, 3) : 0.;
1017  cache.m_combinedEloss.set(0., 0., 0., 0., 0., 0.);
1018  cache.m_combinedThickness = 0.;
1019  }
1020 
1021  if (!Trk::RungeKuttaUtils::transformLocalToGlobal(errorPropagation, *trackParameters, cache.m_P)) {
1022  return nullptr;
1023  }
1024  double path = 0.;
1025 
1026  const Amg::Transform3D& T = targetSurface.transform();
1027  Trk::SurfaceType ty = targetSurface.type();
1028 
1029  if (ty == Trk::SurfaceType::Plane || ty == Trk::SurfaceType::Disc) {
1030  double s[4];
1031  double d = T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
1032 
1033  if (d >= 0.) {
1034  s[0] = T(0, 2);
1035  s[1] = T(1, 2);
1036  s[2] = T(2, 2);
1037  s[3] = d;
1038  } else {
1039  s[0] = -T(0, 2);
1040  s[1] = -T(1, 2);
1041  s[2] = -T(2, 2);
1042  s[3] = -d;
1043  }
1044  if (!propagateWithJacobianImpl(cache, errorPropagation, ty, s, cache.m_P, path)) {
1045  return nullptr;
1046  }
1047  }
1048 
1049  else if (ty == Trk::SurfaceType::Line) {
1050 
1051  double s[6] = { T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2) };
1052  if (!propagateWithJacobianImpl(cache, errorPropagation, ty, s, cache.m_P, path)) {
1053  return nullptr;
1054  }
1055  }
1056 
1057  else if (ty == Trk::SurfaceType::Cylinder) {
1058 
1059  const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(&targetSurface);
1060  double s[9] = { T(0, 3), T(1, 3), T(2, 3), T(0, 2),
1061  T(1, 2), T(2, 2), cyl->bounds().r(), (double)propagationDirection,
1062  0. };
1063  if (!propagateWithJacobianImpl(cache, errorPropagation, ty, s, cache.m_P, path)) {
1064  return nullptr;
1065  }
1066  }
1067 
1068  else if (ty == Trk::SurfaceType::Cone) {
1069 
1070  double k = static_cast<const Trk::ConeSurface*>(&targetSurface)->bounds().tanAlpha();
1071  k = k * k + 1.;
1072  double s[9] = { T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2), k, (double)propagationDirection, 0. };
1073  if (!propagateWithJacobianImpl(cache, errorPropagation, ty, s, cache.m_P, path)) {
1074  return nullptr;
1075  }
1076  }
1077 
1078  else if (ty == Trk::SurfaceType::Perigee) {
1079 
1080  double s[6] = { T(0, 3), T(1, 3), T(2, 3), 0., 0., 1. };
1081  if (!propagateWithJacobianImpl(cache, errorPropagation, ty, s, cache.m_P, path)) {
1082  return nullptr;
1083  }
1084  }
1085 
1086  else { // assume curvilinear
1087 
1088  double s[4];
1089  double d = T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
1090 
1091  if (d >= 0.) {
1092  s[0] = T(0, 2);
1093  s[1] = T(1, 2);
1094  s[2] = T(2, 2);
1095  s[3] = d;
1096  } else {
1097  s[0] = -T(0, 2);
1098  s[1] = -T(1, 2);
1099  s[2] = -T(2, 2);
1100  s[3] = -d;
1101  }
1102  if (!propagateWithJacobianImpl(cache, errorPropagation, ty, s, cache.m_P, path)) {
1103  return nullptr;
1104  }
1105  }
1106 
1107  if (propagationDirection * path < 0.) {
1108  return nullptr;
1109  }
1110 
1111  double localp[5];
1112  // output in curvilinear parameters
1113  if (returnCurv || ty == Trk::SurfaceType::Cone) {
1114 
1116  Amg::Vector3D gp(cache.m_P[0], cache.m_P[1], cache.m_P[2]);
1117 
1118  if (boundaryCheck && !targetSurface.isOnSurface(gp)) {
1119  return nullptr;
1120  }
1121 
1122  if (!errorPropagation || !trackParameters->covariance()) {
1123  return std::make_unique<Trk::CurvilinearParameters>(gp, localp[2], localp[3], localp[4]);
1124  }
1125 
1126  double useless[2];
1127  Trk::RungeKuttaUtils::transformGlobalToCurvilinear(true, cache.m_P, useless, Jacobian);
1128  AmgSymMatrix(5) measurementCovariance =
1129  Trk::RungeKuttaUtils::newCovarianceMatrix(Jacobian, *trackParameters->covariance());
1130 
1131  if (cache.m_matPropOK && (cache.m_multipleScattering || cache.m_straggling) && std::abs(path) > 0.)
1132  covarianceContribution(cache, trackParameters.get(), path, std::abs(1. / cache.m_P[6]), measurementCovariance);
1133 
1134  return std::make_unique<Trk::CurvilinearParameters>(
1135  gp, localp[2], localp[3], localp[4], std::move(measurementCovariance));
1136  }
1137 
1138  // Common transformation for all surfaces
1139  Trk::RungeKuttaUtils::transformGlobalToLocal(&targetSurface, errorPropagation, cache.m_P, localp, Jacobian);
1140 
1141  if (boundaryCheck) {
1142  Amg::Vector2D localPosition(localp[0], localp[1]);
1143  if (!targetSurface.insideBounds(localPosition, 0.)) {
1144  return nullptr;
1145  }
1146  }
1147 
1148  auto onTargetSurf =
1149  targetSurface.createUniqueTrackParameters(localp[0], localp[1], localp[2], localp[3], localp[4], std::nullopt);
1150 
1151  if (!errorPropagation || !trackParameters->covariance()) {
1152  return onTargetSurf;
1153  }
1154 
1155  // Errormatrix is included. Use Jacobian to calculate new covariance
1156  AmgSymMatrix(5) measurementCovariance =
1157  Trk::RungeKuttaUtils::newCovarianceMatrix(Jacobian, *trackParameters->covariance());
1158 
1159  // Calculate multiple scattering and straggling covariance contribution.
1160  if (cache.m_matPropOK && (cache.m_multipleScattering || cache.m_straggling) && std::abs(path) > 0.)
1161  covarianceContribution(cache, trackParameters.get(), path, onTargetSurf.get(), measurementCovariance);
1162 
1163  return targetSurface.createUniqueTrackParameters(
1164  localp[0], localp[1], localp[2], localp[3], localp[4], std::move(measurementCovariance));
1165 }
1166 
1167 } // end of anonymous namespace
1168 
1170 // Constructor
1172 
1173 Trk::STEP_Propagator::STEP_Propagator(const std::string& p, const std::string& n, const IInterface* t)
1174  : AthAlgTool(p, n, t)
1175 {
1176  declareInterface<Trk::IPropagator>(this);
1177 }
1178 
1180 // Destructor
1182 
1184 
1185 // Athena standard methods
1186 // initialize
1188 
1189  // Read handle for AtlasFieldCacheCondObj
1191 
1192  if (!m_materialEffects) { // override all material interactions
1193  m_multipleScattering = false;
1194  m_energyLoss = false;
1195  m_straggling = false;
1196  } else if (!m_energyLoss) { // override straggling
1197  m_straggling = false;
1198  }
1199 
1200  if (m_simulation && m_simMatUpdator.retrieve().isFailure()) {
1201  ATH_MSG_WARNING("Simulation mode requested but material updator not found - no brem photon emission.");
1202  }
1203 
1204  if (m_simulation) {
1205  // get the random generator serice
1206  ATH_CHECK(m_rndGenSvc.retrieve());
1207  m_rngWrapper = m_rndGenSvc->getEngine(this, m_randomEngineName);
1208  }
1209 
1210  return StatusCode::SUCCESS;
1211 }
1212 
1213 // finalize
1215  return StatusCode::SUCCESS;
1216 }
1217 
1219 std::unique_ptr<Trk::NeutralParameters> Trk::STEP_Propagator::propagate(const Trk::NeutralParameters&,
1220  const Trk::Surface&,
1222  const Trk::BoundaryCheck&,
1223  bool) const {
1224  ATH_MSG_WARNING("[STEP_Propagator] STEP_Propagator does not handle neutral track parameters."
1225  << "Use the StraightLinePropagator instead.");
1226  return nullptr;
1227 }
1228 
1230 // Main function for track parameters and covariance matrix propagation
1232 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagate(
1233  const EventContext& ctx, const Trk::TrackParameters& trackParameters, const Trk::Surface& targetSurface,
1234  Trk::PropDirection propagationDirection, const Trk::BoundaryCheck& boundaryCheck,
1235  const MagneticFieldProperties& magneticFieldProperties, ParticleHypothesis particle, bool returnCurv,
1236  const Trk::TrackingVolume* tVol) const {
1237 
1238  // ATH_MSG_WARNING( "[STEP_Propagator] enter 1");
1239 
1240  double Jacobian[25];
1241  Cache cache(ctx);
1242 
1243  // Get field cache object
1244  getFieldCacheObject(cache, ctx);
1245  setCacheFromProperties(cache);
1246  clearMaterialEffects(cache);
1247 
1248  // Check for tracking volume (materialproperties)
1249  cache.m_trackingVolume = tVol;
1250  cache.m_material = tVol;
1251  cache.m_matPropOK = tVol != nullptr;
1252 
1253  cache.m_matupd_lastmom = trackParameters.momentum().mag();
1254  cache.m_matupd_lastpath = 0.;
1255  cache.m_matdump_lastpath = 0.;
1256 
1257  // no identified intersections needed/ no material dump / no path cache
1258  cache.m_identifiedParameters = nullptr;
1259  cache.m_matstates = nullptr;
1260  cache.m_extrapolationCache = nullptr;
1261  cache.m_hitVector = nullptr;
1262 
1263  return propagateRungeKuttaImpl(cache, true, trackParameters, targetSurface, propagationDirection,
1264  magneticFieldProperties, particle, boundaryCheck, Jacobian, returnCurv);
1265 }
1266 
1268 // Main function for track parameters and covariance matrix propagation
1269 // with search of closest surface (ST)
1271 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagate(
1272  const EventContext& ctx, const Trk::TrackParameters& trackParameters,
1273  std::vector<DestSurf>& targetSurfaces, Trk::PropDirection propagationDirection,
1274  const Trk::MagneticFieldProperties& magneticFieldProperties, ParticleHypothesis particle,
1275  std::vector<unsigned int>& solutions, double& path, bool usePathLimit, bool returnCurv,
1276  const Trk::TrackingVolume* tVol) const {
1277 
1278  Cache cache(ctx);
1279 
1280  // Get field cache object
1281  getFieldCacheObject(cache, ctx);
1282  setCacheFromProperties(cache);
1283  clearMaterialEffects(cache);
1284 
1285  // Check for tracking volume (materialproperties)
1286  cache.m_trackingVolume = tVol;
1287  cache.m_material = tVol;
1288  cache.m_matPropOK = tVol != nullptr;
1289 
1290  // no identified intersections needed/ no material dump
1291  cache.m_identifiedParameters = nullptr;
1292  cache.m_matstates = nullptr;
1293  cache.m_extrapolationCache = nullptr;
1294  cache.m_hitVector = nullptr;
1295 
1296  cache.m_matupd_lastmom = trackParameters.momentum().mag();
1297  cache.m_matupd_lastpath = 0.;
1298  cache.m_matdump_lastpath = 0.;
1299 
1300  // resolve path limit input
1301  if (path > 0.) {
1302  cache.m_propagateWithPathLimit = usePathLimit ? 1 : 0;
1303  cache.m_pathLimit = path;
1304  path = 0.;
1305  } else {
1306  cache.m_propagateWithPathLimit = 0;
1307  cache.m_pathLimit = -1.;
1308  path = 0.;
1309  }
1311  return propagateNeutral(trackParameters, targetSurfaces, propagationDirection, solutions, path,
1312  usePathLimit, returnCurv);
1313 
1314  return propagateRungeKutta(cache, true, trackParameters, targetSurfaces, propagationDirection,
1315  magneticFieldProperties, particle, solutions, path, returnCurv);
1316 }
1317 
1319 // Main function for track parameters and covariance matrix propagation
1320 // with search of closest surface and time info (ST)
1322 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagateT(
1323  const EventContext& ctx, const Trk::TrackParameters& trackParameters,
1324  std::vector<DestSurf>& targetSurfaces, Trk::PropDirection propagationDirection,
1325  const Trk::MagneticFieldProperties& magneticFieldProperties, ParticleHypothesis particle,
1326  std::vector<unsigned int>& solutions, PathLimit& pathLim, TimeLimit& timeLim, bool returnCurv,
1327  const Trk::TrackingVolume* tVol, std::vector<Trk::HitInfo>*& hitVector) const {
1328 
1329  Cache cache(ctx);
1330 
1331  // Get field cache object
1332  getFieldCacheObject(cache, ctx);
1333  setCacheFromProperties(cache);
1334  clearMaterialEffects(cache);
1335 
1336  // cache particle mass
1337  cache.m_particleMass = Trk::ParticleMasses::mass[particle]; // Get particle mass from ParticleHypothesis
1338 
1339  // cache input timing - for secondary track emission
1340  cache.m_timeIn = timeLim.time;
1341 
1342  // Check for tracking volume (materialproperties)
1343  cache.m_trackingVolume = tVol;
1344  cache.m_material = tVol;
1345  cache.m_matPropOK = tVol != nullptr;
1346 
1347  // no identified intersections needed/ no material dump
1348  cache.m_identifiedParameters = nullptr;
1349  cache.m_matstates = nullptr;
1350  cache.m_extrapolationCache = nullptr;
1351  cache.m_hitVector = hitVector;
1352 
1353  cache.m_matupd_lastmom = trackParameters.momentum().mag();
1354  cache.m_matupd_lastpath = 0.;
1355  cache.m_matdump_lastpath = 0.;
1356 
1357  // convert time/path limits into trajectory limit (in mm)
1358  double dMat = pathLim.x0Max - pathLim.x0Collected;
1359  double path =
1360  dMat > 0 && cache.m_matPropOK && cache.m_material->x0() > 0. ? dMat * cache.m_material->x0() : -1.;
1361 
1362  double dTim = timeLim.tMax - timeLim.time;
1363  double beta = 1.;
1364  if (dTim > 0.) {
1365  double mom = trackParameters.momentum().mag();
1366  beta = mom / std::sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass);
1367  }
1368  double timMax = dTim > 0 ? dTim * beta * Gaudi::Units::c_light : -1.;
1369 
1370  if (timMax > 0. && timMax < path)
1371  path = timMax;
1372  bool usePathLimit = (path > 0.);
1373 
1374  // resolve path limit input
1375  if (path > 0.) {
1376  cache.m_propagateWithPathLimit = usePathLimit ? 1 : 0;
1377  cache.m_pathLimit = path;
1378  path = 0.;
1379  } else {
1380  cache.m_propagateWithPathLimit = 0;
1381  cache.m_pathLimit = -1.;
1382  path = 0.;
1383  }
1384 
1385  std::unique_ptr<Trk::TrackParameters> nextPar{};
1386 
1388  nextPar = propagateNeutral(trackParameters, targetSurfaces, propagationDirection, solutions, path,
1389  usePathLimit, returnCurv);
1390  } else {
1391  nextPar = propagateRungeKutta(cache, true, trackParameters, targetSurfaces, propagationDirection,
1392  magneticFieldProperties, particle, solutions, path, returnCurv);
1393  }
1394  // update material path
1395  if (cache.m_matPropOK && cache.m_material->x0() > 0. && path > 0.) {
1396  pathLim.updateMat(path / cache.m_material->x0(), cache.m_material->averageZ(), 0.);
1397  }
1398  // return value
1399  timeLim.time += cache.m_timeOfFlight;
1400  return nextPar;
1401 }
1402 
1404 // Main function for track parameters and covariance matrix propagation
1405 // with search of closest surface and material collection (ST)
1407 
1408 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagateM(
1409  const EventContext& ctx, const Trk::TrackParameters& trackParameters,
1410  std::vector<DestSurf>& targetSurfaces, Trk::PropDirection propagationDirection,
1411  const Trk::MagneticFieldProperties& magneticFieldProperties, ParticleHypothesis particle,
1412  std::vector<unsigned int>& solutions, std::vector<const Trk::TrackStateOnSurface*>* matstates,
1413  std::vector<std::pair<std::unique_ptr<Trk::TrackParameters>, int>>* intersections, double& path,
1414  bool usePathLimit, bool returnCurv, const Trk::TrackingVolume* tVol,
1415  Trk::ExtrapolationCache* extrapCache) const {
1416 
1417  Cache cache(ctx);
1418 
1419  // Get field cache object
1420  getFieldCacheObject(cache, ctx);
1421  setCacheFromProperties(cache);
1422  clearMaterialEffects(cache);
1423 
1424  // Check for tracking volume (materialproperties)
1425  cache.m_trackingVolume = tVol;
1426  cache.m_material = tVol;
1427  cache.m_matPropOK = tVol != nullptr;
1428 
1429  cache.m_matstates = matstates;
1430  cache.m_identifiedParameters = intersections;
1431  cache.m_extrapolationCache = extrapCache;
1432  cache.m_hitVector = nullptr;
1433 
1434  cache.m_matupd_lastmom = trackParameters.momentum().mag();
1435  cache.m_matupd_lastpath = 0.;
1436  cache.m_matdump_lastpath = 0.;
1437  cache.m_extrapolationCache = extrapCache;
1438 
1439  // switch on the detailed energy loss
1440  if (cache.m_extrapolationCache) {
1441  cache.m_detailedElossFlag = true;
1442  }
1443  // resolve path limit input
1444  if (path > 0.) {
1445  cache.m_propagateWithPathLimit = usePathLimit ? 1 : 0;
1446  cache.m_pathLimit = path;
1447  path = 0.;
1448  } else {
1449  cache.m_propagateWithPathLimit = 0;
1450  cache.m_pathLimit = -1.;
1451  path = 0.;
1452  }
1454  return propagateNeutral(trackParameters, targetSurfaces, propagationDirection, solutions, path,
1455  usePathLimit, returnCurv);
1456  }
1457  return propagateRungeKutta(cache, true, trackParameters, targetSurfaces, propagationDirection,
1458  magneticFieldProperties, particle, solutions, path, returnCurv);
1459 }
1460 
1462 // Main function for track parameters and covariance matrix propagation.
1464 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagate(
1465  const EventContext& ctx,
1466  const Trk::TrackParameters& trackParameters,
1467  const Trk::Surface& targetSurface,
1468  Trk::PropDirection propagationDirection,
1469  const Trk::BoundaryCheck& boundaryCheck,
1470  const Trk::MagneticFieldProperties& magneticFieldProperties,
1471  std::optional<Trk::TransportJacobian>& jacobian,
1472  double&,
1474  bool returnCurv,
1475  const Trk::TrackingVolume* tVol) const {
1476 
1477  double Jacobian[25];
1478  Cache cache(ctx);
1479 
1480  // Get field cache object
1481  getFieldCacheObject(cache, ctx);
1482  setCacheFromProperties(cache);
1483  clearMaterialEffects(cache);
1484 
1485  // Check for tracking volume (materialproperties)
1486  cache.m_trackingVolume = tVol;
1487  cache.m_material = tVol;
1488  cache.m_matPropOK = tVol != nullptr;
1489 
1490  // no identified intersections needed/ no material dump
1491  cache.m_identifiedParameters = nullptr;
1492  cache.m_matstates = nullptr;
1493  cache.m_extrapolationCache = nullptr;
1494  cache.m_hitVector = nullptr;
1495 
1496  cache.m_matupd_lastmom = trackParameters.momentum().mag();
1497  cache.m_matupd_lastpath = 0.;
1498  cache.m_matdump_lastpath = 0.;
1499 
1500  std::unique_ptr<Trk::TrackParameters> parameters =
1501  propagateRungeKuttaImpl(cache, true, trackParameters, targetSurface, propagationDirection,
1502  magneticFieldProperties, particle, boundaryCheck, Jacobian, returnCurv);
1503 
1504  if (parameters) {
1505  Jacobian[24] = Jacobian[20];
1506  Jacobian[23] = 0.;
1507  Jacobian[22] = 0.;
1508  Jacobian[21] = 0.;
1509  Jacobian[20] = 0.;
1510  jacobian = std::make_optional<Trk::TransportJacobian>(Jacobian);
1511  } else {
1512  jacobian.reset();
1513  }
1514 
1515  return parameters;
1516 }
1517 
1519 // Main function for track parameters propagation without covariance matrix
1521 
1522 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagateParameters(
1523  const EventContext& ctx, const Trk::TrackParameters& trackParameters, const Trk::Surface& targetSurface,
1524  Trk::PropDirection propagationDirection, const Trk::BoundaryCheck& boundaryCheck,
1525  const Trk::MagneticFieldProperties& magneticFieldProperties, ParticleHypothesis particle, bool returnCurv,
1526  const Trk::TrackingVolume* tVol) const {
1527 
1528  double Jacobian[25];
1529  Cache cache(ctx);
1530 
1531  // Get field cache object
1532  getFieldCacheObject(cache, ctx);
1533  setCacheFromProperties(cache);
1534  clearMaterialEffects(cache);
1535 
1536  // Check for tracking volume (materialproperties)
1537  cache.m_trackingVolume = tVol;
1538  cache.m_material = tVol;
1539  cache.m_matPropOK = tVol != nullptr;
1540 
1541  // no identified intersections needed/ no material dump
1542  cache.m_identifiedParameters = nullptr;
1543  cache.m_matstates = nullptr;
1544  cache.m_hitVector = nullptr;
1545 
1546  return propagateRungeKuttaImpl(cache, false, trackParameters, targetSurface, propagationDirection,
1547  magneticFieldProperties, particle, boundaryCheck, Jacobian, returnCurv);
1548 }
1549 
1551 // Main function for track parameters propagation without covariance matrix.
1553 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagateParameters(
1554  const EventContext& ctx,
1555  const Trk::TrackParameters& trackParameters,
1556  const Trk::Surface& targetSurface,
1557  Trk::PropDirection propagationDirection,
1558  const Trk::BoundaryCheck& boundaryCheck,
1559  const Trk::MagneticFieldProperties& magneticFieldProperties,
1560  std::optional<Trk::TransportJacobian>& jacobian,
1562  bool returnCurv,
1563  const Trk::TrackingVolume* tVol) const {
1564 
1565  double Jacobian[25];
1566 
1567  Cache cache(ctx);
1568 
1569  // Get field cache object
1570  getFieldCacheObject(cache, ctx);
1571  setCacheFromProperties(cache);
1572  clearMaterialEffects(cache);
1573 
1574  // Check for tracking volume (materialproperties)
1575  cache.m_trackingVolume = tVol;
1576  cache.m_material = tVol;
1577  cache.m_matPropOK = tVol != nullptr;
1578 
1579  // no identified intersections needed/ no material dump
1580  cache.m_identifiedParameters = nullptr;
1581  cache.m_matstates = nullptr;
1582  cache.m_extrapolationCache = nullptr;
1583  cache.m_hitVector = nullptr;
1584 
1585  std::unique_ptr<Trk::TrackParameters> parameters =
1586  propagateRungeKuttaImpl(cache, true, trackParameters, targetSurface, propagationDirection,
1587  magneticFieldProperties, particle, boundaryCheck, Jacobian, returnCurv);
1588 
1589  if (parameters) {
1590  Jacobian[24] = Jacobian[20];
1591  Jacobian[23] = 0.;
1592  Jacobian[22] = 0.;
1593  Jacobian[21] = 0.;
1594  Jacobian[20] = 0.;
1595  jacobian = std::make_optional<Trk::TransportJacobian>(Jacobian);
1596  } else {
1597  jacobian.reset();
1598  }
1599 
1600  return parameters;
1601 }
1602 
1604 // Function for finding the intersection point with a surface
1606 std::optional<Trk::TrackSurfaceIntersection> Trk::STEP_Propagator::intersect(
1607  const EventContext& ctx,
1608  const Trk::TrackParameters& trackParameters,
1609  const Trk::Surface& targetSurface,
1610  const Trk::MagneticFieldProperties& mft,
1612  const Trk::TrackingVolume* tVol) const {
1613 
1614  Cache cache(ctx);
1615 
1616  // Get field cache object
1617  getFieldCacheObject(cache, ctx);
1618  setCacheFromProperties(cache);
1619  clearMaterialEffects(cache);
1620 
1621  cache.m_particle = particle; // Store for later use
1622 
1623  // Check for tracking volume (materialproperties)
1624  cache.m_trackingVolume = tVol;
1625  cache.m_material = tVol;
1626  cache.m_matPropOK = tVol != nullptr;
1627 
1628  // no identified intersections needed/ no material dump
1629  cache.m_identifiedParameters = nullptr;
1630  cache.m_matstates = nullptr;
1631  cache.m_extrapolationCache = nullptr;
1632  cache.m_hitVector = nullptr;
1633 
1634  // Bfield mode
1635  mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1636 
1637  // Check inputvalues
1638  if (cache.m_tolerance <= 0.) {
1639  return std::nullopt;
1640  }
1641  if (cache.m_momentumCutOff < 0.) {
1642  return std::nullopt;
1643  }
1644  if (std::abs(1. / trackParameters.parameters()[Trk::qOverP]) <= cache.m_momentumCutOff) {
1645  return std::nullopt;
1646  }
1647 
1648  // Check for empty volumes. If x != x then x is not a number.
1649  if (cache.m_matPropOK && ((cache.m_material->zOverAtimesRho() == 0.) || (cache.m_material->x0() == 0.) ||
1650  (cache.m_material->zOverAtimesRho() != cache.m_material->zOverAtimesRho()))) {
1651  cache.m_matPropOK = false;
1652  }
1653 
1654  // double P[45];
1655  if (!Trk::RungeKuttaUtils::transformLocalToGlobal(false, trackParameters, cache.m_P)) {
1656  return std::nullopt;
1657  }
1658  double path = 0.;
1659 
1660  const Amg::Transform3D& T = targetSurface.transform();
1661  Trk::SurfaceType ty = targetSurface.type();
1662 
1663  if (ty == Trk::SurfaceType::Plane || ty == Trk::SurfaceType::Disc) {
1664  double s[4];
1665  double d = T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
1666 
1667  if (d >= 0.) {
1668  s[0] = T(0, 2);
1669  s[1] = T(1, 2);
1670  s[2] = T(2, 2);
1671  s[3] = d;
1672  } else {
1673  s[0] = -T(0, 2);
1674  s[1] = -T(1, 2);
1675  s[2] = -T(2, 2);
1676  s[3] = -d;
1677  }
1678  if (!propagateWithJacobianImpl(cache, false, ty, s, cache.m_P, path)){
1679  return std::nullopt;
1680  }
1681  }
1682 
1683  else if (ty == Trk::SurfaceType::Line) {
1684 
1685  double s[6] = {T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2)};
1686  if (!propagateWithJacobianImpl(cache, false, ty, s, cache.m_P, path)) {
1687  return std::nullopt;
1688  }
1689  }
1690 
1691  else if (ty == Trk::SurfaceType::Cylinder) {
1692 
1693  const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(&targetSurface);
1694  double s[9] = {
1695  T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2), cyl->bounds().r(), Trk::alongMomentum, 0.};
1696  if (!propagateWithJacobianImpl(cache, false, ty, s, cache.m_P, path)) {
1697  return std::nullopt;
1698  }
1699  }
1700 
1701  else if (ty == Trk::SurfaceType::Cone) {
1702 
1703  double k = static_cast<const Trk::ConeSurface*>(&targetSurface)->bounds().tanAlpha();
1704  k = k * k + 1.;
1705  double s[9] = {T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2), k, Trk::alongMomentum, 0.};
1706  if (!propagateWithJacobianImpl(cache, false, ty, s, cache.m_P, path)) {
1707  return std::nullopt;
1708  }
1709  }
1710 
1711  else if (ty == Trk::SurfaceType::Perigee) {
1712 
1713  double s[6] = {T(0, 3), T(1, 3), T(2, 3), 0., 0., 1.};
1714  if (!propagateWithJacobianImpl(cache, false, ty, s, cache.m_P, path)) {
1715  return std::nullopt;
1716  }
1717  }
1718 
1719  else { // presumably curvilinear
1720 
1721  double s[4];
1722  double d = T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
1723 
1724  if (d >= 0.) {
1725  s[0] = T(0, 2);
1726  s[1] = T(1, 2);
1727  s[2] = T(2, 2);
1728  s[3] = d;
1729  } else {
1730  s[0] = -T(0, 2);
1731  s[1] = -T(1, 2);
1732  s[2] = -T(2, 2);
1733  s[3] = -d;
1734  }
1735  if (!propagateWithJacobianImpl(cache, false, ty, s, cache.m_P, path)) {
1736  return std::nullopt;
1737  }
1738  }
1739 
1740  Amg::Vector3D globalPosition(cache.m_P[0], cache.m_P[1], cache.m_P[2]);
1741  Amg::Vector3D direction(cache.m_P[3], cache.m_P[4], cache.m_P[5]);
1742  return std::make_optional<Trk::TrackSurfaceIntersection>(globalPosition, direction, path);
1743 }
1744 
1745 std::optional<Trk::TrackSurfaceIntersection>
1747  const EventContext& ctx, const Trk::Surface& surface,
1748  const Trk::TrackSurfaceIntersection& trackIntersection, const double qOverP,
1749  const Trk::MagneticFieldProperties& mft,
1750  ParticleHypothesis particle) const {
1751 
1752  const Amg::Vector3D& origin = trackIntersection.position();
1753  const Amg::Vector3D& direction = trackIntersection.direction();
1754 
1755  auto perigeeSurface = PerigeeSurface(origin);
1756  perigeeSurface.setOwner(Trk::userOwn); // tmp ones
1757 
1758  auto tmpTrackParameters =
1759  Trk::Perigee(0., 0., direction.phi(), direction.theta(), qOverP, perigeeSurface, std::nullopt);
1760 
1761  std::optional<Trk::TrackSurfaceIntersection> solution =
1762  qOverP == 0
1763  ? intersect(ctx, tmpTrackParameters, surface,
1765  : intersect(ctx, tmpTrackParameters, surface, mft, particle, nullptr);
1766 
1767  return solution;
1768 }
1769 
1771 // Global positions calculation inside CylinderBounds
1772 // if max step allowed > 0 -> propagate along momentum else propagate opposite momentum
1774 void Trk::STEP_Propagator::globalPositions(const EventContext& ctx, std::deque<Amg::Vector3D>& positionsList,
1775  const Trk::TrackParameters& trackParameters,
1776  const Trk::MagneticFieldProperties& mft,
1777  const Trk::CylinderBounds& cylinderBounds, double maxStepSize,
1779  const Trk::TrackingVolume* tVol) const {
1780  Cache cache(ctx);
1781 
1782  // Get field cache object
1783  getFieldCacheObject(cache, ctx);
1784  setCacheFromProperties(cache);
1785  clearMaterialEffects(cache);
1786 
1787  cache.m_particle = particle; // Store for later use
1788 
1789  // Check for tracking volume (materialproperties)
1790  cache.m_trackingVolume = tVol;
1791  cache.m_material = tVol;
1792  cache.m_matPropOK = tVol != nullptr;
1793 
1794  // Check for empty volumes. If x != x then x is not a number.
1795  if (cache.m_matPropOK && ((cache.m_material->zOverAtimesRho() == 0.) || (cache.m_material->x0() == 0.) ||
1796  (cache.m_material->zOverAtimesRho() != cache.m_material->zOverAtimesRho()))) {
1797  cache.m_matPropOK = false;
1798  }
1799 
1800  mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1801 
1802  // Check inputvalues
1803  if (cache.m_tolerance <= 0.)
1804  return;
1805 
1806  double PP[7];
1807  if (!Trk::RungeKuttaUtils::transformLocalToGlobal(false, trackParameters, PP))
1808  return;
1809 
1810  double maxPath = cache.m_maxPath; // Max path allowed
1811  double dDir[3] = {0., 0., 0.}; // Start directions derivs. Zero in case of no RK steps
1812  double distanceStepped = 0.;
1813  double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
1814  // dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz
1815  bool firstStep = true; // Poll B1, else recycle B4
1816  double path = 0.; // path of the trajectory
1817  double radius2Max = cylinderBounds.r() * cylinderBounds.r(); // max. radius**2 of region
1818  double zMax = cylinderBounds.halflengthZ(); // max. Z of region
1819  double radius2 = PP[0] * PP[0] + PP[1] * PP[1]; // Start radius**2
1820  double direction = PP[0] * PP[3] + PP[1] * PP[4]; // Direction
1821  double h = maxStepSize; // max step allowed
1822 
1823  // Test position of the track
1824  if ((std::abs(PP[2]) > zMax) || (radius2 > radius2Max))
1825  return;
1826 
1827  // Store initial position
1828  Amg::Vector3D initialPosition(PP[0], PP[1], PP[2]);
1829  positionsList.push_back(initialPosition);
1830 
1831  bool perigee = false;
1832  if (std::abs(direction) < 0.00001) {
1833  perigee = true;
1834  }
1835 
1836  for (int i = 0; i != 2; ++i) {
1837  if (i) {
1838  if (perigee)
1839  return;
1840  h = -h;
1841  }
1842  double p[7] = {PP[0], PP[1], PP[2], PP[3], PP[4], PP[5], PP[6]};
1843 
1844  while (std::abs(path) < maxPath) {
1845  // Do the step.
1846  if (!rungeKuttaStep(cache, false, h, p, dDir, BG1, firstStep, distanceStepped))
1847  break;
1848  path = path + distanceStepped;
1849 
1850  // Keep h within max stepsize
1851  if (h > maxStepSize) {
1852  h = maxStepSize;
1853  } else if (h < -maxStepSize) {
1854  h = -maxStepSize;
1855  }
1856 
1857  // store current step
1858  Amg::Vector3D globalPosition(p[0], p[1], p[2]);
1859  if (!i) {
1860  positionsList.push_back(globalPosition);
1861  } else {
1862  positionsList.push_front(globalPosition);
1863  }
1864 
1865  // Test position of the track
1866  radius2 = p[0] * p[0] + p[1] * p[1];
1867  if ((std::abs(p[2]) > zMax) || (radius2 > radius2Max))
1868  break;
1869 
1870  // Test perigee
1871  if ((p[0] * p[3] + p[1] * p[4]) * direction < 0.) {
1872  if (i)
1873  return;
1874  perigee = true;
1875  }
1876  }
1877  }
1878 }
1879 
1881 // Main function for track propagation with or without jacobian
1882 // with search of closest surface
1884 std::unique_ptr<Trk::TrackParameters> Trk::STEP_Propagator::propagateRungeKutta(
1885  Cache& cache, bool errorPropagation, const Trk::TrackParameters& inputTrackParameters,
1886  std::vector<DestSurf>& targetSurfaces, Trk::PropDirection propagationDirection,
1888  std::vector<unsigned int>& solutions, double& totalPath, bool returnCurv) const {
1889  // Store for later use
1890  cache.m_particle = particle;
1891  cache.m_charge = inputTrackParameters.charge();
1892  cache.m_inputThetaVariance = 0.;
1893 
1894  std::unique_ptr<Trk::TrackParameters> trackParameters{};
1895 
1896  // Bfield mode
1897  mft.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1898 
1899  // Check inputvalues
1900  if (cache.m_tolerance <= 0.)
1901  return nullptr;
1902  if (cache.m_momentumCutOff < 0.)
1903  return nullptr;
1904 
1905  // Set momentum to 1e10 (straight line) and charge to + if q/p is zero
1906  if (inputTrackParameters.parameters()[Trk::qOverP] == 0) {
1907  trackParameters = createStraightLine(&inputTrackParameters);
1908  if (!trackParameters) {
1909  return nullptr;
1910  }
1911  } else {
1912  trackParameters.reset(inputTrackParameters.clone());
1913  }
1914 
1915  if (std::abs(1. / trackParameters->parameters()[Trk::qOverP]) <= cache.m_momentumCutOff) {
1916  return nullptr;
1917  }
1918 
1919  // Check for empty volumes. If x != x then x is not a number.
1920  if (cache.m_matPropOK && ((cache.m_material->zOverAtimesRho() == 0.) || (cache.m_material->x0() == 0.) ||
1921  (cache.m_material->zOverAtimesRho() != cache.m_material->zOverAtimesRho()))) {
1922  cache.m_matPropOK = false;
1923  }
1924 
1925  if (errorPropagation && !trackParameters->covariance()) {
1926  errorPropagation = false;
1927  }
1928 
1929  if (cache.m_matPropOK && errorPropagation && cache.m_straggling)
1930  cache.m_stragglingVariance = 0.;
1931  cache.m_combinedCovariance.setZero();
1932  cache.m_covariance.setZero();
1933 
1934  if (errorPropagation || cache.m_matstates) {
1935  // this needs debugging
1936  cache.m_inputThetaVariance = trackParameters->covariance() ? (*trackParameters->covariance())(3, 3) : 0.;
1937  cache.m_combinedEloss.set(0., 0., 0., 0., 0., 0.);
1938  cache.m_combinedThickness = 0.;
1939  }
1940 
1941  // double P[45]; // Track parameters and jacobian
1942  if (!Trk::RungeKuttaUtils::transformLocalToGlobal(errorPropagation, *trackParameters, cache.m_P)) {
1943  return nullptr;
1944  }
1945 
1946  double path = 0.;
1947 
1948  // activate brem photon emission if required
1949  cache.m_brem = m_simulation && particle == Trk::electron && m_simMatUpdator;
1950 
1951  // loop while valid solutions
1952  bool validStep = true;
1953  totalPath = 0.;
1954  cache.m_timeOfFlight = 0.;
1955  // Common transformation for all surfaces (angles and momentum)
1956  double localp[5];
1957  double Jacobian[21];
1958  while (validStep) {
1959  // propagation to next surface
1960  validStep = propagateWithJacobian(cache, errorPropagation, targetSurfaces, cache.m_P,
1961  propagationDirection, solutions, path, totalPath);
1962  if (!validStep) {
1963  return nullptr;
1964  }
1965  if (propagationDirection * path <= 0.) {
1966  return nullptr;
1967  }
1968  totalPath += path;
1969  cache.m_timeOfFlight += cache.m_timeStep;
1970  if (cache.m_propagateWithPathLimit > 1 || cache.m_binMat) {
1971  // make sure that for sliding surfaces the result does not get distorted
1972  // return curvilinear parameters
1973  std::unique_ptr<Trk::CurvilinearParameters> cPar = nullptr;
1975  if (!errorPropagation) {
1976  cPar = std::make_unique<Trk::CurvilinearParameters>(
1977  Amg::Vector3D(cache.m_P[0], cache.m_P[1], cache.m_P[2]), localp[2], localp[3], localp[4]);
1978  } else {
1979  double useless[2];
1980  Trk::RungeKuttaUtils::transformGlobalToCurvilinear(true, cache.m_P, useless, Jacobian);
1981  AmgSymMatrix(5) measurementCovariance =
1982  Trk::RungeKuttaUtils::newCovarianceMatrix(Jacobian, *trackParameters->covariance());
1983  // Calculate multiple scattering and straggling covariance contribution.
1984  if (cache.m_matPropOK && (cache.m_multipleScattering || cache.m_straggling) &&
1985  std::abs(totalPath) > 0.) {
1986  covarianceContribution(cache, trackParameters.get(), totalPath, std::abs(1. / cache.m_P[6]),
1987  measurementCovariance);
1988  }
1989  cPar = std::make_unique<Trk::CurvilinearParameters>(
1990  Amg::Vector3D(cache.m_P[0], cache.m_P[1], cache.m_P[2]), localp[2], localp[3], localp[4],
1991  std::move(measurementCovariance));
1992  }
1993  // material collection : first iteration, bin material averaged
1994  // collect material
1995  if (cache.m_binMat && (cache.m_matstates || (errorPropagation && cache.m_extrapolationCache)) &&
1996  std::abs(totalPath - cache.m_matdump_lastpath) > 1.) {
1997  dumpMaterialEffects(cache, cPar.get(), totalPath);
1998  }
1999  return cPar;
2000  }
2001  if (cache.m_propagateWithPathLimit > 0)
2002  cache.m_pathLimit -= path;
2003  // boundary check
2004  // take into account that there may be many identical surfaces with different boundaries
2005  Amg::Vector3D gp(cache.m_P[0], cache.m_P[1], cache.m_P[2]);
2006  bool solution = false;
2007  std::vector<unsigned int> valid_solutions;
2008  valid_solutions.reserve(solutions.size());
2009 
2010  std::vector<unsigned int>::iterator iSol = solutions.begin();
2011  while (iSol != solutions.end()) {
2012  if (targetSurfaces[*iSol].first->isOnSurface(gp, targetSurfaces[*iSol].second, 0.001, 0.001)) {
2013  if (!solution) {
2015  if (returnCurv || targetSurfaces[*iSol].first->type() == Trk::SurfaceType::Cone) {
2016  Trk::RungeKuttaUtils::transformGlobalToCurvilinear(errorPropagation, cache.m_P, localp, Jacobian);
2017  } else
2018  Trk::RungeKuttaUtils::transformGlobalToLocal(targetSurfaces[*iSol].first, errorPropagation,
2019  cache.m_P, localp, Jacobian);
2020  solution = true;
2021  }
2022  valid_solutions.push_back(*iSol);
2023  }
2024  ++iSol;
2025  }
2026  solutions = std::move(valid_solutions);
2027  if (solution)
2028  break;
2029  }
2030 
2031  if (solutions.empty()) {
2032  return nullptr;
2033  }
2034 
2035  // simulation mode : smear momentum
2036  if (m_simulation && cache.m_matPropOK) {
2037  double radDist = totalPath / cache.m_material->x0();
2038  smear(cache, localp[2], localp[3], trackParameters.get(), radDist);
2039  }
2040 
2041  std::unique_ptr<Trk::TrackParameters> onTargetSurf =
2042  (returnCurv || targetSurfaces[solutions[0]].first->type() == Trk::SurfaceType::Cone)
2043  ? nullptr
2044  : targetSurfaces[solutions[0]].first->createUniqueTrackParameters(
2045  localp[0], localp[1], localp[2], localp[3], localp[4], std::nullopt);
2046 
2047  if (!errorPropagation) {
2048  if (returnCurv || targetSurfaces[solutions[0]].first->type() == Trk::SurfaceType::Cone) {
2049  Amg::Vector3D gp(cache.m_P[0], cache.m_P[1], cache.m_P[2]);
2050  return std::make_unique<Trk::CurvilinearParameters>(gp, localp[2], localp[3], localp[4]);
2051  }
2052  return onTargetSurf;
2053  }
2054 
2055  // Errormatrix is included. Use Jacobian to calculate new covariance
2057  for (double i : Jacobian) {
2058  if (!Amg::saneCovarianceElement(i)) {
2059  return nullptr;
2060  }
2061  }
2062 
2063  AmgSymMatrix(5) measurementCovariance =
2064  Trk::RungeKuttaUtils::newCovarianceMatrix(Jacobian, *trackParameters->covariance());
2065  if (!Amg::hasPositiveOrZeroDiagElems(measurementCovariance))
2066  return nullptr;
2067 
2068  // Calculate multiple scattering and straggling covariance contribution.
2069  if (cache.m_matPropOK && (cache.m_multipleScattering || cache.m_straggling) && std::abs(totalPath) > 0.) {
2070  if (returnCurv || targetSurfaces[solutions[0]].first->type() == Trk::SurfaceType::Cone) {
2071  covarianceContribution(cache, trackParameters.get(), totalPath, std::abs(1. / cache.m_P[6]),
2072  measurementCovariance);
2073  } else {
2074  covarianceContribution(cache, trackParameters.get(), totalPath, onTargetSurf.get(),
2075  measurementCovariance);
2076  }
2077  }
2078 
2079  if (returnCurv || targetSurfaces[solutions[0]].first->type() == Trk::SurfaceType::Cone) {
2080  Amg::Vector3D gp(cache.m_P[0], cache.m_P[1], cache.m_P[2]);
2081  return std::make_unique<Trk::CurvilinearParameters>(gp, localp[2], localp[3], localp[4],
2082  std::move(measurementCovariance));
2083  }
2084 
2085  // delete onTargetSurf; // the covariance matrix can be just added instead of recreating ?
2086  return targetSurfaces[solutions[0]].first->createUniqueTrackParameters(
2087  localp[0], localp[1], localp[2], localp[3], localp[4], std::move(measurementCovariance));
2088 }
2089 
2091 // Runge Kutta main program for propagation with or without Jacobian
2092 // with search of closest surface (ST)
2094 bool Trk::STEP_Propagator::propagateWithJacobian(Cache& cache, bool errorPropagation,
2095  std::vector<DestSurf>& sfs, double* P,
2096  Trk::PropDirection propDir,
2097  std::vector<unsigned int>& solutions, double& path,
2098  double sumPath) const {
2099  double maxPath = cache.m_maxPath; // Max path allowed
2100  double* pos = &P[0]; // Start coordinates
2101  double* dir = &P[3]; // Start directions
2102  double dDir[3] = {0., 0., 0.}; // Start directions derivs. Zero in case of no RK steps
2103  // int targetPassed = 0; // how many times have we passed the target?
2104  double previousDistance = 0.;
2105  double distanceStepped = 0.;
2106  double BG1[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; // Bx, By, Bz, dBx/dx, dBx/dy, dBx/dz,
2107  // dBy/dx, dBy/dy, dBy/dz, dBz/dx, dBz/dy, dBz/dz at first point
2108  bool firstStep = true; // Poll BG1, else recycle BG4
2109  int steps = 0;
2110  path = 0.; // path of the trajectory
2111  cache.m_timeStep = 0.; // time separation corresponding to the trajectory
2112  double mom = 0.; // need momentum and beta for timing
2113  double beta = 1.; // need momentum and beta for timing
2114 
2115  // factor to stabilize iteration for soft tracks
2116  double helpSoft = 1.;
2117 
2118  // limit number of recovery attempts
2119  int restartLimit = 10;
2120 
2121  Amg::Vector3D position(P[0], P[1], P[2]);
2122  Amg::Vector3D direction0(P[3], P[4], P[5]);
2123 
2124  // binned material ?
2125  cache.m_binMat = nullptr;
2126  if (cache.m_trackingVolume && cache.m_trackingVolume->isAlignable()) {
2127  const Trk::AlignableTrackingVolume* aliTV =
2128  static_cast<const Trk::AlignableTrackingVolume*>(cache.m_trackingVolume);
2129  cache.m_binMat = aliTV->binnedMaterial();
2130  }
2131 
2132  // closest distance estimate
2133  // maintain count of oscilations and previous distance for each surface;
2134  // skip initial trivial solutions (input parameters at surface) - should be treated before call to the
2135  // propagator
2136  // !
2137  double tol = 0.001;
2138  solutions.clear();
2139  double distanceToTarget = propDir * maxPath;
2140  cache.m_currentDist.resize(sfs.size()); // keep size through the call
2141 
2142  int nextSf = sfs.size();
2143  int nextSfCand = nextSf;
2144  std::vector<DestSurf>::iterator sIter = sfs.begin();
2145  std::vector<DestSurf>::iterator sBeg = sfs.begin();
2146  unsigned int numSf = 0;
2147  unsigned int iCurr = 0; // index for m_currentDist
2148  int startSf = -99;
2149  for (; sIter != sfs.end(); ++sIter) {
2150  Trk::DistanceSolution distSol = (*sIter).first->straightLineDistanceEstimate(position, direction0);
2151  double distEst = -propDir * maxPath;
2152  double dist1Est = distEst;
2153  if (distSol.numberOfSolutions() > 0) {
2154  distEst = distSol.first();
2155  dist1Est = distSol.first();
2156  if (distSol.numberOfSolutions() > 1 &&
2157  (std::abs(distEst) < tol || (propDir * distEst < -tol && propDir * distSol.second() > tol)))
2158  distEst = distSol.second();
2159  }
2160  // select input surfaces;
2161  // do not accept trivial solutions (at the surface)
2162  // but include them into further distance estimate (aca return to the same surface)
2163  if (distEst * propDir > -tol) {
2164  if (distSol.currentDistance() > 500.) {
2165  cache.m_currentDist[iCurr] = std::pair<int, std::pair<double, double>>(
2166  0, std::pair<double, double>(distEst, distSol.currentDistance(true)));
2167  } else {
2168  cache.m_currentDist[iCurr] = std::pair<int, std::pair<double, double>>(
2169  1, std::pair<double, double>(distEst, distSol.currentDistance(true)));
2170  }
2171  if (tol < propDir * distEst && propDir * distEst < propDir * distanceToTarget) {
2172  distanceToTarget = distEst;
2173  nextSf = iCurr;
2174  }
2175  ++numSf;
2176  } else {
2177  // save the nearest distance to surface
2178  cache.m_currentDist[iCurr] = std::pair<int, std::pair<double, double>>(
2179  -1, std::pair<double, double>(distSol.currentDistance(), distSol.currentDistance(true)));
2180  }
2181  if (std::abs(dist1Est) < tol)
2182  startSf = (int)iCurr;
2183  ++iCurr;
2184  }
2185 
2186  if (distanceToTarget == maxPath || numSf == 0)
2187  return false;
2188 
2189  // these do not change
2190  std::vector<std::pair<int, std::pair<double, double>>>::iterator vsBeg = cache.m_currentDist.begin();
2191  std::vector<std::pair<int, std::pair<double, double>>>::iterator vsEnd = cache.m_currentDist.end();
2192  const int num_vs_dist = cache.m_currentDist.size();
2193 
2194  // Set initial step length to 100mm or the distance to the target surface.
2195  double h = 0;
2196  double absPath = 0.;
2197  distanceToTarget > 100. ? h = 100. : distanceToTarget < -100. ? h = -100. : h = distanceToTarget;
2198 
2199  const Trk::IdentifiedMaterial* binIDMat = nullptr;
2200  // Adapt step size to the material binning : change of bin layer triggers dump of material effects
2201  if (cache.m_binMat) {
2202  const Trk::BinUtility* lbu = cache.m_binMat->layerBinUtility(position);
2203  if (lbu) {
2204  cache.m_currentLayerBin = cache.m_binMat->layerBin(position);
2205  binIDMat = cache.m_binMat->material(position);
2206  std::pair<size_t, float> dist2next = lbu->distanceToNext(position, propDir * direction0);
2207  if (dist2next.first < lbu->bins() && std::abs(dist2next.second) > 1. &&
2208  std::abs(dist2next.second) < std::abs(h)) {
2209  h = dist2next.second * propDir;
2210  }
2211  if (binIDMat)
2212  cache.m_material = binIDMat->first.get();
2213  }
2214  }
2215 
2216  // Step to within distanceTolerance, then do the rest as a Taylor expansion.
2217  // Keep distanceTolerance within [1 nanometer, 10 microns].
2218  // This means that no final Taylor expansions beyond 10 microns and no
2219  // Runge-Kutta steps less than 1 nanometer are allowed.
2220  double distanceTolerance = std::min(std::max(std::abs(distanceToTarget) * cache.m_tolerance, 1e-6), 1e-2);
2221 
2222  // bremstrahlung : sample if activated
2223  if (cache.m_brem) {
2224  mom = std::abs(1. / P[6]);
2225  sampleBrem(cache, mom);
2226  }
2227 
2228  while (numSf > 0 && (std::abs(distanceToTarget) > distanceTolerance ||
2229  std::abs(path + distanceStepped) < tol)) { // Step until within tolerance
2230  // Do the step. Stop the propagation if the energy goes below m_momentumCutOff
2231  if (!rungeKuttaStep(cache, errorPropagation, h, P, dDir, BG1, firstStep, distanceStepped)) {
2232  // emit brem photon before stopped ?
2233  if (cache.m_brem) {
2234  if (m_momentumCutOff < cache.m_bremEmitThreshold && m_simMatUpdator) {
2235  Amg::Vector3D position(P[0], P[1], P[2]);
2236  Amg::Vector3D direction(P[3], P[4], P[5]);
2237  m_simMatUpdator->recordBremPhoton(cache.m_timeIn + cache.m_timeOfFlight + cache.m_timeStep,
2238  m_momentumCutOff, cache.m_bremMom, position, direction,
2239  cache.m_particle);
2240  // the recoil can be skipped here
2241  for (int i = 0; i < 3; ++i)
2242  P[3 + i] = direction[i];
2243  // end recoil ( momentum not adjusted here ! continuous energy loss maintained for the moment)
2244  }
2245  }
2246  // collect material and update timing
2247  path = path + distanceStepped;
2248  // timing
2249  mom = std::abs(1. / P[6]);
2250  beta = mom / std::sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass);
2251  cache.m_timeStep += distanceStepped / beta / Gaudi::Units::c_light;
2252 
2253  if (std::abs(distanceStepped) > 0.001) {
2254  cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL * log(std::abs(distanceStepped));
2255  }
2256  // update straggling covariance
2257  if (errorPropagation && cache.m_straggling) {
2258  // 15% of the Radition moves the MOP value thus only 85% is accounted for by the Mean-MOP shift
2259  double sigTot2 = cache.m_sigmaIoni * cache.m_sigmaIoni + cache.m_sigmaRad * cache.m_sigmaRad;
2260  // /(beta*beta*p*p*p*p) transforms Var(E) to Var(q/p)
2261  double bp2 = beta * mom * mom;
2262  cache.m_stragglingVariance += sigTot2 / (bp2 * bp2) * distanceStepped * distanceStepped;
2263  }
2264  if (cache.m_matstates || errorPropagation) {
2265  cache.m_combinedEloss.update(
2266  cache.m_delIoni * distanceStepped, cache.m_sigmaIoni * std::abs(distanceStepped),
2267  cache.m_delRad * distanceStepped, cache.m_sigmaRad * std::abs(distanceStepped), cache.m_MPV);
2268  }
2269  if (cache.m_material && cache.m_material->x0() != 0.) {
2270  cache.m_combinedThickness += propDir * distanceStepped / cache.m_material->x0();
2271  }
2272 
2273  return false;
2274  }
2275  path = path + distanceStepped;
2276  absPath += std::abs(distanceStepped);
2277 
2278  // timing
2279  mom = std::abs(1. / P[6]);
2280  beta = mom / std::sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass);
2281  cache.m_timeStep += distanceStepped / beta / Gaudi::Units::c_light;
2282 
2283  if (std::abs(distanceStepped) > 0.001)
2284  cache.m_sigmaIoni = cache.m_sigmaIoni - cache.m_kazL * log(std::abs(distanceStepped));
2285  // update straggling covariance
2286  if (errorPropagation && cache.m_straggling) {
2287  // 15% of the Radition moves the MOP value thus only 85% is accounted for by the Mean-MOP shift
2288  double sigTot2 = cache.m_sigmaIoni * cache.m_sigmaIoni + cache.m_sigmaRad * cache.m_sigmaRad;
2289  // /(beta*beta*p*p*p*p) transforms Var(E) to Var(q/p)
2290  double bp2 = beta * mom * mom;
2291  cache.m_stragglingVariance += sigTot2 / (bp2 * bp2) * distanceStepped * distanceStepped;
2292  }
2293  if (cache.m_matstates || errorPropagation) {
2294  cache.m_combinedEloss.update(
2295  cache.m_delIoni * distanceStepped, cache.m_sigmaIoni * std::abs(distanceStepped),
2296  cache.m_delRad * distanceStepped, cache.m_sigmaRad * std::abs(distanceStepped), cache.m_MPV);
2297  }
2298  if (cache.m_material && cache.m_material->x0() != 0.) {
2299  cache.m_combinedThickness += propDir * distanceStepped / cache.m_material->x0();
2300  }
2301 
2302  if (absPath > maxPath)
2303  return false;
2304 
2305  // path limit implemented
2306  if (cache.m_propagateWithPathLimit > 0 && cache.m_pathLimit <= path) {
2307  ++cache.m_propagateWithPathLimit;
2308  return true;
2309  }
2310 
2311  bool restart = false;
2312  // in case of problems, make shorter steps
2313  if (propDir * path < -tol || absPath - std::abs(path) > 10.) {
2314  helpSoft = std::abs(path) / absPath > 0.5 ? std::abs(path) / absPath : 0.5;
2315  }
2316 
2317  Amg::Vector3D position(P[0], P[1], P[2]);
2318  Amg::Vector3D direction(P[3], P[4], P[5]);
2319 
2320  // Adapt step size to the material binning : change of bin layer triggers dump of material effects
2321  float distanceToNextBin = h; // default
2322  if (cache.m_binMat) {
2323  const Trk::BinUtility* lbu = cache.m_binMat->layerBinUtility(position);
2324  if (lbu) {
2325  size_t layerBin = cache.m_binMat->layerBin(position);
2326  const Trk::IdentifiedMaterial* iMat = cache.m_binMat->material(position);
2327  std::pair<size_t, float> dist2next = lbu->distanceToNext(position, propDir * direction);
2328  distanceToNextBin = dist2next.second;
2329  if (layerBin != cache.m_currentLayerBin) { // step into new bin
2330  // check the overshoot
2331  std::pair<size_t, float> dist2previous = lbu->distanceToNext(position, -propDir * direction);
2332  float stepOver = dist2previous.second;
2333  double localp[5];
2335  auto cPar = std::make_unique<Trk::CurvilinearParameters>(Amg::Vector3D(P[0], P[1], P[2]), localp[2],
2336  localp[3], localp[4]);
2337  if (cache.m_identifiedParameters) {
2338  if (binIDMat && binIDMat->second > 0 && !iMat) { // exit from active layer
2339  cache.m_identifiedParameters->emplace_back(cPar->clone(), -binIDMat->second);
2340  } else if (binIDMat && binIDMat->second > 0 &&
2341  (iMat->second == 0 || iMat->second == binIDMat->second)) { // exit from active layer
2342  cache.m_identifiedParameters->emplace_back(cPar->clone(), -binIDMat->second);
2343  } else if (iMat && iMat->second > 0) { // entry active layer
2344  cache.m_identifiedParameters->emplace_back(cPar->clone(), iMat->second);
2345  }
2346  }
2347  if (cache.m_hitVector) {
2348  double hitTiming = cache.m_timeIn + cache.m_timeOfFlight + cache.m_timeStep;
2349  if (binIDMat && binIDMat->second > 0 && !iMat) { // exit from active layer
2350  cache.m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2351  } else if (binIDMat && binIDMat->second > 0 &&
2352  (iMat->second == 0 || iMat->second == binIDMat->second)) { // exit from active layer
2353  cache.m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2354  } else if (iMat && iMat->second > 0) { // entry active layer
2355  cache.m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, iMat->second, 0.);
2356  }
2357  }
2358 
2359  cache.m_currentLayerBin = layerBin;
2360  binIDMat = iMat;
2361  if (binIDMat) {
2362  // change of material triggers update of the cache
2363  // @TODO Coverity complains about a possible NULL pointer dereferencing here
2364  // because the code above does not explicitly forbid m_material to be NULL and m_material is used
2365  // unchecked inside updateMaterialEffects.
2366  // Can m_material be NULL at this point ?
2367  if (cache.m_material) {
2368  updateMaterialEffects(cache, mom, sin(direction.theta()), sumPath + path - stepOver);
2369  }
2370  cache.m_material = binIDMat->first.get();
2371  }
2372  // recalculate distance to next bin
2373  if (distanceToNextBin < h) {
2374  Amg::Vector3D probe = position + (distanceToNextBin + h) * propDir * direction;
2375  std::pair<size_t, float> d2n = lbu->distanceToNext(probe, propDir * direction);
2376  distanceToNextBin += d2n.second + h;
2377  }
2378  } else if (dist2next.first < lbu->bins() && std::abs(distanceToNextBin) < 0.01 &&
2379  h > 0.01) { // tolerance 10 microns ?
2380  double localp[5];
2382  auto cPar = std::make_unique<Trk::CurvilinearParameters>(Amg::Vector3D(P[0], P[1], P[2]), localp[2],
2383  localp[3], localp[4]);
2384 
2385  const Trk::IdentifiedMaterial* nextMat = binIDMat;
2386  // need to know what comes next
2387  Amg::Vector3D probe = position + (distanceToNextBin + 0.01) * propDir * direction.normalized();
2388  nextMat = cache.m_binMat->material(probe);
2389 
2390  if (cache.m_identifiedParameters) {
2391  if (binIDMat && binIDMat->second > 0 && !nextMat) { // exit from active layer
2392  cache.m_identifiedParameters->emplace_back(cPar->clone(), -binIDMat->second);
2393  } else if (binIDMat && binIDMat->second > 0 &&
2394  (nextMat->second == 0 || nextMat->second == binIDMat->second)) {
2395  // exit from active layer
2396  cache.m_identifiedParameters->emplace_back(cPar->clone(), -binIDMat->second);
2397  } else if (nextMat && nextMat->second > 0) { // entry active layer
2398  cache.m_identifiedParameters->emplace_back(cPar->clone(), nextMat->second);
2399  }
2400  }
2401  if (cache.m_hitVector) {
2402  double hitTiming = cache.m_timeIn + cache.m_timeOfFlight + cache.m_timeStep;
2403  if (binIDMat && binIDMat->second > 0 && !nextMat) { // exit from active layer
2404  cache.m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2405  } else if (binIDMat && binIDMat->second > 0 &&
2406  (nextMat->second == 0 ||
2407  nextMat->second == binIDMat->second)) { // exit from active layer
2408  cache.m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, -binIDMat->second, 0.);
2409  } else if (nextMat && nextMat->second > 0) { // entry active layer
2410  cache.m_hitVector->emplace_back(cPar->uniqueClone(), hitTiming, nextMat->second, 0.);
2411  }
2412  }
2413 
2414  cache.m_currentLayerBin = dist2next.first;
2415  if (binIDMat != nextMat) { // change of material triggers update of the cache
2416  binIDMat = nextMat;
2417  if (binIDMat) {
2418  assert(cache.m_material);
2419  updateMaterialEffects(cache, mom, sin(direction.theta()), sumPath + path);
2420  cache.m_material = binIDMat->first.get();
2421  }
2422  }
2423  // recalculate distance to next bin
2424  std::pair<size_t, float> d2n = lbu->distanceToNext(probe, propDir * direction.normalized());
2425  distanceToNextBin += d2n.second + 0.01;
2426  }
2427  // TODO: trigger the update of material properties and recalculation of distance to the target sliding
2428  // surface
2429  }
2430  }
2431 
2432  // Calculate new distance to targets
2433  bool flipDirection = false;
2434  numSf = 0;
2435  nextSfCand = nextSf;
2436  double dev = direction0.dot(direction);
2437  std::vector<DestSurf>::iterator sIter = sBeg;
2438  std::vector<std::pair<int, std::pair<double, double>>>::iterator vsIter = vsBeg;
2439  int ic = 0;
2440  int numRestart = 0;
2441 
2442  if (cache.m_brem) {
2443  if (mom < cache.m_bremEmitThreshold && m_simMatUpdator) {
2444  // ST : strictly speaking, the emission point should be shifted backwards a bit
2445  // (mom-m_bremEmitThreshold) this seems to be a minor point
2446  m_simMatUpdator->recordBremPhoton(cache.m_timeIn + cache.m_timeOfFlight + cache.m_timeStep, mom,
2447  cache.m_bremMom, position, direction, cache.m_particle);
2448  cache.m_bremEmitThreshold = 0.;
2449  }
2450  if (mom < cache.m_bremSampleThreshold)
2451  sampleBrem(cache, cache.m_bremSampleThreshold);
2452  }
2453 
2454  for (; vsIter != vsEnd; ++vsIter) {
2455  if (restart) {
2456  ++numRestart;
2457  if (numRestart > restartLimit)
2458  return false;
2459 
2460  vsIter = vsBeg;
2461  ic = 0;
2462  sIter = sBeg;
2463  distanceToTarget = propDir * maxPath;
2464  nextSf = -1;
2465  nextSfCand = -1;
2466  restart = false;
2467  helpSoft = 1.;
2468  }
2469  if ((*vsIter).first != -1 &&
2470  (ic == nextSf || (*vsIter).first == 1 || nextSf < 0 || std::abs((*vsIter).second.first) < 500. ||
2471  std::abs(path) > 0.5 * std::abs((*vsIter).second.second))) {
2472  previousDistance = (*vsIter).second.first;
2473  Trk::DistanceSolution distSol =
2474  (*sIter).first->straightLineDistanceEstimate(position, propDir * direction);
2475  double distanceEst = -propDir * maxPath;
2476  if (distSol.numberOfSolutions() > 0) {
2477  distanceEst = distSol.first();
2478  if (distSol.numberOfSolutions() > 1 &&
2479  std::abs(distSol.first() * propDir + distanceStepped - previousDistance) >
2480  std::abs(distSol.second() * propDir + distanceStepped - previousDistance)) {
2481  distanceEst = distSol.second();
2482  }
2483  // Peter Kluit: avoid jumping into other (distSol.first->second) distance solution for start surface
2484  // with negative distance solution
2485  // negative distanceEst will trigger flipDirection = true and will iterate to the start
2486  // surface this will lead to very similar positions for multiple propagator calls and
2487  // many tiny X0 scatterers
2488  if (ic == startSf && distanceEst < 0 && distSol.first() > 0)
2489  distanceEst = distSol.first();
2490  }
2491  // eliminate close surface if path too small
2492  if (ic == nextSf && std::abs(distanceEst) < tol && std::abs(path) < tol) {
2493  (*vsIter).first = -1;
2494  vsIter = vsBeg;
2495  restart = true;
2496  distanceToTarget = maxPath;
2497  nextSf = -1;
2498  continue;
2499  }
2500 
2501  // If h and distance are in opposite directions, target is passed. Flip propagation direction
2502  // Verify if true intersection
2503  // if ( h * propDir * distanceEst < 0. && std::abs(distanceEst)>distanceTolerance ) {
2504  if ((*vsIter).second.first * propDir * distanceEst < 0. &&
2505  std::abs(distanceEst) > distanceTolerance) {
2506  // verify change of sign in signedDistance ( after eliminating situations where this is meaningless
2507  // )
2508  if (!distSol.signedDistance() || std::abs(distSol.currentDistance(true)) < tol ||
2509  std::abs((*vsIter).second.second) < tol ||
2510  (*vsIter).second.second * distSol.currentDistance(true) < 0) { // true intersection
2511  if (ic == nextSf) {
2512  ((*vsIter).first)++;
2513  // eliminate surface if diverging
2514  if ((*vsIter).first > 3) {
2515  helpSoft = fmax(0.05, 1. - 0.05 * (*vsIter).first);
2516  if ((*vsIter).first > 20)
2517  helpSoft = 1. / (*vsIter).first;
2518  }
2519  // take care of eliminating when number of flips even - otherwise it may end up at the start !
2520  if ((*vsIter).first > 50 && h * propDir > 0) {
2521  // std::abs(distanceEst) >= std::abs(previousDistance) ) {
2522  (*vsIter).first = -1;
2523  vsIter = vsBeg;
2524  restart = true;
2525  continue;
2526  }
2527  if ((*vsIter).first != -1)
2528  flipDirection = true;
2529  } else if (std::abs((*vsIter).second.second) > tol &&
2530  std::abs(distSol.currentDistance(true)) > tol) {
2531  // here we need to compare with distance from current closest
2532  if (ic > nextSf && nextSf != -1) { // easy case, already calculated
2533  if (propDir * distanceEst < (cache.m_currentDist.at(nextSf)).second.first - tol) {
2534  if ((*vsIter).first != -1) {
2535  ((*vsIter).first)++;
2536  flipDirection = true;
2537  nextSf = ic;
2538  }
2539  }
2540  } else if (distanceToTarget >
2541  0.) { // set as nearest (if not iterating already), will be rewritten later
2542  if ((*vsIter).first != -1) {
2543  ((*vsIter).first)++;
2544  flipDirection = true;
2545  nextSf = ic;
2546  }
2547  }
2548  }
2549  } else if (ic == nextSf) {
2550  vsIter = vsBeg;
2551  restart = true;
2552  continue;
2553  }
2554  }
2555 
2556  // save current distance to surface
2557  (*vsIter).second.first = propDir * distanceEst;
2558  (*vsIter).second.second = distSol.currentDistance(true);
2559 
2560  // find closest surface: the step may have been beyond several surfaces
2561  // from all surfaces with 'negative' distance, consider only the one currently designed as 'closest'
2562  // mw if ((*vsIter).first!=-1 && ( distanceEst>-tol || ic==nextSf ) ) {
2563  if ((*vsIter).first != -1 && (distanceEst > 0. || ic == nextSf)) {
2564  ++numSf;
2565  if (distanceEst < std::abs(distanceToTarget)) {
2566  distanceToTarget = propDir * distanceEst;
2567  nextSfCand = ic;
2568  }
2569  }
2570  } else if (std::abs(path) > std::abs((*vsIter).second.second) || dev < 0.985 ||
2571  nextSf < 0) { // keep an eye on surfaces with negative distance; tracks are curved !
2572  Trk::DistanceSolution distSol =
2573  (*sIter).first->straightLineDistanceEstimate(position, propDir * direction);
2574  double distanceEst = -propDir * maxPath;
2575  if (distSol.numberOfSolutions() > 0) {
2576  distanceEst = distSol.first();
2577  }
2578  // save current distance to surface
2579  (*vsIter).second.first = propDir * distanceEst;
2580  (*vsIter).second.second = distSol.currentDistance(true);
2581  // reactivate surface
2582  if (distanceEst > tol && distanceEst < maxPath) {
2583  (*vsIter).first = 0;
2584  } else {
2585  (*vsIter).second.first = distSol.currentDistance() + std::abs(path);
2586  }
2587  if ((*vsIter).first != -1 && distanceEst > 0.) {
2588  ++numSf;
2589  if (distanceEst < std::abs(distanceToTarget)) {
2590  distanceToTarget = propDir * distanceEst;
2591  nextSfCand = ic;
2592  }
2593  }
2594  }
2595  // additional protection - return to the same surface
2596  // eliminate the surface and restart the search
2597  // 04/10/10 ST:infinite loop due to distanceTolerance>tol fixed;
2598  if (std::abs(distanceToTarget) <= distanceTolerance && path * propDir < distanceTolerance) {
2599  (*vsIter).first = -1;
2600  vsIter = vsBeg;
2601  restart = true;
2602  continue;
2603  }
2604  ++sIter;
2605  ++ic;
2606  }
2607  // if next closest not found, propagation failed
2608  if (nextSf < 0 && nextSfCand < 0)
2609  return false;
2610  // flip direction
2611  if (flipDirection) {
2612  // Out of bounds protection
2613  if (nextSf < 0 || nextSf >= num_vs_dist)
2614  return false;
2615  distanceToTarget = (*(vsBeg + nextSf)).second.first;
2616  h = -h;
2617  } else if (nextSfCand != nextSf) {
2618  nextSf = nextSfCand;
2619  // Out of bounds protection
2620  if (nextSf < 0 || nextSf >= num_vs_dist)
2621  return false;
2622  if (cache.m_currentDist[nextSf].first < 3)
2623  helpSoft = 1.;
2624  }
2625 
2626  // don't step beyond surfaces - adjust step
2627  if (std::abs(h) > std::abs(distanceToTarget))
2628  h = distanceToTarget;
2629 
2630  // don't step beyond bin boundary - adjust step
2631  if (cache.m_binMat && std::abs(h) > std::abs(distanceToNextBin) + 0.001) {
2632  if (distanceToNextBin > 0) { // TODO : investigate source of negative distance in BinningData
2633  h = distanceToNextBin * propDir;
2634  }
2635  }
2636 
2637  if (helpSoft < 1.)
2638  h *= helpSoft;
2639 
2640  // don't step much beyond path limit
2641  if (cache.m_propagateWithPathLimit > 0 && h > cache.m_pathLimit)
2642  h = cache.m_pathLimit + tol;
2643 
2644  // Abort if maxPath is reached
2645  if (std::abs(path) > maxPath)
2646  return false;
2647 
2648  if (steps++ > cache.m_maxSteps)
2649  return false; // Too many steps, something is wrong
2650  }
2651 
2652  if (!numSf)
2653  return false;
2654 
2655  // Use Taylor expansions to step the remaining distance (typically microns).
2656  path = path + distanceToTarget;
2657 
2658  // timing
2659  mom = std::abs(1. / P[6]);
2660  beta = mom / std::sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass);
2661  cache.m_timeStep += distanceToTarget / beta / Gaudi::Units::c_light;
2662 
2663  // pos = pos + h*dir + 1/2*h*h*dDir. Second order Taylor expansion.
2664  pos[0] = pos[0] + distanceToTarget * (dir[0] + 0.5 * distanceToTarget * dDir[0]);
2665  pos[1] = pos[1] + distanceToTarget * (dir[1] + 0.5 * distanceToTarget * dDir[1]);
2666  pos[2] = pos[2] + distanceToTarget * (dir[2] + 0.5 * distanceToTarget * dDir[2]);
2667 
2668  // dir = dir + h*dDir. First order Taylor expansion (Euler).
2669  dir[0] = dir[0] + distanceToTarget * dDir[0];
2670  dir[1] = dir[1] + distanceToTarget * dDir[1];
2671  dir[2] = dir[2] + distanceToTarget * dDir[2];
2672 
2673  // Normalize dir
2674  double norm = 1. / std::sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
2675  dir[0] = norm * dir[0];
2676  dir[1] = norm * dir[1];
2677  dir[2] = norm * dir[2];
2678  P[42] = dDir[0];
2679  P[43] = dDir[1];
2680  P[44] = dDir[2];
2681 
2682  // collect all surfaces with distance below tolerance
2683  std::vector<std::pair<int, std::pair<double, double>>>::iterator vsIter = vsBeg;
2684 
2685  int index = 0;
2686  for (; vsIter != vsEnd; ++vsIter) {
2687  if ((*vsIter).first != -1 && propDir * (*vsIter).second.first >= propDir * distanceToTarget - tol &&
2688  propDir * (*vsIter).second.first < 0.01 && index != nextSf) {
2689  solutions.push_back(index);
2690  }
2691  if (index == nextSf)
2692  solutions.push_back(index);
2693  ++index;
2694  }
2695 
2696  return true;
2697 }
2698 
2700 // dump material effects
2703  double path) const {
2704 
2705  // kinematics
2706  double mom = parms->momentum().mag();
2707 
2708  // first update to make sure all material counted
2709  updateMaterialEffects(cache, mom, sin(parms->momentum().theta()), path);
2710 
2711  if (cache.m_extrapolationCache) {
2715  cache.m_combinedEloss.sigmaRad());
2716  }
2717  // output
2718  if (cache.m_matstates) {
2719  auto eloss = !m_detailedEloss
2720  ? std::make_unique<Trk::EnergyLoss>(cache.m_combinedEloss.deltaE(),
2721  cache.m_combinedEloss.sigmaDeltaE())
2722  : std::make_unique<Trk::EnergyLoss>(
2727 
2728  auto sa = Trk::ScatteringAngles(0., 0., std::sqrt(cache.m_covariance(2, 2)),
2729  std::sqrt(cache.m_covariance(3, 3)));
2730 
2731  auto cvlTP = parms->uniqueClone();
2732  auto mefot = std::make_unique<Trk::MaterialEffectsOnTrack>(cache.m_combinedThickness, sa,
2733  std::move(eloss), cvlTP->associatedSurface());
2734 
2735  cache.m_matstates->push_back(new TrackStateOnSurface(nullptr, std::move(cvlTP), std::move(mefot)));
2736  }
2737 
2738  cache.m_matdump_lastpath = path;
2739 
2740  // clean-up
2741  cache.m_combinedCovariance += cache.m_covariance;
2742  cache.m_covariance.setZero();
2743  cache.m_combinedThickness = 0.;
2744  cache.m_combinedEloss.set(0., 0., 0., 0., 0., 0.);
2745 }
2746 
2747 // Smear momentum ( simulation mode )
2749 void Trk::STEP_Propagator::smear(Cache& cache, double& phi, double& theta, const Trk::TrackParameters* parms,
2750  double radDist) const {
2751  if (cache.m_particle == Trk::geantino)
2752  return;
2753  if (!parms)
2754  return;
2755 
2756  if (!cache.m_randomEngine) {
2757  cache.m_randomEngine = getRandomEngine(cache.m_ctx);
2758  }
2759 
2760  // Calculate polar angle
2761  double particleMass = Trk::ParticleMasses::mass[cache.m_particle]; // Get particle mass from
2762  // ParticleHypothesis
2763  double momentum = parms->momentum().mag();
2764  double energy = std::sqrt(momentum * momentum + particleMass * particleMass);
2765  double beta = momentum / energy;
2766  double th = std::sqrt(2.) * 15. * std::sqrt(radDist) / (beta * momentum) *
2767  CLHEP::RandGauss::shoot(cache.m_randomEngine); // Moliere
2768  // double th = (sqrt(2.)*13.6*std::sqrt(radDist)/(beta*momentum)) *
2769  // (1.+0.038*log(radDist/(beta*beta))) * m_gaussian->shoot(); //Highland
2770 
2771  // Calculate azimuthal angle
2772  double ph = 2. * M_PI * CLHEP::RandFlat::shoot(cache.m_randomEngine);
2773 
2775  Amg::AngleAxis3D(-phi, Amg::Vector3D(0., 0., 1.)));
2776  Amg::Vector3D dir0(0., 0., 1.);
2777  Amg::Vector3D rotated = rot.inverse() * Amg::AngleAxis3D(ph, Amg::Vector3D(0., 0., 1.)) *
2778  Amg::AngleAxis3D(th, Amg::Vector3D(0., 1., 0.)) * dir0;
2779 
2780  theta = rotated.theta();
2781  phi = rotated.phi();
2782 }
2783 
2784 // Sample momentum of brem photon ( simulation mode )
2786 void Trk::STEP_Propagator::sampleBrem(Cache& cache, double mom) const {
2787  if (!cache.m_randomEngine) {
2788  cache.m_randomEngine = getRandomEngine(cache.m_ctx);
2789  }
2790  double rndx = CLHEP::RandFlat::shoot(cache.m_randomEngine);
2791  double rnde = CLHEP::RandFlat::shoot(cache.m_randomEngine);
2792 
2793  // sample visible fraction of the mother momentum taken according to 1/f
2794  double eps = cache.m_momentumCutOff / mom;
2795  cache.m_bremMom = pow(eps, pow(rndx, exp(1.))) * mom; // adjustment here ?
2796  cache.m_bremSampleThreshold = mom - cache.m_bremMom;
2797  cache.m_bremEmitThreshold = mom - rnde * cache.m_bremMom;
2798 }
2799 
2800 void Trk::STEP_Propagator::getFieldCacheObject(Cache& cache, const EventContext& ctx) const {
2801  SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, ctx};
2802  const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
2803  if (fieldCondObj == nullptr) {
2804  ATH_MSG_ERROR("extrapolate: Failed to retrieve AtlasFieldCacheCondObj with key "
2805  << m_fieldCacheCondObjInputKey.key());
2806  return;
2807  }
2808  fieldCondObj->getInitializedCache(cache.m_fieldCache);
2809 }
2810 
2811 CLHEP::HepRandomEngine*
2812 Trk::STEP_Propagator::getRandomEngine (const EventContext& ctx) const
2813 {
2814  if (!m_simulation || !m_rngWrapper) {
2815  return nullptr;
2816  }
2817  if (m_rngWrapper->evtSeeded(ctx) != ctx.evt()) {
2818  // Ok, the wrappers are unique to this algorithm and a given slot,
2819  // so cannot be accessed concurrently.
2820  ATHRNG::RNGWrapper* wrapper_nc ATLAS_THREAD_SAFE = m_rngWrapper;
2821  wrapper_nc->setSeed (this->name(), ctx);
2822  }
2823  return m_rngWrapper->getEngine (ctx);
2824 }
BinnedMaterial.h
Trk::ExtrapolationCache::updateX0
void updateX0(double x0)
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Trk::Material::averageZ
float averageZ() const
Definition: Material.h:228
Trk::DistanceSolution::currentDistance
double currentDistance(bool signedDist=false) const
Current distance to surface (spatial), signed (along/opposite to surface normal) if input argument tr...
Trk::STEP_Propagator::dumpMaterialEffects
void dumpMaterialEffects(Cache &cache, const Trk::CurvilinearParameters *trackParameters, double path) const
Definition: STEP_Propagator.cxx:2702
Trk::STEP_Propagator::Cache::m_kazL
double m_kazL
Definition: STEP_Propagator.h:365
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::EnergyLoss::update
void update(double ioni, double sigi, double rad, double sigr, bool mpv=false)
Trk::DistanceSolution::signedDistance
bool signedDistance() const
This method indicates availability of signed current distance (false for Perigee and StraighLineSurfa...
Trk::STEP_Propagator::Cache::m_bremEmitThreshold
double m_bremEmitThreshold
Definition: STEP_Propagator.h:385
Trk::RungeKuttaUtils::jacobianTransformCurvilinearToLocal
void jacobianTransformCurvilinearToLocal(const Trk::TrackParameters &, double *ATH_RESTRICT)
Trk::STEP_Propagator::Cache::m_momentumCutOff
double m_momentumCutOff
Definition: STEP_Propagator.h:377
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
ScatteringAngles.h
Trk::MaterialInteraction::dEdl_ionization
static double dEdl_ionization(double p, const Material &mat, ParticleHypothesis particle, double &sigma, double &kazL)
dE/dl ionization energy loss per path unit
Definition: MaterialInteraction.cxx:117
Trk::STEP_Propagator::Cache::m_bremSampleThreshold
double m_bremSampleThreshold
Definition: STEP_Propagator.h:386
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:79
Trk::k0
@ k0
Definition: ParticleHypothesis.h:38
inline_hints.h
BenchmarkConfig.useCache
useCache
Definition: BenchmarkConfig.py:63
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:128
Trk::RungeKuttaUtils::transformGlobalToCurvilinear
void transformGlobalToCurvilinear(bool, double *ATH_RESTRICT, double *ATH_RESTRICT, double *ATH_RESTRICT)
Trk::STEP_Propagator::Cache::m_multipleScattering
bool m_multipleScattering
Definition: STEP_Propagator.h:356
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
Trk::MagneticFieldProperties
Definition: MagneticFieldProperties.h:31
Trk::ParametersBase::charge
double charge() const
Returns the charge.
Surface.h
AtlasFieldCacheCondObj
Definition: AtlasFieldCacheCondObj.h:19
Trk::STEP_Propagator::propagateM
virtual std::unique_ptr< Trk::TrackParameters > propagateM(const EventContext &ctx, const Trk::TrackParameters &trackParameters, std::vector< Trk::DestSurf > &targetSurfaces, Trk::PropDirection propagationDirection, const MagneticFieldProperties &magneticFieldProperties, ParticleHypothesis particle, std::vector< unsigned int > &solutions, std::vector< const Trk::TrackStateOnSurface * > *matstates, std::vector< std::pair< std::unique_ptr< Trk::TrackParameters >, int >> *intersections, double &path, bool usePathLimit=false, bool returnCurv=false, const Trk::TrackingVolume *tVol=nullptr, Trk::ExtrapolationCache *=nullptr) const override final
Propagate parameters and covariance with search of closest surface and material collection.
Definition: STEP_Propagator.cxx:1408
Trk::STEP_Propagator::Cache::m_ctx
const EventContext & m_ctx
Definition: STEP_Propagator.h:413
Trk::STEP_Propagator::Cache::m_matPropOK
bool m_matPropOK
Definition: STEP_Propagator.h:351
Trk::DistanceSolution
Definition: DistanceSolution.h:25
Trk::STEP_Propagator::Cache::m_P
double m_P[45]
Definition: STEP_Propagator.h:387
MagField::AtlasFieldCache::getFieldZR
void getFieldZR(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field valaue on the z-r plane at given position works only inside the solenoid.
Definition: AtlasFieldCache.cxx:86
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
make_unique
std::unique_ptr< T > make_unique(Args &&... args)
Definition: SkimmingToolEXOT5.cxx:23
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
Trk::DistanceSolution::numberOfSolutions
int numberOfSolutions() const
Number of intersection solutions.
index
Definition: index.py:1
Trk::PathLimit::updateMat
void updateMat(float dX0, float Z, float dL0)
collected material update
Definition: HelperStructs.h:51
Trk::ParametersBase::associatedSurface
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
Trk::ParametersBase::uniqueClone
std::unique_ptr< ParametersBase< DIM, T > > uniqueClone() const
clone method for polymorphic deep copy returning unique_ptr; it is not overriden, but uses the existi...
Definition: ParametersBase.h:97
Trk::STEP_Propagator::Cache::m_stragglingVariance
double m_stragglingVariance
Definition: STEP_Propagator.h:369
hist_file_dump.d
d
Definition: hist_file_dump.py:142
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
Trk::pi0
@ pi0
Definition: ParticleHypothesis.h:37
DMTest::P
P_v1 P
Definition: P.h:23
Trk::CurvilinearParameters
CurvilinearParametersT< TrackParametersDim, Charged, PlaneSurface > CurvilinearParameters
Definition: Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:29
Trk::STEP_Propagator::Cache::m_charge
double m_charge
Definition: STEP_Propagator.h:374
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
EventPrimitivesHelpers.h
Trk::STEP_Propagator::Cache::m_bremMom
double m_bremMom
Definition: STEP_Propagator.h:384
DMTest::C
C_v1 C
Definition: C.h:26
hitVector
std::vector< FPGATrackSimHit > hitVector
Definition: FPGATrackSimCluster.h:22
Amg::saneCovarianceElement
bool saneCovarianceElement(double ele)
A covariance matrix formally needs to be positive semi definite.
Definition: EventPrimitivesCovarianceHelpers.h:63
BinUtility.h
Trk::BinUtility::distanceToNext
std::pair< size_t, float > distanceToNext(const Amg::Vector3D &position, const Amg::Vector3D &direction, size_t ba=0) const
Distance estimate to next bin
Definition: BinUtility.h:154
Trk::STEP_Propagator::globalPositions
virtual void globalPositions(const EventContext &ctx, std::deque< Amg::Vector3D > &positionsList, const TrackParameters &trackParameters, const MagneticFieldProperties &magneticFieldProperties, const CylinderBounds &cylinderBounds, double maxStepSize, ParticleHypothesis particle, const Trk::TrackingVolume *tVol=0) const override final
Return a list of positions along the track.
Definition: STEP_Propagator.cxx:1774
Trk::Material::zOverAtimesRho
float zOverAtimesRho() const
access to members
Definition: Material.h:226
Trk::EnergyLoss::sigmaDeltaE
double sigmaDeltaE() const
returns the symmatric error
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::STEP_Propagator::Cache::m_solenoid
bool m_solenoid
Switch for turning off material effects temporarily.
Definition: STEP_Propagator.h:349
Trk::STEP_Propagator::Cache::m_hitVector
std::vector< Trk::HitInfo > * m_hitVector
Definition: STEP_Propagator.h:397
Trk::ScatteringAngles
represents a deflection of the track caused through multiple scattering in material.
Definition: ScatteringAngles.h:26
Trk::CylinderSurface::bounds
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder)
Amg::saneVector
bool saneVector(const AmgVector(N) &vec)
Definition: EventPrimitivesHelpers.h:28
Trk::SurfaceType
SurfaceType
Definition: SurfaceTypes.h:17
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
STEP_Propagator.h
Trk::TrackSurfaceIntersection
Definition: TrackSurfaceIntersection.h:32
Trk::alongMomentum
@ alongMomentum
Definition: PropDirection.h:20
Trk::EnergyLoss::sigmaRad
double sigmaRad() const
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
Trk::DistanceSolution::first
double first() const
Distance to first intersection solution along direction.
TransportJacobian.h
Trk::STEP_Propagator::m_simulation
BooleanProperty m_simulation
Definition: STEP_Propagator.h:518
InDetAccessor::qOverP
@ qOverP
perigee
Definition: InDetAccessor.h:35
const
bool const RAWDATA *ch2 const
Definition: LArRodBlockPhysicsV0.cxx:560
Trk::EnergyLoss::meanIoni
double meanIoni() const
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
Trk::STEP_Propagator::Cache::m_sigmaIoni
double m_sigmaIoni
Definition: STEP_Propagator.h:364
Trk::RungeKuttaUtils::transformLocalToGlobal
bool transformLocalToGlobal(bool, const Trk::TrackParameters &, double *ATH_RESTRICT)
python.copyTCTOutput.dDir
dDir
Definition: copyTCTOutput.py:60
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition: Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
Trk::STEP_Propagator::Cache::m_tolerance
double m_tolerance
Definition: STEP_Propagator.h:376
Trk::AlignableTrackingVolume::binnedMaterial
const BinnedMaterial * binnedMaterial() const
access to binned material
Definition: AlignableTrackingVolume.h:68
ParamDefs.h
Trk::STEP_Propagator::Cache::m_timeStep
double m_timeStep
Definition: STEP_Propagator.h:372
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
Trk::STEP_Propagator::Cache::m_combinedThickness
double m_combinedThickness
Definition: STEP_Propagator.h:375
Trk::RungeKuttaUtils::transformGlobalToLocal
void transformGlobalToLocal(double *ATH_RESTRICT, double *ATH_RESTRICT)
Trk::BinUtility::bins
size_t bins(size_t ba=0) const
Number of bins.
Definition: BinUtility.h:221
Trk::STEP_Propagator::Cache::m_timeIn
double m_timeIn
Definition: STEP_Propagator.h:383
Trk::STEP_Propagator::Cache::m_maxSteps
int m_maxSteps
Definition: STEP_Propagator.h:380
xAOD::TrackParameters
TrackParameters_v1 TrackParameters
Definition: Event/xAOD/xAODTracking/xAODTracking/TrackParameters.h:11
ATH_RESTRICT
#define ATH_RESTRICT
Definition: restrict.h:31
ExtrapolationCache.h
Trk::STEP_Propagator::m_multipleScattering
BooleanProperty m_multipleScattering
Definition: STEP_Propagator.h:496
Trk::STEP_Propagator::m_randomEngineName
StringProperty m_randomEngineName
Definition: STEP_Propagator.h:527
MagneticFieldProperties.h
MaterialEffectsOnTrack.h
Trk::FastField
@ FastField
call the fast field access method of the FieldSvc
Definition: MagneticFieldMode.h:20
Trk::STEP_Propagator::m_energyLoss
BooleanProperty m_energyLoss
Definition: STEP_Propagator.h:498
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:28
Trk::STEP_Propagator::Cache::m_straggling
bool m_straggling
Definition: STEP_Propagator.h:348
Trk::PropDirection
PropDirection
Definition: PropDirection.h:19
Trk::STEP_Propagator
Definition: STEP_Propagator.h:166
Trk::TimeLimit::tMax
float tMax
Definition: HelperStructs.h:59
Trk::STEP_Propagator::Cache::m_particleMass
double m_particleMass
cache
Definition: STEP_Propagator.h:373
H
#define H(x, y, z)
Definition: MD5.cxx:114
Trk::Surface::isOnSurface
virtual bool isOnSurface(const Amg::Vector3D &glopo, const BoundaryCheck &bchk=true, double tol1=0., double tol2=0.) const
This method returns true if the GlobalPosition is on the Surface for both, within or without check of...
Definition: Surface.cxx:123
Trk::Material::x0
float x0() const
Definition: Material.h:227
Trk::STEP_Propagator::sampleBrem
void sampleBrem(Cache &cache, double mom) const
Definition: STEP_Propagator.cxx:2786
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
beamspotman.steps
int steps
Definition: beamspotman.py:503
Trk::STEP_Propagator::Cache::m_matdump_lastpath
double m_matdump_lastpath
Definition: STEP_Propagator.h:361
python.TriggerHandler.th
th
Definition: TriggerHandler.py:295
Trk::PathLimit::x0Collected
float x0Collected
Definition: HelperStructs.h:36
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
EventPrimitivesToStringConverter.h
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
TrigVtx::gamma
@ gamma
Definition: TrigParticleTable.h:27
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::TrackSurfaceIntersection::position
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
Definition: TrackSurfaceIntersection.h:80
Trk::Surface::createUniqueTrackParameters
virtual ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - charged.
Trk::STEP_Propagator::propagateWithJacobian
bool propagateWithJacobian(Cache &cache, bool errorPropagation, std::vector< DestSurf > &sfs, double *P, Trk::PropDirection propDir, std::vector< unsigned int > &solutions, double &path, double sumPath) const
Definition: STEP_Propagator.cxx:2094
ATH_FLATTEN
#define ATH_FLATTEN
Definition: inline_hints.h:52
beamspotman.n
n
Definition: beamspotman.py:729
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::STEP_Propagator::Cache
stuct to pass information to the heavy lifting calculation internal methods
Definition: STEP_Propagator.h:345
Trk::electron
@ electron
Definition: ParticleHypothesis.h:30
extractSporadic.h
list h
Definition: extractSporadic.py:96
Trk::CylinderSurface
Definition: CylinderSurface.h:55
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
Trk::STEP_Propagator::propagateParameters
virtual std::unique_ptr< Trk::TrackParameters > propagateParameters(const EventContext &ctx, const Trk::TrackParameters &trackParameters, const Trk::Surface &targetSurface, Trk::PropDirection propagationDirection, const Trk::BoundaryCheck &boundaryCheck, const MagneticFieldProperties &magneticFieldProperties, ParticleHypothesis particle, bool returnCurv=false, const Trk::TrackingVolume *tVol=nullptr) const override final
Propagate parameters only.
Definition: STEP_Propagator.cxx:1522
Trk::CylinderBounds
Definition: CylinderBounds.h:46
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
Trk::DistanceSolution::second
double second() const
Distance to second intersection solution along direction (for a cylinder surface)
Trk::EnergyLoss::set
void set(double eLoss, double sigde, double ioni, double sigi, double rad, double sigr)
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Trk::STEP_Propagator::Cache::m_inputThetaVariance
double m_inputThetaVariance
Definition: STEP_Propagator.h:368
Trk::geantino
@ geantino
Definition: ParticleHypothesis.h:29
Trk::EnergyLoss::deltaE
double deltaE() const
returns the
Trk::STEP_Propagator::getFieldCacheObject
void getFieldCacheObject(Cache &cache, const EventContext &ctx) const
Definition: STEP_Propagator.cxx:2800
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::STEP_Propagator::Cache::m_pathLimit
double m_pathLimit
Definition: STEP_Propagator.h:370
Trk::STEP_Propagator::propagateT
virtual std::unique_ptr< Trk::TrackParameters > propagateT(const EventContext &ctx, const Trk::TrackParameters &trackParameters, std::vector< Trk::DestSurf > &targetSurfaces, Trk::PropDirection propagationDirection, const MagneticFieldProperties &magneticFieldProperties, ParticleHypothesis particle, std::vector< unsigned int > &solutions, Trk::PathLimit &path, Trk::TimeLimit &time, bool returnCurv, const Trk::TrackingVolume *tVol, std::vector< Trk::HitInfo > *&hitVector) const override final
Propagate parameters and covariance with search of closest surface time included.
Definition: STEP_Propagator.cxx:1322
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition: AnalysisMisc.h:49
Trk::STEP_Propagator::getRandomEngine
CLHEP::HepRandomEngine * getRandomEngine(const EventContext &ctx) const
Definition: STEP_Propagator.cxx:2812
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::CurvilinearParametersT
Definition: CurvilinearParametersT.h:48
Trk::STEP_Propagator::Cache::m_binMat
const Trk::BinnedMaterial * m_binMat
cache of TrackStateOnSurfaces
Definition: STEP_Propagator.h:389
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Trk::SurfaceType::Cone
@ Cone
Trk::STEP_Propagator::m_simMatUpdator
ToolHandle< ITimedMatEffUpdator > m_simMatUpdator
secondary interactions (brem photon emission)
Definition: STEP_Propagator.h:521
Trk::STEP_Propagator::Cache::m_particle
ParticleHypothesis m_particle
Definition: STEP_Propagator.h:399
Trk::TimeLimit
Definition: HelperStructs.h:58
Trk::muon
@ muon
Definition: ParticleHypothesis.h:31
Trk::neutron
@ neutron
Definition: ParticleHypothesis.h:36
ParticleHypothesis.h
Trk::STEP_Propagator::m_rngWrapper
ATHRNG::RNGWrapper * m_rngWrapper
Random engine.
Definition: STEP_Propagator.h:526
Trk::EnergyLoss::meanRad
double meanRad() const
Trk::STEP_Propagator::Cache::m_delIoni
double m_delIoni
Definition: STEP_Propagator.h:363
RungeKuttaUtils
Definition: RungeKuttaUtils.h:30
Trk::STEP_Propagator::intersect
virtual std::optional< TrackSurfaceIntersection > intersect(const EventContext &ctx, const Trk::TrackParameters &trackParameters, const Trk::Surface &targetSurface, const Trk::MagneticFieldProperties &magneticFieldProperties, ParticleHypothesis particle, const Trk::TrackingVolume *tVol=nullptr) const override final
Propagate parameters and return path (Similar to propagateParameters.
Definition: STEP_Propagator.cxx:1606
Trk::STEP_Propagator::intersectSurface
virtual std::optional< TrackSurfaceIntersection > intersectSurface(const EventContext &ctx, const Surface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP, const MagneticFieldProperties &mft, ParticleHypothesis particle) const override final
Intersection and propagation:
Definition: STEP_Propagator.cxx:1746
beamspotman.dir
string dir
Definition: beamspotman.py:621
Trk::BinUtility
Definition: BinUtility.h:39
grepfile.ic
int ic
Definition: grepfile.py:33
Trk::STEP_Propagator::Cache::m_currentLayerBin
size_t m_currentLayerBin
Definition: STEP_Propagator.h:358
Trk::MaterialInteraction::dEdl_radiation
static double dEdl_radiation(double p, const Material &mat, ParticleHypothesis particle, double &sigma)
dE/dl radiation energy loss per path unit
Definition: MaterialInteraction.cxx:229
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:56
Trk::NoField
@ NoField
Field is set to 0., 0., 0.,.
Definition: MagneticFieldMode.h:18
Trk::STEP_Propagator::m_fieldCacheCondObjInputKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCacheCondObjInputKey
Definition: STEP_Propagator.h:531
EventPrimitives.h
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::SurfaceType::Perigee
@ Perigee
Trk::STEP_Propagator::Cache::m_currentDist
std::vector< std::pair< int, std::pair< double, double > > > m_currentDist
Definition: STEP_Propagator.h:410
Trk::STEP_Propagator::propagateRungeKutta
std::unique_ptr< Trk::TrackParameters > propagateRungeKutta(Cache &cache, bool errorPropagation, const Trk::TrackParameters &trackParameters, std::vector< DestSurf > &targetSurfaces, Trk::PropDirection propagationDirection, const MagneticFieldProperties &magneticFieldProperties, ParticleHypothesis particle, std::vector< unsigned int > &solutions, double &path, bool returnCurv=false) const
Definition: STEP_Propagator.cxx:1884
Trk::TrackStateOnSurface
represents the track state (measurement, material, fit parameters and quality) at a surface.
Definition: TrackStateOnSurface.h:71
Trk::STEP_Propagator::Cache::m_MPV
bool m_MPV
Definition: STEP_Propagator.h:355
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
ATHRNG::RNGWrapper
A wrapper class for event-slot-local random engines.
Definition: RNGWrapper.h:56
VP1PartSpect::E
@ E
Definition: VP1PartSpectFlags.h:21
Trk::STEP_Propagator::Cache::m_detailedElossFlag
bool m_detailedElossFlag
Definition: STEP_Propagator.h:347
xAOD::Curvilinear
@ Curvilinear
Definition: TrackingPrimitives.h:559
charge
double charge(const T &p)
Definition: AtlasPID.h:986
Trk::IdentifiedMaterial
std::pair< std::shared_ptr< Material >, int > IdentifiedMaterial
Definition: BinnedMaterial.h:23
Trk::STEP_Propagator::Cache::m_combinedEloss
Trk::EnergyLoss m_combinedEloss
Definition: STEP_Propagator.h:409
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
Trk::PathLimit
Definition: HelperStructs.h:34
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:73
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::ExtrapolationCache::updateEloss
void updateEloss(double ioni, double sigi, double rad, double sigr)
Trk::MagneticFieldProperties::magneticFieldMode
MagneticFieldMode magneticFieldMode() const
Returns the MagneticFieldMode as specified.
Trk::STEP_Propagator::Cache::m_brem
bool m_brem
Definition: STEP_Propagator.h:352
Trk::TrackSurfaceIntersection::direction
const Amg::Vector3D & direction() const
Method to retrieve the direction at the Intersection.
Definition: TrackSurfaceIntersection.h:88
Trk::STEP_Propagator::smear
void smear(Cache &cache, double &phi, double &theta, const Trk::TrackParameters *parms, double radDist) const
Definition: STEP_Propagator.cxx:2749
RNGWrapper.h
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:16
Trk::STEP_Propagator::Cache::m_extrapolationCache
Trk::ExtrapolationCache * m_extrapolationCache
Definition: STEP_Propagator.h:403
Trk::STEP_Propagator::Cache::m_matupd_lastpath
double m_matupd_lastpath
Definition: STEP_Propagator.h:360
EventPrimitivesCovarianceHelpers.h
Trk::STEP_Propagator::Cache::m_material
const Material * m_material
cache for collecting the total X0 ans Elos
Definition: STEP_Propagator.h:401
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::BinnedMaterial::layerBinUtility
const Trk::BinUtility * layerBinUtility(const Amg::Vector3D &position) const
access to layer bin utility
Definition: BinnedMaterial.h:80
Trk::STEP_Propagator::Cache::m_sigmaRad
double m_sigmaRad
Definition: STEP_Propagator.h:366
DeMoScan.index
string index
Definition: DeMoScan.py:362
h
python.CaloAddPedShiftConfig.int
int
Definition: CaloAddPedShiftConfig.py:45
ITimedMatEffUpdator.h
Trk::STEP_Propagator::Cache::m_matstates
std::vector< const Trk::TrackStateOnSurface * > * m_matstates
cache of intersections
Definition: STEP_Propagator.h:391
Trk::STEP_Propagator::Cache::m_identifiedParameters
std::vector< std::pair< std::unique_ptr< Trk::TrackParameters >, int > > * m_identifiedParameters
cache of intersections/hit info
Definition: STEP_Propagator.h:394
mapkey::sf
@ sf
Definition: TElectronEfficiencyCorrectionTool.cxx:38
Trk::STEP_Propagator::Cache::m_fieldCache
MagField::AtlasFieldCache m_fieldCache
Definition: STEP_Propagator.h:411
Amg::intersect
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
Definition: GeoPrimitivesHelpers.h:347
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::STEP_Propagator::Cache::m_delRad
double m_delRad
Definition: STEP_Propagator.h:362
Trk::STEP_Propagator::Cache::m_trackingVolume
const TrackingVolume * m_trackingVolume
Definition: STEP_Propagator.h:400
Trk::userOwn
@ userOwn
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:61
Trk::BoundaryCheck
Definition: BoundaryCheck.h:51
RungeKuttaUtils.h
Trk::BinnedMaterial::material
const IdentifiedMaterial * material(const Amg::Vector3D &position) const
access to material/id per bin
Definition: BinnedMaterial.cxx:58
Trk::STEP_Propagator::m_rndGenSvc
ServiceHandle< IAthRNGSvc > m_rndGenSvc
Random Generator service.
Definition: STEP_Propagator.h:523
DeMoScan.first
bool first
Definition: DeMoScan.py:534
Trk::STEP_Propagator::Cache::m_propagateWithPathLimit
int m_propagateWithPathLimit
Definition: STEP_Propagator.h:357
Trk::STEP_Propagator::Cache::m_randomEngine
CLHEP::HepRandomEngine * m_randomEngine
Definition: STEP_Propagator.h:412
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
Trk::SurfaceType::Disc
@ Disc
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
Trk::RungeKuttaUtils::stepEstimator
double stepEstimator(int kind, double *ATH_RESTRICT Su, const double *ATH_RESTRICT P, bool &Q)
Definition: RungeKuttaUtils.cxx:1212
Trk::SurfaceType::Cylinder
@ Cylinder
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
get
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
Definition: hcg.cxx:127
Trk::TimeLimit::time
float time
Definition: HelperStructs.h:60
python.SystemOfUnits.s
float s
Definition: SystemOfUnits.py:147
Trk::photon
@ photon
Definition: ParticleHypothesis.h:35
Trk::TrackingVolume::isAlignable
virtual bool isAlignable() const
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Trk::ConeSurface
Definition: ConeSurface.h:51
Amg::AngleAxis3D
Eigen::AngleAxisd AngleAxis3D
Definition: GeoPrimitives.h:45
Trk::SurfaceType::Plane
@ Plane
Trk::phi
@ phi
Definition: ParamDefs.h:75
Trk::PathLimit::x0Max
float x0Max
Definition: HelperStructs.h:35
ATLAS_THREAD_SAFE
#define ATLAS_THREAD_SAFE
Definition: checker_macros.h:211
Trk::BinnedMaterial::layerBin
size_t layerBin(const Amg::Vector3D &position) const
layer bin
Definition: BinnedMaterial.h:85
Trk::nonInteractingMuon
@ nonInteractingMuon
Definition: ParticleHypothesis.h:39
I
#define I(x, y, z)
Definition: MD5.cxx:116
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
Trk::SurfaceType::Line
@ Line
Amg::hasPositiveOrZeroDiagElems
bool hasPositiveOrZeroDiagElems(const AmgSymMatrix(N) &mat)
Returns true if all diagonal elements of the covariance matrix are finite aka sane in the above defin...
Definition: EventPrimitivesCovarianceHelpers.h:73
MuonParameters::beta
@ beta
Definition: MuonParamDefs.h:144
Trk::STEP_Propagator::Cache::m_matupd_lastmom
double m_matupd_lastmom
Definition: STEP_Propagator.h:359
TrackSurfaceIntersection.h
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
Trk::STEP_Propagator::~STEP_Propagator
virtual ~STEP_Propagator()
Trk::STEP_Propagator::Cache::m_timeOfFlight
double m_timeOfFlight
Definition: STEP_Propagator.h:371
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Trk::TrackingVolume
Definition: TrackingVolume.h:119
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
AlignableTrackingVolume.h
Trk::ExtrapolationCache
Definition: ExtrapolationCache.h:26
Trk::STEP_Propagator::initialize
virtual StatusCode initialize() override final
AlgTool initialize method.
Definition: STEP_Propagator.cxx:1187
Trk::Surface::type
constexpr virtual SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
Trk::STEP_Propagator::Cache::m_maxPath
double m_maxPath
Definition: STEP_Propagator.h:379
TrackStateOnSurface.h
Trk::STEP_Propagator::m_materialEffects
BooleanProperty m_materialEffects
Definition: STEP_Propagator.h:488
TSU::T
unsigned long long T
Definition: L1TopoDataTypes.h:35
fitman.k
k
Definition: fitman.py:528
Trk::STEP_Propagator::finalize
virtual StatusCode finalize() override final
AlgTool finalize method.
Definition: STEP_Propagator.cxx:1214
Trk::ParametersBase::clone
virtual ParametersBase< DIM, T > * clone() const override=0
clone method for polymorphic deep copy
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
Trk::CylinderBounds::halflengthZ
double halflengthZ() const
This method returns the halflengthZ.
Trk::EnergyLoss::sigmaIoni
double sigmaIoni() const
Trk::CylinderBounds::r
virtual double r() const override final
This method returns the radius.
Trk::STEP_Propagator::m_straggling
BooleanProperty m_straggling
Definition: STEP_Propagator.h:501
Trk::AlignableTrackingVolume
Definition: AlignableTrackingVolume.h:36
Trk::STEP_Propagator::propagate
virtual std::unique_ptr< Trk::NeutralParameters > propagate(const Trk::NeutralParameters &, const Trk::Surface &, Trk::PropDirection, const Trk::BoundaryCheck &, bool rC=false) const override final
Main propagation method NeutralParameters.
Definition: STEP_Propagator.cxx:1219