18 #include "GaudiKernel/MsgStream.h"
19 #include "GaudiKernel/SystemOfUnits.h"
42 const MeasurementBase* measurementBase)
44 m_alignmentEffects(nullptr),
45 m_alignmentParameter(0),
46 m_alignmentParameter2(0),
48 m_derivative(nullptr),
49 m_derivative2(nullptr),
54 m_flippedDriftDistance(false),
56 m_hitOnTrack(hitOnTrack),
58 m_materialEffects(nullptr),
59 m_materialEffectsOwner(false),
60 m_measurementBase(measurementBase),
61 m_minEnergyDeposit(0.),
62 m_minimizationDirection{},
64 m_numberDoF(measurementBase->localCovariance().cols()),
65 m_numericalDerivative(
false),
67 m_particleMassSquared(0.),
70 m_position(measurementBase->associatedSurface().center()),
72 m_radiationThickness(0.),
77 m_scatteringAngle(0.),
78 m_scatteringAngleOffSet(0.),
84 m_signedDriftDistance(0.),
86 m_surface(&measurementBase->associatedSurface()),
99 if (
dynamic_cast<const PlaneSurface*
>(m_surface)) {
100 m_normal =
Amg::Vector3D(m_surface->transform().rotation().col(2));
102 m_surface->localToGlobal(measurementBase->localParameters());
109 if (m_numberDoF == 2 &&
110 dynamic_cast<const TrapezoidBounds*
>(&m_surface->bounds())) {
113 double aa = measurementBase->localCovariance()(
locX,
locX);
114 double ab = measurementBase->localCovariance()(
locX,
locY);
115 double bb = measurementBase->localCovariance()(
locY,
locY);
116 double sum =
aa + bb;
117 double diff = std::sqrt(
sum *
sum - 4. * (
aa * bb - ab * ab));
120 double widthSq = 0.5 * (
sum -
diff);
121 sigma = std::sqrt(widthSq);
122 double term = 0.5 * (
aa - bb) /
diff;
123 double cosStereo = std::sqrt(0.5 - term);
124 double sinStereo = 0.0;
126 sinStereo = std::sqrt(0.5 + term);
127 if (ab * m_normal.z() < 0.)
128 sinStereo = -sinStereo;
140 m_minimizationDirection =
Amg::Vector3D(m_sensorDirection.cross(m_normal));
141 if (m_numberDoF == 2)
143 }
else if (
dynamic_cast<const StraightLineSurface*
>(m_surface)) {
145 if (!measurementBase->localParameters().contains(
locY)) {
154 m_signedDriftDistance = measurementBase->localParameters()[
driftRadius];
158 m_minimizationDirection =
160 double mag = m_surface->center().mag();
163 m_position = measurementBase->globalPosition();
166 m_signedDriftDistance = 0.;
171 const PerigeeSurface* perigee =
172 dynamic_cast<const PerigeeSurface*
>(m_surface);
174 m_position = m_surface->center();
177 if (m_numberDoF == 2)
194 m_weight = 1. /
sigma;
196 m_weight2 = 1. / sigma2;
204 : m_afterCalo(false),
205 m_alignmentEffects(nullptr),
206 m_alignmentParameter(0),
207 m_alignmentParameter2(0),
209 m_derivative(nullptr),
210 m_derivative2(nullptr),
215 m_flippedDriftDistance(false),
217 m_hitOnTrack(nullptr),
219 m_materialEffects(materialEffects),
220 m_materialEffectsOwner(false),
221 m_measurementBase(nullptr),
222 m_minEnergyDeposit(0.),
223 m_minimizationDirection{},
226 m_numericalDerivative(
false),
228 m_particleMassSquared(particleMass * particleMass),
231 m_position(position),
233 m_radiationThickness(materialEffects->thicknessInX0()),
238 m_scatteringAngle(0.),
239 m_scatteringAngleOffSet(0.),
240 m_secondResidual(0.),
245 m_signedDriftDistance(0.),
247 m_surface(&materialEffects->associatedSurface()),
251 if (
dynamic_cast<const CylinderSurface*
>(m_surface) ||
252 std::abs(m_surface->normal()(2)) < 0.5)
259 const EnergyLoss* energyLoss =
nullptr;
260 const ScatteringAngles* scattering =
nullptr;
261 const MaterialEffectsOnTrack* meot =
262 dynamic_cast<const MaterialEffectsOnTrack*
>(materialEffects);
264 energyLoss = meot->energyLoss();
265 scattering = meot->scatteringAngles();
273 m_energyLoss = std::abs(energyLoss->deltaE());
277 if (calo && !scattering && energyLoss->sigmaDeltaE() > 0.) {
278 m_energyLoss = energyLoss->deltaE();
280 m_sigma = energyLoss->sigmaDeltaE();
281 m_sigmaMinus = energyLoss->sigmaMinusDeltaE();
282 m_sigmaPlus = energyLoss->sigmaPlusDeltaE();
284 m_weight = 1. / m_sigma;
285 m_minEnergyDeposit = 0.5 * m_energyLoss;
286 if (m_energyLoss > 0. && m_energyLoss > 5.0 * m_sigmaMinus) {
287 m_minEnergyDeposit = m_energyLoss - 2.5 * m_sigmaMinus;
288 }
else if (m_energyLoss < 0. && m_energyLoss < -5.0 * m_sigmaMinus) {
289 m_minEnergyDeposit = m_energyLoss + 2.5 * m_sigmaMinus;
296 m_scatterPhi = scattering->deltaPhi();
297 m_scatterTheta = scattering->deltaTheta();
299 m_scatteringAngleOffSet = scattering->sigmaDeltaTheta();
312 : m_afterCalo(false),
313 m_alignmentEffects(nullptr),
314 m_alignmentParameter(0),
315 m_alignmentParameter2(0),
317 m_derivative(nullptr),
318 m_derivative2(nullptr),
321 m_energyLoss(-deltaE),
323 m_flippedDriftDistance(false),
325 m_hitOnTrack(nullptr),
327 m_materialEffects(nullptr),
328 m_materialEffectsOwner(true),
329 m_measurementBase(nullptr),
330 m_minEnergyDeposit(0.),
331 m_minimizationDirection{},
334 m_numericalDerivative(
false),
336 m_particleMassSquared(particleMass * particleMass),
339 m_position(position),
341 m_radiationThickness(radiationThickness),
346 m_scatteringAngle(0.),
347 m_scatteringAngleOffSet(0.),
348 m_secondResidual(0.),
353 m_signedDriftDistance(0.),
360 if (m_surface && (
dynamic_cast<const CylinderSurface*
>(m_surface) ||
361 std::abs(m_surface->normal()(2)) < 0.5))
363 else if (std::abs(direction(2)) < 0.5)
366 CurvilinearUVT uvt(direction);
367 m_surface =
new PlaneSurface(position, uvt);
371 std::bitset<MaterialEffectsBase::NumberOfMaterialEffectsTypes> typeMaterial;
374 auto energyLoss = std::make_unique<EnergyLoss>(deltaE, 0., 0., 0.);
375 m_materialEffects =
new MaterialEffectsOnTrack(
376 radiationThickness, std::move(energyLoss), *m_surface, typeMaterial);
379 m_surface = &m_materialEffects->associatedSurface();
382 std::make_optional<TrackSurfaceIntersection>(position, direction, 0.);
389 : m_afterCalo(false),
390 m_alignmentEffects(alignmentEffects),
391 m_alignmentParameter(0),
392 m_alignmentParameter2(0),
394 m_derivative(nullptr),
395 m_derivative2(nullptr),
400 m_flippedDriftDistance(false),
402 m_hitOnTrack(nullptr),
404 m_materialEffects(nullptr),
405 m_materialEffectsOwner(false),
406 m_measurementBase(nullptr),
407 m_minEnergyDeposit(0.),
408 m_minimizationDirection{},
411 m_numericalDerivative(
false),
413 m_particleMassSquared(0.),
416 m_position(position),
418 m_radiationThickness(0.),
421 m_scatterPhi(alignmentEffects->deltaAngle()),
422 m_scatterTheta(alignmentEffects->deltaTranslation()),
423 m_scatteringAngle(0.),
424 m_scatteringAngleOffSet(0.),
425 m_secondResidual(0.),
428 m_sigmaMinus(alignmentEffects->sigmaDeltaAngle()),
429 m_sigmaPlus(alignmentEffects->sigmaDeltaTranslation()),
430 m_signedDriftDistance(0.),
432 m_surface(&alignmentEffects->associatedSurface()),
438 m_weight = 1. / m_sigmaMinus;
440 m_weight2 = 1. / m_sigmaPlus;
443 std::make_optional<TrackSurfaceIntersection>(position, direction, 0.);
449 : m_afterCalo(false),
450 m_alignmentEffects(nullptr),
451 m_alignmentParameter(0),
452 m_alignmentParameter2(0),
454 m_derivative(nullptr),
455 m_derivative2(nullptr),
460 m_flippedDriftDistance(false),
462 m_hitOnTrack(nullptr),
464 m_materialEffects(nullptr),
465 m_materialEffectsOwner(false),
466 m_measurementBase(nullptr),
467 m_minEnergyDeposit(0.),
468 m_minimizationDirection{},
471 m_numericalDerivative(
false),
473 m_particleMassSquared(0.),
478 m_radiationThickness(0.),
483 m_scatteringAngle(0.),
484 m_scatteringAngleOffSet(0.),
485 m_secondResidual(0.),
490 m_signedDriftDistance(0.),
500 m_surface =
new PlaneSurface(m_position, uvt);
502 m_intersection[
FittedTrajectory] = std::make_optional<TrackSurfaceIntersection>(
508 : m_afterCalo(false),
509 m_alignmentEffects(nullptr),
510 m_alignmentParameter(0),
511 m_alignmentParameter2(0),
513 m_derivative(nullptr),
514 m_derivative2(nullptr),
519 m_flippedDriftDistance(false),
521 m_hitOnTrack(nullptr),
523 m_materialEffects(nullptr),
524 m_materialEffectsOwner(false),
525 m_measurementBase(nullptr),
526 m_minEnergyDeposit(0.),
527 m_minimizationDirection{},
530 m_numericalDerivative(
false),
532 m_particleMassSquared(0.),
535 m_position(TSOS.trackParameters()->position()),
537 m_radiationThickness(0.),
542 m_scatteringAngle(0.),
543 m_scatteringAngleOffSet(0.),
544 m_secondResidual(0.),
549 m_signedDriftDistance(0.),
551 m_surface(&TSOS.trackParameters()->associatedSurface()),
568 double sigma2,
double sinStereo,
int status,
570 : m_afterCalo(false),
571 m_alignmentEffects(nullptr),
572 m_alignmentParameter(0),
573 m_alignmentParameter2(0),
575 m_derivative(nullptr),
576 m_derivative2(nullptr),
581 m_flippedDriftDistance(false),
582 m_hitIndex(hitIndex),
583 m_hitOnTrack(hitOnTrack),
585 m_materialEffects(nullptr),
586 m_materialEffectsOwner(false),
587 m_measurementBase(nullptr),
588 m_minEnergyDeposit(0.),
590 m_numericalDerivative(false),
592 m_particleMassSquared(0.),
595 m_position(position),
597 m_radiationThickness(0.),
602 m_scatteringAngle(0.),
603 m_scatteringAngleOffSet(0.),
604 m_secondResidual(0.),
608 m_signedDriftDistance(0.),
622 m_normal =
Amg::Vector3D(m_surface->transform().rotation().col(2));
623 if (m_numberDoF == 1 && std::abs(m_normal.z()) > 0.99 &&
624 std::abs(sinStereo) < 0.5)
626 double cosStereo = std::sqrt(1. - sinStereo * sinStereo);
628 Amg::Vector3D(position(0) * cosStereo + position(1) * sinStereo,
629 -position(0) * sinStereo + position(1) * cosStereo, 0.);
630 (m_sensorDirection) /= m_sensorDirection.perp();
634 m_sensorDirection =
Amg::Vector3D(m_surface->transform().rotation().col(1));
636 m_minimizationDirection =
Amg::Vector3D(m_sensorDirection.cross(m_normal));
640 m_weight = 1. /
sigma;
645 if (m_numberDoF == 2) {
647 m_weight2 = 1. / sigma2;
658 double signedDriftDistance,
double,
660 : m_afterCalo(false),
661 m_alignmentEffects(nullptr),
662 m_alignmentParameter(0),
663 m_alignmentParameter2(0),
665 m_derivative(nullptr),
666 m_derivative2(nullptr),
671 m_flippedDriftDistance(false),
672 m_hitIndex(hitIndex),
673 m_hitOnTrack(hitOnTrack),
675 m_materialEffects(nullptr),
676 m_materialEffectsOwner(false),
677 m_measurementBase(nullptr),
678 m_minEnergyDeposit(0.),
679 m_minimizationDirection{},
682 m_numericalDerivative(
false),
684 m_particleMassSquared(0.),
687 m_position(position),
689 m_radiationThickness(0.),
694 m_scatteringAngle(0.),
695 m_scatteringAngleOffSet(0.),
696 m_secondResidual(0.),
700 m_signedDriftDistance(signedDriftDistance),
706 m_sensorDirection =
Amg::Vector3D(m_surface->transform().rotation().col(2));
710 m_weight = 1. /
sigma;
719 : m_afterCalo(false),
720 m_alignmentEffects(nullptr),
721 m_alignmentParameter(0),
722 m_alignmentParameter2(0),
724 m_derivative(nullptr),
725 m_derivative2(nullptr),
730 m_flippedDriftDistance(false),
732 m_hitOnTrack(nullptr),
734 m_materialEffects(nullptr),
735 m_materialEffectsOwner(false),
736 m_measurementBase(nullptr),
737 m_minEnergyDeposit(0.),
738 m_minimizationDirection{},
741 m_numericalDerivative(
false),
744 m_particleMassSquared(0.),
747 m_position(perigee.associatedSurface().center()),
749 m_radiationThickness(0.),
754 m_scatteringAngle(0.),
755 m_scatteringAngleOffSet(0.),
756 m_secondResidual(0.),
761 m_signedDriftDistance(0.),
763 m_surface(&perigee.associatedSurface()),
768 m_sensorDirection =
Amg::Vector3D(m_surface->transform().rotation().col(2));
771 if (perigee.covariance() && !m_outlier) {
778 double ptInv0 = 1. /
momentum.perp();
779 double cosPhi = ptInv0 *
momentum(0);
780 double sinPhi = ptInv0 *
momentum(1);
782 ptInv0 *= perigee.charge();
799 covariance(
row, 4) = covariance(4,
row);
802 covariance.inverse();
804 JacobianCotThetaPtToThetaP jacobian(
cotTheta, ptInv0);
806 Amg::MatrixX(jacobian * covariance * jacobian.transpose());
813 : m_afterCalo(false),
814 m_alignmentEffects(nullptr),
815 m_alignmentParameter(0),
816 m_alignmentParameter2(0),
818 m_derivative(nullptr),
819 m_derivative2(nullptr),
824 m_flippedDriftDistance(false),
826 m_hitOnTrack(nullptr),
828 m_materialEffects(nullptr),
829 m_materialEffectsOwner(false),
830 m_measurementBase(nullptr),
831 m_minEnergyDeposit(0.),
832 m_minimizationDirection{},
835 m_numericalDerivative(
false),
837 m_particleMassSquared(0.),
840 m_position(position),
842 m_radiationThickness(0.),
847 m_scatteringAngle(0.),
848 m_scatteringAngleOffSet(0.),
849 m_secondResidual(0.),
854 m_signedDriftDistance(0.),
858 m_weight(1. /
sigma),
871 const std::optional<TrackSurfaceIntersection>&
value) {
876 log <<
" residual 1........2 r phi z"
877 <<
" sigma 1.......2 energy energyLoss scatteringAngle "
884 log <<
m_type << std::setiosflags(std::ios::fixed);
888 log << std::setw(9) << std::setprecision(3) << *(
m_residual + 1);
899 log << std::setw(9) << std::setprecision(3) << *
m_residual <<
" outlier ";
906 log << std::setw(10) << std::setprecision(1) <<
position.perp()
907 << std::setw(9) << std::setprecision(4) <<
position.phi() << std::setw(10)
908 << std::setprecision(1) <<
position(2);
911 log << std::setw(13) << std::setprecision(3) << 1. /
m_weight;
913 log << std::setw(8) << std::setprecision(3) << 1. /
m_weight2;
915 log <<
"(" << std::setw(7) << std::setprecision(3)
921 log << std::setw(33) << std::setprecision(3)
924 if (m_type < barrelInert || m_scatteringAngle > 0.) {
928 log << std::setw(16) << std::setprecision(6) << totScat << std::setw(13)
932 log << std::setw(13) << std::setprecision(3) << 1. /
m_weight;
934 log << std::setw(8) << std::setprecision(3) << 1. /
m_weight2;
938 log << std::setw(13) << std::setprecision(3)
942 log << std::setw(33) << std::setprecision(3)
944 log << std::setw(12) << std::setprecision(4)
974 double totalRadiationThickness) {
989 angle_iPat * angle_iPat);
1022 double sigma = std::sqrt(