ATLAS Offline Software
Loading...
Searching...
No Matches
RiddersAlgorithm.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4
6// RiddersAlgorithm.cxx, (c) ATLAS Detector software
8
9#include "TrkExAlgs/RiddersAlgorithm.h"
10// Trk stuff
11#include "TrkEventPrimitives/TransportJacobian.h"
12#include "TrkGeometry/MagneticFieldProperties.h"
13#include "TrkSurfaces/PlaneSurface.h"
14#include "TrkSurfaces/CylinderSurface.h"
15#include "TrkParameters/TrackParameters.h"
16#include "TrkVolumes/CylinderVolumeBounds.h"
17// Validation mode - TTree includes
18#include "TTree.h"
19#include "GaudiKernel/ITHistSvc.h"
20#include <cmath>
21
22//================ Constructor =================================================
23
24Trk::RiddersAlgorithm::RiddersAlgorithm(const std::string& name, ISvcLocator* pSvcLocator) :
25 AthAlgorithm(name, pSvcLocator) {}
26
27//================ Destructor =================================================
28
34
35//================ Initialisation =================================================
36
37StatusCode Trk::RiddersAlgorithm::initialize() {
38 // Code entered here will be executed once at program start.
39 ATH_MSG_INFO(" initialize()");
40
41 ATH_CHECK(m_propagator.retrieve());
42
43 // Prepare the magnetic field properties
44 if (!m_useCustomField) m_magFieldProperties = new Trk::MagneticFieldProperties();
45 else {
46 // the field
47 Amg::Vector3D bField(0., 0., m_fieldValue);
48 // create the custom magnetic field
49 m_magFieldProperties = new Trk::MagneticFieldProperties(bField);
50 }
51
52 // intialize the random number generators
53 m_gaussDist = new Rndm::Numbers(randSvc(), Rndm::Gauss(0., 1.));
54 m_flatDist = new Rndm::Numbers(randSvc(), Rndm::Flat(0., 1.));
55
56
57 // create the new Tree
58 m_validationTree = new TTree(m_validationTreeName.value().c_str(),
59 m_validationTreeDescription.value().c_str());
60
61 // the branches for the start
62 m_validationTree->Branch("RiddersSteps", &m_steps, "steps/I");
63 // loc 1
64 m_validationTree->Branch("Loc1Loc1", m_loc1loc1, "loc1loc1[steps]/F");
65 m_validationTree->Branch("Loc1Loc2", m_loc1loc2, "loc1loc2[steps]/F");
66 m_validationTree->Branch("Loc1Phi", m_loc1phi, "loc1phi[steps]/F");
67 m_validationTree->Branch("Loc1Theta", m_loc1theta, "loc1theta[steps]/F");
68 m_validationTree->Branch("Loc1qOp", m_loc1qop, "loc1qop[steps]/F");
69 m_validationTree->Branch("Loc1Steps", m_loc1steps, "loc1steps[steps]/F");
70 // loc 2
71 m_validationTree->Branch("Loc2Loc1", m_loc2loc1, "loc2loc1[steps]/F");
72 m_validationTree->Branch("Loc2Loc2", m_loc2loc2, "loc2loc2[steps]/F");
73 m_validationTree->Branch("Loc2Phi", m_loc2phi, "loc2phi[steps]/F");
74 m_validationTree->Branch("Loc2Theta", m_loc2theta, "loc2theta[steps]/F");
75 m_validationTree->Branch("Loc2qOp", m_loc2qop, "loc2qop[steps]/F");
76 m_validationTree->Branch("Loc2Steps", m_loc2steps, "loc2steps[steps]/F");
77 // phi
78 m_validationTree->Branch("PhiLoc1", m_philoc1, "philoc1[steps]/F");
79 m_validationTree->Branch("PhiLoc2", m_philoc2, "philoc2[steps]/F");
80 m_validationTree->Branch("PhiPhi", m_phiphi, "phiphi[steps]/F");
81 m_validationTree->Branch("PhiTheta", m_phitheta, "phitheta[steps]/F");
82 m_validationTree->Branch("PhiqOp", m_phiqop, "phiqop[steps]/F");
83 m_validationTree->Branch("PhiSteps", m_phisteps, "phisteps[steps]/F");
84 // Theta
85 m_validationTree->Branch("ThetaLoc1", m_thetaloc1, "thetaloc1[steps]/F");
86 m_validationTree->Branch("ThetaLoc2", m_thetaloc2, "thetaloc2[steps]/F");
87 m_validationTree->Branch("ThetaPhi", m_thetaphi, "thetaphi[steps]/F");
88 m_validationTree->Branch("ThetaTheta", m_thetatheta, "thetatheta[steps]/F");
89 m_validationTree->Branch("ThetaqOp", m_thetaqop, "thetaqop[steps]/F");
90 m_validationTree->Branch("ThetaSteps", m_thetasteps, "thetasteps[steps]/F");
91 // Qop
92 m_validationTree->Branch("QopLoc1", m_qoploc1, "qoploc1[steps]/F");
93 m_validationTree->Branch("QopLoc2", m_qoploc2, "qoploc2[steps]/F");
94 m_validationTree->Branch("QopPhi", m_qopphi, "qopphi[steps]/F");
95 m_validationTree->Branch("QopTheta", m_qoptheta, "qoptheta[steps]/F");
96 m_validationTree->Branch("QopqOp", m_qopqop, "qopqop[steps]/F");
97 m_validationTree->Branch("QopSteps", m_qopsteps, "qopsteps[steps]/F");
98
99 // now register the Tree
100 SmartIF<ITHistSvc> tHistSvc {
101 service("THistSvc")
102 };
103 if (!tHistSvc) {
104 ATH_MSG_ERROR("initialize() Could not find Hist Service -> Switching ValidationMode Off !");
105 delete m_validationTree;
106 m_validationTree = nullptr;
107 }
108 if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
109 ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !");
110 delete m_validationTree;
111 m_validationTree = nullptr;
112 }
113
114 if (m_localVariations.empty()) m_localVariations = {
115 0.01, 0.001, 0.0001
116 };
117
118 if (m_angularVariations.empty()) m_angularVariations = {
119 0.01, 0.001, 0.0001
120 };
121
122 if (m_qOpVariations.empty()) m_qOpVariations = {
123 0.0001, 0.00001, 0.000001
124 };
125
126 ATH_MSG_INFO("initialize() successful in ");
127 return StatusCode::SUCCESS;
128}
129
130//================ Finalisation =================================================
131
132StatusCode Trk::RiddersAlgorithm::finalize() {
133 // Code entered here will be executed once at the end of the program run.
134 return StatusCode::SUCCESS;
135}
136
137//================ Execution ====================================================
138
139StatusCode Trk::RiddersAlgorithm::execute(const EventContext& ctx) {
140 // this is fine
141 double p = m_minP + m_flatDist->shoot() * (m_maxP - m_minP);
142 double charge = (m_flatDist->shoot() > 0.5) ? -1. : 1.;
143 double qOverP = charge / p;
144
145 // for the momentum logging
146 // m_startP = p;
147
148 // the local start start
149 double loc1 = m_sigmaLoc * m_gaussDist->shoot();
150 double loc2 = m_sigmaLoc * m_gaussDist->shoot();
151 // are adopted for planar and straight line surfaces
152 double phi = m_minPhi + m_maxPhi * m_flatDist->shoot();
153 double eta = m_minEta + m_flatDist->shoot() * (m_maxEta - m_minEta);
154 double theta = 2. * atan(exp(-eta));
155
156 // start
157 double startR = std::abs(m_sigmaR * m_gaussDist->shoot());
158 double surfacePhi = M_PI * m_flatDist->shoot();
159
160 surfacePhi *= (m_flatDist->shoot() > 0.5) ? -1. : 1.;
161 double startX = startR * cos(surfacePhi);
162 double startY = startR * sin(surfacePhi);
163 double startZ = m_sigmaLoc * m_gaussDist->shoot();
164
165 // rotate it around Z
166 double alphaZ = M_PI * m_flatDist->shoot();
167 alphaZ *= (m_flatDist->shoot() > 0.5) ? -1. : 1.;
168
169 // create the plane surface
170 Trk::PlaneSurface startSurface(createTransform(startX,
171 startY,
172 startZ,
173 phi, theta,
174 alphaZ),
175 10e3, 10e3);
176
177 // make a covariance Matrix
178 AmgMatrix(5, 5) covMat;
179 covMat.setZero();
180
181 // the initial perigee with random numbers
182 Trk::AtaPlane startParameters(loc1,
183 loc2,
184 phi,
185 theta,
186 qOverP,
187 startSurface,
188 covMat);
189
190 ATH_MSG_VERBOSE("Start Parameters : " << startParameters);
191
192
193 // destination position
194 double estimationR = m_minimumR + (m_maximumR - m_minimumR) * m_flatDist->shoot();
195
196
197 // --------------- propagate to find a first intersection ---------------------
198 Trk::CylinderSurface estimationCylinder(Amg::Transform3D(), estimationR, 10e10);
199
200 ATH_MSG_VERBOSE("Cylinder to be intersected : " << estimationCylinder);
201
202 auto estimationParameters = m_propagator->propagateParameters(ctx,
203 startParameters,
204 estimationCylinder,
205 Trk::alongMomentum,
206 false,
207 *m_magFieldProperties);
208 if (!estimationParameters) {
209 ATH_MSG_VERBOSE("Estimation of intersection did not work - skip event !");
210 return StatusCode::SUCCESS;
211 }
212
213 ATH_MSG_VERBOSE("Estimation Parameters: " << *estimationParameters);
214
215 const Amg::Vector3D& estimatedPosition = estimationParameters->position();
216
217 double estimationX = estimatedPosition.x();
218 double estimationY = estimatedPosition.y();
219 double estimationZ = estimatedPosition.z();
220
221 double estimationPhi = estimatedPosition.phi();
222 double estimationTheta = estimatedPosition.theta();
223
224
225
226 double rotateTrans = M_PI * m_flatDist->shoot();
227 rotateTrans *= (m_flatDist->shoot() > 0.5) ? -1. : 1.;
228
229 Amg::Transform3D surfaceTransform;
230
231 if (m_useAlignedSurfaces) {
232 // create a radial vector
233 Amg::Vector3D radialVector(estimatedPosition.x(), estimatedPosition.y(), 0.);
234 Amg::Vector3D surfaceZdirection(radialVector.unit());
235 Amg::Vector3D surfaceYdirection(0., 0., 1.);
236 // the x direction
237 Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
238 // the rotation
239 Amg::RotationMatrix3D surfaceRotation;
240 surfaceRotation.col(0) = surfaceXdirection;
241 surfaceRotation.col(1) = surfaceYdirection;
242 surfaceRotation.col(2) = surfaceZdirection;
243 Amg::Transform3D nominalTransform(surfaceRotation, estimatedPosition);
244 surfaceTransform = Amg::Transform3D(nominalTransform * Amg::AngleAxis3D(rotateTrans, Amg::Vector3D(0., 0., 1.)));
245 } else
246 surfaceTransform = createTransform(estimationX,
247 estimationY,
248 estimationZ,
249 estimationPhi,
250 estimationTheta,
251 rotateTrans);
252
253 // cleanup for memory reasons
254
255 Trk::PlaneSurface destinationSurface(surfaceTransform, 10e5, 10e5);
256
257
258 // transport the start to the destination surface
259 std::optional<Trk::TransportJacobian> optTransportJacobian {};
260 AmgMatrix(5, 5) testMatrix;
261 testMatrix.setZero();
262 Trk::TransportJacobian currentStepJacobian(testMatrix);
263 double pathLimit = -1.;
264
265 auto trackParameters = m_propagator->propagate(ctx,
266 startParameters,
267 destinationSurface,
268 Trk::alongMomentum,
269 false,
270 *m_magFieldProperties,
271 optTransportJacobian,
272 pathLimit);
273
274 // --------------------- check if test propagation was successful ------------------------------
275 if (trackParameters && optTransportJacobian) {
276 unsigned int recStep = 0;
277 const auto& transportJacobian = (*optTransportJacobian);
278 // [0] Transport Jacobian -----------------------------------------------------
279 ATH_MSG_VERBOSE("TransportJacobian : " << transportJacobian);
280
281
282 // and now fill the variables
283 m_loc1loc1[recStep] = (transportJacobian) (0, 0);
284 m_loc1loc2[recStep] = (transportJacobian) (0, 1);
285 m_loc1phi[recStep] = (transportJacobian) (0, 2);
286 m_loc1theta[recStep] = (transportJacobian) (0, 3);
287 m_loc1qop[recStep] = (transportJacobian) (0, 4);
288 m_loc1steps[recStep] = 0.;
289
290 m_loc2loc1[recStep] = (transportJacobian) (1, 0);
291 m_loc2loc2[recStep] = (transportJacobian) (1, 1);
292 m_loc2phi[recStep] = (transportJacobian) (1, 2);
293 m_loc2theta[recStep] = (transportJacobian) (1, 3);
294 m_loc2qop[recStep] = (transportJacobian) (1, 4);
295 m_loc2steps[recStep] = 0.;
296
297 m_philoc1[recStep] = (transportJacobian) (2, 0);
298 m_philoc2[recStep] = (transportJacobian) (2, 1);
299 m_phiphi[recStep] = (transportJacobian) (2, 2);
300 m_phitheta[recStep] = (transportJacobian) (2, 3);
301 m_phiqop[recStep] = (transportJacobian) (2, 4);
302 m_phisteps[recStep] = 0.;
303
304 m_thetaloc1[recStep] = (transportJacobian) (3, 0);
305 m_thetaloc2[recStep] = (transportJacobian) (3, 1);
306 m_thetaphi[recStep] = (transportJacobian) (3, 2);
307 m_thetatheta[recStep] = (transportJacobian) (3, 3);
308 m_thetaqop[recStep] = (transportJacobian) (3, 4);
309 m_thetasteps[recStep] = 0.;
310
311 m_qoploc1[recStep] = (transportJacobian) (4, 0);
312 m_qoploc2[recStep] = (transportJacobian) (4, 1);
313 m_qopphi[recStep] = (transportJacobian) (4, 2);
314 m_qoptheta[recStep] = (transportJacobian) (4, 3);
315 m_qopqop[recStep] = (transportJacobian) (4, 4);
316 m_qopsteps[recStep] = 0.;
317
318 ++recStep;
319
320 // start the riddlers algorithm
321 // [1-2-3] Riddlers jacobians -----------------------------------------
322 for (unsigned int istep = 0; istep < m_localVariations.size(); ++istep) {
323 // --------------------------
324 ATH_MSG_VERBOSE("Performing step : " << istep);
325 // the initial perigee with random numbers
326 // for the first row
327 Trk::AtaPlane startLoc1Minus(loc1 - m_localVariations[istep], loc2, phi, theta, qOverP, startSurface);
328 Trk::AtaPlane startLoc1Plus(loc1 + m_localVariations[istep], loc2, phi, theta, qOverP, startSurface);
329 // for the second row
330 Trk::AtaPlane startLoc2Minus(loc1, loc2 - m_localVariations[istep], phi, theta, qOverP, startSurface);
331 Trk::AtaPlane startLoc2Plus(loc1, loc2 + m_localVariations[istep], phi, theta, qOverP, startSurface);
332 // for the third row
333 Trk::AtaPlane startPhiMinus(loc1, loc2, phi - m_angularVariations[istep], theta, qOverP, startSurface);
334 Trk::AtaPlane startPhiPlus(loc1, loc2, phi + m_angularVariations[istep], theta, qOverP, startSurface);
335 // for the fourth row
336 Trk::AtaPlane startThetaMinus(loc1, loc2, phi, theta - m_angularVariations[istep], qOverP, startSurface);
337 Trk::AtaPlane startThetaPlus(loc1, loc2, phi, theta + m_angularVariations[istep], qOverP, startSurface);
338 // for the fifth row
339 Trk::AtaPlane startQopMinus(loc1, loc2, phi, theta, qOverP - m_qOpVariations[istep], startSurface);
340 Trk::AtaPlane startQopPlus(loc1, loc2, phi, theta, qOverP + m_qOpVariations[istep], startSurface);
341
342 // the propagations --- 10 times
343 auto endLoc1Minus = m_propagator->propagateParameters(ctx,
344 startLoc1Minus,
345 destinationSurface,
346 Trk::alongMomentum,
347 false,
348 *m_magFieldProperties);
349
350
351 auto endLoc1Plus = m_propagator->propagateParameters(ctx,
352 startLoc1Plus,
353 destinationSurface,
354 Trk::alongMomentum,
355 false,
356 *m_magFieldProperties);
357
358 auto endLoc2Minus = m_propagator->propagateParameters(ctx,
359 startLoc2Minus,
360 destinationSurface,
361 Trk::alongMomentum,
362 false,
363 *m_magFieldProperties);
364
365 auto endLoc2Plus = m_propagator->propagateParameters(ctx,
366 startLoc2Plus,
367 destinationSurface,
368 Trk::alongMomentum,
369 false,
370 *m_magFieldProperties);
371
372 auto endPhiMinus = m_propagator->propagateParameters(ctx,
373 startPhiMinus,
374 destinationSurface,
375 Trk::alongMomentum,
376 false,
377 *m_magFieldProperties);
378
379 auto endPhiPlus = m_propagator->propagateParameters(ctx,
380 startPhiPlus,
381 destinationSurface,
382 Trk::alongMomentum,
383 false,
384 *m_magFieldProperties);
385
386 auto endThetaMinus = m_propagator->propagateParameters(ctx,
387 startThetaMinus,
388 destinationSurface,
389 Trk::alongMomentum,
390 false,
391 *m_magFieldProperties);
392
393 auto endThetaPlus = m_propagator->propagateParameters(ctx,
394 startThetaPlus,
395 destinationSurface,
396 Trk::alongMomentum,
397 false,
398 *m_magFieldProperties);
399
400 auto endQopMinus = m_propagator->propagateParameters(ctx,
401 startQopMinus,
402 destinationSurface,
403 Trk::alongMomentum,
404 false,
405 *m_magFieldProperties);
406
407 auto endQopPlus = m_propagator->propagateParameters(ctx,
408 startQopPlus,
409 destinationSurface,
410 Trk::alongMomentum,
411 false,
412 *m_magFieldProperties);
413 if (endLoc1Minus
414 && endLoc1Plus
415 && endLoc2Minus
416 && endLoc2Plus
417 && endPhiMinus
418 && endPhiPlus
419 && endThetaMinus
420 && endThetaPlus
421 && endQopMinus
422 && endQopPlus) {
423 // Loc 1
424 const Amg::VectorX& endLoc1MinusPar = endLoc1Minus->parameters();
425 const Amg::VectorX& endLoc1PlusPar = endLoc1Plus->parameters();
426 // Loc 2
427 const Amg::VectorX& endLoc2MinusPar = endLoc2Minus->parameters();
428 const Amg::VectorX& endLoc2PlusPar = endLoc2Plus->parameters();
429 // Phi
430 const Amg::VectorX& endPhiMinusPar = endPhiMinus->parameters();
431 const Amg::VectorX& endPhiPlusPar = endPhiPlus->parameters();
432 // Theta
433 const Amg::VectorX& endThetaMinusPar = endThetaMinus->parameters();
434 const Amg::VectorX& endThetaPlusPar = endThetaPlus->parameters();
435 // qop
436 const Amg::VectorX& endQopMinusPar = endQopMinus->parameters();
437 const Amg::VectorX& endQopPlusPar = endQopPlus->parameters();
438
439 // the deltas
440 Amg::VectorX endLoc1Diff(endLoc1PlusPar - endLoc1MinusPar);
441 Amg::VectorX endLoc2Diff(endLoc2PlusPar - endLoc2MinusPar);
442 Amg::VectorX endPhiDiff(endPhiPlusPar - endPhiMinusPar);
443 Amg::VectorX endThetaDiff(endThetaPlusPar - endThetaMinusPar);
444 Amg::VectorX endQopDiff(endQopPlusPar - endQopMinusPar);
445
446 currentStepJacobian(0, 0) = endLoc1Diff[0] / (2. * m_localVariations[istep]);
447 currentStepJacobian(0, 1) = endLoc2Diff[0] / (2. * m_localVariations[istep]);
448 currentStepJacobian(0, 2) = endPhiDiff[0] / (2. * m_angularVariations[istep]);
449 currentStepJacobian(0, 3) = endThetaDiff[0] / (2. * m_angularVariations[istep]);
450 currentStepJacobian(0, 4) = endQopDiff[0] / (2. * m_qOpVariations[istep]);
451
452 m_loc1loc1[recStep] = currentStepJacobian(0, 0);
453 m_loc1loc2[recStep] = currentStepJacobian(0, 1);
454 m_loc1phi[recStep] = currentStepJacobian(0, 2);
455 m_loc1theta[recStep] = currentStepJacobian(0, 3);
456 m_loc1qop[recStep] = currentStepJacobian(0, 4);
457 m_loc1steps[recStep] = m_localVariations[istep];
458
459 currentStepJacobian(1, 0) = endLoc1Diff[1] / (2. * m_localVariations[istep]);
460 currentStepJacobian(1, 1) = endLoc2Diff[1] / (2. * m_localVariations[istep]);
461 currentStepJacobian(1, 2) = endPhiDiff[1] / (2. * m_angularVariations[istep]);
462 currentStepJacobian(1, 3) = endThetaDiff[1] / (2. * m_angularVariations[istep]);
463 currentStepJacobian(1, 4) = endQopDiff[1] / (2. * m_qOpVariations[istep]);
464
465 m_loc2loc1[recStep] = currentStepJacobian(1, 0);
466 m_loc2loc2[recStep] = currentStepJacobian(1, 1);
467 m_loc2phi[recStep] = currentStepJacobian(1, 2);
468 m_loc2theta[recStep] = currentStepJacobian(1, 3);
469 m_loc2qop[recStep] = currentStepJacobian(1, 4);
470 m_loc2steps[recStep] = m_localVariations[istep];
471
472 currentStepJacobian(2, 0) = endLoc1Diff[2] / (2. * m_localVariations[istep]);
473 currentStepJacobian(2, 1) = endLoc2Diff[2] / (2. * m_localVariations[istep]);
474 currentStepJacobian(2, 2) = endPhiDiff[2] / (2. * m_angularVariations[istep]);
475 currentStepJacobian(2, 3) = endThetaDiff[2] / (2. * m_angularVariations[istep]);
476 currentStepJacobian(2, 4) = endQopDiff[2] / (2. * m_qOpVariations[istep]);
477
478 m_philoc1[recStep] = currentStepJacobian(2, 0);
479 m_philoc2[recStep] = currentStepJacobian(2, 1);
480 m_phiphi[recStep] = currentStepJacobian(2, 2);
481 m_phitheta[recStep] = currentStepJacobian(2, 3);
482 m_phiqop[recStep] = currentStepJacobian(2, 4);
483 m_phisteps[recStep] = m_angularVariations[istep];
484
485 currentStepJacobian(3, 0) = endLoc1Diff[3] / (2. * m_localVariations[istep]);
486 currentStepJacobian(3, 1) = endLoc2Diff[3] / (2. * m_localVariations[istep]);
487 currentStepJacobian(3, 2) = endPhiDiff[3] / (2. * m_angularVariations[istep]);
488 currentStepJacobian(3, 3) = endThetaDiff[3] / (2. * m_angularVariations[istep]);
489 currentStepJacobian(3, 4) = endQopDiff[3] / (2. * m_qOpVariations[istep]);
490
491 m_thetaloc1[recStep] = currentStepJacobian(3, 0);
492 m_thetaloc2[recStep] = currentStepJacobian(3, 1);
493 m_thetaphi[recStep] = currentStepJacobian(3, 2);
494 m_thetatheta[recStep] = currentStepJacobian(3, 3);
495 m_thetaqop[recStep] = currentStepJacobian(3, 4);
496 m_thetasteps[recStep] = m_angularVariations[istep];
497
498 currentStepJacobian(4, 0) = endLoc1Diff[4] / (2. * m_localVariations[istep]);
499 currentStepJacobian(4, 1) = endLoc2Diff[4] / (2. * m_localVariations[istep]);
500 currentStepJacobian(4, 2) = endPhiDiff[4] / (2. * m_angularVariations[istep]);
501 currentStepJacobian(4, 3) = endThetaDiff[4] / (2. * m_angularVariations[istep]);
502 currentStepJacobian(4, 4) = endQopDiff[4] / (2. * m_qOpVariations[istep]);
503
504 m_qoploc1[recStep] = currentStepJacobian(4, 0);
505 m_qoploc2[recStep] = currentStepJacobian(4, 1);
506 m_qopphi[recStep] = currentStepJacobian(4, 2);
507 m_qoptheta[recStep] = currentStepJacobian(4, 3);
508 m_qopqop[recStep] = currentStepJacobian(4, 4);
509 m_qopsteps[recStep] = m_qOpVariations[istep];
510
511 ATH_MSG_DEBUG("Current TransportJacobian : " << currentStepJacobian);
512
513 ++recStep;
514 }
515 }
516
517 // -------------------------------------------------------------------------------
518 if (recStep > 2) {
519 // The parabolic interpolation -------------------------------------------------------------------------------
520 // dL1
521 m_loc1loc1[recStep] =
522 parabolicInterpolation(m_loc1loc1[recStep - 1], m_loc1loc1[recStep - 2], m_loc1loc1[recStep - 3],
523 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
524 m_loc1loc2[recStep] =
525 parabolicInterpolation(m_loc1loc2[recStep - 1], m_loc1loc2[recStep - 2], m_loc1loc2[recStep - 3],
526 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
527 m_loc1phi[recStep] = parabolicInterpolation(m_loc1phi[recStep - 1], m_loc1phi[recStep - 2],
528 m_loc1phi[recStep - 3],
529 m_angularVariations[2], m_angularVariations[1],
530 m_angularVariations[0]);
531 m_loc1theta[recStep] =
532 parabolicInterpolation(m_loc1theta[recStep - 1], m_loc1theta[recStep - 2], m_loc1theta[recStep - 3],
533 m_angularVariations[2], m_angularVariations[1],
534 m_angularVariations[0]);
535 m_loc1qop[recStep] = parabolicInterpolation(m_loc1qop[recStep - 1], m_loc1qop[recStep - 2],
536 m_loc1qop[recStep - 3],
537 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
538 m_loc1steps[recStep] = 1;
539 // dL2
540 m_loc2loc1[recStep] =
541 parabolicInterpolation(m_loc2loc1[recStep - 1], m_loc2loc1[recStep - 2], m_loc2loc1[recStep - 3],
542 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
543 m_loc2loc2[recStep] =
544 parabolicInterpolation(m_loc2loc2[recStep - 1], m_loc2loc2[recStep - 2], m_loc2loc2[recStep - 3],
545 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
546 m_loc2phi[recStep] = parabolicInterpolation(m_loc2phi[recStep - 1], m_loc2phi[recStep - 2],
547 m_loc2phi[recStep - 3],
548 m_angularVariations[2], m_angularVariations[1],
549 m_angularVariations[0]);
550 m_loc2theta[recStep] =
551 parabolicInterpolation(m_loc2theta[recStep - 1], m_loc2theta[recStep - 2], m_loc2theta[recStep - 3],
552 m_angularVariations[2], m_angularVariations[1],
553 m_angularVariations[0]);
554 m_loc2qop[recStep] = parabolicInterpolation(m_loc2qop[recStep - 1], m_loc2qop[recStep - 2],
555 m_loc2qop[recStep - 3],
556 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
557 m_loc2steps[recStep] = 1;
558 // dPhi
559 m_philoc1[recStep] = parabolicInterpolation(m_philoc1[recStep - 1], m_philoc1[recStep - 2],
560 m_philoc1[recStep - 3],
561 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
562 m_philoc2[recStep] = parabolicInterpolation(m_philoc2[recStep - 1], m_philoc2[recStep - 2],
563 m_philoc2[recStep - 3],
564 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
565 m_phiphi[recStep] = parabolicInterpolation(m_phiphi[recStep - 1], m_phiphi[recStep - 2], m_phiphi[recStep - 3],
566 m_angularVariations[2], m_angularVariations[1],
567 m_angularVariations[0]);
568 m_phitheta[recStep] =
569 parabolicInterpolation(m_phitheta[recStep - 1], m_phitheta[recStep - 2], m_phitheta[recStep - 3],
570 m_angularVariations[2], m_angularVariations[1],
571 m_angularVariations[0]);
572 m_phiqop[recStep] = parabolicInterpolation(m_phiqop[recStep - 1], m_phiqop[recStep - 2], m_phiqop[recStep - 3],
573 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
574 m_phisteps[recStep] = 1;
575 // dTheta
576 m_thetaloc1[recStep] =
577 parabolicInterpolation(m_thetaloc1[recStep - 1], m_thetaloc1[recStep - 2], m_thetaloc1[recStep - 3],
578 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
579 m_thetaloc2[recStep] =
580 parabolicInterpolation(m_thetaloc2[recStep - 1], m_thetaloc2[recStep - 2], m_thetaloc2[recStep - 3],
581 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
582 m_thetaphi[recStep] =
583 parabolicInterpolation(m_thetaphi[recStep - 1], m_thetaphi[recStep - 2], m_thetaphi[recStep - 3],
584 m_angularVariations[2], m_angularVariations[1],
585 m_angularVariations[0]);
586 m_thetatheta[recStep] =
587 parabolicInterpolation(m_thetatheta[recStep - 1], m_thetatheta[recStep - 2], m_thetatheta[recStep - 3],
588 m_angularVariations[2], m_angularVariations[1],
589 m_angularVariations[0]);
590 m_thetaqop[recStep] =
591 parabolicInterpolation(m_thetaqop[recStep - 1], m_thetaqop[recStep - 2], m_thetaqop[recStep - 3],
592 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
593 m_thetasteps[recStep] = 1;
594 // dTheta
595 m_qoploc1[recStep] = parabolicInterpolation(m_qoploc1[recStep - 1], m_qoploc1[recStep - 2],
596 m_qoploc1[recStep - 3],
597 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
598 m_qoploc2[recStep] = parabolicInterpolation(m_qoploc2[recStep - 1], m_qoploc2[recStep - 2],
599 m_qoploc2[recStep - 3],
600 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
601 m_qopphi[recStep] = parabolicInterpolation(m_qopphi[recStep - 1], m_qopphi[recStep - 2], m_qopphi[recStep - 3],
602 m_angularVariations[2], m_angularVariations[1],
603 m_angularVariations[0]);
604 m_qoptheta[recStep] =
605 parabolicInterpolation(m_qoptheta[recStep - 1], m_qoptheta[recStep - 2], m_qoptheta[recStep - 3],
606 m_angularVariations[2], m_angularVariations[1],
607 m_angularVariations[0]);
608 m_qopqop[recStep] = parabolicInterpolation(m_qopqop[recStep - 1], m_qopqop[recStep - 2], m_qopqop[recStep - 3],
609 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
610 m_qopsteps[recStep] = 1;
611
612 currentStepJacobian(0, 0) = m_loc1loc1[recStep];
613 currentStepJacobian(0, 1) = m_loc1loc2[recStep];
614 currentStepJacobian(0, 2) = m_loc1phi[recStep];
615 currentStepJacobian(0, 3) = m_loc1theta[recStep];
616 currentStepJacobian(0, 4) = m_loc1qop[recStep];
617
618 currentStepJacobian(1, 0) = m_loc2loc1[recStep];
619 currentStepJacobian(1, 1) = m_loc2loc2[recStep];
620 currentStepJacobian(1, 2) = m_loc2phi[recStep];
621 currentStepJacobian(1, 3) = m_loc2theta[recStep];
622 currentStepJacobian(1, 4) = m_loc2qop[recStep];
623
624 currentStepJacobian(2, 0) = m_philoc1[recStep];
625 currentStepJacobian(2, 1) = m_philoc2[recStep];
626 currentStepJacobian(2, 2) = m_phiphi[recStep];
627 currentStepJacobian(2, 3) = m_phitheta[recStep];
628 currentStepJacobian(2, 4) = m_phiqop[recStep];
629
630 currentStepJacobian(3, 0) = m_thetaloc1[recStep];
631 currentStepJacobian(3, 1) = m_thetaloc2[recStep];
632 currentStepJacobian(3, 2) = m_thetaphi[recStep];
633 currentStepJacobian(3, 3) = m_thetatheta[recStep];
634 currentStepJacobian(3, 4) = m_thetaqop[recStep];
635
636 currentStepJacobian(4, 0) = m_qoploc1[recStep];
637 currentStepJacobian(4, 1) = m_qoploc2[recStep];
638 currentStepJacobian(4, 2) = m_qopphi[recStep];
639 currentStepJacobian(4, 3) = m_qoptheta[recStep];
640 currentStepJacobian(4, 4) = m_qopqop[recStep];
641 }
642
643 ATH_MSG_DEBUG("Interpolated TransportJacobian : " << currentStepJacobian);
644 ++recStep;
645
646
647 // now fill the results
648 TransportJacobian diffMatrix(transportJacobian - currentStepJacobian);
649
650 ATH_MSG_VERBOSE("Absolute Differences of the TransportJacobian : " << diffMatrix);
651
652 // Absolute differences ----------------------------------------------------
653 // (A1)
654 // fill the differences into the last one log (loc1)
655 m_loc1loc1[recStep] = diffMatrix(0, 0);
656 m_loc1loc2[recStep] = diffMatrix(0, 1);
657 m_loc1phi[recStep] = diffMatrix(0, 2);
658 m_loc1theta[recStep] = diffMatrix(0, 3);
659 m_loc1qop[recStep] = diffMatrix(0, 4);
660 m_loc1steps[recStep] = 2;
661 // fill the differences into the last one log (loc2)
662 m_loc2loc1[recStep] = diffMatrix(1, 0);
663 m_loc2loc2[recStep] = diffMatrix(1, 1);
664 m_loc2phi[recStep] = diffMatrix(1, 2);
665 m_loc2theta[recStep] = diffMatrix(1, 3);
666 m_loc2qop[recStep] = diffMatrix(1, 4);
667 m_loc2steps[recStep] = 2;
668 // fill the differences into the last one log (phi)
669 m_philoc1[recStep] = diffMatrix(2, 0);
670 m_philoc2[recStep] = diffMatrix(2, 1);
671 m_phiphi[recStep] = diffMatrix(2, 2);
672 m_phitheta[recStep] = diffMatrix(2, 3);
673 m_phiqop[recStep] = diffMatrix(2, 4);
674 m_phisteps[recStep] = 2;
675 // fill the differences into the last one log (theta)
676 m_thetaloc1[recStep] = diffMatrix(3, 0);
677 m_thetaloc2[recStep] = diffMatrix(3, 1);
678 m_thetaphi[recStep] = diffMatrix(3, 2);
679 m_thetatheta[recStep] = diffMatrix(3, 3);
680 m_thetaqop[recStep] = diffMatrix(3, 4);
681 m_thetasteps[recStep] = 2;
682 // fill the differences into the last one log (qop)
683 m_qoploc1[recStep] = diffMatrix(4, 0);
684 m_qoploc2[recStep] = diffMatrix(4, 1);
685 m_qopphi[recStep] = diffMatrix(4, 2);
686 m_qoptheta[recStep] = diffMatrix(4, 3);
687 m_qopqop[recStep] = diffMatrix(4, 4);
688 m_qopsteps[recStep] = 2;
689 ++recStep;
690
691 // (A2)
692 if (recStep > 1) {
693 // fill the differences into the last one log (loc1)
694 m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep - 1]) >
695 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep - 1])) : 0.;
696 m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep - 1]) >
697 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep - 1])) : 0.;
698 m_loc1phi[recStep] = std::abs(m_loc1phi[recStep - 1]) >
699 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep - 1])) : 0.;
700 m_loc1theta[recStep] = std::abs(m_loc1theta[recStep - 1]) > 1e-50 ? -std::log10(std::abs(
701 m_loc1theta[recStep -
702 1])) : 0.;
703 m_loc1qop[recStep] = std::abs(m_loc1qop[recStep - 1]) >
704 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep - 1])) : 0.;
705 m_loc1steps[recStep] = 3;
706 // fill the differences into the last one log (loc2)
707 m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep - 1]) >
708 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep - 1])) : 0.;
709 m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep - 1]) >
710 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep - 1])) : 0.;
711 m_loc2phi[recStep] = std::abs(m_loc2phi[recStep - 1]) >
712 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep - 1])) : 0.;
713 m_loc2theta[recStep] = std::abs(m_loc2theta[recStep - 1]) >
714 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep - 1])) : 0.;
715 m_loc2qop[recStep] = std::abs(m_loc2qop[recStep - 1]) >
716 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep - 1])) : 0.;
717 m_loc2steps[recStep] = 3;
718 // fill the differences into the last one log (phi)
719 m_philoc1[recStep] = std::abs(m_philoc1[recStep - 1]) >
720 1e-50 ? -std::log10(std::abs(m_philoc1[recStep - 1])) : 0.;
721 m_philoc2[recStep] = std::abs(m_philoc2[recStep - 1]) >
722 1e-50 ? -std::log10(std::abs(m_philoc2[recStep - 1])) : 0.;
723 m_phiphi[recStep] = std::abs(m_phiphi[recStep - 1]) >
724 1e-50 ? -std::log10(std::abs(m_phiphi[recStep - 1])) : 0.;
725 m_phitheta[recStep] = std::abs(m_phitheta[recStep - 1]) >
726 1e-50 ? -std::log10(std::abs(m_phitheta[recStep - 1])) : 0.;
727 m_phiqop[recStep] = std::abs(m_phiqop[recStep - 1]) >
728 1e-50 ? -std::log10(std::abs(m_phiqop[recStep - 1])) : 0.;
729 m_phisteps[recStep] = 3;
730 // fill the differences into the last one log (theta)
731 m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep - 1]) >
732 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep - 1])) : 0.;
733 m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep - 1]) >
734 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep - 1])) : 0.;
735 m_thetaphi[recStep] = std::abs(m_thetaphi[recStep - 1]) >
736 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep - 1])) : 0.;
737 m_thetatheta[recStep] = std::abs(m_thetatheta[recStep - 1]) >
738 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep - 1])) : 0.;
739 m_thetaqop[recStep] = std::abs(m_thetaqop[recStep - 1]) >
740 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep - 1])) : 0.;
741 m_thetasteps[recStep] = 3;
742 // fill the differences into the last one log (qop)
743 m_qoploc1[recStep] = std::abs(m_qoploc1[recStep - 1]) >
744 1e-50 ? -std::log10(std::abs(m_qoploc1[recStep - 1])) : 0.;
745 m_qoploc2[recStep] = std::abs(m_qoploc2[recStep - 1]) >
746 1e-50 ? -std::log10(std::abs(m_qoploc2[recStep - 1])) : 0.;
747 m_qopphi[recStep] = std::abs(m_qopphi[recStep - 1]) >
748 1e-50 ? -std::log10(std::abs(m_qopphi[recStep - 1])) : 0.;
749 m_qoptheta[recStep] = std::abs(m_qoptheta[recStep - 1]) > 1e-50 ? -std::log10(std::abs(
750 m_qoptheta[recStep -
751 1])) : 0.;
752 m_qopqop[recStep] = std::abs(m_qopqop[recStep - 1]) >
753 1e-50 ? -std::log10(std::abs(m_qopqop[recStep - 1])) : 0.;
754 m_qopsteps[recStep] = 3;
755 }
756 ++recStep;
757
758
759 // Relative differences) ----------------------------------------------------
760 // (R1)
761 // fill the differences into the last one log (loc1)
762 m_loc1loc1[recStep] =
763 std::abs((transportJacobian) (0, 0)) > 1e-50 ? diffMatrix(0, 0) / ((transportJacobian) (0, 0)) : 0.;
764 m_loc1loc2[recStep] =
765 std::abs((transportJacobian) (0, 1)) > 1e-50 ? diffMatrix(0, 1) / ((transportJacobian) (0, 1)) : 0.;
766 m_loc1phi[recStep] =
767 std::abs((transportJacobian) (0, 2)) > 1e-50 ? diffMatrix(0, 2) / ((transportJacobian) (0, 2)) : 0.;
768 m_loc1theta[recStep] =
769 std::abs((transportJacobian) (0, 3)) > 1e-50 ? diffMatrix(0, 3) / ((transportJacobian) (0, 3)) : 0.;
770 m_loc1qop[recStep] =
771 std::abs((transportJacobian) (0, 4)) > 1e-50 ? diffMatrix(0, 4) / ((transportJacobian) (0, 4)) : 0.;
772 m_loc1steps[recStep] = 4;
773 // fill the differences into the last one log (loc2)
774 m_loc2loc1[recStep] =
775 std::abs((transportJacobian) (1, 0)) > 1e-50 ? diffMatrix(1, 0) / ((transportJacobian) (1, 0)) : 0.;
776 m_loc2loc2[recStep] =
777 std::abs((transportJacobian) (1, 1)) > 1e-50 ? diffMatrix(1, 1) / ((transportJacobian) (1, 1)) : 0.;
778 m_loc2phi[recStep] =
779 std::abs((transportJacobian) (1, 2)) > 1e-50 ? diffMatrix(1, 2) / ((transportJacobian) (1, 2)) : 0.;
780 m_loc2theta[recStep] =
781 std::abs((transportJacobian) (1, 3)) > 1e-50 ? diffMatrix(1, 3) / ((transportJacobian) (1, 3)) : 0.;
782 m_loc2qop[recStep] =
783 std::abs((transportJacobian) (1, 4)) > 1e-50 ? diffMatrix(1, 4) / ((transportJacobian) (1, 4)) : 0.;
784 m_loc2steps[recStep] = 4;
785 // fill the differences into the last one log (phi)
786 m_philoc1[recStep] =
787 std::abs((transportJacobian) (2, 0)) > 1e-50 ? diffMatrix(2, 0) / ((transportJacobian) (2, 0)) : 0.;
788 m_philoc2[recStep] =
789 std::abs((transportJacobian) (2, 1)) > 1e-50 ? diffMatrix(2, 1) / ((transportJacobian) (2, 1)) : 0.;
790 m_phiphi[recStep] =
791 std::abs((transportJacobian) (2, 2)) > 1e-50 ? diffMatrix(2, 2) / ((transportJacobian) (2, 2)) : 0.;
792 m_phitheta[recStep] =
793 std::abs((transportJacobian) (2, 3)) > 1e-50 ? diffMatrix(2, 3) / ((transportJacobian) (2, 3)) : 0.;
794 m_phiqop[recStep] =
795 std::abs((transportJacobian) (2, 4)) > 1e-50 ? diffMatrix(2, 4) / ((transportJacobian) (2, 4)) : 0.;
796 m_phisteps[recStep] = 4;
797 // fill the differences into the last one log (theta)
798 m_thetaloc1[recStep] =
799 std::abs((transportJacobian) (3, 0)) > 1e-50 ? diffMatrix(3, 0) / ((transportJacobian) (3, 0)) : 0.;
800 m_thetaloc2[recStep] =
801 std::abs((transportJacobian) (3, 1)) > 1e-50 ? diffMatrix(3, 1) / ((transportJacobian) (3, 1)) : 0.;
802 m_thetaphi[recStep] =
803 std::abs((transportJacobian) (3, 2)) > 1e-50 ? diffMatrix(3, 2) / ((transportJacobian) (3, 2)) : 0.;
804 m_thetatheta[recStep] =
805 std::abs((transportJacobian) (3, 3)) > 1e-50 ? diffMatrix(3, 3) / ((transportJacobian) (3, 3)) : 0.;
806 m_thetaqop[recStep] =
807 std::abs((transportJacobian) (3, 4)) > 1e-50 ? diffMatrix(3, 4) / ((transportJacobian) (3, 4)) : 0.;
808 m_thetasteps[recStep] = 4;
809 // fill the differences into the last one log (qop)
810 m_qoploc1[recStep] =
811 std::abs((transportJacobian) (4, 0)) > 1e-50 ? diffMatrix(4, 0) / ((transportJacobian) (4, 0)) : 0.;
812 m_qoploc2[recStep] =
813 std::abs((transportJacobian) (4, 1)) > 1e-50 ? diffMatrix(4, 1) / ((transportJacobian) (4, 1)) : 0.;
814 m_qopphi[recStep] =
815 std::abs((transportJacobian) (4, 2)) > 1e-50 ? diffMatrix(4, 2) / ((transportJacobian) (4, 2)) : 0.;
816 m_qoptheta[recStep] =
817 std::abs((transportJacobian) (4, 3)) > 1e-50 ? diffMatrix(4, 3) / ((transportJacobian) (4, 3)) : 0.;
818 m_qopqop[recStep] =
819 std::abs((transportJacobian) (4, 4)) > 1e-50 ? diffMatrix(4, 4) / ((transportJacobian) (4, 4)) : 0.;
820 m_qopsteps[recStep] = 4;
821 ++recStep;
822
823 // (R2)
824 // Relative differences ----------------------------------------------------
825 m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep - 1]) >
826 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep - 1])) : 0.;
827 m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep - 1]) >
828 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep - 1])) : 0.;
829 m_loc1phi[recStep] = std::abs(m_loc1phi[recStep - 1]) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep - 1])) : 0.;
830 m_loc1theta[recStep] = std::abs(m_loc1theta[recStep - 1]) >
831 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep - 1])) : 0.;
832 m_loc1qop[recStep] = std::abs(m_loc1qop[recStep - 1]) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep - 1])) : 0.;
833 m_loc1steps[recStep] = 5;
834 // fill the differences into the last one log (loc2)
835 m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep - 1]) >
836 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep - 1])) : 0.;
837 m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep - 1]) >
838 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep - 1])) : 0.;
839 m_loc2phi[recStep] = std::abs(m_loc2phi[recStep - 1]) >
840 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep - 1])) : 0.;
841 m_loc2theta[recStep] = std::abs(m_loc2theta[recStep - 1]) > 1e-50 ? -std::log10(std::abs(
842 m_loc2theta[recStep -
843 1])) : 0.;
844 m_loc2qop[recStep] = std::abs(m_loc2qop[recStep - 1]) >
845 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep - 1])) : 0.;
846 m_loc2steps[recStep] = 5;
847 // fill the differences into the last one log (phi)
848 m_philoc1[recStep] = std::abs(m_philoc1[recStep - 1]) >
849 1e-50 ? -std::log10(std::abs(m_philoc1[recStep - 1])) : 0.;
850 m_philoc2[recStep] = std::abs(m_philoc2[recStep - 1]) >
851 1e-50 ? -std::log10(std::abs(m_philoc2[recStep - 1])) : 0.;
852 m_phiphi[recStep] = std::abs(m_phiphi[recStep - 1]) >
853 1e-50 ? -std::log10(std::abs(m_phiphi[recStep - 1])) : 0.;
854 m_phitheta[recStep] = std::abs(m_phitheta[recStep - 1]) >
855 1e-50 ? -std::log10(std::abs(m_phitheta[recStep - 1])) : 0.;
856 m_phiqop[recStep] = std::abs(m_phiqop[recStep - 1]) >
857 1e-50 ? -std::log10(std::abs(m_phiqop[recStep - 1])) : 0.;
858 m_phisteps[recStep] = 5;
859 // fill the differences into the last one log (theta)
860 m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep - 1]) > 1e-50 ? -std::log10(std::abs(
861 m_thetaloc1[recStep -
862 1])) : 0.;
863 m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep - 1]) > 1e-50 ? -std::log10(std::abs(
864 m_thetaloc2[recStep -
865 1])) : 0.;
866 m_thetaphi[recStep] = std::abs(m_thetaphi[recStep - 1]) >
867 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep - 1])) : 0.;
868 m_thetatheta[recStep] = std::abs(m_thetatheta[recStep - 1]) >
869 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep - 1])) : 0.;
870 m_thetaqop[recStep] = std::abs(m_thetaqop[recStep - 1]) >
871 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep - 1])) : 0.;
872 m_thetasteps[recStep] = 5;
873 // fill the differences into the last one log (qop)
874 m_qoploc1[recStep] = std::abs(m_qoploc1[recStep - 1]) >
875 1e-50 ? -std::log10(std::abs(m_qoploc1[recStep - 1])) : 0.;
876 m_qoploc2[recStep] = std::abs(m_qoploc2[recStep - 1]) >
877 1e-50 ? -std::log10(std::abs(m_qoploc2[recStep - 1])) : 0.;
878 m_qopphi[recStep] = std::abs(m_qopphi[recStep - 1]) > 1e-50 ? -std::log10(std::abs(m_qopphi[recStep - 1])) : 0.;
879 m_qoptheta[recStep] = std::abs(m_qoptheta[recStep - 1]) >
880 1e-50 ? -std::log10(std::abs(m_qoptheta[recStep - 1])) : 0.;
881 m_qopqop[recStep] = std::abs(m_qopqop[recStep - 1]) > 1e-50 ? -std::log10(std::abs(m_qopqop[recStep - 1])) : 0.;
882 m_qopsteps[recStep] = 5;
883 ++recStep;
884
885 m_steps = recStep;
886 m_validationTree->Fill();
887 }
888
889 // Code entered here will be executed once per event
890 return StatusCode::SUCCESS;
891}
892
893//============================================================================================
894Amg::Transform3D
895Trk::RiddersAlgorithm::createTransform(double x, double y, double z, double phi, double theta, double alphaZ) {
896 if (phi != 0. && theta != 0.) {
897 // create the Start Surface
898 Amg::Vector3D surfacePosition(x, y, z);
899 // z direction
900 Amg::Vector3D surfaceZdirection(cos(phi) * sin(theta),
901 sin(phi) * sin(theta),
902 cos(theta));
903 // the global z axis
904 Amg::Vector3D zAxis(0., 0., 1.);
905 // the y direction
906 Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
907 // the x direction
908 Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
909 // the rotation
910 Amg::RotationMatrix3D surfaceRotation;
911 surfaceRotation.col(0) = surfaceXdirection;
912 surfaceRotation.col(1) = surfaceYdirection;
913 surfaceRotation.col(2) = surfaceZdirection;
914 Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
915 return Amg::Transform3D(nominalTransform * Amg::AngleAxis3D(alphaZ, zAxis));
916 }
917
918 return Amg::Transform3D(Amg::Translation3D(x, y, z));
919}
920
921double Trk::RiddersAlgorithm::parabolicInterpolation(double y0, double y1, double y2,
922 double x0, double x1, double x2) {
923 double Z0 = x1 * x2 * y0;
924 double Z1 = x0 * x2 * y1;
925 double Z2 = x0 * x1 * y2;
926 double N0 = (x0 - x1) * (x0 - x2);
927 double N1 = (x1 - x2) * (x1 - x0);
928 double N2 = (x2 - x0) * (x2 - x1);
929
930 return Z0 / N0 + Z1 / N1 + Z2 / N2;
931}
M_PI
#define M_PI
Definition ActiveFraction.h:14
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x,...)
Definition AthMsgStreamMacros.h:43
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x,...)
Definition AthMsgStreamMacros.h:47
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x,...)
Definition AthMsgStreamMacros.h:42
ATH_MSG_INFO
#define ATH_MSG_INFO(x,...)
Definition AthMsgStreamMacros.h:45
charge
double charge(const T &p)
Definition AtlasPID.h:1003
AmgMatrix
#define AmgMatrix(rows, cols)
Definition EventPrimitives.h:49
AthAlgorithm::AthAlgorithm
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor.
Definition AthAlgorithm.cxx:10
Trk::CylinderSurface
Class for a CylinderSurface in the ATLAS detector.
Definition CylinderSurface.h:55
Trk::MagneticFieldProperties
magnetic field properties to steer the behavior of the extrapolation
Definition MagneticFieldProperties.h:31
Trk::PlaneSurface
Class for a planaer rectangular or trapezoidal surface in the ATLAS detector.
Definition PlaneSurface.h:64
Trk::RiddersAlgorithm::m_validationTreeFolder
StringProperty m_validationTreeFolder
Definition RiddersAlgorithm.h:109
Trk::RiddersAlgorithm::m_phisteps
float m_phisteps[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:132
Trk::RiddersAlgorithm::m_maximumR
DoubleProperty m_maximumR
Definition RiddersAlgorithm.h:93
Trk::RiddersAlgorithm::m_flatDist
Rndm::Numbers * m_flatDist
Definition RiddersAlgorithm.h:151
Trk::RiddersAlgorithm::m_minimumR
DoubleProperty m_minimumR
To create the first extimations.
Definition RiddersAlgorithm.h:92
Trk::RiddersAlgorithm::m_phiqop
float m_phiqop[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:131
Trk::RiddersAlgorithm::m_qopphi
float m_qopphi[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:143
Trk::RiddersAlgorithm::m_loc1theta
float m_loc1theta[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:116
Trk::RiddersAlgorithm::m_loc2theta
float m_loc2theta[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:123
Trk::RiddersAlgorithm::parabolicInterpolation
static double parabolicInterpolation(double y0, double y1, double y2, double x0, double x1, double x2)
Langrange-parabolic interpolation.
Definition RiddersAlgorithm.cxx:921
Trk::RiddersAlgorithm::m_validationTree
TTree * m_validationTree
Root Validation Tree.
Definition RiddersAlgorithm.h:100
Trk::RiddersAlgorithm::m_thetaloc1
float m_thetaloc1[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:134
Trk::RiddersAlgorithm::m_angularVariations
DoubleArrayProperty m_angularVariations
Definition RiddersAlgorithm.h:97
Trk::RiddersAlgorithm::m_gaussDist
Rndm::Numbers * m_gaussDist
Random Number setup.
Definition RiddersAlgorithm.h:150
Trk::RiddersAlgorithm::m_loc1loc1
float m_loc1loc1[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:113
Trk::RiddersAlgorithm::m_propagator
PublicToolHandle< IPropagator > m_propagator
member variables for algorithm properties:
Definition RiddersAlgorithm.h:70
Trk::RiddersAlgorithm::m_loc1qop
float m_loc1qop[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:117
Trk::RiddersAlgorithm::m_qoptheta
float m_qoptheta[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:144
Trk::RiddersAlgorithm::m_qopqop
float m_qopqop[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:145
Trk::RiddersAlgorithm::m_localVariations
DoubleArrayProperty m_localVariations
variations
Definition RiddersAlgorithm.h:96
Trk::RiddersAlgorithm::m_loc2steps
float m_loc2steps[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:125
Trk::RiddersAlgorithm::m_loc1loc2
float m_loc1loc2[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:114
Trk::RiddersAlgorithm::m_thetaloc2
float m_thetaloc2[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:135
Trk::RiddersAlgorithm::m_loc2qop
float m_loc2qop[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:124
Trk::RiddersAlgorithm::m_validationTreeDescription
StringProperty m_validationTreeDescription
Definition RiddersAlgorithm.h:106
Trk::RiddersAlgorithm::initialize
StatusCode initialize()
standard Athena-Algorithm method
Definition RiddersAlgorithm.cxx:37
Trk::RiddersAlgorithm::finalize
StatusCode finalize()
standard Athena-Algorithm method
Definition RiddersAlgorithm.cxx:132
Trk::RiddersAlgorithm::m_phiphi
float m_phiphi[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:129
Trk::RiddersAlgorithm::~RiddersAlgorithm
~RiddersAlgorithm()
Default Destructor.
Definition RiddersAlgorithm.cxx:29
Trk::RiddersAlgorithm::RiddersAlgorithm
RiddersAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Standard Athena-Algorithm Constructor.
Definition RiddersAlgorithm.cxx:24
Trk::RiddersAlgorithm::m_phitheta
float m_phitheta[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:130
Trk::RiddersAlgorithm::m_useCustomField
BooleanProperty m_useCustomField
Definition RiddersAlgorithm.h:71
Trk::RiddersAlgorithm::m_philoc2
float m_philoc2[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:128
Trk::RiddersAlgorithm::m_philoc1
float m_philoc1[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:127
Trk::RiddersAlgorithm::m_loc2loc2
float m_loc2loc2[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:121
Trk::RiddersAlgorithm::m_loc2phi
float m_loc2phi[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:122
Trk::RiddersAlgorithm::m_qoploc1
float m_qoploc1[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:141
Trk::RiddersAlgorithm::createTransform
static Amg::Transform3D createTransform(double x, double y, double z, double phi=0., double theta=0., double alphaZ=0.)
private helper method to create a HepTransform
Definition RiddersAlgorithm.cxx:895
Trk::RiddersAlgorithm::m_qoploc2
float m_qoploc2[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:142
Trk::RiddersAlgorithm::m_thetaqop
float m_thetaqop[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:138
Trk::RiddersAlgorithm::m_sigmaLoc
DoubleProperty m_sigmaLoc
The smearing.
Definition RiddersAlgorithm.h:79
Trk::RiddersAlgorithm::m_thetaphi
float m_thetaphi[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:136
Trk::RiddersAlgorithm::m_qOpVariations
DoubleArrayProperty m_qOpVariations
Definition RiddersAlgorithm.h:98
Trk::RiddersAlgorithm::m_validationTreeName
StringProperty m_validationTreeName
Definition RiddersAlgorithm.h:103
Trk::RiddersAlgorithm::m_qopsteps
float m_qopsteps[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:146
Trk::RiddersAlgorithm::execute
StatusCode execute(const EventContext &ctx)
standard Athena-Algorithm method
Definition RiddersAlgorithm.cxx:139
Trk::RiddersAlgorithm::m_loc1steps
float m_loc1steps[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:118
Trk::RiddersAlgorithm::m_thetasteps
float m_thetasteps[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:139
Trk::RiddersAlgorithm::m_loc1phi
float m_loc1phi[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:115
Trk::RiddersAlgorithm::m_thetatheta
float m_thetatheta[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:137
Trk::RiddersAlgorithm::m_magFieldProperties
MagneticFieldProperties * m_magFieldProperties
Definition RiddersAlgorithm.h:75
Trk::RiddersAlgorithm::m_useAlignedSurfaces
BooleanProperty m_useAlignedSurfaces
Definition RiddersAlgorithm.h:72
Trk::RiddersAlgorithm::m_fieldValue
DoubleProperty m_fieldValue
Definition RiddersAlgorithm.h:74
Trk::RiddersAlgorithm::m_loc2loc1
float m_loc2loc1[RIDDLERSSTEPS]
Definition RiddersAlgorithm.h:120
Trk::TransportJacobian
This class represents the jacobian for transforming initial track parameters to new track parameters ...
Definition TransportJacobian.h:46
Amg::AngleAxis3D
Eigen::AngleAxisd AngleAxis3D
Definition GeoPrimitives.h:45
Amg::RotationMatrix3D
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
Definition GeoPrimitives.h:50
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition GeoPrimitives.h:46
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:48
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition GeoPrimitives.h:44
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition EventPrimitives.h:30
Trk::alongMomentum
@ alongMomentum
Definition PropDirection.h:20
Trk::x
@ x
Definition ParamDefs.h:55
Trk::z
@ z
global position (cartesian)
Definition ParamDefs.h:57
Trk::theta
@ theta
Definition ParamDefs.h:66
Trk::qOverP
@ qOverP
perigee
Definition ParamDefs.h:67
Trk::y
@ y
Definition ParamDefs.h:56
Trk::loc2
@ loc2
generic first and second local coordinate
Definition ParamDefs.h:35
Trk::phi
@ phi
Definition ParamDefs.h:75
Trk::loc1
@ loc1
Definition ParamDefs.h:34
Trk::AtaPlane
ParametersT< TrackParametersDim, Charged, PlaneSurface > AtaPlane
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:34
Generated on for ATLAS Offline Software by doxygen 1.17.0