ATLAS Offline Software
ExtrapolatorComparisonTest.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 // ExtrapolatorComparisonTest.cxx, (c) ATLAS Detector software
8 
9 #include "TrkExAlgs/ExtrapolatorComparisonTest.h"
10 
11 // Tracking
12 #include "TrkSurfaces/PerigeeSurface.h"
13 #include "TrkSurfaces/CylinderSurface.h"
14 #include "TrkSurfaces/DiscSurface.h"
15 #include "TrkEventPrimitives/ParticleHypothesis.h"
16 #include "TrkParametersBase/ParametersT.h"
17 #include "EventPrimitives/EventPrimitives.h"
18 
19 #include "GaudiKernel/ISvcLocator.h"
20 
21 // ACTS
22 #include "Acts/Definitions/Algebra.hpp"
23 #include "Acts/Definitions/TrackParametrization.hpp"
24 #include "Acts/EventData/ParticleHypothesis.hpp"
25 #include "Acts/Geometry/GeometryContext.hpp"
26 #include "Acts/Propagator/detail/SteppingLogger.hpp"
27 #include "Acts/Surfaces/CylinderSurface.hpp"
28 #include "Acts/Surfaces/DiscSurface.hpp"
29 #include "Acts/Surfaces/PerigeeSurface.hpp"
30 #include "Acts/Surfaces/Surface.hpp"
31 #include "Acts/Utilities/Helpers.hpp"
32 #include "Acts/Utilities/Logger.hpp"
33 #include "ActsGeometryInterfaces/ITrackingGeometryTool.h"
34 #include "ActsInterop/Logger.h"
35 
36 // OTHER
37 #include "CLHEP/Random/RandomEngine.h"
38 #include "CLHEP/Random/Randomize.h"
39 
40 // STL
41 #include <string>
42 #include <fstream>
43 #include <chrono>
44 #include <cmath>
45 
46 using xclock = std::chrono::steady_clock;
47 
48 //================ Constructor =================================================
49 
50 Trk::ExtrapolatorComparisonTest::ExtrapolatorComparisonTest(const std::string& name, ISvcLocator* pSvcLocator) :
51  AthReentrantAlgorithm(name,pSvcLocator) {}
52 
53 //================ Destructor =================================================
54 
55 Trk::ExtrapolatorComparisonTest::~ExtrapolatorComparisonTest()
56 {
57  // cleanup of the Trk::Surfaces
58  for (const auto& surfaceTriple : m_atlasReferenceSurfaceTriples) {
59  for (const auto* surface : surfaceTriple) {
60  delete surface;
61  }
62  }
63 }
64 
65 
66 //================ Initialisation =================================================
67 
68 StatusCode Trk::ExtrapolatorComparisonTest::initialize()
69 {
70  // Code entered here will be executed once at program start.
71 
72  ATH_MSG_INFO(" initialize()");
73 
74  ATH_CHECK( m_extrapolationTool.retrieve() );
75  ATH_CHECK( m_atlasExtrapolator.retrieve() );
76  ATH_CHECK( m_trackingGeometryTool.retrieve() );
77 
78  // Create the destination surfaces for extrapolation
79  // --> you need the Trk::Surfaces and the Acts::Surfaces
80  if (m_referenceSurfaceRadius.size() == m_referenceSurfaceHalflength.size()) {
81  // assign the size
82  m_referenceSurfaces = m_referenceSurfaceRadius.size();
83  // loop over it and create the
84  for (unsigned int surface = 0; surface < m_referenceSurfaces; surface++) {
85 
86  double radius = m_referenceSurfaceRadius[surface];
87  double halfZ = m_referenceSurfaceHalflength[surface];
88 
89  // create the Surface triplet
90  std::vector< const Trk::Surface*> trkSurfaceTriplet;
91  trkSurfaceTriplet.push_back(new Trk::DiscSurface (Amg::Transform3D(Amg::Translation3D(0.,0., halfZ)), 0.,radius));
92  trkSurfaceTriplet.push_back(new Trk::CylinderSurface(Amg::Transform3D(Amg::Translation3D(0.,0., 0.)),radius, halfZ));
93  trkSurfaceTriplet.push_back(new Trk::DiscSurface (Amg::Transform3D(Amg::Translation3D(0.,0.,-halfZ)), 0.,radius));
94  ATH_MSG_INFO("Creating Trk::Surface at: R " << radius << " Z " << halfZ);
95  m_atlasReferenceSurfaceTriples.push_back(trkSurfaceTriplet);
96 
97  // create the Surface triplet
98  std::vector<std::shared_ptr<const Acts::Surface>> actsSurfaceTriplet;
99 
100  Acts::Transform3 posTransf(Acts::Transform3::Identity()*Acts::Translation3(Acts::Vector3(0.,0., halfZ)));
101  Acts::Transform3 cTransf(Acts::Transform3::Identity()*Acts::Translation3(Acts::Vector3(0.,0., 0.)));
102  Acts::Transform3 negTransf(Acts::Transform3::Identity()*Acts::Translation3(Acts::Vector3(0.,0.,-halfZ)));
103 
104  auto posSurface = Acts::Surface::makeShared<Acts::DiscSurface> (posTransf , 0.,radius);
105  auto cSurface = Acts::Surface::makeShared<Acts::CylinderSurface>( cTransf ,radius, halfZ);
106  auto negSurface = Acts::Surface::makeShared<Acts::DiscSurface> (negTransf , 0.,radius);
107 
108  actsSurfaceTriplet.push_back(posSurface);
109  actsSurfaceTriplet.push_back(cSurface );
110  actsSurfaceTriplet.push_back(negSurface);
111  ATH_MSG_INFO("Creating Acts::Surface at: R " << radius << " Z " << halfZ);
112  m_actsReferenceSurfaceTriples.push_back(actsSurfaceTriplet);
113 
114  m_referenceSurfaceNegativeBoundary.push_back(atan2(radius,-halfZ));
115  m_referenceSurfacePositiveBoundary.push_back(atan2(radius, halfZ));
116  }
117  } else {
118  ATH_MSG_WARNING("Not compatible size of ReferenceSurfaceRadius and ReferenceSurfaceHalfZ!! Returning FAILURE!");
119  return StatusCode::FAILURE;
120  }
121 
122  ATH_CHECK( m_atlasPropResultWriterSvc.retrieve() );
123  ATH_CHECK( m_actsPropResultWriterSvc.retrieve() );
124  ATH_CHECK( m_rndmSvc.retrieve() );
125  m_randomEngine = m_rndmSvc->getEngine (this, "ExtrapolatorComparisonTest");
126 
127  return StatusCode::SUCCESS;
128 }
129 
130 //================ Finalisation =================================================
131 
132 StatusCode Trk::ExtrapolatorComparisonTest::finalize()
133 {
134  // Code entered here will be executed once at the end of the program run.
135  return StatusCode::SUCCESS;
136 }
137 
138 //================ Execution ====================================================
139 
140 StatusCode Trk::ExtrapolatorComparisonTest::execute(const EventContext& ctx) const {
141 
142  float milliseconds_to_seconds = 1000.;
143 
144  // generate perigees with random number generator
145  CLHEP::HepRandomEngine* engine = m_randomEngine->getEngine(ctx);
146 
147  std::vector<perigeeParameters> parameters = {};
148  for (int ext = 0; ext<m_eventsPerExecute; ext++) {
149  // generate with random number generator
150  double d0 = CLHEP::RandGauss::shoot(engine) * m_sigmaD0;
151  double z0 = CLHEP::RandGauss::shoot(engine) * m_sigmaZ0;
152  double phi = m_minPhi + (m_maxPhi-m_minPhi)* CLHEP::RandFlat::shoot(engine);
153  double eta = m_minEta + CLHEP::RandFlat::shoot(engine)*(m_maxEta-m_minEta);
154  double pt = m_minPt + CLHEP::RandFlat::shoot(engine)*(m_maxPt-m_minPt);
155  double charge = (CLHEP::RandFlat::shoot(engine) > 0.5 ) ? -1. : 1.;
156  parameters.emplace_back(d0, z0, phi, eta, pt, charge);
157  }
158 
159  int n_extraps = 0;
160  auto start = xclock::now();
161  for (auto& perigee : parameters) {
162 
163  Acts::Vector3 momentum(perigee.m_pt * std::cos(perigee.m_phi), perigee.m_pt * std::sin(perigee.m_phi), perigee.m_pt * std::sinh(perigee.m_eta));
164  double theta = Acts::VectorHelpers::theta(momentum);
165  double qOverP = perigee.m_charge / momentum.norm();
166 
167  const Trk::PerigeeSurface atlPerigeeSurface;
168  const Trk::Perigee * atlPerigee = new const Trk::Perigee(perigee.m_d0, perigee.m_z0, perigee.m_phi, theta, qOverP, atlPerigeeSurface);
169 
170  for (unsigned int surface = 0; surface < m_atlasReferenceSurfaceTriples.size(); surface++) {
171  n_extraps++;
172 
173  double negRef = m_referenceSurfaceNegativeBoundary.at(surface);
174  double posRef = m_referenceSurfacePositiveBoundary.at(surface);
175 
176  // decide which reference surface to take
177  int refSurface = theta < posRef ? 2 : 1;
178  refSurface = theta > negRef ? 0 : 1;
179 
180  const Trk::Surface* destinationSurface = m_atlasReferenceSurfaceTriples.at(surface).at(refSurface);
181 
182  ATH_MSG_VERBOSE("Starting extrapolation " << n_extraps << " from : " << *atlPerigee << " to : " << *destinationSurface);
183 
184  auto start_fwd = xclock::now();
185  const Trk::TrackParameters* destParameters =
186  m_atlasExtrapolator->extrapolate(
187  ctx,
188  *atlPerigee,
189  *destinationSurface,
190  Trk::alongMomentum,
191  true,
192  static_cast<Trk::ParticleHypothesis>(m_particleType.value())).release();
193  auto end_fwd = xclock::now();
194  float ms_fwd = std::chrono::duration_cast<std::chrono::milliseconds>(end_fwd-start_fwd).count();
195 
196  if (destParameters) {
197  ATH_MSG_VERBOSE(" ATLAS Extrapolator succeded!! --> Forward" );
198  ATH_MSG_VERBOSE(" [ intersection ] with surface at (x,y,z) = " << destParameters->position().x() << ", " << destParameters->position().y() << ", " << destParameters->position().z() );
199  ATH_MSG_VERBOSE(" [ intersection ] parameters: " << destParameters->parameters() );
200  ATH_MSG_VERBOSE(" [ intersection ] cov matrix: " << destParameters->covariance() );
201 
202  // now try backward extrapolation
203  auto start_bkw = xclock::now();
204  const Trk::TrackParameters* finalperigee =
205  m_atlasExtrapolator->extrapolate(
206  ctx,
207  *destParameters,
208  atlPerigee->associatedSurface(),
209  Trk::oppositeMomentum,
210  true,
211  static_cast<Trk::ParticleHypothesis>(m_particleType.value())).release();
212  auto end_bkw = xclock::now();
213  float ms_bkw = std::chrono::duration_cast<std::chrono::milliseconds>(end_bkw-start_bkw).count();
214 
215  if (finalperigee) {
216  ATH_MSG_VERBOSE(" ATLAS Extrapolator succeded!! --> Backward" );
217  ATH_MSG_VERBOSE(" [extrapolation to perigee] input: " << atlPerigee->parameters() );
218  ATH_MSG_VERBOSE(" [extrapolation to perigee] output: " << finalperigee->parameters() );
219  ATH_MSG_VERBOSE(" [extrapolation to perigee] cov matrix: " << finalperigee->covariance() );
220 
221  } else if (!finalperigee) {
222  ATH_MSG_DEBUG(" ATLAS Extrapolation to perigee failed for input parameters: " << destParameters->parameters());
223  }
224 
225  m_atlasPropResultWriterSvc->write<Trk::TrackParameters>(atlPerigee, destParameters, ms_fwd, finalperigee, ms_bkw);
226  delete finalperigee;
227  } else if (!destParameters) {
228  ATH_MSG_DEBUG(" ATLAS Extrapolation not successful! " );
229  m_atlasPropResultWriterSvc->write<Trk::TrackParameters>(atlPerigee);
230  }
231  delete destParameters;
232  }
233  delete atlPerigee;
234  }
235  auto end = xclock::now();
236  auto secs = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count() / milliseconds_to_seconds;
237  double secs_per_ex = secs / n_extraps;
238 
239  ATH_MSG_INFO("ATLAS : Time for " << n_extraps << " iterations: " << secs << "s (" << secs_per_ex << "s per extrapolation)");
240 
241  n_extraps = 0;
242  start = xclock::now();
243  for (auto& perigee : parameters) {
244  Acts::Vector3 momentum(perigee.m_pt * std::cos(perigee.m_phi), perigee.m_pt * std::sin(perigee.m_phi), perigee.m_pt * std::sinh(perigee.m_eta));
245  double theta = Acts::VectorHelpers::theta(momentum);
246  double qOverP = perigee.m_charge / momentum.norm();
247 
248  std::shared_ptr<Acts::PerigeeSurface> actsPerigeeSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(Acts::Vector3(0, 0, 0));
249  double t = 0.;
250  Acts::BoundVector pars;
251  //cppcheck-suppress constStatement
252  pars << perigee.m_d0, perigee.m_z0, perigee.m_phi, theta, qOverP, t;
253  std::optional<Acts::BoundSquareMatrix> cov = std::nullopt;
254 
255  // Perigee, no alignment -> default geo context
256  ActsGeometryContext gctx = m_trackingGeometryTool->getNominalGeometryContext();
257  auto anygctx = gctx.context();
258  const auto* startParameters = new const Acts::GenericBoundTrackParameters(std::move(actsPerigeeSurface), pars, std::move(cov), Acts::ParticleHypothesis::pion());
259 
260  for (unsigned int surface = 0; surface < m_actsReferenceSurfaceTriples.size(); surface++) {
261  n_extraps++;
262 
263  double negRef = m_referenceSurfaceNegativeBoundary.at(surface);
264  double posRef = m_referenceSurfacePositiveBoundary.at(surface);
265 
266  // decide which reference surface to take
267  int refSurface = theta < posRef ? 2 : 1;
268  refSurface = theta > negRef ? 0 : 1;
269 
270  auto destinationSurface = m_actsReferenceSurfaceTriples.at(surface).at(refSurface);
271 
272  ATH_MSG_VERBOSE("Starting extrapolation " << n_extraps << " from : " << pars << " to : " << destinationSurface);
273 
274  auto start_fwd = xclock::now();
275  auto destParameters = m_extrapolationTool->propagate(ctx, *startParameters, *destinationSurface, Acts::Direction::Forward());
276  auto end_fwd = xclock::now();
277  float ms_fwd = std::chrono::duration_cast<std::chrono::milliseconds>(end_fwd-start_fwd).count();
278 
279  if (destParameters) {
280  ATH_MSG_VERBOSE(" ACTS Extrapolator succeded!! --> Forward" );
281  ATH_MSG_VERBOSE(" [ intersection ] with surface at (x,y,z) = " << destParameters->position(anygctx).x() << ", " << destParameters->position(anygctx).y() << ", " << destParameters->position(anygctx).z() );
282  ATH_MSG_VERBOSE(" [ intersection ] parameters: " << destParameters->parameters() );
283  ATH_MSG_VERBOSE(" [ intersection ] cov matrix: " << *destParameters->covariance() );
284 
285  // now try backward extrapolation
286  auto start_bkw = xclock::now();
287  auto finalperigee = m_extrapolationTool->propagate(ctx, *destParameters, startParameters->referenceSurface(), Acts::Direction::Backward());
288  auto end_bkw = xclock::now();
289  float ms_bkw = std::chrono::duration_cast<std::chrono::milliseconds>(end_bkw-start_bkw).count();
290 
291  if (finalperigee) {
292  ATH_MSG_VERBOSE(" ACTS Extrapolator succeded!! --> Backward" );
293  ATH_MSG_VERBOSE(" [extrapolation to perigee] input: " << startParameters->parameters() );
294  ATH_MSG_VERBOSE(" [extrapolation to perigee] output: " << finalperigee->parameters() );
295  ATH_MSG_VERBOSE(" [extrapolation to perigee] cov matrix: " << *finalperigee->covariance() );
296 
297  } else if (!finalperigee) {
298  ATH_MSG_DEBUG(" ACTS Extrapolation to perigee failed for input parameters: " << destParameters->parameters());
299  }
300 
301  // Construct wrappers for Acts track parameters
302  const ActsTrackWrapper* startWrapper = new ActsTrackWrapper(startParameters, anygctx);
303  const ActsTrackWrapper* destWrapper = new ActsTrackWrapper(&destParameters.value(), anygctx);
304  const ActsTrackWrapper* finalWrapper = new ActsTrackWrapper(&finalperigee.value(), anygctx);
305 
306  m_actsPropResultWriterSvc->write<ActsTrackWrapper>(startWrapper, destWrapper, ms_fwd, finalWrapper, ms_bkw);
307 
308  delete startWrapper;
309  delete destWrapper;
310  delete finalWrapper;
311  } else if (!destParameters) {
312  ATH_MSG_DEBUG(" ACTS Extrapolation not successful! " );
313  const ActsTrackWrapper* startWrapper = new ActsTrackWrapper(startParameters, anygctx);
314  m_actsPropResultWriterSvc->write<ActsTrackWrapper>(startWrapper);
315  delete startWrapper;
316  }
317  }
318  delete startParameters;
319  }
320 
321  end = xclock::now();
322  secs = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count() / milliseconds_to_seconds;
323  secs_per_ex = secs / n_extraps;
324 
325  ATH_MSG_INFO("ACTS : Time for " << n_extraps << " iterations: " << secs << "s (" << secs_per_ex << "s per extrapolation)");
326 
327  return StatusCode::SUCCESS;
328 }
329 
330 //============================================================================================
331 
332 
333 
make_hlt_rep.pars
pars
Definition: make_hlt_rep.py:90
Trk::ExtrapolatorComparisonTest::ActsTrackWrapper
Definition: ExtrapolatorComparisonTest.h:94
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Trk::ExtrapolatorComparisonTest::~ExtrapolatorComparisonTest
~ExtrapolatorComparisonTest()
Default Destructor.
Definition: ExtrapolatorComparisonTest.cxx:55
PerigeeSurface.h
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
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
mergePhysValFiles.start
start
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:13
Trk::ExtrapolatorComparisonTest::ExtrapolatorComparisonTest
ExtrapolatorComparisonTest(const std::string &name, ISvcLocator *pSvcLocator)
Standard Athena-Algorithm Constructor.
Definition: ExtrapolatorComparisonTest.cxx:50
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:200
test_pyathena.pt
pt
Definition: test_pyathena.py:11
Trk::ExtrapolatorComparisonTest::finalize
StatusCode finalize() override
standard Athena-Algorithm method
Definition: ExtrapolatorComparisonTest.cxx:132
xAOD::pion
@ pion
Definition: TrackingPrimitives.h:197
Trk::z0
@ z0
Definition: ParamDefs.h:64
Trk::alongMomentum
@ alongMomentum
Definition: PropDirection.h:20
Trk::DiscSurface
Definition: DiscSurface.h:54
ExtrapolatorComparisonTest.h
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::Perigee
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
Definition: Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:33
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:92
AthReentrantAlgorithm
An algorithm that can be simultaneously executed in multiple threads.
Definition: AthReentrantAlgorithm.h:74
ActsGeometryContext::context
Acts::GeometryContext context() const
Definition: ActsGeometryContext.h:45
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:28
Trk::ParametersT::associatedSurface
virtual const S & associatedSurface() const override final
Access to the Surface method.
python.handimod.now
now
Definition: handimod.py:674
Trk::ExtrapolatorComparisonTest::initialize
StatusCode initialize() override
standard Athena-Algorithm method
Definition: ExtrapolatorComparisonTest.cxx:68
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::CylinderSurface
Definition: CylinderSurface.h:55
ParametersT.h
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
CylinderSurface.h
Trk::ExtrapolatorComparisonTest::execute
StatusCode execute(const EventContext &ctx) const override
standard Athena-Algorithm method
Definition: ExtrapolatorComparisonTest.cxx:140
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
MakeFileForMJB.ext
string ext
Definition: Moriond2016/MakeFileForMJB.py:41
Trk::ParametersBase
Definition: ParametersBase.h:55
ParticleHypothesis.h
ActsGeometryContext
Include the GeoPrimitives which need to be put first.
Definition: ActsGeometryContext.h:27
EventPrimitives.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Trk::d0
@ d0
Definition: ParamDefs.h:63
charge
double charge(const T &p)
Definition: AtlasPID.h:986
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition: GeoPrimitives.h:44
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
ITrackingGeometryTool.h
DiscSurface.h
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Trk::phi
@ phi
Definition: ParamDefs.h:75
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Logger.h
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
xclock
std::chrono::steady_clock xclock
Definition: ExtrapolatorComparisonTest.cxx:46
Generated on Mon Sep 1 2025 21:10:17 for ATLAS Offline Software by   doxygen 1.8.18