9 #include "GaudiKernel/MsgStream.h"
10 #include "GaudiKernel/PhysicalConstants.h"
60 declareInterface<ITimedExtrapolator>(
this);
72 if (m_propagators.empty()) {
73 m_propagators.push_back(
"Trk::RungeKuttaPropagator/DefaultPropagator");
75 if (m_updators.empty()) {
76 m_updators.push_back(
"Trk::MaterialEffectsUpdator/DefaultMaterialEffectsUpdator");
78 if (m_msupdators.empty()) {
79 m_msupdators.push_back(
"Trk::MultipleScatteringUpdator/AtlasMultipleScatteringUpdator");
83 if (!m_propagators.empty()) {
84 if (m_propagators.retrieve().isFailure()) {
86 return StatusCode::FAILURE;
94 unsigned int validprop = m_propagators.size();
98 ATH_MSG_WARNING(
" Extrapolators jumps back in unconfigured mode, only strategy pattern methods can be used.");
100 m_configurationLevel = validprop - 1;
101 ATH_MSG_VERBOSE(
"Configuration level automatically set to " << m_configurationLevel);
105 if (m_navigator.retrieve().isFailure()) {
107 return StatusCode::FAILURE;
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;
123 unsigned int validmeuts = m_updators.size();
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));
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));
142 m_propNames.value().push_back(m_propNames[0]);
145 m_updatNames.value().push_back(m_updatNames[0]);
147 if (validprop && validmeuts) {
150 unsigned int index = 0;
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) {
159 m_subPropagators[isign] = (
index < validprop) ? &(*m_propagators[
index]) : &(*m_propagators[
Trk::Global]);
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) {
173 <<
" -- At least one IPropagator and IMaterialUpdator instance have to be given.! ");
177 m_maxNavigSurf = 1000;
182 return StatusCode::SUCCESS;
189 return StatusCode::SUCCESS;
192 std::unique_ptr<const Trk::TrackParameters>
198 std::vector<Trk::HitInfo> * &hitInfo,
210 "M-[" << ++cache.
m_methodSequence <<
"] extrapolateWithPathLimit(...): resolve active layers? " << m_resolveActive);
212 if (!m_stepPropagator) {
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");
232 if (boundaryVol && !boundaryVol->
inside(parm.
position(), m_tolerance)) {
237 std::unique_ptr<const Trk::TrackParameters> returnParms =
238 extrapolateToVolumeWithPathLimit(
239 cache, parm, timeLim,
dir,
particle, nextGeoID, boundaryVol);
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() <<
","
257 for (; garbageIter != garbageEnd; ++garbageIter)
if (garbageIter->first) {
258 if(garbageIter->first == returnParms.get()) {
261 << positionOutput(garbageIter->first->position())
262 <<
" parm=" << garbageIter->first
263 <<
" is the return param. Cloning to" << ret.get());
264 returnParms = std::move(ret);
271 std::unique_ptr<const Trk::TrackParameters>
285 std::unique_ptr<const Trk::TrackParameters> returnParameters =
nullptr;
289 std::vector<unsigned int> solutions;
291 unsigned int iDest = 0;
292 const EventContext& ctx = Gaudi::Hive::currentContext();
293 ATH_MSG_DEBUG(
" [+] start extrapolateToVolumeWithPathLimit - at " << positionOutput(parm.
position())<<
" parm="<<&parm);
295 if (destVol && m_navigator->atVolumeBoundary(currPar, destVol,
dir, nextVol, m_tolerance) && nextVol != destVol) {
303 emptyGarbageBin(cache,&parm);
313 if (!tgVol || tgVol != destVol) {
315 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
319 iDest = bounds.size();
324 bool updateStatic =
false;
328 cache.
m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
354 return returnParameters;
356 throwIntoGarbageBin(cache,aPar);
358 return extrapolateToVolumeWithPathLimit(cache,*aPar, timeLim,
dir,
particle, nextGeoID, destVol);
366 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
398 if (!detVols.empty()) {
400 for (; iTer != detVols.end(); ++iTer) {
402 const Trk::Layer *layR = (*iTer)->layerRepresentation();
404 const auto& detBounds = (*iTer)->trackingVolume()->boundarySurfaces();
408 for (
unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
409 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
413 !m_useMuonMatApprox ||
414 (*iTer)->name().substr((*iTer)->name().size() - 4, 4) ==
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();
431 if (!(*iTer)->trackingVolume()->confinedDenseVolumes().empty() || (confLays.size() > detBounds.size())) {
433 for (
unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
434 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
437 }
else if (!confLays.empty()) {
438 for (
const Trk::Layer*
const lIt : confLays) {
439 cache.
m_layers.emplace_back(&(lIt->surfaceRepresentation()),
441 cache.
m_navigLays.emplace_back((*iTer)->trackingVolume(), lIt);
477 std::vector<std::pair<const Trk::TrackingVolume *, unsigned int> > navigVols;
480 std::vector<const Trk::DetachedTrackingVolume *> detVols =
481 m_navigator->trackingGeometry(ctx)->lowestDetachedTrackingVolumes(gp);
483 for (; dIter != detVols.end(); ++dIter) {
484 const Trk::Layer *layR = (*dIter)->layerRepresentation();
486 if (
active && !m_resolveActive) {
490 m_useMuonMatApprox && (*dIter)->name().substr((*dIter)->name().size() - 4, 4) !=
"PERM") {
495 bool dExit = m_navigator->atVolumeBoundary(currPar, dVol,
dir, nextVol, m_tolerance) && !nextVol;
503 if (!
active && confinedDense.empty() && confinedLays.empty()) {
507 if (!
active && confinedDense.empty() && confinedLays.size() <= bounds.size()) {
510 if (!confinedDense.empty() || !confinedLays.empty()) {
511 navigVols.emplace_back(dVol, bounds.size());
512 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
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++) {
529 if (!confinedLays.empty()) {
530 for (
const auto *confinedLay : confinedLays) {
531 cache.
m_layers.emplace_back(&(confinedLay->surfaceRepresentation()),
true);
539 for (
const auto *subvol : subvols) {
540 if (subvol->inside(gp, m_tolerance)) {
550 bool vExit = m_navigator->atVolumeBoundary(currPar, detVol,
dir, nextVol, m_tolerance) && nextVol != detVol;
551 if (vExit && nextVol && nextVol->
inside(gp, m_tolerance)) {
557 navigVols.emplace_back(detVol, bounds.size());
558 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
563 cache.
m_denseVols.emplace_back(detVol, bounds.size());
564 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
573 for (
const auto *cLay : cLays) {
574 if (cLay->layerType() > 0 || cLay->layerMaterialProperties()) {
575 cache.
m_layers.emplace_back(&(cLay->surfaceRepresentation()),
true);
591 if (nextLayer && nextLayer != lay) {
599 cache.
m_layers.emplace_back(&(
layer->surfaceRepresentation()),
true);
612 for (
const auto *cLay : cLays) {
613 if (cLay->layerType() > 0 || cLay->layerMaterialProperties()) {
614 cache.
m_layers.emplace_back(&(cLay->surfaceRepresentation()),
631 if (nextLayer && nextLayer != lay) {
638 if (backLayer && backLayer != lay) {
673 if (!m_navigator->atVolumeBoundary(currPar, dVol,
dir, nextVol, m_tolerance) || nextVol == dVol) {
687 overlapSearch(cache,*m_subPropagators[0], *currPar, *currPar, *cache.
m_navigLays[
i].second, timeLim.
time,
dir,
true,
690 ATH_MSG_VERBOSE(
" [o] Collecting intersection with active input layer.");
703 std::vector<unsigned int> solutions;
707 <<
" (current momentum: " << currPar->
momentum().mag() <<
717 || m_navigator->atVolumeBoundary(currPar, cache.
m_currentDense,
dir, assocVol, m_tolerance))) {
736 ATH_MSG_DEBUG(
" [+] Position after propagation - at " << positionOutput(
743 return returnParameters;
746 throwIntoGarbageBin(cache,nextPar);
764 return returnParameters;
767 throwIntoGarbageBin(cache,iPar);
768 return extrapolateToVolumeWithPathLimit(cache,*iPar, timeLim,
dir,
particle, nextGeoID, destVol);
770 return returnParameters;
774 if (timeLim.
tMax > 0. && timeLim.
time >= timeLim.
tMax) {
783 return returnParameters;
785 throwIntoGarbageBin(cache,iPar);
786 return extrapolateToVolumeWithPathLimit(cache,*iPar, timeLim,
dir,
particle, nextGeoID, destVol);
788 return returnParameters;
794 || m_navigator->atVolumeBoundary(nextPar, cache.
m_currentDense,
dir, assocVol, m_tolerance))) {
799 ATH_MSG_DEBUG(
" [+] Number of intersection solutions: " << solutions.size());
801 unsigned int iSol = 0;
802 while (iSol < solutions.size()) {
803 if (solutions[iSol] < iDest) {
808 if (
mb && m_includeMaterialEffects) {
809 if (
mb->layerMaterialProperties() &&
mb->layerMaterialProperties()->fullMaterial(nextPar->
position())) {
811 nextPar = currentUpdator ? currentUpdator
823 ATH_MSG_VERBOSE(
" [+] Update may have killed neutral track - return.");
825 return returnParameters;
827 throwIntoGarbageBin(cache,nextPar);
834 unsigned int index = solutions[iSol] - iDest;
844 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
855 m_tolerance) && assocVol != nextVol) {
860 ATH_MSG_DEBUG(
" [+] World boundary reached - at " << positionOutput(
876 return extrapolateToVolumeWithPathLimit(cache,*nextPar, timeLim,
dir,
particle, nextGeoID, destVol);
890 double pIn = nextPar->
momentum().mag();
891 nextPar = currentUpdator ? currentUpdator->
update(nextPar, *nextLayer, timeLim, cache.
m_path,
897 return returnParameters;
900 " Layer energy loss:" << nextPar->
momentum().mag() - pIn <<
"at position:" << nextPar->
position() <<
", current momentum:" <<
902 throwIntoGarbageBin(cache,nextPar);
910 overlapSearch(cache,*m_subPropagators[0], *currPar, *nextPar, *nextLayer, timeLim.
time,
dir,
true,
particle);
922 if (newLayer && newLayer != nextLayer) {
940 return extrapolateToVolumeWithPathLimit(cache,*nextPar, timeLim,
dir,
particle, nextGeoID, destVol);
950 std::vector< std::pair<const Trk::TrackingVolume *, unsigned int> >
::iterator dIter = cache.
m_denseVols.begin();
952 index -= (*dIter).second;
956 currVol = (*dIter).first;
960 if (!nextVol || !nextVol->
inside(
tp, m_tolerance)) {
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;
989 if (nIter != navigVols.end()) {
990 currVol = (*nIter).first;
991 nextVol = ((*nIter).first->boundarySurfaces())[
index]->attachedVolume(*nextPar,
dir);
993 ATH_MSG_DEBUG(
" [+] Navigation volume boundary, leaving volume '"
1002 return extrapolateToVolumeWithPathLimit(cache,*currPar, timeLim,
dir,
particle, nextGeoID, destVol);
1013 index -= (*dIter).second;
1017 currVol = (*dIter).first->trackingVolume();
1019 ((*dIter).first->trackingVolume()->boundarySurfaces())[
index]->attachedVolume(*nextPar,
dir);
1028 return extrapolateToVolumeWithPathLimit(cache, *currPar, timeLim,
dir,
particle, nextGeoID, destVol);
1037 return returnParameters;
1051 bool startingLayer)
const {
1053 const EventContext& ctx = Gaudi::Hive::currentContext();
1055 bool isDestinationLayer =
false;
1065 ATH_MSG_VERBOSE(
" [o] OverlapSearch called " << (startSurface ?
"with " :
"w/o ") <<
"start, "
1066 << (endSurface ?
"with " :
"w/o ") <<
"end surface.");
1074 ATH_MSG_VERBOSE(
" [o] Detector surface found through subSurface() call");
1079 ATH_MSG_VERBOSE(
" [o] Detector surface found through parameter on layer association");
1083 bool isStartLayer = (detSurface && detSurface == startSurface);
1087 std::vector<const Trk::TrackParameters*> detParametersOnLayer;
1088 bool reorderDetParametersOnLayer =
false;
1092 if (isDestinationLayer) {
1093 detParameters = (&parsOnLayer);
1094 }
else if (isStartLayer) {
1095 detParameters = (&parm);
1096 }
else if (detSurface) {
1098 detParameters = prop.
propagate(ctx,parm, *detSurface,
dir,
false, m_fieldProperties,
particle).release();
1102 bool surfaceHit =
true;
1103 if (detParameters &&
1105 !isDestinationLayer) {
1106 ATH_MSG_VERBOSE(
" [o] First intersection with Detector surface: " << *detParameters);
1108 surfaceHit = detParameters && detSurface ? detSurface->isOnSurface(detParameters->
position()) : 0;
1114 surfaceHit = (surfaceHit && startSurface) ?
1117 surfaceHit = (surfaceHit && endSurface) ?
1128 detSurface != startSurface) {
1131 detParametersOnLayer.push_back(detParameters);
1132 }
else if (detParameters) {
1134 ATH_MSG_VERBOSE(
" [-] Detector surface hit cancelled through bounds check or start/end surface check.");
1135 throwIntoGarbageBin(cache,detParameters);
1140 if (detParameters) {
1142 std::vector<Trk::SurfaceIntersection> cSurfaces;
1147 ATH_MSG_VERBOSE(
"found " << ncSurfaces <<
" candidate sensitive surfaces to test.");
1150 for (
auto &csf : cSurfaces) {
1161 if (overlapParameters) {
1162 ATH_MSG_VERBOSE(
" [+] Overlap surface was hit, checking start/end surface condition.");
1164 surfaceHit = (startSurface) ?
1167 surfaceHit = (surfaceHit && endSurface) ?
1169 parsOnLayer.
momentum().normalized())
1171 if (surfaceHit && csf.object!=detSurface) {
1174 reorderDetParametersOnLayer =
true;
1176 detParametersOnLayer.push_back(overlapParameters);
1179 ATH_MSG_VERBOSE(
" [-] Detector surface hit cancelled through start/end surface check.");
1180 throwIntoGarbageBin(cache,overlapParameters);
1188 std::vector<const Trk::TrackParameters *>::const_iterator parsOnLayerIter = detParametersOnLayer.begin();
1189 std::vector<const Trk::TrackParameters *>::const_iterator parsOnLayerIterEnd = detParametersOnLayer.end();
1192 if (reorderDetParametersOnLayer) {
1195 sort(detParametersOnLayer.begin(), detParametersOnLayer.end(), parameterSorter);
1199 parsOnLayerIter = detParametersOnLayer.begin();
1200 parsOnLayerIterEnd = detParametersOnLayer.end();
1202 for (; parsOnLayerIter != parsOnLayerIterEnd; ++parsOnLayerIter) {
1205 std::unique_ptr<const Trk::TrackParameters>(*parsOnLayerIter),
1215 std::stringstream outStream;
1217 if (m_printRzOutput) {
1218 outStream <<
"[r,phi,z] = [ " <<
pos.perp() <<
", " <<
pos.phi() <<
", " <<
pos.z() <<
" ]";
1220 outStream <<
"[xyz] = [ " <<
pos.x() <<
", " <<
pos.y() <<
", " <<
pos.z() <<
" ]";
1222 return outStream.str();
1227 std::stringstream outStream;
1229 outStream <<
"[eta,phi] = [ " <<
mom.eta() <<
", " <<
mom.phi() <<
" ]";
1230 return outStream.str();
1240 bool throwCurrent =
false;
1242 for (; garbageIter != garbageEnd; ++garbageIter) {
1243 if (garbageIter->first && garbageIter->first != trPar) {
1244 delete (garbageIter->first);
1246 if (garbageIter->first && garbageIter->first == trPar) {
1247 throwCurrent =
true;
1253 throwIntoGarbageBin(cache,trPar);
1261 for (
const auto *subUpdator : m_subUpdators) {
1262 subUpdator->validationAction();
1267 std::unique_ptr<const Trk::TrackParameters>
1272 std::vector<Trk::HitInfo> * &hitInfo,
1284 "M-[" << ++cache.
m_methodSequence <<
"] transportNeutralsWithPathLimit(...) " << pathLim.
x0Max <<
", from " <<
1301 if (boundaryVol && !boundaryVol->
inside(parm.
position(), m_tolerance)) {
1308 std::unique_ptr<const Trk::TrackParameters> returnParms =
1309 transportToVolumeWithPathLimit(
1310 cache, parm, timeLim,
dir,
particle, nextGeoID, boundaryVol);
1322 for (; garbageIter != garbageEnd; ++garbageIter)
if (garbageIter->first) {
1323 if(garbageIter->first == returnParms.get()) {
1326 << positionOutput(garbageIter->first->position())
1327 <<
" parm=" << garbageIter->first
1328 <<
" is the return param. Cloning to" << ret.get());
1329 returnParms=std::move(ret);
1336 std::unique_ptr<const Trk::TrackParameters>
1351 std::unique_ptr<const Trk::TrackParameters> returnParameters =
nullptr;
1357 unsigned int iDest = 0;
1359 const EventContext& ctx = Gaudi::Hive::currentContext();
1361 if (destVol && m_navigator->atVolumeBoundary(currPar, destVol,
dir, nextVol, m_tolerance) && nextVol != destVol) {
1370 emptyGarbageBin(cache,&parm);
1372 if (cache.
m_trSurfs.capacity() > m_maxNavigSurf) {
1373 cache.
m_trSurfs.reserve(m_maxNavigSurf);
1380 if (!tgVol || tgVol != destVol) {
1382 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
1383 const Trk::Surface &surf = (bounds[
ib])->surfaceRepresentation();
1412 cache.
m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
1439 return returnParameters;
1441 throwIntoGarbageBin(cache,aPar);
1442 return transportToVolumeWithPathLimit(cache,*aPar, timeLim,
dir,
particle, nextGeoID, destVol);
1450 for (
unsigned int ib = 0;
ib < bounds.size(); ++
ib) {
1451 const Trk::Surface &surf = (bounds[
ib])->surfaceRepresentation();
1456 distSol.
first() > m_tolerance) {
1457 double dist = distSol.
first();
1464 if (surf.
isOnSurface(gp,
true, m_tolerance, m_tolerance)) {
1469 double dist = distSol.
second();
1472 if (surf.
isOnSurface(gp,
true, m_tolerance, m_tolerance)) {
1473 if (dist > m_tolerance) {
1482 " transportToVolumeWithPathLimit() - at " << currPar->
position() <<
", missing static volume boundary "
1484 ": transport interrupted");
1487 "---> particle R,phi,z, momentum:" << currPar->
position().perp() <<
"," << currPar->
position().phi() <<
"," << currPar->
position().z() <<
"," <<
1499 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
1500 const Trk::Surface &surf = (bounds[
ib])->surfaceRepresentation();
1504 "---> decomposed boundary surface position, normal, estimated distance:" <<
ib <<
"," << surf.
center() <<
"," <<
1507 "---> estimated distance to (first solution):boundary check:" << distSol.
numberOfSolutions() <<
"," << distSol.
first() <<
":" <<
1509 m_tolerance, m_tolerance));
1511 ATH_MSG_DEBUG(
"---> estimated distance to (second solution):boundary check:" << distSol.
second() <<
"," <<
1513 m_tolerance, m_tolerance));
1517 return returnParameters;
1523 cache.
m_currentStatic = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
1526 return transportToVolumeWithPathLimit(cache,parm, timeLim,
dir,
particle, nextGeoID, destVol);
1528 ATH_MSG_DEBUG(
" [+] World boundary reached - at " << positionOutput(
1544 if (!detVols.empty()) {
1546 for (; iTer != detVols.end(); ++iTer) {
1548 const Trk::Layer *layR = (*iTer)->layerRepresentation();
1552 if (!m_resolveMultilayers || (*iTer)->multilayerRepresentation().empty()) {
1561 cache.
m_navigLays.emplace_back((*iTer)->trackingVolume(), layR);
1565 const auto&
multi = (*iTer)->multilayerRepresentation();
1566 for (
const auto *
i :
multi) {
1575 cache.
m_navigLays.emplace_back((*iTer)->trackingVolume(),
i);
1582 (*iTer)->name().substr((*iTer)->name().size() - 4, 4) ==
"PERM") {
1585 if ((*iTer)->trackingVolume()->zOverAtimesRho() != 0. &&
1586 ((*iTer)->trackingVolume()->confinedDenseVolumes().empty())
1587 && ((*iTer)->trackingVolume()->confinedArbitraryLayers().empty())) {
1588 const auto& detBounds = (*iTer)->trackingVolume()->boundarySurfaces();
1590 for (
unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
1591 const Trk::Surface &surf = (detBounds[ibb])->surfaceRepresentation();
1604 cache.
m_denseVols.emplace_back((*iTer)->trackingVolume(), newB);
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();
1621 for (
unsigned int ibb = 0; ibb < bounds.size(); ibb++) {
1622 const Trk::Surface &surf = (bounds[ibb])->surfaceRepresentation();
1637 if (!(*vIter)->confinedArbitraryLayers().empty()) {
1639 " transportToVolumeWithPathLimit() - at " << currPar->
position() <<
", unresolved sublayers/subvolumes for "
1640 << (*vIter)->volumeName());
1647 if (!confLays.empty()) {
1648 for (
const Trk::Layer*
const lIt: confLays) {
1649 const Trk::Surface &surf = lIt->surfaceRepresentation();
1657 cache.
m_navigLays.emplace_back((*iTer)->trackingVolume(), lIt);
1670 cache.
m_trSurfs.emplace_back((*bIter).surface, (*bIter).distance);
1699 for (
const auto *cLay : cLays) {
1700 if (cLay->layerMaterialProperties()) {
1701 const Trk::Surface &surf = cLay->surfaceRepresentation();
1733 if (!m_navigator->atVolumeBoundary(currPar, dVol,
dir, nextVol, m_tolerance) ||
1747 std::vector<unsigned int> sols;
1749 for (
unsigned int i = 0;
i < cache.
m_trSurfs.size(); ++
i) {
1753 if (sols.size() > 1) {
1754 unsigned int itest = 1;
1755 while (itest < sols.size()) {
1757 unsigned int iex = sols[itest - 1];
1758 sols[itest - 1] = sols[itest];
1766 for (
unsigned int is = 1; is < sols.size(); is++) {
1768 std::cout <<
"wrong intersection ordering" << std::endl;
1787 for (
unsigned int sol : sols) {
1788 if (cache.
m_trSurfs[sol].second == 0.) {
1840 double fr = fmin(frT, frM);
1865 if (nextPar &&
process == 121) {
1866 ATH_MSG_DEBUG(
" [!] WARNING: failed hadronic interaction, killing the input particle anyway");
1867 return returnParameters;
1871 return returnParameters;
1874 throwIntoGarbageBin(cache,nextPar);
1877 return returnParameters;
1889 throwIntoGarbageBin(cache,nextPar);
1896 if (
mb && m_includeMaterialEffects) {
1897 if (
mb->layerMaterialProperties() &&
mb->layerMaterialProperties()->fullMaterial(nextPos)) {
1908 ATH_MSG_VERBOSE(
" [+] Update may have killed neutral track - return.");
1910 return returnParameters;
1912 throwIntoGarbageBin(cache,nextPar);
1929 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
1944 ATH_MSG_DEBUG(
" [+] World boundary reached - at " << positionOutput(
1959 return transportToVolumeWithPathLimit(cache,*nextPar, timeLim,
dir,
particle, nextGeoID, destVol);
1962 return transportToVolumeWithPathLimit(cache,*nextPar, timeLim,
dir,
particle, nextGeoID, destVol);
1972 if (matUp && m_includeMaterialEffects) {
1975 nextPar = currentUpdator ? currentUpdator
1987 ATH_MSG_VERBOSE(
" [+] Update may have killed neutral track - return.");
1989 return returnParameters;
1991 throwIntoGarbageBin(cache,nextPar);
1998 std::vector< std::pair<const Trk::TrackingVolume *, unsigned int> >
::iterator dIter = cache.
m_denseVols.begin();
2000 index -= (*dIter).second;
2004 currVol = (*dIter).first;
2010 }
else if (currVol->
inside(nextPos + 0.002 *
dir * nextPar->
momentum().normalized())) {
2015 if (m_useMuonMatApprox && cache.
m_denseVols.empty()) {
2032 throwIntoGarbageBin(cache,nextPar);
2054 const std::string
m = aliTV ? aliTV->
volumeName() :
" NULLPTR!";
2055 ATH_MSG_DEBUG(
" [0] starting transport of neutral particle in alignable volume " <<
m);
2064 std::vector<Trk::IdentifiedIntersection> iis;
2066 emptyGarbageBin(cache,&parm);
2068 const EventContext& ctx = Gaudi::Hive::currentContext();
2070 return {
nullptr,
nullptr,
nullptr};
2096 if (binIDMat->second > 0) {
2103 unsigned int cbin = lbu->
bin(
pos);
2111 if (d2n.first == cbin) {
2122 if (d2n.first == cbin && fabs(d2n.second) < 0.002) {
2128 if (d2n.second > 0.001) {
2129 pot =
pos + 0.5 * d2n.second *
dir * umo;
2131 iis.emplace_back(distTot, binIDMat->second, binIDMat->first.get());
2142 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
2143 const Trk::Surface &surf = (bounds[
ib])->surfaceRepresentation();
2146 double dist = distSol.
first();
2155 if (dist > m_tolerance && surf.
isOnSurface(gp,
true, m_tolerance, m_tolerance)) {
2158 if (attachedVol && !(attachedVol->
inside(gp + 0.01 *
dir * currPar->
momentum().normalized(), m_tolerance))) {
2162 attachedVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
2163 gp + 0.01 *
dir * currPar->
momentum().normalized());
2170 nextVol = attachedVol;
2172 }
else if (dist > 0.001) {
2177 "gluing problem at the exit from alignable volume: " << gp.perp() <<
"," << gp.z() <<
":" <<
2186 "next volume resolved to:" << testVol->
volumeName() <<
" at the position(R,Z):" << gp.perp() <<
"," <<
2200 return {
nullptr,
nullptr,
nullptr};
2204 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
2209 if (!iis.empty() && cache.
m_trStaticBounds[0].distance - iis.back().distance < 0.01) {
2225 for (
unsigned int is = 0; is < iis.size(); is++) {
2230 double step = iis[is].distance - dist;
2232 nextPos = currPar->
position() +
dir * currPar->
momentum().normalized() * iis[is].distance;
2251 double mDeltaL = currMat->
L0 > 0. ?
step / currMat->
L0 : mDelta / 0.37 / currMat->
averageZ();
2258 double fr = fmin(frT, frM);
2266 if (mDelta > 0 && currMat->
averageZ() > 0) {
2273 if (m_caloMsSecondary) {
2278 throwIntoGarbageBin(cache, nextPar);
2284 if (nextPar &&
process == 121) {
2285 ATH_MSG_DEBUG(
" [!] WARNING: failed hadronic interaction, killing the input particle anyway");
2286 return {
nullptr,
nullptr,
nullptr};
2290 return {
nullptr,
nullptr,
nullptr};
2295 return {
nullptr,
nullptr,
nullptr};
2300 dist = iis[is].distance;
2301 if (mDelta > 0 && currMat->
averageZ() > 0) {
2306 if (is < iis.size() - 1) {
2309 currMat = iis[is].material;
2310 currLay = iis[is].identifier;
2312 if (cache.
m_hitVector && iis[is].identifier > 0) {
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.);
2323 ATH_MSG_VERBOSE(
"active layer/volume exit:" << currLay <<
" at R,z:" << nextPos.perp() <<
"," << nextPos.z());
2324 if (binIDMat and(binIDMat->second > 0)) {
2329 throwIntoGarbageBin(cache,nextPar);
2337 ATH_MSG_DEBUG(
" [+] World boundary reached - at " << positionOutput(
2358 const std::string
m = vol ? vol->
volumeName():
"NULLPTR";
2373 std::vector<unsigned int> solutions;
2376 const EventContext& ctx = Gaudi::Hive::currentContext();
2381 emptyGarbageBin(cache,&parm);
2385 if (vol && vol->
inside(gp, m_tolerance)) {
2388 currVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(gp);
2390 if (m_navigator->atVolumeBoundary(currPar, currVol,
dir, nextStatVol, m_tolerance) && nextStatVol != currVol) {
2391 currVol = nextStatVol;
2393 if (currVol && currVol != vol) {
2402 ATH_MSG_DEBUG(
" [!] failing in retrieval of AlignableTV, return 0");
2403 return {
nullptr,
nullptr,
nullptr};
2414 if (binIDMat->second > 0) {
2430 for (
unsigned int ib = 0;
ib < bounds.size(); ++
ib) {
2431 const Trk::Surface &surf = (bounds[
ib])->surfaceRepresentation();
2445 std::vector<unsigned int> solutions;
2448 <<
" (current momentum: " << currPar->
momentum().mag() <<
2475 ATH_MSG_DEBUG(
" [+] Number of intersection solutions: " << solutions.size());
2478 throwIntoGarbageBin(cache,nextPar);
2500 return {
nullptr,
nullptr,
nullptr};
2503 throwIntoGarbageBin(cache,iPar);
2506 ATH_MSG_DEBUG(
" [!] WARNING: failed hadronic interaction, killing the input particle anyway");
2507 return {
nullptr,
nullptr,
nullptr};
2512 return {
nullptr,
nullptr,
nullptr};
2517 if (timeLim.
tMax > 0. && timeLim.
time >= timeLim.
tMax) {
2524 return {
nullptr,
nullptr,
nullptr};
2527 throwIntoGarbageBin(cache,iPar);
2529 return {
nullptr,
nullptr,
nullptr};
2531 return {
nullptr,
nullptr,
nullptr};
2536 unsigned int iSol = 0;
2537 while (iSol < solutions.size()) {
2541 unsigned int index = solutions[iSol];
2550 nextVol = m_navigator->trackingGeometry(ctx)->lowestStaticTrackingVolume(
2563 if (binIDMat && binIDMat->second > 0 && !cache.
m_hitVector->empty() &&
2564 cache.
m_hitVector->back().detID == binIDMat->second) {
2577 ATH_MSG_DEBUG(
" [+] World boundary reached - at " << positionOutput(
2593 return {
nullptr,
nullptr,
nullptr};