ATLAS Offline Software
FitMeasurement.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /***************************************************************************
6  for any measurement type (cluster, drift circle or material)
7  stores the quantities needed during track fitting
8  i.e. position, surface, weights, intersection, derivatives, residual etc
9  ***************************************************************************/
10 
12 
13 #include <cmath>
14 #include <iomanip>
15 #include <iostream>
16 
18 #include "GaudiKernel/MsgStream.h"
19 #include "GaudiKernel/SystemOfUnits.h"
33 #include "TrkSurfaces/Surface.h"
37 
38 namespace Trk {
39 
40 // MeasurementBase
41 FitMeasurement::FitMeasurement(int hitIndex, HitOnTrack* hitOnTrack,
42  const MeasurementBase* measurementBase)
43  : m_afterCalo(false),
44  m_alignmentEffects(nullptr),
45  m_alignmentParameter(0),
46  m_alignmentParameter2(0),
47  m_betaSquared(1.),
48  m_derivative(nullptr),
49  m_derivative2(nullptr),
50  m_derivativeRow(-1),
51  m_d0(0.),
52  m_energyLoss(0.),
53  m_firstParameter(0),
54  m_flippedDriftDistance(false),
55  m_hitIndex(hitIndex),
56  m_hitOnTrack(hitOnTrack),
57  m_lastParameter(0),
58  m_materialEffects(nullptr),
59  m_materialEffectsOwner(false),
60  m_measurementBase(measurementBase),
61  m_minEnergyDeposit(0.),
62  m_minimizationDirection{},
63  m_normal{},
64  m_numberDoF(measurementBase->localCovariance().cols()),
65  m_numericalDerivative(false),
66  m_outlier(false),
67  m_particleMassSquared(0.),
68  m_perigee{},
69  m_perigeeWeight{},
70  m_position(measurementBase->associatedSurface().center()),
71  m_qOverP(0.),
72  m_radiationThickness(0.),
73  m_residual(0),
74  m_scaleFactor(0.),
75  m_scatterPhi(0.),
76  m_scatterTheta(0.),
77  m_scatteringAngle(0.),
78  m_scatteringAngleOffSet(0.),
79  m_secondResidual(0.),
80  m_sensorDirection{},
81  m_sigma(0.),
82  m_sigmaMinus(0.),
83  m_sigmaPlus(0.),
84  m_signedDriftDistance(0.),
85  m_status(0),
86  m_surface(&measurementBase->associatedSurface()),
88  m_weight(1.),
89  m_weight2(1.) {
90  double sigma = 0.;
91  if (m_numberDoF > 0)
92  sigma = Amg::error(measurementBase->localCovariance(), locX);
93  double sigma2 = 0.;
94  if (m_numberDoF > 1)
95  sigma2 = Amg::error(measurementBase->localCovariance(), locY);
96 
97  // remaining data according to surface (normal, sensor, minimization
98  // directions etc)
99  if (dynamic_cast<const PlaneSurface*>(m_surface)) {
100  m_normal = Amg::Vector3D(m_surface->transform().rotation().col(2));
101  const Amg::Vector3D posptr =
102  m_surface->localToGlobal(measurementBase->localParameters());
103  m_position = posptr;
104 
105  // special case to get sensor for endcap trapezoids (discs represented as
106  // planes) thus sensor direction rotates according to localX - rather than
107  // being parallel and there is only 1 'real' degree of freedom for
108  // covariance
109  if (m_numberDoF == 2 &&
110  dynamic_cast<const TrapezoidBounds*>(&m_surface->bounds())) {
111  m_numberDoF = 1;
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));
118  // used by obsolete scaling
119  // double lengthSq = 0.5*(sum + diff);
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;
125  if (term > -0.5) {
126  sinStereo = std::sqrt(0.5 + term);
127  if (ab * m_normal.z() < 0.)
128  sinStereo = -sinStereo;
129  }
130 
131  const Amg::Vector3D& axis = m_surface->transform().rotation().col(1);
132  m_sensorDirection =
133  Amg::Vector3D(axis(0) * cosStereo + axis(1) * sinStereo,
134  axis(1) * cosStereo - axis(0) * sinStereo, axis(2));
135  } else {
136  m_sensorDirection =
137  Amg::Vector3D(m_surface->transform().rotation().col(1));
138  }
139 
140  m_minimizationDirection = Amg::Vector3D(m_sensorDirection.cross(m_normal));
141  if (m_numberDoF == 2)
143  } else if (dynamic_cast<const StraightLineSurface*>(m_surface)) {
144  // StraightLines can be driftCircles or pseudoMeasurements along the wire
145  if (!measurementBase->localParameters().contains(locY)) {
146  // driftCircles will eventually have sagged surfaces.
147  // then may need something like:
148  // get z-along wire from (globPos-center).dot(sensor)
149  // then loc3D = (0,0,zalongwire)
150  // position = surface->transform() * loc3D
151  // but probably just same as using centre with wire direction ...
152  m_sensorDirection =
153  Amg::Vector3D(m_surface->transform().rotation().col(2));
154  m_signedDriftDistance = measurementBase->localParameters()[driftRadius];
156  } else {
157  // pseudomeasurement - minimize along wire direction
158  m_minimizationDirection =
159  Amg::Vector3D(m_surface->transform().rotation().col(2));
160  double mag = m_surface->center().mag();
161  m_normal = mag > 1e-6 ? Amg::Vector3D(m_surface->center() / mag)
162  : Amg::Vector3D(m_surface->normal());
163  m_position = measurementBase->globalPosition();
164  m_sensorDirection =
165  Amg::Vector3D(m_normal.cross(m_minimizationDirection));
166  m_signedDriftDistance = 0.;
168  }
169  }
170 
171  const PerigeeSurface* perigee =
172  dynamic_cast<const PerigeeSurface*>(m_surface);
173  if (perigee) {
174  m_position = m_surface->center();
175  m_sensorDirection = Amg::Vector3D(0., 0., 1.);
177  if (m_numberDoF == 2)
178  m_type = vertex;
179  }
180 
181  // there are no measurements from Atlas detectors of cylinder or disc type
182  // const CylinderSurface* cylinder = dynamic_cast<const
183  // CylinderSurface*>(m_surface); if (cylinder)
184  // {
185  // }
186 
187  // const DiscSurface* disc = dynamic_cast<const
188  // DiscSurface*>(m_surface); if (disc)
189  // {
190  // }
191 
192  // add protection against junk input
193  if (sigma > 0.)
194  m_weight = 1. / sigma;
195  if (sigma2 > 0.)
196  m_weight2 = 1. / sigma2;
197 }
198 
199 // MaterialEffectsBase constructor
200 FitMeasurement::FitMeasurement(const MaterialEffectsBase* materialEffects,
201  double particleMass,
202  const Amg::Vector3D& position, double qOverP,
203  bool calo)
204  : m_afterCalo(false),
205  m_alignmentEffects(nullptr),
206  m_alignmentParameter(0),
207  m_alignmentParameter2(0),
208  m_betaSquared(1.),
209  m_derivative(nullptr),
210  m_derivative2(nullptr),
211  m_derivativeRow(-1),
212  m_d0(0.),
213  m_energyLoss(0.),
214  m_firstParameter(0),
215  m_flippedDriftDistance(false),
216  m_hitIndex(0),
217  m_hitOnTrack(nullptr),
218  m_lastParameter(0),
219  m_materialEffects(materialEffects),
220  m_materialEffectsOwner(false),
221  m_measurementBase(nullptr),
222  m_minEnergyDeposit(0.),
223  m_minimizationDirection{},
224  m_normal{},
225  m_numberDoF(0),
226  m_numericalDerivative(false),
227  m_outlier(false),
228  m_particleMassSquared(particleMass * particleMass),
229  m_perigee{},
230  m_perigeeWeight{},
231  m_position(position),
232  m_qOverP(qOverP),
233  m_radiationThickness(materialEffects->thicknessInX0()),
234  m_residual(0),
235  m_scaleFactor(0.),
236  m_scatterPhi(0.),
237  m_scatterTheta(0.),
238  m_scatteringAngle(0.),
239  m_scatteringAngleOffSet(0.),
240  m_secondResidual(0.),
241  m_sensorDirection{},
242  m_sigma(0.),
243  m_sigmaMinus(0.),
244  m_sigmaPlus(0.),
245  m_signedDriftDistance(0.),
246  m_status(0),
247  m_surface(&materialEffects->associatedSurface()),
249  m_weight(0.),
250  m_weight2(0.) {
251  if (dynamic_cast<const CylinderSurface*>(m_surface) ||
252  std::abs(m_surface->normal()(2)) < 0.5)
254 
255  if (calo)
257 
258  // set any energy loss
259  const EnergyLoss* energyLoss = nullptr;
260  const ScatteringAngles* scattering = nullptr;
261  const MaterialEffectsOnTrack* meot =
262  dynamic_cast<const MaterialEffectsOnTrack*>(materialEffects);
263  if (meot) {
264  energyLoss = meot->energyLoss();
265  scattering = meot->scatteringAngles();
266  }
267 
268  if (energyLoss) {
269  // note: EDM defines energy loss to be negative
270  // but here take opposite as this convention is not accepted by
271  // CaloEnergy clients (exception to take verbatim for CaloEnergy)
272  // m_energyLoss = -energyLoss->deltaE();
273  m_energyLoss = std::abs(energyLoss->deltaE());
274 
275  // calo energy deposit treated as a single pure energy loss measurement,
276  // with fit taking error into account
277  if (calo && !scattering && energyLoss->sigmaDeltaE() > 0.) {
278  m_energyLoss = energyLoss->deltaE();
279  m_numberDoF = 1;
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;
290  }
291  }
292  }
293 
294  // initialize scattering angles
295  if (scattering) {
296  m_scatterPhi = scattering->deltaPhi();
297  m_scatterTheta = scattering->deltaTheta();
298  if (calo)
299  m_scatteringAngleOffSet = scattering->sigmaDeltaTheta();
300  // if(calo) std::cout << " Calorimeter scaterrer with sigmaDeltaPhi "
301  // << scattering->sigmaDeltaPhi() << " sigmaDeltaTheta " <<
302  // scattering->sigmaDeltaTheta() << std::endl;
303  }
304 }
305 
306 // constructor creating MaterialEffects for merged scatterers
307 FitMeasurement::FitMeasurement(double radiationThickness, double deltaE,
308  double particleMass,
309  const Amg::Vector3D& position,
310  const Amg::Vector3D& direction, double qOverP,
311  const Surface* surface)
312  : m_afterCalo(false),
313  m_alignmentEffects(nullptr),
314  m_alignmentParameter(0),
315  m_alignmentParameter2(0),
316  m_betaSquared(1.),
317  m_derivative(nullptr),
318  m_derivative2(nullptr),
319  m_derivativeRow(-1),
320  m_d0(0.),
321  m_energyLoss(-deltaE),
322  m_firstParameter(0),
323  m_flippedDriftDistance(false),
324  m_hitIndex(0),
325  m_hitOnTrack(nullptr),
326  m_lastParameter(0),
327  m_materialEffects(nullptr),
328  m_materialEffectsOwner(true),
329  m_measurementBase(nullptr),
330  m_minEnergyDeposit(0.),
331  m_minimizationDirection{},
332  m_normal{},
333  m_numberDoF(0),
334  m_numericalDerivative(false),
335  m_outlier(false),
336  m_particleMassSquared(particleMass * particleMass),
337  m_perigee{},
338  m_perigeeWeight{},
339  m_position(position),
340  m_qOverP(qOverP),
341  m_radiationThickness(radiationThickness),
342  m_residual(0),
343  m_scaleFactor(0.),
344  m_scatterPhi(0.),
345  m_scatterTheta(0.),
346  m_scatteringAngle(0.),
347  m_scatteringAngleOffSet(0.),
348  m_secondResidual(0.),
349  m_sensorDirection{},
350  m_sigma(0.),
351  m_sigmaMinus(0.),
352  m_sigmaPlus(0.),
353  m_signedDriftDistance(0.),
354  m_status(0),
355  m_surface(surface),
357  m_weight(0.),
358  m_weight2(0.) {
359  // plane surface with normal along input direction
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)
365  if (!m_surface) {
366  CurvilinearUVT uvt(direction);
367  m_surface = new PlaneSurface(position, uvt);
368  }
369 
370  // create MaterialEffects
371  std::bitset<MaterialEffectsBase::NumberOfMaterialEffectsTypes> typeMaterial;
372  if (deltaE != 0.)
373  typeMaterial.set(MaterialEffectsBase::EnergyLossEffects);
374  auto energyLoss = std::make_unique<EnergyLoss>(deltaE, 0., 0., 0.);
375  m_materialEffects = new MaterialEffectsOnTrack(
376  radiationThickness, std::move(energyLoss), *m_surface, typeMaterial);
377  if (!surface)
378  delete m_surface;
379  m_surface = &m_materialEffects->associatedSurface();
380 
381  m_intersection[FittedTrajectory] =
382  std::make_optional<TrackSurfaceIntersection>(position, direction, 0.);
383 }
384 
385 // constructor for adding (mis-)alignment effects
386 FitMeasurement::FitMeasurement(const AlignmentEffectsOnTrack* alignmentEffects,
387  const Amg::Vector3D& direction,
388  const Amg::Vector3D& position)
389  : m_afterCalo(false),
390  m_alignmentEffects(alignmentEffects),
391  m_alignmentParameter(0),
392  m_alignmentParameter2(0),
393  m_betaSquared(1.),
394  m_derivative(nullptr),
395  m_derivative2(nullptr),
396  m_derivativeRow(-1),
397  m_d0(0.),
398  m_energyLoss(0.),
399  m_firstParameter(0),
400  m_flippedDriftDistance(false),
401  m_hitIndex(0),
402  m_hitOnTrack(nullptr),
403  m_lastParameter(0),
404  m_materialEffects(nullptr),
405  m_materialEffectsOwner(false),
406  m_measurementBase(nullptr),
407  m_minEnergyDeposit(0.),
408  m_minimizationDirection{},
409  m_normal{},
410  m_numberDoF(2),
411  m_numericalDerivative(false),
412  m_outlier(false),
413  m_particleMassSquared(0.),
414  m_perigee{},
415  m_perigeeWeight{},
416  m_position(position),
417  m_qOverP(0.),
418  m_radiationThickness(0.),
419  m_residual(0),
420  m_scaleFactor(0.),
421  m_scatterPhi(alignmentEffects->deltaAngle()),
422  m_scatterTheta(alignmentEffects->deltaTranslation()),
423  m_scatteringAngle(0.),
424  m_scatteringAngleOffSet(0.),
425  m_secondResidual(0.),
426  m_sensorDirection{},
427  m_sigma(0.),
428  m_sigmaMinus(alignmentEffects->sigmaDeltaAngle()),
429  m_sigmaPlus(alignmentEffects->sigmaDeltaTranslation()),
430  m_signedDriftDistance(0.),
431  m_status(0),
432  m_surface(&alignmentEffects->associatedSurface()),
433  m_type(alignment),
434  m_weight(1.),
435  m_weight2(1.) {
436  // set weights
437  if (m_sigmaMinus)
438  m_weight = 1. / m_sigmaMinus;
439  if (m_sigmaPlus)
440  m_weight2 = 1. / m_sigmaPlus;
441 
442  m_intersection[FittedTrajectory] =
443  std::make_optional<TrackSurfaceIntersection>(position, direction, 0.);
444 }
445 
446 // constructor creating placeholder Surface for delimiting material aggregation
447 FitMeasurement::FitMeasurement(const TrackSurfaceIntersection& intersection,
448  double shift)
449  : m_afterCalo(false),
450  m_alignmentEffects(nullptr),
451  m_alignmentParameter(0),
452  m_alignmentParameter2(0),
453  m_betaSquared(1.),
454  m_derivative(nullptr),
455  m_derivative2(nullptr),
456  m_derivativeRow(-1),
457  m_d0(0.),
458  m_energyLoss(0.),
459  m_firstParameter(0),
460  m_flippedDriftDistance(false),
461  m_hitIndex(0),
462  m_hitOnTrack(nullptr),
463  m_lastParameter(0),
464  m_materialEffects(nullptr),
465  m_materialEffectsOwner(false),
466  m_measurementBase(nullptr),
467  m_minEnergyDeposit(0.),
468  m_minimizationDirection{},
469  m_normal{},
470  m_numberDoF(0),
471  m_numericalDerivative(false),
472  m_outlier(false),
473  m_particleMassSquared(0.),
474  m_perigee{},
475  m_perigeeWeight{},
476  m_position(intersection.position()),
477  m_qOverP(0.),
478  m_radiationThickness(0.),
479  m_residual(0),
480  m_scaleFactor(0.),
481  m_scatterPhi(0.),
482  m_scatterTheta(0.),
483  m_scatteringAngle(0.),
484  m_scatteringAngleOffSet(0.),
485  m_secondResidual(0.),
486  m_sensorDirection{},
487  m_sigma(0.),
488  m_sigmaMinus(0.),
489  m_sigmaPlus(0.),
490  m_signedDriftDistance(0.),
491  m_status(0),
492  m_surface(nullptr),
494  m_weight(0.),
495  m_weight2(0.) {
496  // plane surface with normal along input direction and shift wrt position
497  Amg::Vector3D offset = intersection.direction() * shift;
498  m_position += offset;
499  CurvilinearUVT uvt(intersection.direction());
500  m_surface = new PlaneSurface(m_position, uvt);
501 
502  m_intersection[FittedTrajectory] = std::make_optional<TrackSurfaceIntersection>(
503  m_position, intersection.direction(), 0.);
504 }
505 
506 // other TrackStateOnSurface types
507 FitMeasurement::FitMeasurement(const TrackStateOnSurface& TSOS)
508  : m_afterCalo(false),
509  m_alignmentEffects(nullptr),
510  m_alignmentParameter(0),
511  m_alignmentParameter2(0),
512  m_betaSquared(1.),
513  m_derivative(nullptr),
514  m_derivative2(nullptr),
515  m_derivativeRow(-1),
516  m_d0(0.),
517  m_energyLoss(0.),
518  m_firstParameter(0),
519  m_flippedDriftDistance(false),
520  m_hitIndex(0),
521  m_hitOnTrack(nullptr),
522  m_lastParameter(0),
523  m_materialEffects(nullptr),
524  m_materialEffectsOwner(false),
525  m_measurementBase(nullptr),
526  m_minEnergyDeposit(0.),
527  m_minimizationDirection{},
528  m_normal{},
529  m_numberDoF(0),
530  m_numericalDerivative(false),
531  m_outlier(false),
532  m_particleMassSquared(0.),
533  m_perigee{},
534  m_perigeeWeight{},
535  m_position(TSOS.trackParameters()->position()),
536  m_qOverP(0.),
537  m_radiationThickness(0.),
538  m_residual(0),
539  m_scaleFactor(0.),
540  m_scatterPhi(0.),
541  m_scatterTheta(0.),
542  m_scatteringAngle(0.),
543  m_scatteringAngleOffSet(0.),
544  m_secondResidual(0.),
545  m_sensorDirection{},
546  m_sigma(0.),
547  m_sigmaMinus(0.),
548  m_sigmaPlus(0.),
549  m_signedDriftDistance(0.),
550  m_status(0),
551  m_surface(&TSOS.trackParameters()->associatedSurface()),
552  m_type(hole),
553  m_weight(0.),
554  m_weight2(0.) {
555  // set type if necessary (hole = default)
556  if (TSOS.type(TrackStateOnSurface::Outlier)) {
557  m_outlier = true;
558  }
559  // else if (TSOS.type(TrackStateOnSurface::Perigee))
560  // {
561  // m_type = MSperigee;
562  // }
563 }
564 
565 // SiliconTracker constructor (from iPatRec)
566 FitMeasurement::FitMeasurement(int hitIndex, HitOnTrack* hitOnTrack,
567  const Amg::Vector3D& position, double sigma,
568  double sigma2, double sinStereo, int status,
569  const Surface* surface, MeasurementType type)
570  : m_afterCalo(false),
571  m_alignmentEffects(nullptr),
572  m_alignmentParameter(0),
573  m_alignmentParameter2(0),
574  m_betaSquared(1.),
575  m_derivative(nullptr),
576  m_derivative2(nullptr),
577  m_derivativeRow(-1),
578  m_d0(0.),
579  m_energyLoss(0.),
580  m_firstParameter(0),
581  m_flippedDriftDistance(false),
582  m_hitIndex(hitIndex),
583  m_hitOnTrack(hitOnTrack),
584  m_lastParameter(0),
585  m_materialEffects(nullptr),
586  m_materialEffectsOwner(false),
587  m_measurementBase(nullptr),
588  m_minEnergyDeposit(0.),
589  m_numberDoF(1),
590  m_numericalDerivative(false),
591  m_outlier(false),
592  m_particleMassSquared(0.),
593  m_perigee{},
594  m_perigeeWeight{},
595  m_position(position),
596  m_qOverP(0.),
597  m_radiationThickness(0.),
598  m_residual(0),
599  m_scaleFactor(0.),
600  m_scatterPhi(0.),
601  m_scatterTheta(0.),
602  m_scatteringAngle(0.),
603  m_scatteringAngleOffSet(0.),
604  m_secondResidual(0.),
605  m_sigma(0.),
606  m_sigmaMinus(0.),
607  m_sigmaPlus(0.),
608  m_signedDriftDistance(0.),
609  m_status(status),
610  m_surface(surface),
611  m_type(type),
612  m_weight(1.),
613  m_weight2(1.) {
614  // pixel has 2-D measurement
615  if (type == pixelCluster) {
616  m_numberDoF = 2;
617  }
618 
619  // special treatment for projective trapezoidal chambers in the endcap
620  // take sensorDirection from stereo angle as I don't understand Surface axes
621  // in this case
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) // end-cap projective geometry
625  {
626  double cosStereo = std::sqrt(1. - sinStereo * sinStereo);
627  m_sensorDirection =
628  Amg::Vector3D(position(0) * cosStereo + position(1) * sinStereo,
629  -position(0) * sinStereo + position(1) * cosStereo, 0.);
630  (m_sensorDirection) /= m_sensorDirection.perp();
631  } else // otherwise chambers have parallel strips with sensor direction =
632  // appropriate module axis
633  {
634  m_sensorDirection = Amg::Vector3D(m_surface->transform().rotation().col(1));
635  }
636  m_minimizationDirection = Amg::Vector3D(m_sensorDirection.cross(m_normal));
637 
638  // add protection against junk input
639  if (sigma > 0.) {
640  m_weight = 1. / sigma;
641  } else {
642  m_numberDoF = 0;
643  m_outlier = true;
644  }
645  if (m_numberDoF == 2) {
646  if (sigma2 > 0.) {
647  m_weight2 = 1. / sigma2;
648  } else {
649  m_numberDoF = 0;
650  m_outlier = true;
651  }
652  }
653 }
654 
655 // drift constructor (TRT from iPatRec)
656 FitMeasurement::FitMeasurement(int hitIndex, HitOnTrack* hitOnTrack,
657  const Amg::Vector3D& position, double sigma,
658  double signedDriftDistance, double,
659  const Surface* surface)
660  : m_afterCalo(false),
661  m_alignmentEffects(nullptr),
662  m_alignmentParameter(0),
663  m_alignmentParameter2(0),
664  m_betaSquared(1.),
665  m_derivative(nullptr),
666  m_derivative2(nullptr),
667  m_derivativeRow(-1),
668  m_d0(0.),
669  m_energyLoss(0.),
670  m_firstParameter(0),
671  m_flippedDriftDistance(false),
672  m_hitIndex(hitIndex),
673  m_hitOnTrack(hitOnTrack),
674  m_lastParameter(0),
675  m_materialEffects(nullptr),
676  m_materialEffectsOwner(false),
677  m_measurementBase(nullptr),
678  m_minEnergyDeposit(0.),
679  m_minimizationDirection{},
680  m_normal{},
681  m_numberDoF(1),
682  m_numericalDerivative(false),
683  m_outlier(false),
684  m_particleMassSquared(0.),
685  m_perigee{},
686  m_perigeeWeight{},
687  m_position(position),
688  m_qOverP(0.),
689  m_radiationThickness(0.),
690  m_residual(0),
691  m_scaleFactor(0.),
692  m_scatterPhi(0.),
693  m_scatterTheta(0.),
694  m_scatteringAngle(0.),
695  m_scatteringAngleOffSet(0.),
696  m_secondResidual(0.),
697  m_sigma(0.),
698  m_sigmaMinus(0.),
699  m_sigmaPlus(0.),
700  m_signedDriftDistance(signedDriftDistance),
701  m_status(0),
702  m_surface(surface),
704  m_weight(1.),
705  m_weight2(1.) {
706  m_sensorDirection = Amg::Vector3D(m_surface->transform().rotation().col(2));
707 
708  // add protection against junk input
709  if (sigma > 0.) {
710  m_weight = 1. / sigma;
711  } else {
712  m_numberDoF = 0;
713  m_outlier = true;
714  }
715 }
716 
717 // perigee (TrackParameters) constructor
719  : m_afterCalo(false),
720  m_alignmentEffects(nullptr),
721  m_alignmentParameter(0),
722  m_alignmentParameter2(0),
723  m_betaSquared(1.),
724  m_derivative(nullptr),
725  m_derivative2(nullptr),
726  m_derivativeRow(-1),
727  m_d0(0.),
728  m_energyLoss(0.),
729  m_firstParameter(0),
730  m_flippedDriftDistance(false),
731  m_hitIndex(0),
732  m_hitOnTrack(nullptr),
733  m_lastParameter(0),
734  m_materialEffects(nullptr),
735  m_materialEffectsOwner(false),
736  m_measurementBase(nullptr),
737  m_minEnergyDeposit(0.),
738  m_minimizationDirection{},
739  m_normal{},
740  m_numberDoF(0),
741  m_numericalDerivative(false),
742  m_outlier(true), // use base class for additional trackParameters at
743  // detector boundary
744  m_particleMassSquared(0.),
745  m_perigee{},
746  m_perigeeWeight{},
747  m_position(perigee.associatedSurface().center()),
748  m_qOverP(0.),
749  m_radiationThickness(0.),
750  m_residual(0),
751  m_scaleFactor(0.),
752  m_scatterPhi(0.),
753  m_scatterTheta(0.),
754  m_scatteringAngle(0.),
755  m_scatteringAngleOffSet(0.),
756  m_secondResidual(0.),
757  m_sensorDirection{},
758  m_sigma(0.),
759  m_sigmaMinus(0.),
760  m_sigmaPlus(0.),
761  m_signedDriftDistance(0.),
762  m_status(0),
763  m_surface(&perigee.associatedSurface()),
765  m_weight(1.),
766  m_weight2(1.) {
767  // perigee axis needed for propagation of fitted parameters
768  m_sensorDirection = Amg::Vector3D(m_surface->transform().rotation().col(2));
769 
770  // is this perigee to be used as a measurement?
771  if (perigee.covariance() && !m_outlier) {
772  // use as measurement
773  m_numberDoF = 5;
774  m_outlier = false;
775 
776  // perigeeParameters as HepVector
777  Amg::Vector3D momentum = perigee.momentum();
778  double ptInv0 = 1. / momentum.perp();
779  double cosPhi = ptInv0 * momentum(0);
780  double sinPhi = ptInv0 * momentum(1);
781  double cotTheta = ptInv0 * momentum(2);
782  ptInv0 *= perigee.charge();
783 
785  parameters(0) = perigee.parameters()[Trk::d0];
786  parameters(1) = perigee.parameters()[Trk::z0];
787  parameters(2) = cosPhi;
788  parameters(3) = sinPhi;
789  parameters(4) = cotTheta;
790  parameters(5) = ptInv0;
791  m_perigee = Amg::VectorX(parameters);
792 
793  // weight = inverse covariance
794  AmgSymMatrix(5) covariance(*perigee.covariance());
795 
796  // convert to internal units (TeV) to avoid rounding
797  for (int row = 0; row < 5; ++row) {
798  covariance(4, row) *= Gaudi::Units::TeV;
799  covariance(row, 4) = covariance(4, row);
800  }
801  covariance(4, 4) *= Gaudi::Units::TeV;
802  covariance.inverse();
803 
804  JacobianCotThetaPtToThetaP jacobian(cotTheta, ptInv0);
805  m_perigeeWeight =
806  Amg::MatrixX(jacobian * covariance * jacobian.transpose());
807  }
808 }
809 
810 // transverseVertex constructor
811 FitMeasurement::FitMeasurement(double d0, const Amg::Vector3D& position,
812  double sigma)
813  : m_afterCalo(false),
814  m_alignmentEffects(nullptr),
815  m_alignmentParameter(0),
816  m_alignmentParameter2(0),
817  m_betaSquared(1.),
818  m_derivative(nullptr),
819  m_derivative2(nullptr),
820  m_derivativeRow(-1),
821  m_d0(d0), // FIXME:: kept for cache tag as d0 is never used anywhere
822  m_energyLoss(0.),
823  m_firstParameter(0),
824  m_flippedDriftDistance(false),
825  m_hitIndex(0),
826  m_hitOnTrack(nullptr),
827  m_lastParameter(0),
828  m_materialEffects(nullptr),
829  m_materialEffectsOwner(false),
830  m_measurementBase(nullptr),
831  m_minEnergyDeposit(0.),
832  m_minimizationDirection{},
833  m_normal{},
834  m_numberDoF(1),
835  m_numericalDerivative(false),
836  m_outlier(false),
837  m_particleMassSquared(0.),
838  m_perigee{},
839  m_perigeeWeight{},
840  m_position(position),
841  m_qOverP(0.),
842  m_radiationThickness(0.),
843  m_residual(0),
844  m_scaleFactor(0.),
845  m_scatterPhi(0.),
846  m_scatterTheta(0.),
847  m_scatteringAngle(0.),
848  m_scatteringAngleOffSet(0.),
849  m_secondResidual(0.),
850  m_sensorDirection(Amg::Vector3D(0., 0., 1.)),
851  m_sigma(0.),
852  m_sigmaMinus(0.),
853  m_sigmaPlus(0.),
854  m_signedDriftDistance(0.),
855  m_status(0),
856  m_surface(new const Trk::PerigeeSurface(position)),
858  m_weight(1. / sigma),
859  m_weight2(1.) {}
860 
861 // destructor
865  delete m_surface;
867  delete m_materialEffects;
868 }
869 
871  const std::optional<TrackSurfaceIntersection>& value) {
873 }
874 
875 void FitMeasurement::printHeading(MsgStream& log) {
876  log << " residual 1........2 r phi z"
877  << " sigma 1.......2 energy energyLoss scatteringAngle "
878  "integral X0"
879  << std::endl;
880 }
881 
882 void FitMeasurement::print(MsgStream& log) const {
884  log << m_type << std::setiosflags(std::ios::fixed);
885  if (numberDoF()) {
886  log << std::setw(9) << std::setprecision(3) << *m_residual;
887  if (m_numberDoF > 1) {
888  log << std::setw(9) << std::setprecision(3) << *(m_residual + 1);
889  } else if (m_alignmentParameter2) {
890  log << " A" << std::setw(1) << m_alignmentParameter << " A"
891  << std::setw(1) << m_alignmentParameter2 << " ";
892  } else if (m_alignmentParameter) {
893  log << " A" << std::setw(1) << m_alignmentParameter << " ";
894  } else {
895  log << " ";
896  }
897  } else {
898  if (isPositionMeasurement()) {
899  log << std::setw(9) << std::setprecision(3) << *m_residual << " outlier ";
900  } else {
901  log << " ";
902  }
905  }
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);
909 
910  if (isPositionMeasurement()) {
911  log << std::setw(13) << std::setprecision(3) << 1. / m_weight;
912  if (m_numberDoF == 2) {
913  log << std::setw(8) << std::setprecision(3) << 1. / m_weight2;
914  } else if (isDrift()) {
915  log << "(" << std::setw(7) << std::setprecision(3)
916  << m_signedDriftDistance << ")";
917  } else {
918  log << " ";
919  }
920  } else if (isScatterer()) {
921  log << std::setw(33) << std::setprecision(3)
922  << 1. / std::abs(m_qOverP * Gaudi::Units::GeV) << std::setw(12)
923  << std::setprecision(4) << m_energyLoss / Gaudi::Units::GeV;
924  if (m_type < barrelInert || m_scatteringAngle > 0.) {
925  double totScat = sqrt(m_scatteringAngle * std::abs(m_qOverP) *
926  m_scatteringAngle * std::abs(m_qOverP) +
928  log << std::setw(16) << std::setprecision(6) << totScat << std::setw(13)
929  << std::setprecision(3) << m_radiationThickness;
930  }
931  } else if (isAlignment()) {
932  log << std::setw(13) << std::setprecision(3) << 1. / m_weight;
933  if (m_numberDoF == 2) {
934  log << std::setw(8) << std::setprecision(3) << 1. / m_weight2;
935  }
936  } else if (isEnergyDeposit()) {
937  if (m_numberDoF) {
938  log << std::setw(13) << std::setprecision(3)
939  << 1. / (m_weight * Gaudi::Units::GeV) << std::setw(32)
940  << std::setprecision(4) << m_energyLoss / Gaudi::Units::GeV;
941  } else {
942  log << std::setw(33) << std::setprecision(3)
943  << 1. / std::abs(m_qOverP * Gaudi::Units::GeV);
944  log << std::setw(12) << std::setprecision(4)
946  }
947  }
948  log << std::endl;
949 }
950 
951 void FitMeasurement::qOverP(double value) {
952  m_qOverP = value;
953 
954  // for scatterer measurements: correct the weight for the local momentum value
957  m_scatteringAngle > 0.) {
958  double pSquare = 1. / (value * value);
959  m_betaSquared = pSquare / (pSquare + m_particleMassSquared);
960  m_weight = std::sqrt(m_betaSquared * pSquare) / m_scatteringAngle;
961  }
962 
964  if (m_weight > 0) {
965  m_weight = sqrt(1. / (1. / m_weight / m_weight +
967  } else {
969  }
970  }
971 }
972 
974  double totalRadiationThickness) {
975  // update the m_scatteringAngleOffSet at initialisation
976 
978  m_qOverP != 0.) {
979  //
980  double angle_iPat = angle * std::abs(m_qOverP);
981 
982  // std::cout << " scatteringAngle type " << m_type << " angle_iPat "
983  // << angle_iPat << " m_scatteringAngleOffSet " <<
984  // m_scatteringAngleOffSet;
985 
986  if (angle_iPat < m_scatteringAngleOffSet) {
989  angle_iPat * angle_iPat);
990  } else {
992  }
993  // std::cout << " corrected m_scatteringAngleOffSet " <<
994  // m_scatteringAngleOffSet << std::endl;
995  }
996 
998  m_radiationThickness = totalRadiationThickness;
999  if (m_type == barrelInert) {
1001  } else if (m_type == endcapInert) {
1003  }
1004 
1005  if (m_qOverP != 0.) {
1006  double pSquare = 1. / (m_qOverP * m_qOverP);
1007  m_betaSquared = pSquare / (pSquare + m_particleMassSquared);
1008  m_weight = std::sqrt(m_betaSquared * pSquare) / m_scatteringAngle;
1010  if (m_weight > 0) {
1011  m_weight =
1012  sqrt(1. / (1. / m_weight / m_weight +
1014  } else {
1016  }
1017  }
1018  }
1019 }
1020 
1022  double sigma = std::sqrt(
1024  m_weight = 1. / sigma;
1026 
1027 } // namespace Trk
TrapezoidBounds.h
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
query_example.row
row
Definition: query_example.py:24
Trk::FitMeasurement::intersection
const TrackSurfaceIntersection & intersection(ExtrapolationType type) const
Definition: FitMeasurement.h:359
EnergyLoss.h
Trk::FitMeasurement::m_betaSquared
double m_betaSquared
Definition: FitMeasurement.h:194
ScatteringAngles.h
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:17
AlignmentEffectsOnTrack.h
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:30
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
Trk::FitMeasurement::numberDoF
int numberDoF(void) const
Definition: FitMeasurement.h:460
perigeeParameters
Definition: ExtrapolatorComparisonTest.h:42
StraightLineSurface.h
Trk::FitMeasurement::printHeading
static void printHeading(MsgStream &log)
Definition: FitMeasurement.cxx:879
MeasurementBase.h
PerigeeSurface.h
Trk::locX
@ locX
Definition: ParamDefs.h:37
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:38
Trk::FitMeasurement::m_intersection
std::array< std::optional< TrackSurfaceIntersection >, ExtrapolationTypes > m_intersection
Definition: FitMeasurement.h:205
Surface.h
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
Trk::materialDelimiter
@ materialDelimiter
Definition: MeasurementType.h:37
Trk::MeasurementType
MeasurementType
Definition: MeasurementType.h:18
Trk::FitMeasurement::isDrift
bool isDrift(void) const
Definition: FitMeasurement.h:373
EventPrimitivesHelpers.h
yodamerge_tmp.axis
list axis
Definition: yodamerge_tmp.py:241
Trk::calorimeterScatterer
@ calorimeterScatterer
Definition: MeasurementType.h:29
Surface
Definition: Trigger/TrigAccel/TrigCudaFitter/src/Surface.h:8
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
Trk::z0
@ z0
Definition: ParamDefs.h:64
athena.value
value
Definition: athena.py:124
Trk::FitMeasurement::~FitMeasurement
~FitMeasurement(void)
Definition: FitMeasurement.cxx:866
python.SystemOfUnits.TeV
int TeV
Definition: SystemOfUnits.py:158
Trk::FitMeasurement::m_materialEffects
const MaterialEffectsBase * m_materialEffects
Definition: FitMeasurement.h:207
Trk::FitMeasurement::m_particleMassSquared
double m_particleMassSquared
Definition: FitMeasurement.h:216
Trk::FitMeasurement::position
const Amg::Vector3D & position(void) const
Definition: FitMeasurement.h:482
Trk::barrelScatterer
@ barrelScatterer
Definition: MeasurementType.h:27
Trk::FitMeasurement::m_sigmaPlus
double m_sigmaPlus
Definition: FitMeasurement.h:232
Trk::FitMeasurement::setSigmaSymmetric
void setSigmaSymmetric(void)
Definition: FitMeasurement.cxx:1025
Trk::FitMeasurement::m_alignmentParameter
unsigned m_alignmentParameter
Definition: FitMeasurement.h:192
intersection
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Definition: compareFlatTrees.cxx:25
Trk::energyDeposit
@ energyDeposit
Definition: MeasurementType.h:32
Trk::FitMeasurement::m_alignmentParameter2
unsigned m_alignmentParameter2
Definition: FitMeasurement.h:193
Trk::stripCluster
@ stripCluster
Definition: MeasurementType.h:23
Trk::transverseVertex
@ transverseVertex
Definition: MeasurementType.h:20
Trk::FitMeasurement::type
MeasurementType type(void) const
Definition: FitMeasurement.h:589
Trk::TrackStateOnSurface::Outlier
@ Outlier
This TSoS contains an outlier, that is, it contains a MeasurementBase/RIO_OnTrack which was not used ...
Definition: TrackStateOnSurface.h:122
MaterialEffectsOnTrack.h
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
m_type
TokenType m_type
the type
Definition: TProperty.cxx:44
GeoPrimitives.h
Trk::FitMeasurement::m_radiationThickness
double m_radiationThickness
Definition: FitMeasurement.h:221
Trk::trapezoidCluster
@ trapezoidCluster
Definition: MeasurementType.h:24
Trk::FitMeasurement::scatteringAngle
void scatteringAngle(double angle, double totalRadiationThickness)
Definition: FitMeasurement.cxx:977
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
Trk::barrelInert
@ barrelInert
Definition: MeasurementType.h:30
Trk::FitMeasurement::m_weight
double m_weight
Definition: FitMeasurement.h:237
Trk::FitMeasurement::m_residual
std::vector< double >::iterator m_residual
Definition: FitMeasurement.h:222
Trk::endcapInert
@ endcapInert
Definition: MeasurementType.h:31
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition: TRTDetectorFactory_Full.cxx:73
Trk::driftRadius
@ driftRadius
trt, straws
Definition: ParamDefs.h:53
CylinderSurface.h
Trk::FitMeasurement::m_weight2
double m_weight2
Definition: FitMeasurement.h:238
Trk::FitMeasurement::print
void print(MsgStream &log) const
Definition: FitMeasurement.cxx:886
Trk::FitMeasurement::m_position
Amg::Vector3D m_position
Definition: FitMeasurement.h:219
Trk::FitMeasurement::m_sigmaMinus
double m_sigmaMinus
Definition: FitMeasurement.h:231
Trk::FitMeasurement::FitMeasurement
FitMeasurement(int hitIndex, HitOnTrack *hitOnTrack, const MeasurementBase *measurementBase)
Definition: FitMeasurement.cxx:45
Trk::FitMeasurement::m_energyLoss
double m_energyLoss
Definition: FitMeasurement.h:199
Trk::FitMeasurement::sigma
double sigma(void) const
Definition: FitMeasurement.h:565
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::FitMeasurement::m_signedDriftDistance
double m_signedDriftDistance
Definition: FitMeasurement.h:233
Trk::ExtrapolationType
ExtrapolationType
Definition: ExtrapolationType.h:18
Trk::FitMeasurement::isEnergyDeposit
bool isEnergyDeposit(void) const
Definition: FitMeasurement.h:377
Amg
Definition of ATLAS Math & Geometry primitives (Amg)
Definition: AmgStringHelpers.h:19
Trk::d0
@ d0
Definition: ParamDefs.h:63
Amg::error
double error(const Amg::MatrixX &mat, int index)
return diagonal error of the matrix caller should ensure the matrix is symmetric and the index is in ...
Definition: EventPrimitivesHelpers.h:40
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition: Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
Trk::FitMeasurement::m_numberDoF
int m_numberDoF
Definition: FitMeasurement.h:213
Trk::FittedTrajectory
@ FittedTrajectory
Definition: ExtrapolationType.h:19
Trk::FitMeasurement::isAlignment
bool isAlignment(void) const
Definition: FitMeasurement.h:364
Trk::FitMeasurement::m_scatteringAngle
double m_scatteringAngle
Definition: FitMeasurement.h:226
Trk::pseudoMeasurement
@ pseudoMeasurement
Definition: MeasurementType.h:26
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::FitMeasurement::m_qOverP
double m_qOverP
Definition: FitMeasurement.h:220
Trk::MaterialEffectsBase::EnergyLossEffects
@ EnergyLossEffects
contains energy loss corrections
Definition: MaterialEffectsBase.h:48
FitMeasurement.h
Trk::pixelCluster
@ pixelCluster
Definition: MeasurementType.h:22
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
Trk::FitMeasurement::m_type
MeasurementType m_type
Definition: FitMeasurement.h:236
Trk::FitMeasurement::m_materialEffectsOwner
bool m_materialEffectsOwner
Definition: FitMeasurement.h:208
PlaneSurface.h
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
JacobianCotThetaPtToThetaP.h
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
convertTimingResiduals.offset
offset
Definition: convertTimingResiduals.py:71
DiscSurface.h
TRT::Track::cotTheta
@ cotTheta
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:65
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Trk::FitMeasurement::m_scatteringAngleOffSet
double m_scatteringAngleOffSet
Definition: FitMeasurement.h:227
Trk::hole
@ hole
Definition: MeasurementType.h:36
Trk::FitMeasurement::isPositionMeasurement
bool isPositionMeasurement(void) const
Definition: FitMeasurement.h:401
merge.status
status
Definition: merge.py:17
RotatedTrapezoidBounds.h
TrackSurfaceIntersection.h
Trk::FitMeasurement::qOverP
double qOverP(void) const
Definition: FitMeasurement.h:486
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
Trk::FitMeasurement::m_measurementBase
const MeasurementBase * m_measurementBase
Definition: FitMeasurement.h:209
Trk::FitMeasurement::isScatterer
bool isScatterer(void) const
Definition: FitMeasurement.h:409
test_AnalysisBaseEventLoopJob.aa
aa
Definition: test_AnalysisBaseEventLoopJob.py:37
Trk::FitMeasurement::m_surface
const Surface * m_surface
Definition: FitMeasurement.h:235
Trk::alignment
@ alignment
Definition: MeasurementType.h:33
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
Trk::endcapScatterer
@ endcapScatterer
Definition: MeasurementType.h:28
Trk::driftCircle
@ driftCircle
Definition: MeasurementType.h:25
TrackStateOnSurface.h