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