ATLAS Offline Software
Loading...
Searching...
No Matches
Trk::GeantFollowerMSHelper Class Reference

#include <GeantFollowerMSHelper.h>

Inheritance diagram for Trk::GeantFollowerMSHelper:
Collaboration diagram for Trk::GeantFollowerMSHelper:

Classes

struct  TreeData
 Ntuple variables : initial parameters Split this out into a separate, dynamically-allocated block. More...

Public Member Functions

 GeantFollowerMSHelper (const std::string &, const std::string &, const IInterface *)
virtual ~GeantFollowerMSHelper ()
virtual StatusCode initialize () override
virtual StatusCode finalize () override
virtual void beginEvent () override
virtual void trackParticle (const G4ThreeVector &pos, const G4ThreeVector &mom, int pdg, double charge, float t, float X0) override
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
virtual void endEvent () override

Static Public Attributes

static constexpr int MAXPROBES {50000}

Private Attributes

ToolHandle< IExtrapolatorm_extrapolator
PublicToolHandle< IEnergyLossUpdatorm_elossupdator {this,"EnergyLossUpdator","Trk::EnergyLossUpdator/AtlasEnergyLossUpdator",""}
bool m_extrapolateDirectly
bool m_extrapolateIncrementally
bool m_speedup
bool m_useCovMatrix
bool m_useIDExit
const TrackParametersm_parameterCache
const TrackParametersm_parameterCacheCov
const TrackParametersm_parameterCacheEntry
const TrackParametersm_parameterCacheEntryCov
float m_tX0Cache
bool m_crossedEntry
bool m_exitLayer
PlaneSurface m_destinationSurface
std::string m_validationTreeName
 validation tree name - to be acessed by this from root
std::string m_validationTreeDescription
 validation tree description - second argument in TTree
std::string m_validationTreeFolder
 stream/folder to for the TTree to be written out
TTree * m_validationTree
 Root Validation Tree.
std::unique_ptr< TreeDatam_treeData

Detailed Description

Definition at line 32 of file GeantFollowerMSHelper.h.

Constructor & Destructor Documentation

◆ GeantFollowerMSHelper()

Trk::GeantFollowerMSHelper::GeantFollowerMSHelper ( const std::string & t,
const std::string & n,
const IInterface * p )

Definition at line 27 of file GeantFollowerMSHelper.cxx.

30 : base_class(t, n, p),
34 m_speedup(false),
35 m_useCovMatrix(true),
36 m_useIDExit(false),
37 m_parameterCache(nullptr),
38 m_parameterCacheCov(nullptr),
39 m_parameterCacheEntry(nullptr),
41 m_tX0Cache(0.),
42 m_crossedEntry(false),
43 m_exitLayer(false),
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}
const TrackParameters * m_parameterCacheCov
std::string m_validationTreeFolder
stream/folder to for the TTree to be written out
const TrackParameters * m_parameterCache
TTree * m_validationTree
Root Validation Tree.
const TrackParameters * m_parameterCacheEntry
std::string m_validationTreeName
validation tree name - to be acessed by this from root
std::string m_validationTreeDescription
validation tree description - second argument in TTree
ToolHandle< IExtrapolator > m_extrapolator
const TrackParameters * m_parameterCacheEntryCov

◆ ~GeantFollowerMSHelper()

Trk::GeantFollowerMSHelper::~GeantFollowerMSHelper ( )
virtualdefault

Member Function Documentation

◆ beginEvent()

void Trk::GeantFollowerMSHelper::beginEvent ( )
overridevirtual

Definition at line 246 of file GeantFollowerMSHelper.cxx.

246 {
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}
std::unique_ptr< TreeData > m_treeData

◆ endEvent()

void Trk::GeantFollowerMSHelper::endEvent ( )
overridevirtual

Definition at line 1545 of file GeantFollowerMSHelper.cxx.

1545 {
1546 // fill the validation tree
1547 m_validationTree->Fill();
1548 delete m_parameterCache;
1549 delete m_parameterCacheCov;
1550
1551 if (m_crossedEntry) {
1554 }
1555}

◆ finalize()

StatusCode Trk::GeantFollowerMSHelper::finalize ( )
overridevirtual

Definition at line 242 of file GeantFollowerMSHelper.cxx.

242 {
243 return StatusCode::SUCCESS;
244}

◆ initialize()

StatusCode Trk::GeantFollowerMSHelper::initialize ( )
overridevirtual

Definition at line 66 of file GeantFollowerMSHelper.cxx.

66 {
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(),
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 ) {
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()) {
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}
#define ATH_MSG_ERROR(x)
#define ATH_MSG_INFO(x)
PublicToolHandle< IEnergyLossUpdator > m_elossupdator

◆ modifyTSOSvector()

std::vector< const Trk::TrackStateOnSurface * > Trk::GeantFollowerMSHelper::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 at line 1052 of file GeantFollowerMSHelper.cxx.

1056 {
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);
1126 // meotPattern.set(Trk::MaterialEffectsBase::FittedMaterialEffects);
1127
1128 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1129 typePattern(0);
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}
Scalar deltaPhi(const MatrixBase< Derived > &vec) const
Scalar mag() const
mag method
Scalar mag2() const
mag2 method - forward to squaredNorm()
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
double meanRad() const
double length() const
double sigmaPlusDeltaE() const
returns the positive side
double sigmaMinusDeltaE() const
returns the negative side
double sigmaIoni() const
double meanIoni() const
double sigmaDeltaE() const
returns the symmatric error
double sigmaRad() const
double deltaE() const
returns the
@ ScatteringEffects
contains material effects due to multiple scattering
@ EnergyLossEffects
contains energy loss corrections
const Surface & associatedSurface() const
returns the surface to which these m.eff. are associated.
const EnergyLoss * energyLoss() const
returns the energy loss object.
const ScatteringAngles * scatteringAngles() const
returns the MCS-angles object.
const Amg::Vector3D & momentum() const
Access method for the momentum.
double charge() const
Returns the charge.
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.
double sigmaDeltaPhi() const
returns the
double sigmaDeltaTheta() const
returns the
const TrackParameters * trackParameters() const
return ptr to trackparameters const overload
@ InertMaterial
This represents inert material, and so will contain MaterialEffectsBase.
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
@ CaloDeposit
This TSOS contains a CaloEnergy object.
std::string depth
tag string for intendation
Definition fastadd.cxx:46
Eigen::Affine3d Transform3D
Eigen::Matrix< double, 3, 1 > Vector3D

◆ trackParticle()

void Trk::GeantFollowerMSHelper::trackParticle ( const G4ThreeVector & pos,
const G4ThreeVector & mom,
int pdg,
double charge,
float t,
float X0 )
overridevirtual

Definition at line 284 of file GeantFollowerMSHelper.cxx.

287 {
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
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
368 //m_parameterCache = new Trk::CurvilinearParameters(npos, nmom, charge);
369 AmgSymMatrix(5) covMatrix;
370 covMatrix.setZero();
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
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 =
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) {
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,
491 : m_extrapolator->extrapolate(
492 ctx, *m_parameterCacheCov, destinationSurface,
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,
503 : m_extrapolator->extrapolate(
504 ctx, *m_parameterCache, destinationSurface,
506 Trk::addNoise, &extrapolationCache);
507
508 } else {
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,
554 : m_extrapolator->extrapolate(
555 ctx, *m_parameterCacheEntryCov, destinationSurface,
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 ");
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,
570 : m_extrapolator->extrapolate(
571 ctx, *m_parameterCacheEntry, destinationSurface,
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,
584 : m_extrapolator->extrapolate(
585 ctx, *m_parameterCacheEntry, destinationSurface,
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 =
597 ? m_extrapolator->extrapolateDirectly(
598 ctx, *trkParameters, m_destinationSurface,
600 : m_extrapolator->extrapolate(
601 ctx, *trkParameters, m_destinationSurface,
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.)
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,
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) {
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()) {
733 " G4 extrapolateM failed with covariance matrix to Muon Entry or ID Exit ");
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 {
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()) {
760 " G4 extrapolateM failed with covariance matrix to Muon or Calo Exit ");
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 {
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}
Scalar perp() const
perp method - perpendicular length
double charge(const T &p)
Definition AtlasPID.h:997
#define AmgSymMatrix(dim)
if(febId1==febId2)
virtual const S & associatedSurface() const override final
Access to the Surface method.
const EnergyLoss * eloss() const
double x0tot() const
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
double thicknessInX0() const
returns the actually traversed material .
@ oppositeMomentum
@ alongMomentum
CurvilinearParametersT< TrackParametersDim, Charged, PlaneSurface > CurvilinearParameters
@ locY
local cartesian
Definition ParamDefs.h:38
@ x
Definition ParamDefs.h:55
@ locX
Definition ParamDefs.h:37
@ z
global position (cartesian)
Definition ParamDefs.h:57
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ y
Definition ParamDefs.h:56
@ phi
Definition ParamDefs.h:75
@ d0
Definition ParamDefs.h:63
@ z0
Definition ParamDefs.h:64

Member Data Documentation

◆ m_crossedEntry

bool Trk::GeantFollowerMSHelper::m_crossedEntry
private

Definition at line 82 of file GeantFollowerMSHelper.h.

◆ m_destinationSurface

PlaneSurface Trk::GeantFollowerMSHelper::m_destinationSurface
private

Definition at line 84 of file GeantFollowerMSHelper.h.

◆ m_elossupdator

PublicToolHandle<IEnergyLossUpdator> Trk::GeantFollowerMSHelper::m_elossupdator {this,"EnergyLossUpdator","Trk::EnergyLossUpdator/AtlasEnergyLossUpdator",""}
private

Definition at line 67 of file GeantFollowerMSHelper.h.

68{this,"EnergyLossUpdator","Trk::EnergyLossUpdator/AtlasEnergyLossUpdator",""};

◆ m_exitLayer

bool Trk::GeantFollowerMSHelper::m_exitLayer
private

Definition at line 83 of file GeantFollowerMSHelper.h.

◆ m_extrapolateDirectly

bool Trk::GeantFollowerMSHelper::m_extrapolateDirectly
private

Definition at line 70 of file GeantFollowerMSHelper.h.

◆ m_extrapolateIncrementally

bool Trk::GeantFollowerMSHelper::m_extrapolateIncrementally
private

Definition at line 71 of file GeantFollowerMSHelper.h.

◆ m_extrapolator

ToolHandle<IExtrapolator> Trk::GeantFollowerMSHelper::m_extrapolator
private

Definition at line 66 of file GeantFollowerMSHelper.h.

◆ m_parameterCache

const TrackParameters* Trk::GeantFollowerMSHelper::m_parameterCache
private

Definition at line 76 of file GeantFollowerMSHelper.h.

◆ m_parameterCacheCov

const TrackParameters* Trk::GeantFollowerMSHelper::m_parameterCacheCov
private

Definition at line 77 of file GeantFollowerMSHelper.h.

◆ m_parameterCacheEntry

const TrackParameters* Trk::GeantFollowerMSHelper::m_parameterCacheEntry
private

Definition at line 78 of file GeantFollowerMSHelper.h.

◆ m_parameterCacheEntryCov

const TrackParameters* Trk::GeantFollowerMSHelper::m_parameterCacheEntryCov
private

Definition at line 79 of file GeantFollowerMSHelper.h.

◆ m_speedup

bool Trk::GeantFollowerMSHelper::m_speedup
private

Definition at line 72 of file GeantFollowerMSHelper.h.

◆ m_treeData

std::unique_ptr<TreeData> Trk::GeantFollowerMSHelper::m_treeData
private

Definition at line 174 of file GeantFollowerMSHelper.h.

◆ m_tX0Cache

float Trk::GeantFollowerMSHelper::m_tX0Cache
private

Definition at line 80 of file GeantFollowerMSHelper.h.

◆ m_useCovMatrix

bool Trk::GeantFollowerMSHelper::m_useCovMatrix
private

Definition at line 73 of file GeantFollowerMSHelper.h.

◆ m_useIDExit

bool Trk::GeantFollowerMSHelper::m_useIDExit
private

Definition at line 74 of file GeantFollowerMSHelper.h.

◆ m_validationTree

TTree* Trk::GeantFollowerMSHelper::m_validationTree
private

Root Validation Tree.

Definition at line 91 of file GeantFollowerMSHelper.h.

◆ m_validationTreeDescription

std::string Trk::GeantFollowerMSHelper::m_validationTreeDescription
private

validation tree description - second argument in TTree

Definition at line 88 of file GeantFollowerMSHelper.h.

◆ m_validationTreeFolder

std::string Trk::GeantFollowerMSHelper::m_validationTreeFolder
private

stream/folder to for the TTree to be written out

Definition at line 89 of file GeantFollowerMSHelper.h.

◆ m_validationTreeName

std::string Trk::GeantFollowerMSHelper::m_validationTreeName
private

validation tree name - to be acessed by this from root

Definition at line 87 of file GeantFollowerMSHelper.h.

◆ MAXPROBES

int Trk::GeantFollowerMSHelper::MAXPROBES {50000}
staticconstexpr

Definition at line 37 of file GeantFollowerMSHelper.h.

37{50000};

The documentation for this class was generated from the following files: