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