ATLAS Offline Software
Loading...
Searching...
No Matches
TimedExtrapolator.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// TimedExtrapolator.cxx, (c) ATLAS Detector software
8
9#include "GaudiKernel/MsgStream.h"
10#include "GaudiKernel/PhysicalConstants.h"
11// Trk include
12#include "TrkExTools/TimedExtrapolator.h"
13#include "TrkExInterfaces/IPropagator.h"
14#include "TrkExInterfaces/IMultipleScatteringUpdator.h"
15#include "TrkExInterfaces/IEnergyLossUpdator.h"
16#include "TrkSurfaces/SurfaceBounds.h"
17#include "TrkSurfaces/DiscBounds.h"
18#include "TrkSurfaces/PerigeeSurface.h"
19#include "TrkSurfaces/StraightLineSurface.h"
20#include "TrkSurfaces/CylinderSurface.h"
21#include "TrkTrack/Track.h"
22#include "TrkGeometry/DetachedTrackingVolume.h"
23#include "TrkGeometry/AlignableTrackingVolume.h"
24#include "TrkGeometry/Layer.h"
25#include "TrkGeometry/MaterialLayer.h"
26#include "TrkGeometry/CylinderLayer.h"
27#include "TrkGeometry/SubtractedCylinderLayer.h"
28#include "TrkGeometry/TrackingGeometry.h"
29#include "TrkGeometry/MaterialProperties.h"
30#include "TrkVolumes/BoundarySurface.h"
31#include "TrkVolumes/BoundarySurfaceFace.h"
32#include "TrkVolumes/CylinderVolumeBounds.h"
33#include "TrkVolumes/Volume.h"
34#include "TrkParticleBase/TrackParticleBase.h"
35#include "TrkEventUtils/TrkParametersComparisonFunction.h"
36#include "TrkDetDescrInterfaces/IDynamicLayerCreator.h"
37#include "TrkDetDescrUtils/GeometrySignature.h"
38#include "TrkMaterialOnTrack/EnergyLoss.h"
39#include "TrkMaterialOnTrack/ScatteringAngles.h"
40#include "TrkParameters/TrackParameters.h"
41#include "TrkGeometry/MagneticFieldProperties.h"
42#include "TrkEventPrimitives/ParticleHypothesis.h"
43
44// for the comparison with a pointer
45#include <cstdint>
46
47// Amg
48#include "EventPrimitives/EventPrimitives.h"
49#include "GeoPrimitives/GeoPrimitives.h"
50
51// Trk
52#include "TrkSurfaces/PlaneSurface.h"
53
54
55// constructor
56Trk::TimedExtrapolator::TimedExtrapolator(const std::string &t, const std::string &n, const IInterface *p) :
57 AthAlgTool(t, n, p),
58 m_subPropagators(Trk::NumberOfSignatures),
59 m_subUpdators(Trk::NumberOfSignatures) {
60 declareInterface<ITimedExtrapolator>(this);
61}
62
63// destructor
64Trk::TimedExtrapolator::~TimedExtrapolator() = default;
65
66// Athena standard methods
67// initialize
68StatusCode
69Trk::TimedExtrapolator::initialize() {
70 m_fieldProperties = m_fastField ? Trk::MagneticFieldProperties(Trk::FastField) : Trk::MagneticFieldProperties(
71 Trk::FullField);
72 if (m_propagators.empty()) {
73 m_propagators.push_back("Trk::RungeKuttaPropagator/DefaultPropagator");
74 }
75 if (m_updators.empty()) {
76 m_updators.push_back("Trk::MaterialEffectsUpdator/DefaultMaterialEffectsUpdator");
77 }
78 if (m_msupdators.empty()) {
79 m_msupdators.push_back("Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator");
80 }
81
82
83 if (!m_propagators.empty()) {
84 if (m_propagators.retrieve().isFailure()) {
85 ATH_MSG_FATAL("Failed to retrieve tool " << m_propagators);
86 return StatusCode::FAILURE;
87 }
88 ATH_MSG_INFO("Retrieved tools " << m_propagators);
89
90 }
91
92
93 // from the number of retrieved propagators set the configurationLevel
94 unsigned int validprop = m_propagators.size();
95
96 if (!validprop) {
97 ATH_MSG_WARNING("None of the defined propagators could be retrieved!");
98 ATH_MSG_WARNING(" Extrapolators jumps back in unconfigured mode, only strategy pattern methods can be used.");
99 } else {
100 m_configurationLevel = validprop - 1;
101 ATH_MSG_VERBOSE("Configuration level automatically set to " << m_configurationLevel);
102 }
103
104 // Get the Navigation AlgTools
105 if (m_navigator.retrieve().isFailure()) {
106 ATH_MSG_FATAL("Failed to retrieve tool " << m_navigator);
107 return StatusCode::FAILURE;
108 }
109 ATH_MSG_INFO("Retrieved tool " << m_navigator);
110
111 // Get the Material Updator
112 if (m_includeMaterialEffects && !m_updators.empty()) {
113 if (m_updators.retrieve().isFailure()) {
114 ATH_MSG_FATAL("None of the defined material updatros could be retrieved!");
115 ATH_MSG_FATAL("No multiple scattering and energy loss material update will be done.");
116 return StatusCode::FAILURE;
117 }
118 ATH_MSG_INFO("Retrieved tools: " << m_updators);
119
120 }
121
122 // from the number of retrieved propagators set the configurationLevel
123 unsigned int validmeuts = m_updators.size();
124
125 // -----------------------------------------------------------
126 // Sanity check 1
127
128 if (m_propNames.empty() && !m_propagators.empty()) {
129 ATH_MSG_DEBUG("Inconsistent setup of Extrapolator, no sub-propagators configured, doing it for you. ");
130 m_propNames.value().push_back(m_propagators[0]->name().substr(8, m_propagators[0]->name().size() - 8));
131 }
132
133 if (m_updatNames.empty() && !m_updators.empty()) {
134 ATH_MSG_DEBUG("Inconsistent setup of Extrapolator, no sub-materialupdators configured, doing it for you. ");
135 m_updatNames.value().push_back(m_updators[0]->name().substr(8, m_updators[0]->name().size() - 8));
136 }
137
138 // -----------------------------------------------------------
139 // Sanity check 2
140 // fill the number of propagator names and updator names up with first one
141 while (int(m_propNames.size()) < int(Trk::NumberOfSignatures)) {
142 m_propNames.value().push_back(m_propNames[0]);
143 }
144 while (int(m_updatNames.size()) < int(Trk::NumberOfSignatures)) {
145 m_updatNames.value().push_back(m_updatNames[0]);
146 }
147 if (validprop && validmeuts) {
148 // Per definition: if configured not found, take the lowest one
149 for (unsigned int isign = 0; int(isign) < int(Trk::NumberOfSignatures); ++isign) {
150 unsigned int index = 0;
151
152 for (unsigned int iProp = 0; iProp < m_propagators.size(); iProp++) {
153 std::string pname = m_propagators[iProp]->name().substr(8, m_propagators[iProp]->name().size() - 8);
154 if (m_propNames[isign] == pname) {
155 index = iProp;
156 }
157 }
158 ATH_MSG_DEBUG(" subPropagator:" << isign << " pointing to propagator: " << m_propagators[index]->name());
159 m_subPropagators[isign] = (index < validprop) ? &(*m_propagators[index]) : &(*m_propagators[Trk::Global]);
160
161 index = 0;
162 for (unsigned int iUp = 0; iUp < m_updators.size(); iUp++) {
163 std::string uname = m_updators[iUp]->name().substr(8, m_updators[iUp]->name().size() - 8);
164 if (m_updatNames[isign] == uname) {
165 index = iUp;
166 }
167 }
168 ATH_MSG_DEBUG(" subMEUpdator:" << isign << " pointing to updator: " << m_updators[index]->name());
169 m_subUpdators[isign] = (index < validmeuts) ? &(*m_updators[index]) : &(*m_updators[Trk::Global]);
170 }
171 } else {
172 ATH_MSG_FATAL("Configuration Problem of Extrapolator: "
173 << " -- At least one IPropagator and IMaterialUpdator instance have to be given.! ");
174 }
175
176
177 m_maxNavigSurf = 1000;
178 m_maxNavigVol = 50;
179
180
181 ATH_MSG_INFO("initialize() successful");
182 return StatusCode::SUCCESS;
183}
184
185// finalize
186StatusCode
187Trk::TimedExtrapolator::finalize() {
188 ATH_MSG_INFO("finalize() successful");
189 return StatusCode::SUCCESS;
190}
191
192std::unique_ptr<const Trk::TrackParameters>
193Trk::TimedExtrapolator::extrapolateWithPathLimit(
194 const Trk::TrackParameters &parm,
195 Trk::PathLimit &pathLim, Trk::TimeLimit &timeLim,
196 Trk::PropDirection dir,
197 Trk::ParticleHypothesis particle,
198 std::vector<Trk::HitInfo> * &hitInfo,
199 Trk::GeometrySignature &nextGeoID,
200 const Trk::TrackingVolume *boundaryVol) const {
201// extrapolation method intended for simulation of particle decay; collects intersections with active layers
202// possible outcomes:1/ returns curvilinear parameters after reaching the maximal path
203// 2/ returns parameters at destination volume boundary
204// 3/ returns 0 ( particle stopped ) but keeps vector of hits
205
206 Trk::TimedExtrapolator::Cache cache(m_maxNavigSurf);
207 ATH_MSG_DEBUG(
208 "M-[" << ++cache.m_methodSequence << "] extrapolateWithPathLimit(...) " << pathLim.x0Max << ", from " << parm.position());
209 ATH_MSG_DEBUG(
210 "M-[" << ++cache.m_methodSequence << "] extrapolateWithPathLimit(...): resolve active layers? " << m_resolveActive);
211
212 if (!m_stepPropagator) {
213 // Get the STEP_Propagator AlgTool
214 if (m_stepPropagator.retrieve().isFailure()) {
215 ATH_MSG_ERROR("Failed to retrieve tool " << m_stepPropagator);
216 ATH_MSG_ERROR("Configure STEP Propagator for extrapolation with path limit");
217 return nullptr;
218 }
219 ATH_MSG_INFO("Retrieved tool " << m_stepPropagator);
220
221 }
222
223 // reset the path ( in x0 !!)
224 cache.m_path = PathLimit(pathLim.x0Max - pathLim.x0Collected, pathLim.process); // collect material locally
225
226 // initialize hit vector
227 cache.m_hitVector = hitInfo;
228
229 // if no input volume, define as highest volume
230 // const Trk::TrackingVolume* destVolume = boundaryVol ? boundaryVol : m_navigator->highestVolume();
231 cache.m_currentStatic = nullptr;
232 if (boundaryVol && !boundaryVol->inside(parm.position(), m_tolerance)) {
233 return nullptr;
234 }
235
236 // extrapolate to destination volume boundary with path limit
237 std::unique_ptr<const Trk::TrackParameters> returnParms =
238 extrapolateToVolumeWithPathLimit(
239 cache, parm, timeLim, dir, particle, nextGeoID, boundaryVol);
240
241 // save actual path on output
242 if (cache.m_path.x0Collected > 0.) {
243 pathLim.updateMat(cache.m_path.x0Collected, cache.m_path.weightedZ / cache.m_path.x0Collected, cache.m_path.l0Collected);
244 }
245
246 if (hitInfo) {
247 ATH_MSG_DEBUG(hitInfo->size() << " identified intersections found");
248 for (auto & ih : *hitInfo) {
249 ATH_MSG_DEBUG("R,z,ID:" << ih.trackParms->position().perp() << ","
250 << ih.trackParms->position().z() << ","
251 << ih.detID);
252 }
253 }
254
255 std::map<const Trk::TrackParameters *, bool>::iterator garbageIter = cache.m_garbageBin.begin();
256 std::map<const Trk::TrackParameters *, bool>::iterator garbageEnd = cache.m_garbageBin.end();
257 for (; garbageIter != garbageEnd; ++garbageIter) if (garbageIter->first) {
258 if(garbageIter->first == returnParms.get()) {
259 auto ret=returnParms->uniqueClone();
260 ATH_MSG_DEBUG(" [+] garbage - at "
261 << positionOutput(garbageIter->first->position())
262 << " parm=" << garbageIter->first
263 << " is the return param. Cloning to" << ret.get());
264 returnParms = std::move(ret);
265 }
266 }
267
268 return returnParms;
269}
270
271std::unique_ptr<const Trk::TrackParameters>
272Trk::TimedExtrapolator::extrapolateToVolumeWithPathLimit(
273 Trk::TimedExtrapolator::Cache &cache,
274 const Trk::TrackParameters &parm,
275 Trk::TimeLimit &timeLim,
276 Trk::PropDirection dir,
277 Trk::ParticleHypothesis particle,
278 Trk::GeometrySignature &nextGeoID,
279 const Trk::TrackingVolume *destVol) const {
280 // returns:
281 // A) curvilinear track parameters if path limit reached
282 // B) boundary parameters (at destination volume boundary)
283
284 // initialize the return parameters vector
285 std::unique_ptr<const Trk::TrackParameters> returnParameters = nullptr;
286 const Trk::TrackParameters *currPar = &parm;
287 const Trk::TrackingVolume *currVol = nullptr;
288 const Trk::TrackingVolume *nextVol = nullptr;
289 std::vector<unsigned int> solutions;
290 const Trk::TrackingVolume *assocVol = nullptr;
291 unsigned int iDest = 0;
292 const EventContext& ctx = Gaudi::Hive::currentContext();
293 ATH_MSG_DEBUG(" [+] start extrapolateToVolumeWithPathLimit - at " << positionOutput(parm.position())<<" parm="<<&parm);
294 // destination volume boundary ?
295 if (destVol && m_navigator->atVolumeBoundary(currPar, destVol, dir, nextVol, m_tolerance) && nextVol != destVol) {
296 return parm.uniqueClone();
297 }
298
299 if (!cache.m_highestVolume) {
300 cache.m_highestVolume = m_navigator->highestVolume(ctx);
301 }
302
303 emptyGarbageBin(cache,&parm);
304 // navigation surfaces
305 if (cache.m_navigSurfs.capacity() > m_maxNavigSurf) {
306 cache.m_navigSurfs.reserve(m_maxNavigSurf);
307 }
308 cache.m_navigSurfs.clear();
309
310 // target volume may not be part of tracking geometry
311 if (destVol) {
312 const Trk::TrackingVolume *tgVol = m_navigator->trackingGeometry(ctx)->trackingVolume(destVol->volumeName());
313 if (!tgVol || tgVol != destVol) {
314 const auto & bounds = destVol->boundarySurfaces();
315 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
316 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
317 cache.m_navigSurfs.emplace_back(&surf, true);
318 }
319 iDest = bounds.size();
320 }
321 }
322
323 // resolve current position
324 bool updateStatic = false;
325 Amg::Vector3D gp = parm.position();
326
327 if (!cache.m_currentStatic || !cache.m_currentStatic->inside(gp, m_tolerance)) {
328 cache.m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
329 updateStatic = true;
330 }
331 if (m_navigator->atVolumeBoundary(currPar, cache.m_currentStatic, dir, nextVol,
332 m_tolerance) && nextVol != cache.m_currentStatic) {
333 // no next volume found --- end of the world
334 if (!nextVol) {
335 ATH_MSG_DEBUG(" [+] Word boundary reached - at " << positionOutput(currPar->position()));
336 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
337 return currPar->uniqueClone();
338 }
339 cache.m_currentStatic = nextVol;
340 updateStatic = true;
341 }
342
343 // current frame volume known-retrieve geoID
344 nextGeoID = cache.m_currentStatic->geometrySignature();
345
346 // resolve active Calo volumes if hit info required
347 if (cache.m_hitVector && nextGeoID == Trk::Calo) {
348 const Trk::AlignableTrackingVolume *alignTV = dynamic_cast<const Trk::AlignableTrackingVolume *> (cache.m_currentStatic);
349 if (alignTV) {
350 Trk::BoundaryTrackParameters boundPar = extrapolateInAlignableTV(cache,*currPar, timeLim, dir, particle, nextGeoID,
351 alignTV);
352 const Trk::TrackParameters *aPar = boundPar.trPar;
353 if (!aPar) {
354 return returnParameters;
355 }
356 throwIntoGarbageBin(cache,aPar);
357 // cache.m_currentStatic = boundPar.exitVol;
358 return extrapolateToVolumeWithPathLimit(cache,*aPar, timeLim, dir, particle, nextGeoID, destVol);
359 }
360 }
361
362 // update if new static volume
363 if (updateStatic) { // retrieve boundaries
364 cache.m_staticBoundaries.clear();
365 const auto& bounds = cache.m_currentStatic->boundarySurfaces();
366 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
367 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
368 cache.m_staticBoundaries.emplace_back(&surf, true);
369 }
370
371 cache.m_detachedVols.clear();
372 cache.m_detachedBoundaries.clear();
373 cache.m_denseVols.clear();
374 cache.m_denseBoundaries.clear();
375 cache.m_layers.clear();
376 cache.m_navigLays.clear();
377
378 // new: ID volumes may have special material layers ( entry layers ) - add them here
379 // if (cache.m_currentStatic->entryLayerProvider()) {
380 // const std::vector<const Trk::Layer*>& entryLays = cache.m_currentStatic->entryLayerProvider()->layers();
381 // for (unsigned int i=0; i < entryLays.size(); i++) {
382 // if (entryLays[i]->layerType()>0 || entryLays[i]->layerMaterialProperties()) {
383 // cache.m_layers.push_back(std::pair<const
384 // Trk::Surface*,Trk::BoundaryCheck>(&(entryLays[i]->surfaceRepresentation()),true));
385 // cache.m_navigLays.push_back(std::pair<const Trk::TrackingVolume*,const Trk::Layer*> (cache.m_currentStatic,entryLays[i])
386 // );
387 // Trk::DistanceSolution distSol = cache.m_layers.back().first->straightLineDistanceEstimate(currPar->position(),
388 //
389 //
390 //
391 // currPar->momentum().normalized());
392 // }
393 // }
394 // }
395
396 // detached volume boundaries
397 Trk::ArraySpan<const Trk::DetachedTrackingVolume* const> detVols = cache.m_currentStatic->confinedDetachedVolumes();
398 if (!detVols.empty()) {
399 Trk::ArraySpan<const Trk::DetachedTrackingVolume* const>::iterator iTer = detVols.begin();
400 for (; iTer != detVols.end(); ++iTer) {
401 // active station ?
402 const Trk::Layer *layR = (*iTer)->layerRepresentation();
403 bool active = layR && layR->layerType();
404 const auto& detBounds = (*iTer)->trackingVolume()->boundarySurfaces();
405 if (active) {
406 cache.m_detachedVols.emplace_back(*iTer,
407 detBounds.size());
408 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
409 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
410 cache.m_detachedBoundaries.emplace_back(&surf, true);
411 }
412 } else if (cache.m_currentStatic->geometrySignature() != Trk::MS ||
413 !m_useMuonMatApprox ||
414 (*iTer)->name().substr((*iTer)->name().size() - 4, 4) ==
415 "PERM") { // retrieve
416 // inert
417 // detached
418 // objects
419 // only if
420 // needed
421 if ((*iTer)->trackingVolume()->zOverAtimesRho() != 0. &&
422 ((*iTer)->trackingVolume()->confinedDenseVolumes().empty())
423 && (*iTer)->trackingVolume()->confinedArbitraryLayers().empty()) {
424 cache.m_denseVols.emplace_back((*iTer)->trackingVolume(), detBounds.size());
425 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
426 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
427 cache.m_denseBoundaries.emplace_back(&surf, true);
428 }
429 }
430 Trk::ArraySpan<const Trk::Layer* const> confLays = (*iTer)->trackingVolume()->confinedArbitraryLayers();
431 if (!(*iTer)->trackingVolume()->confinedDenseVolumes().empty() || (confLays.size() > detBounds.size())) {
432 cache.m_detachedVols.emplace_back(*iTer, detBounds.size());
433 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
434 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
435 cache.m_detachedBoundaries.emplace_back(&surf, true);
436 }
437 } else if (!confLays.empty()) {
438 for (const Trk::Layer* const lIt : confLays) {
439 cache.m_layers.emplace_back(&(lIt->surfaceRepresentation()),
440 true);
441 cache.m_navigLays.emplace_back((*iTer)->trackingVolume(), lIt);
442 }
443 }
444 }
445 }
446 }
447 cache.m_denseResolved = std::pair<unsigned int, unsigned int> (cache.m_denseVols.size(), cache.m_denseBoundaries.size());
448 cache.m_layerResolved = cache.m_layers.size();
449 }
450
451 cache.m_navigSurfs.insert(cache.m_navigSurfs.end(), cache.m_staticBoundaries.begin(), cache.m_staticBoundaries.end());
452
453 // resolve the use of dense volumes
454 cache.m_dense = (cache.m_currentStatic->geometrySignature() == Trk::MS && m_useMuonMatApprox) ||
455 (cache.m_currentStatic->geometrySignature() != Trk::MS && m_useDenseVolumeDescription);
456
457 // reset remaining counters
458 cache.m_currentDense = cache.m_dense ? cache.m_currentStatic : cache.m_highestVolume;
459 cache.m_navigBoundaries.clear();
460 if (cache.m_denseVols.size() > cache.m_denseResolved.first) {
461 cache.m_denseVols.resize(cache.m_denseResolved.first);
462 }
463 while (cache.m_denseBoundaries.size() > cache.m_denseResolved.second) {
464 cache.m_denseBoundaries.pop_back();
465 }
466 if (cache.m_layers.size() > cache.m_layerResolved) {
467 cache.m_navigLays.resize(cache.m_layerResolved);
468 }
469 while (cache.m_layers.size() > cache.m_layerResolved) {
470 cache.m_layers.pop_back();
471 }
472
473 // current detached volumes
474 // collect : subvolume boundaries, ordered/unordered layers, confined dense volumes
476 // const Trk::DetachedTrackingVolume* currentActive = 0;
477 std::vector<std::pair<const Trk::TrackingVolume *, unsigned int> > navigVols;
478
479 gp = currPar->position();
480 std::vector<const Trk::DetachedTrackingVolume *> detVols =
481 m_navigator->trackingGeometry(ctx)->lowestDetachedTrackingVolumes(gp);
482 std::vector<const Trk::DetachedTrackingVolume *>::iterator dIter = detVols.begin();
483 for (; dIter != detVols.end(); ++dIter) {
484 const Trk::Layer *layR = (*dIter)->layerRepresentation();
485 bool active = layR && layR->layerType();
486 if (active && !m_resolveActive) {
487 continue;
488 }
489 if (!active && cache.m_currentStatic->geometrySignature() == Trk::MS &&
490 m_useMuonMatApprox && (*dIter)->name().substr((*dIter)->name().size() - 4, 4) != "PERM") {
491 continue;
492 }
493 const Trk::TrackingVolume *dVol = (*dIter)->trackingVolume();
494 // detached volume exit ?
495 bool dExit = m_navigator->atVolumeBoundary(currPar, dVol, dir, nextVol, m_tolerance) && !nextVol;
496 if (dExit) {
497 continue;
498 }
499 // inert material
500 const auto& confinedDense = dVol->confinedDenseVolumes();
501 const auto& confinedLays = dVol->confinedArbitraryLayers();
502
503 if (!active && confinedDense.empty() && confinedLays.empty()) {
504 continue;
505 }
506 const auto &bounds = dVol->boundarySurfaces();
507 if (!active && confinedDense.empty() && confinedLays.size() <= bounds.size()) {
508 continue;
509 }
510 if (!confinedDense.empty() || !confinedLays.empty()) {
511 navigVols.emplace_back(dVol, bounds.size());
512 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
513 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
514 cache.m_navigBoundaries.emplace_back(&surf, true);
515 }
516 // collect dense volume boundary
517 if (!confinedDense.empty()) {
518 auto vIter = confinedDense.begin();
519 for (; vIter != confinedDense.end(); ++vIter) {
520 const auto& bounds = (*vIter)->boundarySurfaces();
521 cache.m_denseVols.emplace_back(*vIter, bounds.size());
522 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
523 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
524 cache.m_denseBoundaries.emplace_back(&surf, true);
525 }
526 }
527 }
528 // collect unordered layers
529 if (!confinedLays.empty()) {
530 for (const auto *confinedLay : confinedLays) {
531 cache.m_layers.emplace_back(&(confinedLay->surfaceRepresentation()), true);
532 cache.m_navigLays.emplace_back(dVol, confinedLay);
533 }
534 }
535 } else { // active material
536 const Trk::TrackingVolume *detVol = dVol->associatedSubVolume(gp);
537 if (!detVol && dVol->confinedVolumes()) {
538 std::span<Trk::TrackingVolume const * const> subvols = dVol->confinedVolumes()->arrayObjects();
539 for (const auto *subvol : subvols) {
540 if (subvol->inside(gp, m_tolerance)) {
541 detVol = subvol;
542 break;
543 }
544 }
545 }
546
547 if (!detVol) {
548 detVol = dVol;
549 }
550 bool vExit = m_navigator->atVolumeBoundary(currPar, detVol, dir, nextVol, m_tolerance) && nextVol != detVol;
551 if (vExit && nextVol && nextVol->inside(gp, m_tolerance)) {
552 detVol = nextVol;
553 vExit = false;
554 }
555 if (!vExit) {
556 const auto &bounds = detVol->boundarySurfaces();
557 navigVols.emplace_back(detVol, bounds.size());
558 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
559 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
560 cache.m_navigBoundaries.emplace_back(&surf, true);
561 }
562 if (detVol->zOverAtimesRho() != 0.) {
563 cache.m_denseVols.emplace_back(detVol, bounds.size());
564 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
565 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
566 cache.m_denseBoundaries.emplace_back(&surf, true);
567 }
568 }
569 // layers ?
570 if (detVol->confinedLayers()) {
571 if (m_robustSampling || cache.m_currentStatic->geometrySignature() == Trk::MS) {
572 std::span<Trk::Layer const * const> cLays = detVol->confinedLayers()->arrayObjects();
573 for (const auto *cLay : cLays) {
574 if (cLay->layerType() > 0 || cLay->layerMaterialProperties()) {
575 cache.m_layers.emplace_back(&(cLay->surfaceRepresentation()), true);
576 cache.m_navigLays.emplace_back(cache.m_currentStatic,
577 cLay);
578 }
579 }
580 } else {
581 const Trk::Layer *lay = detVol->associatedLayer(gp);
582 // if (lay && ( (*dIter)->layerRepresentation()
583 // &&(*dIter)->layerRepresentation()->layerType()>0 ) ) currentActive=(*dIter);
584 if (lay) {
585 cache.m_layers.emplace_back(&(lay->surfaceRepresentation()),
586 true);
587 cache.m_navigLays.emplace_back(detVol, lay);
588 }
589 const Trk::Layer *nextLayer = detVol->nextLayer(currPar->position(),
590 dir * currPar->momentum().normalized(), true);
591 if (nextLayer && nextLayer != lay) {
592 cache.m_layers.emplace_back(&(nextLayer->surfaceRepresentation()), true);
593 cache.m_navigLays.emplace_back(detVol, nextLayer);
594 }
595 }
596 } else if (!detVol->confinedArbitraryLayers().empty()) {
597 Trk::ArraySpan<const Trk::Layer* const> layers = detVol->confinedArbitraryLayers();
598 for (const auto *layer : layers) {
599 cache.m_layers.emplace_back(&(layer->surfaceRepresentation()), true);
600 cache.m_navigLays.emplace_back(detVol, layer);
601 }
602 }
603 }
604 }
605 }
606
607 // confined layers
608 if (cache.m_currentStatic->confinedLayers() && updateStatic) {
609 // if ( cache.m_currentStatic->confinedLayers() ) {
610 if (m_robustSampling || cache.m_currentStatic->geometrySignature() == Trk::MS) {
611 std::span<Trk::Layer const * const> cLays = cache.m_currentStatic->confinedLayers()->arrayObjects();
612 for (const auto *cLay : cLays) {
613 if (cLay->layerType() > 0 || cLay->layerMaterialProperties()) {
614 cache.m_layers.emplace_back(&(cLay->surfaceRepresentation()),
615 true);
616 cache.m_navigLays.emplace_back(cache.m_currentStatic, cLay);
617 }
618 }
619 } else {
620 // * this does not work - debug !
621 const Trk::Layer *lay = cache.m_currentStatic->associatedLayer(gp);
622 // if (!lay) {
623 // lay = cache.m_currentStatic->associatedLayer(gp+m_tolerance*parm.momentum().normalized());
624 // std::cout<<" find input associated layer, second attempt:"<< lay<< std::endl;
625 // }
626 if (lay) {
627 cache.m_layers.emplace_back(&(lay->surfaceRepresentation()),
628 Trk::BoundaryCheck(false));
629 cache.m_navigLays.emplace_back(cache.m_currentStatic, lay);
630 const Trk::Layer *nextLayer = lay->nextLayer(currPar->position(), dir * currPar->momentum().normalized());
631 if (nextLayer && nextLayer != lay) {
632 cache.m_layers.emplace_back(&(nextLayer->surfaceRepresentation()),
633 Trk::BoundaryCheck(false));
634 cache.m_navigLays.emplace_back(cache.m_currentStatic,
635 nextLayer);
636 }
637 const Trk::Layer *backLayer = lay->nextLayer(currPar->position(), -dir * currPar->momentum().normalized());
638 if (backLayer && backLayer != lay) {
639 cache.m_layers.emplace_back(&(backLayer->surfaceRepresentation()),
640 Trk::BoundaryCheck(false));
641 cache.m_navigLays.emplace_back(cache.m_currentStatic,
642 backLayer);
643 }
644 }
645 }
646 }
647
648 // cache.m_navigSurfs contains destination surface (if it exists), static volume boundaries
649 // complete with TG cache.m_layers/dynamic layers, cache.m_denseBoundaries, cache.m_navigBoundaries, m_detachedBoundaries
650
651 if (!cache.m_layers.empty()) {
652 cache.m_navigSurfs.insert(cache.m_navigSurfs.end(), cache.m_layers.begin(), cache.m_layers.end());
653 }
654 if (!cache.m_denseBoundaries.empty()) {
655 cache.m_navigSurfs.insert(cache.m_navigSurfs.end(), cache.m_denseBoundaries.begin(), cache.m_denseBoundaries.end());
656 }
657 if (!cache.m_navigBoundaries.empty()) {
658 cache.m_navigSurfs.insert(cache.m_navigSurfs.end(), cache.m_navigBoundaries.begin(), cache.m_navigBoundaries.end());
659 }
660 if (!cache.m_detachedBoundaries.empty()) {
661 cache.m_navigSurfs.insert(cache.m_navigSurfs.end(), cache.m_detachedBoundaries.begin(), cache.m_detachedBoundaries.end());
662 }
663
664
665 // current dense
666 cache.m_currentDense = cache.m_highestVolume;
667 if (cache.m_dense && cache.m_denseVols.empty()) {
668 cache.m_currentDense = cache.m_currentStatic;
669 } else {
670 for (unsigned int i = 0; i < cache.m_denseVols.size(); i++) {
671 const Trk::TrackingVolume *dVol = cache.m_denseVols[i].first;
672 if (dVol->inside(currPar->position(), m_tolerance) && dVol->zOverAtimesRho() != 0.) {
673 if (!m_navigator->atVolumeBoundary(currPar, dVol, dir, nextVol, m_tolerance) || nextVol == dVol) {
674 cache.m_currentDense = dVol;
675 }
676 }
677 }
678 }
679
680 // before propagation, loop over layers and collect hits
681 if (cache.m_hitVector) {
682 for (unsigned int i = 0; i < cache.m_navigLays.size(); i++) {
683 if (cache.m_navigLays[i].second->layerType() > 0 && cache.m_navigLays[i].second->isOnLayer(currPar->position())) {
684 if (cache.m_navigLays[i].second->surfaceArray()) {
685 // perform the overlap Search on this layer
686 ATH_MSG_VERBOSE(" [o] Calling overlapSearch() on input layer.");
687 overlapSearch(cache,*m_subPropagators[0], *currPar, *currPar, *cache.m_navigLays[i].second, timeLim.time, dir, true,
688 particle);
689 } else {
690 ATH_MSG_VERBOSE(" [o] Collecting intersection with active input layer.");
691 cache.m_hitVector->emplace_back(currPar->uniqueClone(), timeLim.time, cache.m_navigLays[i].second->layerType(), 0.);
692 }
693 } // ------------------------------------------------- Fatras mode off -----------------------------------
694 }
695 }
696
697 // ready to propagate
698 // till: A/ static volume boundary(bcheck=true) , B/ material layer(bcheck=true), C/ destination surface(bcheck=false)
699 // update of cache.m_navigSurfs required if I/ entry into new navig volume, II/ exit from currentActive without overlaps
700
701 nextVol = nullptr;
702 while (currPar) {
703 std::vector<unsigned int> solutions;
704 // double time_backup = timeLim.time;
705 // double path_backup = cache.m_path.x0Collected;
706 ATH_MSG_DEBUG(" [+] Starting propagation at position " << positionOutput(currPar->position())
707 << " (current momentum: " << currPar->momentum().mag() <<
708 ")");
709 ATH_MSG_DEBUG(" [+] " << cache.m_navigSurfs.size() << " target surfaces in '" << cache.m_currentDense->volumeName() << "'."); //
710 // verify
711 // that
712 // material
713 // input
714 // makes
715 // sense
716 if (!(cache.m_currentDense->inside(currPar->position(), m_tolerance)
717 || m_navigator->atVolumeBoundary(currPar, cache.m_currentDense, dir, assocVol, m_tolerance))) {
718 cache.m_currentDense = cache.m_highestVolume;
719 }
720 const Trk::TrackParameters* nextPar = m_stepPropagator
721 ->propagateT(ctx,
722 *currPar,
723 cache.m_navigSurfs,
724 dir,
725 m_fieldProperties,
726 particle,
727 solutions,
728 cache.m_path,
729 timeLim,
730 true,
731 cache.m_currentDense,
732 cache.m_hitVector)
733 .release();
734 ATH_MSG_VERBOSE(" [+] Propagation done. ");
735 if (nextPar) {
736 ATH_MSG_DEBUG(" [+] Position after propagation - at " << positionOutput(
737 nextPar->position()) << ", timed at " << timeLim.time);
738 }
739
740 if (!nextPar) {
741 ATH_MSG_DEBUG(" [!] Propagation failed, return 0");
742 cache.m_parametersAtBoundary.boundaryInformation(cache.m_currentStatic, nextPar, nextPar);
743 return returnParameters;
744 }
745
746 throwIntoGarbageBin(cache,nextPar);
747
748 // material update has been done already by the propagator
749 if (cache.m_path.x0Max > 0. &&
750 ((cache.m_path.process < 100 && cache.m_path.x0Collected >= cache.m_path.x0Max) ||
751 (cache.m_path.process > 100 && cache.m_path.l0Collected >= cache.m_path.x0Max))) {
752 // trigger presampled interaction, provide material properties if needed
753 // process interaction only if creation of secondaries allowed
754 if (cache.m_currentStatic->geometrySignature() == Trk::ID || m_caloMsSecondary) {
755 const Trk::Material *extMprop = cache.m_path.process > 100 ? cache.m_currentDense : nullptr;
756
757 const Trk::TrackParameters* iPar =
758 m_updators[0]
759 ->interact(
760 timeLim.time, nextPar->position(), nextPar->momentum(), particle, cache.m_path.process, extMprop)
761 .release();
762
763 if (!iPar) {
764 return returnParameters;
765 }
766
767 throwIntoGarbageBin(cache,iPar);
768 return extrapolateToVolumeWithPathLimit(cache,*iPar, timeLim, dir, particle, nextGeoID, destVol);
769 } // kill the particle without trace ( some validation info can be included here eventually )
770 return returnParameters;
771
772 }
773 // decay ?
774 if (timeLim.tMax > 0. && timeLim.time >= timeLim.tMax) {
775 // process interaction only if creation of secondaries allowed
776 if (cache.m_currentStatic->geometrySignature() == Trk::ID || m_caloMsSecondary) {
777 // trigger presampled interaction
778 const Trk::TrackParameters* iPar =
779 m_updators[0]
780 ->interact(timeLim.time, nextPar->position(), nextPar->momentum(), particle, timeLim.process)
781 .release();
782 if (!iPar) {
783 return returnParameters;
784 }
785 throwIntoGarbageBin(cache,iPar);
786 return extrapolateToVolumeWithPathLimit(cache,*iPar, timeLim, dir, particle, nextGeoID, destVol);
787 } // kill the particle without trace ( some validation info can be included here eventually )
788 return returnParameters;
789
790 }
791
792 // check missing volume boundary
793 if (nextPar && !(cache.m_currentDense->inside(nextPar->position(), m_tolerance)
794 || m_navigator->atVolumeBoundary(nextPar, cache.m_currentDense, dir, assocVol, m_tolerance))) {
795 ATH_MSG_DEBUG(" [!] ERROR: missing volume boundary for volume" << cache.m_currentDense->volumeName());
796 }
797
798
799 ATH_MSG_DEBUG(" [+] Number of intersection solutions: " << solutions.size());
800
801 unsigned int iSol = 0;
802 while (iSol < solutions.size()) {
803 if (solutions[iSol] < iDest) {
804 return nextPar->uniqueClone();
805 } if (solutions[iSol] < iDest + cache.m_staticBoundaries.size()) {
806 // material attached ?
807 const Trk::Layer *mb = cache.m_navigSurfs[solutions[iSol]].first->materialLayer();
808 if (mb && m_includeMaterialEffects) {
809 if (mb->layerMaterialProperties() && mb->layerMaterialProperties()->fullMaterial(nextPar->position())) {
810 const ITimedMatEffUpdator *currentUpdator = subMaterialEffectsUpdator(*cache.m_currentStatic);
811 nextPar = currentUpdator ? currentUpdator
812 ->update(nextPar,
813 *mb,
814 timeLim,
815 cache.m_path,
816 cache.m_currentStatic->geometrySignature(),
817 dir,
818 particle)
819 .release()
820 : nextPar;
821
822 if (!nextPar) {
823 ATH_MSG_VERBOSE(" [+] Update may have killed neutral track - return.");
824 cache.m_parametersAtBoundary.resetBoundaryInformation();
825 return returnParameters;
826 }
827 throwIntoGarbageBin(cache,nextPar);
828 } else { // material layer without material ?
829 ATH_MSG_VERBOSE(" boundary layer without material:" << mb->layerIndex());
830 }
831 }
832
833 // static volume boundary; return to the main loop
834 unsigned int index = solutions[iSol] - iDest;
835
836 // use global coordinates to retrieve attached volume (just for static!)
837 nextVol = (cache.m_currentStatic->boundarySurfaces())[index]->attachedVolume(
838 nextPar->position(), nextPar->momentum(), dir);
839 // double check the next volume
840 if (nextVol && !(nextVol->inside(nextPar->position() + 0.01 * dir * nextPar->momentum().normalized(), 0.))) {
841 ATH_MSG_DEBUG(
842 " [!] WARNING: wrongly assigned static volume ?" << cache.m_currentStatic->volumeName() << "->" <<
843 nextVol->volumeName());
844 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
845 nextPar->position() + 0.01 * dir * nextPar->momentum().normalized());
846 if (nextVol) {
847 ATH_MSG_DEBUG(" new search yields: " << nextVol->volumeName());
848 }
849 }
850 // end double check - to be removed after validation of the geometry gluing
851 if (nextVol != cache.m_currentStatic) {
852 cache.m_parametersAtBoundary.boundaryInformation(nextVol, nextPar, nextPar);
853 ATH_MSG_DEBUG(" [+] StaticVol boundary reached of '" << cache.m_currentStatic->volumeName() << "'.");
854 if (m_navigator->atVolumeBoundary(nextPar, cache.m_currentStatic, dir, assocVol,
855 m_tolerance) && assocVol != nextVol) {
856 cache.m_currentDense = cache.m_dense ? nextVol : cache.m_highestVolume;
857 }
858 // no next volume found --- end of the world
859 if (!nextVol) {
860 ATH_MSG_DEBUG(" [+] World boundary reached - at " << positionOutput(
861 nextPar->position()) << ", timed at " << timeLim.time);
862 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
863 if (!destVol) {
864 return nextPar->uniqueClone();
865 }
866 }
867 // next volume found and parameters are at boundary
868 if (nextVol /*&& nextPar nextPar is dereferenced anyway */) {
869 ATH_MSG_DEBUG(" [+] Crossing to next volume '" << nextVol->volumeName() << "'");
870 ATH_MSG_DEBUG(" [+] Crossing position is - at " << positionOutput(nextPar->position()));
871 if (!destVol && cache.m_currentStatic->geometrySignature() != nextVol->geometrySignature()) {
872 nextGeoID = nextVol->geometrySignature();
873 return nextPar->uniqueClone();
874 }
875 }
876 return extrapolateToVolumeWithPathLimit(cache,*nextPar, timeLim, dir, particle, nextGeoID, destVol);
877 }
878 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() + cache.m_layers.size()) {
879 // next layer; don't return passive material layers unless required
880 unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size();
881 const Trk::Layer *nextLayer = cache.m_navigLays[index].second;
882 // material update ?
883 // bool matUp = nextLayer->layerMaterialProperties() && m_includeMaterialEffects &&
884 // nextLayer->isOnLayer(nextPar->position());
885 bool matUp = nextLayer->fullUpdateMaterialProperties(*nextPar) && m_includeMaterialEffects &&
886 nextLayer->isOnLayer(nextPar->position());
887 // material update
888 const ITimedMatEffUpdator *currentUpdator = subMaterialEffectsUpdator(*cache.m_currentStatic);
889 if (matUp) {
890 double pIn = nextPar->momentum().mag();
891 nextPar = currentUpdator ? currentUpdator->update(nextPar, *nextLayer, timeLim, cache.m_path,
892 cache.m_currentStatic->geometrySignature(), dir,
893 particle).release() : nextPar;
894 if (!nextPar) {
895 ATH_MSG_VERBOSE(" [+] Update may have killed track - return.");
896 cache.m_parametersAtBoundary.resetBoundaryInformation();
897 return returnParameters;
898 }
899 ATH_MSG_VERBOSE(
900 " Layer energy loss:" << nextPar->momentum().mag() - pIn << "at position:" << nextPar->position() << ", current momentum:" <<
901 nextPar->momentum());
902 throwIntoGarbageBin(cache,nextPar);
903
904 }
905 // active surface intersections ( Fatras hits ...)
906 if (cache.m_hitVector && particle != Trk::neutron) {
907 if (nextLayer->surfaceArray()) {
908 // perform the overlap Search on this layer
909 ATH_MSG_VERBOSE(" [o] Calling overlapSearch() on layer.");
910 overlapSearch(cache,*m_subPropagators[0], *currPar, *nextPar, *nextLayer, timeLim.time, dir, true, particle);
911 } else if (nextLayer->layerType() > 0 && nextLayer->isOnLayer(nextPar->position())) {
912 ATH_MSG_VERBOSE(" [o] Collecting intersection with active layer.");
913 cache.m_hitVector->emplace_back(nextPar->uniqueClone(), timeLim.time, nextLayer->layerType(), 0.);
914 }
915 } // ------------------------------------------------- Fatras mode off -----------------------------------
916
917 // TODO : debug the retrieval of next layer
918 if (!m_robustSampling && cache.m_currentStatic->geometrySignature() != Trk::MS) {
919 if (cache.m_navigLays[index].first && cache.m_navigLays[index].first->confinedLayers()) {
920 const Trk::Layer *newLayer = nextLayer->nextLayer(nextPar->position(),
921 dir * nextPar->momentum().normalized());
922 if (newLayer && newLayer != nextLayer) {
923 bool found = false;
924 int replace = -1;
925 for (unsigned int i = 0; i < cache.m_navigLays.size(); i++) {
926 if (cache.m_navigLays[i].second == newLayer) {
927 found = true;
928 break;
929 }
930 if (cache.m_navigLays[i].second != nextLayer) {
931 replace = i;
932 }
933 }
934 if (!found) {
935 if (replace > -1) {
936 cache.m_navigLays[replace].second = newLayer;
937 cache.m_navigSurfs[solutions[iSol] + replace - index].first = &(newLayer->surfaceRepresentation());
938 } else {
939 // can't insert a surface in middle
940 return extrapolateToVolumeWithPathLimit(cache,*nextPar, timeLim, dir, particle, nextGeoID, destVol);
941 }
942 }
943 }
944 }
945 }
946 currPar = nextPar;
947 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() + cache.m_layers.size() + cache.m_denseBoundaries.size()) {
948 // dense volume boundary
949 unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size() - cache.m_layers.size();
950 std::vector< std::pair<const Trk::TrackingVolume *, unsigned int> >::iterator dIter = cache.m_denseVols.begin();
951 while (dIter != cache.m_denseVols.end() && index >= (*dIter).second) {
952 index -= (*dIter).second;
953 ++dIter;
954 }
955 if (dIter != cache.m_denseVols.end()) {
956 currVol = (*dIter).first;
957 nextVol = (currVol->boundarySurfaces())[index]->attachedVolume(*nextPar, dir);
958 // the boundary orientation is not reliable
959 Amg::Vector3D tp = nextPar->position() + 2 * m_tolerance * dir * nextPar->momentum().normalized();
960 if (!nextVol || !nextVol->inside(tp, m_tolerance)) { // search for dense volumes
961 cache.m_currentDense = cache.m_highestVolume;
962 if (cache.m_dense && cache.m_denseVols.empty()) {
963 cache.m_currentDense = cache.m_currentStatic;
964 } else {
965 for (unsigned int i = 0; i < cache.m_denseVols.size(); i++) {
966 const Trk::TrackingVolume *dVol = cache.m_denseVols[i].first;
967 if (dVol->inside(tp, m_tolerance) && dVol->zOverAtimesRho() != 0.) {
968 cache.m_currentDense = dVol;
969 ATH_MSG_DEBUG(" [+] Next dense volume found: '" << cache.m_currentDense->volumeName() << "'.");
970 break;
971 }
972 } // loop over dense volumes
973 }
974 } else {
975 cache.m_currentDense = nextVol;
976 ATH_MSG_DEBUG(" [+] Next dense volume: '" << cache.m_currentDense->volumeName() << "'.");
977 }
978 }
979 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() + cache.m_layers.size() + cache.m_denseBoundaries.size()
980 + cache.m_navigBoundaries.size()) {
981 // navig volume boundary
982 unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size() - cache.m_layers.size() -
983 cache.m_denseBoundaries.size();
984 std::vector< std::pair<const Trk::TrackingVolume *, unsigned int> >::iterator nIter = navigVols.begin();
985 while (nIter != navigVols.end() && index >= (*nIter).second) {
986 index -= (*nIter).second;
987 ++nIter;
988 }
989 if (nIter != navigVols.end()) {
990 currVol = (*nIter).first;
991 nextVol = ((*nIter).first->boundarySurfaces())[index]->attachedVolume(*nextPar, dir);
992 if (!nextVol) {
993 ATH_MSG_DEBUG(" [+] Navigation volume boundary, leaving volume '"
994 << currVol->volumeName() << "'.");
995 } else {
996 ATH_MSG_DEBUG(" [+] Navigation volume boundary, entering volume '" << nextVol->volumeName() << "'.");
997 }
998 currPar = nextPar;
999 // return only if detached volume boundaries not collected
1000 // if ( nextVol || !detachedBoundariesIncluded )
1001 if (nextVol) {
1002 return extrapolateToVolumeWithPathLimit(cache,*currPar, timeLim, dir, particle, nextGeoID, destVol);
1003 }
1004 }
1005 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() + cache.m_layers.size() + cache.m_denseBoundaries.size()
1006 + cache.m_navigBoundaries.size() + cache.m_detachedBoundaries.size()) {
1007 // detached volume boundary
1008 unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size() - cache.m_layers.size()
1009 - cache.m_denseBoundaries.size() - cache.m_navigBoundaries.size();
1010 std::vector< std::pair<const Trk::DetachedTrackingVolume *,
1011 unsigned int> >::iterator dIter = cache.m_detachedVols.begin();
1012 while (dIter != cache.m_detachedVols.end() && index >= (*dIter).second) {
1013 index -= (*dIter).second;
1014 ++dIter;
1015 }
1016 if (dIter != cache.m_detachedVols.end()) {
1017 currVol = (*dIter).first->trackingVolume();
1018 nextVol =
1019 ((*dIter).first->trackingVolume()->boundarySurfaces())[index]->attachedVolume(*nextPar, dir);
1020 if (!nextVol) {
1021 ATH_MSG_DEBUG(" [+] Detached volume boundary, leaving volume '" << currVol->volumeName() << "'.");
1022 } else {
1023 ATH_MSG_DEBUG(" [+] Detached volume boundary, entering volume '" << nextVol->volumeName() << "'.");
1024 }
1025 currPar = nextPar;
1026 // if ( nextVol || !detachedBoundariesIncluded)
1027 if (nextVol) {
1028 return extrapolateToVolumeWithPathLimit(cache, *currPar, timeLim, dir, particle, nextGeoID, destVol);
1029 }
1030 }
1031 }
1032 iSol++;
1033 }
1034 currPar = nextPar;
1035 }
1036
1037 return returnParameters;
1038}
1039
1040void
1041Trk::TimedExtrapolator::overlapSearch(Trk::TimedExtrapolator::Cache &cache,
1042 const IPropagator &prop,
1043 const TrackParameters &parm,
1044 const TrackParameters &parsOnLayer,
1045 const Layer &lay,
1046 // const TrackingVolume& tvol,
1047 float time,
1048 PropDirection dir,
1049 const BoundaryCheck& bcheck, // bcheck
1050 ParticleHypothesis particle,
1051 bool startingLayer) const {
1052
1053 const EventContext& ctx = Gaudi::Hive::currentContext();
1054 // indicate destination layer
1055 bool isDestinationLayer = false;
1056 // start and end surface for on-layer navigation
1057 // -> take the start surface if ther parameter surface is owned by detector element
1058 const Trk::Surface *startSurface = ((parm.associatedSurface()).associatedDetectorElement() && startingLayer) ?
1059 &parm.associatedSurface() : nullptr;
1060 const Trk::Surface *endSurface = nullptr;
1061 // - the best detSurface to start from is the one associated to the detector element
1062 const Trk::Surface *detSurface = (parsOnLayer.associatedSurface()).associatedDetectorElement() ?
1063 &parsOnLayer.associatedSurface() : nullptr;
1064
1065 ATH_MSG_VERBOSE(" [o] OverlapSearch called " << (startSurface ? "with " : "w/o ") << "start, "
1066 << (endSurface ? "with " : "w/o ") << "end surface.");
1067
1068 if (!detSurface) {
1069 // of parsOnLayer are different from parm, then local position is safe, because the extrapolation
1070 // to the detector surface has been done !
1071 detSurface = isDestinationLayer ? lay.subSurface(parsOnLayer.localPosition()) : lay.subSurface(
1072 parsOnLayer.position());
1073 if (detSurface) {
1074 ATH_MSG_VERBOSE(" [o] Detector surface found through subSurface() call");
1075 } else {
1076 ATH_MSG_VERBOSE(" [o] No Detector surface found on this layer.");
1077 }
1078 } else {
1079 ATH_MSG_VERBOSE(" [o] Detector surface found through parameter on layer association");
1080 }
1081
1082 // indicate the start layer
1083 bool isStartLayer = (detSurface && detSurface == startSurface);
1084
1085 const Trk::TrackParameters *detParameters = nullptr;
1086 // the temporary vector (might have to be ordered)
1087 std::vector<const Trk::TrackParameters*> detParametersOnLayer;
1088 bool reorderDetParametersOnLayer = false;
1089 // the first test for the detector surface to be hit (false test)
1090 // - only do this if the parameters aren't on the surface
1091 // (i.e. search on the start layer or end layer)
1092 if (isDestinationLayer) {
1093 detParameters = (&parsOnLayer);
1094 } else if (isStartLayer) {
1095 detParameters = (&parm);
1096 } else if (detSurface) {
1097 // detParameters = prop.propagate(parm, *detSurface, dir, false, tvol, particle);
1098 detParameters = prop.propagate(ctx,parm, *detSurface, dir, false, m_fieldProperties, particle).release();
1099 }
1100
1101 // set the surface hit to true, it is anyway overruled
1102 bool surfaceHit = true;
1103 if (detParameters &&
1104 !isStartLayer &&
1105 !isDestinationLayer) {
1106 ATH_MSG_VERBOSE(" [o] First intersection with Detector surface: " << *detParameters);
1107 // for the later use in the overlapSearch
1108 surfaceHit = detParameters && detSurface ? detSurface->isOnSurface(detParameters->position()) : 0; // ,bcheck) -
1109 // creates
1110 // problems on
1111 // start layer;
1112 // check also for start/endSurface on this level
1113
1114 surfaceHit = (surfaceHit && startSurface) ?
1115 ((detParameters->position() - parm.position()).dot(dir * parm.momentum().normalized()) >
1116 0) : surfaceHit;
1117 surfaceHit = (surfaceHit && endSurface) ?
1118 ((detParameters->position() - parsOnLayer.position()).dot(dir * parsOnLayer.momentum().normalized()) <
1119 0) : surfaceHit;
1120
1121 // surface is hit within bounds (or at least with given boundary check directive) -> it counts
1122 // surface hit also survived start/endsurface search
1123 //
1124 // Convention for Fatras: always apply the full update on the last parameters
1125 // of the gathered vector (no pre/post schema)
1126 // don't record a hit on the destination surface
1127 if (surfaceHit &&
1128 detSurface != startSurface) {
1129 ATH_MSG_VERBOSE(" [H] Hit with detector surface recorded ! ");
1130 // push into the temporary vector
1131 detParametersOnLayer.push_back(detParameters);
1132 } else if (detParameters) {
1133 // no hit -> fill into the garbage bin
1134 ATH_MSG_VERBOSE(" [-] Detector surface hit cancelled through bounds check or start/end surface check.");
1135 throwIntoGarbageBin(cache,detParameters);
1136 }
1137 }
1138
1139 // search for the overlap ------------------------------------------------------------------------
1140 if (detParameters) {
1141 // retrieve compatible subsurfaces
1142 std::vector<Trk::SurfaceIntersection> cSurfaces;
1143 size_t ncSurfaces = lay.compatibleSurfaces(cSurfaces, *detParameters, Trk::anyDirection, bcheck, false);
1144
1145 // import from StaticEngine.icc
1146 if (ncSurfaces) {
1147 ATH_MSG_VERBOSE("found " << ncSurfaces << " candidate sensitive surfaces to test.");
1148 // now loop over the surfaces:
1149 // the surfaces will be sorted @TODO integrate pathLength propagation into this
1150 for (auto &csf : cSurfaces) {
1151 // propagate to the compatible surface, return types are (pathLimit failure is excluded by Trk::anyDirection for
1152 // the moment):
1153 const Trk::TrackParameters *overlapParameters = prop.propagate(ctx,
1154 parm,
1155 *(csf.object),
1156 Trk::anyDirection,
1157 true,
1158 m_fieldProperties,
1159 particle).release();
1160
1161 if (overlapParameters) {
1162 ATH_MSG_VERBOSE(" [+] Overlap surface was hit, checking start/end surface condition.");
1163 // check on start / end surface for on-layer navigaiton action
1164 surfaceHit = (startSurface) ?
1165 ((overlapParameters->position() - parm.position()).dot(dir * parm.momentum().normalized()) >
1166 0) : true;
1167 surfaceHit = (surfaceHit && endSurface) ?
1168 ((overlapParameters->position() - parsOnLayer.position()).dot(dir *
1169 parsOnLayer.momentum().normalized())
1170 < 0) : surfaceHit;
1171 if (surfaceHit && csf.object!=detSurface) { //skipping the initial surface on which a hit has already been created
1172 ATH_MSG_VERBOSE(" [H] Hit with detector surface recorded !");
1173 // distinguish whether sorting is needed or not
1174 reorderDetParametersOnLayer = true;
1175 // push back into the temporary vector
1176 detParametersOnLayer.push_back(overlapParameters);
1177 } else { // the parameters have been cancelled by start/end surface
1178 // no hit -> fill into the garbage bin
1179 ATH_MSG_VERBOSE(" [-] Detector surface hit cancelled through start/end surface check.");
1180 throwIntoGarbageBin(cache,overlapParameters);
1181 }
1182 }
1183 } // loop over test surfaces done
1184 } // there are compatible surfaces
1185 } // ---------------------------------------------------------------------------------------------
1186
1187 // push them into the parameters vector
1188 std::vector<const Trk::TrackParameters *>::const_iterator parsOnLayerIter = detParametersOnLayer.begin();
1189 std::vector<const Trk::TrackParameters *>::const_iterator parsOnLayerIterEnd = detParametersOnLayer.end();
1190
1191 // reorder the track parameters if neccessary, the overlap descriptor did not provide the ordered surfaces
1192 if (reorderDetParametersOnLayer) {
1193 // sort to reference of incoming parameters
1194 Trk::TrkParametersComparisonFunction parameterSorter(parm.position());
1195 sort(detParametersOnLayer.begin(), detParametersOnLayer.end(), parameterSorter);
1196 }
1197
1198 // after sorting : reset the iterators
1199 parsOnLayerIter = detParametersOnLayer.begin();
1200 parsOnLayerIterEnd = detParametersOnLayer.end();
1201 // now fill them into the parameter vector -------> hit creation done <----------------------
1202 for (; parsOnLayerIter != parsOnLayerIterEnd; ++parsOnLayerIter) {
1203 if (cache.m_hitVector) {
1204 cache.m_hitVector->emplace_back(
1205 std::unique_ptr<const Trk::TrackParameters>(*parsOnLayerIter),
1206 time,
1207 0,
1208 0.);
1209 }
1210 }
1211}
1212
1213std::string
1214Trk::TimedExtrapolator::positionOutput(const Amg::Vector3D &pos) const {
1215 std::stringstream outStream;
1216
1217 if (m_printRzOutput) {
1218 outStream << "[r,phi,z] = [ " << pos.perp() << ", " << pos.phi() << ", " << pos.z() << " ]";
1219 } else {
1220 outStream << "[xyz] = [ " << pos.x() << ", " << pos.y() << ", " << pos.z() << " ]";
1221 }
1222 return outStream.str();
1223}
1224
1225std::string
1226Trk::TimedExtrapolator::momentumOutput(const Amg::Vector3D &mom) {
1227 std::stringstream outStream;
1228
1229 outStream << "[eta,phi] = [ " << mom.eta() << ", " << mom.phi() << " ]";
1230 return outStream.str();
1231}
1232
1233void
1234Trk::TimedExtrapolator::emptyGarbageBin(Trk::TimedExtrapolator::Cache &cache,
1235 const Trk::TrackParameters *trPar) const {
1236 // empty the garbage
1237 std::map<const Trk::TrackParameters *, bool>::iterator garbageIter = cache.m_garbageBin.begin();
1238 std::map<const Trk::TrackParameters *, bool>::iterator garbageEnd = cache.m_garbageBin.end();
1239
1240 bool throwCurrent = false;
1241
1242 for (; garbageIter != garbageEnd; ++garbageIter) {
1243 if (garbageIter->first && garbageIter->first != trPar) {
1244 delete (garbageIter->first);
1245 }
1246 if (garbageIter->first && garbageIter->first == trPar) {
1247 throwCurrent = true;
1248 }
1249 }
1250
1251 cache.m_garbageBin.clear();
1252 if (throwCurrent) {
1253 throwIntoGarbageBin(cache,trPar);
1254 }
1255}
1256
1257// the validation action -> propagated to the SubTools
1258void
1259Trk::TimedExtrapolator::validationAction() const {
1260 // record the updator validation information
1261 for (const auto *subUpdator : m_subUpdators) {
1262 subUpdator->validationAction();
1263 }
1264 // record the navigator validation information
1265}
1266
1267std::unique_ptr<const Trk::TrackParameters>
1268Trk::TimedExtrapolator::transportNeutralsWithPathLimit(const Trk::TrackParameters &parm,
1269 Trk::PathLimit &pathLim, Trk::TimeLimit &timeLim,
1270 Trk::PropDirection dir,
1271 Trk::ParticleHypothesis particle,
1272 std::vector<Trk::HitInfo> * &hitInfo,
1273 Trk::GeometrySignature &nextGeoID,
1274 const Trk::TrackingVolume *boundaryVol) const {
1275 Trk::TimedExtrapolator::Cache cache(m_maxNavigSurf);
1276// extrapolation method intended for simulation of particle decay; collects the material up to pre-defined limit and
1277// triggers
1278// material interaction
1279// possible outcomes:1/ returns curvilinear parameters after reaching the maximal path (if to be destroyed)
1280// 2/ returns parameters at destination volume boundary
1281// 3/ returns 0 ( particle stopped ) but keeps material and timing info
1282
1283 ATH_MSG_DEBUG(
1284 "M-[" << ++cache.m_methodSequence << "] transportNeutralsWithPathLimit(...) " << pathLim.x0Max << ", from " <<
1285 parm.position());
1286
1287 // reset the path ( in x0 !!)
1288 cache.m_path = PathLimit(pathLim.x0Max - pathLim.x0Collected, pathLim.process); // collect material locally
1289
1290 // initialize time info
1291 cache.m_time = timeLim.time;
1292
1293 // initialize hit vector
1294 cache.m_hitVector = hitInfo;
1295
1296 cache.m_parametersAtBoundary.resetBoundaryInformation();
1297
1298 // if no input volume, define as highest volume
1299 // const Trk::TrackingVolume* destVolume = boundaryVol ? boundaryVol : m_navigator->highestVolume();
1300 cache.m_currentStatic = nullptr;
1301 if (boundaryVol && !boundaryVol->inside(parm.position(), m_tolerance)) {
1302 return nullptr;
1303 }
1304
1305 cache.m_particleMass = Trk::ParticleMasses::mass[particle];
1306
1307 // extrapolate to destination volume boundary with path limit
1308 std::unique_ptr<const Trk::TrackParameters> returnParms =
1309 transportToVolumeWithPathLimit(
1310 cache, parm, timeLim, dir, particle, nextGeoID, boundaryVol);
1311
1312 // save actual path on output
1313 if (cache.m_path.x0Collected > 0.) {
1314 pathLim.updateMat(cache.m_path.x0Collected, cache.m_path.weightedZ / cache.m_path.x0Collected, cache.m_path.l0Collected);
1315 }
1316
1317 // return timing
1318 timeLim.time = cache.m_time;
1319
1320 std::map<const Trk::TrackParameters *, bool>::iterator garbageIter = cache.m_garbageBin.begin();
1321 std::map<const Trk::TrackParameters *, bool>::iterator garbageEnd = cache.m_garbageBin.end();
1322 for (; garbageIter != garbageEnd; ++garbageIter) if (garbageIter->first) {
1323 if(garbageIter->first == returnParms.get()) {
1324 auto ret=returnParms->uniqueClone();
1325 ATH_MSG_DEBUG(" [+] garbage - at "
1326 << positionOutput(garbageIter->first->position())
1327 << " parm=" << garbageIter->first
1328 << " is the return param. Cloning to" << ret.get());
1329 returnParms=std::move(ret);
1330 }
1331 }
1332
1333 return returnParms;
1334}
1335
1336std::unique_ptr<const Trk::TrackParameters>
1337Trk::TimedExtrapolator::transportToVolumeWithPathLimit(
1338 Trk::TimedExtrapolator::Cache& cache,
1339 const Trk::TrackParameters& parm,
1340 Trk::TimeLimit& timeLim,
1341 Trk::PropDirection dir,
1342 Trk::ParticleHypothesis particle,
1343 Trk::GeometrySignature& nextGeoID,
1344 const Trk::TrackingVolume* destVol) const
1345{
1346 // returns:
1347 // A) curvilinear track parameters if path or time limit reached
1348 // B) boundary parameters (at destination volume boundary)
1349
1350 // initialize the return parameters vector
1351 std::unique_ptr<const Trk::TrackParameters> returnParameters = nullptr;
1352 const Trk::TrackParameters *currPar = &parm;
1353 const Trk::TrackingVolume *currVol = nullptr;
1354 const Trk::TrackingVolume *nextVol = nullptr;
1355 const Trk::TrackingVolume *assocVol = nullptr;
1356 // int nEntryLays = 0;
1357 unsigned int iDest = 0;
1358
1359 const EventContext& ctx = Gaudi::Hive::currentContext();
1360 // destination volume boundary ?
1361 if (destVol && m_navigator->atVolumeBoundary(currPar, destVol, dir, nextVol, m_tolerance) && nextVol != destVol) {
1362 return parm.uniqueClone();
1363 }
1364
1365 // bool resolveActive = m_resolveActive;
1366 if (!cache.m_highestVolume) {
1367 cache.m_highestVolume = m_navigator->highestVolume(ctx);
1368 }
1369
1370 emptyGarbageBin(cache,&parm);
1371 // transport surfaces: collect only those with valid intersection (easy to calculate for neutrals)
1372 if (cache.m_trSurfs.capacity() > m_maxNavigSurf) {
1373 cache.m_trSurfs.reserve(m_maxNavigSurf);
1374 }
1375 cache.m_trSurfs.clear();
1376
1377 // target volume may not be part of tracking geometry
1378 if (destVol) {
1379 const Trk::TrackingVolume *tgVol = m_navigator->trackingGeometry(ctx)->trackingVolume(destVol->volumeName());
1380 if (!tgVol || tgVol != destVol) {
1381 const auto& bounds = destVol->boundarySurfaces();
1382 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
1383 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
1384 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1385 dir * currPar->momentum().normalized());
1386 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1387 // boundary check
1388 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1389 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1390 iDest++;
1391 cache.m_trSurfs.emplace_back(&surf, distSol.first());
1392 } // valid intersection
1393 } // along path
1394 if (distSol.numberOfSolutions() > 1 && distSol.second() > 0.) {
1395 // boundary check
1396 Amg::Vector3D gp = currPar->position() + distSol.second() * dir * currPar->momentum().normalized();
1397 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1398 iDest++;
1399 cache.m_trSurfs.emplace_back(&surf, distSol.second());
1400 } // valid intersection
1401 }
1402 } // end loop over boundaries
1403 } // end process external volume
1404 }
1405
1406 // resolve current position
1407 if (cache.m_parametersAtBoundary.nextParameters == currPar) {
1408 cache.m_currentStatic = cache.m_parametersAtBoundary.nextVolume;
1409 } else {
1410 const Amg::Vector3D& gp = parm.position();
1411 if (!cache.m_currentStatic || !cache.m_currentStatic->inside(gp, m_tolerance)) {
1412 cache.m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
1413
1414 if (!cache.m_currentStatic ||
1415 !cache.m_currentStatic->inside(currPar->position() + 0.01 * dir * currPar->momentum().normalized(), 0.)) {
1416 cache.m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(currPar->position()
1417 + 0.01 * dir *
1418 currPar->momentum().normalized());
1419 }
1420 }
1421
1422 if (!cache.m_currentStatic) {
1423 // no next volume found --- end of the world
1424 ATH_MSG_DEBUG(" [+] Word boundary reached - at " << positionOutput(currPar->position()));
1425 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
1426 return currPar->uniqueClone();
1427 }
1428 }
1429
1430 // current frame volume known-retrieve geoID
1431 nextGeoID = cache.m_currentStatic->geometrySignature();
1432
1433 // resolve active Calo volumes if hit info required
1434 if (cache.m_hitVector && nextGeoID == Trk::Calo) {
1435 const Trk::AlignableTrackingVolume *alignTV = dynamic_cast<const Trk::AlignableTrackingVolume *> (cache.m_currentStatic);
1436 if (alignTV) {
1437 const Trk::TrackParameters *aPar = transportInAlignableTV(cache,parm, timeLim, dir, particle, nextGeoID, alignTV).trPar;
1438 if (!aPar) {
1439 return returnParameters;
1440 }
1441 throwIntoGarbageBin(cache,aPar);
1442 return transportToVolumeWithPathLimit(cache,*aPar, timeLim, dir, particle, nextGeoID, destVol);
1443 }
1444 }
1445
1446 // distance to static volume boundaries recalculated
1447 // retrieve boundaries along path
1448 cache.m_trStaticBounds.clear();
1449 const auto& bounds = cache.m_currentStatic->boundarySurfaces();
1450 for (unsigned int ib = 0; ib < bounds.size(); ++ib) {
1451 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
1452 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1453 dir * currPar->momentum().normalized());
1454 if (distSol.numberOfSolutions() > 0 &&
1455 (distSol.currentDistance(false) > m_tolerance || distSol.numberOfSolutions() > 1) &&
1456 distSol.first() > m_tolerance) {
1457 double dist = distSol.first();
1458 // resolve multiple intersection solutions
1459 if (distSol.numberOfSolutions() > 1 && dist < m_tolerance && distSol.second() > dist) {
1460 dist = distSol.second();
1461 }
1462 // boundary check
1463 Amg::Vector3D gp = currPar->position() + dist * dir * currPar->momentum().normalized();
1464 if (surf.isOnSurface(gp, true, m_tolerance, m_tolerance)) {
1465 cache.m_trStaticBounds.insert(cache.m_trStaticBounds.begin(), Trk::DestBound(&surf, dist, ib));
1466 }
1467 } // along path
1468 if (distSol.numberOfSolutions() > 1 && distSol.second() > m_tolerance) {
1469 double dist = distSol.second();
1470 // boundary check
1471 Amg::Vector3D gp = currPar->position() + dist * dir * currPar->momentum().unit();
1472 if (surf.isOnSurface(gp, true, m_tolerance, m_tolerance)) {
1473 if (dist > m_tolerance) { // valid intersection
1474 cache.m_trStaticBounds.insert(cache.m_trStaticBounds.begin(), Trk::DestBound(&surf, dist, ib));
1475 }
1476 }
1477 } // along path
1478 } // end loop over boundaries
1479
1480 if (cache.m_trStaticBounds.empty()) {
1481 ATH_MSG_WARNING(
1482 " transportToVolumeWithPathLimit() - at " << currPar->position() << ", missing static volume boundary "
1483 << cache.m_currentStatic->volumeName() <<
1484 ": transport interrupted");
1485
1486 ATH_MSG_DEBUG(
1487 "---> particle R,phi,z, momentum:" << currPar->position().perp() << "," << currPar->position().phi() << "," << currPar->position().z() << "," <<
1488 currPar->momentum());
1489 ATH_MSG_DEBUG("---> static volume position:" << cache.m_currentStatic->center());
1490 const Trk::CylinderVolumeBounds *cyl =
1491 dynamic_cast<const Trk::CylinderVolumeBounds *> (&(cache.m_currentStatic->volumeBounds()));
1492 if (cyl) {
1493 ATH_MSG_DEBUG(
1494 "---> cylinder volume dimensions:" << cyl->innerRadius() << "," << cyl->outerRadius() << "," <<
1495 cyl->halflengthZ());
1496 }
1497
1498
1499 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
1500 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
1501 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1502 dir * currPar->momentum().unit());
1503 ATH_MSG_DEBUG(
1504 "---> decomposed boundary surface position, normal, estimated distance:" << ib << "," << surf.center() << "," <<
1505 surf.normal());
1506 ATH_MSG_DEBUG(
1507 "---> estimated distance to (first solution):boundary check:" << distSol.numberOfSolutions() << "," << distSol.first() << ":" <<
1508 surf.isOnSurface(currPar->position() + distSol.first() * dir * currPar->momentum().unit(), true,
1509 m_tolerance, m_tolerance));
1510 if (distSol.numberOfSolutions() > 1) {
1511 ATH_MSG_DEBUG("---> estimated distance to (second solution):boundary check:" << distSol.second() << "," <<
1512 surf.isOnSurface(currPar->position() + distSol.second() * dir * currPar->momentum().unit(), true,
1513 m_tolerance, m_tolerance));
1514 }
1515 }
1516
1517 return returnParameters;
1518 } if (cache.m_trStaticBounds[0].distance < m_tolerance) {
1519 // TODO find out why this case (=exit from volume) haven't been handled by Navigator
1520 // ATH_MSG_WARNING( " recovering from glitch at the static volume boundary:"<<cache.m_trStaticBounds[0].distance );
1521
1522 Amg::Vector3D gp = currPar->position() + m_tolerance * dir * currPar->momentum().unit();
1523 cache.m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
1524
1525 if (cache.m_currentStatic) {
1526 return transportToVolumeWithPathLimit(cache,parm, timeLim, dir, particle, nextGeoID, destVol);
1527 }
1528 ATH_MSG_DEBUG(" [+] World boundary reached - at " << positionOutput(
1529 currPar->position()) << ", timed at " << cache.m_time);
1530 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
1531 // if (!destVol) { return currPar;}
1532 return currPar->uniqueClone();
1533
1534 }
1535
1536 cache.m_detachedVols.clear();
1537 cache.m_denseVols.clear();
1538 cache.m_trDenseBounds.clear();
1539 cache.m_trLays.clear();
1540 cache.m_navigLays.clear();
1541
1542 // detached volume boundaries
1543 Trk::ArraySpan<const Trk::DetachedTrackingVolume* const> detVols = cache.m_currentStatic->confinedDetachedVolumes();
1544 if (!detVols.empty()) {
1545 Trk::ArraySpan<const Trk::DetachedTrackingVolume* const>::iterator iTer = detVols.begin();
1546 for (; iTer != detVols.end(); ++iTer) {
1547 // active station ?
1548 const Trk::Layer *layR = (*iTer)->layerRepresentation();
1549 bool active = layR && layR->layerType();
1550
1551 if (active) {
1552 if (!m_resolveMultilayers || (*iTer)->multilayerRepresentation().empty()) {
1553 const Trk::Surface &surf = layR->surfaceRepresentation();
1554 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1555 dir * currPar->momentum().normalized());
1556 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1557 // boundary check
1558 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1559 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1560 cache.m_trLays.emplace_back(&surf, distSol.first());
1561 cache.m_navigLays.emplace_back((*iTer)->trackingVolume(), layR);
1562 }
1563 }
1564 } else {
1565 const auto& multi = (*iTer)->multilayerRepresentation();
1566 for (const auto *i : multi) {
1567 const Trk::Surface &surf = i->surfaceRepresentation();
1568 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1569 dir * currPar->momentum().normalized());
1570 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1571 // boundary check
1572 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1573 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1574 cache.m_trLays.emplace_back(&surf, distSol.first());
1575 cache.m_navigLays.emplace_back((*iTer)->trackingVolume(), i);
1576 }
1577 }
1578 } // end loop over multilayers
1579 } // end unresolved active
1580 } // active done
1581 else if (cache.m_currentStatic->geometrySignature() != Trk::MS || !m_useMuonMatApprox ||
1582 (*iTer)->name().substr((*iTer)->name().size() - 4, 4) == "PERM") { // retrieve inert detached objects
1583 // only if needed
1584 // dense volume boundaries
1585 if ((*iTer)->trackingVolume()->zOverAtimesRho() != 0. &&
1586 ((*iTer)->trackingVolume()->confinedDenseVolumes().empty())
1587 && ((*iTer)->trackingVolume()->confinedArbitraryLayers().empty())) {
1588 const auto& detBounds = (*iTer)->trackingVolume()->boundarySurfaces();
1589 int newB = 0;
1590 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
1591 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
1592 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1593 dir * currPar->momentum().normalized());
1594 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1595 // boundary check
1596 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1597 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1598 cache.m_trDenseBounds.emplace_back(&surf, distSol.first());
1599 newB++;
1600 } // valid intersection
1601 } // along path
1602 } // end loop over boundaries
1603 if (newB > 0) {
1604 cache.m_denseVols.emplace_back((*iTer)->trackingVolume(), newB);
1605 }
1606 }
1607 // subvolumes ?
1608 // if ((*iTer)->trackingVolume()->confinedDenseVolumes() &&
1609 // (*iTer)->trackingVolume()->confinedDenseVolumes()->size())
1610 // ATH_MSG_WARNING( " transportToVolumeWithPathLimit() - at " << currPar->position() <<", unresolved
1611 // subvolumes for "
1612 // << (*iTer)->trackingVolume()->volumeName() );
1613
1614 const auto confinedDense =
1615 (*iTer)->trackingVolume()->confinedDenseVolumes();
1616 if (!confinedDense.empty()) {
1617 auto vIter = confinedDense.begin();
1618 for (; vIter != confinedDense.end(); ++vIter) {
1619 const auto& bounds = (*vIter)->boundarySurfaces();
1620 int newB = 0;
1621 for (unsigned int ibb = 0; ibb < bounds.size(); ibb++) {
1622 const Trk::Surface &surf = (bounds[ibb])->surfaceRepresentation();
1623 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1624 dir * currPar->momentum().normalized());
1625 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1626 // boundary check
1627 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1628 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1629 cache.m_trDenseBounds.emplace_back(&surf, distSol.first());
1630 newB++;
1631 } // valid intersection
1632 } // along path
1633 } // end loop over boundaries
1634 if (newB > 0) {
1635 cache.m_denseVols.emplace_back((*vIter), newB);
1636 }
1637 if (!(*vIter)->confinedArbitraryLayers().empty()) {
1638 ATH_MSG_DEBUG(
1639 " transportToVolumeWithPathLimit() - at " << currPar->position() << ", unresolved sublayers/subvolumes for "
1640 << (*vIter)->volumeName());
1641 }
1642 }
1643 }
1644
1645 // confined layers
1646 Trk::ArraySpan<const Trk::Layer* const>confLays = (*iTer)->trackingVolume()->confinedArbitraryLayers();
1647 if (!confLays.empty()) {
1648 for (const Trk::Layer* const lIt: confLays) {
1649 const Trk::Surface &surf = lIt->surfaceRepresentation();
1650 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1651 dir * currPar->momentum().normalized());
1652 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1653 // boundary check
1654 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1655 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1656 cache.m_trLays.emplace_back(&surf, distSol.first());
1657 cache.m_navigLays.emplace_back((*iTer)->trackingVolume(), lIt);
1658 } // valid intersection
1659 } // along path
1660 }
1661 } // end confined layers
1662 } // end inert material
1663 }
1664 } // end detached volumes
1665 cache.m_denseResolved = std::pair<unsigned int, unsigned int> (cache.m_denseVols.size(), cache.m_trDenseBounds.size());
1666 cache.m_layerResolved = cache.m_trLays.size();
1667
1668 std::vector< Trk::DestBound >::iterator bIter = cache.m_trStaticBounds.begin();
1669 while (bIter != cache.m_trStaticBounds.end()) {
1670 cache.m_trSurfs.emplace_back((*bIter).surface, (*bIter).distance);
1671 ++bIter;
1672 }
1673
1674 // std::cout <<"navigation in current static:"<< cache.m_trSurfs.size()<<","<<cache.m_trStaticBounds.size()<< std::endl;
1675 // for (unsigned int ib=0; ib<cache.m_trSurfs.size(); ib++) std::cout <<"distance to static:"<<
1676 // ib<<","<<cache.m_trSurfs[ib].second<<std::endl;
1677
1678 // resolve the use of dense volumes
1679 cache.m_dense = (cache.m_currentStatic->geometrySignature() == Trk::MS && m_useMuonMatApprox) ||
1680 (cache.m_currentStatic->geometrySignature() != Trk::MS && m_useDenseVolumeDescription);
1681
1682 // reset remaining counters
1683 cache.m_currentDense = cache.m_dense ? cache.m_currentStatic : cache.m_highestVolume;
1684 cache.m_navigBoundaries.clear();
1685 if (cache.m_denseVols.size() > cache.m_denseResolved.first) {
1686 cache.m_denseVols.resize(cache.m_denseResolved.first);
1687 cache.m_trDenseBounds.resize(cache.m_denseResolved.second);
1688 }
1689 if (cache.m_layers.size() > cache.m_layerResolved) {
1690 cache.m_trLays.resize(cache.m_layerResolved);
1691 cache.m_navigLays.resize(cache.m_layerResolved);
1692 }
1693
1694 // if (cache.m_currentStatic->entryLayerProvider()) nEntryLays = cache.m_currentStatic->entryLayerProvider()->layers().size();
1695
1696 // confined layers
1697 if (cache.m_currentStatic->confinedLayers()) {
1698 std::span <Trk::Layer const * const> cLays = cache.m_currentStatic->confinedLayers()->arrayObjects();
1699 for (const auto *cLay : cLays) {
1700 if (cLay->layerMaterialProperties()) {
1701 const Trk::Surface &surf = cLay->surfaceRepresentation();
1702 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
1703 dir * currPar->momentum().normalized());
1704 if (distSol.numberOfSolutions() > 0 && distSol.first() > 0.) {
1705 // boundary check
1706 Amg::Vector3D gp = currPar->position() + distSol.first() * dir * currPar->momentum().normalized();
1707 if (surf.isOnSurface(gp, true, 0.001, 0.001)) {
1708 cache.m_trLays.emplace_back(&surf, distSol.first());
1709 cache.m_navigLays.emplace_back(cache.m_currentStatic,
1710 cLay);
1711 } // valid intersection
1712 } // along path
1713 }
1714 }
1715 }
1716
1717 // cache.m_trSurfs contains destination surface (if it exists), static volume boundaries
1718 // complete with TG cache.m_layers/dynamic layers, cache.m_denseBoundaries, cache.m_navigBoundaries, m_detachedBoundaries
1719
1720 if (!cache.m_trLays.empty()) {
1721 cache.m_trSurfs.insert(cache.m_trSurfs.end(), cache.m_trLays.begin(), cache.m_trLays.end());
1722 }
1723 if (!cache.m_trDenseBounds.empty()) {
1724 cache.m_trSurfs.insert(cache.m_trSurfs.end(), cache.m_trDenseBounds.begin(), cache.m_trDenseBounds.end());
1725 }
1726
1727 // current dense
1728 cache.m_currentDense = cache.m_highestVolume;
1729
1730 for (unsigned int i = 0; i < cache.m_denseVols.size(); i++) {
1731 const Trk::TrackingVolume *dVol = cache.m_denseVols[i].first;
1732 if (dVol->inside(currPar->position(), m_tolerance) && dVol->zOverAtimesRho() != 0.) {
1733 if (!m_navigator->atVolumeBoundary(currPar, dVol, dir, nextVol, m_tolerance) ||
1734 dVol->inside(currPar->position() + 2 * m_tolerance * currPar->momentum().unit(), m_tolerance)) {
1735 cache.m_currentDense = dVol;
1736 }
1737 }
1738 }
1739
1740 if (cache.m_dense && cache.m_currentDense == cache.m_highestVolume) {
1741 cache.m_currentDense = cache.m_currentStatic;
1742 }
1743
1744 // ready to process
1745 // 1/ order valid intersections ( already in trSurfs )
1746
1747 std::vector<unsigned int> sols;
1748 sols.reserve(cache.m_trSurfs.size());
1749 for (unsigned int i = 0; i < cache.m_trSurfs.size(); ++i) {
1750 sols.push_back(i);
1751 }
1752
1753 if (sols.size() > 1) {
1754 unsigned int itest = 1;
1755 while (itest < sols.size()) {
1756 if (cache.m_trSurfs[sols[itest]].second < cache.m_trSurfs[sols[itest - 1]].second) {
1757 unsigned int iex = sols[itest - 1];
1758 sols[itest - 1] = sols[itest];
1759 sols[itest] = iex;
1760 itest = 1;
1761 } else {
1762 itest++;
1763 }
1764 }
1765 // check ordering
1766 for (unsigned int is = 1; is < sols.size(); is++) {
1767 if (cache.m_trSurfs[sols[is]].second < cache.m_trSurfs[sols[is - 1]].second) {
1768 std::cout << "wrong intersection ordering" << std::endl;
1769 }
1770 }
1771 }
1772
1773
1774 // 2/ check time/material/boundary limit
1775
1776 // update of cache.m_navigSurfs required if I/ entry into new navig volume, II/ exit from currentActive without overlaps
1777
1778 nextVol = nullptr;
1779 const Trk::TrackParameters *nextPar = nullptr;
1780
1781 double dist = 0.;
1782 double mom = currPar->momentum().mag();
1783 double beta = mom / sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass) * Gaudi::Units::c_light;
1784
1785 ATH_MSG_DEBUG(" [0] starting transport of neutral particle in (dense) volume " << cache.m_currentDense->volumeName());
1786
1787 for (unsigned int sol : sols) {
1788 if (cache.m_trSurfs[sol].second == 0.) {
1789 continue;
1790 }
1791
1792 double step = cache.m_trSurfs[sol].second - dist;
1793
1794 Amg::Vector3D nextPos = currPar->position() + dir * currPar->momentum().normalized() * cache.m_trSurfs[sol].second;
1795 // Amg::Vector3D halfStep = nextPos - 0.5*step*dir*currPar->momentum().normalized();
1796
1797 // check missing volume boundary
1798 if (!(cache.m_currentDense->inside(nextPos, m_tolerance))) {
1799 ATH_MSG_DEBUG(" [!] WARNING: missing volume boundary for volume" << cache.m_currentDense->volumeName());
1800 // new search
1801 cache.m_currentDense = cache.m_highestVolume;
1802 for (unsigned int i = 0; i < cache.m_denseVols.size(); i++) {
1803 const Trk::TrackingVolume *dVol = cache.m_denseVols[i].first;
1804 if (dVol->inside(nextPos, m_tolerance) && dVol->zOverAtimesRho() != 0.) {
1805 cache.m_currentDense = dVol;
1806 }
1807 }
1808 if (cache.m_dense && cache.m_currentDense == cache.m_highestVolume) {
1809 cache.m_currentDense = cache.m_currentStatic;
1810 }
1811
1812 ATH_MSG_DEBUG(" [!] new search for dense volume : " << cache.m_currentDense->volumeName());
1813 }
1814
1815 double tDelta = step / beta;
1816
1817 double mDelta = (cache.m_currentDense->zOverAtimesRho() != 0.) ? step / cache.m_currentDense->x0() : 0.;
1818
1819 // in case of hadronic interaction retrieve nuclear interaction properties, too
1820
1821 double frT = 1.;
1822 if (step > 0 && timeLim.tMax > cache.m_time && cache.m_time + tDelta >= timeLim.tMax) {
1823 frT = (timeLim.tMax - cache.m_time) * beta / step;
1824 }
1825
1826 // TODO : compare x0 or l0 according to the process type
1827 double frM = 1.;
1828 if (mDelta > 0 && cache.m_path.x0Max > 0.) {
1829 if (cache.m_path.process < 100 && cache.m_path.x0Collected + mDelta > cache.m_path.x0Max) {
1830 frM = (cache.m_path.x0Max - cache.m_path.x0Collected) / mDelta;
1831 } else { // waiting for hadronic interaction, retrieve nuclear interaction properties
1832 double mDeltaL = cache.m_currentDense->L0 >
1833 0. ? step / cache.m_currentDense->L0 : mDelta / 0.37 / cache.m_currentDense->averageZ();
1834 if (cache.m_path.l0Collected + mDeltaL > cache.m_path.x0Max) {
1835 frM = (cache.m_path.x0Max - cache.m_path.l0Collected) / mDeltaL;
1836 }
1837 }
1838 }
1839
1840 double fr = fmin(frT, frM);
1841
1842 // std::cout << "looping over intersections:"<<is<<","<< cache.m_trSurfs[sols[is]].second<<","<<step << ","<<
1843 // tDelta<<","<<mDelta << std::endl;
1844
1845 if (fr < 1.) { // decay or material interaction during the step
1846 int process = frT < frM ? timeLim.process : cache.m_path.process;
1847 cache.m_time += fr * step / beta;
1848 if (mDelta > 0 && cache.m_currentDense->averageZ() > 0) {
1849 cache.m_path.updateMat(fr * mDelta, cache.m_currentDense->averageZ(), 0.);
1850 }
1851
1852 nextPos = currPar->position() + dir * currPar->momentum().normalized() * (dist + fr * step);
1853
1854 // process interaction only if creation of secondaries allowed
1855 if (cache.m_currentStatic->geometrySignature() == Trk::ID || m_caloMsSecondary) {
1856 const Trk::TrackParameters* nextPar =
1857 m_updators[0]
1858 ->interact(cache.m_time, nextPos, currPar->momentum(), particle, process, cache.m_currentDense)
1859 .release();
1860
1861 if (nextPar) {
1862 ATH_MSG_DEBUG(" [!] WARNING: particle survives the interaction " << process);
1863 }
1864
1865 if (nextPar && process == 121) {
1866 ATH_MSG_DEBUG(" [!] WARNING: failed hadronic interaction, killing the input particle anyway");
1867 return returnParameters;
1868 }
1869
1870 if (!nextPar) {
1871 return returnParameters;
1872 }
1873
1874 throwIntoGarbageBin(cache,nextPar);
1875 // return transportToVolumeWithPathLimit(cache,*nextPar, timeLim, dir, particle, nextGeoID, destVol);
1876 } else { // kill particle without trace
1877 return returnParameters;
1878 }
1879 } // end decay or material interaction durign the step
1880
1881 // update
1882 dist = cache.m_trSurfs[sol].second;
1883 if (mDelta > 0 && cache.m_currentDense->averageZ() > 0) {
1884 cache.m_path.updateMat(mDelta, cache.m_currentDense->averageZ(), 0.);
1885 }
1886 cache.m_time += tDelta;
1887
1888 nextPar = new Trk::CurvilinearParameters(nextPos, currPar->momentum(), 1.); // fake charge
1889 throwIntoGarbageBin(cache,nextPar);
1890
1891 if (sol < iDest) { // destination volume (most often, subdetector boundary)
1892 return nextPar->uniqueClone();
1893 } if (sol < iDest + cache.m_trStaticBounds.size()) { // tracking geometry frame
1894 // material attached ?
1895 const Trk::Layer *mb = cache.m_trStaticBounds[sol - iDest].surface->materialLayer();
1896 if (mb && m_includeMaterialEffects) {
1897 if (mb->layerMaterialProperties() && mb->layerMaterialProperties()->fullMaterial(nextPos)) {
1898 const ITimedMatEffUpdator *currentUpdator = subMaterialEffectsUpdator(*cache.m_currentStatic);
1899 nextPar =
1900 currentUpdator
1901 ? currentUpdator
1902 ->update(
1903 nextPar, *mb, timeLim, cache.m_path, cache.m_currentStatic->geometrySignature(), dir, particle)
1904 .release()
1905 : nextPar;
1906
1907 if (!nextPar) {
1908 ATH_MSG_VERBOSE(" [+] Update may have killed neutral track - return.");
1909 cache.m_parametersAtBoundary.resetBoundaryInformation();
1910 return returnParameters;
1911 }
1912 throwIntoGarbageBin(cache,nextPar);
1913
1914 } else { // material layer without material ?
1915 ATH_MSG_VERBOSE(" boundary layer without material:" << mb->layerIndex());
1916 }
1917 }
1918
1919 // static volume boundary; return to the main loop
1920 unsigned int index = cache.m_trStaticBounds[sol - iDest].bIndex;
1921 // use global coordinates to retrieve attached volume (just for static!)
1922 nextVol = (cache.m_currentStatic->boundarySurfaces())[index]->attachedVolume(
1923 nextPar->position(), nextPar->momentum(), dir);
1924 // double check the next volume
1925 if (nextVol && !(nextVol->inside(nextPar->position() + 0.01 * dir * nextPar->momentum().normalized(), 0.))) {
1926 ATH_MSG_DEBUG(
1927 " [!] WARNING: wrongly assigned static volume ?" << cache.m_currentStatic->volumeName() << "->" <<
1928 nextVol->volumeName());
1929 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
1930 nextPar->position() + 0.01 * dir * nextPar->momentum().normalized());
1931 if (nextVol) {
1932 ATH_MSG_DEBUG(" new search yields: " << nextVol->volumeName());
1933 }
1934 }
1935 // end double check - to be removed after validation of the geometry gluing
1936 if (nextVol != cache.m_currentStatic) {
1937 ATH_MSG_DEBUG(" [+] StaticVol boundary reached of '" << cache.m_currentStatic->volumeName() << "'.");
1938 if (m_navigator->atVolumeBoundary(nextPar, cache.m_currentStatic, dir, assocVol,
1939 m_tolerance) && assocVol != cache.m_currentStatic) {
1940 cache.m_currentDense = cache.m_dense ? nextVol : cache.m_highestVolume;
1941 }
1942 // no next volume found --- end of the world
1943 if (!nextVol) {
1944 ATH_MSG_DEBUG(" [+] World boundary reached - at " << positionOutput(
1945 nextPar->position()) << ", timed at " << cache.m_time);
1946 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
1947 return nextPar->uniqueClone();
1948 }
1949 // next volume found and parameters are at boundary
1950 if (nextVol /*&& nextPar nextPar is dereferenced anyway*/) {
1951 ATH_MSG_DEBUG(" [+] Crossing to next volume '" << nextVol->volumeName() << "'");
1952 ATH_MSG_DEBUG(" [+] Crossing position is - at " << positionOutput(nextPar->position()));
1953 if (!destVol && cache.m_currentStatic->geometrySignature() != nextVol->geometrySignature()) {
1954 nextGeoID = nextVol->geometrySignature();
1955 return nextPar->uniqueClone();
1956 }
1957 }
1958 cache.m_parametersAtBoundary.boundaryInformation(nextVol, nextPar, nextPar);
1959 return transportToVolumeWithPathLimit(cache,*nextPar, timeLim, dir, particle, nextGeoID, destVol);
1960 }
1961 if (dist > 0.) {
1962 return transportToVolumeWithPathLimit(cache,*nextPar, timeLim, dir, particle, nextGeoID, destVol);
1963 }
1964 } else if (sol < iDest + cache.m_trStaticBounds.size() + cache.m_trLays.size()) { // layer
1965 // material thickness - simple approach
1966 unsigned int index = sol - iDest - cache.m_trStaticBounds.size();
1967 const Trk::Layer *nextLayer = cache.m_navigLays[index].second;
1968
1969 bool matUp = nextLayer->layerMaterialProperties()->fullMaterial(nextPos) && m_includeMaterialEffects;
1970
1971 // material update
1972 if (matUp && m_includeMaterialEffects) {
1973 const ITimedMatEffUpdator *currentUpdator = subMaterialEffectsUpdator(*cache.m_currentStatic);
1974
1975 nextPar = currentUpdator ? currentUpdator
1976 ->update(nextPar,
1977 *nextLayer,
1978 timeLim,
1979 cache.m_path,
1980 cache.m_currentStatic->geometrySignature(),
1981 dir,
1982 particle)
1983 .release()
1984 : nextPar;
1985
1986 if (!nextPar) {
1987 ATH_MSG_VERBOSE(" [+] Update may have killed neutral track - return.");
1988 cache.m_parametersAtBoundary.resetBoundaryInformation();
1989 return returnParameters;
1990 }
1991 throwIntoGarbageBin(cache,nextPar);
1992
1993 }
1994 } else if (sol < iDest + cache.m_trStaticBounds.size() + cache.m_trLays.size() + cache.m_trDenseBounds.size()) {
1995 // dense volume boundary : no material update here, navigation only ( set cache.m_currentDense for next step )
1996
1997 unsigned int index = sol - iDest - cache.m_trStaticBounds.size() - cache.m_trLays.size();
1998 std::vector< std::pair<const Trk::TrackingVolume *, unsigned int> >::iterator dIter = cache.m_denseVols.begin();
1999 while (dIter != cache.m_denseVols.end() && index >= (*dIter).second) {
2000 index -= (*dIter).second;
2001 ++dIter;
2002 }
2003 if (dIter != cache.m_denseVols.end()) {
2004 currVol = (*dIter).first;
2005
2006 if (Trk::TrackingGeometry::atVolumeBoundary(nextPos, nextPar->momentum(), currVol, assocVol, dir,
2007 m_tolerance)) {
2008 if (assocVol && assocVol->zOverAtimesRho() != 0.) {
2009 cache.m_currentDense = assocVol;
2010 } else if (currVol->inside(nextPos + 0.002 * dir * nextPar->momentum().normalized())) {
2011 cache.m_currentDense = currVol;
2012 } else {
2013 // new search
2014 cache.m_currentDense = cache.m_highestVolume;
2015 if (m_useMuonMatApprox && cache.m_denseVols.empty()) {
2016 cache.m_currentDense = cache.m_currentStatic;
2017 } else {
2018 for (unsigned int i = 0; i < cache.m_denseVols.size(); i++) {
2019 const Trk::TrackingVolume *dVol = cache.m_denseVols[i].first;
2020 if (dVol->inside(nextPos + 0.002 * dir * nextPar->momentum().normalized(),
2021 m_tolerance) && dVol->zOverAtimesRho() != 0.) {
2022 cache.m_currentDense = dVol;
2023 }
2024 }
2025 }
2026 }
2027 }
2028 }
2029 } else { // detached volume bounds - not relevant ?
2030 }
2031
2032 throwIntoGarbageBin(cache,nextPar);
2033 }
2034
2035
2036
2037 if (nextPar) {
2038 ATH_MSG_DEBUG(
2039 " transportToVolumeWithPathLimit() - return from volume " << cache.m_currentStatic->volumeName() << " at position:" <<
2040 nextPar->position());
2041 return nextPar->uniqueClone();
2042 }
2043 return nullptr;
2044}
2045
2046Trk::BoundaryTrackParameters
2047Trk::TimedExtrapolator::transportInAlignableTV(Trk::TimedExtrapolator::Cache &cache,
2048 const Trk::TrackParameters &parm,
2049 Trk::TimeLimit &timeLim,
2050 Trk::PropDirection dir,
2051 Trk::ParticleHypothesis particle,
2052 Trk::GeometrySignature &nextGeoID,
2053 const Trk::AlignableTrackingVolume *aliTV) const {
2054 const std::string m = aliTV ? aliTV->volumeName() : " NULLPTR!";
2055 ATH_MSG_DEBUG(" [0] starting transport of neutral particle in alignable volume " << m);
2056
2057 // material loop in sensitive Calo volumes
2058 // returns: boundary parameters (static volume boundary)
2059 // material collection / intersection with active layers ( binned material used )
2060
2061 // initialize the return parameters vector
2062 const Trk::TrackParameters *currPar = &parm;
2063 const Trk::TrackingVolume *nextVol = nullptr;
2064 std::vector<Trk::IdentifiedIntersection> iis;
2065
2066 emptyGarbageBin(cache,&parm);
2067
2068 const EventContext& ctx = Gaudi::Hive::currentContext();
2069 if (!aliTV) {
2070 return {nullptr, nullptr, nullptr};
2071 }
2072
2073 // TODO if volume entry go to entry of misaligned volume
2074
2075 // save volume entry if collection present
2076
2077 const Trk::BinnedMaterial *binMat = aliTV->binnedMaterial();
2078
2079 const Trk::IdentifiedMaterial *binIDMat = nullptr;
2080
2081 const Trk::Material *currMat = aliTV; // material to be used
2082
2083
2084 // loop through binned material : save identifier, material, distance
2085
2086 // binned material
2087 if (binMat) {
2088 Amg::Vector3D pos = currPar->position();
2089 Amg::Vector3D pot = currPar->position();
2090 Amg::Vector3D umo = currPar->momentum().normalized();
2091
2092 binIDMat = binMat->material(pos);
2093
2094 if (cache.m_hitVector && binIDMat) {
2095 // std::cout <<"id info at the alignable volume entry:"<<binIDMat->second<<std::endl;
2096 if (binIDMat->second > 0) {
2097 cache.m_hitVector->emplace_back(currPar->uniqueClone(), timeLim.time, binIDMat->second, 0.);
2098 }
2099 }
2100
2101 const Trk::BinUtility *lbu = binMat->layerBinUtility(pos);
2102 if (lbu) {
2103 unsigned int cbin = lbu->bin(pos);
2104 // std::cout <<"layerBinUtility retrieved:"<<lbu->bins()<< std::endl;
2105 std::pair<size_t, float> d2n = lbu->distanceToNext(pos, dir * umo);
2106 // std::cout<<"estimated distance to the next bin:"<<d2n.first<<","<<d2n.second<< std::endl;
2107 float dTot = 0.;
2108 float distTot = 0.;
2109 // std::cout <<"input bin:"<<cbin<<", next: "<<d2n.first<<", at distance:"<<d2n.second<< std::endl;
2110 while (true) {
2111 if (d2n.first == cbin) {
2112 break;
2113 }
2114 dTot += d2n.second;
2115 distTot = dTot;
2116 pos = pos + d2n.second * dir * umo;
2117 if (!aliTV->inside(pos)) {
2118 break; // step outside volume
2119 }
2120 cbin = d2n.first;
2121 d2n = lbu->distanceToNext(pos, dir * umo);
2122 if (d2n.first == cbin && fabs(d2n.second) < 0.002) { // move ahead
2123 pos = pos + 0.002 * dir * umo;
2124 dTot += 0.002;
2125 d2n = lbu->distanceToNext(pos, dir * umo);
2126 }
2127 // std::cout <<"finding next bin?:"<<d2n.first<<","<<dTot<<"+"<<d2n.second<< std::endl;
2128 if (d2n.second > 0.001) { // retrieve material and save bin entry
2129 pot = pos + 0.5 * d2n.second * dir * umo;
2130 binIDMat = binMat->material(pot);
2131 iis.emplace_back(distTot, binIDMat->second, binIDMat->first.get());
2132 // std::cout <<"saving next bin entry:"<< distTot<<","<<binIDMat->second<<std::endl;
2133 }
2134 }
2135 }
2136 }
2137
2138 // resolve exit from the volume
2139
2140 cache.m_trStaticBounds.clear();
2141 const auto &bounds = aliTV->boundarySurfaces();
2142 for (unsigned int ib = 0; ib < bounds.size(); ib++) {
2143 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
2144 Trk::DistanceSolution distSol = surf.straightLineDistanceEstimate(currPar->position(),
2145 dir * currPar->momentum().normalized());
2146 double dist = distSol.first();
2147 // resolve multiple intersection solutions
2148 if (distSol.numberOfSolutions() > 1 && dist < m_tolerance && distSol.second() > dist) {
2149 dist = distSol.second();
2150 }
2151 // boundary check
2152 Amg::Vector3D gp = currPar->position() + dist * dir * currPar->momentum().normalized();
2153 // std::cout<<"alignable volume boundary:"<< ib<<","<<dist<<","<<
2154 // surf.isOnSurface(gp,true,m_tolerance,m_tolerance)<<std::endl;
2155 if (dist > m_tolerance && surf.isOnSurface(gp, true, m_tolerance, m_tolerance)) {
2156 const Trk::TrackingVolume *attachedVol = (bounds[ib])->attachedVolume(gp, currPar->momentum(), dir);
2157
2158 if (attachedVol && !(attachedVol->inside(gp + 0.01 * dir * currPar->momentum().normalized(), m_tolerance))) {
2159 ATH_MSG_DEBUG(
2160 " [!] WARNING: wrongly assigned exit volume ?" << cache.m_currentStatic->volumeName() << "->" <<
2161 attachedVol->volumeName());
2162 attachedVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
2163 gp + 0.01 * dir * currPar->momentum().normalized());
2164 if (attachedVol) {
2165 ATH_MSG_DEBUG(" new search yields: " << attachedVol->volumeName());
2166 }
2167 }
2168
2169 if (attachedVol != cache.m_currentStatic) { // exit
2170 nextVol = attachedVol;
2171 cache.m_trStaticBounds.insert(cache.m_trStaticBounds.begin(), Trk::DestBound(&surf, dist, ib));
2172 } else if (dist > 0.001) {
2173 const Trk::TrackingVolume *testVol = (bounds[ib])->attachedVolume(gp,
2174 currPar->momentum(),
2175 Trk::oppositeMomentum);
2176 ATH_MSG_WARNING(
2177 "gluing problem at the exit from alignable volume: " << gp.perp() << "," << gp.z() << ":" <<
2178 cache.m_currentStatic->volumeName());
2179 if (testVol) {
2180 ATH_MSG_DEBUG("inverted direction:" << testVol->volumeName());
2181 }
2182 if (testVol &&
2183 testVol->inside(gp + 0.01 * dir * currPar->momentum().normalized(),
2184 m_tolerance) && testVol != cache.m_currentStatic) {
2185 ATH_MSG_DEBUG(
2186 "next volume resolved to:" << testVol->volumeName() << " at the position(R,Z):" << gp.perp() << "," <<
2187 gp.z());
2188 nextVol = testVol;
2189 cache.m_trStaticBounds.insert(cache.m_trStaticBounds.begin(), Trk::DestBound(&surf, dist, ib));
2190 }
2191 }
2192 }
2193 } // end loop over boundaries
2194
2195 // if (nextVol) std::cout <<"nextVol, number of exit solutions:"<<
2196 // nextVol->volumeName()<<","<<cache.m_trStaticBounds.size()<< std::endl;
2197
2198 if (cache.m_trStaticBounds.empty()) {
2199 ATH_MSG_WARNING("exit from alignable volume " << aliTV->volumeName() << " not resolved, aborting");
2200 return {nullptr, nullptr, nullptr};
2201 } if (cache.m_trStaticBounds.size() > 1) { // hit edge ?
2202 Amg::Vector3D gp = currPar->position() + (cache.m_trStaticBounds[0].distance + 1.) * dir *
2203 currPar->momentum().normalized();
2204 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
2205 ATH_MSG_DEBUG("exit volume reassigned:" << nextVol->volumeName());
2206 }
2207
2208 // exit from the volume may coincide with the last bin boundary - leave 10 microns marge
2209 if (!iis.empty() && cache.m_trStaticBounds[0].distance - iis.back().distance < 0.01) {
2210 iis.pop_back();
2211 }
2212
2213 // add volume exit
2214 iis.emplace_back(cache.m_trStaticBounds[0].distance, 0, nullptr);
2215
2216 // loop over intersection taking into account the material effects
2217
2218 double dist = 0.;
2219 double mom = currPar->momentum().mag();
2220 double beta = mom / sqrt(mom * mom + cache.m_particleMass * cache.m_particleMass) * Gaudi::Units::c_light;
2221 Amg::Vector3D nextPos = currPar->position();
2222
2223 int currLay = 0;
2224
2225 for (unsigned int is = 0; is < iis.size(); is++) {
2226 if (iis[is].distance == 0.) {
2227 continue;
2228 }
2229
2230 double step = iis[is].distance - dist;
2231
2232 nextPos = currPar->position() + dir * currPar->momentum().normalized() * iis[is].distance;
2233
2234 double tDelta = step / beta;
2235
2236 double mDelta = (currMat->zOverAtimesRho() != 0.) ? step / currMat->x0() : 0.;
2237
2238 // in case of hadronic interaction retrieve nuclear interaction properties, too
2239
2240 double frT = 1.;
2241 if (step > 0 && timeLim.tMax > cache.m_time && cache.m_time + tDelta >= timeLim.tMax) {
2242 frT = (timeLim.tMax - cache.m_time) * beta / step;
2243 }
2244
2245 // TODO : compare x0 or l0 according to the process type
2246 double frM = 1.;
2247 if (mDelta > 0 && cache.m_path.x0Max > 0.) {
2248 if (cache.m_path.process < 100 && cache.m_path.x0Collected + mDelta > cache.m_path.x0Max) {
2249 frM = (cache.m_path.x0Max - cache.m_path.x0Collected) / mDelta;
2250 } else { // waiting for hadronic interaction, retrieve nuclear interaction properties
2251 double mDeltaL = currMat->L0 > 0. ? step / currMat->L0 : mDelta / 0.37 / currMat->averageZ();
2252 if (cache.m_path.l0Collected + mDeltaL > cache.m_path.x0Max) {
2253 frM = (cache.m_path.x0Max - cache.m_path.l0Collected) / mDeltaL;
2254 }
2255 }
2256 }
2257
2258 double fr = fmin(frT, frM);
2259
2260 // std::cout << "looping over intersections:"<<is<<","<< cache.m_trSurfs[sols[is]].second<<","<<step << ","<<
2261 // tDelta<<","<<mDelta << std::endl;
2262
2263 if (fr < 1.) { // decay or material interaction during the step
2264 int process = frT < frM ? timeLim.process : cache.m_path.process;
2265 cache.m_time += fr * step / beta;
2266 if (mDelta > 0 && currMat->averageZ() > 0) {
2267 cache.m_path.updateMat(fr * mDelta, currMat->averageZ(), 0.);
2268 }
2269
2270 nextPos = currPar->position() + dir * currPar->momentum().normalized() * (dist + fr * step);
2271
2272 // process interaction only if creation of secondaries allowed
2273 if (m_caloMsSecondary) {
2274 const Trk::TrackParameters* nextPar =
2275 m_updators[0]
2276 ->interact(cache.m_time, nextPos, currPar->momentum(), particle, process, currMat)
2277 .release();
2278 throwIntoGarbageBin(cache, nextPar);
2279
2280 if (nextPar) {
2281 ATH_MSG_DEBUG(" [!] WARNING: particle survives the interaction " << process);
2282 }
2283
2284 if (nextPar && process == 121) {
2285 ATH_MSG_DEBUG(" [!] WARNING: failed hadronic interaction, killing the input particle anyway");
2286 return {nullptr, nullptr, nullptr};
2287 }
2288
2289 if (!nextPar) {
2290 return {nullptr, nullptr, nullptr};
2291 }
2292
2293 // return transportToVolumeWithPathLimit(*nextPar, timeLim, dir, particle, nextGeoID, destVol);
2294 } else { // kill particle without trace ?
2295 return {nullptr, nullptr, nullptr};
2296 }
2297 } // end decay or material interaction during the step
2298
2299 // update
2300 dist = iis[is].distance;
2301 if (mDelta > 0 && currMat->averageZ() > 0) {
2302 cache.m_path.updateMat(mDelta, currMat->averageZ(), 0.);
2303 }
2304 cache.m_time += tDelta;
2305
2306 if (is < iis.size() - 1) { // update bin material info
2307 // binIDMat = binMat->material(nextPos);
2308 // currMat = binIDMat->first;
2309 currMat = iis[is].material;
2310 currLay = iis[is].identifier;
2311
2312 if (cache.m_hitVector && iis[is].identifier > 0) { // save entry to the next layer
2313 ATH_MSG_VERBOSE("active layer entry:" << currLay << " at R,z:" << nextPos.perp() << "," << nextPos.z());
2314 auto nextPar = std::make_unique<Trk::CurvilinearParameters>(nextPos, currPar->momentum(), 0.);
2315 cache.m_hitVector->emplace_back(std::move(nextPar), timeLim.time, iis[is].identifier, 0.);
2316 }
2317 }
2318 } // end loop over intersections
2319
2320 Trk::CurvilinearParameters *nextPar = new Trk::CurvilinearParameters(nextPos, currPar->momentum(), 0.);
2321
2322 if (cache.m_hitVector) { // save volume exit /active layer only ?
2323 ATH_MSG_VERBOSE("active layer/volume exit:" << currLay << " at R,z:" << nextPos.perp() << "," << nextPos.z());
2324 if (binIDMat and(binIDMat->second > 0)) {
2325 cache.m_hitVector->emplace_back(nextPar->uniqueClone(), timeLim.time, currLay, 0.);
2326 }
2327 }
2328
2329 throwIntoGarbageBin(cache,nextPar);
2330
2331
2332
2333 ATH_MSG_DEBUG(" [+] StaticVol boundary reached of '" << cache.m_currentStatic->volumeName() << "'.");
2334
2335 // no next volume found --- end of the world
2336 if (!nextVol) {
2337 ATH_MSG_DEBUG(" [+] World boundary reached - at " << positionOutput(
2338 nextPar->position()) << ", timed at " << cache.m_time);
2339 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
2340 } else {
2341 ATH_MSG_DEBUG(" [+] Crossing to next volume '" << nextVol->volumeName() << "'");
2342 ATH_MSG_DEBUG(" [+] Crossing position is - at " << positionOutput(nextPar->position()));
2343 }
2344
2345 cache.m_parametersAtBoundary.boundaryInformation(nextVol, nextPar, nextPar);
2346
2347 return {nextPar, nextVol, cache.m_currentStatic};
2348}
2349
2350Trk::BoundaryTrackParameters
2351Trk::TimedExtrapolator::extrapolateInAlignableTV(Trk::TimedExtrapolator::Cache &cache,
2352 const Trk::TrackParameters &parm,
2353 Trk::TimeLimit &timeLim,
2354 Trk::PropDirection dir,
2355 Trk::ParticleHypothesis particle,
2356 Trk::GeometrySignature &nextGeoID,
2357 const Trk::AlignableTrackingVolume *vol) const {
2358 const std::string m = vol ? vol->volumeName():"NULLPTR";
2359 ATH_MSG_DEBUG("M-[" << ++cache.m_methodSequence << "] extrapolateInAlignableTV(...) " << m);
2360
2361 // material loop in sensitive Calo volumes
2362 // extrapolation without target surface returns:
2363 // A) boundary parameters (static volume boundary)
2364 // if target surface:
2365 // B) trPar at target surface
2366 // material collection done by the propagator ( binned material used )
2367
2368 // initialize the return parameters vector
2369 const Trk::TrackParameters *currPar = &parm;
2370 const Trk::AlignableTrackingVolume *staticVol = nullptr;
2371 const Trk::TrackingVolume *currVol = nullptr;
2372 const Trk::TrackingVolume *nextVol = nullptr;
2373 std::vector<unsigned int> solutions;
2374 // double tol = 0.001;
2375 // double path = 0.;
2376 const EventContext& ctx = Gaudi::Hive::currentContext();
2377 if (!cache.m_highestVolume) {
2378 cache.m_highestVolume = m_navigator->highestVolume(ctx);
2379 }
2380
2381 emptyGarbageBin(cache,&parm);
2382
2383 // verify current position
2384 const Amg::Vector3D& gp = parm.position();
2385 if (vol && vol->inside(gp, m_tolerance)) {
2386 staticVol = vol;
2387 } else {
2388 currVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
2389 const Trk::TrackingVolume *nextStatVol = nullptr;
2390 if (m_navigator->atVolumeBoundary(currPar, currVol, dir, nextStatVol, m_tolerance) && nextStatVol != currVol) {
2391 currVol = nextStatVol;
2392 }
2393 if (currVol && currVol != vol) {
2394 const Trk::AlignableTrackingVolume *aliTG = dynamic_cast<const Trk::AlignableTrackingVolume *> (currVol);
2395 if (aliTG) {
2396 staticVol = aliTG;
2397 }
2398 }
2399 }
2400
2401 if (!staticVol) {
2402 ATH_MSG_DEBUG(" [!] failing in retrieval of AlignableTV, return 0");
2403 return {nullptr, nullptr, nullptr};
2404 }
2405
2406 // TODO if volume entry go to entry of misaligned volume
2407
2408 // save volume entry if collection present
2409
2410 if (cache.m_hitVector) {
2411 const Trk::BinnedMaterial *binMat = staticVol->binnedMaterial();
2412 if (binMat) {
2413 const Trk::IdentifiedMaterial *binIDMat = binMat->material(currPar->position());
2414 if (binIDMat->second > 0) {
2415 cache.m_hitVector->emplace_back(currPar->uniqueClone(), timeLim.time, binIDMat->second, 0.);
2416 }
2417 }
2418 }
2419
2420 // navigation surfaces
2421 if (cache.m_navigSurfs.capacity() > m_maxNavigSurf) {
2422 cache.m_navigSurfs.reserve(m_maxNavigSurf);
2423 }
2424 cache.m_navigSurfs.clear();
2425
2426 // assume new static volume, retrieve boundaries
2427 cache.m_currentStatic = staticVol;
2428 cache.m_staticBoundaries.clear();
2429 const auto &bounds = staticVol->boundarySurfaces();
2430 for (unsigned int ib = 0; ib < bounds.size(); ++ib) {
2431 const Trk::Surface &surf = (bounds[ib])->surfaceRepresentation();
2432 cache.m_staticBoundaries.emplace_back(&surf, true);
2433 }
2434
2435 cache.m_navigSurfs.insert(cache.m_navigSurfs.end(), cache.m_staticBoundaries.begin(), cache.m_staticBoundaries.end());
2436
2437 // current dense
2438 cache.m_currentDense = staticVol;
2439
2440 // ready to propagate
2441 // till: A/ static volume boundary(bcheck=true) , B/ destination surface(bcheck=false)
2442
2443 nextVol = nullptr;
2444 while (currPar) {
2445 std::vector<unsigned int> solutions;
2446 // propagate now
2447 ATH_MSG_DEBUG(" [+] Starting propagation at position " << positionOutput(currPar->position())
2448 << " (current momentum: " << currPar->momentum().mag() <<
2449 ")");
2450 ATH_MSG_DEBUG(" [+] " << cache.m_navigSurfs.size() << " target surfaces in '" << cache.m_currentDense->volumeName() << "'.");
2451 // arguments : inputParameters, vector of navigation surfaces, propagation direction, b field service, particle
2452 // type, result,
2453 // material collection, intersection collection, path limit, switch for use of path limit, switch for
2454 // curvilinear on return, current TG volume
2455 const Trk::TrackParameters* nextPar = m_stepPropagator
2456 ->propagateT(ctx,
2457 *currPar,
2458 cache.m_navigSurfs,
2459 dir,
2460 m_fieldProperties,
2461 particle,
2462 solutions,
2463 cache.m_path,
2464 timeLim,
2465 true,
2466 cache.m_currentDense,
2467 cache.m_hitVector)
2468 .release();
2469 ATH_MSG_VERBOSE(" [+] Propagation done. ");
2470 if (nextPar) {
2471 ATH_MSG_DEBUG(" [+] Position after propagation - at " << positionOutput(nextPar->position()));
2472 }
2473
2474 if (nextPar) {
2475 ATH_MSG_DEBUG(" [+] Number of intersection solutions: " << solutions.size());
2476 }
2477 if (nextPar) {
2478 throwIntoGarbageBin(cache,nextPar);
2479 }
2480
2481 // material update has been done already by the propagator
2482 if (cache.m_path.x0Max > 0. &&
2483 ((cache.m_path.process < 100 && cache.m_path.x0Collected >= cache.m_path.x0Max) ||
2484 (cache.m_path.process > 100 && cache.m_path.l0Collected >= cache.m_path.x0Max))) {
2485 // trigger presampled interaction, provide material properties if needed
2486 // process interaction only if creation of secondaries allowed
2487 if (cache.m_currentStatic->geometrySignature() == Trk::ID || m_caloMsSecondary) {
2488 const Trk::Material *extMprop = cache.m_path.process > 100 ? cache.m_currentDense : nullptr;
2489
2490 const Trk::TrackParameters *iPar = nullptr;
2491 if (nextPar) {
2492 iPar =
2493 m_updators[0]
2494 ->interact(
2495 timeLim.time, nextPar->position(), nextPar->momentum(), particle, cache.m_path.process, extMprop)
2496 .release();
2497 }
2498
2499 if (!iPar) {
2500 return {nullptr, nullptr, nullptr};
2501 }
2502
2503 throwIntoGarbageBin(cache,iPar);
2504
2505 if (iPar && cache.m_path.process == 121) {
2506 ATH_MSG_DEBUG(" [!] WARNING: failed hadronic interaction, killing the input particle anyway");
2507 return {nullptr, nullptr, nullptr};
2508 }
2509
2510 // return transportToVolumeWithPathLimit(*nextPar, timeLim, dir, particle, nextGeoID, destVol);
2511 } else { // kill particle without trace ?
2512 return {nullptr, nullptr, nullptr};
2513 }
2514 }
2515
2516 // decay ?
2517 if (timeLim.tMax > 0. && timeLim.time >= timeLim.tMax) {
2518 // process interaction only if creation of secondaries allowed
2519 if (cache.m_currentStatic->geometrySignature() == Trk::ID || m_caloMsSecondary) {
2520 // trigger presampled interaction
2521 const Trk::TrackParameters* iPar = m_updators[0]->interact(
2522 timeLim.time, nextPar->position(), nextPar->momentum(), particle, timeLim.process).release();
2523 if (!iPar) {
2524 return {nullptr, nullptr, nullptr};
2525 }
2526
2527 throwIntoGarbageBin(cache,iPar);
2528 ATH_MSG_WARNING("particle decay survival?" << particle << "," << timeLim.process);
2529 return {nullptr, nullptr, nullptr};
2530 } // kill the particle without trace ( some validation info can be included here eventually )
2531 return {nullptr, nullptr, nullptr};
2532
2533 }
2534
2535 if (nextPar) {
2536 unsigned int iSol = 0;
2537 while (iSol < solutions.size()) {
2538 if (solutions[iSol] < cache.m_staticBoundaries.size()) {
2539 // TODO if massive boundary coded, add the material effects here
2540 // static volume boundary; return to the main loop : TODO move from misaligned to static
2541 unsigned int index = solutions[iSol];
2542 // use global coordinates to retrieve attached volume (just for static!)
2543 nextVol = (cache.m_currentStatic->boundarySurfaces())[index]->attachedVolume(
2544 nextPar->position(), nextPar->momentum(), dir);
2545 // double check the next volume
2546 if (nextVol && !(nextVol->inside(nextPar->position() + 0.01 * dir * nextPar->momentum().normalized(), 0.))) {
2547 ATH_MSG_DEBUG(
2548 " [!] WARNING: wrongly assigned static volume ?" << cache.m_currentStatic->volumeName() << "->" <<
2549 nextVol->volumeName());
2550 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
2551 nextPar->position() + 0.01 * dir * nextPar->momentum().normalized());
2552 if (nextVol) {
2553 ATH_MSG_DEBUG(" new search yields: " << nextVol->volumeName());
2554 }
2555 }
2556 // end double check - to be removed after validation of the geometry gluing
2557 // lateral exit from calo sample can be handled here
2558 if (cache.m_hitVector) {
2559 const Trk::BinnedMaterial *binMat = staticVol->binnedMaterial();
2560 if (binMat) {
2561 const Trk::IdentifiedMaterial *binIDMat = binMat->material(nextPar->position());
2562 // save only if entry to the sample present, the exit missing and non-zero step in the sample
2563 if (binIDMat && binIDMat->second > 0 && !cache.m_hitVector->empty() &&
2564 cache.m_hitVector->back().detID == binIDMat->second) {
2565 // double s = (nextPar->position()-m_identifiedParameters->back().first->position()).mag();
2566 // if (s>0.001) m_identifiedParameters->push_back(std::pair<const Trk::TrackParameters*,int>
2567 // (nextPar->clone(), -binIDMat->second));
2568 cache.m_hitVector->emplace_back(nextPar->uniqueClone(), timeLim.time, -binIDMat->second, 0.);
2569 }
2570 }
2571 }
2572 // end lateral exit handling
2573
2574 ATH_MSG_DEBUG(" [+] StaticVol boundary reached of '" << cache.m_currentStatic->volumeName() << "'.");
2575 // no next volume found --- end of the world
2576 if (!nextVol) {
2577 ATH_MSG_DEBUG(" [+] World boundary reached - at " << positionOutput(
2578 nextPar->position()) << ", timed at " << cache.m_time);
2579 nextGeoID = Trk::GeometrySignature(Trk::Unsigned);
2580 } else {
2581 ATH_MSG_DEBUG(" [+] Crossing to next volume '" << nextVol->volumeName() << "'");
2582 ATH_MSG_DEBUG(" [+] Crossing position is - at " << positionOutput(nextPar->position()));
2583 }
2584
2585 return {nextPar, nextVol, cache.m_currentStatic};
2586 }
2587 }
2588 }
2589
2590 currPar = nextPar;
2591 }
2592
2593 return {nullptr, nullptr, nullptr};
2594}
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
size
size_t size() const
Number of registered mappings.
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
Trk::AlignableTrackingVolume
Base Class for a navigation object (active) in the Calo realm.
Definition AlignableTrackingVolume.h:36
Trk::AlignableTrackingVolume::binnedMaterial
const BinnedMaterial * binnedMaterial() const
access to binned material
Definition AlignableTrackingVolume.h:68
Trk::BinUtility
A generic symmetric BinUtility, for fully symmetric binning in terms of binning grid and binning type...
Definition BinUtility.h:39
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::BinUtility::bin
size_t bin(const Amg::Vector3D &position, size_t ba=0) const
Bin from a 3D vector (already in binning frame).
Definition BinUtility.h:126
Trk::BinnedArray::arrayObjects
virtual std::span< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
Trk::BinnedMaterial
for description of non-homogenous dense volumes
Definition BinnedMaterial.h:31
Trk::BinnedMaterial::material
const IdentifiedMaterial * material(const Amg::Vector3D &position) const
access to material/id per bin
Definition BinnedMaterial.cxx:58
Trk::BinnedMaterial::layerBinUtility
const Trk::BinUtility * layerBinUtility(const Amg::Vector3D &position) const
access to layer bin utility
Definition BinnedMaterial.h:80
Trk::BoundaryCheck
The BoundaryCheck class allows to steer the way surface boundaries are used for inside/outside checks...
Definition BoundaryCheck.h:51
Trk::CylinderVolumeBounds
Bounds for a cylindrical Volume, the decomposeToSurfaces method creates a vector of up to 6 surfaces:
Definition CylinderVolumeBounds.h:70
Trk::CylinderVolumeBounds::innerRadius
double innerRadius() const
This method returns the inner radius.
Definition CylinderVolumeBounds.h:187
Trk::CylinderVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflengthZ.
Definition CylinderVolumeBounds.h:207
Trk::CylinderVolumeBounds::outerRadius
double outerRadius() const
This method returns the outer radius.
Definition CylinderVolumeBounds.h:191
Trk::DetachedTrackingVolume
Base Class for a navigation object (active/passive) in the Tracking realm.
Definition DetachedTrackingVolume.h:37
Trk::DistanceSolution
Access to distance solutions.
Definition DistanceSolution.h:25
Trk::DistanceSolution::second
double second() const
Distance to second intersection solution along direction (for a cylinder surface).
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::DistanceSolution::numberOfSolutions
int numberOfSolutions() const
Number of intersection solutions.
Trk::DistanceSolution::first
double first() const
Distance to first intersection solution along direction.
Trk::IPropagator
Interface class IPropagators It inherits from IAlgTool.
Definition IPropagator.h:54
Trk::IPropagator::propagate
virtual std::unique_ptr< NeutralParameters > propagate(const NeutralParameters &parameters, const Surface &sf, PropDirection dir, const BoundaryCheck &bcheck, bool returnCurv=false) const =0
Main propagation method for NeutralParameters.
Trk::ITimedMatEffUpdator
Interface class for the updater AlgTool, it inherits from IAlgTool Detailed information about private...
Definition ITimedMatEffUpdator.h:43
Trk::ITimedMatEffUpdator::update
virtual std::unique_ptr< TrackParameters > update(const TrackParameters *parm, const Layer &sf, TimeLimit &time, PathLimit &path, Trk::GeometrySignature geoID, PropDirection dir=alongMomentum, ParticleHypothesis particle=pion) const =0
Updator interface (full update for a layer): The parmeters are given as a const pointer.
Trk::LayerMaterialProperties::fullMaterial
virtual const MaterialProperties * fullMaterial(const Amg::Vector3D &gp) const =0
Return method for full material description of the Layer.
Trk::Layer
Base Class for a Detector Layer in the Tracking realm.
Definition Layer.h:72
Trk::Layer::compatibleSurfaces
size_t compatibleSurfaces(std::vector< SurfaceIntersection > &cSurfaces, const TrackParameters &pars, PropDirection pdir, const BoundaryCheck &bcheck, bool materialSurfacesOnly=true, const Surface *startSurface=nullptr, const Surface *endSurface=nullptr, const ICompatibilityEstimator *ice=nullptr) const
get compatible surfaces starting from charged parameters
Trk::Layer::layerType
int layerType() const
get the Layer coding
Trk::Layer::surfaceRepresentation
virtual const Surface & surfaceRepresentation() const =0
Transforms the layer into a Surface representation for extrapolation.
Trk::Layer::fullUpdateMaterialProperties
const MaterialProperties * fullUpdateMaterialProperties(const TrackParameters &par) const
getting the MaterialProperties back - for full update
Definition Layer.cxx:169
Trk::Layer::isOnLayer
virtual bool isOnLayer(const Amg::Vector3D &gp, const BoundaryCheck &bcheck=BoundaryCheck(true)) const
isOnLayer() method, using isOnSurface() with Layer specific tolerance
Definition Layer.cxx:135
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
Trk::Layer::surfaceArray
const SurfaceArray * surfaceArray() const
Return the entire SurfaceArray, returns nullptr if no SurfaceArray.
Trk::Layer::layerMaterialProperties
const LayerMaterialProperties * layerMaterialProperties() const
getting the LayerMaterialProperties including full/pre/post update
Trk::Layer::subSurface
const Surface * subSurface(const Amg::Vector3D &gp) const
If no subSurface array is defined or no subSurface can be found to the given Amg::Vector3D,...
Definition Layer.cxx:107
Trk::MagneticFieldProperties
magnetic field properties to steer the behavior of the extrapolation
Definition MagneticFieldProperties.h:31
Trk::Material
A common object to be contained by.
Definition Material.h:117
Trk::Material::x0
float x0() const
Definition Material.h:227
Trk::Material::averageZ
float averageZ() const
Definition Material.h:228
Trk::Material::zOverAtimesRho
float zOverAtimesRho() const
access to members
Definition Material.h:226
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
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::ParametersCommon::localPosition
Amg::Vector2D localPosition() const
Access method for the local coordinates, local parameter definitions differ for each surface type.
Trk::Surface
Abstract Base Class for tracking surfaces.
Definition Surface.h:79
Trk::Surface::associatedDetectorElement
const TrkDetElementBase * associatedDetectorElement() const
return associated Detector Element
Trk::Surface::straightLineDistanceEstimate
virtual DistanceSolution straightLineDistanceEstimate(const Amg::Vector3D &pos, const Amg::Vector3D &dir) const =0
fast straight line distance evaluation to Surface
Trk::Surface::normal
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
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::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Trk::TimedExtrapolator::transportInAlignableTV
BoundaryTrackParameters transportInAlignableTV(Trk::TimedExtrapolator::Cache &cache, const Trk::TrackParameters &parm, Trk::TimeLimit &time, Trk::PropDirection dir, Trk::ParticleHypothesis particle, Trk::GeometrySignature &nextGeoId, const Trk::AlignableTrackingVolume *aliTV) const
Definition TimedExtrapolator.cxx:2047
Trk::TimedExtrapolator::m_propagators
ToolHandleArray< IPropagator > m_propagators
Definition TimedExtrapolator.h:293
Trk::TimedExtrapolator::TimedExtrapolator
TimedExtrapolator(const std::string &, const std::string &, const IInterface *)
Constructor.
Definition TimedExtrapolator.cxx:56
Trk::TimedExtrapolator::m_resolveActive
BooleanProperty m_resolveActive
Definition TimedExtrapolator.h:357
Trk::TimedExtrapolator::m_printRzOutput
BooleanProperty m_printRzOutput
Definition TimedExtrapolator.h:363
Trk::TimedExtrapolator::m_useDenseVolumeDescription
BooleanProperty m_useDenseVolumeDescription
Definition TimedExtrapolator.h:352
Trk::TimedExtrapolator::finalize
virtual StatusCode finalize() override
AlgTool finalize method.
Definition TimedExtrapolator.cxx:187
Trk::TimedExtrapolator::transportNeutralsWithPathLimit
virtual std::unique_ptr< const Trk::TrackParameters > transportNeutralsWithPathLimit(const Trk::TrackParameters &parm, Trk::PathLimit &pathLim, Trk::TimeLimit &time, Trk::PropDirection dir, Trk::ParticleHypothesis particle, std::vector< Trk::HitInfo > *&hitVector, Trk::GeometrySignature &nextGeoId, const Trk::TrackingVolume *boundaryVol=nullptr) const override
Transport method for neutral, possibly unstable particles.
Definition TimedExtrapolator.cxx:1268
Trk::TimedExtrapolator::m_resolveMultilayers
BooleanProperty m_resolveMultilayers
Definition TimedExtrapolator.h:358
Trk::TimedExtrapolator::extrapolateToVolumeWithPathLimit
std::unique_ptr< const Trk::TrackParameters > extrapolateToVolumeWithPathLimit(Trk::TimedExtrapolator::Cache &cache, const Trk::TrackParameters &parm, Trk::TimeLimit &time, Trk::PropDirection dir, Trk::ParticleHypothesis particle, Trk::GeometrySignature &nextGeoID, const Trk::TrackingVolume *destVol) const
Definition TimedExtrapolator.cxx:272
Trk::TimedExtrapolator::m_navigator
ToolHandle< INavigator > m_navigator
Definition TimedExtrapolator.h:297
Trk::TimedExtrapolator::m_updators
ToolHandleArray< ITimedMatEffUpdator > m_updators
Definition TimedExtrapolator.h:300
Trk::TimedExtrapolator::m_propNames
StringArrayProperty m_propNames
Definition TimedExtrapolator.h:316
Trk::TimedExtrapolator::m_msupdators
ToolHandleArray< IMultipleScatteringUpdator > m_msupdators
Definition TimedExtrapolator.h:302
Trk::TimedExtrapolator::m_fieldProperties
Trk::MagneticFieldProperties m_fieldProperties
Definition TimedExtrapolator.h:436
Trk::TimedExtrapolator::m_stepPropagator
ToolHandle< IPropagator > m_stepPropagator
Definition TimedExtrapolator.h:295
Trk::TimedExtrapolator::m_useMuonMatApprox
BooleanProperty m_useMuonMatApprox
Definition TimedExtrapolator.h:355
Trk::TimedExtrapolator::m_includeMaterialEffects
BooleanProperty m_includeMaterialEffects
Definition TimedExtrapolator.h:328
Trk::TimedExtrapolator::extrapolateInAlignableTV
BoundaryTrackParameters extrapolateInAlignableTV(Trk::TimedExtrapolator::Cache &cache, const Trk::TrackParameters &parm, Trk::TimeLimit &time, Trk::PropDirection dir, Trk::ParticleHypothesis particle, Trk::GeometrySignature &nextGeoId, const Trk::AlignableTrackingVolume *aliTV) const
Definition TimedExtrapolator.cxx:2351
Trk::TimedExtrapolator::validationAction
virtual void validationAction() const override
Validation Action: Can be implemented optionally, outside access to internal validation steps.
Definition TimedExtrapolator.cxx:1259
Trk::TimedExtrapolator::momentumOutput
static std::string momentumOutput(const Amg::Vector3D &mom)
For the output - global momentum.
Definition TimedExtrapolator.cxx:1226
Trk::TimedExtrapolator::~TimedExtrapolator
virtual ~TimedExtrapolator()
Destructor.
Trk::TimedExtrapolator::positionOutput
std::string positionOutput(const Amg::Vector3D &pos) const
Private method for conversion of the synchronized geometry signature to the natural subdetector order...
Definition TimedExtrapolator.cxx:1214
Trk::TimedExtrapolator::extrapolateWithPathLimit
virtual std::unique_ptr< const Trk::TrackParameters > extrapolateWithPathLimit(const Trk::TrackParameters &parm, Trk::PathLimit &pathLim, Trk::TimeLimit &time, Trk::PropDirection dir, Trk::ParticleHypothesis particle, std::vector< Trk::HitInfo > *&hitVector, Trk::GeometrySignature &nextGeoID, const Trk::TrackingVolume *boundaryVol=nullptr) const override
Extrapolation method for charged, possibly unstable particles.
Definition TimedExtrapolator.cxx:193
Trk::TimedExtrapolator::m_subPropagators
std::vector< const IPropagator * > m_subPropagators
Propagators to chose from (steered by signature).
Definition TimedExtrapolator.h:310
Trk::TimedExtrapolator::m_configurationLevel
unsigned int m_configurationLevel
see the supported levels of configuration above
Definition TimedExtrapolator.h:326
Trk::TimedExtrapolator::m_subUpdators
std::vector< const ITimedMatEffUpdator * > m_subUpdators
Updators to chose from (steered by signature).
Definition TimedExtrapolator.h:312
Trk::TimedExtrapolator::emptyGarbageBin
void emptyGarbageBin(Trk::TimedExtrapolator::Cache &cache, const Trk::TrackParameters *) const
Private method for emptying the GarbageBin.
Definition TimedExtrapolator.cxx:1234
Trk::TimedExtrapolator::m_updatNames
StringArrayProperty m_updatNames
Definition TimedExtrapolator.h:318
Trk::TimedExtrapolator::transportToVolumeWithPathLimit
std::unique_ptr< const Trk::TrackParameters > transportToVolumeWithPathLimit(Trk::TimedExtrapolator::Cache &cache, const Trk::TrackParameters &parm, Trk::TimeLimit &time, Trk::PropDirection dir, Trk::ParticleHypothesis particle, Trk::GeometrySignature &nextGeoId, const Trk::TrackingVolume *boundaryVol) const
Definition TimedExtrapolator.cxx:1337
Trk::TimedExtrapolator::m_caloMsSecondary
BooleanProperty m_caloMsSecondary
Definition TimedExtrapolator.h:346
Trk::TimedExtrapolator::subMaterialEffectsUpdator
const ITimedMatEffUpdator * subMaterialEffectsUpdator(const TrackingVolume &tvol) const
Access the subPropagator to the given volume.
Definition TimedExtrapolator.h:464
Trk::TimedExtrapolator::initialize
virtual StatusCode initialize() override
AlgTool initailize method.
Definition TimedExtrapolator.cxx:69
Trk::TimedExtrapolator::throwIntoGarbageBin
void throwIntoGarbageBin(Trk::TimedExtrapolator::Cache &cache, const Trk::TrackParameters *garbage) const
Private method for throwing into the GarbageBin.
Definition TimedExtrapolator.h:471
Trk::TimedExtrapolator::m_robustSampling
BooleanProperty m_robustSampling
Definition TimedExtrapolator.h:351
Trk::TimedExtrapolator::overlapSearch
void overlapSearch(Trk::TimedExtrapolator::Cache &cache, const IPropagator &prop, const TrackParameters &parm, const TrackParameters &parsOnLayer, const Layer &lay, float time, PropDirection dir=anyDirection, const BoundaryCheck &bcheck=true, ParticleHypothesis particle=pion, bool startingLayer=false) const
Private to search for overlap surfaces.
Definition TimedExtrapolator.cxx:1041
Trk::TrackingGeometry::atVolumeBoundary
static bool atVolumeBoundary(const Amg::Vector3D &gp, const TrackingVolume *vol, double tol)
check position at volume boundary
Definition TrackingGeometry.cxx:244
Trk::TrackingVolume
Full Volume description used in Tracking, it inherits from Volume to get the geometrical structure,...
Definition TrackingVolume.h:119
Trk::TrackingVolume::confinedLayers
const LayerArray * confinedLayers() const
Return the subLayer array.
Trk::TrackingVolume::geometrySignature
GeometrySignature geometrySignature() const
return the Signature
Trk::TrackingVolume::associatedLayer
const Layer * associatedLayer(const Amg::Vector3D &gp) const
Return the associated Layer.
Definition TrackingVolume.cxx:382
Trk::TrackingVolume::confinedVolumes
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
Trk::TrackingVolume::associatedSubVolume
const TrackingVolume * associatedSubVolume(const Amg::Vector3D &gp) const
Return the associated sub Volume, returns THIS if no subVolume exists.
Definition TrackingVolume.cxx:530
Trk::TrackingVolume::boundarySurfaces
std::vector< std::shared_ptr< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
Definition TrackingVolume.cxx:809
Trk::TrackingVolume::volumeName
const std::string & volumeName() const
Returns the VolumeName - for debug reason, might be depreciated later.
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::TrackingVolume::confinedArbitraryLayers
ArraySpan< Layer const *const > confinedArbitraryLayers() const
Return the confined subLayer array.
Trk::TrackingVolume::confinedDenseVolumes
ArraySpan< TrackingVolume const *const > confinedDenseVolumes() const
Return unordered subVolumes - not the ownership.
Trk::Volume::inside
bool inside(const Amg::Vector3D &gp, double tol=0.) const
Inside() method for checks.
Definition Volume.cxx:72
process
const std::string process
Definition fbtTestBasics.cxx:76
replace
std::string replace(std::string s, const std::string &s2, const std::string &s3)
Definition hcg.cxx:312
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition ParticleHypothesis.h:56
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::PropDirection
PropDirection
PropDirection, enum for direction of the propagation.
Definition PropDirection.h:19
Trk::oppositeMomentum
@ oppositeMomentum
Definition PropDirection.h:21
Trk::anyDirection
@ anyDirection
Definition PropDirection.h:22
Trk::ArraySpan
std::span< T > ArraySpan
Definition TrackingVolume.h:59
Trk::CurvilinearParameters
CurvilinearParametersT< TrackParametersDim, Charged, PlaneSurface > CurvilinearParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:29
Trk::TrkParametersComparisonFunction
ComparisonFunction< TrackParameters > TrkParametersComparisonFunction
Definition TrkParametersComparisonFunction.h:26
Trk::FastField
@ FastField
call the fast field access method of the FieldSvc
Definition MagneticFieldMode.h:20
Trk::FullField
@ FullField
Field is set to be realistic, but within a given Volume.
Definition MagneticFieldMode.h:21
Trk::ParticleHypothesis
ParticleHypothesis
Enumeration for Particle hypothesis respecting the interaction with material.
Definition ParticleHypothesis.h:28
Trk::IdentifiedMaterial
std::pair< std::shared_ptr< Material >, int > IdentifiedMaterial
Definition BinnedMaterial.h:23
Trk::NumberOfSignatures
@ NumberOfSignatures
Definition GeometrySignature.h:32
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
Trk::cache
const IIntersectionCache * cache() const
Retrieve the associated cache block, if it exists.
Trk::active
@ active
Definition Layer.h:47
index
Definition index.py:1
Trk::BoundaryTrackParameters::trPar
const TrackParameters * trPar
Definition TimedExtrapolator.h:75
Trk::PathLimit::updateMat
void updateMat(float dX0, float Z, float dL0)
collected material update
Definition HelperStructs.h:51
Generated on for ATLAS Offline Software by doxygen 1.16.1