736{
737 RecursionCounter
counter(cache.m_recursionCount);
739 ATH_MSG_WARNING(
"Too many recursive calls of extrapolateToNextMaterialLayer: "
742 return {};
743 }
744 ++cache.m_methodSequence;
745 ATH_MSG_DEBUG(
"M-[" << cache.m_methodSequence <<
"] extrapolateToNextMaterialLayer(...) ");
747 ATH_MSG_WARNING(
"Too many method sequence calls of extrapolateToNextMaterialLayer: "
748 << cache.m_methodSequence);
750 return {};
751 }
752
753
754
755
756
757
758
759
760
761
763 const Trk::TrackingVolume* staticVol = nullptr;
764 const Trk::TrackingVolume* currVol = nullptr;
765 const Trk::TrackingVolume* nextVol = nullptr;
766 std::vector<unsigned int> solutions;
767 const Trk::TrackingVolume* assocVol = nullptr;
768
770 bool resolveActive = destSurf == nullptr;
773 }
774 if (cache.m_lastMaterialLayer && !cache.m_lastMaterialLayer->isOnLayer(parm->position())) {
775 cache.m_lastMaterialLayer = nullptr;
776 }
777
779 if (!cache.m_highestVolume) {
780 cache.m_highestVolume = cache.m_trackingGeometry->highestTrackingVolume();
781 }
782
785 staticVol = vol;
786 } else {
787 staticVol = cache.m_trackingGeometry->lowestStaticTrackingVolume(gp);
788 const Trk::TrackingVolume* nextStatVol = nullptr;
790 nextStatVol != staticVol) {
791 staticVol = nextStatVol;
792 }
793 }
794
795
798 }
799 cache.m_navigSurfs.clear();
800 if (destSurf) {
801 cache.m_navigSurfs.emplace_back(destSurf, false);
802 }
803
806 const Trk::AlignableTrackingVolume* alignTV =
807 static_cast<const Trk::AlignableTrackingVolume*>(staticVol);
808 cache.m_identifiedParameters.reset();
810 alignTV, dir, particle);
811 }
812 }
813
814
815 if (staticVol && (staticVol != cache.m_currentStatic || resolveActive !=
m_resolveActive)) {
816
817 cache.m_currentStatic = staticVol;
818 cache.retrieveBoundaries();
819
820 cache.m_detachedVols.clear();
821 cache.m_detachedBoundaries.clear();
822 cache.m_denseVols.clear();
823 cache.m_denseBoundaries.clear();
824 cache.m_layers.clear();
825 cache.m_navigLays.clear();
826
829 if (!detVols.empty()) {
830 Trk::ArraySpan<const Trk::DetachedTrackingVolume* const>::iterator iTer = detVols.begin();
831 for (; iTer != detVols.end(); ++iTer) {
832
833 const Trk::Layer* layR = (*iTer)->layerRepresentation();
835 const auto detBounds = (*iTer)->trackingVolume()->boundarySurfaces();
837 if (resolveActive) {
838 cache.m_detachedVols.emplace_back(*iTer, detBounds.size());
839 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
840 const Trk::Surface& surf = (detBounds[ibb])->surfaceRepresentation();
841 cache.m_detachedBoundaries.emplace_back(&surf, true);
842 }
843 } else {
845 cache.addOneNavigationLayer((*iTer)->trackingVolume(), layR);
846
847 } else {
848 const auto&
multi = (*iTer)->multilayerRepresentation();
849 for (const auto *i : multi) {
850 cache.addOneNavigationLayer((*iTer)->trackingVolume(), i);
851 }
852 }
853 }
856 (*iTer)->name().compare((*iTer)->name().size() - 4, 4, "PERM") ==
857 0) {
858
859
860
861
862
863 if ((*iTer)->trackingVolume()->zOverAtimesRho() != 0. &&
864 ((*iTer)->trackingVolume()->confinedDenseVolumes().empty()) &&
865 ((*iTer)->trackingVolume()->confinedArbitraryLayers().empty())) {
866 cache.m_denseVols.emplace_back((*iTer)->trackingVolume(), detBounds.size());
867
868 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
869 const Trk::Surface& surf = (detBounds[ibb])->surfaceRepresentation();
870 cache.m_denseBoundaries.emplace_back(&surf, true);
871 }
872 }
874 (*iTer)->trackingVolume()->confinedArbitraryLayers();
875 if (!(*iTer)->trackingVolume()->confinedDenseVolumes().empty() ||
876 (confLays.size() > detBounds.size())) {
877 cache.m_detachedVols.emplace_back(*iTer, detBounds.size());
878 for (unsigned int ibb = 0; ibb < detBounds.size(); ibb++) {
879 const Trk::Surface& surf =
880 (detBounds[ibb])->surfaceRepresentation();
881 cache.m_detachedBoundaries.emplace_back(&surf, true);
882 }
883 } else if (!confLays.empty()) {
884 for (const Trk::Layer* const lIt : confLays) {
885 cache.addOneNavigationLayer((*iTer)->trackingVolume(), (lIt));
886 }
887 }
888 }
889 }
890 }
891 cache.m_denseResolved = std::pair<unsigned int, unsigned int>(cache.m_denseVols.size(),
892 cache.m_denseBoundaries.size());
893 cache.m_layerResolved = cache.m_layers.size();
894 }
895
896 cache.m_navigSurfs.insert(
897 cache.m_navigSurfs.end(), cache.m_staticBoundaries.begin(), cache.m_staticBoundaries.end());
898
899
900 if (staticVol) {
903 }
905
906 solutions.resize(0);
907 const Trk::TrackingVolume* propagVol = cache.m_dense ? staticVol : cache.m_highestVolume;
910 << "'");
911
912
913 std::unique_ptr<Trk::TrackParameters> pNextPar = prop.propagate(ctx,*currPar,cache.m_navigSurfs,
915 false,propagVol);
917 if (nextPar) {
918 ++cache.m_nPropagations;
921 }
922 if (!nextPar) {
923 cache.m_parametersAtBoundary.resetBoundaryInformation();
924 return {};
925 }
926 if (nextPar) {
927
928 if (propagVol->
zOverAtimesRho() != 0. && !cache.m_matstates && cache.m_extrapolationCache) {
929 if (not cache.elossPointerOverwritten()) {
931 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer"));
932 }
933 const double dInX0 = std::abs(path) / propagVol->
x0();
935 cache.m_extrapolationCache->updateX0(dInX0);
936 const Trk::MaterialProperties materialProperties(*propagVol, std::abs(path));
937 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
939 materialProperties, std::abs(1. / currentqoverp), 1., dir, particle);
941 << nextPar->momentum().mag() - currPar->momentum().mag() << ","
942 << eloss.deltaE());
943 cache.m_extrapolationCache->updateEloss(
944 eloss.meanIoni(), eloss.sigmaIoni(), eloss.meanRad(), eloss.sigmaRad());
947 }
948 } else {
950 }
951 }
953 const double dInX0 = std::abs(path) / propagVol->
x0();
954 const Trk::MaterialProperties materialProperties(*propagVol, std::abs(path));
955 const double scatsigma = std::sqrt(
m_msupdater->sigmaSquare(
956 materialProperties, 1. / std::abs(nextPar->parameters()[
qOverP]), 1., particle));
957 auto newsa = Trk::ScatteringAngles(
958 0, 0, scatsigma / std::sin(nextPar->parameters()[
Trk::theta]), scatsigma);
959
960 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
962 materialProperties, std::abs(1. / currentqoverp), 1., dir, particle);
963
965 << nextPar->momentum().mag() - currPar->momentum().mag() << ","
966 << eloss.deltaE());
967
969 nextPar->position(), nextPar->momentum(), nextPar->charge()));
970
971 if (cache.m_extrapolationCache) {
973 ATH_MSG_DEBUG(cache.to_string(
" mat states extrapolateToNextMaterialLayer"));
974 }
975 cache.m_extrapolationCache->updateX0(dInX0);
976 cache.m_extrapolationCache->updateEloss(
977 eloss.meanIoni(), eloss.sigmaIoni(), eloss.meanRad(), eloss.sigmaRad());
980 }
981 }
982 auto mefot = std::make_unique<Trk::MaterialEffectsOnTrack>(
983 dInX0, newsa, std::make_unique<Trk::EnergyLoss>(std::move(eloss)),
984 cvlTP->associatedSurface());
985 cache.m_matstates->push_back(new TrackStateOnSurface(
986 nullptr, std::move(cvlTP), std::move(mefot)));
988 << "', t/X0 = " << dInX0);
989 }
990 }
991 currPar = nextPar;
992 const unsigned int isurf = destSurf ? 1 : 0;
993 if (destSurf && solutions[0] == 0) {
994 return currPar;
995 }
996 if (destSurf && solutions.size() > 1 && solutions[1] == 0) {
997 return currPar;
998 }
999 if (solutions[0] <= isurf + cache.m_staticBoundaries.size()) {
1000
1001 const Trk::TrackingVolume* nextVol =
1002 cache.m_currentStatic->
boundarySurfaces()[solutions[0] - isurf]->attachedVolume(
1003 currPar->position(), currPar->momentum(), dir);
1004 cache.m_parametersAtBoundary.boundaryInformation(nextVol, currPar, currPar);
1005 if (!nextVol) {
1006 ATH_MSG_DEBUG(
" [!] World boundary at position R,z: " << currPar->position().perp() <<
","
1007 << currPar->position().z());
1008 } else {
1010 }
1011 }
1012 return {};
1013 }
1014
1016 return {};
1017 }
1018
1019
1020 cache.m_currentDense = cache.m_dense ? cache.m_currentStatic : cache.m_highestVolume;
1021 cache.m_navigBoundaries.clear();
1022 if (cache.m_denseVols.size() > cache.m_denseResolved.first) {
1023 cache.m_denseVols.resize(cache.m_denseResolved.first);
1024 }
1025 while (cache.m_denseBoundaries.size() > cache.m_denseResolved.second) {
1026 cache.m_denseBoundaries.pop_back();
1027 }
1028 if (cache.m_layers.size() > cache.m_layerResolved) {
1029 cache.m_navigLays.resize(cache.m_layerResolved);
1030 }
1031 while (cache.m_layers.size() > cache.m_layerResolved) {
1032 cache.m_layers.pop_back();
1033 }
1034
1035
1036
1038
1041 }
1042 cache.m_navigVolsInt.clear();
1043
1044 gp = currPar->position();
1045 std::vector<const Trk::DetachedTrackingVolume*> detVols =
1046 cache.m_trackingGeometry->lowestDetachedTrackingVolumes(gp);
1047 std::vector<const Trk::DetachedTrackingVolume*>::iterator dIter = detVols.begin();
1048 for (; dIter != detVols.end(); ++dIter) {
1049 const Trk::Layer* layR = (*dIter)->layerRepresentation();
1051 if (
active && !resolveActive) {
1052 continue;
1053 }
1055 (*dIter)->name().compare((*dIter)->name().size() - 4, 4, "PERM") != 0) {
1056 continue;
1057 }
1058 const Trk::TrackingVolume* dVol = (*dIter)->trackingVolume();
1059
1060 const bool dExit =
1062 if (dExit) {
1063 continue;
1064 }
1065
1068
1069 if (!
active && confinedDense.empty() && confinedLays.empty()) {
1070 continue;
1071 }
1073 if (!
active && confinedDense.empty() && confinedLays.size() <= bounds.size()) {
1074 continue;
1075 }
1076 if (!confinedDense.empty() || !confinedLays.empty()) {
1077 cache.m_navigVolsInt.emplace_back(dVol, bounds.size());
1078 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
1079 const Trk::Surface& surf = (bounds[
ib])->surfaceRepresentation();
1080 cache.m_navigBoundaries.emplace_back(&surf, true);
1081 }
1082
1083 if (!confinedDense.empty()) {
1084 auto vIter = confinedDense.begin();
1085 for (; vIter != confinedDense.end(); ++vIter) {
1086 const auto bounds = (*vIter)->boundarySurfaces();
1087 cache.m_denseVols.emplace_back(*vIter, bounds.size());
1088 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
1089 const Trk::Surface& surf = (bounds[
ib])->surfaceRepresentation();
1090 cache.m_denseBoundaries.emplace_back(&surf, true);
1091 }
1092 }
1093 }
1094
1095 if (!confinedLays.empty()) {
1096 for (const auto *confinedLay : confinedLays) {
1097 cache.addOneNavigationLayer(dVol, confinedLay);
1098 }
1099 }
1100 } else {
1104 for (const auto *subvol : subvols) {
1106 detVol = subvol;
1107 break;
1108 }
1109 }
1110 }
1111
1112 if (!detVol) {
1113 detVol = dVol;
1114 }
1115 bool vExit =
1117 nextVol != detVol;
1119 detVol = nextVol;
1120 vExit = false;
1121 }
1122 if (!vExit) {
1124 cache.m_navigVolsInt.emplace_back(detVol, bounds.size());
1125 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
1126 const Trk::Surface& surf = (bounds[
ib])->surfaceRepresentation();
1127 cache.m_navigBoundaries.emplace_back(&surf, true);
1128 }
1130 cache.m_denseVols.emplace_back(detVol, bounds.size());
1131 for (
unsigned int ib = 0;
ib < bounds.size();
ib++) {
1132 const Trk::Surface& surf = (bounds[
ib])->surfaceRepresentation();
1133 cache.m_denseBoundaries.emplace_back(&surf, true);
1134 }
1135 }
1136
1139
1140 if (lay) {
1141 cache.addOneNavigationLayer(detVol, lay);
1142
1143 }
1144 const Trk::Layer* nextLayer =
1145 detVol->
nextLayer(currPar->position(), dir * currPar->momentum().unit(),
true);
1146 if (nextLayer && nextLayer != lay) {
1147 cache.addOneNavigationLayer(detVol, nextLayer);
1148 }
1151 for (const auto *layer : layers) {
1152 cache.addOneNavigationLayer(detVol, layer);
1153 }
1154 }
1155 }
1156 }
1157 }
1158 cache.copyToNavigationSurfaces();
1159
1160
1161 cache.m_currentDense = cache.m_highestVolume;
1162 if (cache.m_dense && cache.m_denseVols.empty()) {
1163 cache.m_currentDense = cache.m_currentStatic;
1164 } else {
1165 for (
unsigned int i = 0;
i < cache.m_denseVols.size();
i++) {
1166 const Trk::TrackingVolume* dVol = cache.m_denseVols[
i].first;
1169 nextVol == dVol) {
1170 cache.m_currentDense = dVol;
1171 }
1172 }
1173 }
1174 }
1175
1176
1177
1178
1179 nextVol = nullptr;
1180 while (currPar) {
1182 std::vector<unsigned int> solutions;
1183
1185 currPar->position() + 2 *
m_tolerance *
dir * currPar->momentum().normalized();
1186 if (!(cache.m_currentDense->inside(tp, 0.))) {
1187 cache.m_currentDense = cache.m_highestVolume;
1188 if (cache.m_dense && cache.m_denseVols.empty()) {
1189 cache.m_currentDense = cache.m_currentStatic;
1190 } else {
1191 for (
unsigned int i = 0;
i < cache.m_denseVols.size();
i++) {
1192 const Trk::TrackingVolume* dVol = cache.m_denseVols[
i].first;
1194 cache.m_currentDense = dVol;
1195 }
1196 }
1197 }
1198 }
1199
1202 << " (current momentum: " << currPar->momentum().mag() << ")");
1203 ATH_MSG_DEBUG(
" [+] " << cache.m_navigSurfs.size() <<
" target surfaces in '"
1204 << cache.m_currentDense->volumeName() << "'.");
1207 solutions, path, false, false, cache.m_currentDense));
1208 if (nextPar) {
1209 ++cache.m_nPropagations;
1212 }
1213
1214 if (nextPar && !(cache.m_currentDense->inside(nextPar->position(),
m_tolerance) ||
1216 nextPar, cache.m_currentDense, dir, assocVol,
m_tolerance))) {
1217 ATH_MSG_DEBUG(
" [!] ERROR: missing volume boundary for volume"
1218 << cache.m_currentDense->volumeName());
1219 if (cache.m_currentDense->zOverAtimesRho() != 0.) {
1220 ATH_MSG_DEBUG(
" [!] ERROR: trying to recover: repeat the propagation step in"
1221 << cache.m_highestVolume->volumeName());
1222 cache.m_currentDense = cache.m_highestVolume;
1223 continue;
1224 }
1225 }
1226 if (nextPar) {
1227 ATH_MSG_DEBUG(
" [+] Number of intersection solutions: " << solutions.size());
1228 if (cache.m_currentDense->zOverAtimesRho() != 0. && !cache.m_matstates &&
1229 cache.m_extrapolationCache) {
1230 if (not cache.elossPointerOverwritten()) {
1232 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer dense "));
1233 }
1234 const double dInX0 = std::abs(path) / cache.m_currentDense->x0();
1235 cache.m_extrapolationCache->updateX0(dInX0);
1236 const Trk::MaterialProperties materialProperties(*cache.m_currentDense, std::abs(path));
1237 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
1239 materialProperties, std::abs(1. / currentqoverp), 1., dir, particle);
1240 cache.m_extrapolationCache->updateEloss(
1244 }
1245 } else {
1247 }
1248 }
1249
1250 if (cache.m_currentDense->zOverAtimesRho() != 0. && cache.m_matstates) {
1251 const double dInX0 = std::abs(path) / cache.m_currentDense->x0();
1252 if (path * dir < 0.) {
1254 }
1255 const Trk::MaterialProperties materialProperties(*cache.m_currentDense, std::abs(path));
1256 const double scatsigma = std::sqrt(
m_msupdater->sigmaSquare(
1257 materialProperties, 1. / std::abs(nextPar->parameters()[
qOverP]), 1., particle));
1258 auto newsa = Trk::ScatteringAngles(
1259 0, 0, scatsigma / std::sin(nextPar->parameters()[
Trk::theta]), scatsigma);
1260
1261 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
1263 materialProperties, std::abs(1. / currentqoverp), 1., dir, particle);
1264
1266 << nextPar->momentum().mag() - currPar->momentum().mag() << ","
1267 << eloss.deltaE());
1268
1269
1271 nextPar->position(), nextPar->momentum(), nextPar->charge()));
1272
1273 if (cache.m_extrapolationCache) {
1275 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer dense"));
1276 }
1277 cache.m_extrapolationCache->updateX0(dInX0);
1278 cache.m_extrapolationCache->updateEloss(eloss.meanIoni(),
1279 eloss.sigmaIoni(),
1280 eloss.meanRad(),
1281 eloss.sigmaRad());
1284 }
1285 }
1286 auto mefot = std::make_unique<Trk::MaterialEffectsOnTrack>(
1287 dInX0, newsa, std::make_unique<Trk::EnergyLoss>(std::move(eloss)), cvlTP->associatedSurface());
1288
1289 cache.m_matstates->push_back(new TrackStateOnSurface(
1290 nullptr, std::move(cvlTP), std::move(mefot)));
1291
1292 ATH_MSG_DEBUG(
" [M] Collecting material from dense volume '"
1293 << cache.m_currentDense->volumeName() << "', t/X0 = " << dInX0);
1294 }
1295
1296 if (destSurf && solutions[0] == 0) {
1297 return nextPar;
1298 }
1299 if (destSurf && solutions.size() > 1 && solutions[1] == 0) {
1300 return nextPar;
1301 }
1302
1303 if (destSurf) {
1304 double dist = 0.;
1306 nextPar->position(), nextPar->momentum().normalized());
1308 dist = distSol.
first();
1311 }
1314 }
1315 } else {
1317 }
1318 if (dist * dir < 0.) {
1319 ATH_MSG_DEBUG(
" [+] Destination surface missed ? " << dist <<
"," << dir);
1320 cache.m_parametersAtBoundary.resetBoundaryInformation();
1321 return {};
1322 }
1323 ATH_MSG_DEBUG(
" [+] New 3D-distance to destinatiion - d3 = " << dist * dir);
1324 }
1325
1326 int const iDest = destSurf ? 1 : 0;
1327 unsigned int iSol = 0;
1328 while (iSol < solutions.size()) {
1329 if (solutions[iSol] < iDest + cache.m_staticBoundaries.size()) {
1330
1331 const Trk::Layer*
mb = cache.m_navigSurfs[solutions[iSol]].first->materialLayer();
1332 if (mb) {
1333 if (
mb->layerMaterialProperties() &&
1334 mb->layerMaterialProperties()->fullMaterial(nextPar->position())) {
1335
1336 const IMaterialEffectsUpdator* currentUpdator =
1338 IMaterialEffectsUpdator::ICache& currentUpdatorCache =
1339 cache.subMaterialEffectsUpdatorCache();
1340
1341 if (currentUpdator) {
1342 nextPar = cache.m_ownedPtrs.push(
1343 currentUpdator->update(currentUpdatorCache, nextPar,
1344 *mb, dir, particle, matupmode));
1345 }
1346 if (!nextPar) {
1347 cache.m_parametersAtBoundary.resetBoundaryInformation();
1348 return {};
1349 }
1350
1351
1352 const Trk::MaterialProperties* lmat =
mb->fullUpdateMaterialProperties(*nextPar);
1353 const double lx0 = lmat->
x0();
1354 const double layThick =
mb->thickness();
1355
1356 double thick = 0.;
1357 const double costr =
1358 std::abs(nextPar->momentum().normalized().dot(
mb->surfaceRepresentation().normal()));
1359
1360 if (
mb->surfaceRepresentation().isOnSurface(
1361 mb->surfaceRepresentation().center(),
false, 0., 0.)) {
1362 thick = fmin(
mb->surfaceRepresentation().bounds().r(),
1363 layThick / std::abs(nextPar->momentum().normalized().dot(
1364 mb->surfaceRepresentation().normal())));
1365 } else {
1366 thick = fmin(2 *
mb->thickness(), layThick / (1 - costr));
1367 }
1368
1369 if (!cache.m_matstates && cache.m_extrapolationCache) {
1370 if (not cache.elossPointerOverwritten()) {
1371 const double dInX0 = thick / lx0;
1373 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer thin "));
1374 }
1375 cache.m_extrapolationCache->updateX0(dInX0);
1376 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
1378 *lmat, std::abs(1. / currentqoverp), 1. / costr, dir, particle);
1379 cache.m_extrapolationCache->updateEloss(
1380 eloss.meanIoni(), eloss.sigmaIoni(), eloss.meanRad(), eloss.sigmaRad());
1383 }
1384 } else {
1386 }
1387
1388 }
1389
1390 if (cache.m_matstates) {
1391 const double dInX0 = thick / lx0;
1392 const double scatsigma = std::sqrt(
m_msupdater->sigmaSquare(
1393 *lmat, 1. / std::abs(nextPar->parameters()[
qOverP]), 1., particle));
1394 auto newsa = Trk::ScatteringAngles(
1395 0, 0, scatsigma / std::sin(nextPar->parameters()[
Trk::theta]), scatsigma);
1396
1397 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
1399 *lmat, std::abs(1. / currentqoverp), 1. / costr, dir, particle);
1400
1401
1403 nextPar->position(), nextPar->momentum(), nextPar->charge()));
1404 if (cache.m_extrapolationCache) {
1405 if (not cache.elossPointerOverwritten()) {
1407 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer thin "));
1408 }
1409 cache.m_extrapolationCache->updateX0(dInX0);
1410 cache.m_extrapolationCache->updateEloss(
1411 eloss.meanIoni(), eloss.sigmaIoni(), eloss.meanRad(), eloss.sigmaRad());
1414 }
1415 } else {
1417 }
1418 }
1419 auto mefot =
1420 std::make_unique<Trk::MaterialEffectsOnTrack>(
1421 dInX0, newsa,
1422 std::make_unique<Trk::EnergyLoss>(std::move(eloss)),
1423 cvlTP->associatedSurface());
1424 cache.m_matstates->push_back(new TrackStateOnSurface(
1425 nullptr, std::move(cvlTP), std::move(mefot)));
1426 }
1427 }
1428 }
1429
1430
1431 const unsigned int index = solutions[iSol] - iDest;
1432
1433 nextVol = (cache.m_currentStatic->boundarySurfaces())[
index]->attachedVolume(
1434 nextPar->position(), nextPar->momentum(), dir);
1435
1436 if (nextVol &&
1437 !(nextVol->
inside(nextPar->position() + 0.01 * dir * nextPar->momentum().normalized(),
1439 ATH_MSG_DEBUG(
" [!] WARNING: wrongly assigned static volume ?"
1440 << cache.m_currentStatic->volumeName() <<
"->" << nextVol->
volumeName());
1441 nextVol = cache.m_trackingGeometry->lowestStaticTrackingVolume(
1442 nextPar->position() + 0.01 * nextPar->momentum().normalized());
1443 if (nextVol) {
1445 }
1446 }
1447
1448 if (nextVol != cache.m_currentStatic) {
1449 cache.m_parametersAtBoundary.boundaryInformation(nextVol, nextPar, nextPar);
1451 << cache.m_currentStatic->volumeName() << "'.");
1453 nextPar, cache.m_currentStatic, dir, assocVol,
m_tolerance) &&
1454 assocVol != cache.m_currentStatic) {
1456 }
1457
1458 if (!nextVol) {
1461 }
1462
1463 if (nextVol && nextPar) {
1467 }
1468 return {};
1469 }
1470 } else if (solutions[iSol] <
1471 iDest + cache.m_staticBoundaries.size() + cache.m_layers.size()) {
1472
1473 const unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size();
1474 const Trk::Layer* nextLayer = cache.m_navigLays[
index].second;
1475
1476
1478 if (nextLayer == cache.m_lastMaterialLayer &&
1481 " [!] This layer is identical to the one with last material update, return layer "
1482 "without repeating the update");
1484 if (!destSurf && (nextLayer->
layerType() > 0)) {
1485 return nextPar;
1486 }
1487 }
1488 const double layThick = nextLayer->
thickness();
1489 if (collect && layThick > 0.) {
1490
1491
1492 const IMaterialEffectsUpdator* currentUpdator =
1494 IMaterialEffectsUpdator::ICache& currentUpdatorCache =
1495 cache.subMaterialEffectsUpdatorCache();
1496
1497 if (currentUpdator) {
1498 nextPar = cache.m_ownedPtrs.push(
1499 currentUpdator->update(currentUpdatorCache, nextPar,
1500 *nextLayer, dir, particle, matupmode));
1501 }
1502 if (!nextPar) {
1503 cache.m_parametersAtBoundary.resetBoundaryInformation();
1504 return {};
1505 }
1506
1507
1509
1510 double thick = 0.;
1511 const double costr = std::abs(
1513
1517 layThick / std::abs(nextPar->momentum().normalized().dot(
1519 } else {
1520 thick = fmin(2 * nextLayer->
thickness(), layThick / (1 - costr));
1521 }
1522
1523 if (!cache.m_matstates && cache.m_extrapolationCache) {
1524 if (not cache.elossPointerOverwritten()) {
1525 const double dInX0 = thick / lx0;
1527 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer thin "));
1528 }
1529 cache.m_extrapolationCache->updateX0(dInX0);
1530 const Trk::MaterialProperties materialProperties(
1532 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
1534 materialProperties, std::abs(1. / currentqoverp), 1. / costr, dir, particle);
1535 cache.m_extrapolationCache->updateEloss(
1536 eloss.meanIoni(), eloss.sigmaIoni(), eloss.meanRad(), eloss.sigmaRad());
1539 }
1540 } else {
1542 }
1543 }
1544
1545 if (cache.m_matstates) {
1546 const double dInX0 = thick / lx0;
1547 const Trk::MaterialProperties materialProperties(
1549 const double scatsigma = std::sqrt(
m_msupdater->sigmaSquare(
1550 materialProperties, 1. / std::abs(nextPar->parameters()[
qOverP]), 1., particle));
1551 const double par_theta = std::abs(nextPar->parameters()[
Trk::theta]) > FLT_EPSILON
1553 : FLT_EPSILON;
1554 Trk::ScatteringAngles newsa(0, 0, scatsigma / std::sin(par_theta), scatsigma);
1555
1556 const double currentqoverp = nextPar->parameters()[
Trk::qOverP];
1558 materialProperties, std::abs(1. / currentqoverp), 1. / costr, dir, particle);
1559
1560
1561 auto cvlTP = std::make_unique<Trk::CurvilinearParameters>(
1562 nextPar->position(), nextPar->momentum(), nextPar->charge());
1563 if (cache.m_extrapolationCache) {
1564 if (not cache.elossPointerOverwritten()) {
1566 ATH_MSG_DEBUG(cache.to_string(
" extrapolateToNextMaterialLayer thin "));
1567 }
1568 cache.m_extrapolationCache->updateX0(dInX0);
1569 cache.m_extrapolationCache->updateEloss(
1570 eloss.meanIoni(), eloss.sigmaIoni(), eloss.meanRad(), eloss.sigmaRad());
1573 }
1574 } else {
1576 }
1577 }
1578 auto mefot = std::make_unique<Trk::MaterialEffectsOnTrack>(
1579 dInX0, newsa, std::make_unique<Trk::EnergyLoss>(std::move(eloss)), cvlTP->associatedSurface());
1580 cache.m_matstates->push_back(new TrackStateOnSurface(
1581 nullptr, std::move(cvlTP), std::move(mefot)));
1582 }
1583
1584 if (cache.m_cacheLastMatLayer) {
1585 cache.m_lastMaterialLayer = nextLayer;
1586 }
1587 if (!destSurf && nextLayer->
layerType() > 0) {
1588 return nextPar;
1589 }
1590 }
1591 if (resolveActive) {
1592
1593 if (cache.m_navigLays[index].first &&
1594 cache.m_navigLays[index].first->confinedLayers()) {
1595 const Trk::Layer* newLayer = cache.m_navigLays[
index].first->nextLayer(
1596 nextPar->position(), dir * nextPar->momentum().normalized(), true);
1597 if (newLayer) {
1598 cache.m_navigLays[
index].second = newLayer;
1600 }
1601 }
1602 }
1603
1604
1605 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() +
1606 cache.m_layers.size() + cache.m_denseBoundaries.size()) {
1607
1608 unsigned int index =
1609 solutions[iSol] - iDest - cache.m_staticBoundaries.size() - cache.m_layers.size();
1610 std::vector<std::pair<const Trk::TrackingVolume*, unsigned int>>::iterator dIter =
1611 cache.m_denseVols.begin();
1612 while (dIter != cache.m_denseVols.end() && index >= (*dIter).second) {
1613 index -= (*dIter).second;
1614 ++dIter;
1615 }
1616 if (dIter != cache.m_denseVols.end()) {
1617 currVol = (*dIter).first;
1618 nextVol =
1619 ((*dIter).first->boundarySurfaces())[
index]->attachedVolume(*nextPar, dir);
1620
1622 nextPar->position() + 2 *
m_tolerance *
dir * nextPar->momentum().normalized();
1624 cache.m_currentDense = currVol;
1626 cache.m_currentDense = cache.m_highestVolume;
1627 if (cache.m_dense && cache.m_denseVols.empty()) {
1628 cache.m_currentDense = cache.m_currentStatic;
1629 } else {
1630 for (
unsigned int i = 0;
i < cache.m_denseVols.size();
i++) {
1631 const Trk::TrackingVolume* dVol = cache.m_denseVols[
i].first;
1633 cache.m_currentDense = dVol;
1635 << cache.m_currentDense->volumeName() << "'.");
1636 break;
1637 }
1638 }
1639 }
1640 } else {
1641 cache.m_currentDense = nextVol;
1642 ATH_MSG_DEBUG(
" [+] Next dense volume: '" << cache.m_currentDense->volumeName()
1643 << "'.");
1644 }
1645 }
1646 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() +
1647 cache.m_layers.size() + cache.m_denseBoundaries.size() +
1648 cache.m_navigBoundaries.size()) {
1649
1650 unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size() -
1651 cache.m_layers.size() - cache.m_denseBoundaries.size();
1652 std::vector<std::pair<const Trk::TrackingVolume*, unsigned int>>::iterator nIter =
1653 cache.m_navigVolsInt.begin();
1654 while (nIter != cache.m_navigVolsInt.end() && index >= (*nIter).second) {
1655 index -= (*nIter).second;
1656 ++nIter;
1657 }
1658 if (nIter != cache.m_navigVolsInt.end()) {
1659 currVol = (*nIter).first;
1660 nextVol =
1661 ((*nIter).first->boundarySurfaces())[
index]->attachedVolume(*nextPar, dir);
1662
1664 nextPar->position() + 2 *
m_tolerance *
dir * nextPar->momentum().normalized();
1665 if (nextVol && nextVol->
inside(tp, 0.)) {
1666 ATH_MSG_DEBUG(
" [+] Navigation volume boundary, entering volume '"
1668 }
else if (currVol->
inside(tp, 0.)) {
1669 nextVol = currVol;
1670 ATH_MSG_DEBUG(
" [+] Navigation volume boundary, entering volume '"
1672 } else {
1673 nextVol = nullptr;
1674 ATH_MSG_DEBUG(
" [+] Navigation volume boundary, leaving volume '"
1676 }
1677
1678
1679 if (nextVol) {
1681 ctx, cache, prop, nextPar, destSurf, cache.m_currentStatic,
1682 dir, bcheck, particle, matupmode);
1683 }
1684 }
1685 } else if (solutions[iSol] < iDest + cache.m_staticBoundaries.size() +
1686 cache.m_layers.size() + cache.m_denseBoundaries.size() +
1687 cache.m_navigBoundaries.size() +
1688 cache.m_detachedBoundaries.size()) {
1689
1690 unsigned int index = solutions[iSol] - iDest - cache.m_staticBoundaries.size() -
1691 cache.m_layers.size() - cache.m_denseBoundaries.size() -
1692 cache.m_navigBoundaries.size();
1693 std::vector<std::pair<const Trk::DetachedTrackingVolume*, unsigned int>>::iterator dIter =
1694 cache.m_detachedVols.begin();
1695 while (dIter != cache.m_detachedVols.end() && index >= (*dIter).second) {
1696 index -= (*dIter).second;
1697 ++dIter;
1698 }
1699 if (dIter != cache.m_detachedVols.end()) {
1700 currVol = (*dIter).first->trackingVolume();
1701
1702 nextVol =
1703 ((*dIter).first->trackingVolume()->boundarySurfaces())[
index]->attachedVolume(
1704 *nextPar, dir);
1706 nextPar->position() + 2 *
m_tolerance *
dir * nextPar->momentum().normalized();
1707 if (nextVol && nextVol->
inside(tp, 0.)) {
1708 ATH_MSG_DEBUG(
" [+] Detached volume boundary, entering volume '"
1710 }
else if (currVol->
inside(tp, 0.)) {
1711 nextVol = currVol;
1712 ATH_MSG_DEBUG(
" [+] Detached volume boundary, entering volume '"
1714 } else {
1715 nextVol = nullptr;
1716 ATH_MSG_DEBUG(
" [+] Detached volume boundary, leaving volume '"
1718 }
1719
1720
1721 if (nextVol) {
1723 ctx, cache, prop, nextPar, destSurf, cache.m_currentStatic,
1724 dir, bcheck, particle, matupmode);
1725 }
1726 }
1727 }
1728 iSol++;
1729 }
1730 } else {
1732 cache.m_parametersAtBoundary.boundaryInformation(cache.m_currentStatic, nextPar, nextPar);
1733 return {};
1734 }
1735 currPar = nextPar;
1736 }
1737
1738 return {};
1739}
virtual std::span< T *const > arrayObjects()=0
Return all objects of the Array non-const we can still modify the T.
double thickness() const
Return the Thickness of the Layer.
float x0() const
Return the radiation length.
float zOverAtimesRho() const
access to members
virtual double r() const =0
Interface method for the maximal extension or the radius.
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
virtual constexpr SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
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...
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
const Amg::Vector3D & center() const
Returns the center position of the Surface.
const LayerArray * confinedLayers() const
Return the subLayer array.
const Layer * associatedLayer(const Amg::Vector3D &gp) const
Return the associated Layer.
const TrackingVolumeArray * confinedVolumes() const
Return the subLayer array.
const TrackingVolume * associatedSubVolume(const Amg::Vector3D &gp) const
Return the associated sub Volume, returns THIS if no subVolume exists.
std::vector< std::shared_ptr< BoundarySurface< TrackingVolume > > > & boundarySurfaces()
Method to return the BoundarySurfaces.
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.
ArraySpan< DetachedTrackingVolume const *const > confinedDetachedVolumes() const
Return detached subVolumes - not the ownership.
ArraySpan< Layer const *const > confinedArbitraryLayers() const
Return the confined subLayer array.
ArraySpan< TrackingVolume const *const > confinedDenseVolumes() const
Return unordered subVolumes - not the ownership.
layers(flags, cells_name, *args, **kw)
Here we define wrapper functions to set up all of the standard corrections.
const double mb
1mb to cm2
void collect(const HLT::TriggerElement *te, std::vector< Trig::Feature< T > > &data, const std::string &label, unsigned int condition, const std::string &teName, const HLT::TrigNavStructure *navstructure)
actual feature acceess implementation It has (thanks to the ClassTraits) functionality to flatten con...
CurvilinearParametersT< TrackParametersDim, Charged, PlaneSurface > CurvilinearParameters
@ kRecursionCountExceeded
1.14.0