Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 
24 Trk::RiddersAlgorithm::RiddersAlgorithm(const std::string& name, ISvcLocator* pSvcLocator) :
25  AthAlgorithm(name,pSvcLocator) {}
26 
27 //================ Destructor =================================================
28 
29 Trk::RiddersAlgorithm::~RiddersAlgorithm()
30 {
31  delete m_gaussDist;
32  delete m_flatDist;
33  delete m_magFieldProperties;
34 }
35 
36 
37 //================ Initialisation =================================================
38 
39 StatusCode Trk::RiddersAlgorithm::initialize()
40 {
41  // Code entered here will be executed once at program start.
42  ATH_MSG_INFO( " initialize()" );
43 
44  ATH_CHECK( m_propagator.retrieve() );
45 
46  // Prepare the magnetic field properties
47  if (!m_useCustomField)
48  m_magFieldProperties = new Trk::MagneticFieldProperties();
49  else {
50  // the field
51  Amg::Vector3D bField(0.,0.,m_fieldValue);
52  // create the custom magnetic field
53  m_magFieldProperties = new Trk::MagneticFieldProperties(bField);
54  }
55 
56  // intialize the random number generators
57  m_gaussDist = new Rndm::Numbers(randSvc(), Rndm::Gauss(0.,1.));
58  m_flatDist = new Rndm::Numbers(randSvc(), Rndm::Flat(0.,1.));
59 
60 
61  // create the new Tree
62  m_validationTree = new TTree(m_validationTreeName.value().c_str(),
63  m_validationTreeDescription.value().c_str());
64 
65  // the branches for the start
66  m_validationTree->Branch("RiddersSteps", &m_steps, "steps/I");
67  // loc 1
68  m_validationTree->Branch("Loc1Loc1", m_loc1loc1, "loc1loc1[steps]/F");
69  m_validationTree->Branch("Loc1Loc2", m_loc1loc2, "loc1loc2[steps]/F");
70  m_validationTree->Branch("Loc1Phi", m_loc1phi, "loc1phi[steps]/F");
71  m_validationTree->Branch("Loc1Theta", m_loc1theta, "loc1theta[steps]/F");
72  m_validationTree->Branch("Loc1qOp", m_loc1qop, "loc1qop[steps]/F");
73  m_validationTree->Branch("Loc1Steps", m_loc1steps, "loc1steps[steps]/F");
74  // loc 2
75  m_validationTree->Branch("Loc2Loc1", m_loc2loc1, "loc2loc1[steps]/F");
76  m_validationTree->Branch("Loc2Loc2", m_loc2loc2, "loc2loc2[steps]/F");
77  m_validationTree->Branch("Loc2Phi", m_loc2phi, "loc2phi[steps]/F");
78  m_validationTree->Branch("Loc2Theta", m_loc2theta, "loc2theta[steps]/F");
79  m_validationTree->Branch("Loc2qOp", m_loc2qop, "loc2qop[steps]/F");
80  m_validationTree->Branch("Loc2Steps", m_loc2steps, "loc2steps[steps]/F");
81  // phi
82  m_validationTree->Branch("PhiLoc1", m_philoc1, "philoc1[steps]/F");
83  m_validationTree->Branch("PhiLoc2", m_philoc2 , "philoc2[steps]/F");
84  m_validationTree->Branch("PhiPhi", m_phiphi, "phiphi[steps]/F");
85  m_validationTree->Branch("PhiTheta", m_phitheta, "phitheta[steps]/F");
86  m_validationTree->Branch("PhiqOp", m_phiqop, "phiqop[steps]/F");
87  m_validationTree->Branch("PhiSteps", m_phisteps, "phisteps[steps]/F");
88  // Theta
89  m_validationTree->Branch("ThetaLoc1", m_thetaloc1, "thetaloc1[steps]/F");
90  m_validationTree->Branch("ThetaLoc2", m_thetaloc2, "thetaloc2[steps]/F");
91  m_validationTree->Branch("ThetaPhi", m_thetaphi, "thetaphi[steps]/F");
92  m_validationTree->Branch("ThetaTheta", m_thetatheta, "thetatheta[steps]/F");
93  m_validationTree->Branch("ThetaqOp", m_thetaqop, "thetaqop[steps]/F");
94  m_validationTree->Branch("ThetaSteps", m_thetasteps, "thetasteps[steps]/F");
95  // Qop
96  m_validationTree->Branch("QopLoc1", m_qoploc1, "qoploc1[steps]/F");
97  m_validationTree->Branch("QopLoc2", m_qoploc2, "qoploc2[steps]/F");
98  m_validationTree->Branch("QopPhi", m_qopphi, "qopphi[steps]/F");
99  m_validationTree->Branch("QopTheta", m_qoptheta, "qoptheta[steps]/F");
100  m_validationTree->Branch("QopqOp", m_qopqop, "qopqop[steps]/F");
101  m_validationTree->Branch("QopSteps", m_qopsteps, "qopsteps[steps]/F");
102 
103  // now register the Tree
104  SmartIF<ITHistSvc> tHistSvc{service("THistSvc")};
105  if (!tHistSvc){
106  ATH_MSG_ERROR( "initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
107  delete m_validationTree; m_validationTree = nullptr;
108  }
109  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
110  ATH_MSG_ERROR( "initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
111  delete m_validationTree; m_validationTree = nullptr;
112  }
113 
114  if (m_localVariations.empty())
115  m_localVariations = {0.01, 0.001, 0.0001};
116 
117  if (m_angularVariations.empty())
118  m_angularVariations = {0.01, 0.001, 0.0001};
119 
120  if (m_qOpVariations.empty())
121  m_qOpVariations = {0.0001, 0.00001, 0.000001};
122 
123  ATH_MSG_INFO( "initialize() successful in " );
124  return StatusCode::SUCCESS;
125 }
126 
127 //================ Finalisation =================================================
128 
129 StatusCode Trk::RiddersAlgorithm::finalize()
130 {
131  // Code entered here will be executed once at the end of the program run.
132  return StatusCode::SUCCESS;
133 }
134 
135 //================ Execution ====================================================
136 
137 StatusCode Trk::RiddersAlgorithm::execute()
138 {
139  const EventContext& ctx = Gaudi::Hive::currentContext();
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  surfacePhi *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
160  double startX = startR*cos(surfacePhi);
161  double startY = startR*sin(surfacePhi);
162  double startZ = m_sigmaLoc * m_gaussDist->shoot();
163 
164  // rotate it around Z
165  double alphaZ = M_PI * m_flatDist->shoot();
166  alphaZ *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
167 
168  // create the plane surface
169  Trk::PlaneSurface startSurface(createTransform(startX,
170  startY,
171  startZ,
172  phi, theta,
173  alphaZ),
174  10e3,10e3);
175 
176  // make a covariance Matrix
177  AmgMatrix(5,5) covMat;
178  covMat.setZero();
179 
180  // the initial perigee with random numbers
181  Trk::AtaPlane startParameters(loc1,
182  loc2,
183  phi,
184  theta,
185  qOverP,
186  startSurface,
187  covMat);
188 
189  ATH_MSG_VERBOSE( "Start Parameters : " << startParameters );
190 
191 
192  // destination position
193  double estimationR = m_minimumR + (m_maximumR-m_minimumR) * m_flatDist->shoot();
194 
195 
196  // --------------- propagate to find a first intersection ---------------------
197  Trk::CylinderSurface estimationCylinder(Amg::Transform3D(), estimationR, 10e10);
198 
199  ATH_MSG_VERBOSE( "Cylinder to be intersected : " << estimationCylinder );
200 
201  auto estimationParameters = m_propagator->propagateParameters(ctx,
202  startParameters,
203  estimationCylinder,
204  Trk::alongMomentum,
205  false,
206  *m_magFieldProperties);
207  if (!estimationParameters) {
208  ATH_MSG_VERBOSE( "Estimation of intersection did not work - skip event !" );
209  return StatusCode::SUCCESS;
210  }
211 
212  ATH_MSG_VERBOSE( "Estimation Parameters: " << *estimationParameters );
213 
214  const Amg::Vector3D& estimatedPosition = estimationParameters->position();
215 
216  double estimationX = estimatedPosition.x();
217  double estimationY = estimatedPosition.y();
218  double estimationZ = estimatedPosition.z();
219 
220  double estimationPhi = estimatedPosition.phi();
221  double estimationTheta = estimatedPosition.theta();
222 
223 
224 
225  double rotateTrans = M_PI * m_flatDist->shoot();
226  rotateTrans *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
227 
228  Amg::Transform3D surfaceTransform;
229 
230  if (m_useAlignedSurfaces){
231  // create a radial vector
232  Amg::Vector3D radialVector(estimatedPosition.x(), estimatedPosition.y(), 0.);
233  Amg::Vector3D surfaceZdirection(radialVector.unit());
234  Amg::Vector3D surfaceYdirection(0.,0.,1.);
235  // the x direction
236  Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
237  // the rotation
238  Amg::RotationMatrix3D surfaceRotation;
239  surfaceRotation.col(0) = surfaceXdirection;
240  surfaceRotation.col(1) = surfaceYdirection;
241  surfaceRotation.col(2) = surfaceZdirection;
242  Amg::Transform3D nominalTransform(surfaceRotation, estimatedPosition);
243  surfaceTransform = Amg::Transform3D(nominalTransform*Amg::AngleAxis3D(rotateTrans,Amg::Vector3D(0.,0.,1.)));
244  } else
245  surfaceTransform = createTransform(estimationX,
246  estimationY,
247  estimationZ,
248  estimationPhi,
249  estimationTheta,
250  rotateTrans);
251 
252  // cleanup for memory reasons
253 
254  Trk::PlaneSurface destinationSurface(surfaceTransform,10e5 , 10e5);
255 
256 
257  // transport the start to the destination surface
258  std::optional<Trk::TransportJacobian> optTransportJacobian{};
259  AmgMatrix(5,5) testMatrix; testMatrix.setZero();
260  Trk::TransportJacobian currentStepJacobian(testMatrix);
261  double pathLimit = -1.;
262 
263  auto trackParameters = m_propagator->propagate(ctx,
264  startParameters,
265  destinationSurface,
266  Trk::alongMomentum,
267  false,
268  *m_magFieldProperties,
269  optTransportJacobian,
270  pathLimit);
271 
272  // --------------------- check if test propagation was successful ------------------------------
273  if (trackParameters && optTransportJacobian){
274 
275  unsigned int recStep = 0;
276  const auto& transportJacobian = (*optTransportJacobian);
277  // [0] Transport Jacobian -----------------------------------------------------
278  ATH_MSG_VERBOSE( "TransportJacobian : " << transportJacobian );
279 
280 
281  // and now fill the variables
282  m_loc1loc1[recStep] = (transportJacobian)(0,0);
283  m_loc1loc2[recStep] = (transportJacobian)(0,1);
284  m_loc1phi[recStep] = (transportJacobian)(0,2);
285  m_loc1theta[recStep] = (transportJacobian)(0,3);
286  m_loc1qop[recStep] = (transportJacobian)(0,4);
287  m_loc1steps[recStep] = 0.;
288 
289  m_loc2loc1[recStep] = (transportJacobian)(1,0);
290  m_loc2loc2[recStep] = (transportJacobian)(1,1);
291  m_loc2phi[recStep] = (transportJacobian)(1,2);
292  m_loc2theta[recStep] = (transportJacobian)(1,3);
293  m_loc2qop[recStep] = (transportJacobian)(1,4);
294  m_loc2steps[recStep] = 0.;
295 
296  m_philoc1[recStep] = (transportJacobian)(2,0);
297  m_philoc2[recStep] = (transportJacobian)(2,1);
298  m_phiphi[recStep] = (transportJacobian)(2,2);
299  m_phitheta[recStep] = (transportJacobian)(2,3);
300  m_phiqop[recStep] = (transportJacobian)(2,4);
301  m_phisteps[recStep] = 0.;
302 
303  m_thetaloc1[recStep] = (transportJacobian)(3,0);
304  m_thetaloc2[recStep] = (transportJacobian)(3,1);
305  m_thetaphi[recStep] = (transportJacobian)(3,2);
306  m_thetatheta[recStep] = (transportJacobian)(3,3);
307  m_thetaqop[recStep] = (transportJacobian)(3,4);
308  m_thetasteps[recStep] = 0.;
309 
310  m_qoploc1[recStep] = (transportJacobian)(4,0);
311  m_qoploc2[recStep] = (transportJacobian)(4,1);
312  m_qopphi[recStep] = (transportJacobian)(4,2);
313  m_qoptheta[recStep] = (transportJacobian)(4,3);
314  m_qopqop[recStep] = (transportJacobian)(4,4);
315  m_qopsteps[recStep] = 0.;
316 
317  ++recStep;
318 
319  // start the riddlers algorithm
320  // [1-2-3] Riddlers jacobians -----------------------------------------
321  for (unsigned int istep = 0; istep < m_localVariations.size(); ++istep){
322  // --------------------------
323  ATH_MSG_VERBOSE( "Performing step : " << istep );
324  // the initial perigee with random numbers
325  // for the first row
326  Trk::AtaPlane startLoc1Minus(loc1-m_localVariations[istep],loc2,phi,theta,qOverP,startSurface);
327  Trk::AtaPlane startLoc1Plus(loc1+m_localVariations[istep],loc2,phi,theta,qOverP,startSurface);
328  // for the second row
329  Trk::AtaPlane startLoc2Minus(loc1,loc2-m_localVariations[istep],phi,theta,qOverP,startSurface);
330  Trk::AtaPlane startLoc2Plus(loc1,loc2+m_localVariations[istep],phi,theta,qOverP,startSurface);
331  // for the third row
332  Trk::AtaPlane startPhiMinus(loc1,loc2,phi-m_angularVariations[istep],theta,qOverP,startSurface);
333  Trk::AtaPlane startPhiPlus(loc1,loc2,phi+m_angularVariations[istep],theta,qOverP,startSurface);
334  // for the fourth row
335  Trk::AtaPlane startThetaMinus(loc1,loc2,phi,theta-m_angularVariations[istep],qOverP,startSurface);
336  Trk::AtaPlane startThetaPlus(loc1,loc2,phi,theta+m_angularVariations[istep],qOverP,startSurface);
337  // for the fifth row
338  Trk::AtaPlane startQopMinus(loc1,loc2,phi,theta,qOverP-m_qOpVariations[istep],startSurface);
339  Trk::AtaPlane startQopPlus(loc1,loc2,phi,theta,qOverP+m_qOpVariations[istep],startSurface);
340 
341  // the propagations --- 10 times
342  auto endLoc1Minus = m_propagator->propagateParameters(ctx,
343  startLoc1Minus,
344  destinationSurface,
345  Trk::alongMomentum,
346  false,
347  *m_magFieldProperties);
348 
349 
350  auto endLoc1Plus = m_propagator->propagateParameters(ctx,
351  startLoc1Plus,
352  destinationSurface,
353  Trk::alongMomentum,
354  false,
355  *m_magFieldProperties);
356 
357  auto endLoc2Minus = m_propagator->propagateParameters(ctx,
358  startLoc2Minus,
359  destinationSurface,
360  Trk::alongMomentum,
361  false,
362  *m_magFieldProperties);
363 
364  auto endLoc2Plus = m_propagator->propagateParameters(ctx,
365  startLoc2Plus,
366  destinationSurface,
367  Trk::alongMomentum,
368  false,
369  *m_magFieldProperties);
370 
371  auto endPhiMinus = m_propagator->propagateParameters(ctx,
372  startPhiMinus,
373  destinationSurface,
374  Trk::alongMomentum,
375  false,
376  *m_magFieldProperties);
377 
378  auto endPhiPlus = m_propagator->propagateParameters(ctx,
379  startPhiPlus,
380  destinationSurface,
381  Trk::alongMomentum,
382  false,
383  *m_magFieldProperties);
384 
385  auto endThetaMinus = m_propagator->propagateParameters(ctx,
386  startThetaMinus,
387  destinationSurface,
388  Trk::alongMomentum,
389  false,
390  *m_magFieldProperties);
391 
392  auto endThetaPlus = m_propagator->propagateParameters(ctx,
393  startThetaPlus,
394  destinationSurface,
395  Trk::alongMomentum,
396  false,
397  *m_magFieldProperties);
398 
399  auto endQopMinus = m_propagator->propagateParameters(ctx,
400  startQopMinus,
401  destinationSurface,
402  Trk::alongMomentum,
403  false,
404  *m_magFieldProperties);
405 
406  auto endQopPlus = m_propagator->propagateParameters(ctx,
407  startQopPlus,
408  destinationSurface,
409  Trk::alongMomentum,
410  false,
411  *m_magFieldProperties);
412  if (endLoc1Minus
413  && endLoc1Plus
414  && endLoc2Minus
415  && endLoc2Plus
416  && endPhiMinus
417  && endPhiPlus
418  && endThetaMinus
419  && endThetaPlus
420  && endQopMinus
421  && endQopPlus){
422 
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] = parabolicInterpolation(m_loc1loc1[recStep-1],m_loc1loc1[recStep-2],m_loc1loc1[recStep-3],
522  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
523  m_loc1loc2[recStep] = parabolicInterpolation(m_loc1loc2[recStep-1],m_loc1loc2[recStep-2],m_loc1loc2[recStep-3],
524  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
525  m_loc1phi[recStep] = parabolicInterpolation(m_loc1phi[recStep-1],m_loc1phi[recStep-2],m_loc1phi[recStep-3],
526  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
527  m_loc1theta[recStep] = parabolicInterpolation(m_loc1theta[recStep-1],m_loc1theta[recStep-2],m_loc1theta[recStep-3],
528  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
529  m_loc1qop[recStep] = parabolicInterpolation(m_loc1qop[recStep-1],m_loc1qop[recStep-2],m_loc1qop[recStep-3],
530  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
531  m_loc1steps[recStep] = 1;
532  // dL2
533  m_loc2loc1[recStep] = parabolicInterpolation(m_loc2loc1[recStep-1],m_loc2loc1[recStep-2],m_loc2loc1[recStep-3],
534  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
535  m_loc2loc2[recStep] = parabolicInterpolation(m_loc2loc2[recStep-1],m_loc2loc2[recStep-2],m_loc2loc2[recStep-3],
536  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
537  m_loc2phi[recStep] = parabolicInterpolation(m_loc2phi[recStep-1],m_loc2phi[recStep-2],m_loc2phi[recStep-3],
538  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
539  m_loc2theta[recStep] = parabolicInterpolation(m_loc2theta[recStep-1],m_loc2theta[recStep-2],m_loc2theta[recStep-3],
540  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
541  m_loc2qop[recStep] = parabolicInterpolation(m_loc2qop[recStep-1],m_loc2qop[recStep-2],m_loc2qop[recStep-3],
542  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
543  m_loc2steps[recStep] = 1;
544  // dPhi
545  m_philoc1[recStep] = parabolicInterpolation(m_philoc1[recStep-1],m_philoc1[recStep-2],m_philoc1[recStep-3],
546  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
547  m_philoc2[recStep] = parabolicInterpolation(m_philoc2[recStep-1],m_philoc2[recStep-2],m_philoc2[recStep-3],
548  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
549  m_phiphi[recStep] = parabolicInterpolation(m_phiphi[recStep-1],m_phiphi[recStep-2],m_phiphi[recStep-3],
550  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
551  m_phitheta[recStep] = parabolicInterpolation(m_phitheta[recStep-1],m_phitheta[recStep-2],m_phitheta[recStep-3],
552  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
553  m_phiqop[recStep] = parabolicInterpolation(m_phiqop[recStep-1],m_phiqop[recStep-2],m_phiqop[recStep-3],
554  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
555  m_phisteps[recStep] = 1;
556  // dTheta
557  m_thetaloc1[recStep] = parabolicInterpolation(m_thetaloc1[recStep-1],m_thetaloc1[recStep-2],m_thetaloc1[recStep-3],
558  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
559  m_thetaloc2[recStep] = parabolicInterpolation(m_thetaloc2[recStep-1],m_thetaloc2[recStep-2],m_thetaloc2[recStep-3],
560  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
561  m_thetaphi[recStep] = parabolicInterpolation(m_thetaphi[recStep-1],m_thetaphi[recStep-2],m_thetaphi[recStep-3],
562  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
563  m_thetatheta[recStep] = parabolicInterpolation(m_thetatheta[recStep-1],m_thetatheta[recStep-2],m_thetatheta[recStep-3],
564  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
565  m_thetaqop[recStep] = parabolicInterpolation(m_thetaqop[recStep-1],m_thetaqop[recStep-2],m_thetaqop[recStep-3],
566  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
567  m_thetasteps[recStep] = 1;
568  // dTheta
569  m_qoploc1[recStep] = parabolicInterpolation(m_qoploc1[recStep-1],m_qoploc1[recStep-2],m_qoploc1[recStep-3],
570  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
571  m_qoploc2[recStep] = parabolicInterpolation(m_qoploc2[recStep-1],m_qoploc2[recStep-2],m_qoploc2[recStep-3],
572  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
573  m_qopphi[recStep] = parabolicInterpolation(m_qopphi[recStep-1],m_qopphi[recStep-2],m_qopphi[recStep-3],
574  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
575  m_qoptheta[recStep] = parabolicInterpolation(m_qoptheta[recStep-1],m_qoptheta[recStep-2],m_qoptheta[recStep-3],
576  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
577  m_qopqop[recStep] = parabolicInterpolation(m_qopqop[recStep-1],m_qopqop[recStep-2],m_qopqop[recStep-3],
578  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
579  m_qopsteps[recStep] = 1;
580 
581  currentStepJacobian(0,0)= m_loc1loc1[recStep];
582  currentStepJacobian(0,1)= m_loc1loc2[recStep];
583  currentStepJacobian(0,2)= m_loc1phi[recStep];
584  currentStepJacobian(0,3)= m_loc1theta[recStep];
585  currentStepJacobian(0,4)= m_loc1qop[recStep];
586 
587  currentStepJacobian(1,0)= m_loc2loc1[recStep];
588  currentStepJacobian(1,1)= m_loc2loc2[recStep];
589  currentStepJacobian(1,2)= m_loc2phi[recStep];
590  currentStepJacobian(1,3)= m_loc2theta[recStep];
591  currentStepJacobian(1,4)= m_loc2qop[recStep];
592 
593  currentStepJacobian(2,0)= m_philoc1[recStep];
594  currentStepJacobian(2,1)= m_philoc2[recStep];
595  currentStepJacobian(2,2)= m_phiphi[recStep];
596  currentStepJacobian(2,3)= m_phitheta[recStep];
597  currentStepJacobian(2,4)= m_phiqop[recStep];
598 
599  currentStepJacobian(3,0)= m_thetaloc1[recStep];
600  currentStepJacobian(3,1)= m_thetaloc2[recStep];
601  currentStepJacobian(3,2)= m_thetaphi[recStep];
602  currentStepJacobian(3,3)= m_thetatheta[recStep];
603  currentStepJacobian(3,4)= m_thetaqop[recStep];
604 
605  currentStepJacobian(4,0)= m_qoploc1[recStep];
606  currentStepJacobian(4,1)= m_qoploc2[recStep];
607  currentStepJacobian(4,2)= m_qopphi[recStep];
608  currentStepJacobian(4,3)= m_qoptheta[recStep];
609  currentStepJacobian(4,4)= m_qopqop[recStep];
610 
611  }
612 
613  ATH_MSG_DEBUG( "Interpolated TransportJacobian : " << currentStepJacobian );
614  ++recStep;
615 
616 
617  // now fill the results
618  TransportJacobian diffMatrix(transportJacobian-currentStepJacobian);
619 
620  ATH_MSG_VERBOSE( "Absolute Differences of the TransportJacobian : " << diffMatrix );
621 
622  // Absolute differences ----------------------------------------------------
623  // (A1)
624  // fill the differences into the last one log (loc1)
625  m_loc1loc1[recStep] = diffMatrix(0,0);
626  m_loc1loc2[recStep] = diffMatrix(0,1);
627  m_loc1phi[recStep] = diffMatrix(0,2);
628  m_loc1theta[recStep] = diffMatrix(0,3);
629  m_loc1qop[recStep] = diffMatrix(0,4);
630  m_loc1steps[recStep] = 2;
631  // fill the differences into the last one log (loc2)
632  m_loc2loc1[recStep] = diffMatrix(1,0);
633  m_loc2loc2[recStep] = diffMatrix(1,1);
634  m_loc2phi[recStep] = diffMatrix(1,2);
635  m_loc2theta[recStep] = diffMatrix(1,3);
636  m_loc2qop[recStep] = diffMatrix(1,4);
637  m_loc2steps[recStep] = 2;
638  // fill the differences into the last one log (phi)
639  m_philoc1[recStep] = diffMatrix(2,0);
640  m_philoc2[recStep] = diffMatrix(2,1);
641  m_phiphi[recStep] = diffMatrix(2,2);
642  m_phitheta[recStep] = diffMatrix(2,3);
643  m_phiqop[recStep] = diffMatrix(2,4);
644  m_phisteps[recStep] = 2;
645  // fill the differences into the last one log (theta)
646  m_thetaloc1[recStep] = diffMatrix(3,0);
647  m_thetaloc2[recStep] = diffMatrix(3,1);
648  m_thetaphi[recStep] = diffMatrix(3,2);
649  m_thetatheta[recStep] = diffMatrix(3,3);
650  m_thetaqop[recStep] = diffMatrix(3,4);
651  m_thetasteps[recStep] = 2;
652  // fill the differences into the last one log (qop)
653  m_qoploc1[recStep] = diffMatrix(4,0);
654  m_qoploc2[recStep] = diffMatrix(4,1);
655  m_qopphi[recStep] = diffMatrix(4,2);
656  m_qoptheta[recStep] = diffMatrix(4,3);
657  m_qopqop[recStep] = diffMatrix(4,4);
658  m_qopsteps[recStep] = 2;
659  ++recStep;
660 
661  // (A2)
662  if (recStep > 1){
663  // fill the differences into the last one log (loc1)
664  m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep-1] )) : 0.;
665  m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep-1] )) : 0.;
666  m_loc1phi[recStep] = std::abs(m_loc1phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep-1] )) : 0.;
667  m_loc1theta[recStep] = std::abs(m_loc1theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep-1])) : 0.;
668  m_loc1qop[recStep] = std::abs(m_loc1qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep-1] )) : 0.;
669  m_loc1steps[recStep] = 3;
670  // fill the differences into the last one log (loc2)
671  m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep-1] )) : 0.;
672  m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep-1] )) : 0.;
673  m_loc2phi[recStep] = std::abs(m_loc2phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep-1] )) : 0.;
674  m_loc2theta[recStep] = std::abs(m_loc2theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep-1])) : 0.;
675  m_loc2qop[recStep] = std::abs(m_loc2qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep-1] )) : 0.;
676  m_loc2steps[recStep] = 3;
677  // fill the differences into the last one log (phi)
678  m_philoc1[recStep] = std::abs(m_philoc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc1[recStep-1] )) : 0.;
679  m_philoc2[recStep] = std::abs(m_philoc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc2[recStep-1] )) : 0.;
680  m_phiphi[recStep] = std::abs(m_phiphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiphi[recStep-1] )) : 0.;
681  m_phitheta[recStep] = std::abs(m_phitheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phitheta[recStep-1])) : 0.;
682  m_phiqop[recStep] = std::abs(m_phiqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiqop[recStep-1] )) : 0.;
683  m_phisteps[recStep] = 3;
684  // fill the differences into the last one log (theta)
685  m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep-1] )) : 0.;
686  m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep-1] )) : 0.;
687  m_thetaphi[recStep] = std::abs(m_thetaphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep-1] )) : 0.;
688  m_thetatheta[recStep] = std::abs(m_thetatheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep-1])) : 0.;
689  m_thetaqop[recStep] = std::abs(m_thetaqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep-1] )) : 0.;
690  m_thetasteps[recStep] = 3;
691  // fill the differences into the last one log (qop)
692  m_qoploc1[recStep] = std::abs(m_qoploc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoploc1[recStep-1] )) : 0.;
693  m_qoploc2[recStep] = std::abs(m_qoploc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoploc2[recStep-1] )) : 0.;
694  m_qopphi[recStep] = std::abs(m_qopphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qopphi[recStep-1] )) : 0.;
695  m_qoptheta[recStep] = std::abs(m_qoptheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoptheta[recStep-1])) : 0.;
696  m_qopqop[recStep] = std::abs(m_qopqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qopqop[recStep-1] )) : 0.;
697  m_qopsteps[recStep] = 3;
698  }
699  ++recStep;
700 
701 
702  // Relative differences) ----------------------------------------------------
703  // (R1)
704  // fill the differences into the last one log (loc1)
705  m_loc1loc1[recStep] = std::abs((transportJacobian)(0,0)) > 1e-50 ? diffMatrix(0,0)/((transportJacobian)(0,0)) : 0.;
706  m_loc1loc2[recStep] = std::abs((transportJacobian)(0,1)) > 1e-50 ? diffMatrix(0,1)/((transportJacobian)(0,1)) : 0.;
707  m_loc1phi[recStep] = std::abs((transportJacobian)(0,2)) > 1e-50 ? diffMatrix(0,2)/((transportJacobian)(0,2)) : 0.;
708  m_loc1theta[recStep] = std::abs((transportJacobian)(0,3)) > 1e-50 ? diffMatrix(0,3)/((transportJacobian)(0,3)) : 0.;
709  m_loc1qop[recStep] = std::abs((transportJacobian)(0,4)) > 1e-50 ? diffMatrix(0,4)/((transportJacobian)(0,4)) : 0.;
710  m_loc1steps[recStep] = 4;
711  // fill the differences into the last one log (loc2)
712  m_loc2loc1[recStep] = std::abs((transportJacobian)(1,0)) > 1e-50 ? diffMatrix(1,0)/((transportJacobian)(1,0)) : 0.;
713  m_loc2loc2[recStep] = std::abs((transportJacobian)(1,1)) > 1e-50 ? diffMatrix(1,1)/((transportJacobian)(1,1)) : 0.;
714  m_loc2phi[recStep] = std::abs((transportJacobian)(1,2)) > 1e-50 ? diffMatrix(1,2)/((transportJacobian)(1,2)) : 0.;
715  m_loc2theta[recStep] = std::abs((transportJacobian)(1,3)) > 1e-50 ? diffMatrix(1,3)/((transportJacobian)(1,3)) : 0.;
716  m_loc2qop[recStep] = std::abs((transportJacobian)(1,4)) > 1e-50 ? diffMatrix(1,4)/((transportJacobian)(1,4)) : 0.;
717  m_loc2steps[recStep] = 4;
718  // fill the differences into the last one log (phi)
719  m_philoc1[recStep] = std::abs((transportJacobian)(2,0)) > 1e-50 ? diffMatrix(2,0)/((transportJacobian)(2,0)) : 0.;
720  m_philoc2[recStep] = std::abs((transportJacobian)(2,1)) > 1e-50 ? diffMatrix(2,1)/((transportJacobian)(2,1)) : 0.;
721  m_phiphi[recStep] = std::abs((transportJacobian)(2,2)) > 1e-50 ? diffMatrix(2,2)/((transportJacobian)(2,2)) : 0.;
722  m_phitheta[recStep] = std::abs((transportJacobian)(2,3)) > 1e-50 ? diffMatrix(2,3)/((transportJacobian)(2,3)) : 0.;
723  m_phiqop[recStep] = std::abs((transportJacobian)(2,4)) > 1e-50 ? diffMatrix(2,4)/((transportJacobian)(2,4)) : 0.;
724  m_phisteps[recStep] = 4;
725  // fill the differences into the last one log (theta)
726  m_thetaloc1[recStep] = std::abs((transportJacobian)(3,0)) > 1e-50 ? diffMatrix(3,0)/((transportJacobian)(3,0)) : 0.;
727  m_thetaloc2[recStep] = std::abs((transportJacobian)(3,1)) > 1e-50 ? diffMatrix(3,1)/((transportJacobian)(3,1)) : 0.;
728  m_thetaphi[recStep] = std::abs((transportJacobian)(3,2)) > 1e-50 ? diffMatrix(3,2)/((transportJacobian)(3,2)) : 0.;
729  m_thetatheta[recStep] = std::abs((transportJacobian)(3,3)) > 1e-50 ? diffMatrix(3,3)/((transportJacobian)(3,3)) : 0.;
730  m_thetaqop[recStep] = std::abs((transportJacobian)(3,4)) > 1e-50 ? diffMatrix(3,4)/((transportJacobian)(3,4)) : 0.;
731  m_thetasteps[recStep] = 4;
732  // fill the differences into the last one log (qop)
733  m_qoploc1[recStep] = std::abs((transportJacobian)(4,0)) > 1e-50 ? diffMatrix(4,0)/((transportJacobian)(4,0)) : 0.;
734  m_qoploc2[recStep] = std::abs((transportJacobian)(4,1)) > 1e-50 ? diffMatrix(4,1)/((transportJacobian)(4,1)) : 0.;
735  m_qopphi[recStep] = std::abs((transportJacobian)(4,2)) > 1e-50 ? diffMatrix(4,2)/((transportJacobian)(4,2)) : 0.;
736  m_qoptheta[recStep] = std::abs((transportJacobian)(4,3)) > 1e-50 ? diffMatrix(4,3)/((transportJacobian)(4,3)) : 0.;
737  m_qopqop[recStep] = std::abs((transportJacobian)(4,4)) > 1e-50 ? diffMatrix(4,4)/((transportJacobian)(4,4)) : 0.;
738  m_qopsteps[recStep] = 4;
739  ++recStep;
740 
741  // (R2)
742  // Relative differences ----------------------------------------------------
743  m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep-1] )) : 0.;
744  m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep-1] )) : 0.;
745  m_loc1phi[recStep] = std::abs(m_loc1phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep-1] )) : 0.;
746  m_loc1theta[recStep] = std::abs(m_loc1theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep-1])) : 0.;
747  m_loc1qop[recStep] = std::abs(m_loc1qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep-1] )) : 0.;
748  m_loc1steps[recStep] = 5;
749  // fill the differences into the last one log (loc2)
750  m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep-1] )) : 0.;
751  m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep-1] )) : 0.;
752  m_loc2phi[recStep] = std::abs(m_loc2phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep-1] )) : 0.;
753  m_loc2theta[recStep] = std::abs(m_loc2theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep-1])) : 0.;
754  m_loc2qop[recStep] = std::abs(m_loc2qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep-1] )) : 0.;
755  m_loc2steps[recStep] = 5;
756  // fill the differences into the last one log (phi)
757  m_philoc1[recStep] = std::abs(m_philoc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc1[recStep-1] )) : 0.;
758  m_philoc2[recStep] = std::abs(m_philoc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc2[recStep-1] )) : 0.;
759  m_phiphi[recStep] = std::abs(m_phiphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiphi[recStep-1] )) : 0.;
760  m_phitheta[recStep] = std::abs(m_phitheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phitheta[recStep-1])) : 0.;
761  m_phiqop[recStep] = std::abs(m_phiqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiqop[recStep-1] )) : 0.;
762  m_phisteps[recStep] = 5;
763  // fill the differences into the last one log (theta)
764  m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep-1] )) : 0.;
765  m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep-1] )) : 0.;
766  m_thetaphi[recStep] = std::abs(m_thetaphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep-1] )) : 0.;
767  m_thetatheta[recStep] = std::abs(m_thetatheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep-1])) : 0.;
768  m_thetaqop[recStep] = std::abs(m_thetaqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep-1] )) : 0.;
769  m_thetasteps[recStep] = 5;
770  // fill the differences into the last one log (qop)
771  m_qoploc1[recStep] = std::abs(m_qoploc1[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoploc1[recStep-1] ) ): 0.;
772  m_qoploc2[recStep] = std::abs(m_qoploc2[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoploc2[recStep-1] ) ): 0.;
773  m_qopphi[recStep] = std::abs(m_qopphi[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qopphi[recStep-1] ) ): 0.;
774  m_qoptheta[recStep] = std::abs(m_qoptheta[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoptheta[recStep-1]) ): 0.;
775  m_qopqop[recStep] = std::abs(m_qopqop[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qopqop[recStep-1] ) ): 0.;
776  m_qopsteps[recStep] = 5;
777  ++recStep;
778 
779  m_steps = recStep;
780  m_validationTree->Fill();
781  }
782 
783  // Code entered here will be executed once per event
784  return StatusCode::SUCCESS;
785 }
786 
787 //============================================================================================
788 Amg::Transform3D
789 Trk::RiddersAlgorithm::createTransform(double x, double y, double z, double phi, double theta, double alphaZ)
790 {
791 
792  if (phi!=0. && theta != 0.){
793  // create the Start Surface
794  Amg::Vector3D surfacePosition(x,y,z);
795  // z direction
796  Amg::Vector3D surfaceZdirection(cos(phi)*sin(theta),
797  sin(phi)*sin(theta),
798  cos(theta));
799  // the global z axis
800  Amg::Vector3D zAxis(0.,0.,1.);
801  // the y direction
802  Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
803  // the x direction
804  Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
805  // the rotation
806  Amg::RotationMatrix3D surfaceRotation;
807  surfaceRotation.col(0) = surfaceXdirection;
808  surfaceRotation.col(1) = surfaceYdirection;
809  surfaceRotation.col(2) = surfaceZdirection;
810  Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
811  return Amg::Transform3D(nominalTransform*Amg::AngleAxis3D(alphaZ,zAxis));
812 
813  }
814 
815  return Amg::Transform3D(Amg::Translation3D(x,y,z));
816 }
817 
818 double Trk::RiddersAlgorithm::parabolicInterpolation(double y0, double y1, double y2,
819  double x0, double x1, double x2)
820 {
821  double Z0 = x1*x2*y0;
822  double Z1 = x0*x2*y1;
823  double Z2 = x0*x1*y2;
824  double N0 = (x0-x1)*(x0-x2);
825  double N1 = (x1-x2)*(x1-x0);
826  double N2 = (x2-x0)*(x2-x1);
827  return Z0/N0 + Z1/N1 + Z2/N2;
828 }
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::y
@ y
Definition: ParamDefs.h:56
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
Trk::TransportJacobian
Definition: TransportJacobian.h:46
Trk::AmgMatrix
AmgMatrix(3, 3) NeutralParticleParameterCalculator
Definition: NeutralParticleParameterCalculator.cxx:233
dqt_zlumi_pandas.N2
N2
Definition: dqt_zlumi_pandas.py:325
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:30
TrackParameters.h
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:57
Trk::MagneticFieldProperties
Definition: MagneticFieldProperties.h:31
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
PlotCalibFromCool.zAxis
zAxis
Definition: PlotCalibFromCool.py:76
Trk::RiddersAlgorithm::parabolicInterpolation
static double parabolicInterpolation(double y0, double y1, double y2, double x0, double x1, double x2)
Langrange-parabolic interpolation.
Definition: RiddersAlgorithm.cxx:818
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
Trk::RiddersAlgorithm::RiddersAlgorithm
RiddersAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Standard Athena-Algorithm Constructor.
Definition: RiddersAlgorithm.cxx:24
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::loc2
@ loc2
generic first and second local coordinate
Definition: ParamDefs.h:35
Trk::alongMomentum
@ alongMomentum
Definition: PropDirection.h:20
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
TransportJacobian.h
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
RiddersAlgorithm.h
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
Trk::Z0
@ Z0
Definition: ParameterType.h:18
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
MagneticFieldProperties.h
CylinderVolumeBounds.h
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::CylinderSurface
Definition: CylinderSurface.h:55
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
CylinderSurface.h
dqt_zlumi_pandas.N1
int N1
Definition: dqt_zlumi_pandas.py:322
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::RiddersAlgorithm::~RiddersAlgorithm
~RiddersAlgorithm()
Default Destructor.
Definition: RiddersAlgorithm.cxx:29
AthAlgorithm
Definition: AthAlgorithm.h:47
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
charge
double charge(const T &p)
Definition: AtlasPID.h:931
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::RiddersAlgorithm::initialize
StatusCode initialize()
standard Athena-Algorithm method
Definition: RiddersAlgorithm.cxx:39
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:789
Trk::PlaneSurface
Definition: PlaneSurface.h:64
Amg::RotationMatrix3D
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
Definition: GeoPrimitives.h:49
PlaneSurface.h
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition: GeoPrimitives.h:44
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
Amg::AngleAxis3D
Eigen::AngleAxisd AngleAxis3D
Definition: GeoPrimitives.h:45
Trk::phi
@ phi
Definition: ParamDefs.h:75
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Trk::x
@ x
Definition: ParamDefs.h:55
Trk::loc1
@ loc1
Definition: ParamDefs.h:34
Trk::RiddersAlgorithm::execute
StatusCode execute()
standard Athena-Algorithm method
Definition: RiddersAlgorithm.cxx:137
Trk::RiddersAlgorithm::finalize
StatusCode finalize()
standard Athena-Algorithm method
Definition: RiddersAlgorithm.cxx:129
Generated on Mon Mar 31 2025 21:17:32 for ATLAS Offline Software by   doxygen 1.8.18