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