ATLAS Offline Software
RiddersAlgorithm.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 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  SmartIF<ITHistSvc> tHistSvc{service("THistSvc")};
206  if (!tHistSvc){
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  unsigned int recStep = 0;
386  const auto& transportJacobian = (*optTransportJacobian);
387  // [0] Transport Jacobian -----------------------------------------------------
388  ATH_MSG_VERBOSE( "TransportJacobian : " << transportJacobian );
389 
390 
391  // and now fill the variables
392  m_loc1loc1[recStep] = (transportJacobian)(0,0);
393  m_loc1loc2[recStep] = (transportJacobian)(0,1);
394  m_loc1phi[recStep] = (transportJacobian)(0,2);
395  m_loc1theta[recStep] = (transportJacobian)(0,3);
396  m_loc1qop[recStep] = (transportJacobian)(0,4);
397  m_loc1steps[recStep] = 0.;
398 
399  m_loc2loc1[recStep] = (transportJacobian)(1,0);
400  m_loc2loc2[recStep] = (transportJacobian)(1,1);
401  m_loc2phi[recStep] = (transportJacobian)(1,2);
402  m_loc2theta[recStep] = (transportJacobian)(1,3);
403  m_loc2qop[recStep] = (transportJacobian)(1,4);
404  m_loc2steps[recStep] = 0.;
405 
406  m_philoc1[recStep] = (transportJacobian)(2,0);
407  m_philoc2[recStep] = (transportJacobian)(2,1);
408  m_phiphi[recStep] = (transportJacobian)(2,2);
409  m_phitheta[recStep] = (transportJacobian)(2,3);
410  m_phiqop[recStep] = (transportJacobian)(2,4);
411  m_phisteps[recStep] = 0.;
412 
413  m_thetaloc1[recStep] = (transportJacobian)(3,0);
414  m_thetaloc2[recStep] = (transportJacobian)(3,1);
415  m_thetaphi[recStep] = (transportJacobian)(3,2);
416  m_thetatheta[recStep] = (transportJacobian)(3,3);
417  m_thetaqop[recStep] = (transportJacobian)(3,4);
418  m_thetasteps[recStep] = 0.;
419 
420  m_qoploc1[recStep] = (transportJacobian)(4,0);
421  m_qoploc2[recStep] = (transportJacobian)(4,1);
422  m_qopphi[recStep] = (transportJacobian)(4,2);
423  m_qoptheta[recStep] = (transportJacobian)(4,3);
424  m_qopqop[recStep] = (transportJacobian)(4,4);
425  m_qopsteps[recStep] = 0.;
426 
427  ++recStep;
428 
429  // start the riddlers algorithm
430  // [1-2-3] Riddlers jacobians -----------------------------------------
431  for (unsigned int istep = 0; istep < m_localVariations.size(); ++istep){
432  // --------------------------
433  ATH_MSG_VERBOSE( "Performing step : " << istep );
434  // the initial perigee with random numbers
435  // for the first row
436  Trk::AtaPlane startLoc1Minus(loc1-m_localVariations[istep],loc2,phi,theta,qOverP,startSurface);
437  Trk::AtaPlane startLoc1Plus(loc1+m_localVariations[istep],loc2,phi,theta,qOverP,startSurface);
438  // for the second row
439  Trk::AtaPlane startLoc2Minus(loc1,loc2-m_localVariations[istep],phi,theta,qOverP,startSurface);
440  Trk::AtaPlane startLoc2Plus(loc1,loc2+m_localVariations[istep],phi,theta,qOverP,startSurface);
441  // for the third row
442  Trk::AtaPlane startPhiMinus(loc1,loc2,phi-m_angularVariations[istep],theta,qOverP,startSurface);
443  Trk::AtaPlane startPhiPlus(loc1,loc2,phi+m_angularVariations[istep],theta,qOverP,startSurface);
444  // for the fourth row
445  Trk::AtaPlane startThetaMinus(loc1,loc2,phi,theta-m_angularVariations[istep],qOverP,startSurface);
446  Trk::AtaPlane startThetaPlus(loc1,loc2,phi,theta+m_angularVariations[istep],qOverP,startSurface);
447  // for the fifth row
448  Trk::AtaPlane startQopMinus(loc1,loc2,phi,theta,qOverP-m_qOpVariations[istep],startSurface);
449  Trk::AtaPlane startQopPlus(loc1,loc2,phi,theta,qOverP+m_qOpVariations[istep],startSurface);
450 
451  // the propagations --- 10 times
452  auto endLoc1Minus = m_propagator->propagateParameters(ctx,
453  startLoc1Minus,
454  destinationSurface,
455  Trk::alongMomentum,
456  false,
457  *m_magFieldProperties);
458 
459 
460  auto endLoc1Plus = m_propagator->propagateParameters(ctx,
461  startLoc1Plus,
462  destinationSurface,
463  Trk::alongMomentum,
464  false,
465  *m_magFieldProperties);
466 
467  auto endLoc2Minus = m_propagator->propagateParameters(ctx,
468  startLoc2Minus,
469  destinationSurface,
470  Trk::alongMomentum,
471  false,
472  *m_magFieldProperties);
473 
474  auto endLoc2Plus = m_propagator->propagateParameters(ctx,
475  startLoc2Plus,
476  destinationSurface,
477  Trk::alongMomentum,
478  false,
479  *m_magFieldProperties);
480 
481  auto endPhiMinus = m_propagator->propagateParameters(ctx,
482  startPhiMinus,
483  destinationSurface,
484  Trk::alongMomentum,
485  false,
486  *m_magFieldProperties);
487 
488  auto endPhiPlus = m_propagator->propagateParameters(ctx,
489  startPhiPlus,
490  destinationSurface,
491  Trk::alongMomentum,
492  false,
493  *m_magFieldProperties);
494 
495  auto endThetaMinus = m_propagator->propagateParameters(ctx,
496  startThetaMinus,
497  destinationSurface,
498  Trk::alongMomentum,
499  false,
500  *m_magFieldProperties);
501 
502  auto endThetaPlus = m_propagator->propagateParameters(ctx,
503  startThetaPlus,
504  destinationSurface,
505  Trk::alongMomentum,
506  false,
507  *m_magFieldProperties);
508 
509  auto endQopMinus = m_propagator->propagateParameters(ctx,
510  startQopMinus,
511  destinationSurface,
512  Trk::alongMomentum,
513  false,
514  *m_magFieldProperties);
515 
516  auto endQopPlus = m_propagator->propagateParameters(ctx,
517  startQopPlus,
518  destinationSurface,
519  Trk::alongMomentum,
520  false,
521  *m_magFieldProperties);
522  if (endLoc1Minus
523  && endLoc1Plus
524  && endLoc2Minus
525  && endLoc2Plus
526  && endPhiMinus
527  && endPhiPlus
528  && endThetaMinus
529  && endThetaPlus
530  && endQopMinus
531  && endQopPlus){
532 
533  // Loc 1
534  const Amg::VectorX& endLoc1MinusPar = endLoc1Minus->parameters();
535  const Amg::VectorX& endLoc1PlusPar = endLoc1Plus->parameters();
536  // Loc 2
537  const Amg::VectorX& endLoc2MinusPar = endLoc2Minus->parameters();
538  const Amg::VectorX& endLoc2PlusPar = endLoc2Plus->parameters();
539  // Phi
540  const Amg::VectorX& endPhiMinusPar = endPhiMinus->parameters();
541  const Amg::VectorX& endPhiPlusPar = endPhiPlus->parameters();
542  // Theta
543  const Amg::VectorX& endThetaMinusPar = endThetaMinus->parameters();
544  const Amg::VectorX& endThetaPlusPar = endThetaPlus->parameters();
545  // qop
546  const Amg::VectorX& endQopMinusPar = endQopMinus->parameters();
547  const Amg::VectorX& endQopPlusPar = endQopPlus->parameters();
548 
549  // the deltas
550  Amg::VectorX endLoc1Diff(endLoc1PlusPar-endLoc1MinusPar);
551  Amg::VectorX endLoc2Diff(endLoc2PlusPar-endLoc2MinusPar);
552  Amg::VectorX endPhiDiff(endPhiPlusPar-endPhiMinusPar);
553  Amg::VectorX endThetaDiff(endThetaPlusPar-endThetaMinusPar);
554  Amg::VectorX endQopDiff(endQopPlusPar-endQopMinusPar);
555 
556  currentStepJacobian(0,0) = endLoc1Diff[0]/(2.*m_localVariations[istep]);
557  currentStepJacobian(0,1) = endLoc2Diff[0]/(2.*m_localVariations[istep]);
558  currentStepJacobian(0,2) = endPhiDiff[0]/(2.*m_angularVariations[istep]);
559  currentStepJacobian(0,3) = endThetaDiff[0]/(2.*m_angularVariations[istep]);
560  currentStepJacobian(0,4) = endQopDiff[0]/(2.*m_qOpVariations[istep]);
561 
562  m_loc1loc1[recStep] = currentStepJacobian(0,0);
563  m_loc1loc2[recStep] = currentStepJacobian(0,1);
564  m_loc1phi[recStep] = currentStepJacobian(0,2);
565  m_loc1theta[recStep] = currentStepJacobian(0,3);
566  m_loc1qop[recStep] = currentStepJacobian(0,4);
567  m_loc1steps[recStep] = m_localVariations[istep];
568 
569  currentStepJacobian(1,0) = endLoc1Diff[1]/(2.*m_localVariations[istep]);
570  currentStepJacobian(1,1) = endLoc2Diff[1]/(2.*m_localVariations[istep]);
571  currentStepJacobian(1,2) = endPhiDiff[1]/(2.*m_angularVariations[istep]);
572  currentStepJacobian(1,3) = endThetaDiff[1]/(2.*m_angularVariations[istep]);
573  currentStepJacobian(1,4) = endQopDiff[1]/(2.*m_qOpVariations[istep]);
574 
575  m_loc2loc1[recStep] = currentStepJacobian(1,0);
576  m_loc2loc2[recStep] = currentStepJacobian(1,1);
577  m_loc2phi[recStep] = currentStepJacobian(1,2);
578  m_loc2theta[recStep] = currentStepJacobian(1,3);
579  m_loc2qop[recStep] = currentStepJacobian(1,4);
580  m_loc2steps[recStep] = m_localVariations[istep];
581 
582  currentStepJacobian(2,0) = endLoc1Diff[2]/(2.*m_localVariations[istep]);
583  currentStepJacobian(2,1) = endLoc2Diff[2]/(2.*m_localVariations[istep]);
584  currentStepJacobian(2,2) = endPhiDiff[2]/(2.*m_angularVariations[istep]);
585  currentStepJacobian(2,3) = endThetaDiff[2]/(2.*m_angularVariations[istep]);
586  currentStepJacobian(2,4) = endQopDiff[2]/(2.*m_qOpVariations[istep]);
587 
588  m_philoc1[recStep] = currentStepJacobian(2,0);
589  m_philoc2[recStep] = currentStepJacobian(2,1);
590  m_phiphi[recStep] = currentStepJacobian(2,2);
591  m_phitheta[recStep] = currentStepJacobian(2,3);
592  m_phiqop[recStep] = currentStepJacobian(2,4);
593  m_phisteps[recStep] = m_angularVariations[istep];
594 
595  currentStepJacobian(3,0) = endLoc1Diff[3]/(2.*m_localVariations[istep]);
596  currentStepJacobian(3,1) = endLoc2Diff[3]/(2.*m_localVariations[istep]);
597  currentStepJacobian(3,2) = endPhiDiff[3]/(2.*m_angularVariations[istep]);
598  currentStepJacobian(3,3) = endThetaDiff[3]/(2.*m_angularVariations[istep]);
599  currentStepJacobian(3,4) = endQopDiff[3]/(2.*m_qOpVariations[istep]);
600 
601  m_thetaloc1[recStep] = currentStepJacobian(3,0);
602  m_thetaloc2[recStep] = currentStepJacobian(3,1);
603  m_thetaphi[recStep] = currentStepJacobian(3,2);
604  m_thetatheta[recStep] = currentStepJacobian(3,3);
605  m_thetaqop[recStep] = currentStepJacobian(3,4);
606  m_thetasteps[recStep] = m_angularVariations[istep];
607 
608  currentStepJacobian(4,0) = endLoc1Diff[4]/(2.*m_localVariations[istep]);
609  currentStepJacobian(4,1) = endLoc2Diff[4]/(2.*m_localVariations[istep]);
610  currentStepJacobian(4,2) = endPhiDiff[4]/(2.*m_angularVariations[istep]);
611  currentStepJacobian(4,3) = endThetaDiff[4]/(2.*m_angularVariations[istep]);
612  currentStepJacobian(4,4) = endQopDiff[4]/(2.*m_qOpVariations[istep]);
613 
614  m_qoploc1[recStep] = currentStepJacobian(4,0);
615  m_qoploc2[recStep] = currentStepJacobian(4,1);
616  m_qopphi[recStep] = currentStepJacobian(4,2);
617  m_qoptheta[recStep] = currentStepJacobian(4,3);
618  m_qopqop[recStep] = currentStepJacobian(4,4);
619  m_qopsteps[recStep] = m_qOpVariations[istep];
620 
621  ATH_MSG_DEBUG( "Current TransportJacobian : " << currentStepJacobian );
622 
623  ++recStep;
624  }
625  }
626 
627  // -------------------------------------------------------------------------------
628  if (recStep > 2){
629  // The parabolic interpolation -------------------------------------------------------------------------------
630  // dL1
631  m_loc1loc1[recStep] = parabolicInterpolation(m_loc1loc1[recStep-1],m_loc1loc1[recStep-2],m_loc1loc1[recStep-3],
632  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
633  m_loc1loc2[recStep] = parabolicInterpolation(m_loc1loc2[recStep-1],m_loc1loc2[recStep-2],m_loc1loc2[recStep-3],
634  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
635  m_loc1phi[recStep] = parabolicInterpolation(m_loc1phi[recStep-1],m_loc1phi[recStep-2],m_loc1phi[recStep-3],
636  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
637  m_loc1theta[recStep] = parabolicInterpolation(m_loc1theta[recStep-1],m_loc1theta[recStep-2],m_loc1theta[recStep-3],
638  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
639  m_loc1qop[recStep] = parabolicInterpolation(m_loc1qop[recStep-1],m_loc1qop[recStep-2],m_loc1qop[recStep-3],
640  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
641  m_loc1steps[recStep] = 1;
642  // dL2
643  m_loc2loc1[recStep] = parabolicInterpolation(m_loc2loc1[recStep-1],m_loc2loc1[recStep-2],m_loc2loc1[recStep-3],
644  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
645  m_loc2loc2[recStep] = parabolicInterpolation(m_loc2loc2[recStep-1],m_loc2loc2[recStep-2],m_loc2loc2[recStep-3],
646  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
647  m_loc2phi[recStep] = parabolicInterpolation(m_loc2phi[recStep-1],m_loc2phi[recStep-2],m_loc2phi[recStep-3],
648  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
649  m_loc2theta[recStep] = parabolicInterpolation(m_loc2theta[recStep-1],m_loc2theta[recStep-2],m_loc2theta[recStep-3],
650  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
651  m_loc2qop[recStep] = parabolicInterpolation(m_loc2qop[recStep-1],m_loc2qop[recStep-2],m_loc2qop[recStep-3],
652  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
653  m_loc2steps[recStep] = 1;
654  // dPhi
655  m_philoc1[recStep] = parabolicInterpolation(m_philoc1[recStep-1],m_philoc1[recStep-2],m_philoc1[recStep-3],
656  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
657  m_philoc2[recStep] = parabolicInterpolation(m_philoc2[recStep-1],m_philoc2[recStep-2],m_philoc2[recStep-3],
658  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
659  m_phiphi[recStep] = parabolicInterpolation(m_phiphi[recStep-1],m_phiphi[recStep-2],m_phiphi[recStep-3],
660  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
661  m_phitheta[recStep] = parabolicInterpolation(m_phitheta[recStep-1],m_phitheta[recStep-2],m_phitheta[recStep-3],
662  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
663  m_phiqop[recStep] = parabolicInterpolation(m_phiqop[recStep-1],m_phiqop[recStep-2],m_phiqop[recStep-3],
664  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
665  m_phisteps[recStep] = 1;
666  // dTheta
667  m_thetaloc1[recStep] = parabolicInterpolation(m_thetaloc1[recStep-1],m_thetaloc1[recStep-2],m_thetaloc1[recStep-3],
668  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
669  m_thetaloc2[recStep] = parabolicInterpolation(m_thetaloc2[recStep-1],m_thetaloc2[recStep-2],m_thetaloc2[recStep-3],
670  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
671  m_thetaphi[recStep] = parabolicInterpolation(m_thetaphi[recStep-1],m_thetaphi[recStep-2],m_thetaphi[recStep-3],
672  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
673  m_thetatheta[recStep] = parabolicInterpolation(m_thetatheta[recStep-1],m_thetatheta[recStep-2],m_thetatheta[recStep-3],
674  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
675  m_thetaqop[recStep] = parabolicInterpolation(m_thetaqop[recStep-1],m_thetaqop[recStep-2],m_thetaqop[recStep-3],
676  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
677  m_thetasteps[recStep] = 1;
678  // dTheta
679  m_qoploc1[recStep] = parabolicInterpolation(m_qoploc1[recStep-1],m_qoploc1[recStep-2],m_qoploc1[recStep-3],
680  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
681  m_qoploc2[recStep] = parabolicInterpolation(m_qoploc2[recStep-1],m_qoploc2[recStep-2],m_qoploc2[recStep-3],
682  m_localVariations[2], m_localVariations[1], m_localVariations[0]);
683  m_qopphi[recStep] = parabolicInterpolation(m_qopphi[recStep-1],m_qopphi[recStep-2],m_qopphi[recStep-3],
684  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
685  m_qoptheta[recStep] = parabolicInterpolation(m_qoptheta[recStep-1],m_qoptheta[recStep-2],m_qoptheta[recStep-3],
686  m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
687  m_qopqop[recStep] = parabolicInterpolation(m_qopqop[recStep-1],m_qopqop[recStep-2],m_qopqop[recStep-3],
688  m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
689  m_qopsteps[recStep] = 1;
690 
691  currentStepJacobian(0,0)= m_loc1loc1[recStep];
692  currentStepJacobian(0,1)= m_loc1loc2[recStep];
693  currentStepJacobian(0,2)= m_loc1phi[recStep];
694  currentStepJacobian(0,3)= m_loc1theta[recStep];
695  currentStepJacobian(0,4)= m_loc1qop[recStep];
696 
697  currentStepJacobian(1,0)= m_loc2loc1[recStep];
698  currentStepJacobian(1,1)= m_loc2loc2[recStep];
699  currentStepJacobian(1,2)= m_loc2phi[recStep];
700  currentStepJacobian(1,3)= m_loc2theta[recStep];
701  currentStepJacobian(1,4)= m_loc2qop[recStep];
702 
703  currentStepJacobian(2,0)= m_philoc1[recStep];
704  currentStepJacobian(2,1)= m_philoc2[recStep];
705  currentStepJacobian(2,2)= m_phiphi[recStep];
706  currentStepJacobian(2,3)= m_phitheta[recStep];
707  currentStepJacobian(2,4)= m_phiqop[recStep];
708 
709  currentStepJacobian(3,0)= m_thetaloc1[recStep];
710  currentStepJacobian(3,1)= m_thetaloc2[recStep];
711  currentStepJacobian(3,2)= m_thetaphi[recStep];
712  currentStepJacobian(3,3)= m_thetatheta[recStep];
713  currentStepJacobian(3,4)= m_thetaqop[recStep];
714 
715  currentStepJacobian(4,0)= m_qoploc1[recStep];
716  currentStepJacobian(4,1)= m_qoploc2[recStep];
717  currentStepJacobian(4,2)= m_qopphi[recStep];
718  currentStepJacobian(4,3)= m_qoptheta[recStep];
719  currentStepJacobian(4,4)= m_qopqop[recStep];
720 
721  }
722 
723  ATH_MSG_DEBUG( "Interpolated TransportJacobian : " << currentStepJacobian );
724  ++recStep;
725 
726 
727  // now fill the results
728  TransportJacobian diffMatrix(transportJacobian-currentStepJacobian);
729 
730  ATH_MSG_VERBOSE( "Absolute Differences of the TransportJacobian : " << diffMatrix );
731 
732  // Absolute differences ----------------------------------------------------
733  // (A1)
734  // fill the differences into the last one log (loc1)
735  m_loc1loc1[recStep] = diffMatrix(0,0);
736  m_loc1loc2[recStep] = diffMatrix(0,1);
737  m_loc1phi[recStep] = diffMatrix(0,2);
738  m_loc1theta[recStep] = diffMatrix(0,3);
739  m_loc1qop[recStep] = diffMatrix(0,4);
740  m_loc1steps[recStep] = 2;
741  // fill the differences into the last one log (loc2)
742  m_loc2loc1[recStep] = diffMatrix(1,0);
743  m_loc2loc2[recStep] = diffMatrix(1,1);
744  m_loc2phi[recStep] = diffMatrix(1,2);
745  m_loc2theta[recStep] = diffMatrix(1,3);
746  m_loc2qop[recStep] = diffMatrix(1,4);
747  m_loc2steps[recStep] = 2;
748  // fill the differences into the last one log (phi)
749  m_philoc1[recStep] = diffMatrix(2,0);
750  m_philoc2[recStep] = diffMatrix(2,1);
751  m_phiphi[recStep] = diffMatrix(2,2);
752  m_phitheta[recStep] = diffMatrix(2,3);
753  m_phiqop[recStep] = diffMatrix(2,4);
754  m_phisteps[recStep] = 2;
755  // fill the differences into the last one log (theta)
756  m_thetaloc1[recStep] = diffMatrix(3,0);
757  m_thetaloc2[recStep] = diffMatrix(3,1);
758  m_thetaphi[recStep] = diffMatrix(3,2);
759  m_thetatheta[recStep] = diffMatrix(3,3);
760  m_thetaqop[recStep] = diffMatrix(3,4);
761  m_thetasteps[recStep] = 2;
762  // fill the differences into the last one log (qop)
763  m_qoploc1[recStep] = diffMatrix(4,0);
764  m_qoploc2[recStep] = diffMatrix(4,1);
765  m_qopphi[recStep] = diffMatrix(4,2);
766  m_qoptheta[recStep] = diffMatrix(4,3);
767  m_qopqop[recStep] = diffMatrix(4,4);
768  m_qopsteps[recStep] = 2;
769  ++recStep;
770 
771  // (A2)
772  if (recStep > 1){
773  // fill the differences into the last one log (loc1)
774  m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep-1] )) : 0.;
775  m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep-1] )) : 0.;
776  m_loc1phi[recStep] = std::abs(m_loc1phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep-1] )) : 0.;
777  m_loc1theta[recStep] = std::abs(m_loc1theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep-1])) : 0.;
778  m_loc1qop[recStep] = std::abs(m_loc1qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep-1] )) : 0.;
779  m_loc1steps[recStep] = 3;
780  // fill the differences into the last one log (loc2)
781  m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep-1] )) : 0.;
782  m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep-1] )) : 0.;
783  m_loc2phi[recStep] = std::abs(m_loc2phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep-1] )) : 0.;
784  m_loc2theta[recStep] = std::abs(m_loc2theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep-1])) : 0.;
785  m_loc2qop[recStep] = std::abs(m_loc2qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep-1] )) : 0.;
786  m_loc2steps[recStep] = 3;
787  // fill the differences into the last one log (phi)
788  m_philoc1[recStep] = std::abs(m_philoc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc1[recStep-1] )) : 0.;
789  m_philoc2[recStep] = std::abs(m_philoc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc2[recStep-1] )) : 0.;
790  m_phiphi[recStep] = std::abs(m_phiphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiphi[recStep-1] )) : 0.;
791  m_phitheta[recStep] = std::abs(m_phitheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phitheta[recStep-1])) : 0.;
792  m_phiqop[recStep] = std::abs(m_phiqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiqop[recStep-1] )) : 0.;
793  m_phisteps[recStep] = 3;
794  // fill the differences into the last one log (theta)
795  m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep-1] )) : 0.;
796  m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep-1] )) : 0.;
797  m_thetaphi[recStep] = std::abs(m_thetaphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep-1] )) : 0.;
798  m_thetatheta[recStep] = std::abs(m_thetatheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep-1])) : 0.;
799  m_thetaqop[recStep] = std::abs(m_thetaqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep-1] )) : 0.;
800  m_thetasteps[recStep] = 3;
801  // fill the differences into the last one log (qop)
802  m_qoploc1[recStep] = std::abs(m_qoploc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoploc1[recStep-1] )) : 0.;
803  m_qoploc2[recStep] = std::abs(m_qoploc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoploc2[recStep-1] )) : 0.;
804  m_qopphi[recStep] = std::abs(m_qopphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qopphi[recStep-1] )) : 0.;
805  m_qoptheta[recStep] = std::abs(m_qoptheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoptheta[recStep-1])) : 0.;
806  m_qopqop[recStep] = std::abs(m_qopqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qopqop[recStep-1] )) : 0.;
807  m_qopsteps[recStep] = 3;
808  }
809  ++recStep;
810 
811 
812  // Relative differences) ----------------------------------------------------
813  // (R1)
814  // fill the differences into the last one log (loc1)
815  m_loc1loc1[recStep] = std::abs((transportJacobian)(0,0)) > 1e-50 ? diffMatrix(0,0)/((transportJacobian)(0,0)) : 0.;
816  m_loc1loc2[recStep] = std::abs((transportJacobian)(0,1)) > 1e-50 ? diffMatrix(0,1)/((transportJacobian)(0,1)) : 0.;
817  m_loc1phi[recStep] = std::abs((transportJacobian)(0,2)) > 1e-50 ? diffMatrix(0,2)/((transportJacobian)(0,2)) : 0.;
818  m_loc1theta[recStep] = std::abs((transportJacobian)(0,3)) > 1e-50 ? diffMatrix(0,3)/((transportJacobian)(0,3)) : 0.;
819  m_loc1qop[recStep] = std::abs((transportJacobian)(0,4)) > 1e-50 ? diffMatrix(0,4)/((transportJacobian)(0,4)) : 0.;
820  m_loc1steps[recStep] = 4;
821  // fill the differences into the last one log (loc2)
822  m_loc2loc1[recStep] = std::abs((transportJacobian)(1,0)) > 1e-50 ? diffMatrix(1,0)/((transportJacobian)(1,0)) : 0.;
823  m_loc2loc2[recStep] = std::abs((transportJacobian)(1,1)) > 1e-50 ? diffMatrix(1,1)/((transportJacobian)(1,1)) : 0.;
824  m_loc2phi[recStep] = std::abs((transportJacobian)(1,2)) > 1e-50 ? diffMatrix(1,2)/((transportJacobian)(1,2)) : 0.;
825  m_loc2theta[recStep] = std::abs((transportJacobian)(1,3)) > 1e-50 ? diffMatrix(1,3)/((transportJacobian)(1,3)) : 0.;
826  m_loc2qop[recStep] = std::abs((transportJacobian)(1,4)) > 1e-50 ? diffMatrix(1,4)/((transportJacobian)(1,4)) : 0.;
827  m_loc2steps[recStep] = 4;
828  // fill the differences into the last one log (phi)
829  m_philoc1[recStep] = std::abs((transportJacobian)(2,0)) > 1e-50 ? diffMatrix(2,0)/((transportJacobian)(2,0)) : 0.;
830  m_philoc2[recStep] = std::abs((transportJacobian)(2,1)) > 1e-50 ? diffMatrix(2,1)/((transportJacobian)(2,1)) : 0.;
831  m_phiphi[recStep] = std::abs((transportJacobian)(2,2)) > 1e-50 ? diffMatrix(2,2)/((transportJacobian)(2,2)) : 0.;
832  m_phitheta[recStep] = std::abs((transportJacobian)(2,3)) > 1e-50 ? diffMatrix(2,3)/((transportJacobian)(2,3)) : 0.;
833  m_phiqop[recStep] = std::abs((transportJacobian)(2,4)) > 1e-50 ? diffMatrix(2,4)/((transportJacobian)(2,4)) : 0.;
834  m_phisteps[recStep] = 4;
835  // fill the differences into the last one log (theta)
836  m_thetaloc1[recStep] = std::abs((transportJacobian)(3,0)) > 1e-50 ? diffMatrix(3,0)/((transportJacobian)(3,0)) : 0.;
837  m_thetaloc2[recStep] = std::abs((transportJacobian)(3,1)) > 1e-50 ? diffMatrix(3,1)/((transportJacobian)(3,1)) : 0.;
838  m_thetaphi[recStep] = std::abs((transportJacobian)(3,2)) > 1e-50 ? diffMatrix(3,2)/((transportJacobian)(3,2)) : 0.;
839  m_thetatheta[recStep] = std::abs((transportJacobian)(3,3)) > 1e-50 ? diffMatrix(3,3)/((transportJacobian)(3,3)) : 0.;
840  m_thetaqop[recStep] = std::abs((transportJacobian)(3,4)) > 1e-50 ? diffMatrix(3,4)/((transportJacobian)(3,4)) : 0.;
841  m_thetasteps[recStep] = 4;
842  // fill the differences into the last one log (qop)
843  m_qoploc1[recStep] = std::abs((transportJacobian)(4,0)) > 1e-50 ? diffMatrix(4,0)/((transportJacobian)(4,0)) : 0.;
844  m_qoploc2[recStep] = std::abs((transportJacobian)(4,1)) > 1e-50 ? diffMatrix(4,1)/((transportJacobian)(4,1)) : 0.;
845  m_qopphi[recStep] = std::abs((transportJacobian)(4,2)) > 1e-50 ? diffMatrix(4,2)/((transportJacobian)(4,2)) : 0.;
846  m_qoptheta[recStep] = std::abs((transportJacobian)(4,3)) > 1e-50 ? diffMatrix(4,3)/((transportJacobian)(4,3)) : 0.;
847  m_qopqop[recStep] = std::abs((transportJacobian)(4,4)) > 1e-50 ? diffMatrix(4,4)/((transportJacobian)(4,4)) : 0.;
848  m_qopsteps[recStep] = 4;
849  ++recStep;
850 
851  // (R2)
852  // Relative differences ----------------------------------------------------
853  m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep-1] )) : 0.;
854  m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep-1] )) : 0.;
855  m_loc1phi[recStep] = std::abs(m_loc1phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep-1] )) : 0.;
856  m_loc1theta[recStep] = std::abs(m_loc1theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep-1])) : 0.;
857  m_loc1qop[recStep] = std::abs(m_loc1qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep-1] )) : 0.;
858  m_loc1steps[recStep] = 5;
859  // fill the differences into the last one log (loc2)
860  m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep-1] )) : 0.;
861  m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep-1] )) : 0.;
862  m_loc2phi[recStep] = std::abs(m_loc2phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep-1] )) : 0.;
863  m_loc2theta[recStep] = std::abs(m_loc2theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep-1])) : 0.;
864  m_loc2qop[recStep] = std::abs(m_loc2qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep-1] )) : 0.;
865  m_loc2steps[recStep] = 5;
866  // fill the differences into the last one log (phi)
867  m_philoc1[recStep] = std::abs(m_philoc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc1[recStep-1] )) : 0.;
868  m_philoc2[recStep] = std::abs(m_philoc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc2[recStep-1] )) : 0.;
869  m_phiphi[recStep] = std::abs(m_phiphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiphi[recStep-1] )) : 0.;
870  m_phitheta[recStep] = std::abs(m_phitheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phitheta[recStep-1])) : 0.;
871  m_phiqop[recStep] = std::abs(m_phiqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiqop[recStep-1] )) : 0.;
872  m_phisteps[recStep] = 5;
873  // fill the differences into the last one log (theta)
874  m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep-1] )) : 0.;
875  m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep-1] )) : 0.;
876  m_thetaphi[recStep] = std::abs(m_thetaphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep-1] )) : 0.;
877  m_thetatheta[recStep] = std::abs(m_thetatheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep-1])) : 0.;
878  m_thetaqop[recStep] = std::abs(m_thetaqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep-1] )) : 0.;
879  m_thetasteps[recStep] = 5;
880  // fill the differences into the last one log (qop)
881  m_qoploc1[recStep] = std::abs(m_qoploc1[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoploc1[recStep-1] ) ): 0.;
882  m_qoploc2[recStep] = std::abs(m_qoploc2[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoploc2[recStep-1] ) ): 0.;
883  m_qopphi[recStep] = std::abs(m_qopphi[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qopphi[recStep-1] ) ): 0.;
884  m_qoptheta[recStep] = std::abs(m_qoptheta[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoptheta[recStep-1]) ): 0.;
885  m_qopqop[recStep] = std::abs(m_qopqop[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qopqop[recStep-1] ) ): 0.;
886  m_qopsteps[recStep] = 5;
887  ++recStep;
888 
889  m_steps = recStep;
890  m_validationTree->Fill();
891  }
892 
893  // Code entered here will be executed once per event
894  return StatusCode::SUCCESS;
895 }
896 
897 //============================================================================================
898 Amg::Transform3D
899 Trk::RiddersAlgorithm::createTransform(double x, double y, double z, double phi, double theta, double alphaZ)
900 {
901 
902  if (phi!=0. && theta != 0.){
903  // create the Start Surface
904  Amg::Vector3D surfacePosition(x,y,z);
905  // z direction
906  Amg::Vector3D surfaceZdirection(cos(phi)*sin(theta),
907  sin(phi)*sin(theta),
908  cos(theta));
909  // the global z axis
910  Amg::Vector3D zAxis(0.,0.,1.);
911  // the y direction
912  Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
913  // the x direction
914  Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
915  // the rotation
916  Amg::RotationMatrix3D surfaceRotation;
917  surfaceRotation.col(0) = surfaceXdirection;
918  surfaceRotation.col(1) = surfaceYdirection;
919  surfaceRotation.col(2) = surfaceZdirection;
920  Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
921  return Amg::Transform3D(nominalTransform*Amg::AngleAxis3D(alphaZ,zAxis));
922 
923  }
924 
925  return Amg::Transform3D(Amg::Translation3D(x,y,z));
926 }
927 
928 double Trk::RiddersAlgorithm::parabolicInterpolation(double y0, double y1, double y2,
929  double x0, double x1, double x2)
930 {
931  double Z0 = x1*x2*y0;
932  double Z1 = x0*x2*y1;
933  double Z2 = x0*x1*y2;
934  double N0 = (x0-x1)*(x0-x2);
935  double N1 = (x1-x2)*(x1-x0);
936  double N2 = (x2-x0)*(x2-x1);
937  return Z0/N0 + Z1/N1 + Z2/N2;
938 }
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
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:17
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
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:928
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:35
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
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
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:228
python.SystemOfUnits.tesla
int tesla
Definition: SystemOfUnits.py:228
charge
double charge(const T &p)
Definition: AtlasPID.h:756
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
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:899
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
Gaudi
=============================================================================
Definition: CaloGPUClusterAndCellDataMonitorOptions.h:273
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:247
Trk::RiddersAlgorithm::finalize
StatusCode finalize()
standard Athena-Algorithm method
Definition: RiddersAlgorithm.cxx:239
Generated on Mon Dec 23 2024 21:17:16 for ATLAS Offline Software by   doxygen 1.8.18