ATLAS Offline Software
RungeKuttaIntersector.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // Runge-Kutta method for tracking a particle through a magnetic field
7 // (c) ATLAS Detector software
9 
10 #include <cmath>
11 #include <iomanip>
20 #include "TrkSurfaces/Surface.h"
21 
22 namespace{
23 
24 bool
25 isTrapped(const double distance,
26  double& previousDistance,
27  unsigned long long& stepsUntilTrapped)
28 {
29  // catch loopers
30  if (distance > previousDistance) {
31  return true;
32  }
33 
34  if (stepsUntilTrapped <= 0) {
35  return true;
36  }
37 
38  --stepsUntilTrapped;
39 
40  // store previous value to check steps are converging (with 10% safety factor)
41  previousDistance = 1.1 * distance;
42 
43  return false;
44 }
45 
46 }//end anonymous namespace
47 namespace Trk
48 {
49 
51  const std::string& name,
52  const IInterface* parent)
53  : base_class (type, name, parent),
54  m_productionMode (true),
55  m_caloR0 (1900.*Gaudi::Units::mm), // r min for calo high field gradient region
56  m_caloR1 (2500.*Gaudi::Units::mm), // r max for calo high field gradient region
57  m_caloR2 (3500.*Gaudi::Units::mm), // r min for calo medium field gradient region
58  m_caloR3 (3700.*Gaudi::Units::mm), // r min for calo outer flux return region
59  m_caloR4 (3800.*Gaudi::Units::mm), // r max for calo medium field gradient region
60  m_caloZ0 (2350.*Gaudi::Units::mm), // z min for calo medium field gradient region
61  m_caloZ1 (2600.*Gaudi::Units::mm), // z min for calo high field gradient region
62  m_caloZ2 (3600.*Gaudi::Units::mm), // z max for calo high field gradient region
63  m_caloZ3 (6000.*Gaudi::Units::mm), // z min for calo outer flux return region
64  m_inDetR0 (400.*Gaudi::Units::mm), // end of central barrel near constant field region
65  m_inDetR1 (350.*Gaudi::Units::mm), // inner radius of middle/outer transition region
66  m_inDetR2 (800.*Gaudi::Units::mm), // outer radius of low field gradient field region
67  m_inDetZ0 (350.*Gaudi::Units::mm), // end of central barrel near constant field region
68  m_inDetZ1 (420.*Gaudi::Units::mm), // start of well behaved transition region
69  m_inDetZ2 (700.*Gaudi::Units::mm), // start of endcap region
70  m_momentumThreshold (1./20.*Gaudi::Units::MeV), // protection against loopers
71  m_momentumWarnThreshold (1./450.*Gaudi::Units::MeV), // warning threshold for intersection failure
72  m_muonR0 (4300.*Gaudi::Units::mm), // inner radius of barrel toroid region
73  m_muonZ0 (6600.*Gaudi::Units::mm), // start of endcap toroid region
74  m_shortStepMax (3.0*Gaudi::Units::mm),
75  m_shortStepMin (10.*Gaudi::Units::nanometer),
76  m_solenoidR (1300.*Gaudi::Units::mm), // r max after coil (will take small steps near coil)
77  m_solenoidZ (3500.*Gaudi::Units::mm), // z end of InDet region
78  m_stepMax0 (8.0*Gaudi::Units::mm),
79  m_stepMax1 (40.0*Gaudi::Units::mm),
80  m_stepMax2 (80.0*Gaudi::Units::mm),
81  m_stepMax3 (160.0*Gaudi::Units::mm),
82  m_stepMax4 (320.0*Gaudi::Units::mm),
83  m_stepsUntilTrapped (2000),
84  m_third (1./3.),
85  m_toroidR0 (1850.0*Gaudi::Units::mm), // endcap toroid central field - inner radius
86  m_toroidR1 (3500.0*Gaudi::Units::mm), // endcap toroid central field - outer radius
87  m_toroidR2 (6000.0*Gaudi::Units::mm), // after inner barrel or outer endcap coil
88  m_toroidR3 (6500.0*Gaudi::Units::mm), // toroid region - radius above which long steps OK
89  m_toroidZ0 (7000.0*Gaudi::Units::mm), // endcap - near iron structure
90  m_toroidZ1 (8000.0*Gaudi::Units::mm), // endcap - high field gradient
91  m_toroidZ2 (8700.0*Gaudi::Units::mm), // barrel - before coil end loop
92  m_toroidZ3 (9100.0*Gaudi::Units::mm), // barrel - nearing coil end loop
93  m_toroidZ4 (9500.0*Gaudi::Units::mm), // endcap toroid full field and barrel coil regions
94  m_toroidZ5 (9800.0*Gaudi::Units::mm), // endcap toroid central field - inner z
95  m_toroidZ6 (11400.0*Gaudi::Units::mm), // endcap toroid central field - outer z
96  m_toroidZ7 (12900.0*Gaudi::Units::mm), // toroid exit fringe fields
97  m_toroidZ8 (14000.0*Gaudi::Units::mm), // essentially out of any toroid fields
98  m_countExtrapolations (0),
99  m_countShortStep (0),
100  m_countStep (0),
101  m_countStepReduction (0)
102 {
103  declareProperty("ProductionMode", m_productionMode);
104 }
105 
108  // print name and package version
109  ATH_MSG_DEBUG( "RungeKuttaIntersector::initialize()" );
110 
111  // initialize base class
112  if (StatusCode::SUCCESS != AlgTool::initialize()) return StatusCode::FAILURE;
113 
114  // retrieve MagneticFieldTool and StraightLineIntersector
116 
117  // productionMode gives steplengths tuned for use in production (adequate precision),
118  // otherwise take very small steps for maximum precision but with a heavy execution penalty
119  if (! m_productionMode){
126  m_stepsUntilTrapped = 20000;
127  }
128 
129  return StatusCode::SUCCESS;
130 }
131 
134 {
136  {
137  msg(MSG::INFO) << "finalized after " << m_countExtrapolations << " extrapolations,";
138  double norm = 1./static_cast<double>(m_countExtrapolations);
139  msg(MSG::INFO) << std::setiosflags(std::ios::fixed)
140  << " taking an average" << std::setw(7) << std::setprecision(1)
141  << norm*static_cast<double>(m_countStep)
142  << " full steps," << std::setw(5) << std::setprecision(2)
143  << norm*static_cast<double>(m_countStepReduction)
144  << " step reductions and" << std::setw(5) << std::setprecision(2)
145  << norm*static_cast<double>(m_countShortStep)
146  << " short final steps";
147  msg(MSG::INFO) << endmsg;
148  }
149 
150  return StatusCode::SUCCESS;
151 }
152 
154 std::optional<Trk::TrackSurfaceIntersection>
156  const TrackSurfaceIntersection& trackIntersection,
157  const double qOverP) const
158 {
159  // trap low momentum
160  if (std::abs(qOverP) > m_momentumThreshold)
161  {
162  ATH_MSG_DEBUG(" trapped as below momentum threshold" );
163  return std::nullopt;
164  }
165 
166  const auto surfaceType = surface.type();
167  if (surfaceType == Trk::SurfaceType::Plane) {
168  return intersectPlaneSurface(static_cast<const PlaneSurface&>(surface),
169  trackIntersection, qOverP);
170  }
171  if (surfaceType == Trk::SurfaceType::Line) {
173  static_cast<const StraightLineSurface&>(surface), trackIntersection,
174  qOverP);
175  }
176  if (surfaceType == Trk::SurfaceType::Cylinder) {
178  static_cast<const CylinderSurface&>(surface), trackIntersection,
179  qOverP);
180  }
181  if (surfaceType == Trk::SurfaceType::Disc) {
182  return intersectDiscSurface(static_cast<const DiscSurface&>(surface),
183  trackIntersection, qOverP);
184  }
185  if (surfaceType == Trk::SurfaceType::Perigee) {
186  return approachPerigeeSurface(static_cast<const PerigeeSurface&>(surface),
187  trackIntersection, qOverP);
188  }
189 
190  ATH_MSG_WARNING( " unrecognized Surface" );
191  return std::nullopt;
192 }
193 
195 std::optional<Trk::TrackSurfaceIntersection>
197  const TrackSurfaceIntersection& trackIntersection,
198  const double qOverP) const
199 {
201  TrackSurfaceIntersection isect = trackIntersection;
202  const Amg::Vector3D& pos = isect.position();
203  const double rStart = pos.perp();
204  const double zStart = pos.z();
205  MagField::AtlasFieldCache fieldCache;
206  initializeFieldCache(fieldCache);
207  Amg::Vector3D fieldValue = field(pos, fieldCache);
208 
209  // straight line distance along track to closest approach to line
210  const Amg::Vector3D& lineDirection = (surface.transform().rotation()).col(2);
211  double stepLength = 0;
212  double distance = distanceToLine (isect, surface.center(),lineDirection, stepLength);
213  unsigned long long stepsUntilTrapped = m_stepsUntilTrapped;
214  double previousDistance = 1.1*distance;
215 
216  // integration loop (step)
217  while (distance > m_shortStepMax)
218  {
219  if (isTrapped(distance, previousDistance, stepsUntilTrapped)) {
220  // trapped
221  if (msgLvl(MSG::DEBUG)) debugFailure (std::move(isect), surface, qOverP,
222  rStart, zStart, true);
223  return std::nullopt;
224  }
225  assignStepLength(isect, stepLength);
226  step(isect, fieldValue, stepLength, qOverP, fieldCache);
227  distance = distanceToLine (isect, surface.center(),lineDirection, stepLength);
228  }
229 
230  // if intersection OK: make short final step to surface - assuming constant field
231  if (distance > m_shortStepMin) shortStep(isect, fieldValue, stepLength, qOverP);
232  return newIntersection(std::move(isect), surface, qOverP,
233  rStart, zStart);
234 }
235 
237 std::optional<Trk::TrackSurfaceIntersection>
239  const TrackSurfaceIntersection& trackIntersection,
240  const double qOverP) const
241 {
243  TrackSurfaceIntersection isect = trackIntersection;
244  const Amg::Vector3D& pos = isect.position();
245  const double rStart = pos.perp();
246  const double zStart = pos.z();
247  MagField::AtlasFieldCache fieldCache;
248  initializeFieldCache(fieldCache);
249  Amg::Vector3D fieldValue = field(pos, fieldCache);
250 
251  // straight line distance along track to closest approach to line
252  const Amg::Vector3D& lineDirection = (surface.transform().rotation()).col(2);
253  double stepLength = 0;
254  double distance = distanceToLine (isect, surface.center(),lineDirection, stepLength);
255  unsigned long long stepsUntilTrapped = m_stepsUntilTrapped;
256  double previousDistance = 1.1*distance;
257 
258  // integration loop (step)
259  while (distance > m_shortStepMax)
260  {
261  if (isTrapped(distance, previousDistance, stepsUntilTrapped)) {
262  // trapped
263  if (msgLvl(MSG::DEBUG)) debugFailure (std::move(isect), surface, qOverP,
264  rStart, zStart, true);
265  return std::nullopt;
266  }
267  assignStepLength(isect, stepLength);
268  step(isect, fieldValue, stepLength, qOverP,fieldCache);
269  distance = distanceToLine (isect, surface.center(),lineDirection, stepLength);
270  }
271 
272  // if intersection OK: make short final step assuming constant field
273  if (distance > m_shortStepMin) shortStep(isect, fieldValue, stepLength, qOverP);
274  return newIntersection(std::move(isect), surface, qOverP,
275  rStart, zStart);
276 }
277 
279 std::optional<Trk::TrackSurfaceIntersection>
281  const TrackSurfaceIntersection& trackIntersection,
282  const double qOverP) const
283 {
285  TrackSurfaceIntersection isect = trackIntersection;
286  const Amg::Vector3D& pos = isect.position();
287  const double rStart = pos.perp();
288  const double zStart = pos.z();
289  MagField::AtlasFieldCache fieldCache;
290  initializeFieldCache(fieldCache);
291  Amg::Vector3D fieldValue = field(pos, fieldCache);
292 
293  // calculate straight line distance along track to intersect with cylinder radius
294  double cylinderRadius = (surface.globalReferencePoint() - surface.center()).perp();
295  Amg::Vector3D offset = surface.center() - isect.position();
296  double rCurrent = offset.perp();
297  double stepLength = 0;
298  double distance = distanceToCylinder (isect, cylinderRadius,rCurrent,offset, stepLength);
299  unsigned long long stepsUntilTrapped = m_stepsUntilTrapped;
300  double previousDistance = 1.1*distance;
301  bool trapped = false;
302 
303  // integration loop (step)
304  while (distance > m_shortStepMax)
305  {
306  if (isTrapped(distance, previousDistance, stepsUntilTrapped)) {
307  trapped = true;
308  break;
309  }
310  assignStepLength(isect, stepLength);
311  double rPrevious= rCurrent;
312  step(isect, fieldValue, stepLength, qOverP,fieldCache);
313  offset = surface.center() - isect.position();
314  rCurrent = offset.perp();
315  double deltaR1 = rCurrent - rPrevious;
316  double deltaR2 = cylinderRadius - rCurrent;
317  double scale = 1.;
318  if (deltaR1 != 0.) scale = deltaR2 / deltaR1;
319  if (scale < 1.)
320  {
321  stepLength *= scale;
322  if (scale < -1.) {
323  trapped = true;
324  break;
325  }
326  }
327  distance = std::abs(stepLength);
328  }
329 
330  // if intersection OK: make short final step assuming constant field
331  if (! trapped
332  && std::abs(cylinderRadius - rCurrent) < m_shortStepMax)
333  {
334  distance = distanceToCylinder (isect, cylinderRadius,rCurrent,offset, stepLength);
335  // protect against divergence (looping)
336  if (distance < m_shortStepMax)
337  {
338  if (distance > m_shortStepMin) shortStep(isect, fieldValue, stepLength, qOverP);
339  return newIntersection(std::move(isect), surface, qOverP,
340  rStart, zStart);
341  }
342  }
343 
344  // trapped
345  if (msgLvl(MSG::DEBUG)) debugFailure (std::move(isect), surface, qOverP,
346  rStart, zStart, true);
347  return std::nullopt;
348 }
349 
351 std::optional<Trk::TrackSurfaceIntersection>
353  const TrackSurfaceIntersection& trackIntersection,
354  const double qOverP) const
355 {
357  TrackSurfaceIntersection isect = trackIntersection;
358  const Amg::Vector3D& pos = isect.position();
359  const double rStart = pos.perp();
360  const double zStart = pos.z();
361  MagField::AtlasFieldCache fieldCache;
362  initializeFieldCache(fieldCache);
363  Amg::Vector3D fieldValue = field(pos, fieldCache);
364 
365  // straight line distance along track to intersect with disc
366  double stepLength = 0;
367  double distance = distanceToDisc (isect, surface.center().z(), stepLength);
368  unsigned long long stepsUntilTrapped = m_stepsUntilTrapped;
369  double previousDistance = 1.1*distance;
370 
371  // integration loop
372  while (std::abs(stepLength) > m_shortStepMax)
373  {
374  if (isTrapped(distance, previousDistance, stepsUntilTrapped)) {
375  // trapped
376  if (msgLvl(MSG::DEBUG)) debugFailure (std::move(isect), surface, qOverP,
377  rStart, zStart, true);
378  return std::nullopt;
379  }
380  assignStepLength(isect, stepLength);
381  step(isect, fieldValue, stepLength, qOverP,fieldCache);
382  distance = distanceToDisc (isect, surface.center().z(), stepLength);
383  }
384 
385  // if intersection OK: make short final step assuming constant field
386  if (std::abs(stepLength) > m_shortStepMin) shortStep(isect, fieldValue, stepLength, qOverP);
387  return newIntersection(std::move(isect), surface, qOverP,
388  rStart, zStart);
389 }
390 
392 std::optional<Trk::TrackSurfaceIntersection>
394  const TrackSurfaceIntersection& trackIntersection,
395  const double qOverP) const
396 {
398  TrackSurfaceIntersection isect = trackIntersection;
399  const Amg::Vector3D& pos = isect.position();
400  const double rStart = pos.perp();
401  const double zStart = pos.z();
402  MagField::AtlasFieldCache fieldCache;
403  initializeFieldCache(fieldCache);
404  Amg::Vector3D fieldValue = field(pos, fieldCache);
405 
406  // straight line distance along track to intersect with plane
407  double stepLength = 0;
408  double distance = distanceToPlane (isect, surface.center(),surface.normal(), stepLength);
409  unsigned long long stepsUntilTrapped = m_stepsUntilTrapped;
410  double previousDistance = 1.1*distance;
411 
412  // integration loop (step)
413  while (distance > m_shortStepMax)
414  {
415  if (isTrapped(distance, previousDistance, stepsUntilTrapped)) {
416  // trapped
417  if (msgLvl(MSG::DEBUG)) debugFailure (std::move(isect), surface, qOverP,
418  rStart, zStart, true);
419  return std::nullopt;
420  }
421  assignStepLength(isect, stepLength);
422  step(isect, fieldValue, stepLength, qOverP, fieldCache);
423  distance = distanceToPlane (isect, surface.center(),surface.normal(), stepLength);
424  }
425 
426  // if intersection OK: make short final step assuming constant field
427  if (distance > m_shortStepMin) shortStep(isect, fieldValue, stepLength, qOverP);
428  return newIntersection(std::move(isect), surface, qOverP,
429  rStart, zStart);
430 }
431 
432 
433 //<<<<<< PRIVATE FUNCTION DEFINITIONS >>>>>>
434 
435 void
437  double& stepLength) const
438 {
439  const Amg::Vector3D& pos = isect.position();
440  const double sinTheta = isect.direction().perp();
441 
442  // step length assigned according to abs value of current rz position
443  // use default (large) value for most regions of InDet and Calo
444  // TODO: consider bounded steps acc region => solves in/outwards problem
445  const double zCurrent = std::abs(pos.z());
446  const double rCurrent = pos.perp();
447  double stepMax = m_stepMax3;
448  // m_stepFlag = 0;
449 
450  // first treat InDet and Solenoid regions
451  // TODO: so far validated out to R = 1.0m
452  if (rCurrent < m_solenoidR && zCurrent < m_solenoidZ) {
453  if (rCurrent < m_inDetR2) {
454  // inner barrel, forward tracks and endcap regions generally happy
455  // with long steps double number of steps for middle barrel and much
456  // of the transition region
457  if (zCurrent < m_inDetZ0 && rCurrent < m_inDetR0) {
458  // inner barrel
459  // stepMax = m_stepMax3;
460  } else if (zCurrent > m_inDetZ2 || sinTheta < 0.35) {
461  // endcap and far forward
462  // stepMax = m_stepMax3;
463  }
464 
465  // double number of steps for most of the remainder
466  else if (zCurrent > m_inDetZ1 && sinTheta < 0.60) {
467  // barrel:endcap transition
468  stepMax = m_stepMax2;
469  } else if (rCurrent > m_inDetR0) {
470  // middle barrel
471  stepMax = m_stepMax2;
472  }
473 
474  // the remaining transition region requires much shorter steps
475  else {
476  if (rCurrent > m_inDetR1) {
477  // middle/outer transition region
478  stepMax = m_stepMax1;
479  } else {
480  // inner transition (anomalous behaviour causing low pt track
481  // rotation)
482  stepMax = m_stepMax0;
483  }
484  }
485  }
486  // take care when moving to region with shorter steps
487 
488  // field less uniform (approaching coil structure)
489  else if (rCurrent < 1000.) // 1000
490  {
491  if (zCurrent > 700.) {
492  // stepMax = m_stepMax3;
493  // exception at higher radii (influence of coil structure)
494  if (rCurrent > 850.)
495  stepMax = m_stepMax2;
496  }
497  // // double number of steps for most of the remainder
498  else if (zCurrent > 420. && sinTheta < 0.60) {
499  // barrel:endcap transition
500  stepMax = m_stepMax2;
501  } else {
502  // decrease step length with increasing radius
503  if (rCurrent > 900.) {
504  // coil structure
505  stepMax = m_stepMax0;
506  } else {
507  // middle
508  stepMax = m_stepMax1;
509  }
510  }
511  } else {
512  // field less uniform in vicinity of solenoid coils
513  stepMax = m_stepMax1;
514  if (zCurrent < 3000.)
515  stepMax = m_stepMax0;
516  }
517  }
518 
519  // secondly treat MuonSpectrometer regions
520  // start with central barrel
521  else if (rCurrent > m_muonR0 && zCurrent < m_toroidZ4) // Z3
522  {
523  // ++m_count0;
524  // optimize - only refresh phi occasionally
525  double period = M_PI / 4.;
526  double phi = pos.phi() + 0.5 * period;
527  if (phi < 0.)
528  phi += 2. * M_PI;
529  int n = static_cast<int>(phi / period);
530  phi -= static_cast<double>(n) * period;
531  if (phi > 0.5 * period)
532  phi = period - phi;
533 
534  // sol 1: 12.1 full steps, 0.34 reduc, 0 short
535  // sol 2: 12.0 full steps, 0.34 reduc, 0 short
536  // sol 3: 11.9 full steps, 0.34 reduc, 0 short (tidy-up)
537  double radius = rCurrent;
538  if (stepLength < 0.)
539  radius -= m_stepMax2;
540  if (zCurrent > m_toroidZ3 && phi < 0.04) {
541  // fringe and coil regions
542  stepMax = m_stepMax1;
543  // m_stepFlag = 22;
544  // ++m_count5;
545  } else if (radius < 5300.) {
546  // fringe and coil regions
547  stepMax = m_stepMax1;
548  // m_stepFlag = 21;
549  // ++m_count1;
550  }
551  // note: long sector (phi>0.175) has max 1/3 outer steps as finishes at
552  // lower r try: phi-extended long sector with stepMax3 for all r > 5500
553  else if (phi > 0.04 ||
554  (radius > 5500. &&
555  (radius < 9000. ||
556  (phi > 0.028 && radius < 9350.)))) // (extended) long sector
557  {
558  // central 'aircore' between coils
559  stepMax = m_stepMax3;
560  // m_stepFlag = 27;
561  // ++m_count2;
562  } else if (radius > 5500. && radius < 9350.) {
563  // outer coil fringe region
564  stepMax = m_stepMax2;
565  // m_stepFlag = 28;
566  // ++m_count3;
567  } else {
568  // outer coil region
569  stepMax = m_stepMax1;
570  // m_stepFlag = 29;
571  // ++m_count4;
572  }
573  } else if (rCurrent > m_muonR0 || zCurrent > m_muonZ0) {
574  if (zCurrent > m_toroidZ4) {
575  if (zCurrent > m_toroidZ8) {
576  // essentially out of any toroid fields
577  stepMax = m_stepMax4;
578  } else if (zCurrent > m_toroidZ7) {
579  // toroid exit fringe fields
580  stepMax = m_stepMax3;
581  } else if (rCurrent > m_toroidR0 && rCurrent < m_toroidR1 &&
582  zCurrent > m_toroidZ5 && zCurrent < m_toroidZ6) {
583  // endcap toroid central field
584  stepMax = m_stepMax3;
585  } else {
586  // endcap toroid high field gradient and barrel coil regions
587  stepMax = m_stepMax1;
588  }
589  } else if (rCurrent > m_toroidR3) {
590  // main region of barrel toroid :
591  // generally take long steps but with care near coil end loop
592  if (zCurrent < m_toroidZ2) {
593  stepMax = m_stepMax4;
594  // if (rCurrent > 9000.) stepMax = m_stepMax3;
595  // if (rCurrent > 9400.) stepMax = m_stepMax2;
596  } else if (zCurrent < m_toroidZ3) {
597  stepMax = m_stepMax3;
598  } else {
599  stepMax = m_stepMax1;
600  }
601  } else if (rCurrent > m_toroidR2) {
602  // toroid field after inner barrel or outer endcap coil
603  stepMax = m_stepMax3;
604  } else if (zCurrent < m_toroidZ0 || zCurrent > m_toroidZ1 ||
605  rCurrent > m_muonR0) {
606  // small steps in toroid high field regions and near iron structures
607  stepMax = m_stepMax1;
608  } else {
609  // endcap toroid entrance fringe fields
610  stepMax = m_stepMax3;
611  }
612  }
613 
614  // finally assign Calo regions
615  else {
616  if ((rCurrent < m_solenoidR) ||
617  (rCurrent > m_caloR0 && rCurrent < m_caloR1) ||
618  rCurrent > m_caloR3 ||
619  (zCurrent > m_caloZ1 && zCurrent < m_caloZ2) ||
620  zCurrent > m_caloZ3) {
621  stepMax = m_stepMax1;
622  } else if ((rCurrent > m_caloR2 && rCurrent < m_caloR4) ||
623  zCurrent > m_caloZ0) {
624  stepMax = m_stepMax2;
625  }
626  }
627 
628  // finally reset stepLength as necessary
629  if (std::abs(stepLength) < stepMax)
630  return;
631  if (stepLength > 0.) {
632  stepLength = stepMax;
633  } else {
634  stepLength = -stepMax;
635  }
636 }
637 
638 void
640  const Surface& surface,
641  const double qOverP,
642  const double rStart,
643  const double zStart,
644  const bool trapped) const
645 {
646  // forget low momentum (prone to looping...)
647  if (std::abs(qOverP) > m_momentumWarnThreshold) return;
648 
649  double pt = 1.E+8;
650  const double sinTheta = isect.direction().perp();
651  if (qOverP != 0.) pt = sinTheta/(qOverP*Gaudi::Units::GeV);
652 
653  const double rCurrent = isect.position().perp();
654 
655  MsgStream log(msgSvc(), name());
656  if (rCurrent > rStart)
657  {
658  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right)
659  << " fail to intersect surface when extrapolating outwards from R,Z"
660  << std::setw(8) << std::setprecision(1) << rStart << ","
661  << std::setw(7) << std::setprecision(0) << zStart << " mm, with pt"
662  << std::setw(7) << std::setprecision(2) << pt << " GeV, direction eta"
663  << std::setw(5) << std::setprecision(2) << isect.direction().eta();
664  }
665  else
666  {
667  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right)
668  << " fail to intersect surface when extrapolating inwards from R,Z"
669  << std::setw(8) << std::setprecision(1) << rStart << ","
670  << std::setw(7) << std::setprecision(0) << zStart << " mm, with pt"
671  << std::setw(7) << std::setprecision(2) << pt << " GeV, direction eta"
672  << std::setw(5) << std::setprecision(2) << isect.direction().eta();
673  }
674  // if (trapped) log << MSG::DEBUG << " looping in mag field ";
675  log << MSG::DEBUG << endmsg;
676 
677  if (dynamic_cast<const PlaneSurface*>(&surface))
678  {
679  double stepLength = 0;
680  (void)distanceToPlane (isect, surface.center(),surface.normal(), stepLength);
681  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right) << " PlaneSurface"
682  << " at R,Z" << std::setw(8) << std::setprecision(1) << surface.center().perp() << ","
683  << std::setw(7) << std::setprecision(0) << surface.center().z()
684  << " at line distance " << std::setw(9) << std::setprecision(1) << stepLength;
685  // if (trapped) log << MSG::DEBUG << " looping in mag field ";
686  log << MSG::DEBUG << endmsg;
687  }
688  else if (dynamic_cast<const CylinderSurface*>(&surface))
689  {
690  double cylinderRadius = (surface.globalReferencePoint() - surface.center()).perp();
691  Amg::Vector3D offset = surface.center() - isect.position();
692  double rCurrent = offset.perp();
693  double stepLength = 0;
694  double distance = distanceToCylinder (isect, cylinderRadius,rCurrent,offset, stepLength);
695  if (distance < m_shortStepMin)
696  {
697  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right)
698  << " closest approach to CylinderSurface at radius "
699  << std::setw(9) << std::setprecision(4) << rCurrent
700  << " mm. Cylinder radius " << std::setw(9) << std::setprecision(4) << cylinderRadius << " mm"
701  << endmsg;
702  }
703  else
704  {
705  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right) << " CylinderSurface"
706  << " radius " << std::setw(6) << std::setprecision(1) << cylinderRadius
707  << " rCurrent " << std::setw(6) << std::setprecision(1) << rCurrent
708  << " distance " << std::setw(6) << std::setprecision(1) << stepLength;
709  if (trapped) log << MSG::DEBUG << " looping in mag field ";
710  log << MSG::DEBUG << endmsg;
711  }
712  }
713  else if (dynamic_cast<const DiscSurface*>(&surface))
714  {
715  double stepLength = 0;
716  (void)distanceToDisc (isect, surface.center().z(), stepLength);
717  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right) << " DiscSurface"
718  << " at R,Z" << std::setw(8) << std::setprecision(1) << surface.center().perp() << ","
719  << std::setw(7) << std::setprecision(0) << surface.center().z()
720  << " at line distance " << std::setw(9) << std::setprecision(1) << stepLength;
721  if (trapped) log << MSG::DEBUG << " looping in mag field ";
722  log << MSG::DEBUG << endmsg;
723  }
724  else if (dynamic_cast<const PerigeeSurface*>(&surface))
725  {
726  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right) << " PerigeeSurface "
727  << endmsg;
728  }
729  else if (dynamic_cast<const StraightLineSurface*>(&surface))
730  {
731  log << MSG::DEBUG << std::setiosflags(std::ios::fixed|std::ios::right) << " StraightLineSurface "
732  << endmsg;
733  }
734  log << MSG::DEBUG << endmsg;
735 }
736 
737 void
739  const Amg::Vector3D& fieldValue,
740  const double stepLength,
741  const double qOverP) const
742 {
743  Amg::Vector3D& pos = isect.position();
744  Amg::Vector3D& dir = isect.direction();
745  const double cOverP = Gaudi::Units::c_light*qOverP;
746 
747  // as step except for const field assumption
748  double stepOverP = 0.5*stepLength*cOverP;
749  Amg::Vector3D product0 = stepOverP*dir.cross(fieldValue);
750  // intermediate point (half way through step)
751  Amg::Vector3D direction1 = dir + product0;
752  Amg::Vector3D product1 = stepOverP*direction1.cross(fieldValue);
753  Amg::Vector3D direction2 = dir + product1;
754  Amg::Vector3D product2 = stepOverP*direction2.cross(fieldValue);
755  // step end point
756  pos += stepLength*(dir + m_third*(product0+product1+product2));
757  Amg::Vector3D direction3 = dir + 2.*product2;
758  Amg::Vector3D product3 = stepOverP*direction3.cross(fieldValue);
759  dir += m_third*(product0+product3 + 2.*(product1+product2));
760  isect.pathlength() += stepLength;
762 }
763 
764 void
767  double& stepLength,
768  const double qOverP,
769  MagField::AtlasFieldCache &fieldCache) const
770 {
771  Amg::Vector3D& pos = isect.position();
772  Amg::Vector3D& dir = isect.direction();
773  const double cOverP = Gaudi::Units::c_light*qOverP;
774 
775  double stepOverP = 0.5*stepLength*cOverP;
776  Amg::Vector3D product0 = stepOverP*dir.cross(fieldValue);
777 
778  // intermediate field look-up point (half way through step)
779  Amg::Vector3D position1 = pos + 0.5*stepLength*(dir + 0.5*product0);
780  Amg::Vector3D fieldValue1 = field(position1, fieldCache);
781  Amg::Vector3D direction1 = dir + product0;
782  Amg::Vector3D product1 = stepOverP*direction1.cross(fieldValue1);
783  Amg::Vector3D direction2 = dir + product1;
784  Amg::Vector3D product2 = stepOverP*direction2.cross(fieldValue1);
785 
786  // field look-up at step end point
787  Amg::Vector3D offsetAtEnd = stepLength*(dir + m_third*(product0+product1+product2));
788  Amg::Vector3D fieldAtEnd = field(pos+offsetAtEnd, fieldCache);
789 
790  // redo with reduced stepLength if non-uniform field derivative detected
791  if ((fieldValue1 - 0.5 * (fieldValue + fieldAtEnd)).mag() > 0.00001 &&
792  std::abs(stepLength) > m_stepMax0) {
793  if (stepLength > 0.) {
794  stepLength = m_stepMax0;
795  } else {
796  stepLength = -m_stepMax0;
797  }
799  stepOverP = 0.5 * stepLength * cOverP;
800  product0 = stepOverP * dir.cross(fieldValue);
801  // intermediate point (half way through step)
802  Amg::Vector3D position1p =
803  pos + 0.5 * stepLength * (dir + 0.5 * product0);
804  Amg::Vector3D fieldValue1p = field(position1p, fieldCache);
805  Amg::Vector3D direction1p = dir + product0;
806  product1 = stepOverP * direction1p.cross(fieldValue1p);
807  Amg::Vector3D direction2p = dir + product1;
808  product2 = stepOverP * direction2p.cross(fieldValue1p);
809  // step end point
810  offsetAtEnd =
811  stepLength * (dir + m_third * (product0 + product1 + product2));
812  fieldAtEnd = field(pos + offsetAtEnd, fieldCache);
813  }
814 
815  Amg::Vector3D direction3 = dir + 2. * product2;
816  Amg::Vector3D product3 = stepOverP * direction3.cross(fieldAtEnd);
817  dir += m_third * (product0 + product3 + 2. * (product1 + product2));
818  fieldValue = fieldAtEnd;
819  pos += offsetAtEnd;
820  isect.pathlength() += stepLength;
821  ++m_countStep;
822 }
823 
824 } // end of namespace
RungeKuttaIntersector.h
Trk::RungeKuttaIntersector::m_muonZ0
const double m_muonZ0
Definition: RungeKuttaIntersector.h:151
Trk::RungeKuttaIntersector::distanceToPlane
double distanceToPlane(const TrackSurfaceIntersection &isect, const Amg::Vector3D &planePosition, const Amg::Vector3D &planeNormal, double &stepLength) const
Definition: RungeKuttaIntersector.h:232
python.AtlRunQueryAMI.period
period
Definition: AtlRunQueryAMI.py:225
Trk::RungeKuttaIntersector::m_toroidR2
const double m_toroidR2
Definition: RungeKuttaIntersector.h:165
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
Trk::RungeKuttaIntersector::m_momentumWarnThreshold
const double m_momentumWarnThreshold
Definition: RungeKuttaIntersector.h:149
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:17
Trk::RungeKuttaIntersector::m_countStepReduction
std::atomic< unsigned long long > m_countStepReduction
Definition: RungeKuttaIntersector.h:181
StraightLineSurface.h
TrackParameters.h
Trk::RungeKuttaIntersector::m_inDetR1
const double m_inDetR1
Definition: RungeKuttaIntersector.h:143
PerigeeSurface.h
Surface.h
Trk::RungeKuttaIntersector::m_third
const double m_third
Definition: RungeKuttaIntersector.h:162
perp
Scalar perp() const
perp method - perpenticular length
Definition: AmgMatrixBasePlugin.h:44
Trk::RungeKuttaIntersector::m_stepMax1
double m_stepMax1
Definition: RungeKuttaIntersector.h:157
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
Trk::RungeKuttaIntersector::m_caloR2
const double m_caloR2
Definition: RungeKuttaIntersector.h:135
Trk::RungeKuttaIntersector::m_countExtrapolations
std::atomic< unsigned long long > m_countExtrapolations
Definition: RungeKuttaIntersector.h:178
Trk::RungeKuttaIntersector::finalize
virtual StatusCode finalize() override
Definition: RungeKuttaIntersector.cxx:133
Trk::RungeKuttaIntersector::m_toroidZ6
const double m_toroidZ6
Definition: RungeKuttaIntersector.h:173
initialize
void initialize()
Definition: run_EoverP.cxx:894
Trk::RungeKuttaIntersector::intersectPlaneSurface
virtual std::optional< TrackSurfaceIntersection > intersectPlaneSurface(const PlaneSurface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP) const override
IIntersector interface method for specific Surface type : PlaneSurface.
Definition: RungeKuttaIntersector.cxx:393
Trk::RungeKuttaIntersector::m_inDetZ2
const double m_inDetZ2
Definition: RungeKuttaIntersector.h:147
Trk::RungeKuttaIntersector::m_toroidZ8
const double m_toroidZ8
Definition: RungeKuttaIntersector.h:175
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
Trk::RungeKuttaIntersector::initialize
virtual StatusCode initialize() override
Definition: RungeKuttaIntersector.cxx:107
Trk::RungeKuttaIntersector::distanceToCylinder
double distanceToCylinder(const TrackSurfaceIntersection &isect, const double cylinderRadius, const double offsetRadius, const Amg::Vector3D &offset, double &stepLength) const
Definition: RungeKuttaIntersector.h:186
Trk::RungeKuttaIntersector::m_stepsUntilTrapped
int m_stepsUntilTrapped
Definition: RungeKuttaIntersector.h:161
Trk::RungeKuttaIntersector::m_stepMax4
double m_stepMax4
Definition: RungeKuttaIntersector.h:160
test_pyathena.pt
pt
Definition: test_pyathena.py:11
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::RungeKuttaIntersector::m_toroidR3
const double m_toroidR3
Definition: RungeKuttaIntersector.h:166
Trk::RungeKuttaIntersector::m_productionMode
bool m_productionMode
Definition: RungeKuttaIntersector.h:130
Trk::RungeKuttaIntersector::RungeKuttaIntersector
RungeKuttaIntersector(const std::string &type, const std::string &name, const IInterface *parent)
Definition: RungeKuttaIntersector.cxx:50
Trk::RungeKuttaIntersector::m_shortStepMax
double m_shortStepMax
Definition: RungeKuttaIntersector.h:152
Trk::TrackSurfaceIntersection
Definition: TrackSurfaceIntersection.h:32
Trk::RungeKuttaIntersector::field
Amg::Vector3D field(const Amg::Vector3D &point, MagField::AtlasFieldCache &fieldCache) const
Definition: RungeKuttaIntersector.h:257
Trk::DiscSurface
Definition: DiscSurface.h:54
Trk::Surface::globalReferencePoint
virtual const Amg::Vector3D & globalReferencePoint() const
Returns a global reference point on the surface, for PlaneSurface, StraightLineSurface,...
Trk::TrackSurfaceIntersection::pathlength
double pathlength() const
Method to retrieve the pathlength propagated till the Intersection.
Definition: TrackSurfaceIntersection.h:96
xAOD::P4Helpers::deltaR2
double deltaR2(double rapidity1, double phi1, double rapidity2, double phi2)
from bare rapidity,phi
Definition: xAODP4Helpers.h:111
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
Trk::RungeKuttaIntersector::m_caloR4
const double m_caloR4
Definition: RungeKuttaIntersector.h:137
Trk::RungeKuttaIntersector::m_caloR0
const double m_caloR0
Definition: RungeKuttaIntersector.h:133
Trk::RungeKuttaIntersector::m_stepMax2
double m_stepMax2
Definition: RungeKuttaIntersector.h:158
Trk::RungeKuttaIntersector::m_momentumThreshold
const double m_momentumThreshold
Definition: RungeKuttaIntersector.h:148
Trk::RungeKuttaIntersector::m_caloR1
const double m_caloR1
Definition: RungeKuttaIntersector.h:134
Trk::RungeKuttaIntersector::m_stepMax3
double m_stepMax3
Definition: RungeKuttaIntersector.h:159
Trk::RungeKuttaIntersector::m_toroidZ3
const double m_toroidZ3
Definition: RungeKuttaIntersector.h:170
Trk::RungeKuttaIntersector::m_toroidZ7
const double m_toroidZ7
Definition: RungeKuttaIntersector.h:174
Trk::RungeKuttaIntersector::m_solenoidZ
const double m_solenoidZ
Definition: RungeKuttaIntersector.h:155
StdJOSetup.msgSvc
msgSvc
Provide convenience handles for various services.
Definition: StdJOSetup.py:36
Trk::RungeKuttaIntersector::m_caloR3
const double m_caloR3
Definition: RungeKuttaIntersector.h:136
Trk::RungeKuttaIntersector::distanceToDisc
double distanceToDisc(const TrackSurfaceIntersection &isect, const double discZ, double &stepLength) const
Definition: RungeKuttaIntersector.h:202
Trk::TrackSurfaceIntersection::position
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
Definition: TrackSurfaceIntersection.h:80
Trk::RungeKuttaIntersector::m_toroidR1
const double m_toroidR1
Definition: RungeKuttaIntersector.h:164
Trk::RungeKuttaIntersector::m_caloZ2
const double m_caloZ2
Definition: RungeKuttaIntersector.h:140
beamspotman.n
n
Definition: beamspotman.py:731
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::CylinderSurface
Definition: CylinderSurface.h:55
Trk::RungeKuttaIntersector::step
void step(TrackSurfaceIntersection &isect, Amg::Vector3D &fieldValue, double &stepLength, const double qOverP, MagField::AtlasFieldCache &fieldCache) const
Definition: RungeKuttaIntersector.cxx:765
Trk::RungeKuttaIntersector::shortStep
void shortStep(TrackSurfaceIntersection &isect, const Amg::Vector3D &fieldValue, const double stepLength, const double qOverP) const
Definition: RungeKuttaIntersector.cxx:738
CylinderSurface.h
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Trk::RungeKuttaIntersector::distanceToLine
double distanceToLine(const TrackSurfaceIntersection &isect, const Amg::Vector3D &linePosition, const Amg::Vector3D &lineDirection, double &stepLength) const
Definition: RungeKuttaIntersector.h:213
Trk::Surface::normal
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::RungeKuttaIntersector::initializeFieldCache
void initializeFieldCache(MagField::AtlasFieldCache &fieldCache) const
Definition: RungeKuttaIntersector.h:246
Trk::RungeKuttaIntersector::m_caloZ1
const double m_caloZ1
Definition: RungeKuttaIntersector.h:139
Trk::RungeKuttaIntersector::m_countShortStep
std::atomic< unsigned long long > m_countShortStep
Definition: RungeKuttaIntersector.h:179
python.SystemOfUnits.nanometer
int nanometer
Definition: SystemOfUnits.py:72
Trk::RungeKuttaIntersector::m_toroidZ5
const double m_toroidZ5
Definition: RungeKuttaIntersector.h:172
beamspotman.dir
string dir
Definition: beamspotman.py:623
Athena::Units
Definition: Units.h:45
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::SurfaceType::Perigee
@ Perigee
Trk::RungeKuttaIntersector::m_inDetZ1
const double m_inDetZ1
Definition: RungeKuttaIntersector.h:146
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
Trk::RungeKuttaIntersector::m_muonR0
const double m_muonR0
Definition: RungeKuttaIntersector.h:150
Trk::RungeKuttaIntersector::assignStepLength
void assignStepLength(const TrackSurfaceIntersection &isect, double &stepLength) const
Definition: RungeKuttaIntersector.cxx:436
Trk::RungeKuttaIntersector::m_fieldCacheCondObjInputKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCacheCondObjInputKey
Definition: RungeKuttaIntersector.h:125
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
Trk::RungeKuttaIntersector::m_solenoidR
const double m_solenoidR
Definition: RungeKuttaIntersector.h:154
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
query_example.col
col
Definition: query_example.py:7
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
Trk::TrackSurfaceIntersection::direction
const Amg::Vector3D & direction() const
Method to retrieve the direction at the Intersection.
Definition: TrackSurfaceIntersection.h:88
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
python.SystemOfUnits.mm
int mm
Definition: SystemOfUnits.py:83
Trk::RungeKuttaIntersector::m_inDetZ0
const double m_inDetZ0
Definition: RungeKuttaIntersector.h:145
taskman.fieldValue
fieldValue
Definition: taskman.py:493
Trk::RungeKuttaIntersector::intersectCylinderSurface
virtual std::optional< TrackSurfaceIntersection > intersectCylinderSurface(const CylinderSurface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP) const override
IIntersector interface method for specific Surface type : CylinderSurface.
Definition: RungeKuttaIntersector.cxx:280
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
Trk::RungeKuttaIntersector::m_toroidR0
const double m_toroidR0
Definition: RungeKuttaIntersector.h:163
PlaneSurface.h
Trk::RungeKuttaIntersector::m_caloZ3
const double m_caloZ3
Definition: RungeKuttaIntersector.h:141
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
Definition: AtlasFieldCache.h:43
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
Trk::SurfaceType::Disc
@ Disc
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
Trk::RungeKuttaIntersector::m_inDetR0
const double m_inDetR0
Definition: RungeKuttaIntersector.h:142
Trk::SurfaceType::Cylinder
@ Cylinder
convertTimingResiduals.offset
offset
Definition: convertTimingResiduals.py:71
DiscSurface.h
Trk::RungeKuttaIntersector::m_caloZ0
const double m_caloZ0
Definition: RungeKuttaIntersector.h:138
Gaudi
=============================================================================
Definition: CaloGPUClusterAndCellDataMonitorOptions.h:273
Trk::RungeKuttaIntersector::m_toroidZ4
const double m_toroidZ4
Definition: RungeKuttaIntersector.h:171
Trk::SurfaceType::Plane
@ Plane
Trk::phi
@ phi
Definition: ParamDefs.h:75
Trk::RungeKuttaIntersector::m_stepMax0
double m_stepMax0
Definition: RungeKuttaIntersector.h:156
Trk::RungeKuttaIntersector::intersectSurface
virtual std::optional< TrackSurfaceIntersection > intersectSurface(const Surface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP) const override
IIntersector interface method for general Surface type.
Definition: RungeKuttaIntersector.cxx:155
Trk::RungeKuttaIntersector::approachPerigeeSurface
virtual std::optional< TrackSurfaceIntersection > approachPerigeeSurface(const PerigeeSurface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP) const override
IIntersector interface method for specific Surface type : PerigeeSurface.
Definition: RungeKuttaIntersector.cxx:196
Trk::RungeKuttaIntersector::m_countStep
std::atomic< unsigned long long > m_countStep
Definition: RungeKuttaIntersector.h:180
Trk::RungeKuttaIntersector::approachStraightLineSurface
virtual std::optional< TrackSurfaceIntersection > approachStraightLineSurface(const StraightLineSurface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP) const override
IIntersector interface method for specific Surface type : StraightLineSurface.
Definition: RungeKuttaIntersector.cxx:238
Trk::SurfaceType::Line
@ Line
Trk::RungeKuttaIntersector::m_inDetR2
const double m_inDetR2
Definition: RungeKuttaIntersector.h:144
Trk::CylinderSurface::globalReferencePoint
virtual const Amg::Vector3D & globalReferencePoint() const override final
Returns a global reference point: For the Cylinder this is Where denotes the averagePhi() of the cy...
Definition: CylinderSurface.cxx:167
Trk::RungeKuttaIntersector::intersectDiscSurface
virtual std::optional< TrackSurfaceIntersection > intersectDiscSurface(const DiscSurface &surface, const TrackSurfaceIntersection &trackIntersection, const double qOverP) const override
IIntersector interface method for specific Surface type : DiscSurface.
Definition: RungeKuttaIntersector.cxx:352
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
Trk::RungeKuttaIntersector::newIntersection
std::optional< TrackSurfaceIntersection > newIntersection(TrackSurfaceIntersection &&isect, const Surface &surface, const double qOverP, const double rStart, const double zStart) const
Definition: RungeKuttaIntersector.h:265
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
Trk::RungeKuttaIntersector::m_toroidZ2
const double m_toroidZ2
Definition: RungeKuttaIntersector.h:169
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
Intersection.h
Trk::RungeKuttaIntersector::m_shortStepMin
const double m_shortStepMin
Definition: RungeKuttaIntersector.h:153
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
Trk::Surface::type
constexpr virtual SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
Trk::RungeKuttaIntersector::m_toroidZ1
const double m_toroidZ1
Definition: RungeKuttaIntersector.h:168
Trk::StraightLineSurface
Definition: StraightLineSurface.h:51
Trk::RungeKuttaIntersector::debugFailure
void debugFailure(TrackSurfaceIntersection &&isect, const Surface &surface, const double qOverP, const double rStart, const double zStart, const bool trapped) const
Definition: RungeKuttaIntersector.cxx:639