1056 {
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082 std::vector<const Trk::TrackStateOnSurface*> newTSOSvector;
1083 int maxsize = 2 * matvec.size();
1084 if (aggregate) maxsize = 2;
1085 newTSOSvector.reserve(maxsize);
1086
1087
1088
1089
1090 Eloss_tot = 0.;
1091
1092 double X0_tot = 0.;
1093
1094 double sigmaDeltaPhi2_tot = 0.;
1095 double sigmaDeltaTheta2_tot = 0.;
1096 double deltaE_tot = 0.;
1097 double sigmaDeltaE_tot = 0.;
1098 double sigmaPlusDeltaE_tot = 0.;
1099 double sigmaMinusDeltaE_tot = 0.;
1100 double deltaE_ioni_tot = 0.;
1101 double sigmaDeltaE_ioni_tot = 0.;
1102 double deltaE_rad_tot = 0.;
1103 double sigmaDeltaE_rad_tot = 0.;
1104
1105 const Trk::TrackStateOnSurface* mprevious = nullptr;
1106 const Trk::TrackStateOnSurface* mfirst = nullptr;
1107 const Trk::TrackStateOnSurface* mlast = nullptr;
1109
1110 double deltaEFirst = 0.;
1111
1113 double deltaTheta = 0.;
1114
1116
1117 double w_tot = 0.;
1118 double wdist2 = 0.;
1121
1122 std::bitset<Trk::MaterialEffectsBase::NumberOfMaterialEffectsTypes>
1123 meotPattern(0);
1126
1127
1128 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1129 typePattern(0);
1132
1133 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1134 typePatternDeposit(0);
1138
1139 for (const auto* m : matvec) {
1140 if (!
m->trackParameters()) {
1142 continue;
1143 }
1144 if (
m->materialEffectsOnTrack()) {
1145 double X0 =
m->materialEffectsOnTrack()->thicknessInX0();
1146 const Trk::MaterialEffectsOnTrack* meot =
1147 dynamic_cast<const Trk::MaterialEffectsOnTrack*>(
1148 m->materialEffectsOnTrack());
1149 const Trk::EnergyLoss* energyLoss = nullptr;
1150 const Trk::ScatteringAngles* scat = nullptr;
1151 if (meot) {
1153 if (energyLoss) {
1154
1155 } else {
1157 continue;
1158 }
1160 if (scat) {
1161
1162 } else {
1164 continue;
1165 }
1166 } else {
1168 continue;
1169 }
1170
1174 <<
m->dumpType() <<
" TSOS surface "
1175 <<
m->trackParameters()->associatedSurface()
1176 <<
" position x " <<
m->trackParameters()->position().x()
1177 <<
" y " <<
m->trackParameters()->position().y() <<
" z "
1178 <<
m->trackParameters()->position().z() <<
" direction x "
1179 <<
m->trackParameters()->momentum().unit().x() <<
" y "
1180 <<
m->trackParameters()->momentum().unit().y() <<
" z "
1181 <<
m->trackParameters()->momentum().unit().z() <<
" p "
1182 <<
m->trackParameters()->momentum().mag() <<
" X0 " << X0
1183 <<
" deltaE " << energyLoss->
deltaE()
1185 <<
" depth " <<
depth);
1186
1187 X0_tot += scaleX0 *
X0;
1188
1189 sigmaDeltaTheta2_tot +=
1191 sigmaDeltaPhi2_tot +=
1193
1194
1195
1196
1197 deltaE_tot += scaleEloss * energyLoss->
deltaE();
1198 sigmaDeltaE_tot += scaleEloss * energyLoss->
sigmaDeltaE();
1201 deltaE_ioni_tot += scaleEloss * energyLoss->
meanIoni();
1202 sigmaDeltaE_ioni_tot += scaleEloss * energyLoss->
sigmaIoni();
1203 deltaE_rad_tot += scaleEloss * energyLoss->
meanRad();
1204 sigmaDeltaE_rad_tot += scaleEloss * energyLoss->
sigmaRad();
1205
1207
1210 if (mprevious) {
1212 }
1213
1216 <<
pos.z() <<
" perp " <<
pos.perp());
1221
1223 << pos0.x() << " y " << pos0.y() << " z " << pos0.z()
1224 << " perp " << pos0.perp());
1226 << posNew.x() << " y " << posNew.y() << " z " << posNew.z()
1227 << " perp " << posNew.perp() << " distance "
1228 << (pos0 - posNew).mag() <<
" depth " <<
depth);
1229 if (!mfirst) {
1231 posFirst = pos0;
1232 deltaEFirst = energyLoss->
deltaE();
1233 }
1235
1238 wpos +=
w * pos0 / 2.;
1239 wpos +=
w * posNew / 2.;
1241
1242 wdist2 +=
w * (pos0 - posFirst).
mag2() / 2.;
1243 wdist2 +=
w * (posNew - posFirst).
mag2() / 2.;
1244
1245 if (!aggregate && !reposition) {
1246 auto scatNew = ScatteringAngles(
deltaPhi, deltaTheta,
1247 std::sqrt(sigmaDeltaPhi2_tot),
1248 std::sqrt(sigmaDeltaTheta2_tot));
1249 auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1250 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1251 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1252 deltaE_rad_tot, sigmaDeltaE_rad_tot,
depth);
1253 Eloss_tot += energyLossNew->deltaE();
1255 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1256 X0_tot, scatNew, std::move(energyLossNew), surf, meotPattern);
1257 auto pars =
m->trackParameters()->uniqueClone();
1258
1259
1260 const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1261 nullptr, std::move(pars), std::move(meotLast), typePattern);
1262 newTSOSvector.push_back(newTSOS);
1263
1264 X0_tot = 0.;
1265 sigmaDeltaTheta2_tot = 0.;
1266 sigmaDeltaPhi2_tot = 0.;
1267 deltaE_tot = 0.;
1268 sigmaDeltaE_tot = 0;
1269 sigmaPlusDeltaE_tot = 0.;
1270 sigmaMinusDeltaE_tot = 0.;
1271 deltaE_ioni_tot = 0.;
1272 sigmaDeltaE_ioni_tot = 0.;
1273 deltaE_rad_tot = 0.;
1274 sigmaDeltaE_rad_tot = 0.;
1275
1276 } else if (!aggregate && reposition) {
1277 if (std::abs(
depth) < 10.) {
1278 auto scatNew = ScatteringAngles(
deltaPhi, deltaTheta,
1279 std::sqrt(sigmaDeltaPhi2_tot),
1280 std::sqrt(sigmaDeltaTheta2_tot));
1281 auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1282 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1283 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1284 deltaE_rad_tot, sigmaDeltaE_rad_tot,
depth);
1286 Eloss_tot += energyLossNew->deltaE();
1287 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1288 X0_tot, scatNew, std::move(energyLossNew), surf, meotPattern);
1289 std::unique_ptr<Trk::TrackParameters>
pars =
1290 m->trackParameters()->uniqueClone();
1291
1292 const Trk::TrackStateOnSurface* newTSOS =
1293 new Trk::TrackStateOnSurface(nullptr, std::move(pars),
1294 std::move(meotLast), typePattern);
1295 newTSOSvector.push_back(newTSOS);
1296 X0_tot = 0.;
1297 sigmaDeltaTheta2_tot = 0.;
1298 sigmaDeltaPhi2_tot = 0.;
1299 deltaE_tot = 0.;
1300 sigmaDeltaE_tot = 0;
1301 sigmaPlusDeltaE_tot = 0.;
1302 sigmaMinusDeltaE_tot = 0.;
1303 deltaE_ioni_tot = 0.;
1304 sigmaDeltaE_ioni_tot = 0.;
1305 deltaE_rad_tot = 0.;
1306 sigmaDeltaE_rad_tot = 0.;
1307
1308 } else {
1309
1310
1311
1312
1313
1314 auto energyLoss0 = std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.);
1315 auto scatFirst = ScatteringAngles(
deltaPhi, deltaTheta,
1316 sqrt(sigmaDeltaPhi2_tot / 2.),
1317 sqrt(sigmaDeltaTheta2_tot / 2.));
1318
1319
1320
1321 auto scatNew = ScatteringAngles(
deltaPhi, deltaTheta,
1322 sqrt(sigmaDeltaPhi2_tot / 2.),
1323 sqrt(sigmaDeltaTheta2_tot / 2.));
1324 auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1325 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1326 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1327 deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.);
1329
1332 -
dir.y() *
dir.z() / norm, norm);
1334
1337 Trk::PlaneSurface* surfFirst =
1338 new Trk::PlaneSurface(surfaceTransformFirst);
1339 Trk::PlaneSurface* surfLast =
1340 new Trk::PlaneSurface(surfaceTransformLast);
1341 Eloss_tot += energyLossNew->deltaE();
1342
1343 auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1344 X0_tot / 2., scatFirst, std::move(energyLoss0), *surfFirst,
1345 meotPattern);
1346 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1347 X0_tot / 2., scatNew, std::move(energyLossNew), *surfLast,
1348 meotPattern);
1349
1350
1351 double qOverP0 =
m->trackParameters()->charge() /
1353 std::fabs(energyLoss->
deltaE()));
1354 if (mprevious)
1357 std::unique_ptr<Trk::TrackParameters> parsFirst =
1359 0., 0.,
dir.phi(),
dir.theta(), qOverP0);
1360
1361 double qOverPNew =
m->trackParameters()->charge() /
1362 m->trackParameters()->momentum().mag();
1363 std::unique_ptr<Trk::TrackParameters> parsLast =
1365 0., 0.,
dir.phi(),
dir.theta(), qOverPNew);
1366
1367
1368 const Trk::TrackStateOnSurface* newTSOSFirst =
1369 new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1370 std::move(meotFirst), typePattern);
1371 const Trk::TrackStateOnSurface* newTSOS =
1372 new Trk::TrackStateOnSurface(nullptr, std::move(parsLast),
1373 std::move(meotLast), typePattern);
1374
1375 newTSOSvector.push_back(newTSOSFirst);
1376 newTSOSvector.push_back(newTSOS);
1377
1378 X0_tot = 0.;
1379 sigmaDeltaTheta2_tot = 0.;
1380 sigmaDeltaPhi2_tot = 0.;
1381 deltaE_tot = 0.;
1382 sigmaDeltaE_tot = 0;
1383 sigmaPlusDeltaE_tot = 0.;
1384 sigmaMinusDeltaE_tot = 0.;
1385 deltaE_ioni_tot = 0.;
1386 sigmaDeltaE_ioni_tot = 0.;
1387 deltaE_rad_tot = 0.;
1388 sigmaDeltaE_rad_tot = 0.;
1389 }
1390 }
1391
1393 }
1394 }
1395 if (aggregate && reposition) {
1396 if (n_tot > 0) {
1397
1398
1399
1401 bool threePlanes = false;
1402 if (X0_tot > 50 && std::fabs(
pos.z()) < 6700 &&
pos.perp() < 4200)
1403 threePlanes = true;
1404
1405 auto energyLoss0 = std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.);
1406 auto scatFirst =
1407 ScatteringAngles(
deltaPhi, deltaTheta, sqrt(sigmaDeltaPhi2_tot / 2.),
1408 sqrt(sigmaDeltaTheta2_tot / 2.));
1409
1410 auto scatNew =
1411 ScatteringAngles(
deltaPhi, deltaTheta, sqrt(sigmaDeltaPhi2_tot / 2.),
1412 sqrt(sigmaDeltaTheta2_tot / 2.));
1413 auto energyLoss2 = Trk::EnergyLoss(
1414 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1415 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1416 deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.);
1417
1418 int elossFlag = 0;
1419
1420 auto energyLossNew =
1421 (updateEloss
1423 caloEnergyError, pCaloEntry,
1424 momentumError, elossFlag)
1425 : Trk::EnergyLoss(deltaE_tot, sigmaDeltaE_tot,
1426 sigmaPlusDeltaE_tot, sigmaMinusDeltaE_tot,
1427 deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1428 deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.));
1429
1430
1434
1437 norm);
1439
1440 double halflength2 =
1441 wdist2 / w_tot - (
pos - posFirst).
mag() * (
pos - posFirst).
mag();
1442 double halflength = 0.;
1443 if (halflength2 > 0) halflength = sqrt(halflength2);
1447
1448 ATH_MSG_DEBUG(
" WITH aggregation and WITH reposition center planes x "
1449 <<
pos.x() <<
" y " <<
pos.y() <<
" z " <<
pos.z()
1450 << " halflength " << halflength << " w_tot " << w_tot
1451 << " X0_tot " << X0_tot);
1452
1455 Trk::PlaneSurface* surfFirst =
1456 new Trk::PlaneSurface(surfaceTransformFirst);
1457 Trk::PlaneSurface* surfLast = new Trk::PlaneSurface(surfaceTransformLast);
1458
1461 std::fabs(deltaEFirst));
1462
1465 std::unique_ptr<Trk::TrackParameters> parsFirst =
1467 0., 0.,
dir.phi(),
dir.theta(), qOverP0);
1468 std::unique_ptr<Trk::TrackParameters> parsLast =
1470 0., 0.,
dir.phi(),
dir.theta(), qOverPNew);
1471
1472 Eloss_tot += energyLossNew.deltaE();
1473 if (!threePlanes) {
1474
1475
1476
1477
1478
1479 auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1480 X0_tot / 2., scatFirst, std::move(energyLoss0), *surfFirst,
1481 meotPattern);
1482
1483
1484
1485 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1486 X0_tot / 2., scatNew,
1487 std::make_unique<Trk::EnergyLoss>(std::move(energyLossNew)),
1488 *surfLast, meotPattern);
1489
1490 const Trk::TrackStateOnSurface* newTSOSFirst =
1491 new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1492 std::move(meotFirst), typePattern);
1493 auto whichType = (elossFlag != 0) ? typePatternDeposit : typePattern;
1494 const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1495 nullptr, std::move(parsLast),
1496 std::move(meotLast), whichType);
1497
1498 newTSOSvector.push_back(newTSOSFirst);
1499 newTSOSvector.push_back(newTSOS);
1500 } else {
1501
1502
1503
1504 auto scatZero = ScatteringAngles(0., 0., 0., 0.);
1506 Trk::PlaneSurface* surf = new Trk::PlaneSurface(surfaceTransform);
1507 std::unique_ptr<Trk::TrackParameters>
pars =
1509 0., 0.,
dir.phi(),
dir.theta(), qOverPNew);
1510
1511
1512 auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1513 X0_tot / 2., scatFirst,
1514 std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.), *surfFirst,
1515 meotPattern);
1516
1517
1518 auto meot = std::make_unique<Trk::MaterialEffectsOnTrack>(
1519 0., scatZero,
1520 std::make_unique<Trk::EnergyLoss>(std::move(energyLossNew)), *surf,
1521 meotPattern);
1522
1523
1524 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1525 X0_tot / 2., scatNew,
1526 std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.), *surfLast,
1527 meotPattern);
1528 const Trk::TrackStateOnSurface* newTSOSFirst =
1529 new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1530 std::move(meotFirst), typePattern);
1531 const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1532 nullptr, std::move(pars), std::move(meot), typePatternDeposit);
1533 const Trk::TrackStateOnSurface* newTSOSLast =
1534 new Trk::TrackStateOnSurface(nullptr, std::move(parsLast),
1535 std::move(meotLast), typePattern);
1536 newTSOSvector.push_back(newTSOSFirst);
1537 newTSOSvector.push_back(newTSOS);
1538 newTSOSvector.push_back(newTSOSLast);
1539 }
1540 }
1541 }
1542
1543 return newTSOSvector;
1544}
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar mag() const
mag method
Scalar mag2() const
mag2 method - forward to squaredNorm()
#define ATH_MSG_WARNING(x)
double sigmaPlusDeltaE() const
returns the positive side
double sigmaMinusDeltaE() const
returns the negative side
double sigmaDeltaE() const
returns the symmatric error
double deltaE() const
returns the
@ ScatteringEffects
contains material effects due to multiple scattering
@ EnergyLossEffects
contains energy loss corrections
const Surface & associatedSurface() const
returns the surface to which these m.eff. are associated.
const EnergyLoss * energyLoss() const
returns the energy loss object.
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
const Amg::Vector3D & momentum() const
Access method for the momentum.
double charge() const
Returns the charge.
std::unique_ptr< ParametersT< DIM, T, PlaneSurface > > createUniqueParameters(double l1, double l2, double phi, double theta, double qop, std::optional< AmgSymMatrix(DIM)> cov=std::nullopt) const
Use the Surface as a ParametersBase constructor, from local parameters.
double sigmaDeltaPhi() const
returns the
double sigmaDeltaTheta() const
returns the
const TrackParameters * trackParameters() const
return ptr to trackparameters const overload
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
@ CaloDeposit
This TSOS contains a CaloEnergy object.
std::string depth
tag string for intendation
Eigen::Affine3d Transform3D
Eigen::Matrix< double, 3, 1 > Vector3D