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()
138{
139 const EventContext& ctx = Gaudi::Hive::currentContext();
140 // this is fine
141 double p = m_minP + m_flatDist->shoot()*(m_maxP-m_minP);
142 double charge = (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
143 double qOverP = charge/p;
144
145 // for the momentum logging
146 // m_startP = p;
147
148 // the local start start
149 double loc1 = m_sigmaLoc * m_gaussDist->shoot();
150 double loc2 = m_sigmaLoc * m_gaussDist->shoot();
151 // are adopted for planar and straight line surfaces
152 double phi = m_minPhi + m_maxPhi * m_flatDist->shoot();
153 double eta = m_minEta + m_flatDist->shoot()*(m_maxEta-m_minEta);
154 double theta = 2.*atan(exp(-eta));
155
156 // start
157 double startR = std::abs(m_sigmaR * m_gaussDist->shoot());
158 double surfacePhi = M_PI * m_flatDist->shoot();
159 surfacePhi *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
160 double startX = startR*cos(surfacePhi);
161 double startY = startR*sin(surfacePhi);
162 double startZ = m_sigmaLoc * m_gaussDist->shoot();
163
164 // rotate it around Z
165 double alphaZ = M_PI * m_flatDist->shoot();
166 alphaZ *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
167
168 // create the plane surface
169 Trk::PlaneSurface startSurface(createTransform(startX,
170 startY,
171 startZ,
172 phi, theta,
173 alphaZ),
174 10e3,10e3);
175
176 // make a covariance Matrix
177 AmgMatrix(5,5) covMat;
178 covMat.setZero();
179
180 // the initial perigee with random numbers
181 Trk::AtaPlane startParameters(loc1,
182 loc2,
183 phi,
184 theta,
185 qOverP,
186 startSurface,
187 covMat);
188
189 ATH_MSG_VERBOSE( "Start Parameters : " << startParameters );
190
191
192 // destination position
193 double estimationR = m_minimumR + (m_maximumR-m_minimumR) * m_flatDist->shoot();
194
195
196 // --------------- propagate to find a first intersection ---------------------
197 Trk::CylinderSurface estimationCylinder(Amg::Transform3D(), estimationR, 10e10);
198
199 ATH_MSG_VERBOSE( "Cylinder to be intersected : " << estimationCylinder );
200
201 auto estimationParameters = m_propagator->propagateParameters(ctx,
202 startParameters,
203 estimationCylinder,
204 Trk::alongMomentum,
205 false,
206 *m_magFieldProperties);
207 if (!estimationParameters) {
208 ATH_MSG_VERBOSE( "Estimation of intersection did not work - skip event !" );
209 return StatusCode::SUCCESS;
210 }
211
212 ATH_MSG_VERBOSE( "Estimation Parameters: " << *estimationParameters );
213
214 const Amg::Vector3D& estimatedPosition = estimationParameters->position();
215
216 double estimationX = estimatedPosition.x();
217 double estimationY = estimatedPosition.y();
218 double estimationZ = estimatedPosition.z();
219
220 double estimationPhi = estimatedPosition.phi();
221 double estimationTheta = estimatedPosition.theta();
222
223
224
225 double rotateTrans = M_PI * m_flatDist->shoot();
226 rotateTrans *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
227
228 Amg::Transform3D surfaceTransform;
229
230 if (m_useAlignedSurfaces){
231 // create a radial vector
232 Amg::Vector3D radialVector(estimatedPosition.x(), estimatedPosition.y(), 0.);
233 Amg::Vector3D surfaceZdirection(radialVector.unit());
234 Amg::Vector3D surfaceYdirection(0.,0.,1.);
235 // the x direction
236 Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
237 // the rotation
238 Amg::RotationMatrix3D surfaceRotation;
239 surfaceRotation.col(0) = surfaceXdirection;
240 surfaceRotation.col(1) = surfaceYdirection;
241 surfaceRotation.col(2) = surfaceZdirection;
242 Amg::Transform3D nominalTransform(surfaceRotation, estimatedPosition);
243 surfaceTransform = Amg::Transform3D(nominalTransform*Amg::AngleAxis3D(rotateTrans,Amg::Vector3D(0.,0.,1.)));
244 } else
245 surfaceTransform = createTransform(estimationX,
246 estimationY,
247 estimationZ,
248 estimationPhi,
249 estimationTheta,
250 rotateTrans);
251
252 // cleanup for memory reasons
253
254 Trk::PlaneSurface destinationSurface(surfaceTransform,10e5 , 10e5);
255
256
257 // transport the start to the destination surface
258 std::optional<Trk::TransportJacobian> optTransportJacobian{};
259 AmgMatrix(5,5) testMatrix; testMatrix.setZero();
260 Trk::TransportJacobian currentStepJacobian(testMatrix);
261 double pathLimit = -1.;
262
263 auto trackParameters = m_propagator->propagate(ctx,
264 startParameters,
265 destinationSurface,
266 Trk::alongMomentum,
267 false,
268 *m_magFieldProperties,
269 optTransportJacobian,
270 pathLimit);
271
272 // --------------------- check if test propagation was successful ------------------------------
273 if (trackParameters && optTransportJacobian){
274
275 unsigned int recStep = 0;
276 const auto& transportJacobian = (*optTransportJacobian);
277 // [0] Transport Jacobian -----------------------------------------------------
278 ATH_MSG_VERBOSE( "TransportJacobian : " << transportJacobian );
279
280
281 // and now fill the variables
282 m_loc1loc1[recStep] = (transportJacobian)(0,0);
283 m_loc1loc2[recStep] = (transportJacobian)(0,1);
284 m_loc1phi[recStep] = (transportJacobian)(0,2);
285 m_loc1theta[recStep] = (transportJacobian)(0,3);
286 m_loc1qop[recStep] = (transportJacobian)(0,4);
287 m_loc1steps[recStep] = 0.;
288
289 m_loc2loc1[recStep] = (transportJacobian)(1,0);
290 m_loc2loc2[recStep] = (transportJacobian)(1,1);
291 m_loc2phi[recStep] = (transportJacobian)(1,2);
292 m_loc2theta[recStep] = (transportJacobian)(1,3);
293 m_loc2qop[recStep] = (transportJacobian)(1,4);
294 m_loc2steps[recStep] = 0.;
295
296 m_philoc1[recStep] = (transportJacobian)(2,0);
297 m_philoc2[recStep] = (transportJacobian)(2,1);
298 m_phiphi[recStep] = (transportJacobian)(2,2);
299 m_phitheta[recStep] = (transportJacobian)(2,3);
300 m_phiqop[recStep] = (transportJacobian)(2,4);
301 m_phisteps[recStep] = 0.;
302
303 m_thetaloc1[recStep] = (transportJacobian)(3,0);
304 m_thetaloc2[recStep] = (transportJacobian)(3,1);
305 m_thetaphi[recStep] = (transportJacobian)(3,2);
306 m_thetatheta[recStep] = (transportJacobian)(3,3);
307 m_thetaqop[recStep] = (transportJacobian)(3,4);
308 m_thetasteps[recStep] = 0.;
309
310 m_qoploc1[recStep] = (transportJacobian)(4,0);
311 m_qoploc2[recStep] = (transportJacobian)(4,1);
312 m_qopphi[recStep] = (transportJacobian)(4,2);
313 m_qoptheta[recStep] = (transportJacobian)(4,3);
314 m_qopqop[recStep] = (transportJacobian)(4,4);
315 m_qopsteps[recStep] = 0.;
316
317 ++recStep;
318
319 // start the riddlers algorithm
320 // [1-2-3] Riddlers jacobians -----------------------------------------
321 for (unsigned int istep = 0; istep < m_localVariations.size(); ++istep){
322 // --------------------------
323 ATH_MSG_VERBOSE( "Performing step : " << istep );
324 // the initial perigee with random numbers
325 // for the first row
326 Trk::AtaPlane startLoc1Minus(loc1-m_localVariations[istep],loc2,phi,theta,qOverP,startSurface);
327 Trk::AtaPlane startLoc1Plus(loc1+m_localVariations[istep],loc2,phi,theta,qOverP,startSurface);
328 // for the second row
329 Trk::AtaPlane startLoc2Minus(loc1,loc2-m_localVariations[istep],phi,theta,qOverP,startSurface);
330 Trk::AtaPlane startLoc2Plus(loc1,loc2+m_localVariations[istep],phi,theta,qOverP,startSurface);
331 // for the third row
332 Trk::AtaPlane startPhiMinus(loc1,loc2,phi-m_angularVariations[istep],theta,qOverP,startSurface);
333 Trk::AtaPlane startPhiPlus(loc1,loc2,phi+m_angularVariations[istep],theta,qOverP,startSurface);
334 // for the fourth row
335 Trk::AtaPlane startThetaMinus(loc1,loc2,phi,theta-m_angularVariations[istep],qOverP,startSurface);
336 Trk::AtaPlane startThetaPlus(loc1,loc2,phi,theta+m_angularVariations[istep],qOverP,startSurface);
337 // for the fifth row
338 Trk::AtaPlane startQopMinus(loc1,loc2,phi,theta,qOverP-m_qOpVariations[istep],startSurface);
339 Trk::AtaPlane startQopPlus(loc1,loc2,phi,theta,qOverP+m_qOpVariations[istep],startSurface);
340
341 // the propagations --- 10 times
342 auto endLoc1Minus = m_propagator->propagateParameters(ctx,
343 startLoc1Minus,
344 destinationSurface,
345 Trk::alongMomentum,
346 false,
347 *m_magFieldProperties);
348
349
350 auto endLoc1Plus = m_propagator->propagateParameters(ctx,
351 startLoc1Plus,
352 destinationSurface,
353 Trk::alongMomentum,
354 false,
355 *m_magFieldProperties);
356
357 auto endLoc2Minus = m_propagator->propagateParameters(ctx,
358 startLoc2Minus,
359 destinationSurface,
360 Trk::alongMomentum,
361 false,
362 *m_magFieldProperties);
363
364 auto endLoc2Plus = m_propagator->propagateParameters(ctx,
365 startLoc2Plus,
366 destinationSurface,
367 Trk::alongMomentum,
368 false,
369 *m_magFieldProperties);
370
371 auto endPhiMinus = m_propagator->propagateParameters(ctx,
372 startPhiMinus,
373 destinationSurface,
374 Trk::alongMomentum,
375 false,
376 *m_magFieldProperties);
377
378 auto endPhiPlus = m_propagator->propagateParameters(ctx,
379 startPhiPlus,
380 destinationSurface,
381 Trk::alongMomentum,
382 false,
383 *m_magFieldProperties);
384
385 auto endThetaMinus = m_propagator->propagateParameters(ctx,
386 startThetaMinus,
387 destinationSurface,
388 Trk::alongMomentum,
389 false,
390 *m_magFieldProperties);
391
392 auto endThetaPlus = m_propagator->propagateParameters(ctx,
393 startThetaPlus,
394 destinationSurface,
395 Trk::alongMomentum,
396 false,
397 *m_magFieldProperties);
398
399 auto endQopMinus = m_propagator->propagateParameters(ctx,
400 startQopMinus,
401 destinationSurface,
402 Trk::alongMomentum,
403 false,
404 *m_magFieldProperties);
405
406 auto endQopPlus = m_propagator->propagateParameters(ctx,
407 startQopPlus,
408 destinationSurface,
409 Trk::alongMomentum,
410 false,
411 *m_magFieldProperties);
412 if (endLoc1Minus
413 && endLoc1Plus
414 && endLoc2Minus
415 && endLoc2Plus
416 && endPhiMinus
417 && endPhiPlus
418 && endThetaMinus
419 && endThetaPlus
420 && endQopMinus
421 && endQopPlus){
422
423 // Loc 1
424 const Amg::VectorX& endLoc1MinusPar = endLoc1Minus->parameters();
425 const Amg::VectorX& endLoc1PlusPar = endLoc1Plus->parameters();
426 // Loc 2
427 const Amg::VectorX& endLoc2MinusPar = endLoc2Minus->parameters();
428 const Amg::VectorX& endLoc2PlusPar = endLoc2Plus->parameters();
429 // Phi
430 const Amg::VectorX& endPhiMinusPar = endPhiMinus->parameters();
431 const Amg::VectorX& endPhiPlusPar = endPhiPlus->parameters();
432 // Theta
433 const Amg::VectorX& endThetaMinusPar = endThetaMinus->parameters();
434 const Amg::VectorX& endThetaPlusPar = endThetaPlus->parameters();
435 // qop
436 const Amg::VectorX& endQopMinusPar = endQopMinus->parameters();
437 const Amg::VectorX& endQopPlusPar = endQopPlus->parameters();
438
439 // the deltas
440 Amg::VectorX endLoc1Diff(endLoc1PlusPar-endLoc1MinusPar);
441 Amg::VectorX endLoc2Diff(endLoc2PlusPar-endLoc2MinusPar);
442 Amg::VectorX endPhiDiff(endPhiPlusPar-endPhiMinusPar);
443 Amg::VectorX endThetaDiff(endThetaPlusPar-endThetaMinusPar);
444 Amg::VectorX endQopDiff(endQopPlusPar-endQopMinusPar);
445
446 currentStepJacobian(0,0) = endLoc1Diff[0]/(2.*m_localVariations[istep]);
447 currentStepJacobian(0,1) = endLoc2Diff[0]/(2.*m_localVariations[istep]);
448 currentStepJacobian(0,2) = endPhiDiff[0]/(2.*m_angularVariations[istep]);
449 currentStepJacobian(0,3) = endThetaDiff[0]/(2.*m_angularVariations[istep]);
450 currentStepJacobian(0,4) = endQopDiff[0]/(2.*m_qOpVariations[istep]);
451
452 m_loc1loc1[recStep] = currentStepJacobian(0,0);
453 m_loc1loc2[recStep] = currentStepJacobian(0,1);
454 m_loc1phi[recStep] = currentStepJacobian(0,2);
455 m_loc1theta[recStep] = currentStepJacobian(0,3);
456 m_loc1qop[recStep] = currentStepJacobian(0,4);
457 m_loc1steps[recStep] = m_localVariations[istep];
458
459 currentStepJacobian(1,0) = endLoc1Diff[1]/(2.*m_localVariations[istep]);
460 currentStepJacobian(1,1) = endLoc2Diff[1]/(2.*m_localVariations[istep]);
461 currentStepJacobian(1,2) = endPhiDiff[1]/(2.*m_angularVariations[istep]);
462 currentStepJacobian(1,3) = endThetaDiff[1]/(2.*m_angularVariations[istep]);
463 currentStepJacobian(1,4) = endQopDiff[1]/(2.*m_qOpVariations[istep]);
464
465 m_loc2loc1[recStep] = currentStepJacobian(1,0);
466 m_loc2loc2[recStep] = currentStepJacobian(1,1);
467 m_loc2phi[recStep] = currentStepJacobian(1,2);
468 m_loc2theta[recStep] = currentStepJacobian(1,3);
469 m_loc2qop[recStep] = currentStepJacobian(1,4);
470 m_loc2steps[recStep] = m_localVariations[istep];
471
472 currentStepJacobian(2,0) = endLoc1Diff[2]/(2.*m_localVariations[istep]);
473 currentStepJacobian(2,1) = endLoc2Diff[2]/(2.*m_localVariations[istep]);
474 currentStepJacobian(2,2) = endPhiDiff[2]/(2.*m_angularVariations[istep]);
475 currentStepJacobian(2,3) = endThetaDiff[2]/(2.*m_angularVariations[istep]);
476 currentStepJacobian(2,4) = endQopDiff[2]/(2.*m_qOpVariations[istep]);
477
478 m_philoc1[recStep] = currentStepJacobian(2,0);
479 m_philoc2[recStep] = currentStepJacobian(2,1);
480 m_phiphi[recStep] = currentStepJacobian(2,2);
481 m_phitheta[recStep] = currentStepJacobian(2,3);
482 m_phiqop[recStep] = currentStepJacobian(2,4);
483 m_phisteps[recStep] = m_angularVariations[istep];
484
485 currentStepJacobian(3,0) = endLoc1Diff[3]/(2.*m_localVariations[istep]);
486 currentStepJacobian(3,1) = endLoc2Diff[3]/(2.*m_localVariations[istep]);
487 currentStepJacobian(3,2) = endPhiDiff[3]/(2.*m_angularVariations[istep]);
488 currentStepJacobian(3,3) = endThetaDiff[3]/(2.*m_angularVariations[istep]);
489 currentStepJacobian(3,4) = endQopDiff[3]/(2.*m_qOpVariations[istep]);
490
491 m_thetaloc1[recStep] = currentStepJacobian(3,0);
492 m_thetaloc2[recStep] = currentStepJacobian(3,1);
493 m_thetaphi[recStep] = currentStepJacobian(3,2);
494 m_thetatheta[recStep] = currentStepJacobian(3,3);
495 m_thetaqop[recStep] = currentStepJacobian(3,4);
496 m_thetasteps[recStep] = m_angularVariations[istep];
497
498 currentStepJacobian(4,0) = endLoc1Diff[4]/(2.*m_localVariations[istep]);
499 currentStepJacobian(4,1) = endLoc2Diff[4]/(2.*m_localVariations[istep]);
500 currentStepJacobian(4,2) = endPhiDiff[4]/(2.*m_angularVariations[istep]);
501 currentStepJacobian(4,3) = endThetaDiff[4]/(2.*m_angularVariations[istep]);
502 currentStepJacobian(4,4) = endQopDiff[4]/(2.*m_qOpVariations[istep]);
503
504 m_qoploc1[recStep] = currentStepJacobian(4,0);
505 m_qoploc2[recStep] = currentStepJacobian(4,1);
506 m_qopphi[recStep] = currentStepJacobian(4,2);
507 m_qoptheta[recStep] = currentStepJacobian(4,3);
508 m_qopqop[recStep] = currentStepJacobian(4,4);
509 m_qopsteps[recStep] = m_qOpVariations[istep];
510
511 ATH_MSG_DEBUG( "Current TransportJacobian : " << currentStepJacobian );
512
513 ++recStep;
514 }
515 }
516
517 // -------------------------------------------------------------------------------
518 if (recStep > 2){
519 // The parabolic interpolation -------------------------------------------------------------------------------
520 // dL1
521 m_loc1loc1[recStep] = parabolicInterpolation(m_loc1loc1[recStep-1],m_loc1loc1[recStep-2],m_loc1loc1[recStep-3],
522 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
523 m_loc1loc2[recStep] = parabolicInterpolation(m_loc1loc2[recStep-1],m_loc1loc2[recStep-2],m_loc1loc2[recStep-3],
524 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
525 m_loc1phi[recStep] = parabolicInterpolation(m_loc1phi[recStep-1],m_loc1phi[recStep-2],m_loc1phi[recStep-3],
526 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
527 m_loc1theta[recStep] = parabolicInterpolation(m_loc1theta[recStep-1],m_loc1theta[recStep-2],m_loc1theta[recStep-3],
528 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
529 m_loc1qop[recStep] = parabolicInterpolation(m_loc1qop[recStep-1],m_loc1qop[recStep-2],m_loc1qop[recStep-3],
530 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
531 m_loc1steps[recStep] = 1;
532 // dL2
533 m_loc2loc1[recStep] = parabolicInterpolation(m_loc2loc1[recStep-1],m_loc2loc1[recStep-2],m_loc2loc1[recStep-3],
534 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
535 m_loc2loc2[recStep] = parabolicInterpolation(m_loc2loc2[recStep-1],m_loc2loc2[recStep-2],m_loc2loc2[recStep-3],
536 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
537 m_loc2phi[recStep] = parabolicInterpolation(m_loc2phi[recStep-1],m_loc2phi[recStep-2],m_loc2phi[recStep-3],
538 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
539 m_loc2theta[recStep] = parabolicInterpolation(m_loc2theta[recStep-1],m_loc2theta[recStep-2],m_loc2theta[recStep-3],
540 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
541 m_loc2qop[recStep] = parabolicInterpolation(m_loc2qop[recStep-1],m_loc2qop[recStep-2],m_loc2qop[recStep-3],
542 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
543 m_loc2steps[recStep] = 1;
544 // dPhi
545 m_philoc1[recStep] = parabolicInterpolation(m_philoc1[recStep-1],m_philoc1[recStep-2],m_philoc1[recStep-3],
546 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
547 m_philoc2[recStep] = parabolicInterpolation(m_philoc2[recStep-1],m_philoc2[recStep-2],m_philoc2[recStep-3],
548 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
549 m_phiphi[recStep] = parabolicInterpolation(m_phiphi[recStep-1],m_phiphi[recStep-2],m_phiphi[recStep-3],
550 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
551 m_phitheta[recStep] = parabolicInterpolation(m_phitheta[recStep-1],m_phitheta[recStep-2],m_phitheta[recStep-3],
552 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
553 m_phiqop[recStep] = parabolicInterpolation(m_phiqop[recStep-1],m_phiqop[recStep-2],m_phiqop[recStep-3],
554 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
555 m_phisteps[recStep] = 1;
556 // dTheta
557 m_thetaloc1[recStep] = parabolicInterpolation(m_thetaloc1[recStep-1],m_thetaloc1[recStep-2],m_thetaloc1[recStep-3],
558 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
559 m_thetaloc2[recStep] = parabolicInterpolation(m_thetaloc2[recStep-1],m_thetaloc2[recStep-2],m_thetaloc2[recStep-3],
560 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
561 m_thetaphi[recStep] = parabolicInterpolation(m_thetaphi[recStep-1],m_thetaphi[recStep-2],m_thetaphi[recStep-3],
562 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
563 m_thetatheta[recStep] = parabolicInterpolation(m_thetatheta[recStep-1],m_thetatheta[recStep-2],m_thetatheta[recStep-3],
564 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
565 m_thetaqop[recStep] = parabolicInterpolation(m_thetaqop[recStep-1],m_thetaqop[recStep-2],m_thetaqop[recStep-3],
566 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
567 m_thetasteps[recStep] = 1;
568 // dTheta
569 m_qoploc1[recStep] = parabolicInterpolation(m_qoploc1[recStep-1],m_qoploc1[recStep-2],m_qoploc1[recStep-3],
570 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
571 m_qoploc2[recStep] = parabolicInterpolation(m_qoploc2[recStep-1],m_qoploc2[recStep-2],m_qoploc2[recStep-3],
572 m_localVariations[2], m_localVariations[1], m_localVariations[0]);
573 m_qopphi[recStep] = parabolicInterpolation(m_qopphi[recStep-1],m_qopphi[recStep-2],m_qopphi[recStep-3],
574 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
575 m_qoptheta[recStep] = parabolicInterpolation(m_qoptheta[recStep-1],m_qoptheta[recStep-2],m_qoptheta[recStep-3],
576 m_angularVariations[2], m_angularVariations[1], m_angularVariations[0]);
577 m_qopqop[recStep] = parabolicInterpolation(m_qopqop[recStep-1],m_qopqop[recStep-2],m_qopqop[recStep-3],
578 m_qOpVariations[2], m_qOpVariations[1], m_qOpVariations[0]);
579 m_qopsteps[recStep] = 1;
580
581 currentStepJacobian(0,0)= m_loc1loc1[recStep];
582 currentStepJacobian(0,1)= m_loc1loc2[recStep];
583 currentStepJacobian(0,2)= m_loc1phi[recStep];
584 currentStepJacobian(0,3)= m_loc1theta[recStep];
585 currentStepJacobian(0,4)= m_loc1qop[recStep];
586
587 currentStepJacobian(1,0)= m_loc2loc1[recStep];
588 currentStepJacobian(1,1)= m_loc2loc2[recStep];
589 currentStepJacobian(1,2)= m_loc2phi[recStep];
590 currentStepJacobian(1,3)= m_loc2theta[recStep];
591 currentStepJacobian(1,4)= m_loc2qop[recStep];
592
593 currentStepJacobian(2,0)= m_philoc1[recStep];
594 currentStepJacobian(2,1)= m_philoc2[recStep];
595 currentStepJacobian(2,2)= m_phiphi[recStep];
596 currentStepJacobian(2,3)= m_phitheta[recStep];
597 currentStepJacobian(2,4)= m_phiqop[recStep];
598
599 currentStepJacobian(3,0)= m_thetaloc1[recStep];
600 currentStepJacobian(3,1)= m_thetaloc2[recStep];
601 currentStepJacobian(3,2)= m_thetaphi[recStep];
602 currentStepJacobian(3,3)= m_thetatheta[recStep];
603 currentStepJacobian(3,4)= m_thetaqop[recStep];
604
605 currentStepJacobian(4,0)= m_qoploc1[recStep];
606 currentStepJacobian(4,1)= m_qoploc2[recStep];
607 currentStepJacobian(4,2)= m_qopphi[recStep];
608 currentStepJacobian(4,3)= m_qoptheta[recStep];
609 currentStepJacobian(4,4)= m_qopqop[recStep];
610
611 }
612
613 ATH_MSG_DEBUG( "Interpolated TransportJacobian : " << currentStepJacobian );
614 ++recStep;
615
616
617 // now fill the results
618 TransportJacobian diffMatrix(transportJacobian-currentStepJacobian);
619
620 ATH_MSG_VERBOSE( "Absolute Differences of the TransportJacobian : " << diffMatrix );
621
622 // Absolute differences ----------------------------------------------------
623 // (A1)
624 // fill the differences into the last one log (loc1)
625 m_loc1loc1[recStep] = diffMatrix(0,0);
626 m_loc1loc2[recStep] = diffMatrix(0,1);
627 m_loc1phi[recStep] = diffMatrix(0,2);
628 m_loc1theta[recStep] = diffMatrix(0,3);
629 m_loc1qop[recStep] = diffMatrix(0,4);
630 m_loc1steps[recStep] = 2;
631 // fill the differences into the last one log (loc2)
632 m_loc2loc1[recStep] = diffMatrix(1,0);
633 m_loc2loc2[recStep] = diffMatrix(1,1);
634 m_loc2phi[recStep] = diffMatrix(1,2);
635 m_loc2theta[recStep] = diffMatrix(1,3);
636 m_loc2qop[recStep] = diffMatrix(1,4);
637 m_loc2steps[recStep] = 2;
638 // fill the differences into the last one log (phi)
639 m_philoc1[recStep] = diffMatrix(2,0);
640 m_philoc2[recStep] = diffMatrix(2,1);
641 m_phiphi[recStep] = diffMatrix(2,2);
642 m_phitheta[recStep] = diffMatrix(2,3);
643 m_phiqop[recStep] = diffMatrix(2,4);
644 m_phisteps[recStep] = 2;
645 // fill the differences into the last one log (theta)
646 m_thetaloc1[recStep] = diffMatrix(3,0);
647 m_thetaloc2[recStep] = diffMatrix(3,1);
648 m_thetaphi[recStep] = diffMatrix(3,2);
649 m_thetatheta[recStep] = diffMatrix(3,3);
650 m_thetaqop[recStep] = diffMatrix(3,4);
651 m_thetasteps[recStep] = 2;
652 // fill the differences into the last one log (qop)
653 m_qoploc1[recStep] = diffMatrix(4,0);
654 m_qoploc2[recStep] = diffMatrix(4,1);
655 m_qopphi[recStep] = diffMatrix(4,2);
656 m_qoptheta[recStep] = diffMatrix(4,3);
657 m_qopqop[recStep] = diffMatrix(4,4);
658 m_qopsteps[recStep] = 2;
659 ++recStep;
660
661 // (A2)
662 if (recStep > 1){
663 // fill the differences into the last one log (loc1)
664 m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep-1] )) : 0.;
665 m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep-1] )) : 0.;
666 m_loc1phi[recStep] = std::abs(m_loc1phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep-1] )) : 0.;
667 m_loc1theta[recStep] = std::abs(m_loc1theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep-1])) : 0.;
668 m_loc1qop[recStep] = std::abs(m_loc1qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep-1] )) : 0.;
669 m_loc1steps[recStep] = 3;
670 // fill the differences into the last one log (loc2)
671 m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep-1] )) : 0.;
672 m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep-1] )) : 0.;
673 m_loc2phi[recStep] = std::abs(m_loc2phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep-1] )) : 0.;
674 m_loc2theta[recStep] = std::abs(m_loc2theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep-1])) : 0.;
675 m_loc2qop[recStep] = std::abs(m_loc2qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep-1] )) : 0.;
676 m_loc2steps[recStep] = 3;
677 // fill the differences into the last one log (phi)
678 m_philoc1[recStep] = std::abs(m_philoc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc1[recStep-1] )) : 0.;
679 m_philoc2[recStep] = std::abs(m_philoc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc2[recStep-1] )) : 0.;
680 m_phiphi[recStep] = std::abs(m_phiphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiphi[recStep-1] )) : 0.;
681 m_phitheta[recStep] = std::abs(m_phitheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phitheta[recStep-1])) : 0.;
682 m_phiqop[recStep] = std::abs(m_phiqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiqop[recStep-1] )) : 0.;
683 m_phisteps[recStep] = 3;
684 // fill the differences into the last one log (theta)
685 m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep-1] )) : 0.;
686 m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep-1] )) : 0.;
687 m_thetaphi[recStep] = std::abs(m_thetaphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep-1] )) : 0.;
688 m_thetatheta[recStep] = std::abs(m_thetatheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep-1])) : 0.;
689 m_thetaqop[recStep] = std::abs(m_thetaqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep-1] )) : 0.;
690 m_thetasteps[recStep] = 3;
691 // fill the differences into the last one log (qop)
692 m_qoploc1[recStep] = std::abs(m_qoploc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoploc1[recStep-1] )) : 0.;
693 m_qoploc2[recStep] = std::abs(m_qoploc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoploc2[recStep-1] )) : 0.;
694 m_qopphi[recStep] = std::abs(m_qopphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qopphi[recStep-1] )) : 0.;
695 m_qoptheta[recStep] = std::abs(m_qoptheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qoptheta[recStep-1])) : 0.;
696 m_qopqop[recStep] = std::abs(m_qopqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_qopqop[recStep-1] )) : 0.;
697 m_qopsteps[recStep] = 3;
698 }
699 ++recStep;
700
701
702 // Relative differences) ----------------------------------------------------
703 // (R1)
704 // fill the differences into the last one log (loc1)
705 m_loc1loc1[recStep] = std::abs((transportJacobian)(0,0)) > 1e-50 ? diffMatrix(0,0)/((transportJacobian)(0,0)) : 0.;
706 m_loc1loc2[recStep] = std::abs((transportJacobian)(0,1)) > 1e-50 ? diffMatrix(0,1)/((transportJacobian)(0,1)) : 0.;
707 m_loc1phi[recStep] = std::abs((transportJacobian)(0,2)) > 1e-50 ? diffMatrix(0,2)/((transportJacobian)(0,2)) : 0.;
708 m_loc1theta[recStep] = std::abs((transportJacobian)(0,3)) > 1e-50 ? diffMatrix(0,3)/((transportJacobian)(0,3)) : 0.;
709 m_loc1qop[recStep] = std::abs((transportJacobian)(0,4)) > 1e-50 ? diffMatrix(0,4)/((transportJacobian)(0,4)) : 0.;
710 m_loc1steps[recStep] = 4;
711 // fill the differences into the last one log (loc2)
712 m_loc2loc1[recStep] = std::abs((transportJacobian)(1,0)) > 1e-50 ? diffMatrix(1,0)/((transportJacobian)(1,0)) : 0.;
713 m_loc2loc2[recStep] = std::abs((transportJacobian)(1,1)) > 1e-50 ? diffMatrix(1,1)/((transportJacobian)(1,1)) : 0.;
714 m_loc2phi[recStep] = std::abs((transportJacobian)(1,2)) > 1e-50 ? diffMatrix(1,2)/((transportJacobian)(1,2)) : 0.;
715 m_loc2theta[recStep] = std::abs((transportJacobian)(1,3)) > 1e-50 ? diffMatrix(1,3)/((transportJacobian)(1,3)) : 0.;
716 m_loc2qop[recStep] = std::abs((transportJacobian)(1,4)) > 1e-50 ? diffMatrix(1,4)/((transportJacobian)(1,4)) : 0.;
717 m_loc2steps[recStep] = 4;
718 // fill the differences into the last one log (phi)
719 m_philoc1[recStep] = std::abs((transportJacobian)(2,0)) > 1e-50 ? diffMatrix(2,0)/((transportJacobian)(2,0)) : 0.;
720 m_philoc2[recStep] = std::abs((transportJacobian)(2,1)) > 1e-50 ? diffMatrix(2,1)/((transportJacobian)(2,1)) : 0.;
721 m_phiphi[recStep] = std::abs((transportJacobian)(2,2)) > 1e-50 ? diffMatrix(2,2)/((transportJacobian)(2,2)) : 0.;
722 m_phitheta[recStep] = std::abs((transportJacobian)(2,3)) > 1e-50 ? diffMatrix(2,3)/((transportJacobian)(2,3)) : 0.;
723 m_phiqop[recStep] = std::abs((transportJacobian)(2,4)) > 1e-50 ? diffMatrix(2,4)/((transportJacobian)(2,4)) : 0.;
724 m_phisteps[recStep] = 4;
725 // fill the differences into the last one log (theta)
726 m_thetaloc1[recStep] = std::abs((transportJacobian)(3,0)) > 1e-50 ? diffMatrix(3,0)/((transportJacobian)(3,0)) : 0.;
727 m_thetaloc2[recStep] = std::abs((transportJacobian)(3,1)) > 1e-50 ? diffMatrix(3,1)/((transportJacobian)(3,1)) : 0.;
728 m_thetaphi[recStep] = std::abs((transportJacobian)(3,2)) > 1e-50 ? diffMatrix(3,2)/((transportJacobian)(3,2)) : 0.;
729 m_thetatheta[recStep] = std::abs((transportJacobian)(3,3)) > 1e-50 ? diffMatrix(3,3)/((transportJacobian)(3,3)) : 0.;
730 m_thetaqop[recStep] = std::abs((transportJacobian)(3,4)) > 1e-50 ? diffMatrix(3,4)/((transportJacobian)(3,4)) : 0.;
731 m_thetasteps[recStep] = 4;
732 // fill the differences into the last one log (qop)
733 m_qoploc1[recStep] = std::abs((transportJacobian)(4,0)) > 1e-50 ? diffMatrix(4,0)/((transportJacobian)(4,0)) : 0.;
734 m_qoploc2[recStep] = std::abs((transportJacobian)(4,1)) > 1e-50 ? diffMatrix(4,1)/((transportJacobian)(4,1)) : 0.;
735 m_qopphi[recStep] = std::abs((transportJacobian)(4,2)) > 1e-50 ? diffMatrix(4,2)/((transportJacobian)(4,2)) : 0.;
736 m_qoptheta[recStep] = std::abs((transportJacobian)(4,3)) > 1e-50 ? diffMatrix(4,3)/((transportJacobian)(4,3)) : 0.;
737 m_qopqop[recStep] = std::abs((transportJacobian)(4,4)) > 1e-50 ? diffMatrix(4,4)/((transportJacobian)(4,4)) : 0.;
738 m_qopsteps[recStep] = 4;
739 ++recStep;
740
741 // (R2)
742 // Relative differences ----------------------------------------------------
743 m_loc1loc1[recStep] = std::abs(m_loc1loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc1[recStep-1] )) : 0.;
744 m_loc1loc2[recStep] = std::abs(m_loc1loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1loc2[recStep-1] )) : 0.;
745 m_loc1phi[recStep] = std::abs(m_loc1phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1phi[recStep-1] )) : 0.;
746 m_loc1theta[recStep] = std::abs(m_loc1theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1theta[recStep-1])) : 0.;
747 m_loc1qop[recStep] = std::abs(m_loc1qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc1qop[recStep-1] )) : 0.;
748 m_loc1steps[recStep] = 5;
749 // fill the differences into the last one log (loc2)
750 m_loc2loc1[recStep] = std::abs(m_loc2loc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc1[recStep-1] )) : 0.;
751 m_loc2loc2[recStep] = std::abs(m_loc2loc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2loc2[recStep-1] )) : 0.;
752 m_loc2phi[recStep] = std::abs(m_loc2phi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2phi[recStep-1] )) : 0.;
753 m_loc2theta[recStep] = std::abs(m_loc2theta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2theta[recStep-1])) : 0.;
754 m_loc2qop[recStep] = std::abs(m_loc2qop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_loc2qop[recStep-1] )) : 0.;
755 m_loc2steps[recStep] = 5;
756 // fill the differences into the last one log (phi)
757 m_philoc1[recStep] = std::abs(m_philoc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc1[recStep-1] )) : 0.;
758 m_philoc2[recStep] = std::abs(m_philoc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_philoc2[recStep-1] )) : 0.;
759 m_phiphi[recStep] = std::abs(m_phiphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiphi[recStep-1] )) : 0.;
760 m_phitheta[recStep] = std::abs(m_phitheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phitheta[recStep-1])) : 0.;
761 m_phiqop[recStep] = std::abs(m_phiqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_phiqop[recStep-1] )) : 0.;
762 m_phisteps[recStep] = 5;
763 // fill the differences into the last one log (theta)
764 m_thetaloc1[recStep] = std::abs(m_thetaloc1[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc1[recStep-1] )) : 0.;
765 m_thetaloc2[recStep] = std::abs(m_thetaloc2[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaloc2[recStep-1] )) : 0.;
766 m_thetaphi[recStep] = std::abs(m_thetaphi[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaphi[recStep-1] )) : 0.;
767 m_thetatheta[recStep] = std::abs(m_thetatheta[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetatheta[recStep-1])) : 0.;
768 m_thetaqop[recStep] = std::abs(m_thetaqop[recStep-1] ) > 1e-50 ? -std::log10(std::abs(m_thetaqop[recStep-1] )) : 0.;
769 m_thetasteps[recStep] = 5;
770 // fill the differences into the last one log (qop)
771 m_qoploc1[recStep] = std::abs(m_qoploc1[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoploc1[recStep-1] ) ): 0.;
772 m_qoploc2[recStep] = std::abs(m_qoploc2[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoploc2[recStep-1] ) ): 0.;
773 m_qopphi[recStep] = std::abs(m_qopphi[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qopphi[recStep-1] ) ): 0.;
774 m_qoptheta[recStep] = std::abs(m_qoptheta[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qoptheta[recStep-1]) ): 0.;
775 m_qopqop[recStep] = std::abs(m_qopqop[recStep-1] ) > 1e-50 ? -std::log10( std::abs(m_qopqop[recStep-1] ) ): 0.;
776 m_qopsteps[recStep] = 5;
777 ++recStep;
778
779 m_steps = recStep;
780 m_validationTree->Fill();
781 }
782
783 // Code entered here will be executed once per event
784 return StatusCode::SUCCESS;
785}
786
787//============================================================================================
788Amg::Transform3D
789Trk::RiddersAlgorithm::createTransform(double x, double y, double z, double phi, double theta, double alphaZ)
790{
791
792 if (phi!=0. && theta != 0.){
793 // create the Start Surface
794 Amg::Vector3D surfacePosition(x,y,z);
795 // z direction
796 Amg::Vector3D surfaceZdirection(cos(phi)*sin(theta),
797 sin(phi)*sin(theta),
798 cos(theta));
799 // the global z axis
800 Amg::Vector3D zAxis(0.,0.,1.);
801 // the y direction
802 Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
803 // the x direction
804 Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
805 // the rotation
806 Amg::RotationMatrix3D surfaceRotation;
807 surfaceRotation.col(0) = surfaceXdirection;
808 surfaceRotation.col(1) = surfaceYdirection;
809 surfaceRotation.col(2) = surfaceZdirection;
810 Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
811 return Amg::Transform3D(nominalTransform*Amg::AngleAxis3D(alphaZ,zAxis));
812
813 }
814
815 return Amg::Transform3D(Amg::Translation3D(x,y,z));
816}
817
818double Trk::RiddersAlgorithm::parabolicInterpolation(double y0, double y1, double y2,
819 double x0, double x1, double x2)
820{
821 double Z0 = x1*x2*y0;
822 double Z1 = x0*x2*y1;
823 double Z2 = x0*x1*y2;
824 double N0 = (x0-x1)*(x0-x2);
825 double N1 = (x1-x2)*(x1-x0);
826 double N2 = (x2-x0)*(x2-x1);
827 return Z0/N0 + Z1/N1 + Z2/N2;
828}
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:24
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:818
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:789
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()
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.14.0