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
1083
1084 std::vector<const Trk::TrackStateOnSurface*> newTSOSvector;
1085 int maxsize = 2 * matvec.size();
1086 if (aggregate) maxsize = 2;
1087 newTSOSvector.reserve(maxsize);
1088
1089
1090
1091
1092 Eloss_tot = 0.;
1093
1094 double X0_tot = 0.;
1095
1096 double sigmaDeltaPhi2_tot = 0.;
1097 double sigmaDeltaTheta2_tot = 0.;
1098 double deltaE_tot = 0.;
1099 double sigmaDeltaE_tot = 0.;
1100 double sigmaPlusDeltaE_tot = 0.;
1101 double sigmaMinusDeltaE_tot = 0.;
1102 double deltaE_ioni_tot = 0.;
1103 double sigmaDeltaE_ioni_tot = 0.;
1104 double deltaE_rad_tot = 0.;
1105 double sigmaDeltaE_rad_tot = 0.;
1106
1107 const Trk::TrackStateOnSurface* mprevious = nullptr;
1108 const Trk::TrackStateOnSurface* mfirst = nullptr;
1109 const Trk::TrackStateOnSurface* mlast = nullptr;
1111
1112 double deltaEFirst = 0.;
1113
1115 double deltaTheta = 0.;
1116
1118
1119 double w_tot = 0.;
1120 double wdist2 = 0.;
1123
1124 std::bitset<Trk::MaterialEffectsBase::NumberOfMaterialEffectsTypes>
1125 meotPattern(0);
1128
1129
1130 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1131 typePattern(0);
1134
1135 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1136 typePatternDeposit(0);
1140
1141 for (const auto* m : matvec) {
1142 if (!
m->trackParameters()) {
1144 continue;
1145 }
1146 if (
m->materialEffectsOnTrack()) {
1147 double X0 =
m->materialEffectsOnTrack()->thicknessInX0();
1148 const Trk::MaterialEffectsOnTrack* meot =
1149 dynamic_cast<const Trk::MaterialEffectsOnTrack*>(
1150 m->materialEffectsOnTrack());
1151 const Trk::EnergyLoss* energyLoss = nullptr;
1152 const Trk::ScatteringAngles* scat = nullptr;
1153 if (meot) {
1155 if (energyLoss) {
1156
1157 } else {
1159 continue;
1160 }
1162 if (scat) {
1163
1164 } else {
1166 continue;
1167 }
1168 } else {
1170 continue;
1171 }
1172
1176 <<
m->dumpType() <<
" TSOS surface "
1177 <<
m->trackParameters()->associatedSurface()
1178 <<
" position x " <<
m->trackParameters()->position().x()
1179 <<
" y " <<
m->trackParameters()->position().y() <<
" z "
1180 <<
m->trackParameters()->position().z() <<
" direction x "
1181 <<
m->trackParameters()->momentum().unit().x() <<
" y "
1182 <<
m->trackParameters()->momentum().unit().y() <<
" z "
1183 <<
m->trackParameters()->momentum().unit().z() <<
" p "
1184 <<
m->trackParameters()->momentum().mag() <<
" X0 " << X0
1185 <<
" deltaE " << energyLoss->
deltaE()
1187 <<
" depth " <<
depth);
1188
1189 X0_tot += scaleX0 *
X0;
1190
1191 sigmaDeltaTheta2_tot +=
1193 sigmaDeltaPhi2_tot +=
1195
1196
1197
1198
1199 deltaE_tot += scaleEloss * energyLoss->
deltaE();
1200 sigmaDeltaE_tot += scaleEloss * energyLoss->
sigmaDeltaE();
1203 deltaE_ioni_tot += scaleEloss * energyLoss->
meanIoni();
1204 sigmaDeltaE_ioni_tot += scaleEloss * energyLoss->
sigmaIoni();
1205 deltaE_rad_tot += scaleEloss * energyLoss->
meanRad();
1206 sigmaDeltaE_rad_tot += scaleEloss * energyLoss->
sigmaRad();
1207
1209
1212 if (mprevious) {
1214 }
1215
1218 <<
pos.z() <<
" perp " <<
pos.perp());
1223
1225 << pos0.x() << " y " << pos0.y() << " z " << pos0.z()
1226 << " perp " << pos0.perp());
1228 << posNew.x() << " y " << posNew.y() << " z " << posNew.z()
1229 << " perp " << posNew.perp() << " distance "
1230 << (pos0 - posNew).mag() <<
" depth " <<
depth);
1231 if (!mfirst) {
1233 posFirst = pos0;
1234 deltaEFirst = energyLoss->
deltaE();
1235 }
1237
1240 wpos +=
w * pos0 / 2.;
1241 wpos +=
w * posNew / 2.;
1243
1244 wdist2 +=
w * (pos0 - posFirst).
mag2() / 2.;
1245 wdist2 +=
w * (posNew - posFirst).
mag2() / 2.;
1246
1247 if (!aggregate && !reposition) {
1248 auto scatNew = ScatteringAngles(
deltaPhi, deltaTheta,
1249 std::sqrt(sigmaDeltaPhi2_tot),
1250 std::sqrt(sigmaDeltaTheta2_tot));
1251 auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1252 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1253 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1254 deltaE_rad_tot, sigmaDeltaE_rad_tot,
depth);
1255 Eloss_tot += energyLossNew->deltaE();
1257 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1258 X0_tot, scatNew, std::move(energyLossNew), surf, meotPattern);
1259 auto pars =
m->trackParameters()->uniqueClone();
1260
1261
1262 const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1263 nullptr, std::move(pars), std::move(meotLast), typePattern);
1264 newTSOSvector.push_back(newTSOS);
1265
1266 X0_tot = 0.;
1267 sigmaDeltaTheta2_tot = 0.;
1268 sigmaDeltaPhi2_tot = 0.;
1269 deltaE_tot = 0.;
1270 sigmaDeltaE_tot = 0;
1271 sigmaPlusDeltaE_tot = 0.;
1272 sigmaMinusDeltaE_tot = 0.;
1273 deltaE_ioni_tot = 0.;
1274 sigmaDeltaE_ioni_tot = 0.;
1275 deltaE_rad_tot = 0.;
1276 sigmaDeltaE_rad_tot = 0.;
1277
1278 } else if (!aggregate && reposition) {
1279 if (std::abs(
depth) < 10.) {
1280 auto scatNew = ScatteringAngles(
deltaPhi, deltaTheta,
1281 std::sqrt(sigmaDeltaPhi2_tot),
1282 std::sqrt(sigmaDeltaTheta2_tot));
1283 auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1284 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1285 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1286 deltaE_rad_tot, sigmaDeltaE_rad_tot,
depth);
1288 Eloss_tot += energyLossNew->deltaE();
1289 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1290 X0_tot, scatNew, std::move(energyLossNew), surf, meotPattern);
1291 std::unique_ptr<Trk::TrackParameters> pars =
1292 m->trackParameters()->uniqueClone();
1293
1294 const Trk::TrackStateOnSurface* newTSOS =
1295 new Trk::TrackStateOnSurface(nullptr, std::move(pars),
1296 std::move(meotLast), typePattern);
1297 newTSOSvector.push_back(newTSOS);
1298 X0_tot = 0.;
1299 sigmaDeltaTheta2_tot = 0.;
1300 sigmaDeltaPhi2_tot = 0.;
1301 deltaE_tot = 0.;
1302 sigmaDeltaE_tot = 0;
1303 sigmaPlusDeltaE_tot = 0.;
1304 sigmaMinusDeltaE_tot = 0.;
1305 deltaE_ioni_tot = 0.;
1306 sigmaDeltaE_ioni_tot = 0.;
1307 deltaE_rad_tot = 0.;
1308 sigmaDeltaE_rad_tot = 0.;
1309
1310 } else {
1311
1312
1313
1314
1315
1316 auto energyLoss0 = std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.);
1317 auto scatFirst = ScatteringAngles(
deltaPhi, deltaTheta,
1318 sqrt(sigmaDeltaPhi2_tot / 2.),
1319 sqrt(sigmaDeltaTheta2_tot / 2.));
1320
1321
1322
1323 auto scatNew = ScatteringAngles(
deltaPhi, deltaTheta,
1324 sqrt(sigmaDeltaPhi2_tot / 2.),
1325 sqrt(sigmaDeltaTheta2_tot / 2.));
1326 auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1327 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1328 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1329 deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.);
1331
1334 -
dir.y() *
dir.z() / norm, norm);
1336
1339 Trk::PlaneSurface* surfFirst =
1340 new Trk::PlaneSurface(surfaceTransformFirst);
1341 Trk::PlaneSurface* surfLast =
1342 new Trk::PlaneSurface(surfaceTransformLast);
1343 Eloss_tot += energyLossNew->deltaE();
1344
1345 auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1346 X0_tot / 2., scatFirst, std::move(energyLoss0), *surfFirst,
1347 meotPattern);
1348 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1349 X0_tot / 2., scatNew, std::move(energyLossNew), *surfLast,
1350 meotPattern);
1351
1352
1353 double qOverP0 =
m->trackParameters()->charge() /
1355 std::fabs(energyLoss->
deltaE()));
1356 if (mprevious)
1359 std::unique_ptr<Trk::TrackParameters> parsFirst =
1361 0., 0.,
dir.phi(),
dir.theta(), qOverP0);
1362
1363 double qOverPNew =
m->trackParameters()->charge() /
1364 m->trackParameters()->momentum().mag();
1365 std::unique_ptr<Trk::TrackParameters> parsLast =
1367 0., 0.,
dir.phi(),
dir.theta(), qOverPNew);
1368
1369
1370 const Trk::TrackStateOnSurface* newTSOSFirst =
1371 new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1372 std::move(meotFirst), typePattern);
1373 const Trk::TrackStateOnSurface* newTSOS =
1374 new Trk::TrackStateOnSurface(nullptr, std::move(parsLast),
1375 std::move(meotLast), typePattern);
1376
1377 newTSOSvector.push_back(newTSOSFirst);
1378 newTSOSvector.push_back(newTSOS);
1379
1380 X0_tot = 0.;
1381 sigmaDeltaTheta2_tot = 0.;
1382 sigmaDeltaPhi2_tot = 0.;
1383 deltaE_tot = 0.;
1384 sigmaDeltaE_tot = 0;
1385 sigmaPlusDeltaE_tot = 0.;
1386 sigmaMinusDeltaE_tot = 0.;
1387 deltaE_ioni_tot = 0.;
1388 sigmaDeltaE_ioni_tot = 0.;
1389 deltaE_rad_tot = 0.;
1390 sigmaDeltaE_rad_tot = 0.;
1391 }
1392 }
1393
1395 }
1396 }
1397 if (aggregate && reposition) {
1398 if (n_tot > 0) {
1399
1400
1401
1403 bool threePlanes = false;
1404 if (X0_tot > 50 && std::fabs(
pos.z()) < 6700 &&
pos.perp() < 4200)
1405 threePlanes = true;
1406
1407 auto energyLoss0 = std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.);
1408 auto scatFirst =
1409 ScatteringAngles(
deltaPhi, deltaTheta, sqrt(sigmaDeltaPhi2_tot / 2.),
1410 sqrt(sigmaDeltaTheta2_tot / 2.));
1411
1412 auto scatNew =
1413 ScatteringAngles(
deltaPhi, deltaTheta, sqrt(sigmaDeltaPhi2_tot / 2.),
1414 sqrt(sigmaDeltaTheta2_tot / 2.));
1415 auto energyLoss2 = Trk::EnergyLoss(
1416 deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1417 sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1418 deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.);
1419
1420 int elossFlag = 0;
1421
1422 auto energyLossNew =
1423 (updateEloss
1425 caloEnergyError, pCaloEntry,
1426 momentumError, elossFlag)
1427 : Trk::EnergyLoss(deltaE_tot, sigmaDeltaE_tot,
1428 sigmaPlusDeltaE_tot, sigmaMinusDeltaE_tot,
1429 deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1430 deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.));
1431
1432
1436
1439 norm);
1441
1442 double halflength2 =
1443 wdist2 / w_tot - (
pos - posFirst).
mag() * (
pos - posFirst).
mag();
1444 double halflength = 0.;
1445 if (halflength2 > 0) halflength = sqrt(halflength2);
1449
1450 ATH_MSG_DEBUG(
" WITH aggregation and WITH reposition center planes x "
1451 <<
pos.x() <<
" y " <<
pos.y() <<
" z " <<
pos.z()
1452 << " halflength " << halflength << " w_tot " << w_tot
1453 << " X0_tot " << X0_tot);
1454
1457 Trk::PlaneSurface* surfFirst =
1458 new Trk::PlaneSurface(surfaceTransformFirst);
1459 Trk::PlaneSurface* surfLast = new Trk::PlaneSurface(surfaceTransformLast);
1460
1463 std::fabs(deltaEFirst));
1464
1467 std::unique_ptr<Trk::TrackParameters> parsFirst =
1469 0., 0.,
dir.phi(),
dir.theta(), qOverP0);
1470 std::unique_ptr<Trk::TrackParameters> parsLast =
1472 0., 0.,
dir.phi(),
dir.theta(), qOverPNew);
1473
1474 Eloss_tot += energyLossNew.deltaE();
1475 if (!threePlanes) {
1476
1477
1478
1479
1480
1481 auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1482 X0_tot / 2., scatFirst, std::move(energyLoss0), *surfFirst,
1483 meotPattern);
1484
1485
1486
1487 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1488 X0_tot / 2., scatNew,
1489 std::make_unique<Trk::EnergyLoss>(std::move(energyLossNew)),
1490 *surfLast, meotPattern);
1491
1492 const Trk::TrackStateOnSurface* newTSOSFirst =
1493 new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1494 std::move(meotFirst), typePattern);
1495 auto whichType = (elossFlag != 0) ? typePatternDeposit : typePattern;
1496 const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1497 nullptr, std::move(parsLast),
1498 std::move(meotLast), whichType);
1499
1500 newTSOSvector.push_back(newTSOSFirst);
1501 newTSOSvector.push_back(newTSOS);
1502 } else {
1503
1504
1505
1506 auto scatZero = ScatteringAngles(0., 0., 0., 0.);
1508 Trk::PlaneSurface* surf = new Trk::PlaneSurface(surfaceTransform);
1509 std::unique_ptr<Trk::TrackParameters> pars =
1511 0., 0.,
dir.phi(),
dir.theta(), qOverPNew);
1512
1513
1514 auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1515 X0_tot / 2., scatFirst,
1516 std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.), *surfFirst,
1517 meotPattern);
1518
1519
1520 auto meot = std::make_unique<Trk::MaterialEffectsOnTrack>(
1521 0., scatZero,
1522 std::make_unique<Trk::EnergyLoss>(std::move(energyLossNew)), *surf,
1523 meotPattern);
1524
1525
1526 auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1527 X0_tot / 2., scatNew,
1528 std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.), *surfLast,
1529 meotPattern);
1530 const Trk::TrackStateOnSurface* newTSOSFirst =
1531 new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1532 std::move(meotFirst), typePattern);
1533 const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1534 nullptr, std::move(pars), std::move(meot), typePatternDeposit);
1535 const Trk::TrackStateOnSurface* newTSOSLast =
1536 new Trk::TrackStateOnSurface(nullptr, std::move(parsLast),
1537 std::move(meotLast), typePattern);
1538 newTSOSvector.push_back(newTSOSFirst);
1539 newTSOSvector.push_back(newTSOS);
1540 newTSOSvector.push_back(newTSOSLast);
1541 }
1542 }
1543 }
1544
1545 return newTSOSvector;
1546}
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar mag() const
mag method
Scalar mag2() const
mag2 method - forward to squaredNorm()
#define ATH_MSG_DEBUG(x,...)
#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, const 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