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