ATLAS Offline Software
GeantFollowerMSHelper.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 // GeantFollowerMSHelper.cxx, (c) ATLAS Detector Software
8 
9 // StoreGate
10 #include "TrkG4UserActions/GeantFollowerMSHelper.h"
11 
12 #include "GaudiKernel/ITHistSvc.h"
13 #include "StoreGate/StoreGateSvc.h"
14 #include "TTree.h"
15 // Trk
16 #include <cmath>
17 
18 #include "GeoPrimitives/GeoPrimitives.h"
19 #include "TrkExInterfaces/IExtrapolator.h"
20 #include "TrkExUtils/ExtrapolationCache.h"
21 #include "TrkMaterialOnTrack/EnergyLoss.h"
22 #include "TrkMaterialOnTrack/MaterialEffectsOnTrack.h"
23 #include "TrkMaterialOnTrack/ScatteringAngles.h"
24 #include "TrkParameters/TrackParameters.h"
25 #include "TrkSurfaces/PlaneSurface.h"
26 
27 Trk::GeantFollowerMSHelper::GeantFollowerMSHelper(const std::string& t,
28  const std::string& n,
29  const IInterface* p)
30  : base_class(t, n, p),
31  m_extrapolator(""),
32  m_extrapolateDirectly(false),
33  m_extrapolateIncrementally(false),
34  m_speedup(false),
35  m_useCovMatrix(true),
36  m_useIDExit(false),
37  m_parameterCache(nullptr),
38  m_parameterCacheCov(nullptr),
39  m_parameterCacheEntry(nullptr),
40  m_parameterCacheEntryCov(nullptr),
41  m_tX0Cache(0.),
42  m_crossedEntry(false),
43  m_exitLayer(false),
44  m_destinationSurface(),
45  m_validationTreeName("G4Follower"),
46  m_validationTreeDescription("Output of the G4Follower_"),
47  m_validationTreeFolder("/val/G4Follower"),
48  m_validationTree(nullptr) {
49  declareProperty("Extrapolator", m_extrapolator);
50  declareProperty("ExtrapolateDirectly", m_extrapolateDirectly);
51  declareProperty("ExtrapolateIncrementally", m_extrapolateIncrementally);
52 
53  // SpeedUp False takes more CPU because it will stop at each G4 Step in the
54  declareProperty("SpeedUp", m_speedup);
55  declareProperty("UseCovMatrix", m_useCovMatrix);
56  declareProperty("UseIDExit",m_useIDExit); // used for validating the ID and Calo tracking
57 
58 
59 }
60 
61 // destructor
62 Trk::GeantFollowerMSHelper::~GeantFollowerMSHelper() = default;
63 
64 // Athena standard methods
65 // initialize
66 StatusCode Trk::GeantFollowerMSHelper::initialize() {
67  m_treeData = std::make_unique<TreeData>();
68 
69 
70 
71 
72 
73  if (m_extrapolator.retrieve().isFailure()) {
74  ATH_MSG_ERROR("Could not retrieve Extrapolator " << m_extrapolator
75  << " . Abort.");
76  return StatusCode::FAILURE;
77  }
78 
79  if (m_elossupdator.retrieve().isFailure()) {
80  ATH_MSG_ERROR("Could not retrieve ElossUpdator " << m_elossupdator
81  << " . Abort.");
82  return StatusCode::FAILURE;
83  }
84 
85 
86  if (m_speedup) {
87  ATH_MSG_INFO(" SpeedUp GeantFollowerMS ");
88  } else {
89  ATH_MSG_INFO(" NO SpeedUp GeantFollowerMS ");
90  }
91  ATH_MSG_INFO("initialize()");
92 
93  // create the new Tree
94  m_validationTree = new TTree(m_validationTreeName.c_str(),
95  m_validationTreeDescription.c_str());
96 
97  m_validationTree->Branch("InitX", &m_treeData->m_t_x, "initX/F");
98  m_validationTree->Branch("InitY", &m_treeData->m_t_y, "initY/F");
99  m_validationTree->Branch("InitZ", &m_treeData->m_t_z, "initZ/F");
100  m_validationTree->Branch("InitTheta", &m_treeData->m_t_theta, "initTheta/F");
101  m_validationTree->Branch("InitEta", &m_treeData->m_t_eta, "initEta/F");
102  m_validationTree->Branch("InitPhi", &m_treeData->m_t_phi, "initPhi/F");
103  m_validationTree->Branch("InitP", &m_treeData->m_t_p, "initP/F");
104  m_validationTree->Branch("InitPdg", &m_treeData->m_t_pdg, "initPdg/I");
105  m_validationTree->Branch("InitCharge", &m_treeData->m_t_charge, "initQ/F");
106 
107  m_validationTree->Branch("MEntryX", &m_treeData->m_m_x, "mentryX/F");
108  m_validationTree->Branch("MEntryY", &m_treeData->m_m_y, "mentryY/F");
109  m_validationTree->Branch("MEntryZ", &m_treeData->m_m_z, "mentryZ/F");
110  m_validationTree->Branch("MEntryTheta", &m_treeData->m_m_theta,
111  "mentryTheta/F");
112  m_validationTree->Branch("MEntryEta", &m_treeData->m_m_eta, "mentryEta/F");
113  m_validationTree->Branch("MEntryPhi", &m_treeData->m_m_phi, "mentryPhi/F");
114  m_validationTree->Branch("MEntryP", &m_treeData->m_m_p, "mentryP/F");
115 
116  m_validationTree->Branch("BackX", &m_treeData->m_b_x, "backX/F");
117  m_validationTree->Branch("BackY", &m_treeData->m_b_y, "backY/F");
118  m_validationTree->Branch("BackZ", &m_treeData->m_b_z, "backZ/F");
119  m_validationTree->Branch("BackTheta", &m_treeData->m_b_theta, "backTheta/F");
120  m_validationTree->Branch("BackEta", &m_treeData->m_b_eta, "backEta/F");
121  m_validationTree->Branch("BackPhi", &m_treeData->m_b_phi, "backPhi/F");
122  m_validationTree->Branch("BackP", &m_treeData->m_b_p, "backP/F");
123  m_validationTree->Branch("BackX0", &m_treeData->m_b_X0, "backX0/F");
124  m_validationTree->Branch("BackEloss", &m_treeData->m_b_Eloss, "backEloss/F");
125 
126  m_validationTree->Branch("G4Steps", &m_treeData->m_g4_steps, "g4steps/I");
127  m_validationTree->Branch("TrkStepScats", &m_treeData->m_trk_scats,
128  "trkscats/I");
129 
130  m_validationTree->Branch("G4StepP", m_treeData->m_g4_p, "g4stepP[g4steps]/F");
131  m_validationTree->Branch("G4StepEta", m_treeData->m_g4_eta,
132  "g4stepEta[g4steps]/F");
133  m_validationTree->Branch("G4StepTheta", m_treeData->m_g4_theta,
134  "g4stepTheta[g4steps]/F");
135  m_validationTree->Branch("G4StepPhi", m_treeData->m_g4_phi,
136  "g4stepPhi[g4steps]/F");
137  m_validationTree->Branch("G4StepX", m_treeData->m_g4_x, "g4stepX[g4steps]/F");
138  m_validationTree->Branch("G4StepY", m_treeData->m_g4_y, "g4stepY[g4steps]/F");
139  m_validationTree->Branch("G4StepZ", m_treeData->m_g4_z, "g4stepZ[g4steps]/F");
140  m_validationTree->Branch("G4AccumTX0", m_treeData->m_g4_tX0,
141  "g4stepAccTX0[g4steps]/F");
142  m_validationTree->Branch("G4StepT", m_treeData->m_g4_t,
143  "g4stepTX[g4steps]/F");
144  m_validationTree->Branch("G4StepX0", m_treeData->m_g4_X0,
145  "g4stepX0[g4steps]/F");
146 
147  m_validationTree->Branch("TrkStepStatus", m_treeData->m_trk_status,
148  "trkstepStatus[g4steps]/I");
149  m_validationTree->Branch("TrkStepP", m_treeData->m_trk_p,
150  "trkstepP[g4steps]/F");
151  m_validationTree->Branch("TrkStepEta", m_treeData->m_trk_eta,
152  "trkstepEta[g4steps]/F");
153  m_validationTree->Branch("TrkStepTheta", m_treeData->m_trk_theta,
154  "trkstepTheta[g4steps]/F");
155  m_validationTree->Branch("TrkStepPhi", m_treeData->m_trk_phi,
156  "trkstepPhi[g4steps]/F");
157  m_validationTree->Branch("TrkStepX", m_treeData->m_trk_x,
158  "trkstepX[g4steps]/F");
159  m_validationTree->Branch("TrkStepY", m_treeData->m_trk_y,
160  "trkstepY[g4steps]/F");
161  m_validationTree->Branch("TrkStepZ", m_treeData->m_trk_z,
162  "trkstepZ[g4steps]/F");
163  m_validationTree->Branch("TrkStepLocX", m_treeData->m_trk_lx,
164  "trkstepLX[g4steps]/F");
165  m_validationTree->Branch("TrkStepLocY", m_treeData->m_trk_ly,
166  "trkstepLY[g4steps]/F");
167  m_validationTree->Branch("TrkStepEloss", m_treeData->m_trk_eloss,
168  "trkstepEloss[g4steps]/F");
169  m_validationTree->Branch("TrkStepEloss1", m_treeData->m_trk_eloss1,
170  "trkstepEloss1[g4steps]/F");
171  m_validationTree->Branch("TrkStepEloss0", m_treeData->m_trk_eloss0,
172  "trkstepEloss0[g4steps]/F");
173  m_validationTree->Branch("TrkStepEloss5", m_treeData->m_trk_eloss5,
174  "trkstepEloss5[g4steps]/F");
175  m_validationTree->Branch("TrkStepEloss10", m_treeData->m_trk_eloss10,
176  "trkstepEloss10[g4steps]/F");
177  m_validationTree->Branch("TrkStepScaleEloss", m_treeData->m_trk_scaleeloss,
178  "trkstepScaleEloss[g4steps]/F");
179  m_validationTree->Branch("TrkStepScaleX0", m_treeData->m_trk_scalex0,
180  "trkstepScaleX0[g4steps]/F");
181  m_validationTree->Branch("TrkStepX0", m_treeData->m_trk_x0,
182  "trkstepX0[g4steps]/F");
183  m_validationTree->Branch("TrkStepErd0", m_treeData->m_trk_erd0,
184  "trkstepErd0[g4steps]/F");
185  m_validationTree->Branch("TrkStepErz0", m_treeData->m_trk_erz0,
186  "trkstepErz0[g4steps]/F");
187  m_validationTree->Branch("TrkStepErphi", m_treeData->m_trk_erphi,
188  "trkstepErphi[g4steps]/F");
189  m_validationTree->Branch("TrkStepErtheta", m_treeData->m_trk_ertheta,
190  "trkstepErtheta[g4steps]/F");
191  m_validationTree->Branch("TrkStepErqoverp", m_treeData->m_trk_erqoverp,
192  "trkstepErqoverp[g4steps]/F");
193 
194  m_validationTree->Branch("TrkStepScatStatus", m_treeData->m_trk_sstatus,
195  "trkscatStatus[trkscats]/I");
196  m_validationTree->Branch("TrkStepScatX", m_treeData->m_trk_sx,
197  "trkscatX[trkscats]/F");
198  m_validationTree->Branch("TrkStepScatY", m_treeData->m_trk_sy,
199  "trkscatY[trkscats]/F");
200  m_validationTree->Branch("TrkStepScatZ", m_treeData->m_trk_sz,
201  "trkscatZ[trkscats]/F");
202  m_validationTree->Branch("TrkStepScatX0", m_treeData->m_trk_sx0,
203  "trkscatX0[trkscats]/F");
204  m_validationTree->Branch("TrkStepScatEloss", m_treeData->m_trk_seloss,
205  "trkscatEloss[trkscats]/F");
206  m_validationTree->Branch("TrkStepScatMeanIoni", m_treeData->m_trk_smeanIoni,
207  "trkscatMeanIoni[trkscats]/F");
208  m_validationTree->Branch("TrkStepScatSigIoni", m_treeData->m_trk_ssigIoni,
209  "trkscatSigIoni[trkscats]/F");
210  m_validationTree->Branch("TrkStepScatMeanRad", m_treeData->m_trk_smeanRad,
211  "trkscatMeanRad[trkscats]/F");
212  m_validationTree->Branch("TrkStepScatSigRad", m_treeData->m_trk_ssigRad,
213  "trkscatSigRad[trkscats]/F");
214  m_validationTree->Branch("TrkStepScatSigTheta", m_treeData->m_trk_ssigTheta,
215  "trkscatSigTheta[trkscats]/F");
216  m_validationTree->Branch("TrkStepScatSigPhi", m_treeData->m_trk_ssigPhi,
217  "trkscatSigPhi[trkscats]/F");
218 
219  m_crossedEntry = false;
220  m_exitLayer = false;
221  // now register the Tree
222  SmartIF<ITHistSvc> tHistSvc{Gaudi::svcLocator()->service("THistSvc")};
223  if ( !tHistSvc ) {
224  ATH_MSG_ERROR(
225  "Could not find Hist Service -> Switching ValidationMode Off !");
226  delete m_validationTree;
227  m_validationTree = nullptr;
228  }
229  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree))
230  .isFailure()) {
231  ATH_MSG_ERROR(
232  "Could not register the validation Tree -> Switching ValidationMode "
233  "Off !");
234  delete m_validationTree;
235  m_validationTree = nullptr;
236  }
237 
238  ATH_MSG_INFO("initialize() successful");
239  return StatusCode::SUCCESS;
240 }
241 
242 StatusCode Trk::GeantFollowerMSHelper::finalize() {
243  return StatusCode::SUCCESS;
244 }
245 
246 void Trk::GeantFollowerMSHelper::beginEvent() {
247  m_treeData->m_t_x = 0.;
248  m_treeData->m_t_y = 0.;
249  m_treeData->m_t_z = 0.;
250  m_treeData->m_t_theta = 0.;
251  m_treeData->m_t_eta = 0.;
252  m_treeData->m_t_phi = 0.;
253  m_treeData->m_t_p = 0.;
254  m_treeData->m_t_charge = 0.;
255  m_treeData->m_t_pdg = 0;
256 
257  m_treeData->m_m_x = 0.;
258  m_treeData->m_m_y = 0.;
259  m_treeData->m_m_z = 0.;
260  m_treeData->m_m_theta = 0.;
261  m_treeData->m_m_eta = 0.;
262  m_treeData->m_m_phi = 0.;
263  m_treeData->m_m_p = 0.;
264 
265  m_treeData->m_b_x = 0.;
266  m_treeData->m_b_y = 0.;
267  m_treeData->m_b_z = 0.;
268  m_treeData->m_b_theta = 0.;
269  m_treeData->m_b_eta = 0.;
270  m_treeData->m_b_phi = 0.;
271  m_treeData->m_b_p = 0.;
272  m_treeData->m_b_X0 = 0.;
273  m_treeData->m_b_Eloss = 0.;
274 
275  m_treeData->m_g4_steps = -1;
276  m_treeData->m_g4_stepsEntry = -1;
277  m_treeData->m_trk_scats = 0;
278  m_tX0Cache = 0.;
279 
280  m_crossedEntry = false;
281  m_exitLayer = false;
282 }
283 
284 void Trk::GeantFollowerMSHelper::trackParticle(const G4ThreeVector& pos,
285  const G4ThreeVector& mom,
286  int pdg, double charge, float t,
287  float X0) {
288  const EventContext& ctx = Gaudi::Hive::currentContext();
289  // as the MS starts at 6736 in R.07 the cut is just before
290 
291  double zMuonEntry = 6735.;
292  double rMuonEntry = 4254;
293  double zMuonExit = 21800.;
294  double rMuonExit = 12500.;
295  double zIDExit = 2720.;
296  double rIDExit = 1080.;
297 
298  double zEntry = zMuonEntry;
299  double rEntry = rMuonEntry;
300  double zExit = zMuonExit;
301  double rExit = rMuonExit;
302 
303  if(m_useIDExit) {
304  zEntry = zIDExit;
305  rEntry = rIDExit;
306  zExit = zMuonEntry;
307  rExit = rMuonEntry;
308  }
309 
310  double scale = 1.;
311 
312  Amg::Vector3D npos(scale * pos.x(), scale * pos.y(), scale * pos.z());
313  Amg::Vector3D nmom(mom.x(), mom.y(), mom.z());
314 
315  if (m_treeData->m_g4_steps == -1) {
316  ATH_MSG_INFO("Initial step ... preparing event cache.");
317  m_treeData->m_t_x = npos.x();
318  m_treeData->m_t_y = npos.y();
319  m_treeData->m_t_z = npos.z();
320  m_treeData->m_t_theta = nmom.theta();
321  m_treeData->m_t_eta = nmom.eta();
322  m_treeData->m_t_phi = nmom.phi();
323  m_treeData->m_t_p = nmom.mag();
324  m_treeData->m_t_charge = charge;
325  m_treeData->m_t_pdg = pdg;
326  m_treeData->m_g4_steps = 0;
327  m_tX0Cache = 0.;
328 
329  // construct the intial parameters
330  m_parameterCache = new Trk::CurvilinearParameters(npos, nmom, charge);
331  AmgSymMatrix(5) covMatrix;
332  covMatrix.setZero();
333  // covMatrix(0, 0) = 1e-34;
334  // covMatrix(1, 1) = 1e-34;
335  // covMatrix(2, 2) = 1e-34;
336  // covMatrix(3, 4) = 1e-34;
337  // covMatrix(4, 4) = 1e-34;
338  ATH_MSG_DEBUG(" covMatrix " << covMatrix);
339  m_parameterCacheCov = new Trk::CurvilinearParameters(npos, nmom, charge,
340  std::move(covMatrix));
341  ATH_MSG_DEBUG(" Made m_parameterCacheCov with covMatrix "
342  << *m_parameterCacheCov->covariance());
343  return;
344  }
345 
346  float tX0 = X0 > 10e-5 ? t / X0 : 0.;
347  m_tX0Cache += tX0;
348  ATH_MSG_DEBUG(" position R " << npos.perp() << " z " << npos.z() << " X0 "
349  << X0 << " t " << t << " m_tX0Cache "
350  << m_tX0Cache);
351 
352  bool useEntry = true;
353 
354  // Muon Entry or ID exit
355  if (useEntry && !m_crossedEntry &&
356  (std::fabs(npos.z()) > zEntry || npos.perp() > rEntry)) {
357  m_treeData->m_m_x = npos.x();
358  m_treeData->m_m_y = npos.y();
359  m_treeData->m_m_z = npos.z();
360  m_treeData->m_m_theta = nmom.theta();
361  m_treeData->m_m_eta = nmom.eta();
362  m_treeData->m_m_phi = nmom.phi();
363  m_treeData->m_m_p = nmom.mag();
364  // overwrite everything before ME layer
365  m_treeData->m_g4_stepsEntry = 0;
366  // construct the intial parameters
367  m_parameterCacheEntry = new Trk::CurvilinearParameters(npos, nmom, charge);
368  //m_parameterCache = new Trk::CurvilinearParameters(npos, nmom, charge);
369  AmgSymMatrix(5) covMatrix;
370  covMatrix.setZero();
371  m_parameterCacheEntryCov = new Trk::CurvilinearParameters(
372  npos, nmom, charge, std::move(covMatrix));
373  ATH_MSG_DEBUG("m_crossedEntry x "
374  << m_parameterCacheEntry->position().x() << " y "
375  << m_parameterCacheEntry->position().y() << " z "
376  << m_parameterCacheEntry->position().z());
377  m_crossedEntry = true;
378  Trk::CurvilinearParameters g4Parameters =
379  Trk::CurvilinearParameters(npos, nmom, m_treeData->m_t_charge);
380  // Muon Entry
381  m_destinationSurface = g4Parameters.associatedSurface();
382  }
383 
384  // jumping over inital step
385  m_treeData->m_g4_steps =
386  (m_treeData->m_g4_steps == -1) ? 0 : m_treeData->m_g4_steps;
387 
388  if (!m_parameterCache) {
389  ATH_MSG_WARNING("No Parameters available. Bailing out.");
390  return;
391  }
392 
393  if (m_treeData->m_g4_steps >= MAXPROBES) {
394  ATH_MSG_WARNING("Maximum number of " << MAXPROBES
395  << " reached, step is ignored.");
396  return;
397  }
398 
399  // DO NOT store before Entry is crossed (gain CPU)
400  if (!m_crossedEntry) return;
401  if (m_exitLayer) return;
402 
403  // PK 2023
404  // store G4 steps if m_crossedEntry
405  m_treeData->m_g4_p[m_treeData->m_g4_steps] = nmom.mag();
406  m_treeData->m_g4_eta[m_treeData->m_g4_steps] = nmom.eta();
407  m_treeData->m_g4_theta[m_treeData->m_g4_steps] = nmom.theta();
408  m_treeData->m_g4_phi[m_treeData->m_g4_steps] = nmom.phi();
409  m_treeData->m_g4_x[m_treeData->m_g4_steps] = npos.x();
410  m_treeData->m_g4_y[m_treeData->m_g4_steps] = npos.y();
411  m_treeData->m_g4_z[m_treeData->m_g4_steps] = npos.z();
412  m_treeData->m_g4_tX0[m_treeData->m_g4_steps] = m_tX0Cache;
413  m_treeData->m_g4_t[m_treeData->m_g4_steps] = t;
414  m_treeData->m_g4_X0[m_treeData->m_g4_steps] = X0;
415 
416  m_treeData->m_trk_p[m_treeData->m_g4_steps] = 0.;
417  m_treeData->m_trk_eta[m_treeData->m_g4_steps] = 0.;
418  m_treeData->m_trk_theta[m_treeData->m_g4_steps] = 0.;
419  m_treeData->m_trk_phi[m_treeData->m_g4_steps] = 0.;
420  m_treeData->m_trk_x[m_treeData->m_g4_steps] = 0.;
421  m_treeData->m_trk_y[m_treeData->m_g4_steps] = 0.;
422  m_treeData->m_trk_z[m_treeData->m_g4_steps] = 0.;
423  m_treeData->m_trk_lx[m_treeData->m_g4_steps] = 0.;
424  m_treeData->m_trk_ly[m_treeData->m_g4_steps] = 0.;
425  m_treeData->m_trk_eloss[m_treeData->m_g4_steps] = 0.;
426  m_treeData->m_trk_eloss0[m_treeData->m_g4_steps] = 0.;
427  m_treeData->m_trk_eloss1[m_treeData->m_g4_steps] = 0.;
428  m_treeData->m_trk_eloss5[m_treeData->m_g4_steps] = 0.;
429  m_treeData->m_trk_eloss10[m_treeData->m_g4_steps] = 0.;
430  m_treeData->m_trk_scaleeloss[m_treeData->m_g4_steps] = 0.;
431  m_treeData->m_trk_scalex0[m_treeData->m_g4_steps] = 0.;
432  m_treeData->m_trk_x0[m_treeData->m_g4_steps] = 0.;
433  m_treeData->m_trk_status[m_treeData->m_g4_steps] = 0;
434  m_treeData->m_trk_erd0[m_treeData->m_g4_steps] = 0.;
435  m_treeData->m_trk_erz0[m_treeData->m_g4_steps] = 0.;
436  m_treeData->m_trk_erphi[m_treeData->m_g4_steps] = 0.;
437  m_treeData->m_trk_ertheta[m_treeData->m_g4_steps] = 0.;
438  m_treeData->m_trk_erqoverp[m_treeData->m_g4_steps] = 0.;
439  ++m_treeData->m_g4_steps;
440 
441  ATH_MSG_DEBUG("initialise m_treeData->m_g4_steps" << m_treeData->m_g4_steps);
442 
443  bool crossedExitLayer = false;
444  // ID envelope
445  if (std::fabs(npos.z()) > zExit || npos.perp() > rExit)
446  crossedExitLayer = true;
447 
448  ATH_MSG_DEBUG("npos Z: " << npos.z() << "npos prep: " << npos.perp()
449  << "crossedExitLayer: " << crossedExitLayer);
450 
451  if (m_speedup) {
452  ATH_MSG_DEBUG("Starting speed up:"
453  << "m_crossedEntry: " << m_crossedEntry
454  << "m_treeData->m_g4_steps: " << m_treeData->m_g4_steps
455  << "m_treeData->m_g4_stepsEntry: " << m_treeData->m_g4_stepsEntry);
456  if (m_crossedEntry && m_treeData->m_g4_steps >= 2 && !crossedExitLayer)
457  return;
458  }
459 
460  Trk::EnergyLoss eloss = EnergyLoss(0., 0., 0., 0., 0., 0);
461  Trk::ExtrapolationCache extrapolationCache = ExtrapolationCache(0., &eloss);
462 
463  // Cache ONLY used for extrapolateM and extrapolate with covariance Matrix
464 
465  // parameters of the G4 step point
466  Trk::CurvilinearParameters g4Parameters =
467  Trk::CurvilinearParameters(npos, nmom, m_treeData->m_t_charge);
468  // destination surface
469  const Trk::PlaneSurface& destinationSurface =
470  g4Parameters.associatedSurface();
471  // extrapolate to the destination surface
472  std::unique_ptr<Trk::TrackParameters> trkParameters =
473  m_extrapolateDirectly && m_crossedEntry
474  ? m_extrapolator->extrapolateDirectly(
475  ctx, *m_parameterCache, destinationSurface, Trk::alongMomentum,
476  false, Trk::muon)
477  : m_extrapolator->extrapolate(ctx, *m_parameterCache,
478  destinationSurface, Trk::alongMomentum,
479  false, Trk::muon);
480  if (!trkParameters) {
481  ATH_MSG_DEBUG(
482  " G4 extrapolate failed without covariance to destination surface ");
483  }
484  if (m_treeData->m_g4_stepsEntry == 0 && m_useCovMatrix) {
485  ATH_MSG_DEBUG(" Extrapolate m_parameterCacheCov with covMatrix ");
486  extrapolationCache.reset();
487  trkParameters = m_extrapolateDirectly && m_crossedEntry
488  ? m_extrapolator->extrapolateDirectly(
489  ctx, *m_parameterCacheCov, destinationSurface,
490  Trk::alongMomentum, false, Trk::muon)
491  : m_extrapolator->extrapolate(
492  ctx, *m_parameterCacheCov, destinationSurface,
493  Trk::alongMomentum, false, Trk::muon,
494  Trk::addNoise, &extrapolationCache);
495  if (!trkParameters) {
496  ATH_MSG_DEBUG(" G4 extrapolate failed with covariance to Muon Entry or ID Exit");
497  ATH_MSG_DEBUG(" Redo G4 extrapolateM without covariance matrix to Muon Entry or ID Exit ");
498  extrapolationCache.reset();
499  trkParameters = m_extrapolateDirectly && m_crossedEntry
500  ? m_extrapolator->extrapolateDirectly(
501  ctx, *m_parameterCache, destinationSurface,
502  Trk::alongMomentum, false, Trk::muon)
503  : m_extrapolator->extrapolate(
504  ctx, *m_parameterCache, destinationSurface,
505  Trk::alongMomentum, false, Trk::muon,
506  Trk::addNoise, &extrapolationCache);
507 
508  } else {
509  ATH_MSG_DEBUG(
510  " G4 extrapolate succesfull with covariance to Muon Entry or ID exit "
511  << " X0 " << extrapolationCache.x0tot() << " Eloss deltaE "
512  << extrapolationCache.eloss()->deltaE() << " Eloss sigma "
513  << extrapolationCache.eloss()->sigmaDeltaE() << " meanIoni "
514  << extrapolationCache.eloss()->meanIoni() << " sigmaIoni "
515  << extrapolationCache.eloss()->sigmaIoni() << " meanRad "
516  << extrapolationCache.eloss()->meanRad() << " sigmaRad "
517  << extrapolationCache.eloss()->sigmaRad() << " depth "
518  << extrapolationCache.eloss()->length());
519  }
520  }
521 
522  // sroe: coverity 31530
523  m_treeData->m_trk_status[m_treeData->m_g4_steps] = trkParameters ? 1 : 0;
524  ATH_MSG_DEBUG("m_treeData->m_g4_steps: "
525  << m_treeData->m_g4_steps << " exit Layer: " << crossedExitLayer
526  << " track parameters: " << trkParameters.get());
527  ATH_MSG_DEBUG("m_parameterCache: " << m_parameterCache);
528  if (!trkParameters) {
529  return;
530  }
531 
532  ATH_MSG_DEBUG(" exit Layer: " << crossedExitLayer
533  << "track parameters: " << trkParameters.get());
534  if (crossedExitLayer) {
535  ATH_MSG_DEBUG(" exit layer found ");
536  // PK 2023
537  m_exitLayer = true;
538  m_treeData->m_trk_status[m_treeData->m_g4_steps] = 1000;
539  // Get extrapolatio with errors
540  if (m_useCovMatrix) {
541  extrapolationCache.reset();
542  ATH_MSG_DEBUG(" Extrapolate m_parameterCacheEntryCov with covMatrix "
543  << " x " << m_parameterCacheEntryCov->position().x() << " y "
544  << m_parameterCacheEntryCov->position().y() << " z "
545  << m_parameterCacheEntryCov->position().z());
546  ATH_MSG_DEBUG(" m_parameterCacheEntryCov "
547  << "m_extrapolateDirectly: " << m_extrapolateDirectly
548  << "m_crossedEntry: " << m_crossedEntry);
549 
550  trkParameters = m_extrapolateDirectly && m_crossedEntry
551  ? m_extrapolator->extrapolateDirectly(
552  ctx, *m_parameterCacheEntryCov, destinationSurface,
553  Trk::alongMomentum, false, Trk::muon)
554  : m_extrapolator->extrapolate(
555  ctx, *m_parameterCacheEntryCov, destinationSurface,
556  Trk::alongMomentum, false, Trk::muon,
557  Trk::addNoise, &extrapolationCache);
558  if (trkParameters) {
559  ATH_MSG_DEBUG("extrapolation with m_parameterCacheEntryCov succeeded ");
560  } else {
561  ATH_MSG_DEBUG(" extrapolation failed with m_parameterCacheEntryCov ");
562  if (!m_parameterCacheEntryCov) {
563  ATH_MSG_DEBUG(" failed due to m_parameterCacheEntryCov is zero");
564  }
565  extrapolationCache.reset();
566  trkParameters = m_extrapolateDirectly
567  ? m_extrapolator->extrapolateDirectly(
568  ctx, *m_parameterCacheEntry, destinationSurface,
569  Trk::alongMomentum, false, Trk::muon)
570  : m_extrapolator->extrapolate(
571  ctx, *m_parameterCacheEntry, destinationSurface,
572  Trk::alongMomentum, false, Trk::muon,
573  Trk::addNoise, &extrapolationCache);
574  }
575  if (trkParameters)
576  ATH_MSG_DEBUG("extrapolation with m_parameterCacheEntry succeeded");
577  } else {
578  // no covariance matrix
579  extrapolationCache.reset();
580  trkParameters = m_extrapolateDirectly
581  ? m_extrapolator->extrapolateDirectly(
582  ctx, *m_parameterCacheEntry, destinationSurface,
583  Trk::alongMomentum, false, Trk::muon)
584  : m_extrapolator->extrapolate(
585  ctx, *m_parameterCacheEntry, destinationSurface,
586  Trk::alongMomentum, false, Trk::muon,
587  Trk::addNoise, &extrapolationCache);
588  }
589 
590  // Backwards from Exit to ME
591  if (trkParameters) {
592  ATH_MSG_DEBUG(" forward extrapolation succeeded ");
593  bool doBackWard = false;
594  if (doBackWard) {
595  std::unique_ptr<Trk::TrackParameters> trkParameters_BACK =
596  m_extrapolateDirectly
597  ? m_extrapolator->extrapolateDirectly(
598  ctx, *trkParameters, m_destinationSurface,
599  Trk::oppositeMomentum, false, Trk::muon)
600  : m_extrapolator->extrapolate(
601  ctx, *trkParameters, m_destinationSurface,
602  Trk::oppositeMomentum, false, Trk::muon);
603  if (trkParameters_BACK) {
604  ATH_MSG_DEBUG(" back extrapolation succeeded ");
605  m_exitLayer = true;
606  m_treeData->m_b_p = trkParameters_BACK->momentum().mag();
607  m_treeData->m_b_eta = trkParameters_BACK->momentum().eta();
608  m_treeData->m_b_theta = trkParameters_BACK->momentum().theta();
609  m_treeData->m_b_phi = trkParameters_BACK->momentum().phi();
610  m_treeData->m_b_x = trkParameters_BACK->position().x();
611  m_treeData->m_b_y = trkParameters_BACK->position().y();
612  m_treeData->m_b_z = trkParameters_BACK->position().z();
613  if (std::fabs(m_treeData->m_m_p - m_treeData->m_b_p) > 10.)
614  ATH_MSG_DEBUG(
615  " Back extrapolation to Muon Entry or ID Exit finds different "
616  "momentum difference MeV "
617  << m_treeData->m_m_p - m_treeData->m_b_p);
618  // delete trkParameters_BACK;
619  extrapolationCache.reset();
620  const std::vector<const Trk::TrackStateOnSurface*>* matvec_BACK =
621  m_extrapolator->extrapolateM(
622  ctx, *trkParameters, m_destinationSurface,
623  Trk::oppositeMomentum, false, Trk::muon, &extrapolationCache);
624  double Eloss = 0.;
625  double x0 = 0.;
626 
627  int mmat = 0;
628  if (matvec_BACK && !matvec_BACK->empty()) {
629  std::vector<const Trk::TrackStateOnSurface*>::const_iterator it =
630  matvec_BACK->begin();
631  std::vector<const Trk::TrackStateOnSurface*>::const_iterator
632  it_end = matvec_BACK->end();
633  for (; it != it_end; ++it) {
634  const Trk::MaterialEffectsBase* matEf =
635  (*it)->materialEffectsOnTrack();
636  if (matEf) {
637  mmat++;
638  if (m_treeData->m_trk_status[m_treeData->m_g4_steps] == 1000)
639  ATH_MSG_DEBUG(" mmat " << mmat << " matEf->thicknessInX0() "
640  << matEf->thicknessInX0());
641  x0 += matEf->thicknessInX0();
642  const Trk::MaterialEffectsOnTrack* matEfs =
643  dynamic_cast<const Trk::MaterialEffectsOnTrack*>(matEf);
644  if (not matEfs) continue;
645  //
646  double eloss0 = 0.;
647  double meanIoni = 0.;
648  double sigmaIoni = 0.;
649  double meanRad = 0.;
650  double sigmaRad = 0.;
651  double sigmaTheta = 0.;
652  double sigmaPhi = 0.;
653 
654  const Trk::EnergyLoss* eLoss = (matEfs)->energyLoss();
655  if (eLoss) {
656  Eloss += eLoss->deltaE();
657  eloss0 = eLoss->deltaE();
658  meanIoni = eLoss->meanIoni();
659  sigmaIoni = eLoss->sigmaIoni();
660  meanRad = eLoss->meanRad();
661  sigmaRad = eLoss->sigmaRad();
662  if (m_treeData->m_trk_status[m_treeData->m_g4_steps] == 1000)
663  ATH_MSG_DEBUG(" mmat " << mmat << " eLoss->deltaE() "
664  << eLoss->deltaE()
665  << " eLoss->length() "
666  << eLoss->length());
667  }
668 
669  const Trk::ScatteringAngles* scatAng =
670  (matEfs)->scatteringAngles();
671  if (scatAng) {
672  sigmaTheta = scatAng->sigmaDeltaTheta();
673  sigmaPhi = scatAng->sigmaDeltaPhi();
674  }
675  if (m_treeData->m_trk_scats < 500) {
676  // backwards
677  if (m_treeData->m_trk_status[m_treeData->m_g4_steps] == 1000)
678  m_treeData->m_trk_sstatus[m_treeData->m_trk_scats] = -1000;
679  if ((*it)->trackParameters()) {
680  m_treeData->m_trk_sx[m_treeData->m_trk_scats] =
681  (*it)->trackParameters()->position().x();
682  m_treeData->m_trk_sy[m_treeData->m_trk_scats] =
683  (*it)->trackParameters()->position().y();
684  m_treeData->m_trk_sz[m_treeData->m_trk_scats] =
685  (*it)->trackParameters()->position().z();
686  }
687  m_treeData->m_trk_sx0[m_treeData->m_trk_scats] =
688  matEf->thicknessInX0();
689  m_treeData->m_trk_seloss[m_treeData->m_trk_scats] = eloss0;
690  m_treeData->m_trk_smeanIoni[m_treeData->m_trk_scats] =
691  meanIoni;
692  m_treeData->m_trk_ssigIoni[m_treeData->m_trk_scats] =
693  sigmaIoni;
694  m_treeData->m_trk_smeanRad[m_treeData->m_trk_scats] = meanRad;
695  m_treeData->m_trk_ssigRad[m_treeData->m_trk_scats] = sigmaRad;
696  m_treeData->m_trk_ssigTheta[m_treeData->m_trk_scats] =
697  sigmaTheta;
698  m_treeData->m_trk_ssigPhi[m_treeData->m_trk_scats] = sigmaPhi;
699  m_treeData->m_trk_scats++;
700  }
701  }
702  }
703  }
704  m_treeData->m_b_X0 = x0;
705  m_treeData->m_b_Eloss = Eloss;
706  delete matvec_BACK;
707  }
708  }
709  }
710  }
711 
712  extrapolationCache.reset();
713  const std::vector<const Trk::TrackStateOnSurface*>* matvec =
714  m_extrapolator->extrapolateM(ctx, *m_parameterCache, destinationSurface,
715  Trk::alongMomentum, false, Trk::muon,
716  &extrapolationCache);
717 
718  if (matvec)
719  ATH_MSG_DEBUG("MatVec 1: " << matvec->size());
720  else
721  ATH_MSG_DEBUG("MatVec 1: NULL");
722 
723  if (m_useCovMatrix) {
724  ATH_MSG_DEBUG(
725  "m_treeData->m_g4_stepsEntry (debug): " << m_treeData->m_g4_stepsEntry);
726  if (m_treeData->m_g4_stepsEntry <= 0) {
727  extrapolationCache.reset();
728  matvec = m_extrapolator->extrapolateM(
729  ctx, *m_parameterCacheCov, destinationSurface, Trk::alongMomentum,
730  false, Trk::muon, &extrapolationCache);
731  if (!matvec || matvec->empty()) {
732  ATH_MSG_DEBUG(
733  " G4 extrapolateM failed with covariance matrix to Muon Entry or ID Exit ");
734  ATH_MSG_DEBUG(
735  " Redo G4 extrapolateM without covariance matrix to Muon Entry or ID Exit ");
736  extrapolationCache.reset();
737  matvec = m_extrapolator->extrapolateM(
738  ctx, *m_parameterCache, destinationSurface, Trk::alongMomentum,
739  false, Trk::muon, &extrapolationCache);
740  } else {
741  ATH_MSG_DEBUG(
742  " G4 extrapolateM succesfull with covariance matrix to Muon Entry or ID Exit ");
743  }
744  ATH_MSG_DEBUG("From Entry Cache X0 "
745  << extrapolationCache.x0tot() << " Eloss deltaE "
746  << extrapolationCache.eloss()->deltaE() << " Eloss sigma "
747  << extrapolationCache.eloss()->sigmaDeltaE() << " meanIoni "
748  << extrapolationCache.eloss()->meanIoni() << " sigmaIoni "
749  << extrapolationCache.eloss()->sigmaIoni() << " meanRad "
750  << extrapolationCache.eloss()->meanRad() << " sigmaRad "
751  << extrapolationCache.eloss()->sigmaRad());
752  }
753  if (m_treeData->m_g4_stepsEntry == 1) {
754  extrapolationCache.reset();
755  matvec = m_extrapolator->extrapolateM(
756  ctx, *m_parameterCacheEntryCov, destinationSurface, Trk::alongMomentum,
757  false, Trk::muon, &extrapolationCache);
758  if (!matvec || matvec->empty()) {
759  ATH_MSG_DEBUG(
760  " G4 extrapolateM failed with covariance matrix to Muon or Calo Exit ");
761  ATH_MSG_DEBUG(
762  " Redo G4 extrapolateM without covariance matrix to Muon or Calo Exit ");
763  extrapolationCache.reset();
764  matvec = m_extrapolator->extrapolateM(
765  ctx, *m_parameterCacheEntry, destinationSurface, Trk::alongMomentum,
766  false, Trk::muon, &extrapolationCache);
767  } else {
768  ATH_MSG_DEBUG(
769  " G4 extrapolateM succesfull with covariance matrix to Muon or Calo Exit ");
770  }
771  ATH_MSG_DEBUG("From Muon or Calo Exit Cache X0 "
772  << extrapolationCache.x0tot() << " Eloss deltaE "
773  << extrapolationCache.eloss()->deltaE() << " Eloss sigma "
774  << extrapolationCache.eloss()->sigmaDeltaE() << " meanIoni "
775  << extrapolationCache.eloss()->meanIoni() << " sigmaIoni "
776  << extrapolationCache.eloss()->sigmaIoni() << " meanRad "
777  << extrapolationCache.eloss()->meanRad() << " sigmaRad "
778  << extrapolationCache.eloss()->sigmaRad());
779  }
780  } else {
781  if (m_treeData->m_g4_stepsEntry == 1) {
782  extrapolationCache.reset();
783  matvec = m_extrapolator->extrapolateM(
784  ctx, *m_parameterCacheEntry, destinationSurface, Trk::alongMomentum,
785  false, Trk::muon, &extrapolationCache);
786  ATH_MSG_DEBUG(" G4 extrapolateM without covariance matrix to Muon Entry or ID Exit "
787  << " X0 " << extrapolationCache.x0tot() << " Eloss deltaE "
788  << extrapolationCache.eloss()->deltaE() << " Eloss sigma "
789  << extrapolationCache.eloss()->sigmaDeltaE() << " meanIoni "
790  << extrapolationCache.eloss()->meanIoni() << " sigmaIoni "
791  << extrapolationCache.eloss()->sigmaIoni() << " meanRad "
792  << extrapolationCache.eloss()->meanRad() << " sigmaRad "
793  << extrapolationCache.eloss()->sigmaRad());
794  }
795  }
796 
797  double Elosst = 0.;
798  const std::vector<const Trk::TrackStateOnSurface*> matvecNewRepAggrUp =
799  modifyTSOSvector(*matvec, 1.0, 1.0, true, true, true, 0., 0., 10000., 0.,
800  Elosst);
801 
802  double X0Scale = 1.0;
803  double ElossScale = 1.0;
804 
805  double Eloss0 = 0.;
806  double Eloss1 = 0.;
807  double Eloss5 = 0.;
808  double Eloss10 = 0.;
809  bool system1 = false; // ID or Calorimeter
810  bool system2 = false; // Calorimeter or Muon system
811 
812  if (!matvec->empty()) {
813  if (m_crossedEntry && !m_exitLayer) system1 = true;
814  if (m_crossedEntry && m_exitLayer) system2 = true;
815  }
816 
817  ATH_MSG_DEBUG(" ID or Calorimeter system " << system1 << " Calorimeter or Muon system " << system2);
818  if (system2) {
819  //
820  // Calorimeter or Muon system
821  //
822  m_elossupdator->getX0ElossScales(0, m_treeData->m_m_eta,
823  m_treeData->m_m_phi, X0Scale, ElossScale);
824  ATH_MSG_DEBUG(" system2 scales X0 " << X0Scale << " ElossScale "
825  << ElossScale);
826 
827  const std::vector<const Trk::TrackStateOnSurface*> matvecNew1 =
828  modifyTSOSvector(*matvec, X0Scale, 1., true, true, true, 0., 0.,
829  m_treeData->m_m_p, 0., Eloss1);
830  const std::vector<const Trk::TrackStateOnSurface*> matvecNew0 =
831  modifyTSOSvector(*matvec, X0Scale, ElossScale, true, true, true, 0., 0.,
832  m_treeData->m_m_p, 0., Eloss0);
833  ATH_MSG_DEBUG(" muon system modify with 5 percent ");
834  const std::vector<const Trk::TrackStateOnSurface*> matvecNew5 =
835  modifyTSOSvector(*matvec, X0Scale, ElossScale, true, true, true, 0., 0.,
836  m_treeData->m_m_p, 0.05 * m_treeData->m_m_p, Eloss5);
837  ATH_MSG_DEBUG(" muon system modify with 10 percent ");
838  const std::vector<const Trk::TrackStateOnSurface*> matvecNew10 =
839  modifyTSOSvector(*matvec, X0Scale, ElossScale, true, true, true, 0., 0.,
840  m_treeData->m_m_p, 0.10 * m_treeData->m_m_p, Eloss10);
841  }
842  if (system1) {
843  //
844  // ID or Calorimeter system
845  //
846  double phiCaloExit = atan2(m_treeData->m_m_y, m_treeData->m_m_x);
847  m_elossupdator->getX0ElossScales(1, m_treeData->m_t_eta, phiCaloExit,
848  X0Scale, ElossScale);
849  ATH_MSG_DEBUG(" calorimeter scales X0 " << X0Scale << " ElossScale "
850  << ElossScale);
851  const std::vector<const Trk::TrackStateOnSurface*> matvecNew1 =
852  modifyTSOSvector(*matvec, X0Scale, 1., true, true, true, 0., 0.,
853  m_treeData->m_m_p, 0., Eloss1);
854  const std::vector<const Trk::TrackStateOnSurface*> matvecNew0 =
855  modifyTSOSvector(*matvec, X0Scale, ElossScale, true, true, true, 0., 0.,
856  m_treeData->m_t_p, 0., Eloss0);
857  if (std::fabs(Eloss1) > 0)
858  ATH_MSG_DEBUG(" **** Cross Check calorimeter with Eloss Scale1 "
859  << Eloss1 << " Eloss0 " << Eloss0 << " ratio "
860  << Eloss0 / Eloss1);
861 
862  ATH_MSG_DEBUG(" calorimeter modify with 5 percent ");
863  const std::vector<const Trk::TrackStateOnSurface*> matvecNew5 =
864  modifyTSOSvector(*matvec, X0Scale, ElossScale, true, true, true, 0., 0.,
865  m_treeData->m_t_p, 0.05 * m_treeData->m_m_p, Eloss5);
866  ATH_MSG_DEBUG(" calorimeter modify with 10 percent ");
867  const std::vector<const Trk::TrackStateOnSurface*> matvecNew10 =
868  modifyTSOSvector(*matvec, X0Scale, ElossScale, true, true, true, 0., 0.,
869  m_treeData->m_t_p, 0.10 * m_treeData->m_m_p, Eloss10);
870  }
871 
872  ATH_MSG_DEBUG(" status " << m_treeData->m_trk_status[m_treeData->m_g4_steps]
873  << "Eloss1 " << Eloss1 << " Eloss0 " << Eloss0
874  << " Eloss5 " << Eloss5 << " Eloss10 " << Eloss10);
875 
876  double Eloss = 0.;
877  double x0 = 0.;
878 
879  int mmat = 0;
880  // PK 2023 only add scatterers for the ID or Calorimeter
881  if (!(matvec->empty()) && m_treeData->m_g4_stepsEntry <= 1) {
882  std::vector<const Trk::TrackStateOnSurface*>::const_iterator it =
883  matvec->begin();
884  std::vector<const Trk::TrackStateOnSurface*>::const_iterator it_end =
885  matvec->end();
886  for (; it != it_end; ++it) {
887  const Trk::MaterialEffectsBase* matEf = (*it)->materialEffectsOnTrack();
888  if (matEf) {
889  mmat++;
890  if (m_treeData->m_trk_status[m_treeData->m_g4_steps] == 1000)
891  ATH_MSG_DEBUG(" mmat " << mmat << " matEf->thicknessInX0() "
892  << matEf->thicknessInX0());
893  x0 += matEf->thicknessInX0();
894  const Trk::MaterialEffectsOnTrack* matEfs =
895  dynamic_cast<const Trk::MaterialEffectsOnTrack*>(matEf);
896  double eloss0 = 0.;
897  double meanIoni = 0.;
898  double sigmaIoni = 0.;
899  double meanRad = 0.;
900  double sigmaRad = 0.;
901  double sigmaTheta = 0.;
902  double sigmaPhi = 0.;
903  if (matEfs) {
904  const Trk::EnergyLoss* eLoss = (matEfs)->energyLoss();
905  if (eLoss) {
906  Eloss += eLoss->deltaE();
907  eloss0 = eLoss->deltaE();
908  meanIoni = eLoss->meanIoni();
909  sigmaIoni = eLoss->sigmaIoni();
910  meanRad = eLoss->meanRad();
911  sigmaRad = eLoss->sigmaRad();
912  ATH_MSG_DEBUG("m_treeData->m_g4_stepsEntry "
913  << m_treeData->m_g4_stepsEntry << " mmat " << mmat
914  << " X0 " << matEf->thicknessInX0()
915  << " eLoss->deltaE() " << eLoss->deltaE()
916  << " meanIoni " << meanIoni << " Total Eloss "
917  << Eloss << " eLoss->length() " << eLoss->length());
918  }
919  }
920  // sroe: coverity 31532
921  const Trk::ScatteringAngles* scatAng =
922  (matEfs) ? ((matEfs)->scatteringAngles()) : (nullptr);
923 
924  if (scatAng) {
925  sigmaTheta = scatAng->sigmaDeltaTheta();
926  sigmaPhi = scatAng->sigmaDeltaPhi();
927  ATH_MSG_DEBUG("m_treeData->m_g4_stepsEntry "
928  << m_treeData->m_g4_stepsEntry << " mmat " << mmat
929  << " sigmaTheta " << sigmaTheta << " sigmaPhi "
930  << sigmaPhi);
931  }
932 
933  if (m_treeData->m_trk_scats < 500) {
934  if (m_treeData->m_g4_stepsEntry == 0 ||
935  m_treeData->m_trk_status[m_treeData->m_g4_steps] == 1000) {
936  // forwards
937  if (m_treeData->m_g4_stepsEntry == 0)
938  m_treeData->m_trk_sstatus[m_treeData->m_trk_scats] = 10;
939  if (m_treeData->m_trk_status[m_treeData->m_g4_steps] == 1000)
940  m_treeData->m_trk_sstatus[m_treeData->m_trk_scats] = 1000;
941  if ((*it)->trackParameters()) {
942  m_treeData->m_trk_sx[m_treeData->m_trk_scats] =
943  (*it)->trackParameters()->position().x();
944  m_treeData->m_trk_sy[m_treeData->m_trk_scats] =
945  (*it)->trackParameters()->position().y();
946  m_treeData->m_trk_sz[m_treeData->m_trk_scats] =
947  (*it)->trackParameters()->position().z();
948  }
949  m_treeData->m_trk_sx0[m_treeData->m_trk_scats] =
950  matEf->thicknessInX0();
951  m_treeData->m_trk_seloss[m_treeData->m_trk_scats] = eloss0;
952  m_treeData->m_trk_smeanIoni[m_treeData->m_trk_scats] = meanIoni;
953  m_treeData->m_trk_ssigIoni[m_treeData->m_trk_scats] = sigmaIoni;
954  m_treeData->m_trk_smeanRad[m_treeData->m_trk_scats] = meanRad;
955  m_treeData->m_trk_ssigRad[m_treeData->m_trk_scats] = sigmaRad;
956  m_treeData->m_trk_ssigTheta[m_treeData->m_trk_scats] = sigmaTheta;
957  m_treeData->m_trk_ssigPhi[m_treeData->m_trk_scats] = sigmaPhi;
958  m_treeData->m_trk_scats++;
959  }
960  }
961  }
962  }
963  delete matvec;
964  }
965 
966  ATH_MSG_DEBUG(" m_treeData->m_g4_steps "
967  << m_treeData->m_g4_steps << " Radius " << npos.perp() << " z "
968  << npos.z() << " size matvec "
969  << " total X0 " << x0 << " total Eloss " << Eloss);
970 
971  // go back and refill information
972  // PK 2023
973  if (m_treeData->m_g4_steps > 0) --m_treeData->m_g4_steps;
974  // fill the geant information and the trk information
975  m_treeData->m_g4_p[m_treeData->m_g4_steps] = nmom.mag();
976  m_treeData->m_g4_eta[m_treeData->m_g4_steps] = nmom.eta();
977  m_treeData->m_g4_theta[m_treeData->m_g4_steps] = nmom.theta();
978  m_treeData->m_g4_phi[m_treeData->m_g4_steps] = nmom.phi();
979  m_treeData->m_g4_x[m_treeData->m_g4_steps] = npos.x();
980  m_treeData->m_g4_y[m_treeData->m_g4_steps] = npos.y();
981  m_treeData->m_g4_z[m_treeData->m_g4_steps] = npos.z();
982  m_treeData->m_g4_tX0[m_treeData->m_g4_steps] = m_tX0Cache;
983  m_treeData->m_g4_t[m_treeData->m_g4_steps] = t;
984  m_treeData->m_g4_X0[m_treeData->m_g4_steps] = X0;
985 
986  m_treeData->m_trk_p[m_treeData->m_g4_steps] =
987  trkParameters ? trkParameters->momentum().mag() : 0.;
988  m_treeData->m_trk_eta[m_treeData->m_g4_steps] =
989  trkParameters ? trkParameters->momentum().eta() : 0.;
990  m_treeData->m_trk_theta[m_treeData->m_g4_steps] =
991  trkParameters ? trkParameters->momentum().theta() : 0.;
992  m_treeData->m_trk_phi[m_treeData->m_g4_steps] =
993  trkParameters ? trkParameters->momentum().phi() : 0.;
994  m_treeData->m_trk_x[m_treeData->m_g4_steps] =
995  trkParameters ? trkParameters->position().x() : 0.;
996  m_treeData->m_trk_y[m_treeData->m_g4_steps] =
997  trkParameters ? trkParameters->position().y() : 0.;
998  m_treeData->m_trk_z[m_treeData->m_g4_steps] =
999  trkParameters ? trkParameters->position().z() : 0.;
1000  m_treeData->m_trk_lx[m_treeData->m_g4_steps] =
1001  trkParameters ? trkParameters->parameters()[Trk::locX] : 0.;
1002  m_treeData->m_trk_ly[m_treeData->m_g4_steps] =
1003  trkParameters ? trkParameters->parameters()[Trk::locY] : 0.;
1004  m_treeData->m_trk_eloss[m_treeData->m_g4_steps] = Eloss;
1005  m_treeData->m_trk_eloss0[m_treeData->m_g4_steps] = Eloss0;
1006  m_treeData->m_trk_eloss1[m_treeData->m_g4_steps] = Eloss1;
1007  m_treeData->m_trk_eloss5[m_treeData->m_g4_steps] = Eloss5;
1008  m_treeData->m_trk_eloss10[m_treeData->m_g4_steps] = Eloss10;
1009  m_treeData->m_trk_scaleeloss[m_treeData->m_g4_steps] = ElossScale;
1010  m_treeData->m_trk_scalex0[m_treeData->m_g4_steps] = X0Scale;
1011  m_treeData->m_trk_x0[m_treeData->m_g4_steps] = x0;
1012  if (m_treeData->m_g4_stepsEntry == 0)
1013  m_treeData->m_trk_status[m_treeData->m_g4_steps] = 10;
1014  else
1015  m_treeData->m_trk_status[m_treeData->m_g4_steps] = 1000;
1016 
1017  double errord0 = 0.;
1018  double errorz0 = 0.;
1019  double errorphi = 0.;
1020  double errortheta = 0.;
1021  double errorqoverp = 0.;
1022  if (trkParameters && trkParameters->covariance()) {
1023  errord0 = (*trkParameters->covariance())(Trk::d0, Trk::d0);
1024  errorz0 = (*trkParameters->covariance())(Trk::z0, Trk::z0);
1025  errorphi = (*trkParameters->covariance())(Trk::phi, Trk::phi);
1026  errortheta = (*trkParameters->covariance())(Trk::theta, Trk::theta);
1027  errorqoverp = (*trkParameters->covariance())(Trk::qOverP, Trk::qOverP);
1028  ATH_MSG_DEBUG(" Covariance found for m_treeData->m_trk_status "
1029  << m_treeData->m_trk_status[m_treeData->m_g4_steps]);
1030  }
1031 
1032  m_treeData->m_trk_erd0[m_treeData->m_g4_steps] = sqrt(errord0);
1033  m_treeData->m_trk_erz0[m_treeData->m_g4_steps] = sqrt(errorz0);
1034  m_treeData->m_trk_erphi[m_treeData->m_g4_steps] = sqrt(errorphi);
1035  m_treeData->m_trk_ertheta[m_treeData->m_g4_steps] = sqrt(errortheta);
1036  m_treeData->m_trk_erqoverp[m_treeData->m_g4_steps] = sqrt(errorqoverp);
1037 
1038  // reset X0 at Muon Entry or ID Exit
1039  if (m_treeData->m_g4_stepsEntry == 0) m_tX0Cache = 0.;
1040  // update the parameters if needed/configured
1041  if (m_extrapolateIncrementally && trkParameters) {
1042  delete m_parameterCache;
1043  // Unsure what to do here?
1044  // m_parameterCache = trkParameters;
1045  }
1046 
1047  ++m_treeData->m_g4_steps;
1048  if (m_treeData->m_g4_stepsEntry != -1) ++m_treeData->m_g4_stepsEntry;
1049 }
1050 
1051 std::vector<const Trk::TrackStateOnSurface*>
1052 Trk::GeantFollowerMSHelper::modifyTSOSvector(
1053  const std::vector<const Trk::TrackStateOnSurface*>& matvec, double scaleX0,
1054  double scaleEloss, bool reposition, bool aggregate, bool updateEloss,
1055  double caloEnergy, double caloEnergyError, double pCaloEntry,
1056  double momentumError, double& Eloss_tot) const {
1057  //
1058  // inputs: TSOSs for material (matvec) and scale factors for X0 (scaleX0) and
1059  // Eloss (scaleEloss)
1060  //
1061  // returns: new vector of TSOSs including scaling of X0 and Eloss;
1062  //
1063  // options:
1064  // bool reposition correct repositioning of the scattering centers in space
1065  // bool aggregate put scattering centra together in two planes
1066  // bool update Eloss correct energy loss 1) including the measured
1067  // calorimeter Eloss 2) include smearing of the muon momentum
1068  //
1069  // the routine should NOT be called for the ID
1070  // for best use in the Calorimeter: bool reposition = true, bool
1071  // aggregate = true and updateEloss = true (measured caloEnergy and
1072  // caloEnergyError should be passed)
1073  // note that the updateEloss is only
1074  // active with aggregate = true
1075  // for best use in the Muon Specrometer: bool reposition = true, bool
1076  // aggregate = true and updateEloss = false
1077  //
1078  // if one runs with reposition = false the scattering centra are kept at the
1079  // END of the thick/dense material: that is not right for thick material for
1080  // thin it is OK
1081  //
1082  std::vector<const Trk::TrackStateOnSurface*> newTSOSvector;
1083  int maxsize = 2 * matvec.size();
1084  if (aggregate) maxsize = 2;
1085  newTSOSvector.reserve(maxsize);
1086  //
1087  // initialize total sum variables
1088  //
1089  //
1090  Eloss_tot = 0.;
1091 
1092  double X0_tot = 0.;
1093 
1094  double sigmaDeltaPhi2_tot = 0.;
1095  double sigmaDeltaTheta2_tot = 0.;
1096  double deltaE_tot = 0.;
1097  double sigmaDeltaE_tot = 0.;
1098  double sigmaPlusDeltaE_tot = 0.;
1099  double sigmaMinusDeltaE_tot = 0.;
1100  double deltaE_ioni_tot = 0.;
1101  double sigmaDeltaE_ioni_tot = 0.;
1102  double deltaE_rad_tot = 0.;
1103  double sigmaDeltaE_rad_tot = 0.;
1104 
1105  const Trk::TrackStateOnSurface* mprevious = nullptr;
1106  const Trk::TrackStateOnSurface* mfirst = nullptr;
1107  const Trk::TrackStateOnSurface* mlast = nullptr;
1108  Amg::Vector3D posFirst(0., 0., 0.);
1109 
1110  double deltaEFirst = 0.;
1111 
1112  double deltaPhi = 0.;
1113  double deltaTheta = 0.;
1114 
1115  int n_tot = 0;
1116 
1117  double w_tot = 0.;
1118  double wdist2 = 0.;
1119  Amg::Vector3D wdir(0., 0., 0.);
1120  Amg::Vector3D wpos(0., 0., 0.);
1121 
1122  std::bitset<Trk::MaterialEffectsBase::NumberOfMaterialEffectsTypes>
1123  meotPattern(0);
1124  meotPattern.set(Trk::MaterialEffectsBase::EnergyLossEffects);
1125  meotPattern.set(Trk::MaterialEffectsBase::ScatteringEffects);
1126  // meotPattern.set(Trk::MaterialEffectsBase::FittedMaterialEffects);
1127 
1128  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1129  typePattern(0);
1130  typePattern.set(Trk::TrackStateOnSurface::InertMaterial);
1131  typePattern.set(Trk::TrackStateOnSurface::Scatterer);
1132 
1133  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1134  typePatternDeposit(0);
1135  typePatternDeposit.set(Trk::TrackStateOnSurface::InertMaterial);
1136  typePatternDeposit.set(Trk::TrackStateOnSurface::CaloDeposit);
1137  typePatternDeposit.set(Trk::TrackStateOnSurface::Scatterer);
1138 
1139  for (const auto* m : matvec) {
1140  if (!m->trackParameters()) {
1141  ATH_MSG_WARNING("No trackparameters on TrackStateOnSurface ");
1142  continue;
1143  }
1144  if (m->materialEffectsOnTrack()) {
1145  double X0 = m->materialEffectsOnTrack()->thicknessInX0();
1146  const Trk::MaterialEffectsOnTrack* meot =
1147  dynamic_cast<const Trk::MaterialEffectsOnTrack*>(
1148  m->materialEffectsOnTrack());
1149  const Trk::EnergyLoss* energyLoss = nullptr;
1150  const Trk::ScatteringAngles* scat = nullptr;
1151  if (meot) {
1152  energyLoss = meot->energyLoss();
1153  if (energyLoss) {
1154  // double deltaE = energyLoss->deltaE();
1155  } else {
1156  ATH_MSG_WARNING("No energyLoss on TrackStateOnSurface ");
1157  continue;
1158  }
1159  scat = meot->scatteringAngles();
1160  if (scat) {
1161  // double dtheta = scat->sigmaDeltaTheta();
1162  } else {
1163  ATH_MSG_WARNING("No scatteringAngles on TrackStateOnSurface ");
1164  continue;
1165  }
1166  } else {
1167  ATH_MSG_WARNING("No materialEffectsOnTrack on TrackStateOnSurface ");
1168  continue;
1169  }
1170 
1171  double depth = energyLoss->length();
1172  ATH_MSG_DEBUG(" ");
1173  ATH_MSG_DEBUG(" original TSOS type "
1174  << m->dumpType() << " TSOS surface "
1175  << m->trackParameters()->associatedSurface()
1176  << " position x " << m->trackParameters()->position().x()
1177  << " y " << m->trackParameters()->position().y() << " z "
1178  << m->trackParameters()->position().z() << " direction x "
1179  << m->trackParameters()->momentum().unit().x() << " y "
1180  << m->trackParameters()->momentum().unit().y() << " z "
1181  << m->trackParameters()->momentum().unit().z() << " p "
1182  << m->trackParameters()->momentum().mag() << " X0 " << X0
1183  << " deltaE " << energyLoss->deltaE()
1184  << " sigma deltaTheta " << scat->sigmaDeltaTheta()
1185  << " depth " << depth);
1186 
1187  X0_tot += scaleX0 * X0;
1188 
1189  sigmaDeltaTheta2_tot +=
1190  scaleX0 * scat->sigmaDeltaTheta() * scat->sigmaDeltaTheta();
1191  sigmaDeltaPhi2_tot +=
1192  scaleX0 * scat->sigmaDeltaPhi() * scat->sigmaDeltaPhi();
1193 
1194  // Eloss sigma values add up linearly for Landau and exponential
1195  // distributions
1196 
1197  deltaE_tot += scaleEloss * energyLoss->deltaE();
1198  sigmaDeltaE_tot += scaleEloss * energyLoss->sigmaDeltaE();
1199  sigmaPlusDeltaE_tot += scaleEloss * energyLoss->sigmaPlusDeltaE();
1200  sigmaMinusDeltaE_tot += scaleEloss * energyLoss->sigmaMinusDeltaE();
1201  deltaE_ioni_tot += scaleEloss * energyLoss->meanIoni();
1202  sigmaDeltaE_ioni_tot += scaleEloss * energyLoss->sigmaIoni();
1203  deltaE_rad_tot += scaleEloss * energyLoss->meanRad();
1204  sigmaDeltaE_rad_tot += scaleEloss * energyLoss->sigmaRad();
1205 
1206  n_tot++;
1207 
1208  Amg::Vector3D dir = m->trackParameters()->momentum().unit();
1209  Amg::Vector3D pos = m->trackParameters()->position();
1210  if (mprevious) {
1211  dir += mprevious->trackParameters()->momentum().unit();
1212  }
1213 
1214  dir = dir / dir.mag();
1215  ATH_MSG_DEBUG(" position at end " << pos.x() << " y " << pos.y() << " z "
1216  << pos.z() << " perp " << pos.perp());
1217  ATH_MSG_DEBUG(" direction x " << dir.x() << " y " << dir.y() << " z "
1218  << dir.z());
1219  Amg::Vector3D pos0 = pos - (depth / 2. + depth / sqrt(12.)) * dir;
1220  Amg::Vector3D posNew = pos - (depth / 2. - depth / sqrt(12.)) * dir;
1221 
1222  ATH_MSG_DEBUG(" position scattering centre0 x "
1223  << pos0.x() << " y " << pos0.y() << " z " << pos0.z()
1224  << " perp " << pos0.perp());
1225  ATH_MSG_DEBUG(" position scattering centre1 x "
1226  << posNew.x() << " y " << posNew.y() << " z " << posNew.z()
1227  << " perp " << posNew.perp() << " distance "
1228  << (pos0 - posNew).mag() << " depth " << depth);
1229  if (!mfirst) {
1230  mfirst = m;
1231  posFirst = pos0;
1232  deltaEFirst = energyLoss->deltaE();
1233  }
1234  mlast = m;
1235 
1236  double w = scat->sigmaDeltaTheta() * scat->sigmaDeltaTheta();
1237  w_tot += w;
1238  wpos += w * pos0 / 2.;
1239  wpos += w * posNew / 2.;
1240  wdir += w * dir;
1241 
1242  wdist2 += w * (pos0 - posFirst).mag2() / 2.;
1243  wdist2 += w * (posNew - posFirst).mag2() / 2.;
1244 
1245  if (!aggregate && !reposition) {
1246  auto scatNew = ScatteringAngles(deltaPhi, deltaTheta,
1247  std::sqrt(sigmaDeltaPhi2_tot),
1248  std::sqrt(sigmaDeltaTheta2_tot));
1249  auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1250  deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1251  sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1252  deltaE_rad_tot, sigmaDeltaE_rad_tot, depth);
1253  Eloss_tot += energyLossNew->deltaE();
1254  const Trk::Surface& surf = *(meot->associatedSurface().clone());
1255  auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1256  X0_tot, scatNew, std::move(energyLossNew), surf, meotPattern);
1257  auto pars = m->trackParameters()->uniqueClone();
1258 
1259  // make new TSOS
1260  const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1261  nullptr, std::move(pars), std::move(meotLast), typePattern);
1262  newTSOSvector.push_back(newTSOS);
1263 
1264  X0_tot = 0.;
1265  sigmaDeltaTheta2_tot = 0.;
1266  sigmaDeltaPhi2_tot = 0.;
1267  deltaE_tot = 0.;
1268  sigmaDeltaE_tot = 0;
1269  sigmaPlusDeltaE_tot = 0.;
1270  sigmaMinusDeltaE_tot = 0.;
1271  deltaE_ioni_tot = 0.;
1272  sigmaDeltaE_ioni_tot = 0.;
1273  deltaE_rad_tot = 0.;
1274  sigmaDeltaE_rad_tot = 0.;
1275 
1276  } else if (!aggregate && reposition) {
1277  if (std::abs(depth) < 10.) {
1278  auto scatNew = ScatteringAngles(deltaPhi, deltaTheta,
1279  std::sqrt(sigmaDeltaPhi2_tot),
1280  std::sqrt(sigmaDeltaTheta2_tot));
1281  auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1282  deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1283  sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1284  deltaE_rad_tot, sigmaDeltaE_rad_tot, depth);
1285  const Trk::Surface& surf = *(meot->associatedSurface().clone());
1286  Eloss_tot += energyLossNew->deltaE();
1287  auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1288  X0_tot, scatNew, std::move(energyLossNew), surf, meotPattern);
1289  std::unique_ptr<Trk::TrackParameters> pars =
1290  m->trackParameters()->uniqueClone();
1291  // make new TSOS
1292  const Trk::TrackStateOnSurface* newTSOS =
1293  new Trk::TrackStateOnSurface(nullptr, std::move(pars),
1294  std::move(meotLast), typePattern);
1295  newTSOSvector.push_back(newTSOS);
1296  X0_tot = 0.;
1297  sigmaDeltaTheta2_tot = 0.;
1298  sigmaDeltaPhi2_tot = 0.;
1299  deltaE_tot = 0.;
1300  sigmaDeltaE_tot = 0;
1301  sigmaPlusDeltaE_tot = 0.;
1302  sigmaMinusDeltaE_tot = 0.;
1303  deltaE_ioni_tot = 0.;
1304  sigmaDeltaE_ioni_tot = 0.;
1305  deltaE_rad_tot = 0.;
1306  sigmaDeltaE_rad_tot = 0.;
1307 
1308  } else {
1309  //
1310  // Thick scatterer: make two TSOSs
1311  //
1312  // prepare for first MaterialEffectsOnTrack with X0 = X0/2
1313  // Eloss = 0 and scattering2 = total2 / 2. depth = 0
1314  auto energyLoss0 = std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.);
1315  auto scatFirst = ScatteringAngles(deltaPhi, deltaTheta,
1316  sqrt(sigmaDeltaPhi2_tot / 2.),
1317  sqrt(sigmaDeltaTheta2_tot / 2.));
1318 
1319  // prepare for second MaterialEffectsOnTrack with X0 = X0/2
1320  // Eloss = Eloss total and scattering2 = total2 / 2. depth = 0
1321  auto scatNew = ScatteringAngles(deltaPhi, deltaTheta,
1322  sqrt(sigmaDeltaPhi2_tot / 2.),
1323  sqrt(sigmaDeltaTheta2_tot / 2.));
1324  auto energyLossNew = std::make_unique<Trk::EnergyLoss>(
1325  deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1326  sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1327  deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.);
1328  double norm = dir.perp();
1329  // Rotation matrix representation
1330  Amg::Vector3D colx(-dir.y() / norm, dir.x() / norm, 0);
1331  Amg::Vector3D coly(-dir.x() * dir.z() / norm,
1332  -dir.y() * dir.z() / norm, norm);
1333  Amg::Vector3D colz(dir.x(), dir.y(), dir.z());
1334 
1335  Amg::Transform3D surfaceTransformFirst(colx, coly, colz, pos0);
1336  Amg::Transform3D surfaceTransformLast(colx, coly, colz, posNew);
1337  Trk::PlaneSurface* surfFirst =
1338  new Trk::PlaneSurface(surfaceTransformFirst);
1339  Trk::PlaneSurface* surfLast =
1340  new Trk::PlaneSurface(surfaceTransformLast);
1341  Eloss_tot += energyLossNew->deltaE();
1342  // make MaterialEffectsOnTracks
1343  auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1344  X0_tot / 2., scatFirst, std::move(energyLoss0), *surfFirst,
1345  meotPattern);
1346  auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1347  X0_tot / 2., scatNew, std::move(energyLossNew), *surfLast,
1348  meotPattern);
1349 
1350  // calculate TrackParameters at first surface
1351  double qOverP0 = m->trackParameters()->charge() /
1352  (m->trackParameters()->momentum().mag() -
1353  std::fabs(energyLoss->deltaE()));
1354  if (mprevious)
1355  qOverP0 = mprevious->trackParameters()->charge() /
1356  mprevious->trackParameters()->momentum().mag();
1357  std::unique_ptr<Trk::TrackParameters> parsFirst =
1358  surfFirst->createUniqueParameters<5, Trk::Charged>(
1359  0., 0., dir.phi(), dir.theta(), qOverP0);
1360  // calculate TrackParameters at second surface
1361  double qOverPNew = m->trackParameters()->charge() /
1362  m->trackParameters()->momentum().mag();
1363  std::unique_ptr<Trk::TrackParameters> parsLast =
1364  surfLast->createUniqueParameters<5, Trk::Charged>(
1365  0., 0., dir.phi(), dir.theta(), qOverPNew);
1366  // make TSOS
1367  //
1368  const Trk::TrackStateOnSurface* newTSOSFirst =
1369  new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1370  std::move(meotFirst), typePattern);
1371  const Trk::TrackStateOnSurface* newTSOS =
1372  new Trk::TrackStateOnSurface(nullptr, std::move(parsLast),
1373  std::move(meotLast), typePattern);
1374 
1375  newTSOSvector.push_back(newTSOSFirst);
1376  newTSOSvector.push_back(newTSOS);
1377 
1378  X0_tot = 0.;
1379  sigmaDeltaTheta2_tot = 0.;
1380  sigmaDeltaPhi2_tot = 0.;
1381  deltaE_tot = 0.;
1382  sigmaDeltaE_tot = 0;
1383  sigmaPlusDeltaE_tot = 0.;
1384  sigmaMinusDeltaE_tot = 0.;
1385  deltaE_ioni_tot = 0.;
1386  sigmaDeltaE_ioni_tot = 0.;
1387  deltaE_rad_tot = 0.;
1388  sigmaDeltaE_rad_tot = 0.;
1389  }
1390  }
1391 
1392  mprevious = m;
1393  }
1394  }
1395  if (aggregate && reposition) {
1396  if (n_tot > 0) {
1397  //
1398  // Make three scattering planes in Calorimeter else make two
1399  //
1400  Amg::Vector3D pos = wpos / w_tot;
1401  bool threePlanes = false;
1402  if (X0_tot > 50 && std::fabs(pos.z()) < 6700 && pos.perp() < 4200)
1403  threePlanes = true;
1404  //
1405  auto energyLoss0 = std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.);
1406  auto scatFirst =
1407  ScatteringAngles(deltaPhi, deltaTheta, sqrt(sigmaDeltaPhi2_tot / 2.),
1408  sqrt(sigmaDeltaTheta2_tot / 2.));
1409 
1410  auto scatNew =
1411  ScatteringAngles(deltaPhi, deltaTheta, sqrt(sigmaDeltaPhi2_tot / 2.),
1412  sqrt(sigmaDeltaTheta2_tot / 2.));
1413  auto energyLoss2 = Trk::EnergyLoss(
1414  deltaE_tot, sigmaDeltaE_tot, sigmaPlusDeltaE_tot,
1415  sigmaMinusDeltaE_tot, deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1416  deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.);
1417 
1418  int elossFlag = 0; // return Flag for updateEnergyLoss Calorimeter
1419  // energy (0 = not used)
1420  auto energyLossNew =
1421  (updateEloss
1422  ? m_elossupdator->updateEnergyLoss(energyLoss2, caloEnergy,
1423  caloEnergyError, pCaloEntry,
1424  momentumError, elossFlag)
1425  : Trk::EnergyLoss(deltaE_tot, sigmaDeltaE_tot,
1426  sigmaPlusDeltaE_tot, sigmaMinusDeltaE_tot,
1427  deltaE_ioni_tot, sigmaDeltaE_ioni_tot,
1428  deltaE_rad_tot, sigmaDeltaE_rad_tot, 0.));
1429 
1430  // direction of plane
1431  Amg::Vector3D dir = wdir / w_tot;
1432  dir = dir / dir.mag();
1433  double norm = dir.perp();
1434  // Rotation matrix representation
1435  Amg::Vector3D colx(-dir.y() / norm, dir.x() / norm, 0);
1436  Amg::Vector3D coly(-dir.x() * dir.z() / norm, -dir.y() * dir.z() / norm,
1437  norm);
1438  Amg::Vector3D colz(dir.x(), dir.y(), dir.z());
1439  // Centre position of the two planes
1440  double halflength2 =
1441  wdist2 / w_tot - (pos - posFirst).mag() * (pos - posFirst).mag();
1442  double halflength = 0.;
1443  if (halflength2 > 0) halflength = sqrt(halflength2);
1444  Amg::Vector3D pos0 = pos - halflength * dir;
1445  Amg::Vector3D posNew = pos + halflength * dir;
1446  if (updateEloss) ATH_MSG_DEBUG("WITH updateEloss");
1447 
1448  ATH_MSG_DEBUG(" WITH aggregation and WITH reposition center planes x "
1449  << pos.x() << " y " << pos.y() << " z " << pos.z()
1450  << " halflength " << halflength << " w_tot " << w_tot
1451  << " X0_tot " << X0_tot);
1452 
1453  Amg::Transform3D surfaceTransformFirst(colx, coly, colz, pos0);
1454  Amg::Transform3D surfaceTransformLast(colx, coly, colz, posNew);
1455  Trk::PlaneSurface* surfFirst =
1456  new Trk::PlaneSurface(surfaceTransformFirst);
1457  Trk::PlaneSurface* surfLast = new Trk::PlaneSurface(surfaceTransformLast);
1458  // calculate TrackParameters at first surface
1459  double qOverP0 = mfirst->trackParameters()->charge() /
1460  (mfirst->trackParameters()->momentum().mag() +
1461  std::fabs(deltaEFirst));
1462  // calculate TrackParameters at last surface
1463  double qOverPNew = mlast->trackParameters()->charge() /
1464  mlast->trackParameters()->momentum().mag();
1465  std::unique_ptr<Trk::TrackParameters> parsFirst =
1466  surfFirst->createUniqueParameters<5, Trk::Charged>(
1467  0., 0., dir.phi(), dir.theta(), qOverP0);
1468  std::unique_ptr<Trk::TrackParameters> parsLast =
1469  surfLast->createUniqueParameters<5, Trk::Charged>(
1470  0., 0., dir.phi(), dir.theta(), qOverPNew);
1471 
1472  Eloss_tot += energyLossNew.deltaE();
1473  if (!threePlanes) {
1474  //
1475  // make two scattering planes and TSOS
1476  //
1477  // prepare for first MaterialEffectsOnTrack with X0 = X0/2
1478  // Eloss = 0 and scattering2 = total2 / 2. depth = 0
1479  auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1480  X0_tot / 2., scatFirst, std::move(energyLoss0), *surfFirst,
1481  meotPattern);
1482  // prepare for second MaterialEffectsOnTrack with X0 = X0/2
1483  // Eloss = Eloss total and scattering2 = total2 / 2. depth
1484  // = 0
1485  auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1486  X0_tot / 2., scatNew,
1487  std::make_unique<Trk::EnergyLoss>(std::move(energyLossNew)),
1488  *surfLast, meotPattern);
1489 
1490  const Trk::TrackStateOnSurface* newTSOSFirst =
1491  new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1492  std::move(meotFirst), typePattern);
1493  auto whichType = (elossFlag != 0) ? typePatternDeposit : typePattern;
1494  const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1495  nullptr, std::move(parsLast),
1496  std::move(meotLast), whichType);
1497 
1498  newTSOSvector.push_back(newTSOSFirst);
1499  newTSOSvector.push_back(newTSOS);
1500  } else {
1501  //
1502  // make three scattering planes and TSOS in Calorimeter
1503  //
1504  auto scatZero = ScatteringAngles(0., 0., 0., 0.);
1505  Amg::Transform3D surfaceTransform(colx, coly, colz, pos);
1506  Trk::PlaneSurface* surf = new Trk::PlaneSurface(surfaceTransform);
1507  std::unique_ptr<Trk::TrackParameters> pars =
1508  surf->createUniqueParameters<5, Trk::Charged>(
1509  0., 0., dir.phi(), dir.theta(), qOverPNew);
1510  // prepare for first MaterialEffectsOnTrack with X0 = X0/2
1511  // Eloss = 0 and scattering2 = total2 / 2. depth = 0
1512  auto meotFirst = std::make_unique<Trk::MaterialEffectsOnTrack>(
1513  X0_tot / 2., scatFirst,
1514  std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.), *surfFirst,
1515  meotPattern);
1516  // prepare for middle MaterialEffectsOnTrack with X0 = 0
1517  // Eloss = ElossNew and scattering2 = 0. depth = 0
1518  auto meot = std::make_unique<Trk::MaterialEffectsOnTrack>(
1519  0., scatZero,
1520  std::make_unique<Trk::EnergyLoss>(std::move(energyLossNew)), *surf,
1521  meotPattern);
1522  // prepare for last MaterialEffectsOnTrack with X0 = X0/2
1523  // Eloss = 0 total and scattering2 = total2 / 2. depth = 0
1524  auto meotLast = std::make_unique<Trk::MaterialEffectsOnTrack>(
1525  X0_tot / 2., scatNew,
1526  std::make_unique<Trk::EnergyLoss>(0., 0., 0., 0.), *surfLast,
1527  meotPattern);
1528  const Trk::TrackStateOnSurface* newTSOSFirst =
1529  new Trk::TrackStateOnSurface(nullptr, std::move(parsFirst),
1530  std::move(meotFirst), typePattern);
1531  const Trk::TrackStateOnSurface* newTSOS = new Trk::TrackStateOnSurface(
1532  nullptr, std::move(pars), std::move(meot), typePatternDeposit);
1533  const Trk::TrackStateOnSurface* newTSOSLast =
1534  new Trk::TrackStateOnSurface(nullptr, std::move(parsLast),
1535  std::move(meotLast), typePattern);
1536  newTSOSvector.push_back(newTSOSFirst);
1537  newTSOSvector.push_back(newTSOS);
1538  newTSOSvector.push_back(newTSOSLast);
1539  }
1540  }
1541  }
1542 
1543  return newTSOSvector;
1544 }
1545 void Trk::GeantFollowerMSHelper::endEvent() {
1546  // fill the validation tree
1547  m_validationTree->Fill();
1548  delete m_parameterCache;
1549  delete m_parameterCacheCov;
1550 
1551  if (m_crossedEntry) {
1552  if (m_parameterCacheEntry) delete m_parameterCacheEntry;
1553  if (m_parameterCacheEntryCov) delete m_parameterCacheEntryCov;
1554  }
1555 }
Trk::GeantFollowerMSHelper::m_useIDExit
bool m_useIDExit
Definition: GeantFollowerMSHelper.h:74
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
make_hlt_rep.pars
pars
Definition: make_hlt_rep.py:90
Trk::EnergyLoss::sigmaMinusDeltaE
double sigmaMinusDeltaE() const
returns the negative side
Trk::TrackStateOnSurface::trackParameters
const TrackParameters * trackParameters() const
return ptr to trackparameters const overload
GeantFollowerMSHelper.h
Trk::TrackStateOnSurface::CaloDeposit
@ CaloDeposit
This TSOS contains a CaloEnergy object.
Definition: TrackStateOnSurface.h:135
EnergyLoss.h
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
ScatteringAngles.h
egammaParameters::depth
@ depth
pointing depth of the shower as calculated in egammaqgcld
Definition: egammaParamDefs.h:276
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
TrackParameters.h
Trk::GeantFollowerMSHelper::endEvent
virtual void endEvent() override
Definition: GeantFollowerMSHelper.cxx:1545
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Trk::locX
@ locX
Definition: ParamDefs.h:37
Trk::GeantFollowerMSHelper::m_extrapolateIncrementally
bool m_extrapolateIncrementally
Definition: GeantFollowerMSHelper.h:71
Trk::ParametersBase::charge
double charge() const
Returns the charge.
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:38
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Trk::oppositeMomentum
@ oppositeMomentum
Definition: PropDirection.h:21
xAOD::deltaPhi
setSAddress setEtaMS setDirPhiMS setDirZMS setBarrelRadius setEndcapAlpha setEndcapRadius setInterceptInner setEtaMap setEtaBin setIsTgcFailure setDeltaPt deltaPhi
Definition: L2StandAloneMuon_v1.cxx:160
Trk::CurvilinearParameters
CurvilinearParametersT< TrackParametersDim, Charged, PlaneSurface > CurvilinearParameters
Definition: Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:29
skel.it
it
Definition: skel.GENtoEVGEN.py:396
Trk::EnergyLoss::sigmaDeltaE
double sigmaDeltaE() const
returns the symmatric error
Trk::z0
@ z0
Definition: ParamDefs.h:64
Trk::ScatteringAngles
represents a deflection of the track caused through multiple scattering in material.
Definition: ScatteringAngles.h:26
IExtrapolator.h
Trk::MaterialEffectsBase::thicknessInX0
double thicknessInX0() const
returns the actually traversed material .
Trk::alongMomentum
@ alongMomentum
Definition: PropDirection.h:20
Trk::EnergyLoss::sigmaRad
double sigmaRad() const
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
Trk::GeantFollowerMSHelper::GeantFollowerMSHelper
GeantFollowerMSHelper(const std::string &, const std::string &, const IInterface *)
Definition: GeantFollowerMSHelper.cxx:27
Trk::MaterialEffectsBase
base class to integrate material effects on Trk::Track in a flexible way.
Definition: MaterialEffectsBase.h:35
Trk::EnergyLoss::meanIoni
double meanIoni() const
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
Trk::ExtrapolationCache::reset
void reset()
ExtrapolationCache.h
Trk::Charged
Definition: Charged.h:27
MaterialEffectsOnTrack.h
Trk::GeantFollowerMSHelper::~GeantFollowerMSHelper
virtual ~GeantFollowerMSHelper()
pdg_comparison.X0
X0
Definition: pdg_comparison.py:314
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
GeoPrimitives.h
Trk::MaterialEffectsOnTrack
represents the full description of deflection and e-loss of a track in material.
Definition: MaterialEffectsOnTrack.h:40
Trk::EnergyLoss::length
double length() const
Trk::ScatteringAngles::sigmaDeltaTheta
double sigmaDeltaTheta() const
returns the
Definition: ScatteringAngles.h:100
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
Trk::GeantFollowerMSHelper::modifyTSOSvector
std::vector< const Trk::TrackStateOnSurface * > modifyTSOSvector(const std::vector< const Trk::TrackStateOnSurface * > &matvec, double scaleX0, double scaleEloss, bool reposition, bool aggregate, bool updateEloss, double caloEnergy, double caloEnergyError, double pCaloEntry, double momentumError, double &Eloss_tot) const
modify TSOS vector
Definition: GeantFollowerMSHelper.cxx:1052
updateDocumentation.wdir
wdir
Definition: updateDocumentation.py:8
Trk::PlaneSurface::createUniqueParameters
std::unique_ptr< ParametersT< DIM, T, PlaneSurface > > createUniqueParameters(double l1, double l2, double phi, double theta, double qop, std::optional< AmgSymMatrix(DIM)> cov=std::nullopt) const
Use the Surface as a ParametersBase constructor, from local parameters.
beamspotman.n
n
Definition: beamspotman.py:731
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
xAOD::covMatrix
covMatrix
Definition: TrackMeasurement_v1.cxx:19
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
perfmonmt-refit.n_tot
n_tot
Definition: perfmonmt-refit.py:104
Trk::EnergyLoss::deltaE
double deltaE() const
returns the
Trk::MaterialEffectsBase::ScatteringEffects
@ ScatteringEffects
contains material effects due to multiple scattering
Definition: MaterialEffectsBase.h:45
Trk::GeantFollowerMSHelper::m_extrapolateDirectly
bool m_extrapolateDirectly
Definition: GeantFollowerMSHelper.h:70
Trk::GeantFollowerMSHelper::initialize
virtual StatusCode initialize() override
Definition: GeantFollowerMSHelper.cxx:66
Trk::CurvilinearParametersT
Definition: CurvilinearParametersT.h:48
Trk::muon
@ muon
Definition: ParticleHypothesis.h:28
Trk::EnergyLoss::meanRad
double meanRad() const
python.TriggerConfigAccess.maxsize
maxsize
Definition: TriggerConfigAccess.py:51
beamspotman.dir
string dir
Definition: beamspotman.py:623
Trk::CurvilinearParametersT::associatedSurface
virtual const S & associatedSurface() const override final
Access to the Surface method.
Trk::TrackStateOnSurface
represents the track state (measurement, material, fit parameters and quality) at a surface.
Definition: TrackStateOnSurface.h:71
Trk::d0
@ d0
Definition: ParamDefs.h:63
Trk::TrackStateOnSurface::InertMaterial
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
Definition: TrackStateOnSurface.h:105
Trk::EnergyLoss
This class describes energy loss material effects in the ATLAS tracking EDM.
Definition: EnergyLoss.h:34
charge
double charge(const T &p)
Definition: AtlasPID.h:756
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::MaterialEffectsBase::EnergyLossEffects
@ EnergyLossEffects
contains energy loss corrections
Definition: MaterialEffectsBase.h:48
Trk::Surface::clone
virtual Surface * clone() const =0
Implicit constructor - uses the copy constructor.
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::GeantFollowerMSHelper::trackParticle
virtual void trackParticle(const G4ThreeVector &pos, const G4ThreeVector &mom, int pdg, double charge, float t, float X0) override
Definition: GeantFollowerMSHelper.cxx:284
Trk::GeantFollowerMSHelper::m_extrapolator
ToolHandle< IExtrapolator > m_extrapolator
Definition: GeantFollowerMSHelper.h:66
Trk::MaterialEffectsOnTrack::energyLoss
const EnergyLoss * energyLoss() const
returns the energy loss object.
MAXPROBES
#define MAXPROBES
Definition: ActsGeantFollowerHelper.h:22
Trk::ExtrapolationCache::x0tot
double x0tot() const
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
Trk::addNoise
@ addNoise
Definition: MaterialUpdateMode.h:19
PlaneSurface.h
Trk::MaterialEffectsOnTrack::scatteringAngles
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
Trk::GeantFollowerMSHelper::m_useCovMatrix
bool m_useCovMatrix
Definition: GeantFollowerMSHelper.h:73
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
Trk::ScatteringAngles::sigmaDeltaPhi
double sigmaDeltaPhi() const
returns the
Definition: ScatteringAngles.h:94
Trk::TrackStateOnSurface::Scatterer
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
Definition: TrackStateOnSurface.h:113
Trk::phi
@ phi
Definition: ParamDefs.h:75
Trk::GeantFollowerMSHelper::m_speedup
bool m_speedup
Definition: GeantFollowerMSHelper.h:72
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
mag2
Scalar mag2() const
mag2 method - forward to squaredNorm()
Definition: AmgMatrixBasePlugin.h:31
StoreGateSvc.h
Trk::ExtrapolationCache
Definition: ExtrapolationCache.h:26
Trk::MaterialEffectsBase::associatedSurface
const Surface & associatedSurface() const
returns the surface to which these m.eff. are associated.
Trk::EnergyLoss::sigmaPlusDeltaE
double sigmaPlusDeltaE() const
returns the positive side
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
Trk::GeantFollowerMSHelper::beginEvent
virtual void beginEvent() override
Definition: GeantFollowerMSHelper.cxx:246
Trk::ExtrapolationCache::eloss
const EnergyLoss * eloss() const
Trk::EnergyLoss::sigmaIoni
double sigmaIoni() const
Trk::GeantFollowerMSHelper::finalize
virtual StatusCode finalize() override
Definition: GeantFollowerMSHelper.cxx:242
Generated on Mon Dec 23 2024 21:10:46 for ATLAS Offline Software by   doxygen 1.8.18