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