ATLAS Offline Software
Loading...
Searching...
No Matches
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
27Trk::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
62Trk::GeantFollowerMSHelper::~GeantFollowerMSHelper() = default;
63
64// Athena standard methods
65// initialize
66StatusCode 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
242StatusCode Trk::GeantFollowerMSHelper::finalize() {
243 return StatusCode::SUCCESS;
244}
245
246void 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
284void 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
1051std::vector<const Trk::TrackStateOnSurface*>
1052Trk::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}
1545void 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}
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:997
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:242
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:1052
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:27
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:284
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:246
Trk::GeantFollowerMSHelper::m_parameterCacheEntryCov
const TrackParameters * m_parameterCacheEntryCov
Definition GeantFollowerMSHelper.h:79
Trk::GeantFollowerMSHelper::initialize
virtual StatusCode initialize() override
Definition GeantFollowerMSHelper.cxx:66
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, 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 Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/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
Generated on for ATLAS Offline Software by doxygen 1.14.0