ATLAS Offline Software
GsfExtrapolator.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3  */
4 
13 //
16 //
17 #include "TrkGeometry/Layer.h"
22 #include "TrkSurfaces/Surface.h"
23 //
25 //
30 //
31 #include <boost/container/flat_set.hpp>
32 #include <boost/container/small_vector.hpp>
33 #include <utility>
34 
35 namespace {
36 constexpr bool useBoundaryMaterialUpdate(true);
37 
38 // We have cases where the MultiComponentState is not owned
39 // so we need to set just the cache ptr pointing to that.
40 // In other cases we need to add it to the cache
41 // and make it to point to the last element we push to the m_mcsRecycleBin
42 inline void
43 addMultiComponentToCache(Trk::IMultiStateExtrapolator::Cache& cache,
45  cache.m_mcsRecycleBin.emplace_back(std::move(input));
46  cache.m_stateAtBoundary = &(cache.m_mcsRecycleBin.back());
47 }
48 
49 inline void
50 setRecallInformation(Trk::IMultiStateExtrapolator::Cache& cache,
51  const Trk::Surface& recallSurface,
52  const Trk::Layer& recallLayer,
53  const Trk::TrackingVolume& recallTrackingVolume)
54 {
55  cache.m_recallSurface = &recallSurface;
56  cache.m_recallLayer = &recallLayer;
57  cache.m_recallTrackingVolume = &recallTrackingVolume;
58 }
59 
60 inline void
61 resetRecallInformation(Trk::IMultiStateExtrapolator::Cache& cache)
62 {
63  cache.m_recallSurface = nullptr;
64  cache.m_recallLayer = nullptr;
65  cache.m_recallTrackingVolume = nullptr;
66 }
67 
68 inline void
69 emptyRecycleBins(Trk::IMultiStateExtrapolator::Cache& cache)
70 {
71  // Reset the boundary information
72  cache.m_stateAtBoundary = nullptr;
73  cache.m_navigationParameters = nullptr;
74  cache.m_trackingVolume = nullptr;
75  cache.m_mcsRecycleBin.clear();
76 }
77 
78 int
79 radialDirection(const Trk::MultiComponentState& pars, Trk::PropDirection dir)
80 {
81  // safe inbound/outbound estimation
82  const double prePositionR = pars.begin()->params->position().perp();
83  return (prePositionR > (pars.begin()->params->position() +
84  static_cast<int>(dir) * 0.5 * prePositionR *
85  pars.begin()->params->momentum().unit())
86  .perp())
87  ? -1
88  : 1;
89 }
94 inline bool
95 radialDirectionCheck(const EventContext& ctx,
96  const Trk::IPropagator& prop,
97  const Trk::MultiComponentState& startParm,
98  const Trk::MultiComponentState& parsOnLayer,
99  const Trk::TrackingVolume& tvol,
100  const Trk::MagneticFieldProperties& fieldProperties,
101  const Trk::PropDirection dir)
102 {
103  const Amg::Vector3D& startPosition = startParm.begin()->params->position();
104  const Amg::Vector3D& onLayerPosition = parsOnLayer.begin()->params->position();
105 
106  // the 3D distance to the layer intersection
107  const double distToLayer = (startPosition - onLayerPosition).mag();
108  // get the innermost contained surface for crosscheck
109  const auto& boundarySurfaces = tvol.boundarySurfaces();
110  // only for tubes the crossing makes sense to check for validity
111  if (boundarySurfaces.size() == 4) {
112  // propagate to the inside surface and compare the distance:
113  // it can be either the next layer from the initial point, or the inner tube
114  // boundary surface
115  const Trk::Surface& insideSurface =
116  (boundarySurfaces[Trk::tubeInnerCover])->surfaceRepresentation();
117  auto parsOnInsideSurface =
118  prop.propagateParameters(ctx,
119  *(startParm.begin()->params),
120  insideSurface,
121  dir,
122  true,
123  fieldProperties,
125  const double distToInsideSurface =
126  parsOnInsideSurface
127  ? (startPosition - (parsOnInsideSurface->position())).mag()
128  : 10e10;
129  // the intersection with the original layer is valid if it is before the
130  // inside surface
131  return distToLayer < distToInsideSurface;
132  }
133  return true;
134 }
135 
136 } // end of anonymous namespace
137 
138 /*
139  * AlgTool implementations
140  */
142  const std::string& name,
143  const IInterface* parent)
145 {
146  declareInterface<IMultiStateExtrapolator>(this);
147 }
148 
150 
153 {
154 
155  ATH_CHECK(m_propagator.retrieve());
156  ATH_CHECK(m_navigator.retrieve());
157  ATH_CHECK(m_materialUpdator.retrieve());
158  m_fieldProperties = m_fastField
161 
162  return StatusCode::SUCCESS;
163 }
164 
165 /************************************************************/
166 /*
167  * Implement the public extrapolate methods
168  */
169 /************************************************************/
170 /*
171  * Extrapolate Interface
172  */
175  const EventContext& ctx,
176  Cache& cache,
177  const Trk::MultiComponentState& multiComponentState,
178  const Trk::Surface& surface,
179  Trk::PropDirection direction,
180  const Trk::BoundaryCheck& boundaryCheck) const
181 {
182  if (multiComponentState.empty()) {
183  return {};
184  }
185  return extrapolateImpl(ctx,
186  cache,
187  multiComponentState,
188  surface,
189  direction,
190  boundaryCheck);
191 }
192 
193 /*
194  * Extrapolate Directly method. Does not use a cache
195  */
198  const EventContext& ctx,
199  const Trk::MultiComponentState& multiComponentState,
200  const Trk::Surface& surface,
201  Trk::PropDirection direction,
202  const Trk::BoundaryCheck& boundaryCheck) const
203 {
204  if (multiComponentState.empty()) {
205  return {};
206  }
207 
208  const Trk::TrackingVolume* currentVolume = m_navigator->highestVolume(ctx);
209  if (!currentVolume) {
211  "Current tracking volume could not be determined... returning {}");
212  return {};
213  }
214  return extrapolateDirectlyImpl(ctx,
215  multiComponentState,
216  surface,
217  direction,
218  boundaryCheck);
219 }
220 
221 /************************************************************/
222 /*
223  * Now the implementation of all the internal methods
224  * that perform the actual calculations
225  */
226 /************************************************************/
227 
228 /*
229  * This is the actual extrapolation method implementation
230  */
233  const EventContext& ctx,
234  Cache& cache,
235  const Trk::MultiComponentState& multiComponentState,
236  const Trk::Surface& surface,
237  Trk::PropDirection direction,
238  const Trk::BoundaryCheck& boundaryCheck) const
239 {
240 
241  // Empty the garbage bin
242  emptyRecycleBins(cache);
243  const Trk::Layer* associatedLayer = nullptr;
244  const Trk::TrackingVolume* startVolume = nullptr;
245  const Trk::TrackingVolume* destinationVolume = nullptr;
246  std::unique_ptr<Trk::TrackParameters> referenceParameters =
247  initialiseNavigation(ctx, cache, multiComponentState, surface,
248  associatedLayer, startVolume, destinationVolume,
249  direction);
250 
251  // Bail to direct extrapolation if the direction cannot be determined
252  if (direction == Trk::anyDirection) {
253  return extrapolateDirectlyImpl(ctx,
254  multiComponentState,
255  surface,
256  direction,
257  boundaryCheck);
258  }
259 
260  const Trk::TrackParameters* combinedState =
261  multiComponentState.begin()->params.get();
262 
263  const Trk::MultiComponentState* currentState = &multiComponentState;
264 
265  /* Define the initial distance between destination and current position.
266  Destination should be determined from either
267  - reference parameters (prefered if they exist) or
268  - destination surface
269  */
270  const Amg::Vector3D globalSeparation =
271  referenceParameters
272  ? referenceParameters->position() - combinedState->position()
273  : surface.globalReferencePoint() - combinedState->position();
274  const double initialDistance = globalSeparation.mag();
275  // Clean up memory from combiner. It is no longer needed
276  combinedState = nullptr;
277 
278  /* There are two parts to the extrapolation:
279  - Extrapolate from start point to volume boundary
280  - Extrapolate from volume boundary to destination surface
281  */
282  /*
283  * Extrapolation to destination volume boundary
284  */
285  bool foundFinalBoundary(true);
286  int fallbackOscillationCounter(0);
287  const Trk::TrackingVolume* currentVolume = startVolume;
288  const Trk::TrackingVolume* previousVolume = nullptr;
289 
290  while (currentVolume && currentVolume != destinationVolume) {
291  // Extrapolate to volume boundary
292  extrapolateToVolumeBoundary(ctx,
293  cache,
294  *currentState,
295  associatedLayer,
296  *currentVolume,
297  direction);
298 
299  // New current state is the state extrapolated to the tracking volume
300  // boundary.
301  currentState = cache.m_stateAtBoundary;
302  // The volume that the extrapolation is about to enter into is called the
303  // nextVolume
304  const Trk::TrackingVolume* nextVolume = cache.m_trackingVolume;
305  // Break the loop if the next tracking volume is the same as the current one
306  if (!nextVolume || nextVolume == currentVolume) {
307  foundFinalBoundary = false;
308  break;
309  }
310  // Break the loop if an oscillation is detected
311  if (previousVolume == nextVolume) {
312  ++fallbackOscillationCounter;
313  }
314  if (fallbackOscillationCounter > 10) {
315  foundFinalBoundary = false;
316  break;
317  }
318  // Break the loop if the distance between the surface and the track
319  // parameters has increased
320  combinedState = currentState->begin()->params.get();
321 
322  auto parametersAtDestination =
323  m_propagator->propagateParameters(ctx,
324  *combinedState,
325  surface,
326  direction,
327  false,
328  m_fieldProperties,
330  Amg::Vector3D newDestination;
331  if (parametersAtDestination) {
332  newDestination = parametersAtDestination->position();
333  // delete parametersAtDestination;
334  } else {
335  newDestination = surface.center();
336  }
337 
338  const double revisedDistance =
339  (cache.m_navigationParameters->position() - newDestination).mag();
340 
341  const double distanceChange = std::abs(revisedDistance - initialDistance);
342 
343  if (revisedDistance > initialDistance && distanceChange > 0.01) {
344  foundFinalBoundary = false;
345  break;
346  }
347 
348  combinedState = nullptr;
349  // Initialise the oscillation checker
350  previousVolume = currentVolume;
351  // As the extrapolation is moving into the next volume, the next volume ->
352  // current volume
353  currentVolume = nextVolume;
354  // Associated layer now needs to be reset
355  // if(!entryLayerFound)
356  associatedLayer = nullptr;
357  } // end while loop
358 
359  // catch failures now
360  if (!currentState || (currentVolume != destinationVolume)) {
361  currentState = &multiComponentState;
362  foundFinalBoundary = false;
363  }
364 
365  if (!foundFinalBoundary) {
366  Trk::MultiComponentState bailOutState =
367  m_propagator->multiStatePropagate(ctx,
368  *currentState,
369  surface,
370  m_fieldProperties,
372  boundaryCheck,
374 
375  emptyRecycleBins(cache);
376  return bailOutState;
377  }
378 
379  /*
380  * Extrapolation from volume boundary to surface
381  */
382 
383  // extrapolate inside destination volume
384  Trk::MultiComponentState destinationState =
385  extrapolateInsideVolume(ctx,
386  cache,
387  *currentState,
388  surface,
389  associatedLayer,
390  *currentVolume,
391  direction,
392  boundaryCheck);
393 
394  // FALLBACK POINT: Crisis if extrapolation fails here... As per extrapolation
395  // to volume boundary, in emergency revert to extrapolateDirectly
396 
397  // or we failed to reach the target
398  if (!destinationState.empty() &&
399  &((*(destinationState.begin())).params->associatedSurface()) != &surface) {
400  destinationState.clear();
401  }
402 
403  if (destinationState.empty()) {
404  destinationState = m_propagator->multiStatePropagate(ctx,
405  *currentState,
406  surface,
407  m_fieldProperties,
409  boundaryCheck,
411  }
412  emptyRecycleBins(cache);
413  return destinationState;
414 }
415 
416 /*
417  * Extrapolate Directly method. Does not use a cache
418  */
421  const EventContext& ctx,
422  const Trk::MultiComponentState& multiComponentState,
423  const Trk::Surface& surface,
424  Trk::PropDirection direction,
425  const Trk::BoundaryCheck& boundaryCheck) const
426 {
427  return m_propagator->multiStatePropagate(ctx,
428  multiComponentState,
429  surface,
430  m_fieldProperties,
431  direction,
432  boundaryCheck,
434 }
435 
436 /*
437  * Extrapolate to Volume Boundary!
438  */
439 void
441  const EventContext& ctx,
442  Cache& cache,
443  const Trk::MultiComponentState& multiComponentState,
444  const Trk::Layer* layer,
445  const Trk::TrackingVolume& trackingVolume,
446  Trk::PropDirection direction) const
447 {
448 
449  //We 1st point to the input.
450  //As we move on we might change to point to the
451  //last element held by
452  //cache.m_mcsRecycleBin
453  cache.m_stateAtBoundary = &multiComponentState;
454 
455  const Trk::TrackParameters* combinedState =
456  cache.m_stateAtBoundary->begin()->params.get();
457  const Trk::Layer* associatedLayer = layer;
458 
459  if (!associatedLayer) {
460  // Get entry layer but do not use it as it should have already be hit if it
461  // was desired
462  associatedLayer = trackingVolume.associatedLayer(combinedState->position());
463  associatedLayer = associatedLayer
464  ? associatedLayer
465  : trackingVolume.nextLayer(
466  combinedState->position(),
467  direction * combinedState->momentum().unit(),
468  associatedLayer);
469  }
470  // Only loop over layers if they can be found within the tracking volume
471  else if (trackingVolume.confinedLayers() &&
472  associatedLayer->layerMaterialProperties()) {
473  Trk::MultiComponentState updatedState = m_materialUpdator->postUpdate(
475  *(cache.m_stateAtBoundary),
476  *layer,
477  direction);
478 
479  if (!updatedState.empty()) {
480  addMultiComponentToCache(cache,std::move(updatedState));
481  }
482  }
483 
484  // If an associated surface can be found, extrapolation within the tracking
485  // volume is mandatory This will take extrapolate to the last layer in the
486  // volume
487  if (associatedLayer) {
488  Trk::MultiComponentState nextState = extrapolateFromLayerToLayer(
489  ctx,
490  cache,
491  *(cache.m_stateAtBoundary),
492  trackingVolume,
493  associatedLayer,
494  nullptr,
495  direction);
496  // if we have a next State update the currentState
497  if (!nextState.empty()) {
498  addMultiComponentToCache(cache,std::move(nextState));
499  }
500  }
501 
502  /* =============================================
503  Find the boundary surface using the navigator
504  ============================================= */
505 
506  Trk::NavigationCell nextNavigationCell(nullptr, nullptr);
507 
508  combinedState = cache.m_stateAtBoundary->begin()->params.get();
509 
510  const Trk::TrackingVolume* nextVolume = nullptr;
511  const Trk::TrackParameters* navigationParameters =
513  ? cache.m_navigationParameters.get()
514  : combinedState;
515 
516  nextNavigationCell = m_navigator->nextTrackingVolume(
517  ctx, *m_propagator, *navigationParameters, direction, trackingVolume);
518 
519  nextVolume = nextNavigationCell.nextVolume;
520 
521  std::unique_ptr<Trk::TrackParameters> nextNavigationParameters =
522  std::move(nextNavigationCell.parametersOnBoundary);
523 
524  if (!nextVolume) {
525  // Reset the layer recall
526  resetRecallInformation(cache);
527  }
528 
529  if (useBoundaryMaterialUpdate) {
530  // Check for two things:
531  // 1. If the next volume was found
532  // 2. If there is material associated with the boundary layer.
533  // If so, apply material effects update.
534 
535  // Get layer associated with boundary surface.
536  const Trk::MaterialLayer * layerAtBoundary =
537  (nextNavigationCell.parametersOnBoundary)
538  ? (nextNavigationCell.parametersOnBoundary->associatedSurface())
539  .materialLayer()
540  : nullptr;
541  Trk::MultiComponentState matUpdatedState{};
542  if (nextVolume && layerAtBoundary) {
543  if (layerAtBoundary->layerMaterialProperties()) {
544  matUpdatedState = m_materialUpdator->postUpdate(
546  *(cache.m_stateAtBoundary),
547  *layerAtBoundary,
548  direction);
549  }
550  }
551 
552  // If state has changed due to boundary material, modify state, parameters
553  // accordingly.
554  if (!matUpdatedState.empty()) {
555  addMultiComponentToCache(cache, std::move(matUpdatedState));
556  nextNavigationParameters =
557  cache.m_stateAtBoundary->begin()->params->uniqueClone();
558  }
559  }
560  // Update the rest of the boundary information in the cache
561  cache.m_navigationParameters = std::move(nextNavigationParameters);
562  cache.m_trackingVolume = nextVolume;
563 }
564 
565 /*
566  * Extrapolate inside volume to destination surface!
567  */
570  const EventContext& ctx,
571  Cache& cache,
572  const Trk::MultiComponentState& multiComponentState,
573  const Trk::Surface& surface,
574  const Trk::Layer* layer,
575  const Trk::TrackingVolume& trackingVolume,
576  Trk::PropDirection direction,
577  const Trk::BoundaryCheck& boundaryCheck) const
578 {
579  //curent state is a plainn ptr to keep track
580  const Trk::MultiComponentState* currentState = &multiComponentState;
581 
582  // Retrieve the destination layer
583  // 1. Association
584  const Trk::Layer* destinationLayer = surface.associatedLayer();
585 
586  // 2. Recall and Global Search
587  if (!destinationLayer) {
588  destinationLayer =
589  (&surface == cache.m_recallSurface)
590  ? cache.m_recallLayer
591  : trackingVolume.associatedLayer(surface.globalReferencePoint());
592  }
593 
594  // Retrieve the current layer
595  // Produce a combined state
596  const Trk::TrackParameters* combinedState =
597  currentState->begin()->params.get();
598 
599  const Trk::Layer* associatedLayer = layer;
600 
601  Trk::MultiComponentState updatedState{};
602  if (!associatedLayer) {
603  // Get entry layer but do not use it as it should have already be hit if it
604  // was desired
605  associatedLayer = trackingVolume.associatedLayer(combinedState->position());
606  associatedLayer =
607  associatedLayer
608  ? associatedLayer
609  : trackingVolume.nextLayer(combinedState->position(),
610  direction * combinedState->momentum().unit(),
611  associatedLayer);
612  }
613 
614  else if (associatedLayer != destinationLayer &&
615  trackingVolume.confinedLayers() &&
616  associatedLayer->layerMaterialProperties()) {
617 
618  updatedState = m_materialUpdator->postUpdate(cache.m_materialEffectsCaches,
619  *currentState,
620  *associatedLayer,
621  direction);
622 
623  if (!updatedState.empty()) {
624  // Refresh the current state pointer
625  currentState = &updatedState;
626  }
627  }
628 
629  // Reset combined state target
630  combinedState = nullptr;
631  Trk::MultiComponentState nextState{};
632  if (destinationLayer) {
633  // If there are intermediate layers then additional extrapolations need to
634  // be done
635  if (associatedLayer && associatedLayer != destinationLayer) {
636  nextState = extrapolateFromLayerToLayer(ctx,
637  cache,
638  *currentState,
639  trackingVolume,
640  associatedLayer,
641  destinationLayer,
642  direction);
643 
644  // currentState is now the next
645  if (!nextState.empty()) {
646  // Refresh the current state pointer
647  currentState = &nextState;
648  }
649  }
650  // Final extrapolation to destination surface
651  Trk::MultiComponentState returnState =
652  extrapolateToDestinationLayer(ctx,
653  cache,
654  *currentState,
655  surface,
656  *destinationLayer,
657  associatedLayer,
658  direction,
659  boundaryCheck);
660  // Set the information for the current layer, surface, tracking volume
661  setRecallInformation(cache, surface, *destinationLayer, trackingVolume);
662  return returnState;
663  }
664 
665  // FALLBACK POINT: If no destination layer is found fall-back and extrapolate
666  // directly
667  Trk::MultiComponentState returnState =
668  m_propagator->multiStatePropagate(ctx,
669  *currentState,
670  surface,
671  m_fieldProperties,
672  direction,
673  boundaryCheck,
675 
676  // No destination layer exists so layer recall method cannot be used and
677  // should be reset
678  resetRecallInformation(cache);
679 
680  return returnState;
681 }
682 
683 /*
684  * Extrapolate from Layer to Layer
685  */
688  const EventContext& ctx,
689  Cache& cache,
690  const MultiComponentState& multiComponentState,
691  const TrackingVolume& trackingVolume,
692  const Layer* startLayer,
693  const Layer* destinationLayer,
694  PropDirection direction) const
695 {
696 
697  const Trk::Layer* currentLayer = startLayer;
698  Trk::MultiComponentState currentState{};
699 
700  const Trk::TrackParameters* combinedState =
701  multiComponentState.begin()->params.get();
702  Amg::Vector3D currentPosition = combinedState->position();
703  Amg::Vector3D currentDirection = direction * combinedState->momentum().unit();
704 
705  // No need to extrapolate to start layer, find the next one
706  const Trk::Layer* nextLayer =
707  currentLayer->nextLayer(currentPosition, currentDirection);
708 
709  using LayerSet = boost::container::flat_set<
710  const Trk::Layer*,
711  std::less<const Trk::Layer*>,
712  boost::container::small_vector<const Trk::Layer*, 8>>;
713  LayerSet layersHit;
714 
715  layersHit.insert(currentLayer);
716 
717  // Begin while loop over all intermediate layers
718  while (nextLayer && nextLayer != destinationLayer) {
719  layersHit.insert(nextLayer);
720  // Only extrapolate to an intermediate layer if it requires material
721  // update. Otherwise step over it
722  if (nextLayer->layerMaterialProperties()) {
723  currentState = extrapolateToIntermediateLayer(
724  ctx,
725  cache,
726  !currentState.empty() ? currentState : multiComponentState,
727  *nextLayer,
728  trackingVolume,
729  direction);
730  }
731 
732  if (!currentState.empty()) {
733  combinedState = currentState.begin()->params.get();
734  currentPosition = combinedState->position();
735  currentDirection = direction * combinedState->momentum().unit();
736  }
737 
738  // Find the next layer
739  currentLayer = nextLayer;
740  nextLayer = currentLayer->nextLayer(currentPosition, currentDirection);
741  if (layersHit.find(nextLayer) != layersHit.end()) {
742  break;
743  }
744  }
745 
746  if (destinationLayer && nextLayer != destinationLayer &&
747  !currentState.empty()) {
748  currentState.clear();
749  }
750 
751  return currentState;
752 }
753 
754 /*
755  * Extrapolate to Intermediate Layer
756  */
759  const EventContext& ctx,
760  Cache& cache,
761  const Trk::MultiComponentState& multiComponentState,
762  const Trk::Layer& layer,
763  const Trk::TrackingVolume& trackingVolume,
764  Trk::PropDirection direction) const
765 {
766  const Trk::MultiComponentState* initialState = &multiComponentState;
767 
768  // Propagate over all components
769  Trk::MultiComponentState destinationState =
770  m_propagator->multiStatePropagate(ctx,
771  *initialState,
772  layer.surfaceRepresentation(),
773  m_fieldProperties,
774  direction,
775  true,
777 
778  if (destinationState.empty()) {
779  return {};
780  }
781 
782  // the layer has been intersected
783  // ------------------------------------------------------------------------
784  // check for radial direction change
785  // ---------------------------------------------------------------------
786  const int rDirection = radialDirection(multiComponentState, direction);
787  const int newrDirection = radialDirection(destinationState, direction);
788  if (newrDirection != rDirection) {
789  // it is unfortunate that the cancelling could invalidate the material
790  // collection
791  // reset the nextParameters if the radial change is not allowed
792  // resetting is ok - since the parameters are in the garbage bin already
793  if (!radialDirectionCheck(ctx,
794  *m_propagator,
795  multiComponentState,
796  destinationState,
797  trackingVolume,
798  m_fieldProperties,
799  direction)) {
800  return {};
801  }
802  }
803 
804  /* -------------------------------------
805  Material effects
806  ------------------------------------- */
807 
808  Trk::MultiComponentState updatedState =
809  m_materialUpdator->update(cache.m_materialEffectsCaches,
810  destinationState,
811  layer,
812  direction);
813 
814  if (updatedState.empty()) {
815  return destinationState;
816  }
817 
818  return updatedState;
819 }
820 
821 /*
822  Extrapolate to Destination Layer
823 */
826  const EventContext& ctx,
827  Cache& cache,
828  const Trk::MultiComponentState& multiComponentState,
829  const Trk::Surface& surface,
830  const Trk::Layer& layer,
831  const Trk::Layer* startLayer,
832  Trk::PropDirection direction,
833  const Trk::BoundaryCheck& boundaryCheck) const
834 {
835 
836  const Trk::MultiComponentState* initialState = &multiComponentState;
837  const Trk::TrackParameters* combinedState = nullptr;
838 
839  // Propagate over all components
840  Trk::MultiComponentState destinationState =
841  m_propagator->multiStatePropagate(ctx,
842  multiComponentState,
843  surface,
844  m_fieldProperties,
845  direction,
846  boundaryCheck,
848 
849  // Require a fall-back if the initial state is close to the destination
850  // surface then a fall-back solution is required
851 
852  if (destinationState.empty()) {
853  combinedState = initialState->begin()->params.get();
854  if (surface.isOnSurface(
855  combinedState->position(), true, 0.5 * layer.thickness())) {
856  destinationState = m_propagator->multiStatePropagate(ctx,
857  *initialState,
858  surface,
859  m_fieldProperties,
861  boundaryCheck,
863  }
864  combinedState = nullptr;
865  if (destinationState.empty()) {
866  return {};
867  }
868  }
869 
870  /* ----------------------------------------
871  Material effects
872  ---------------------------------------- */
873 
874  Trk::MultiComponentState updatedState{};
875  if (startLayer != &layer) {
876  updatedState = m_materialUpdator->preUpdate(cache.m_materialEffectsCaches,
877  destinationState,
878  layer,
879  direction);
880  }
881 
882  if (updatedState.empty()) {
883  return destinationState;
884  }
885 
886  return updatedState;
887 }
888 
889 /*
890  * Initialise Navigation
891  */
892 std::unique_ptr<Trk::TrackParameters>
894  const EventContext& ctx,
895  Cache& cache,
896  const Trk::MultiComponentState& multiComponentState,
897  const Trk::Surface& surface,
898  const Trk::Layer*& currentLayer,
899  const Trk::TrackingVolume*& currentVolume,
900  const Trk::TrackingVolume*& destinationVolume,
901  Trk::PropDirection& direction) const
902 {
903  //Get the highest weight parameters. We will just use those
904  const Trk::TrackParameters* combinedState = multiComponentState.begin()->params.get();
905  /* =============================================
906  Look for current volume
907  ============================================= */
908  // 1. See if the current layer is associated with a tracking volume
909  const Trk::Surface* associatedSurface = &(combinedState->associatedSurface());
910  currentLayer = associatedSurface ? associatedSurface->associatedLayer() : currentLayer;
911  currentVolume = currentLayer ? currentLayer->enclosingTrackingVolume() : currentVolume;
912  // If the association method failed
913  if (!currentVolume) {
914  //Try the recall in case the associatedSurface is the recall one
915  if (associatedSurface == cache.m_recallSurface) {
916  currentVolume = cache.m_recallTrackingVolume;
917  currentLayer = cache.m_recallLayer;
918  }
919  // Global search method if this fails
920  else {
921  // If the recall method fails reset the cache
922  resetRecallInformation(cache);
923  currentVolume = m_navigator->volume(ctx, combinedState->position());
924  currentLayer = currentVolume
925  ? currentVolume->associatedLayer(combinedState->position())
926  : nullptr;
927  }
928  }
929  /* =============================================
930  Determine the resolved direction
931  ============================================= */
932  std::unique_ptr<Trk::TrackParameters> referenceParameters =
933  currentVolume ? m_propagator->propagateParameters(
934  ctx, *combinedState, surface, direction, false,
935  m_fieldProperties, Trk::nonInteracting)
936  : nullptr;
937  // Find concrete direction based on reference parameters
938  if (direction == Trk::anyDirection && referenceParameters) {
939  const Amg::Vector3D surfaceDirection(referenceParameters->position() -
940  combinedState->position());
941  direction = (surfaceDirection.dot(combinedState->momentum()) > 0.)
944  }
945  /* =============================================
946  Look for destination volume
947  ============================================= */
948  // See if the destination layer is associated with a tracking volume
949  destinationVolume = surface.associatedLayer()
951  : nullptr;
952  // Association failed
953  if (!destinationVolume) {
954  // See if the cached recall surface is the destination one
955  if (&surface == cache.m_recallSurface) {
956  destinationVolume = cache.m_recallTrackingVolume;
957  } else {
958  // Global search of tracking geometry to find the destination volume
959  destinationVolume = m_navigator->volume(
960  ctx, referenceParameters ? referenceParameters->position()
961  : surface.globalReferencePoint());
962  }
963  }
964  return referenceParameters;
965 }
966 
Trk::GsfExtrapolator::extrapolate
virtual MultiComponentState extrapolate(const EventContext &ctx, Cache &, const MultiComponentState &, const Surface &, PropDirection direction, const BoundaryCheck &boundaryCheck) const override final
Extrapolation method applying material affects.
Definition: GsfExtrapolator.cxx:174
Trk::anyDirection
@ anyDirection
Definition: PropDirection.h:22
make_hlt_rep.pars
pars
Definition: make_hlt_rep.py:90
EnergyLoss.h
ScatteringAngles.h
Trk::TrackingVolume::nextLayer
const Layer * nextLayer(const Amg::Vector3D &gp, const Amg::Vector3D &mom, bool asres=true, bool skipNavLayer=false) const
Return the next Layer if existing, NULL if no next layer corresponds.
Definition: TrackingVolume.cxx:418
Trk::NavigationCell::nextVolume
const TrackingVolume * nextVolume
Definition: INavigator.h:40
Trk::GsfExtrapolator::extrapolateToDestinationLayer
MultiComponentState extrapolateToDestinationLayer(const EventContext &ctx, Cache &cache, const MultiComponentState &, const Surface &, const Layer &, const Layer *, PropDirection direction, const BoundaryCheck &boundaryCheck) const
Final extrapolation step to a destination layer.
Definition: GsfExtrapolator.cxx:825
Trk::NavigationCell::parametersOnBoundary
std::unique_ptr< TrackParameters > parametersOnBoundary
Definition: INavigator.h:42
TrackParameters.h
Trk::MagneticFieldProperties
Definition: MagneticFieldProperties.h:31
Surface.h
Trk::IMultiStateExtrapolator::Cache::m_recallSurface
const Surface * m_recallSurface
< Surface for recall (not owning)
Definition: IMultiStateExtrapolator.h:53
MaterialProperties.h
Trk::GsfExtrapolator::extrapolateFromLayerToLayer
MultiComponentState extrapolateFromLayerToLayer(const EventContext &ctx, Cache &cache, const MultiComponentState &, const TrackingVolume &, const Layer *startLayer, const Layer *destinationLayer, PropDirection direction) const
Layer stepping, stopping at the last layer before destination.
Definition: GsfExtrapolator.cxx:687
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Trk::oppositeMomentum
@ oppositeMomentum
Definition: PropDirection.h:21
Trk::ParametersBase::associatedSurface
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
Trk::GsfExtrapolator::extrapolateDirectly
virtual MultiComponentState extrapolateDirectly(const EventContext &ctx, const MultiComponentState &, const Surface &, PropDirection direction, const BoundaryCheck &boundaryCheck) const override final
Extrapolation method without material effects.
Definition: GsfExtrapolator.cxx:197
Trk::IMultiStateExtrapolator::Cache::m_stateAtBoundary
const MultiComponentState * m_stateAtBoundary
Definition: IMultiStateExtrapolator.h:59
Layer.h
Trk::alongMomentum
@ alongMomentum
Definition: PropDirection.h:20
Trk::Surface::globalReferencePoint
virtual const Amg::Vector3D & globalReferencePoint() const
Returns a global reference point on the surface, for PlaneSurface, StraightLineSurface,...
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Trk::GsfExtrapolator::initialiseNavigation
std::unique_ptr< Trk::TrackParameters > initialiseNavigation(const EventContext &ctx, Cache &cache, const MultiComponentState &initialState, const Surface &surface, const Layer *&associatedLayer, const TrackingVolume *&currentVolume, const TrackingVolume *&destinationVolume, PropDirection &direction) const
Method to initialise navigation parameters including starting state, layer and volume,...
Definition: GsfExtrapolator.cxx:893
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
flat_set
boost::container::flat_set< Key > flat_set
Definition: CheckConfig.cxx:14
Trk::TrackingVolume::confinedLayers
const LayerArray * confinedLayers() const
Return the subLayer array.
MagneticFieldProperties.h
Trk::FastField
@ FastField
call the fast field access method of the FieldSvc
Definition: MagneticFieldMode.h:20
Trk::PropDirection
PropDirection
Definition: PropDirection.h:19
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::IPropagator
Definition: IPropagator.h:55
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
PlotPulseshapeFromCool.input
input
Definition: PlotPulseshapeFromCool.py:106
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Trk::MultiComponentState
std::vector< ComponentParameters > MultiComponentState
Definition: ComponentParameters.h:27
Trk::Layer::nextLayer
const Layer * nextLayer(const Amg::Vector3D &gp, const Amg::Vector3D &udir) const
getting the next/previous Layer if registered - unit for direction vector required
Definition: Layer.cxx:161
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::TrackingVolume::boundarySurfaces
std::vector< std::shared_ptr< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
Definition: TrackingVolume.cxx:809
Trk::IMultiStateExtrapolator::Cache::m_recallTrackingVolume
const TrackingVolume * m_recallTrackingVolume
Definition: IMultiStateExtrapolator.h:57
Trk::IMultiStateExtrapolator::Cache
MultiStateExtrapolator cache class.
Definition: IMultiStateExtrapolator.h:51
Trk::GsfExtrapolator::initialize
virtual StatusCode initialize() override final
Definition: GsfExtrapolator.cxx:152
Trk::FullField
@ FullField
Field is set to be realistic, but within a given Volume.
Definition: MagneticFieldMode.h:21
Trk::GsfExtrapolator::~GsfExtrapolator
virtual ~GsfExtrapolator() override final
beamspotman.dir
string dir
Definition: beamspotman.py:621
Trk::IPropagator::propagateParameters
virtual std::unique_ptr< TrackParameters > propagateParameters(const EventContext &ctx, const TrackParameters &parm, const Surface &sf, PropDirection dir, const BoundaryCheck &bcheck, const MagneticFieldProperties &mprop, ParticleHypothesis particle=pion, bool returnCurv=false, const TrackingVolume *tVol=nullptr) const =0
Main propagation method for parameters only.
Trk::IMultiStateExtrapolator::Cache::m_recallLayer
const Layer * m_recallLayer
Tracking volume for recall (not owning)
Definition: IMultiStateExtrapolator.h:55
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Trk::nonInteracting
@ nonInteracting
Definition: ParticleHypothesis.h:28
MaterialLayer.h
Trk::GsfExtrapolator::extrapolateDirectlyImpl
MultiComponentState extrapolateDirectlyImpl(const EventContext &ctx, const MultiComponentState &, const Surface &, PropDirection direction, const BoundaryCheck &boundaryCheck) const
Implementation of extrapolation without material effects.
Definition: GsfExtrapolator.cxx:420
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::NavigationCell
useful struct for a single navigation cell
Definition: INavigator.h:38
MaterialUpdateMode.h
Trk::Surface::associatedLayer
const Trk::Layer * associatedLayer() const
return the associated Layer
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
TrackingVolume.h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::BoundaryCheck
Definition: BoundaryCheck.h:51
Trk::Layer::enclosingTrackingVolume
const TrackingVolume * enclosingTrackingVolume() const
get the confining TrackingVolume
Trk::Layer::layerMaterialProperties
const LayerMaterialProperties * layerMaterialProperties() const
getting the LayerMaterialProperties including full/pre/post update
Trk::GsfExtrapolator::extrapolateInsideVolume
MultiComponentState extrapolateInsideVolume(const EventContext &ctx, Cache &cache, const MultiComponentState &, const Surface &, const Layer *, const TrackingVolume &, PropDirection direction, const BoundaryCheck &boundaryCheck) const
Definition: GsfExtrapolator.cxx:569
Trk::IMultiStateExtrapolator::Cache::m_materialEffectsCaches
std::vector< GsfMaterial::Combined > m_materialEffectsCaches
Definition: IMultiStateExtrapolator.h:67
Trk::IMultiStateExtrapolator::Cache::m_trackingVolume
const TrackingVolume * m_trackingVolume
keep track of the MultiComponentStates
Definition: IMultiStateExtrapolator.h:63
Trk::MaterialLayer
Definition: MaterialLayer.h:42
GsfConstants.h
Trk::tubeInnerCover
@ tubeInnerCover
Definition: BoundarySurfaceFace.h:39
GsfExtrapolator.h
Trk::IMultiStateExtrapolator::Cache::m_mcsRecycleBin
std::vector< MultiComponentState > m_mcsRecycleBin
Definition: IMultiStateExtrapolator.h:65
AthAlgTool
Definition: AthAlgTool.h:26
IMaterialMixtureConvolution.h
Abstract base class for convolution of material effects.
Trk::IMultiStateExtrapolator::Cache::m_navigationParameters
std::unique_ptr< TrackParameters > m_navigationParameters
Definition: IMultiStateExtrapolator.h:61
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
Trk::TrackingVolume
Definition: TrackingVolume.h:119
Trk::TrackingVolume::associatedLayer
const Layer * associatedLayer(const Amg::Vector3D &gp) const
Return the associated Layer.
Definition: TrackingVolume.cxx:382
Trk::GsfExtrapolator::extrapolateImpl
MultiComponentState extrapolateImpl(const EventContext &ctx, Cache &cache, const MultiComponentState &, const Surface &, PropDirection direction, const BoundaryCheck &boundaryCheck) const
Implementation of main extrapolation method.
Definition: GsfExtrapolator.cxx:232
Trk::GsfExtrapolator::extrapolateToVolumeBoundary
void extrapolateToVolumeBoundary(const EventContext &ctx, Cache &cache, const MultiComponentState &, const Layer *, const TrackingVolume &, PropDirection direction) const
Two primary private extrapolation methods.
Definition: GsfExtrapolator.cxx:440
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
TrackStateOnSurface.h
Trk::GsfExtrapolator::GsfExtrapolator
GsfExtrapolator(const std::string &, const std::string &, const IInterface *)
Constructor with AlgTool parameters.
Definition: GsfExtrapolator.cxx:141
Trk::Layer
Definition: Layer.h:72
Trk::GsfExtrapolator::extrapolateToIntermediateLayer
MultiComponentState extrapolateToIntermediateLayer(const EventContext &ctx, Cache &cache, const MultiComponentState &, const Layer &, const TrackingVolume &, PropDirection direction) const
Single extrapolation step to an intermediate layer.
Definition: GsfExtrapolator.cxx:758