ATLAS Offline Software
ExtrapolationValidation.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // ExtrapolationValidation.cxx, (c) ATLAS Detector software
8 
9 // Tracking
10 #include "TrkExAlgs/ExtrapolationValidation.h"
11 #include "TrkExInterfaces/IExtrapolator.h"
12 #include "TrkSurfaces/CylinderSurface.h"
13 #include "TrkSurfaces/DiscSurface.h"
14 #include "TrkSurfaces/PlaneSurface.h"
15 #include "TrkEventPrimitives/ParticleHypothesis.h"
16 #include "TrkGeometry/TrackingGeometry.h"
17 #include "TrkGeometry/TrackingVolume.h"
18 #include "TrkVolumes/CylinderVolumeBounds.h"
19 #include "TrkTrack/TrackStateOnSurface.h"
20 // Validation mode - TTree includes
21 #include "TTree.h"
22 #include "GaudiKernel/ITHistSvc.h"
23 #include "GaudiKernel/SystemOfUnits.h"
24 
25 //================ Constructor =================================================
26 
27 Trk::ExtrapolationValidation::ExtrapolationValidation(const std::string& name, ISvcLocator* pSvcLocator)
28  :
29  AthAlgorithm(name,pSvcLocator),
30  m_highestVolume(nullptr),
31  m_extrapolator("Trk::Extrapolator/AtlasExtrapolator"),
32  m_gaussDist(nullptr),
33  m_flatDist(nullptr),
34  m_materialCollectionValidation(true),
35  m_direct(false),
36  m_validationTree(nullptr),
37  m_validationTreeName("ExtrapolationValidation"),
38  m_validationTreeDescription("Output of the ExtrapolationValidation Algorithm"),
39  m_validationTreeFolder("/val/ExtrapolationValidation"),
40  m_maximumR(0.),
41  m_maximumZ(0.),
42  m_sigmaLoc(10.*Gaudi::Units::micrometer),
43  m_sigmaR(17.*Gaudi::Units::micrometer),
44  m_sigmaZ(50.*Gaudi::Units::millimeter),
45  m_minEta(-3.),
46  m_maxEta(3.),
47  m_minP(0.5*Gaudi::Units::GeV),
48  m_maxP(100.*Gaudi::Units::GeV),
49  m_particleType(2),
50  m_parameters(0),
51  m_parameterLoc1{},
52  m_parameterLoc2{},
53  m_parameterPhi{},
54  m_parameterTheta{},
55  m_parameterEta{},
56  m_parameterQoverP{},
57  m_covarianceLoc1{},
58  m_covarianceLoc2{},
59  m_covariancePhi{},
60  m_covarianceTheta{},
61  m_covarianceQoverP{},
62  m_covarianceDeterminant{},
63  m_destinationSurfaceType(0),
64  m_startX{},
65  m_startY{},
66  m_startR{},
67  m_startZ{},
68  m_startP{},
69  m_estimationX(0.),
70  m_estimationY(0.),
71  m_estimationR(0.),
72  m_estimationZ(0.),
73  m_destinationX(0.),
74  m_destinationY(0.),
75  m_destinationR(0.),
76  m_destinationZ(0.),
77  m_triesFront(0),
78  m_breaksFront(0),
79  m_triesBack(0),
80  m_breaksBack(0),
81  m_collectedLayerFront(0),
82  m_collectedLayerBack(0)
83 {
84  // used algorithms and alg tools
85  declareProperty("Extrapolator", m_extrapolator);
86  declareProperty("ValidateMaterialCollection", m_materialCollectionValidation);
87  declareProperty("ExtrapolateDirectly", m_direct);
88  // TTree handling
89  declareProperty("ValidationTreeName", m_validationTreeName);
90  declareProperty("ValidationTreeDescription", m_validationTreeDescription);
91  declareProperty("ValidationTreeFolder", m_validationTreeFolder);
92  // algorithm steering
93  declareProperty("StartPerigeeSigmaLoc" , m_sigmaLoc);
94  declareProperty("StartPerigeeSigmaR" , m_sigmaR);
95  declareProperty("StartPerigeeSigmaZ" , m_sigmaZ);
96  declareProperty("StartPerigeeMinEta" , m_minEta);
97  declareProperty("StartPerigeeMaxEta" , m_maxEta);
98  declareProperty("StartPerigeeMinP" , m_minP);
99  declareProperty("StartPerigeeMaxP" , m_maxP);
100  declareProperty("ParticleType" , m_particleType);
101 
102 }
103 
104 //================ Destructor =================================================
105 
106 Trk::ExtrapolationValidation::~ExtrapolationValidation()
107 {
108  // clear random number generators
109  delete m_gaussDist;
110  delete m_flatDist;
111 }
112 
113 
114 //================ Initialisation =================================================
115 
116 StatusCode Trk::ExtrapolationValidation::initialize()
117 {
118  // Code entered here will be executed once at program start.
119  ATH_MSG_INFO( " initialize()" );
120 
121  // Get Extrapolator from ToolService
122  ATH_CHECK( m_extrapolator.retrieve());
123 
124  // create the new Tree
125  m_validationTree = new TTree(m_validationTreeName.c_str(), m_validationTreeDescription.c_str());
126 
127  // the branches for the parameters
128  m_validationTree->Branch("Parameters", &m_parameters, "params/I");
129  m_validationTree->Branch("ParametersLoc1", m_parameterLoc1, "paramLoc1[params]/F");
130  m_validationTree->Branch("ParametersLoc2", m_parameterLoc2, "paramLoc2[params]/F");
131  m_validationTree->Branch("ParametersPhi", m_parameterPhi, "paramPhi[params]/F");
132  m_validationTree->Branch("ParametersTheta", m_parameterTheta, "paramTheta[params]/F");
133  m_validationTree->Branch("ParametersEta", m_parameterEta, "paramEta[params]/F");
134  m_validationTree->Branch("ParametersQoverP", m_parameterQoverP, "paramQoverP[params]/F");
135  // for the covariance diagonals
136  m_validationTree->Branch("CovarianceLoc1", m_covarianceLoc1, "covLoc1[params]/F");
137  m_validationTree->Branch("CovarianceLoc2", m_covarianceLoc2, "covLoc2[params]/F");
138  m_validationTree->Branch("CovariancePhi", m_covariancePhi, "covPhi[params]/F");
139  m_validationTree->Branch("CovarianceTheta", m_covarianceTheta, "covTheta[params]/F");
140  m_validationTree->Branch("CovarianceQoverP", m_covarianceQoverP, "covQoverP[params]/F");
141  m_validationTree->Branch("CovarianceDeterminant", m_covarianceDeterminant, "covDet[params]/F");
142  // the start Momentum
143  m_validationTree->Branch("StartMomentum", &m_startP, "startP/F");
144  // for the start surface
145  m_validationTree->Branch("StartSurfaceX", &m_startX, "startX/F");
146  m_validationTree->Branch("StartSurfaceY", &m_startY, "startY/F");
147  m_validationTree->Branch("StartSurfaceR", &m_startR, "startR/F");
148  m_validationTree->Branch("StartSurfaceZ", &m_startZ, "startZ/F");
149  // the estimation of the parameters
150  m_validationTree->Branch("EstimationSurfaceX", &m_estimationX, "estimateX/F");
151  m_validationTree->Branch("EstimationSurfaceY", &m_estimationY, "estimateY/F");
152  m_validationTree->Branch("EstimationSurfaceR", &m_estimationR, "estimateR/F");
153  m_validationTree->Branch("EstimationSurfaceZ", &m_estimationZ, "estimateZ/F");
154  // for the surface type
155  m_validationTree->Branch("DestinationSurfaceType", &m_destinationSurfaceType, "surfaceType/I");
156  m_validationTree->Branch("DestinationSurfaceX", &m_destinationX, "surfaceX/F");
157  m_validationTree->Branch("DestinationSurfaceY", &m_destinationY, "surfaceY/F");
158  m_validationTree->Branch("DestinationSurfaceR", &m_destinationR, "surfaceR/F");
159  m_validationTree->Branch("DestinationSurfaceZ", &m_destinationZ, "surfaceZ/F");
160 
161  // now register the Tree
162  ITHistSvc* tHistSvc = nullptr;
163  if (service("THistSvc",tHistSvc).isFailure()){
164  ATH_MSG_ERROR("initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
165  delete m_validationTree; m_validationTree = nullptr;
166  }
167  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
168  ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
169  delete m_validationTree; m_validationTree = nullptr;
170  }
171 
172  // intialize the random number generators
173  m_gaussDist = new Rndm::Numbers(randSvc(), Rndm::Gauss(0.,1.));
174  m_flatDist = new Rndm::Numbers(randSvc(), Rndm::Flat(0.,1.));
175 
176  ATH_MSG_INFO( "initialize() successful");
177  return StatusCode::SUCCESS;
178 }
179 
180 //================ Finalisation =================================================
181 
182 StatusCode Trk::ExtrapolationValidation::finalize()
183 {
184  // Code entered here will be executed once at the end of the program run.
185  ATH_MSG_INFO("================== Output Statistics =========================");
186  ATH_MSG_INFO("= Navigation : ");
187  ATH_MSG_INFO("= - breaks fwd : " << double(m_breaksFront)/double(m_triesFront)
188  << " (" << m_breaksFront << "/" << m_triesFront << ")");
189  ATH_MSG_INFO("= - breaks bwd : " << double(m_breaksBack)/double(m_triesBack)
190  << " (" << m_breaksBack << "/" << m_triesBack << ")");
191  if (m_materialCollectionValidation){
192  ATH_MSG_INFO("= Material collection : ");
193  ATH_MSG_INFO("= - layer collected fwd : " << m_collectedLayerFront );
194  ATH_MSG_INFO("= - layer collected bwd : " << m_collectedLayerBack );
195  }
196 
197  ATH_MSG_INFO("==============================================================");
198 
199  return StatusCode::SUCCESS;
200 }
201 
202 //================ Execution ====================================================
203 
204 StatusCode Trk::ExtrapolationValidation::execute()
205 {
206  const EventContext& ctx = Gaudi::Hive::currentContext();
207  // get the overall dimensions
208  if (!m_highestVolume){
209  // get TrackingGeometry and highest volume
210  const Trk::TrackingGeometry* trackingGeometry = m_extrapolator->trackingGeometry();
211  m_highestVolume = trackingGeometry ? trackingGeometry->highestTrackingVolume() : nullptr;
212  const Trk::CylinderVolumeBounds* cylBounds = m_highestVolume ?
213  dynamic_cast<const Trk::CylinderVolumeBounds*>(&(m_highestVolume->volumeBounds())) : nullptr;
214  // bail out
215  if (!cylBounds){
216  ATH_MSG_WARNING("No highest TrackingVolume / no VolumeBounds ... pretty useless! ");
217  return StatusCode::SUCCESS;
218  }
219  // get the numbers
220  m_maximumR = cylBounds->outerRadius();
221  m_maximumZ = cylBounds->halflengthZ();
222  }
223 
224  // intialize the values
225  m_parameters = 0;
226  m_destinationSurfaceType = 0;
227  // -----------> start
228  m_startX = 0.;
229  m_startY = 0.;
230  m_startR = 0.;
231  m_startZ = 0.;
232  // -----------> estimation
233  m_estimationX = 0.;
234  m_estimationY = 0.;
235  m_estimationR = 0.;
236  m_estimationZ = 0.;
237  // -----------> destination
238  m_destinationX = 0.;
239  m_destinationY = 0.;
240  m_destinationR = 0.;
241  m_destinationZ = 0.;
242 
243  // the local start parameters
244  m_parameterLoc1[m_parameters] = m_sigmaLoc * m_gaussDist->shoot();
245  m_parameterLoc2[m_parameters] = m_sigmaLoc * m_gaussDist->shoot();
246  // are adopted for planar and straight line surfaces
247  m_parameterPhi[m_parameters] = M_PI * m_flatDist->shoot();
248  m_parameterPhi[m_parameters] *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
249  m_parameterEta[m_parameters] = m_minEta + m_flatDist->shoot()*(m_maxEta-m_minEta);
250  m_parameterTheta[m_parameters] = 2.*atan(exp(-m_parameterEta[m_parameters]));
251 
252 
253  m_covarianceLoc1[m_parameters] = fabs( m_parameterLoc1[m_parameters] * 0.1);
254  m_covarianceLoc2[m_parameters] = fabs( m_parameterLoc2[m_parameters] * 0.1);
255  m_covariancePhi[m_parameters] = fabs( m_parameterPhi[m_parameters] * 0.1);
256  m_covarianceTheta[m_parameters] = fabs(m_parameterTheta[m_parameters] * 0.1);
257 
258  // this is fine
259  double p = m_minP + m_flatDist->shoot()*(m_maxP-m_minP);
260  double charge = (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
261  m_parameterQoverP[m_parameters] = charge/p;
262 
263  m_covarianceQoverP[m_parameters] = fabs(m_parameterQoverP[m_parameters] * 0.1);
264 
265  // for the momentum logging
266  m_startP = p;
267 
268  // start
269  m_startR = fabs(m_sigmaR * m_gaussDist->shoot());
270  double surfacePhi = M_PI * m_flatDist->shoot();
271  surfacePhi *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
272  m_startX = m_startR*cos(surfacePhi);
273  m_startY = m_startR*sin(surfacePhi);
274  m_startZ = m_sigmaZ * m_gaussDist->shoot();
275 
276  // rotate it around Z
277  double alphaZ = M_PI * m_flatDist->shoot();
278  alphaZ *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
279 
280  // create the plane surface
281  Trk::PlaneSurface startSurface(createTransform(m_startX,
282  m_startY,
283  m_startZ,
284  m_parameterPhi[m_parameters],
285  m_parameterTheta[m_parameters],
286  alphaZ),
287  10e3,10e3);
288 
289 
290 
291  AmgSymMatrix(5) covariance;
292  covariance.setZero();
293  ATH_MSG_VERBOSE(m_covarianceLoc1[m_parameters]);
294  covariance(0,0) = m_covarianceLoc1[m_parameters];
295  ATH_MSG_VERBOSE(m_covarianceLoc2[m_parameters]);
296  covariance(1,1) = m_covarianceLoc2[m_parameters];
297  ATH_MSG_VERBOSE(m_covariancePhi[m_parameters]);
298  covariance(2,2) = m_covariancePhi[m_parameters];
299  ATH_MSG_VERBOSE(m_covarianceTheta[m_parameters]);
300  covariance(3,3) = m_covarianceTheta[m_parameters];
301  ATH_MSG_VERBOSE(m_covarianceQoverP[m_parameters]);
302  covariance(4,4) = m_covarianceQoverP[m_parameters];
303  ATH_MSG_VERBOSE("Initial Setting: \n"<<covariance);
304 
305 
306  m_covarianceDeterminant[m_parameters] = covariance.determinant();
307 
308  // the initial perigee with random numbers
309  Trk::AtaPlane startParameters(m_parameterLoc1[m_parameters],
310  m_parameterLoc2[m_parameters],
311  m_parameterPhi[m_parameters],
312  m_parameterTheta[m_parameters],
313  m_parameterQoverP[m_parameters],
314  startSurface,
315  std::move(covariance));
316 
317  ATH_MSG_VERBOSE( "Start Parameters : " << startParameters );
318  if(startParameters.covariance())ATH_MSG_VERBOSE( "Start Covariance : \n" << *startParameters.covariance() );
319 
320 
321  // destination position
322  m_estimationR = m_maximumR * m_flatDist->shoot();
323 
324  // --------------- propagate to find a first intersection ---------------------
325  Amg::Transform3D CylTrf;
326  CylTrf.setIdentity();
327  Trk::CylinderSurface estimationCylinder(CylTrf, m_estimationR, 10e10);
328  const Trk::TrackParameters* estimationParameters = m_extrapolator->extrapolateDirectly(ctx,
329  startParameters,
330  estimationCylinder,
331  Trk::alongMomentum,
332  false).release();
333  if (!estimationParameters) {
334  ATH_MSG_VERBOSE( "Estimation of intersection did not work - skip event !" );
335  return StatusCode::SUCCESS;
336  }
337  else if (m_highestVolume && estimationParameters && !(m_highestVolume->inside(estimationParameters->position()))){
338  ATH_MSG_VERBOSE( "Estimation of intersection is outside the known world - skip event !" );
339  delete estimationParameters;
340  return StatusCode::SUCCESS;
341  }
342 
343  ATH_MSG_VERBOSE( "Estimation Parameters: " << *estimationParameters );
344 
345  // record the estimation parameters
346  ++m_triesFront;
347  ++m_parameters;
348  m_parameterLoc1[m_parameters] = estimationParameters->parameters()[Trk::loc1];
349  m_parameterLoc2[m_parameters] = estimationParameters->parameters()[Trk::loc2];
350  m_parameterPhi[m_parameters] = estimationParameters->parameters()[Trk::phi];
351  m_parameterEta[m_parameters] = estimationParameters->eta();
352  m_parameterTheta[m_parameters] = estimationParameters->parameters()[Trk::theta];
353  m_parameterQoverP[m_parameters] = estimationParameters->parameters()[Trk::qOverP];
354  if(estimationParameters->covariance()){
355  m_covarianceLoc1[m_parameters] = (*estimationParameters->covariance())(0,0);
356  m_covarianceLoc2[m_parameters] = (*estimationParameters->covariance())(1,1);
357  m_covariancePhi[m_parameters] = (*estimationParameters->covariance())(2,2);
358  m_covarianceTheta[m_parameters] = (*estimationParameters->covariance())(3,3);
359  m_covarianceQoverP[m_parameters] = (*estimationParameters->covariance())(4,4);
360  m_covarianceDeterminant[m_parameters] = (estimationParameters->covariance())->determinant();
361  }
362  else{
363  m_covarianceLoc1[m_parameters] = 0;
364  m_covarianceLoc2[m_parameters] = 0;
365  m_covariancePhi[m_parameters] = 0;
366  m_covarianceTheta[m_parameters] = 0;
367  m_covarianceQoverP[m_parameters] = 0;
368  m_covarianceDeterminant[m_parameters] = 0;
369  }
370  // the start Momentum
371 
372  // get the estimated position
373  const Amg::Vector3D& estimatedPosition = estimationParameters->position();
374 
375  m_estimationX = estimatedPosition.x();
376  m_estimationY = estimatedPosition.y();
377  m_estimationZ = estimatedPosition.z();
378 
379  // cleanup for memory reasons
380  delete estimationParameters; estimationParameters = nullptr;
381 
382  // create the radom surface at the destination point
383  Trk::PlaneSurface destinationSurface(createTransform(m_estimationX,
384  m_estimationY,
385  m_estimationZ,
386  m_parameterPhi[m_parameters],
387  m_parameterTheta[m_parameters]), 10e5 , 10e5);
388 
389 
390  ATH_MSG_VERBOSE( "Extrapolation to Destination Surface: " << destinationSurface );
391 
392  // the destination parameters
393  const Trk::TrackParameters* destParameters = nullptr;
394  // the standard validation ...
395  if (!m_materialCollectionValidation && !m_direct)
396  destParameters = m_extrapolator->extrapolate(ctx,
397  startParameters,
398  destinationSurface,
399  Trk::alongMomentum,
400  false,
401  (Trk::ParticleHypothesis)m_particleType,Trk::addNoise).release();
402  else if(!m_direct){ // material collection validation
403  // get the vector of TrackStateOnSurfaces back
404  const std::vector<const Trk::TrackStateOnSurface*>*
405  collectedMaterial = m_extrapolator->extrapolateM(ctx,
406  startParameters,
407  destinationSurface,
408  Trk::alongMomentum,
409  false,
410  (Trk::ParticleHypothesis)m_particleType);
411 
412  // get the last one and clone it
413  if (collectedMaterial && !collectedMaterial->empty()){
414  // get the last track state on surface & clone the destination parameters
415  const Trk::TrackStateOnSurface* destinationState = collectedMaterial->back();
416  destParameters = destinationState->trackParameters() ? destinationState->trackParameters()->clone() : nullptr;
417  m_collectedLayerFront += collectedMaterial->size();
418  // delete the layers / cleanup
419  std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter = collectedMaterial->begin();
420  std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIterEnd = collectedMaterial->end();
421  for ( ; tsosIter != tsosIterEnd; delete(*tsosIter), ++tsosIter);
422  }
423  }
424 
425  else{
426  destParameters = m_extrapolator->extrapolateDirectly(ctx,
427  startParameters,
428  destinationSurface,
429  Trk::alongMomentum,
430  false,
431  (Trk::ParticleHypothesis)m_particleType).release();
432 
433  }
434  // ----------------------- check if forward call was successful and continue then
435  if (destParameters){
436  // successful tries
437  ++m_triesBack;
438  // record the destination parameters
439  ++m_parameters;
440  m_parameterLoc1[m_parameters] = destParameters->parameters()[Trk::loc1];
441  m_parameterLoc2[m_parameters] = destParameters->parameters()[Trk::loc2];
442  m_parameterPhi[m_parameters] = destParameters->parameters()[Trk::phi];
443  m_parameterEta[m_parameters] = destParameters->eta();
444  m_parameterTheta[m_parameters] = destParameters->parameters()[Trk::theta];
445  m_parameterQoverP[m_parameters] = destParameters->parameters()[Trk::qOverP];
446  if(destParameters->covariance()){
447  m_covarianceLoc1[m_parameters] = (*destParameters->covariance())(0,0);
448  m_covarianceLoc2[m_parameters] = (*destParameters->covariance())(1,1);
449  m_covariancePhi[m_parameters] = (*destParameters->covariance())(2,2);
450  m_covarianceTheta[m_parameters] = (*destParameters->covariance())(3,3);
451  m_covarianceQoverP[m_parameters] = (*destParameters->covariance())(4,4);
452  m_covarianceDeterminant[m_parameters] = (destParameters->covariance())->determinant();
453  }
454  else{
455  m_covarianceLoc1[m_parameters] = 0;
456  m_covarianceLoc2[m_parameters] = 0;
457  m_covariancePhi[m_parameters] = 0;
458  m_covarianceTheta[m_parameters] = 0;
459  m_covarianceQoverP[m_parameters] = 0;
460  m_covarianceDeterminant[m_parameters] = 0;
461  }
462  // record the destination parameters
463  const Amg::Vector3D& destinationPosition = destParameters->position();
464  m_destinationX = destinationPosition.x();
465  m_destinationY = destinationPosition.y();
466  m_destinationZ = destinationPosition.z();
467  m_destinationR = destinationPosition.perp();
468 
469  // now simply go backwards
470  const Trk::TrackParameters* backParameters = nullptr;
471  // the standard validation ...
472  if (!m_materialCollectionValidation && !m_direct)
473  backParameters = m_extrapolator->extrapolate(ctx,
474  *destParameters,
475  startSurface,
476  Trk::oppositeMomentum,
477  false,
478  (Trk::ParticleHypothesis)m_particleType,Trk::removeNoise).release();
479  else if(!m_direct){ // material collection validation
480  // get the vector of TrackStateOnSurfaces back
481  const std::vector<const Trk::TrackStateOnSurface*>*
482  collectedBackMaterial = m_extrapolator->extrapolateM(ctx,
483  *destParameters,
484  startSurface,
485  Trk::oppositeMomentum,
486  false,
487  (Trk::ParticleHypothesis)m_particleType);
488  // get the last one and clone it
489  if (collectedBackMaterial && !collectedBackMaterial->empty()){
490  // get the last track state on surface & clone the destination parameters
491  const Trk::TrackStateOnSurface* startState = collectedBackMaterial->back();
492  // assign the last ones of the call
493  backParameters = startState->trackParameters() ? startState->trackParameters()->clone() : nullptr;
494  m_collectedLayerBack += collectedBackMaterial->size();
495  // delete the layers / cleanup
496  std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter = collectedBackMaterial->begin();
497  std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIterEnd = collectedBackMaterial->end();
498  for ( ; tsosIter != tsosIterEnd; delete(*tsosIter), ++tsosIter);
499  }
500  }
501 
502  else{
503  backParameters = m_extrapolator->extrapolateDirectly(ctx,
504  *destParameters,
505  startSurface,
506  Trk::oppositeMomentum,
507  false,
508  (Trk::ParticleHypothesis)m_particleType).release();
509 
510  }
511  // ----------------------- check if backward call was successful and continue then
512  if (backParameters){
513 
514  ATH_MSG_VERBOSE( "Back Parameters : " << *backParameters );
515 
516  // record the back extrapolated ones
517  ++m_parameters;
518  m_parameterLoc1[m_parameters] = backParameters->parameters()[Trk::loc1];
519  m_parameterLoc2[m_parameters] = backParameters->parameters()[Trk::loc2];
520  m_parameterPhi[m_parameters] = backParameters->parameters()[Trk::phi];
521  m_parameterEta[m_parameters] = backParameters->eta();
522  m_parameterTheta[m_parameters] = backParameters->parameters()[Trk::theta];
523  m_parameterQoverP[m_parameters] = backParameters->parameters()[Trk::qOverP];
524  if(backParameters->covariance()){
525  m_covarianceLoc1[m_parameters] = (*backParameters->covariance())(0,0);
526  m_covarianceLoc2[m_parameters] = (*backParameters->covariance())(1,1);
527  m_covariancePhi[m_parameters] = (*backParameters->covariance())(2,2);
528  m_covarianceTheta[m_parameters] = (*backParameters->covariance())(3,3);
529  m_covarianceQoverP[m_parameters] = (*backParameters->covariance())(4,4);
530  m_covarianceDeterminant[m_parameters] = (backParameters->covariance())->determinant();
531  }
532  else{
533  m_covarianceLoc1[m_parameters] = 0;
534  m_covarianceLoc2[m_parameters] = 0;
535  m_covariancePhi[m_parameters] = 0;
536  m_covarianceTheta[m_parameters] = 0;
537  m_covarianceQoverP[m_parameters] = 0;
538  m_covarianceDeterminant[m_parameters] = 0;
539  }
540  // memory cleanup
541  delete backParameters;
542  } else
543  ++m_breaksBack;
544  // memory cleanup
545  delete destParameters;
546  } else
547  ++m_breaksFront;
548  // increase ones more
549  ++m_parameters;
550  // memory cleanup
551 
552  if (m_validationTree)
553  m_validationTree->Fill();
554 
555 
556  //std::cout<<"Cleaning up..."<<std::endl;
557  //delete covariance;
558 
559  return StatusCode::SUCCESS;
560 }
561 
562 //============================================================================================
563 Amg::Transform3D
564 Trk::ExtrapolationValidation::createTransform(double x, double y, double z, double phi, double theta, double alphaZ)
565 {
566 
567  if (phi!=0. && theta != 0.){
568  // create the Start Surface
569  Amg::Vector3D surfacePosition(x,y,z);
570  // z direction
571  Amg::Vector3D surfaceZdirection(cos(phi)*sin(theta),
572  sin(phi)*sin(theta),
573  cos(theta));
574  // the global z axis
575  Amg::Vector3D zAxis(0.,0.,1.);
576  // the y direction
577  Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
578  // the x direction
579  Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
580 
581  double nx = 1./sqrt(surfaceXdirection[0]*surfaceXdirection[0]+surfaceXdirection[1]*surfaceXdirection[1]+surfaceXdirection[2]*surfaceXdirection[2]);
582  double ny = 1./sqrt(surfaceYdirection[0]*surfaceYdirection[0]+surfaceYdirection[1]*surfaceYdirection[1]+surfaceYdirection[2]*surfaceYdirection[2]);
583  surfaceXdirection[0]*=nx;
584  surfaceXdirection[1]*=nx;
585  surfaceXdirection[2]*=nx;
586 
587  surfaceYdirection[0]*=ny;
588  surfaceYdirection[1]*=ny;
589  surfaceYdirection[2]*=ny;
590  // the rotation
591  Amg::RotationMatrix3D surfaceRotation;
592  surfaceRotation.col(0) = surfaceXdirection;
593  surfaceRotation.col(1) = surfaceYdirection;
594  surfaceRotation.col(2) = surfaceZdirection;
595  // return it
596  if (alphaZ==0.)
597  return Amg::Transform3D(surfaceRotation, surfacePosition);
598  Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
599  return Amg::Transform3D(nominalTransform*Amg::AngleAxis3D(alphaZ,zAxis));
600 
601  }
602 
603  return Amg::Transform3D(Amg::Translation3D(x,y,z));
604 }
605 
Trk::y
@ y
Definition: ParamDefs.h:62
Trk::TrackStateOnSurface::trackParameters
const TrackParameters * trackParameters() const
return ptr to trackparameters const overload
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:63
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
PlotCalibFromCool.zAxis
zAxis
Definition: PlotCalibFromCool.py:76
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Trk::oppositeMomentum
@ oppositeMomentum
Definition: PropDirection.h:21
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::loc2
@ loc2
generic first and second local coordinate
Definition: ParamDefs.h:41
IExtrapolator.h
ExtrapolationValidation.h
Trk::alongMomentum
@ alongMomentum
Definition: PropDirection.h:20
Trk::ExtrapolationValidation::~ExtrapolationValidation
~ExtrapolationValidation()
Default Destructor.
Definition: ExtrapolationValidation.cxx:106
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
drawFromPickle.atan
atan
Definition: drawFromPickle.py:36
Trk::ExtrapolationValidation::finalize
StatusCode finalize()
standard Athena-Algorithm method
Definition: ExtrapolationValidation.cxx:182
Trk::ExtrapolationValidation::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: ExtrapolationValidation.cxx:564
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:25
CylinderVolumeBounds.h
Trk::TrackingGeometry::highestTrackingVolume
const TrackingVolume * highestTrackingVolume() const
return the world
python.SystemOfUnits.millimeter
int millimeter
Definition: SystemOfUnits.py:53
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
Trk::theta
@ theta
Definition: ParamDefs.h:72
Trk::TrackingGeometry
Definition: TrackingGeometry.h:67
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::CylinderSurface
Definition: CylinderSurface.h:55
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
CylinderSurface.h
Trk::CylinderVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflengthZ.
Definition: CylinderVolumeBounds.h:207
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ParametersBase
Definition: ParametersBase.h:55
python.SystemOfUnits.micrometer
int micrometer
Definition: SystemOfUnits.py:71
Trk::ExtrapolationValidation::ExtrapolationValidation
ExtrapolationValidation(const std::string &name, ISvcLocator *pSvcLocator)
Standard Athena-Algorithm Constructor.
Definition: ExtrapolationValidation.cxx:27
ParticleHypothesis.h
AthAlgorithm
Definition: AthAlgorithm.h:47
Athena::Units
Definition: Units.h:45
Trk::CylinderVolumeBounds
Definition: CylinderVolumeBounds.h:70
Trk::TrackStateOnSurface
represents the track state (measurement, material, fit parameters and quality) at a surface.
Definition: TrackStateOnSurface.h:71
python.EventInfoMgtInit.release
release
Definition: EventInfoMgtInit.py:24
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
Trk::CylinderVolumeBounds::outerRadius
double outerRadius() const
This method returns the outer radius.
Definition: CylinderVolumeBounds.h:191
charge
double charge(const T &p)
Definition: AtlasPID.h:494
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::ExtrapolationValidation::initialize
StatusCode initialize()
standard Athena-Algorithm method
Definition: ExtrapolationValidation.cxx:116
TrackingVolume.h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
Amg::RotationMatrix3D
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
Definition: GeoPrimitives.h:49
Trk::addNoise
@ addNoise
Definition: MaterialUpdateMode.h:19
PlaneSurface.h
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition: GeoPrimitives.h:44
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:73
DiscSurface.h
Gaudi
=============================================================================
Definition: CaloGPUClusterAndCellDataMonitorOptions.h:273
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
Trk::ParametersBase::eta
double eta() const
Access method for pseudorapidity - from momentum.
TrackingGeometry.h
Amg::AngleAxis3D
Eigen::AngleAxisd AngleAxis3D
Definition: GeoPrimitives.h:45
Trk::phi
@ phi
Definition: ParamDefs.h:81
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Trk::x
@ x
Definition: ParamDefs.h:61
Trk::ExtrapolationValidation::execute
StatusCode execute()
standard Athena-Algorithm method
Definition: ExtrapolationValidation.cxx:204
Trk::loc1
@ loc1
Definition: ParamDefs.h:40
Trk::removeNoise
@ removeNoise
Definition: MaterialUpdateMode.h:20
GeV
#define GeV
Definition: CaloTransverseBalanceVecMon.cxx:30
TrackStateOnSurface.h
Trk::ParametersBase::clone
virtual ParametersBase< DIM, T > * clone() const override=0
clone method for polymorphic deep copy
Generated on Tue Jul 2 2024 21:15:33 for ATLAS Offline Software by   doxygen 1.8.18