ATLAS Offline Software
MeasurementProcessor.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /***************************************************************************
6  process measurements i.e. extrapolate to surface,
7  compute derivatives and residuals
8  ***************************************************************************/
9 
10 #include "TrkiPatFitterUtils/MeasurementProcessor.h"
11 
12 #include <cmath>
13 #include <iomanip>
14 #include <iostream>
15 
16 #include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
17 #include "GaudiKernel/MsgStream.h"
18 #include "GaudiKernel/SystemOfUnits.h"
19 #include "TrkExInterfaces/IIntersector.h"
20 #include "TrkExInterfaces/IPropagator.h"
21 #include "TrkExUtils/TrackSurfaceIntersection.h"
22 #include "TrkSurfaces/Surface.h"
23 #include "TrkiPatFitterUtils/FitMeasurement.h"
24 #include "TrkiPatFitterUtils/FitParameters.h"
25 #include "TrkiPatFitterUtils/ParameterType.h"
26 namespace Trk {
27 
28 // constructor
29 MeasurementProcessor::MeasurementProcessor(
30  bool asymmetricCaloEnergy, Amg::MatrixX& /*derivativeMatrix*/,
31  const ToolHandle<IIntersector>& intersector,
32  std::vector<FitMeasurement*>& measurements, FitParameters* parameters,
33  const ToolHandle<IIntersector>& rungeKuttaIntersector,
34  const ToolHandle<IPropagator>& stepPropagator, int useStepPropagator)
35  : m_asymmetricCaloEnergy(asymmetricCaloEnergy),
36  m_caloEnergyMeasurement(nullptr),
37  m_cosPhi0(parameters->cosPhi()),
38  m_cosTheta0(parameters->cosTheta()),
39  m_derivQOverP0(0.),
40  m_derivQOverP1(0.),
41  m_energyResidual(0.),
42  m_firstScatteringParameter(parameters->firstScatteringParameter()),
43  // m_havePhiPseudo (false),
44  m_intersector(intersector),
45  m_largeDeltaD0(50. * Gaudi::Units::mm),
46  m_largeDeltaPhi0(0.05),
47  // m_minDistanceForAngle (2.*Gaudi::Units::mm),
48  m_measurements(measurements),
49  m_numericDerivatives(false),
50  m_parameters(parameters),
51  m_phiInstability(false),
52  m_qOverPbeforeCalo(0.),
53  m_qOverPafterCalo(0.),
54  m_rungeKuttaIntersector(rungeKuttaIntersector),
55  m_stepPropagator(stepPropagator),
56  m_useStepPropagator(useStepPropagator),
57  m_sinPhi0(parameters->sinPhi()),
58  m_sinTheta0(parameters->sinTheta()),
59  // m_toroidTurn (0.1),
60  m_vertexIntersect(),
61  m_intersectStartingValue(),
62  m_x0(parameters->position().x()),
63  m_y0(parameters->position().y()),
64  m_z0(parameters->position().z())
65 // m_zInstability (false)
66 {
67  m_alignments.reserve(10);
68  m_intersectStartingValue = m_parameters->intersection();
69  // bool haveDrift = false;
70  int numberPositionMeas = 0;
71  // int numberPseudo = 0;
72  m_scatterers.reserve(100);
73 
74  // Field for StepPropagator
75  m_stepField = Trk::MagneticFieldProperties(Trk::FullField);
76  if (m_useStepPropagator == 2)
77  m_stepField = Trk::MagneticFieldProperties(Trk::FastField);
78 
79  for (auto* m : m_measurements) {
80  if (!m->numberDoF())
81  continue;
82  if (m->isAlignment())
83  m_alignments.push_back(m);
84  // if (m->isDrift()) haveDrift = true;
85  if (m->isPositionMeasurement())
86  ++numberPositionMeas;
87  // if (m->isPseudo()) ++numberPseudo;
88  if (m->isScatterer())
89  m_scatterers.push_back(m);
90 
91  // cache some variables when we fit calo energy deposit
92  if (!m->isEnergyDeposit() || !m_parameters->fitMomentum())
93  continue;
94  m_caloEnergyMeasurement = m;
95  if (numberPositionMeas < 5) {
96  m_asymmetricCaloEnergy = false;
97  m_caloEnergyMeasurement->setSigmaSymmetric();
98  }
99  m_qOverPbeforeCalo = m_parameters->qOverP();
100  m_qOverPafterCalo = m->qOverP();
101  if (m_qOverPbeforeCalo * m_qOverPafterCalo < 0.)
102  m_qOverPafterCalo = m_qOverPbeforeCalo;
103  m_parameters->qOverP1(m_qOverPafterCalo);
104  }
105 
106  // // pseudoMeasurement projectivity constraint needed for some MS tracks
107  // if (numberPseudo && (numberPseudo > 1 ||
108  // m_measurements.front()->isVertex())) m_phiInstability = true;
109 
110  // for numerical derivs
111  for (int typ = 0; typ != ExtrapolationTypes; ++typ) {
112  m_delta[typ] = 0.0001;
113  m_qOverP[typ] = m_parameters->qOverP();
114  }
115 }
116 
117 // destructor
118 MeasurementProcessor::~MeasurementProcessor(void) = default;
119 
120 bool MeasurementProcessor::calculateDerivatives(void) {
121  // extrapolate to all measurement surfaces to compute intersects for momentum
122  // derivatives
123  if (m_parameters->fitMomentum() && !extrapolateToMeasurements(DeltaQOverP0))
124  return false;
125  if (m_parameters->fitEnergyDeposit() &&
126  !extrapolateToMeasurements(DeltaQOverP1))
127  return false;
128 
129  // extrapolate for any other numeric derivatives
130  if (m_numericDerivatives) {
131  double ptInv = std::abs(m_parameters->ptInv0());
132  m_delta[DeltaD0] = 0.010 + 10.0 * ptInv;
133  m_delta[DeltaZ0] = 0.010 + 10.0 * ptInv;
134  m_delta[DeltaPhi0] = 0.0001 + 2.0 * ptInv;
135  m_delta[DeltaTheta0] = 0.0005 + 2.0 * ptInv;
136 
137  TrackSurfaceIntersection intersection = m_vertexIntersect;
138  m_vertexIntersect = TrackSurfaceIntersection(
139  Amg::Vector3D(m_x0 - m_delta[DeltaD0] * m_sinPhi0,
140  m_y0 + m_delta[DeltaD0] * m_cosPhi0, m_z0),
141  intersection.direction(), 0.);
142  if (!extrapolateToMeasurements(DeltaD0)) {
143  return false;
144  }
145  m_vertexIntersect = TrackSurfaceIntersection(
146  Amg::Vector3D(m_x0, m_y0, m_z0 + m_delta[DeltaZ0]),
147  intersection.direction(), 0.);
148  if (!extrapolateToMeasurements(DeltaZ0)) {
149  return false;
150  }
151  double sinTheta = intersection.direction().perp();
152  double delCF = 1. - (0.5 * m_delta[DeltaPhi0] * m_delta[DeltaPhi0]);
153  Amg::Vector3D direction(
154  sinTheta * (m_cosPhi0 * delCF - m_sinPhi0 * m_delta[DeltaPhi0]),
155  sinTheta * (m_sinPhi0 * delCF + m_cosPhi0 * m_delta[DeltaPhi0]),
156  intersection.direction().z());
157  m_vertexIntersect =
158  TrackSurfaceIntersection(intersection.position(), direction, 0.);
159  if (!extrapolateToMeasurements(DeltaPhi0)) {
160  return false;
161  }
162  double cotTheta =
163  m_vertexIntersect.direction().z() / sinTheta + m_delta[DeltaTheta0];
164  sinTheta = 1. / std::sqrt(1. + cotTheta * cotTheta);
165  Amg::Vector3D directionTheta(sinTheta * m_cosPhi0, sinTheta * m_sinPhi0,
166  sinTheta * cotTheta);
167  m_vertexIntersect =
168  TrackSurfaceIntersection(intersection.position(), directionTheta, 0.);
169  if (!extrapolateToMeasurements(DeltaTheta0)) {
170  return false;
171  }
172  m_vertexIntersect = intersection;
173  }
174 
175  // loop over measurements to compute derivatives:
176  for (auto* m : m_measurements) {
177  // strip detector types
178  if (m->isCluster()) {
179  clusterDerivatives(m->numberDoF(), *m);
180  } else if (m->isDrift() || m->isPerigee() ||
181  m->isPseudo()) // else drift circles
182  {
183  const TrackSurfaceIntersection& intersection =
184  m->intersection(FittedTrajectory);
185  Amg::Vector3D driftDirection =
186  Amg::Vector3D(m->sensorDirection().cross(intersection.direction()));
187  driftDirection = driftDirection.unit();
188  m->minimizationDirection(driftDirection);
189  driftDerivatives(m->numberDoF(), *m);
190  } else if (m->isVertex()) {
191  const TrackSurfaceIntersection& intersection =
192  m->intersection(FittedTrajectory);
193  Amg::Vector3D minimizationDirection =
194  Amg::Vector3D(m->sensorDirection().cross(intersection.direction()));
195  minimizationDirection = minimizationDirection.unit();
196  m->minimizationDirection(minimizationDirection);
197  m->derivative(D0,
198  m->weight() * (m_cosPhi0 * m->minimizationDirection().y() -
199  m_sinPhi0 * m->minimizationDirection().x()));
200  if (m->is2Dimensional())
201  m->derivative2(Z0, m->weight2() * m->sensorDirection().z());
202  } else if (!m->numberDoF()) {
203  continue;
204  } else if (m->isScatterer()) {
205  unsigned param = m->lastParameter() - 2;
206  m->derivative(
207  param,
208  m->weight() * m->intersection(FittedTrajectory).direction().perp());
209  m->derivative2(++param, m->weight());
210  } else if (m->isAlignment()) {
211  unsigned param = m->firstParameter();
212  m->derivative(param, m->weight());
213  m->derivative2(++param, m->weight2());
214  } else if (m->isEnergyDeposit()) {
215  double E0 = 1. / std::abs(m_qOverPbeforeCalo);
216  double E1 = 1. / std::abs(m_qOverPafterCalo);
217  double weight = m->weight();
218  m->derivative(QOverP0,
219  weight * E0 / (m_qOverPbeforeCalo * Gaudi::Units::TeV));
220  m->derivative(QOverP1,
221  -weight * E1 / (m_qOverPafterCalo * Gaudi::Units::TeV));
222  } else {
223  // report any mysteries ...
224  continue;
225  }
226  }
227 
228  // flip step sign of numerical derivatives for next iteration
229  if (m_numericDerivatives) {
230  for (int typ = 1; typ != ExtrapolationTypes; ++typ)
231  m_delta[typ] = -m_delta[typ];
232  }
233 
234  return true;
235 }
236 
237 bool MeasurementProcessor::calculateFittedTrajectory(int /*iteration*/) {
238  // cache frequently accessed parameters
239  m_cosPhi0 = m_parameters->cosPhi();
240  m_cosTheta0 = m_parameters->cosTheta();
241  m_sinPhi0 = m_parameters->sinPhi();
242  m_sinTheta0 = m_parameters->sinTheta();
243  m_x0 = m_parameters->position().x();
244  m_y0 = m_parameters->position().y();
245  m_z0 = m_parameters->position().z();
246 
247  // start with intersection placed at perigee
248  m_vertexIntersect = m_parameters->intersection();
249 
250  // increments for momentum derivatives (Mev^-1 : X-over at 30GeV)
251  double floor = 0.000000030;
252  double fraction = 0.0001;
253  if (m_parameters->qOverP() < 0.) {
254  m_delta[DeltaQOverP0] = floor - fraction * m_parameters->qOverP();
255  m_delta[DeltaQOverP1] = m_delta[DeltaQOverP0];
256  } else {
257  m_delta[DeltaQOverP0] = -floor - fraction * m_parameters->qOverP();
258  m_delta[DeltaQOverP1] = m_delta[DeltaQOverP0];
259  }
260  if (m_parameters->fitEnergyDeposit())
261  m_delta[DeltaQOverP1] *= m_qOverPafterCalo / m_qOverPbeforeCalo;
262 
263  // TeV momentum scale to avoid magnitude mismatch (aka rounding)
264  m_derivQOverP0 = 1. / (m_delta[DeltaQOverP0] * Gaudi::Units::TeV);
265  m_derivQOverP1 = 1. / (m_delta[DeltaQOverP1] * Gaudi::Units::TeV);
266 
267  // and set qOverP
268  for (double& typ : m_qOverP) {
269  typ = m_parameters->qOverP();
270  }
271  m_qOverP[DeltaQOverP0] += m_delta[DeltaQOverP0];
272 
273  // use asymmetric error in case of significant energy deposit residual
274  if (m_parameters->fitEnergyDeposit() && m_asymmetricCaloEnergy) {
275  if (m_energyResidual * m_caloEnergyMeasurement->residual() < 0.) {
276  m_caloEnergyMeasurement->setSigma();
277  m_energyResidual = m_caloEnergyMeasurement->residual();
278  } else if (m_energyResidual < 1. &&
279  m_caloEnergyMeasurement->residual() > 1.) {
280  m_caloEnergyMeasurement->setSigmaPlus();
281  m_energyResidual = m_caloEnergyMeasurement->residual();
282  } else if (m_energyResidual > -1. &&
283  m_caloEnergyMeasurement->residual() < -1.) {
284  m_caloEnergyMeasurement->setSigmaMinus();
285  m_energyResidual = m_caloEnergyMeasurement->residual();
286  }
287  }
288 
289  // extrapolate to all measurement surfaces and compute intersects
290  return extrapolateToMeasurements(FittedTrajectory);
291 }
292 
293 void MeasurementProcessor::calculateResiduals(void) {
294  int nAlign = 0;
295  int nScat = 0;
296  for (auto* m : m_measurements) {
297  if (!(*m).numberDoF())
298  continue;
299 
300  // special measurements to constrain additional parameters
301  if (!m->isPositionMeasurement()) {
302  if (m->isScatterer()) // scattering centres
303  {
304  double phiResidual =
305  -m->weight() * m_parameters->scattererPhi(nScat) *
306  m->intersection(FittedTrajectory).direction().perp();
307  m->residual(phiResidual);
308  double thetaResidual =
309  -m->weight() * m_parameters->scattererTheta(nScat);
310  (*m).residual2(thetaResidual);
311  ++nScat;
312  } else if ((*m).isAlignment()) // alignment uncertainties
313  {
314  (*m).residual(-m->weight() *
315  (m_parameters->alignmentAngle(nAlign) -
316  m_parameters->alignmentAngleConstraint(nAlign)));
317  (*m).residual2(-m->weight2() *
318  (m_parameters->alignmentOffset(nAlign) -
319  m_parameters->alignmentOffsetConstraint(nAlign)));
320  ++nAlign;
321  } else if (m->isEnergyDeposit()) {
322  // Add the energy loss as a measurement
323  double E0 = 1. / std::abs(m_qOverPbeforeCalo);
324  double E1 = 1. / std::abs(m_qOverPafterCalo);
325  double residual = m->weight() * (E0 - E1 - m->energyLoss());
326  m->residual(residual);
327  }
328  continue;
329  }
330 
331  // position measurements
332  const TrackSurfaceIntersection& intersection =
333  m->intersection(FittedTrajectory);
334  const Amg::Vector3D& minimizationDirection = m->minimizationDirection();
335  Amg::Vector3D offset = m->position() - intersection.position();
336  if (m->isCluster()) // strip detector types
337  {
338  double residual = m->weight() * minimizationDirection.dot(offset);
339  if (m->alignmentParameter()) {
340  // propagate to residual using derivatives
341  unsigned param = m->alignmentParameter() - 1;
342  unsigned deriv = m_parameters->numberParameters() -
343  2 * m_parameters->numberAlignments() + 2 * param;
344  residual -=
345  m_parameters->alignmentAngle(param) * m->derivative(deriv) +
346  m_parameters->alignmentOffset(param) * m->derivative(deriv + 1);
347 
348  if (m->alignmentParameter2()) {
349  // propagate to residual using derivatives
350  param = m->alignmentParameter2() - 1;
351  deriv = m_parameters->numberParameters() -
352  2 * m_parameters->numberAlignments() + 2 * param;
353  residual -=
354  m_parameters->alignmentAngle(param) * m->derivative(deriv) +
355  m_parameters->alignmentOffset(param) * m->derivative(deriv + 1);
356  }
357  }
358  m->residual(residual);
359  if (!m->is2Dimensional())
360  continue;
361  double residual2 = m->weight2() * m->sensorDirection().dot(offset);
362  m->residual2(residual2);
363  } else if (m->isDrift() ||
364  m->isPerigee()) // else drift circles (perigee is similar)
365  {
366  double residual = m->weight() * (minimizationDirection.dot(offset) +
367  m->signedDriftDistance());
368  if (m->alignmentParameter()) {
369  // propagate to residual using derivatives
370  unsigned param = m->alignmentParameter() - 1;
371  int deriv = m_parameters->numberParameters() -
372  2 * m_parameters->numberAlignments() + 2 * param;
373  residual -=
374  m_parameters->alignmentAngle(param) * m->derivative(deriv) +
375  m_parameters->alignmentOffset(param) * m->derivative(deriv + 1);
376  if (deriv != m_parameters->firstAlignmentParameter())
377 
378  if (m->alignmentParameter2()) {
379  // propagate to residual using derivatives
380  param = m->alignmentParameter2() - 1;
381  deriv = m_parameters->numberParameters() -
382  2 * m_parameters->numberAlignments() + 2 * param;
383  residual -=
384  m_parameters->alignmentAngle(param) * m->derivative(deriv) +
385  m_parameters->alignmentOffset(param) * m->derivative(deriv + 1);
386  }
387  }
388  m->residual(residual);
389  // std::cout << " residual " << residual*m->sigma()
390  // << " drift distance " << m->signedDriftDistance() <<
391  // std::endl;
392  } else if (m->isPseudo()) // else pseudo measurement
393  {
394  double residual = m->weight() * minimizationDirection.dot(offset);
395  m->residual(residual);
396  if (!m->is2Dimensional())
397  continue;
398  double residual2 = m->weight2() * m->sensorDirection().dot(offset);
399  m->residual2(residual2);
400  } else if (m->isVertex()) {
401  double residual = m->weight() *
402  (minimizationDirection.x() * offset.x() +
403  minimizationDirection.y() * offset.y()) /
404  minimizationDirection.perp();
405  m->residual(residual);
406  if (!m->is2Dimensional())
407  continue;
408  double residual2 = m->weight2() * m->sensorDirection().dot(offset);
409  m->residual2(residual2);
410  }
411  }
412 }
413 
414 void MeasurementProcessor::fieldIntegralUncertainty(MsgStream& log,
415  Amg::MatrixX& covariance) {
416  // check field integral stability if there is a large error on the start
417  // position/direction but only meaningful when material taken into account
418  if (!m_parameters->numberScatterers() ||
419  !Amg::hasPositiveOrZeroDiagElems(covariance))
420  return;
421 
422  if (covariance(0, 0) + covariance(1, 1) < m_largeDeltaD0 * m_largeDeltaD0 &&
423  covariance(2, 2) < m_largeDeltaPhi0 * m_largeDeltaPhi0)
424  return;
425 
426  // FIXME: must also account for asymmetry in case of large momentum error
427  // when momentum correlation terms can be significantly overestimated
428 
429  // assume fieldIntegral proportional to angular deflection from vertex to last
430  // measurement as problems are always in the toroid take the theta deflection
431  // contribution to qOverP
432  double sinTheta =
433  1. / std::sqrt(1. + m_parameters->cotTheta() * m_parameters->cotTheta());
434  Amg::Vector3D startDirection(sinTheta * m_parameters->cosPhi(),
435  sinTheta * m_parameters->sinPhi(),
436  sinTheta * m_parameters->cotTheta());
437  const FitMeasurement& lastMeas = **m_measurements.rbegin();
438  Amg::Vector3D endDirection =
439  lastMeas.intersection(FittedTrajectory).direction();
440  if (startDirection.z() == 0. || endDirection.z() == 0.)
441  return;
442 
443  double deflectionPhi = startDirection.x() * endDirection.y() -
444  startDirection.y() * endDirection.x();
445  double deflectionTheta = startDirection.perp() * endDirection.z() -
446  startDirection.z() * endDirection.perp();
447 
448  // poorly measured phi
449  double shiftPhi0 = std::sqrt(covariance(2, 2));
450  if (shiftPhi0 > m_largeDeltaPhi0) {
451  Amg::Vector3D vertex(m_parameters->position());
452  double cosPhi = m_parameters->cosPhi() - m_parameters->sinPhi() * shiftPhi0;
453  double sinPhi = m_parameters->sinPhi() + m_parameters->cosPhi() * shiftPhi0;
454  double cotTheta = m_parameters->cotTheta();
455  sinTheta = 1. / std::sqrt(1. + cotTheta * cotTheta);
456  startDirection = Amg::Vector3D(sinTheta * cosPhi, sinTheta * sinPhi,
457  sinTheta * cotTheta);
458  m_vertexIntersect = TrackSurfaceIntersection(vertex, startDirection, 0.);
459  if (!extrapolateToMeasurements(DeltaPhi0))
460  return;
461 
462  endDirection = lastMeas.intersection(DeltaPhi0).direction();
463  double deltaPhi = startDirection.x() * endDirection.y() -
464  startDirection.y() * endDirection.x() - deflectionPhi;
465  double deltaTheta = startDirection.perp() * endDirection.z() -
466  startDirection.z() * endDirection.perp() -
467  deflectionTheta;
468  covariance(3, 3) += deltaTheta * deltaTheta;
469  if (std::abs(deflectionTheta) > 0.0001) {
470  double deltaQOverP =
471  (deltaTheta / deflectionTheta) * m_parameters->qOverP();
472  covariance(4, 4) += deltaQOverP * deltaQOverP;
473  }
474  if (log.level() < MSG::DEBUG) {
475  log << MSG::VERBOSE << std::setiosflags(std::ios::fixed)
476  << " fieldIntegralUncertainty: "
477  << " shiftPhi0 " << std::setw(8) << std::setprecision(3) << shiftPhi0
478  << " deflection (phi,theta) " << std::setw(8) << std::setprecision(4)
479  << deflectionPhi << std::setw(8) << std::setprecision(4)
480  << deflectionTheta << " diff " << std::setw(8)
481  << std::setprecision(4) << deltaPhi << std::setw(8)
482  << std::setprecision(4) << deltaTheta << endmsg;
483  return;
484  }
485  } else {
486  // compute average deflection for a one sigma move towards origin
487  double shiftD0 = std::sqrt(covariance(0, 0));
488  double shiftZ0 = std::sqrt(covariance(1, 1));
489  if (m_parameters->d0() > 0.)
490  shiftD0 = -shiftD0;
491  if (m_parameters->position().z() > 0.)
492  shiftZ0 = -shiftZ0;
493 
494  Amg::Vector3D vertex(-shiftD0 * m_parameters->sinPhi(),
495  shiftD0 * m_parameters->cosPhi(), shiftZ0);
496  vertex += m_parameters->position();
497  double dPhi = covariance(0, 2) / shiftD0;
498  double cosPhi = m_parameters->cosPhi() - m_parameters->sinPhi() * dPhi;
499  double sinPhi = m_parameters->sinPhi() + m_parameters->cosPhi() * dPhi;
500  double cotTheta = m_parameters->cotTheta();
501  sinTheta = 1. / std::sqrt(1. + cotTheta * cotTheta);
502  cotTheta -= covariance(1, 3) / (shiftZ0 * sinTheta * sinTheta);
503  startDirection = Amg::Vector3D(sinTheta * cosPhi, sinTheta * sinPhi,
504  sinTheta * cotTheta);
505  m_vertexIntersect = TrackSurfaceIntersection(vertex, startDirection, 0.);
506  if (!extrapolateToMeasurements(DeltaD0))
507  return;
508 
509  endDirection = lastMeas.intersection(DeltaD0).direction();
510  double deltaPhi = startDirection.x() * endDirection.y() -
511  startDirection.y() * endDirection.x() - deflectionPhi;
512  double deltaTheta = startDirection.perp() * endDirection.z() -
513  startDirection.z() * endDirection.perp() -
514  deflectionTheta;
515 
516  covariance(2, 2) += deltaPhi * deltaPhi;
517  covariance(3, 3) += deltaTheta * deltaTheta;
518  if (std::abs(deflectionTheta) > 0.0001) {
519  double deltaQOverP =
520  (deltaTheta / deflectionTheta) * m_parameters->qOverP();
521  covariance(4, 4) += deltaQOverP * deltaQOverP;
522  }
523 
524  if (log.level() < MSG::DEBUG) {
525  log << MSG::VERBOSE << std::setiosflags(std::ios::fixed)
526  << " fieldIntegralUncertainty: "
527  << " shiftD0 " << std::setw(8) << std::setprecision(1) << shiftD0
528  << " shiftZ0 " << std::setw(8) << std::setprecision(1) << shiftZ0
529  << " deflection (phi,theta) " << std::setw(8) << std::setprecision(4)
530  << deflectionPhi << std::setw(8) << std::setprecision(4)
531  << deflectionTheta << " diff " << std::setw(8)
532  << std::setprecision(4) << deltaPhi << std::setw(8)
533  << std::setprecision(4) << deltaTheta << endmsg;
534  }
535  }
536 }
537 
538 void MeasurementProcessor::propagationDerivatives(void) {
539  // compute additional derivatives when needed for covariance propagation.
540  // loop over measurements:
541  for (auto* m : m_measurements) {
542  // compute the D0 and Z0 derivs that don't already exist
543  if (!m->isPositionMeasurement() || m->numberDoF() > 1)
544  continue;
545  int derivativeFlag = 0;
546  if (m->numberDoF()) {
547  derivativeFlag = 0;
548  } else {
549  derivativeFlag = 2;
550  }
551 
552  if (m->isCluster()) {
553  clusterDerivatives(derivativeFlag, *m);
554  } else if (m->isDrift() || m->isPerigee() || m->isPseudo()) {
555  driftDerivatives(derivativeFlag, *m);
556  }
557  }
558 }
559 
560 void MeasurementProcessor::clusterDerivatives(int derivativeFlag,
561  FitMeasurement& measurement) {
562  // derivativeFlag definition: 0 take wrt Z0, 1 take wrt D0, 2 take wrt D0 and
563  // Z0
564  double weight = measurement.weight();
565  const TrackSurfaceIntersection& intersection =
566  measurement.intersection(FittedTrajectory);
567  const Amg::Vector3D& minimizationDirection =
568  measurement.minimizationDirection();
569  const Amg::Vector3D& sensorDirection = measurement.sensorDirection();
570 
571  // protect against grazing incidence
572  if (std::abs(measurement.normal().dot(intersection.direction())) < 0.001)
573  return;
574 
575  // transverse distance to measurement
576  double xDistance = intersection.position().x() - m_x0;
577  double yDistance = intersection.position().y() - m_y0;
578  double rDistance = -(m_cosPhi0 * xDistance + m_sinPhi0 * yDistance) /
579  (m_sinTheta0 * m_sinTheta0);
580 
581  if (derivativeFlag != 0) {
582  // momentum derivative - always numeric
583  if (m_parameters->fitEnergyDeposit() && measurement.afterCalo()) {
584  Amg::Vector3D offset = measurement.intersection(DeltaQOverP1).position() -
585  intersection.position();
586  measurement.derivative(
587  QOverP1, weight * minimizationDirection.dot(offset) * m_derivQOverP1);
588  } else if (m_parameters->fitMomentum()) {
589  Amg::Vector3D offset = measurement.intersection(DeltaQOverP0).position() -
590  intersection.position();
591  measurement.derivative(
592  QOverP0, weight * minimizationDirection.dot(offset) * m_derivQOverP0);
593  }
594 
595  // analytic derivatives in longitudinal direction
596  Amg::Vector3D weightedProjection =
597  weight * sensorDirection.cross(intersection.direction()) /
598  measurement.normal().dot(intersection.direction());
599  measurement.derivative(Z0, weightedProjection.z());
600  measurement.derivative(Theta0, weightedProjection.z() * rDistance);
601 
602  // numeric or analytic d0/phi derivative
603  if (measurement.numericalDerivative()) {
604  Amg::Vector3D offset = measurement.intersection(DeltaD0).position() -
605  intersection.position();
606  measurement.derivative(
607  D0, weight * minimizationDirection.dot(offset) / m_delta[DeltaD0]);
608  // offset = measurement.intersection(DeltaZ0).position() -
609  // intersection.position(); measurement.derivative(
610  // Z0,weight*minimizationDirection.dot(offset)/m_delta[DeltaZ0]);
611  offset = measurement.intersection(DeltaPhi0).position() -
612  intersection.position();
613  measurement.derivative(Phi0, weight * minimizationDirection.dot(offset) /
614  m_delta[DeltaPhi0]);
615  // offset = measurement.intersection(DeltaTheta0).position() -
616  // intersection.position(); measurement.derivative(Theta0,
617  // weight*minimizationDirection.dot(offset)/m_delta[DeltaTheta0]);
618  } else {
619  // transverse derivatives
620  measurement.derivative(D0, m_cosPhi0 * weightedProjection.y() -
621  m_sinPhi0 * weightedProjection.x());
622  measurement.derivative(Phi0, xDistance * weightedProjection.y() -
623  yDistance * weightedProjection.x());
624  }
625 
626  // derivatives wrt multiple scattering parameters
627  // similar to phi/theta
628  unsigned param = m_parameters->firstScatteringParameter();
629  std::vector<FitMeasurement*>::const_iterator s = m_scatterers.begin();
630  while (++param < measurement.lastParameter()) {
631  const TrackSurfaceIntersection& scatteringCentre =
632  (**s).intersection(FittedTrajectory);
633  double xDistScat =
634  intersection.position().x() - scatteringCentre.position().x();
635  double yDistScat =
636  intersection.position().y() - scatteringCentre.position().y();
637  double rDistScat = -(scatteringCentre.direction().x() * xDistScat +
638  scatteringCentre.direction().y() * yDistScat) /
639  (scatteringCentre.direction().perp2() *
640  scatteringCentre.direction().perp());
641  measurement.derivative(param, xDistScat * weightedProjection.y() -
642  yDistScat * weightedProjection.x());
643  measurement.derivative(++param, weightedProjection.z() * rDistScat);
644  ++s;
645  }
646 
647  if (measurement.alignmentParameter()) {
648  // TODO: fix projection factors wrt principal measurement axis from
649  // alignment surface
650  param = measurement.alignmentParameter();
651  FitMeasurement* fm = m_alignments[param - 1];
652  param = m_parameters->numberParameters() -
653  2 * m_parameters->numberAlignments() + 2 * (param - 1);
654  double distance = fm->surface()->normal().dot(
655  measurement.intersection(FittedTrajectory).position() -
656  fm->intersection(FittedTrajectory).position());
657  double projection = fm->surface()->normal().dot(
658  measurement.intersection(FittedTrajectory).direction());
659 
660  measurement.derivative(param, weight * distance);
661  measurement.derivative(++param, weight * projection);
662  if (measurement.alignmentParameter2()) {
663  param = measurement.alignmentParameter2();
664  fm = m_alignments[param - 1];
665  param = m_parameters->numberParameters() -
666  2 * m_parameters->numberAlignments() + 2 * (param - 1);
667  distance = fm->surface()->normal().dot(
668  measurement.intersection(FittedTrajectory).position() -
669  fm->intersection(FittedTrajectory).position());
670  projection = fm->surface()->normal().dot(
671  measurement.intersection(FittedTrajectory).direction());
672  measurement.derivative(param, weight * distance);
673  measurement.derivative(++param, weight * projection);
674  }
675  }
676  }
677 
678  // similar derivatives for the 2nd dimension,
679  // with minimization along sensor direction
680  if (derivativeFlag == 1)
681  return;
682  weight = measurement.weight2();
683 
684  // momentum derivative - always numeric
685  if (m_parameters->fitEnergyDeposit() && measurement.afterCalo()) {
686  Amg::Vector3D offset = measurement.intersection(DeltaQOverP1).position() -
687  intersection.position();
688  measurement.derivative2(
689  QOverP1, weight * sensorDirection.dot(offset) * m_derivQOverP1);
690  } else if (m_parameters->fitMomentum()) {
691  Amg::Vector3D offset = measurement.intersection(DeltaQOverP0).position() -
692  intersection.position();
693  measurement.derivative2(
694  QOverP0, weight * sensorDirection.dot(offset) * m_derivQOverP0);
695  }
696 
697  // analytic derivatives in longitudinal direction
698  // FIXME: orthogonal to 1st coord for DOF != 2
699  Amg::Vector3D weightedProjection =
700  -weight * minimizationDirection.cross(intersection.direction()) /
701  measurement.normal().dot(intersection.direction());
702  measurement.derivative2(Z0, weightedProjection.z());
703  measurement.derivative2(Theta0, weightedProjection.z() * rDistance);
704 
705  // numeric or analytic d0/phi derivative
706  if (measurement.numericalDerivative()) {
707  Amg::Vector3D offset =
708  measurement.intersection(DeltaD0).position() - intersection.position();
709  measurement.derivative2(
710  D0, weight * sensorDirection.dot(offset) / m_delta[DeltaD0]);
711  // offset = measurement.intersection(DeltaZ0).position() -
712  // intersection.position(); measurement.derivative2(
713  // Z0,weight*sensorDirection.dot(offset)/m_delta[DeltaZ0]);
714  offset = measurement.intersection(DeltaPhi0).position() -
715  intersection.position();
716  measurement.derivative2(
717  Phi0, weight * sensorDirection.dot(offset) / m_delta[DeltaPhi0]);
718  // offset = measurement.intersection(DeltaTheta0).position()
719  // - intersection.position(); measurement.derivative2(Theta0,
720  // weight*sensorDirection.dot(offset)/m_delta[DeltaTheta0]);
721  } else {
722  measurement.derivative2(D0, m_cosPhi0 * weightedProjection.y() -
723  m_sinPhi0 * weightedProjection.x());
724  measurement.derivative2(Phi0, xDistance * weightedProjection.y() -
725  yDistance * weightedProjection.x());
726  }
727 
728  // derivatives wrt multiple scattering parameters
729  // similar to phi/theta
730  unsigned param = m_parameters->firstScatteringParameter();
731  std::vector<FitMeasurement*>::const_iterator s = m_scatterers.begin();
732  while (++param < measurement.lastParameter()) {
733  const TrackSurfaceIntersection& scatteringCentre =
734  (**s).intersection(FittedTrajectory);
735  double xDistScat =
736  intersection.position().x() - scatteringCentre.position().x();
737  double yDistScat =
738  intersection.position().y() - scatteringCentre.position().y();
739  double rDistScat = -(scatteringCentre.direction().x() * xDistScat +
740  scatteringCentre.direction().y() * yDistScat) /
741  (scatteringCentre.direction().perp2() *
742  scatteringCentre.direction().perp());
743  measurement.derivative2(param, xDistScat * weightedProjection.y() -
744  yDistScat * weightedProjection.x());
745  measurement.derivative2(++param, weightedProjection.z() * rDistScat);
746  ++s;
747  }
748 }
749 
750 void MeasurementProcessor::driftDerivatives(int derivativeFlag,
751  FitMeasurement& measurement) {
752  // transverse distance to measurement
753  const TrackSurfaceIntersection& intersection =
754  measurement.intersection(FittedTrajectory);
755  double xDistance = intersection.position().x() - m_x0;
756  double yDistance = intersection.position().y() - m_y0;
757  double rDistance = -(m_cosPhi0 * xDistance + m_sinPhi0 * yDistance) /
758  (m_sinTheta0 * m_sinTheta0);
759 
760  // derivativeFlag definition: 0 take wrt Z0, 1 take wrt D0, 2 take wrt D0 and
761  // Z0
762  if (derivativeFlag != 0) {
763  double weight = measurement.weight();
764  const Amg::Vector3D& driftDirection = measurement.minimizationDirection();
765  double xDrift = driftDirection.x();
766  double yDrift = driftDirection.y();
767 
768  // momentum derivative - always numeric
769  if (m_parameters->fitEnergyDeposit() && measurement.afterCalo()) {
770  Amg::Vector3D offset = measurement.intersection(DeltaQOverP1).position() -
771  intersection.position();
772  measurement.derivative(
773  QOverP1, weight * driftDirection.dot(offset) * m_derivQOverP1);
774  } else if (m_parameters->fitMomentum()) {
775  Amg::Vector3D offset = measurement.intersection(DeltaQOverP0).position() -
776  intersection.position();
777  measurement.derivative(
778  QOverP0, weight * driftDirection.dot(offset) * m_derivQOverP0);
779  }
780 
781  if (measurement.numericalDerivative()) {
782  Amg::Vector3D offset = measurement.intersection(DeltaD0).position() -
783  intersection.position();
784  measurement.derivative(
785  D0, weight * driftDirection.dot(offset) / m_delta[DeltaD0]);
786  // offset =
787  // measurement.intersection(DeltaZ0).position() - intersection.position();
788  // measurement.derivative(
789  // Z0,weight*driftDirection.dot(offset)/m_delta[DeltaZ0]);
790  offset = measurement.intersection(DeltaPhi0).position() -
791  intersection.position();
792  measurement.derivative(
793  Phi0, weight * driftDirection.dot(offset) / m_delta[DeltaPhi0]);
794  // offset =
795  // measurement.intersection(DeltaTheta0).position() -
796  // intersection.position();
797  // measurement.derivative(Theta0,weight*driftDirection.dot(offset)/m_delta[DeltaTheta0]);
798  } else if (
799  m_phiInstability
800  // &&
801  // std::abs(intersection.direction().z()*m_vertexIntersect->direction().perp()
802  // -
803  // m_vertexIntersect->direction().z()*intersection.direction().perp())
804  // > m_toroidTurn
805  && !measurement.isPseudo()) {
806  measurement.derivative(D0, 0.);
807  measurement.derivative(Phi0, 0.);
808  } else {
809  measurement.derivative(
810  D0, weight * (m_cosPhi0 * yDrift - m_sinPhi0 * xDrift));
811  measurement.derivative(
812  Phi0, weight * (xDistance * yDrift - yDistance * xDrift));
813  }
814 
815  measurement.derivative(Z0, weight * driftDirection.z());
816  measurement.derivative(Theta0, weight * driftDirection.z() * rDistance);
817 
818  // derivatives wrt multiple scattering parameters
819  // similar to phi/theta
820  unsigned param = m_parameters->firstScatteringParameter();
821  std::vector<FitMeasurement*>::const_iterator s = m_scatterers.begin();
822  while (++param < measurement.lastParameter()) {
823  const TrackSurfaceIntersection& scatteringCentre =
824  (**s).intersection(FittedTrajectory);
825  double xDistScat =
826  intersection.position().x() - scatteringCentre.position().x();
827  double yDistScat =
828  intersection.position().y() - scatteringCentre.position().y();
829  double rDistScat = -(scatteringCentre.direction().x() * xDistScat +
830  scatteringCentre.direction().y() * yDistScat) /
831  (scatteringCentre.direction().perp2() *
832  scatteringCentre.direction().perp());
833  measurement.derivative(
834  param, weight * (xDistScat * yDrift - yDistScat * xDrift));
835  measurement.derivative(++param, weight * driftDirection.z() * rDistScat);
836  ++s;
837  }
838 
839  // derivatives wrt alignment parameters
840  if (measurement.alignmentParameter()) {
841  param = measurement.alignmentParameter();
842  FitMeasurement* fm = m_alignments[param - 1];
843  param = m_parameters->numberParameters() -
844  2 * m_parameters->numberAlignments() + 2 * (param - 1);
845 
846  // angle derivative (delta_theta) similar to scatterer delta_theta
847  const TrackSurfaceIntersection& alignmentCentre =
848  fm->intersection(FittedTrajectory);
849  double xDistance =
850  intersection.position().x() - alignmentCentre.position().x();
851  double yDistance =
852  intersection.position().y() - alignmentCentre.position().y();
853  double rDistance = -(alignmentCentre.direction().x() * xDistance +
854  alignmentCentre.direction().y() * yDistance) /
855  (alignmentCentre.direction().perp2() *
856  alignmentCentre.direction().perp());
857  measurement.derivative(param, weight * driftDirection.z() * rDistance);
858 
859  // offset derivative: barrel (endcap) projection factor onto the surface
860  // plane in the z (r) direction
861  double projection = 0;
862  const Surface& surface = *fm->surface();
863  if (std::abs(surface.normal().z()) > 0.5) {
864  projection = (driftDirection.x() * surface.center().x() +
865  driftDirection.y() * surface.center().y()) /
866  surface.center().perp();
867  } else {
868  projection = driftDirection.z();
869  }
870  measurement.derivative(++param, weight * projection);
871 
872  if (measurement.alignmentParameter2()) {
873  param = measurement.alignmentParameter2();
874  fm = m_alignments[param - 1];
875  param = m_parameters->numberParameters() -
876  2 * m_parameters->numberAlignments() + 2 * (param - 1);
877 
878  // angle derivative (delta_theta) similar to scatterer delta_theta
879  const TrackSurfaceIntersection& alignmentCentre =
880  fm->intersection(FittedTrajectory);
881  xDistance =
882  intersection.position().x() - alignmentCentre.position().x();
883  yDistance =
884  intersection.position().y() - alignmentCentre.position().y();
885  rDistance = -(alignmentCentre.direction().x() * xDistance +
886  alignmentCentre.direction().y() * yDistance) /
887  (alignmentCentre.direction().perp2() *
888  alignmentCentre.direction().perp());
889  measurement.derivative(param, weight * driftDirection.z() * rDistance);
890 
891  // offset derivative: barrel (endcap) projection factor onto the surface
892  // plane in the z (r) direction
893  const Surface& surface = *fm->surface();
894  if (surface.normal().dot(surface.center().unit()) < 0.5) {
895  projection = driftDirection.z();
896  } else {
897  projection = (driftDirection.x() * surface.center().x() +
898  driftDirection.y() * surface.center().y()) /
899  surface.center().perp();
900  }
901  measurement.derivative(++param, weight * projection);
902  }
903  }
904  }
905 
906  // similar for derivatives along the wire direction
907  if (derivativeFlag == 1)
908  return;
909  double weight = measurement.weight2();
910  const Amg::Vector3D& wireDirection = measurement.sensorDirection();
911  double xWire = wireDirection.x();
912  double yWire = wireDirection.y();
913 
914  // momentum derivative - always numeric
915  if (m_parameters->fitEnergyDeposit() && measurement.afterCalo()) {
916  Amg::Vector3D offset = measurement.intersection(DeltaQOverP1).position() -
917  intersection.position();
918  measurement.derivative2(
919  QOverP1, weight * wireDirection.dot(offset) * m_derivQOverP1);
920  } else if (m_parameters->fitMomentum()) {
921  Amg::Vector3D offset = measurement.intersection(DeltaQOverP0).position() -
922  intersection.position();
923  measurement.derivative2(
924  QOverP0, weight * wireDirection.dot(offset) * m_derivQOverP0);
925  }
926 
927  if (measurement.numericalDerivative()) {
928  Amg::Vector3D offset =
929  measurement.intersection(DeltaD0).position() - intersection.position();
930  measurement.derivative2(
931  D0, weight * wireDirection.dot(offset) / m_delta[DeltaD0]);
932  // offset = measurement.intersection(DeltaZ0).position() -
933  // intersection.position();
934  // measurement.derivative(Z0,weight*wireDirection.dot(offset)/m_delta[DeltaZ0]);
935  offset = measurement.intersection(DeltaPhi0).position() -
936  intersection.position();
937  measurement.derivative2(
938  Phi0, weight * wireDirection.dot(offset) / m_delta[DeltaPhi0]);
939  // offset = measurement.intersection(DeltaTheta0).position() -
940  // intersection.position();
941  // measurement.derivative(Theta0,weight*wireDirection.dot(offset)/m_delta[DeltaTheta0]);
942  } else {
943  measurement.derivative2(D0,
944  weight * (m_cosPhi0 * yWire - m_sinPhi0 * xWire));
945  measurement.derivative2(Phi0,
946  weight * (xDistance * yWire - yDistance * xWire));
947  }
948 
949  measurement.derivative2(Z0, weight * wireDirection.z());
950  measurement.derivative2(Theta0, weight * wireDirection.z() * rDistance);
951 
952  // derivatives wrt multiple scattering parameters
953  // similar to phi/theta
954  unsigned param = m_parameters->firstScatteringParameter();
955  std::vector<FitMeasurement*>::const_iterator s = m_scatterers.begin();
956  while (++param < measurement.lastParameter()) {
957  const TrackSurfaceIntersection& scatteringCentre =
958  (**s).intersection(FittedTrajectory);
959  double xDistScat =
960  intersection.position().x() - scatteringCentre.position().x();
961  double yDistScat =
962  intersection.position().y() - scatteringCentre.position().y();
963  double rDistScat = -(scatteringCentre.direction().x() * xDistScat +
964  scatteringCentre.direction().y() * yDistScat) /
965  (scatteringCentre.direction().perp2() *
966  scatteringCentre.direction().perp());
967  measurement.derivative2(param,
968  weight * (xDistScat * yWire - yDistScat * xWire));
969  measurement.derivative2(++param, weight * wireDirection.z() * rDistScat);
970  ++s;
971  }
972 }
973 
974 bool MeasurementProcessor::extrapolateToMeasurements(ExtrapolationType type) {
975  // extrapolate to all measurement surfaces
976  int nScat = 0;
977  double qOverP = m_qOverP[type];
978  std::optional<TrackSurfaceIntersection> intersection = m_vertexIntersect;
979  const Surface* surface = nullptr;
980  const EventContext& ctx = Gaudi::Hive::currentContext();
981  // careful: use RungeKutta for extrapolation to vertex measurement
982  std::vector<FitMeasurement*>::iterator m = m_measurements.begin();
983  if ((**m).isVertex()) {
984  if (m_useStepPropagator == 99) {
985  std::optional<TrackSurfaceIntersection> newIntersectionSTEP =
986  m_stepPropagator->intersectSurface(ctx, *(**m).surface(),
987  *intersection, qOverP, m_stepField,
988  Trk::muon);
989  double exdist = 0;
990  if (newIntersectionSTEP)
991  exdist =
992  (newIntersectionSTEP->position() - intersection->position()).mag();
993  intersection = m_rungeKuttaIntersector->intersectSurface(
994  *(**m).surface(), *intersection, qOverP);
995  if (newIntersectionSTEP) {
996  double dist =
997  1000. *
998  (newIntersectionSTEP->position() - intersection->position()).mag();
999  std::cout << " iMeasProc 1 distance STEP and Intersector " << dist
1000  << " ex dist " << exdist << std::endl;
1001  if (dist > 10.)
1002  std::cout << " iMeasProc 1 ALARM distance STEP and Intersector "
1003  << dist << " ex dist " << exdist << std::endl;
1004  } else {
1005  if (intersection)
1006  std::cout << " iMeasProc 1 ALARM STEP did not intersect! "
1007  << std::endl;
1008  }
1009  } else {
1010  intersection = m_useStepPropagator == 0
1011  ? m_rungeKuttaIntersector->intersectSurface(
1012  *(**m).surface(), *intersection, qOverP)
1013  : m_stepPropagator->intersectSurface(
1014  ctx, *(**m).surface(), *intersection, qOverP,
1015  m_stepField, Trk::muon);
1016  }
1017  surface = (**m).surface();
1018  if (!intersection) {
1019  return false;
1020  }
1021  (**m).intersection(type, intersection);
1022  if (type == FittedTrajectory)
1023  (**m).qOverP(qOverP);
1024  ++m;
1025  }
1026 
1027  for (; m != m_measurements.end(); ++m) {
1028  if ((**m).afterCalo()) {
1029  if (type == DeltaQOverP0)
1030  continue;
1031  } else if (type == DeltaQOverP1) {
1032  intersection = (**m).intersection(FittedTrajectory);
1033  qOverP = (**m).qOverP();
1034  if ((**m).numberDoF() && (**m).isScatterer())
1035  ++nScat;
1036  continue;
1037  }
1038 
1039  // to avoid rounding: copy intersection if at previous surface
1040  if ((**m).surface() != surface) {
1041  if (m_useStepPropagator == 99) {
1042  std::optional<TrackSurfaceIntersection> newIntersectionSTEP =
1043  m_stepPropagator->intersectSurface(ctx, *(**m).surface(),
1044  *intersection, qOverP,
1045  m_stepField, Trk::muon);
1046  double exdist = 0;
1047  if (newIntersectionSTEP)
1048  exdist = (newIntersectionSTEP->position() - intersection->position())
1049  .mag();
1050 
1051  intersection = m_intersector->intersectSurface(*(**m).surface(),
1052  *intersection, qOverP);
1053  if (newIntersectionSTEP) {
1054  double dist = 1000. * (newIntersectionSTEP->position() -
1055  intersection->position())
1056  .mag();
1057  std::cout << " iMeasProc 2 distance STEP and Intersector " << dist
1058  << " ex dist " << exdist << std::endl;
1059  if (dist > 10.)
1060  std::cout << " iMeasProc 2 ALARM distance STEP and Intersector "
1061  << dist << " ex dist " << exdist << std::endl;
1062  } else {
1063  if (intersection)
1064  std::cout << " iMeasProc 2 ALARM STEP did not intersect! "
1065  << std::endl;
1066  }
1067  } else {
1068  intersection = m_useStepPropagator == 0
1069  ? m_intersector->intersectSurface(
1070  *(**m).surface(), *intersection, qOverP)
1071  : m_stepPropagator->intersectSurface(
1072  ctx, *(**m).surface(), *intersection, qOverP,
1073  m_stepField, Trk::muon);
1074  }
1075  surface = (**m).surface();
1076  } else {
1077  intersection = TrackSurfaceIntersection(*intersection);
1078  }
1079  if (!intersection) {
1080  return false;
1081  }
1082 
1083  // alignment and material effects (energy loss and trajectory deviations)
1084  if ((**m).isEnergyDeposit() || (**m).isScatterer()) {
1085  // apply fitted parameters
1086  if ((**m).numberDoF()) {
1087  if ((**m).isEnergyDeposit()) // reserved for calorimeter
1088  {
1089  if (type == FittedTrajectory) {
1090  m_qOverPbeforeCalo = qOverP;
1091  m_qOverPafterCalo = m_parameters->qOverP1();
1092  qOverP = m_qOverPafterCalo;
1093  (**m).qOverP(qOverP);
1094  } else if (type == DeltaQOverP1) {
1095  intersection =
1096  TrackSurfaceIntersection((**m).intersection(FittedTrajectory));
1097  qOverP = m_qOverPafterCalo + m_delta[DeltaQOverP1];
1098  } else {
1099  qOverP = m_qOverPafterCalo;
1100  }
1101  (**m).intersection(type, intersection);
1102  continue;
1103  }
1104 
1105  if ((**m).isScatterer()) // scattering centre
1106  {
1107  // update track direction for fitted scattering
1108  double sinDeltaPhi = std::sin(m_parameters->scattererPhi(nScat));
1109  double cosDeltaPhi = std::sqrt(1. - sinDeltaPhi * sinDeltaPhi);
1110  double sinDeltaTheta = std::sin(m_parameters->scattererTheta(nScat));
1111  double tanDeltaTheta = sinDeltaTheta;
1112  if (std::abs(sinDeltaTheta) < 1.)
1113  tanDeltaTheta /= std::sqrt(1. - sinDeltaTheta * sinDeltaTheta);
1114  Amg::Vector3D trackDirection(intersection->direction());
1115  trackDirection /= trackDirection.perp();
1116  double cosPhi = trackDirection.x() * cosDeltaPhi -
1117  trackDirection.y() * sinDeltaPhi;
1118  double sinPhi = trackDirection.y() * cosDeltaPhi +
1119  trackDirection.x() * sinDeltaPhi;
1120  double cotTheta = (trackDirection.z() - tanDeltaTheta) /
1121  (1. + trackDirection.z() * tanDeltaTheta);
1122  trackDirection = Amg::Vector3D(cosPhi, sinPhi, cotTheta);
1123  trackDirection = trackDirection.unit();
1124  intersection =
1125  TrackSurfaceIntersection(intersection->position(), trackDirection,
1126  intersection->pathlength());
1127  ++nScat;
1128  }
1129  }
1130  }
1131 
1132  // apply unfitted energy loss
1133  if (m_parameters->fitMomentum()) {
1134  // at extreme momenta use series expansion to avoid rounding (or FPE)
1135  if (m_parameters->extremeMomentum()) {
1136  double loss = (**m).energyLoss() * std::abs(qOverP);
1137  qOverP *= 1. + loss + loss * loss + loss * loss * loss;
1138  } else {
1139  double momentum = 1. / qOverP;
1140  if (momentum > 0.) {
1141  momentum -= (**m).energyLoss();
1142  } else {
1143  momentum += (**m).energyLoss();
1144  }
1145  qOverP = 1. / momentum;
1146  }
1147  }
1148 
1149  // store intersection and momentum vector (leaving surface)
1150  if (type == FittedTrajectory)
1151  (**m).qOverP(qOverP);
1152  (**m).intersection(type, intersection);
1153  }
1154 
1155  return true;
1156 }
1157 
1158 } // namespace Trk
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Trk::y
@ y
Definition: ParamDefs.h:56
Trk::FitMeasurement::intersection
const TrackSurfaceIntersection & intersection(ExtrapolationType type) const
Definition: FitMeasurement.h:359
Trk::FitParameters::alignmentAngleConstraint
double alignmentAngleConstraint(int alignment) const
Definition: FitParameters.h:160
Trk::QOverP1
@ QOverP1
Definition: ParameterType.h:18
Trk::MeasurementProcessor::MeasurementProcessor
MeasurementProcessor(bool asymmetricCaloEnergy, Amg::MatrixX &derivativeMatrix, const ToolHandle< IIntersector > &intersector, std::vector< FitMeasurement * > &measurements, FitParameters *parameters, const ToolHandle< IIntersector > &rungeKuttaIntersector, const ToolHandle< IPropagator > &stepPropagator, int useStepPropagator)
Definition: MeasurementProcessor.cxx:32
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
Trk::MeasurementProcessor::calculateFittedTrajectory
bool calculateFittedTrajectory(int iteration)
Definition: MeasurementProcessor.cxx:240
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:57
Trk::DeltaZ0
@ DeltaZ0
Definition: ExtrapolationType.h:21
Trk::MagneticFieldProperties
Definition: MagneticFieldProperties.h:31
Trk::MeasurementProcessor::m_qOverPbeforeCalo
double m_qOverPbeforeCalo
Definition: MeasurementProcessor.h:81
Surface.h
Trk::MeasurementProcessor::extrapolateToMeasurements
bool extrapolateToMeasurements(ExtrapolationType type)
Definition: MeasurementProcessor.cxx:977
Trk::FitParameters::alignmentAngle
double alignmentAngle(int alignment) const
Definition: FitParameters.h:156
Trk::MeasurementProcessor::m_derivQOverP1
double m_derivQOverP1
Definition: MeasurementProcessor.h:68
Trk::FitParameters::fitMomentum
bool fitMomentum(void) const
Definition: FitParameters.h:225
Trk::FitParameters::firstScatteringParameter
int firstScatteringParameter(void) const
Definition: FitParameters.h:217
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
Trk::MeasurementProcessor::m_intersectStartingValue
TrackSurfaceIntersection m_intersectStartingValue
Definition: MeasurementProcessor.h:91
Trk::FitMeasurement::setSigmaMinus
void setSigmaMinus(void)
Definition: FitMeasurement.h:557
xAOD::deltaPhi
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
Definition: L2StandAloneMuon_v1.cxx:160
Trk::QOverP0
@ QOverP0
Definition: ParameterType.h:18
Trk::FitParameters::fitEnergyDeposit
bool fitEnergyDeposit(void) const
Definition: FitParameters.h:221
Trk::MeasurementProcessor::m_intersector
const ToolHandle< IIntersector > & m_intersector
Definition: MeasurementProcessor.h:72
Trk::FitMeasurement::surface
const Surface * surface(void) const
Definition: FitMeasurement.h:585
Trk::MeasurementProcessor::m_energyResidual
double m_energyResidual
Definition: MeasurementProcessor.h:69
Trk::MeasurementProcessor::m_stepField
Trk::MagneticFieldProperties m_stepField
Definition: MeasurementProcessor.h:96
Trk::MeasurementProcessor::m_useStepPropagator
int m_useStepPropagator
Definition: MeasurementProcessor.h:85
Trk::FitParameters::alignmentOffsetConstraint
double alignmentOffsetConstraint(int alignment) const
Definition: FitParameters.h:168
Trk::MeasurementProcessor::m_delta
double m_delta[ExtrapolationTypes]
Definition: MeasurementProcessor.h:66
Trk::FitParameters::sinTheta
double sinTheta(void) const
Definition: FitParameters.h:277
Trk::MeasurementProcessor::m_asymmetricCaloEnergy
bool m_asymmetricCaloEnergy
Definition: MeasurementProcessor.h:62
Trk::MeasurementProcessor::m_vertexIntersect
TrackSurfaceIntersection m_vertexIntersect
Definition: MeasurementProcessor.h:90
Trk::MeasurementProcessor::m_y0
double m_y0
Definition: MeasurementProcessor.h:93
Trk::MeasurementProcessor::m_measurements
std::vector< FitMeasurement * > & m_measurements
Definition: MeasurementProcessor.h:75
IPropagator.h
python.SystemOfUnits.TeV
int TeV
Definition: SystemOfUnits.py:158
Trk::TrackSurfaceIntersection
Definition: TrackSurfaceIntersection.h:32
Trk::MeasurementProcessor::m_z0
double m_z0
Definition: MeasurementProcessor.h:94
Trk::FitParameters::scattererTheta
double scattererTheta(int scatterer) const
Definition: FitParameters.h:269
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Trk::MeasurementProcessor::m_sinTheta0
double m_sinTheta0
Definition: MeasurementProcessor.h:88
Trk::Z0
@ Z0
Definition: ParameterType.h:18
Trk::FitMeasurement::setSigmaSymmetric
void setSigmaSymmetric(void)
Definition: FitMeasurement.cxx:1025
Trk::MeasurementProcessor::m_largeDeltaD0
double m_largeDeltaD0
Definition: MeasurementProcessor.h:73
intersection
std::vector< std::string > intersection(std::vector< std::string > &v1, std::vector< std::string > &v2)
Definition: compareFlatTrees.cxx:25
Trk::FitParameters::ptInv0
double ptInv0(void) const
Definition: FitParameters.h:249
FitParameters.h
Trk::MeasurementProcessor::m_rungeKuttaIntersector
const ToolHandle< IIntersector > & m_rungeKuttaIntersector
Definition: MeasurementProcessor.h:83
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
Trk::FitParameters::alignmentOffset
double alignmentOffset(int alignment) const
Definition: FitParameters.h:164
Trk::FastField
@ FastField
call the fast field access method of the FieldSvc
Definition: MagneticFieldMode.h:20
Trk::FitParameters::qOverP1
double qOverP1(void) const
Definition: FitParameters.h:261
Trk::FitParameters::numberScatterers
int numberScatterers(void) const
Definition: FitParameters.h:245
Trk::MeasurementProcessor::fieldIntegralUncertainty
void fieldIntegralUncertainty(MsgStream &log, Amg::MatrixX &covariance)
Definition: MeasurementProcessor.cxx:417
Trk::FitParameters::extremeMomentum
bool extremeMomentum(void) const
Definition: FitParameters.h:205
Trk::MeasurementProcessor::clusterDerivatives
void clusterDerivatives(int derivativeFlag, FitMeasurement &measurement)
Definition: MeasurementProcessor.cxx:563
Trk::FitParameters::numberParameters
int numberParameters(void) const
Definition: FitParameters.h:241
Trk::MeasurementProcessor::m_phiInstability
bool m_phiInstability
Definition: MeasurementProcessor.h:79
Trk::MeasurementProcessor::driftDerivatives
void driftDerivatives(int derivativeFlag, FitMeasurement &measurement)
Definition: MeasurementProcessor.cxx:753
Trk::FitParameters::intersection
TrackSurfaceIntersection intersection(void) const
Definition: FitParameters.cxx:154
Trk::MeasurementProcessor::m_scatterers
std::vector< FitMeasurement * > m_scatterers
Definition: MeasurementProcessor.h:86
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
Trk::DeltaTheta0
@ DeltaTheta0
Definition: ExtrapolationType.h:23
Trk::TrackSurfaceIntersection::position
const Amg::Vector3D & position() const
Method to retrieve the position of the Intersection.
Definition: TrackSurfaceIntersection.h:80
Trk::FitMeasurement
Definition: FitMeasurement.h:40
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:538
Trk::MeasurementProcessor::m_qOverPafterCalo
double m_qOverPafterCalo
Definition: MeasurementProcessor.h:82
Trk::DeltaD0
@ DeltaD0
Definition: ExtrapolationType.h:20
Trk::Surface::normal
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
Trk::FitParameters::position
const Amg::Vector3D & position(void) const
Definition: FitParameters.h:253
Trk::MeasurementProcessor::calculateResiduals
void calculateResiduals(void)
Definition: MeasurementProcessor.cxx:296
Trk::MeasurementProcessor::m_derivQOverP0
double m_derivQOverP0
Definition: MeasurementProcessor.h:67
Trk::FitParameters::qOverP
double qOverP(void) const
Definition: FitParameters.h:257
Trk::MeasurementProcessor::m_parameters
FitParameters * m_parameters
Definition: MeasurementProcessor.h:78
Trk::MeasurementProcessor::m_qOverP
double m_qOverP[ExtrapolationTypes]
Definition: MeasurementProcessor.h:80
Trk::muon
@ muon
Definition: ParticleHypothesis.h:28
Trk::Theta0
@ Theta0
Definition: ParameterType.h:18
Trk::FitParameters::cosPhi
double cosPhi(void) const
Definition: FitParameters.h:172
Trk::DeltaQOverP1
@ DeltaQOverP1
Definition: ExtrapolationType.h:25
Trk::MeasurementProcessor::m_sinPhi0
double m_sinPhi0
Definition: MeasurementProcessor.h:87
Trk::FitParameters::cosTheta
double cosTheta(void) const
Definition: FitParameters.h:176
Trk::FullField
@ FullField
Field is set to be realistic, but within a given Volume.
Definition: MagneticFieldMode.h:21
Trk::FitParameters::firstAlignmentParameter
int firstAlignmentParameter(void) const
Definition: FitParameters.h:213
Athena::Units
Definition: Units.h:45
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::ExtrapolationType
ExtrapolationType
Definition: ExtrapolationType.h:18
Trk::MeasurementProcessor::calculateDerivatives
bool calculateDerivatives(void)
Definition: MeasurementProcessor.cxx:123
Trk::MeasurementProcessor::m_largeDeltaPhi0
double m_largeDeltaPhi0
Definition: MeasurementProcessor.h:74
Trk::FittedTrajectory
@ FittedTrajectory
Definition: ExtrapolationType.h:19
Trk::MeasurementProcessor::m_alignments
std::vector< FitMeasurement * > m_alignments
Definition: MeasurementProcessor.h:61
Trk::MeasurementProcessor::~MeasurementProcessor
~MeasurementProcessor(void)
Trk::FitParameters::d0
double d0(void) const
Definition: FitParameters.h:184
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
IIntersector.h
Trk::FitParameters::cotTheta
double cotTheta(void) const
Definition: FitParameters.h:180
Trk::TrackSurfaceIntersection::direction
const Amg::Vector3D & direction() const
Method to retrieve the direction at the Intersection.
Definition: TrackSurfaceIntersection.h:88
Trk::D0
@ D0
Definition: ParameterType.h:18
FitMeasurement.h
python.SystemOfUnits.mm
int mm
Definition: SystemOfUnits.py:83
EventPrimitivesCovarianceHelpers.h
Trk::DeltaPhi0
@ DeltaPhi0
Definition: ExtrapolationType.h:22
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
Trk::FitMeasurement::residual
double residual(void) const
Definition: FitMeasurement.h:494
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
Trk::MeasurementProcessor::m_numericDerivatives
bool m_numericDerivatives
Definition: MeasurementProcessor.h:77
MeasurementProcessor.h
Trk::DeltaQOverP0
@ DeltaQOverP0
Definition: ExtrapolationType.h:24
Trk::ExtrapolationTypes
@ ExtrapolationTypes
Definition: ExtrapolationType.h:26
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
Trk::MeasurementProcessor::m_cosPhi0
double m_cosPhi0
Definition: MeasurementProcessor.h:64
convertTimingResiduals.offset
offset
Definition: convertTimingResiduals.py:71
TRT::Track::cotTheta
@ cotTheta
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:65
Trk::FitParameters::sinPhi
double sinPhi(void) const
Definition: FitParameters.h:273
Gaudi
=============================================================================
Definition: CaloGPUClusterAndCellDataMonitorOptions.h:273
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Trk::FitParameters::numberAlignments
int numberAlignments(void) const
Definition: FitParameters.h:233
Trk::MeasurementProcessor::m_cosTheta0
double m_cosTheta0
Definition: MeasurementProcessor.h:65
Trk::Phi0
@ Phi0
Definition: ParameterType.h:18
Trk::FitParameters::scattererPhi
double scattererPhi(int scatterer) const
Definition: FitParameters.h:265
Trk::FitParameters
Definition: FitParameters.h:29
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
Trk::FitMeasurement::setSigma
void setSigma(void)
Definition: FitMeasurement.h:553
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Trk::x
@ x
Definition: ParamDefs.h:55
Amg::hasPositiveOrZeroDiagElems
bool hasPositiveOrZeroDiagElems(const AmgSymMatrix(N) &mat)
Returns true if all diagonal elements of the covariance matrix are finite aka sane in the above defin...
Definition: EventPrimitivesCovarianceHelpers.h:73
TrackSurfaceIntersection.h
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
Trk::MeasurementProcessor::propagationDerivatives
void propagationDerivatives(void)
Definition: MeasurementProcessor.cxx:541
Trk::MeasurementProcessor::m_caloEnergyMeasurement
FitMeasurement * m_caloEnergyMeasurement
Definition: MeasurementProcessor.h:63
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
Trk::MeasurementProcessor::m_x0
double m_x0
Definition: MeasurementProcessor.h:92
Trk::FitMeasurement::setSigmaPlus
void setSigmaPlus(void)
Definition: FitMeasurement.h:561
Trk::MeasurementProcessor::m_stepPropagator
const ToolHandle< IPropagator > & m_stepPropagator
Definition: MeasurementProcessor.h:84
ParameterType.h
Generated on Thu Nov 7 2024 21:20:50 for ATLAS Offline Software by   doxygen 1.8.18