2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
9 #include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
11 #include "TrkParameters/TrackParameters.h"
12 #include "TrkNeutralParameters/NeutralParameters.h"
13 #include "TrkValInterfaces/IPositionMomentumWriter.h"
14 #include "TrkEventPrimitives/ParticleHypothesis.h"
15 #include "TrkEventPrimitives/PdgToParticleHypothesis.h"
16 #include "TrkExUtils/MaterialInteraction.h"
17 #include "TrkGeometry/Layer.h"
18 #include "InDetIdentifier/PixelID.h"
19 #include "InDetIdentifier/SCT_ID.h"
20 #include "HGTD_Identifier/HGTD_ID.h"
21 #include "AtlasDetDescr/AtlasDetectorID.h"
22 #include "Identifier/Identifier.h"
23 #include "Identifier/IdentifierHash.h"
24 #include "InDetReadoutGeometry/SiDetectorElement.h"
26 template < class T, class P > StatusCode Trk::ExtrapolationEngineTest::runTestT() {
27 // ----------------- creation of the surfaces & track parameters -------------
28 for (size_t it = 0; it < TrkExUnitTestBase::m_numTests; ++it) {
30 ATH_MSG_DEBUG("===> starting test " << it << " <<===");
31 // create the curvilinear parameters
32 double eta = (m_scanMode) ? m_currentEta : m_etaMin + (m_etaMax - m_etaMin) *
33 Trk::TrkExUnitTestBase::m_flatDist->shoot();
34 double theta = 2. * atan(exp(-eta));
35 double phi = (m_scanMode) ? m_currentPhi : m_phiMin + (m_phiMax - m_phiMin) *
36 Trk::TrkExUnitTestBase::m_flatDist->shoot();
37 double pt = m_ptMin + (m_ptMax - m_ptMin) * Trk::TrkExUnitTestBase::m_flatDist->shoot();
38 double p = pt / sin(theta);
39 double q = m_splitCharge ? m_charge *
40 -1. : (m_parametersMode ? (Trk::TrkExUnitTestBase::m_flatDist->shoot() > 0.5 ? 1. : -1) : 1.); // charge
44 // initializa the validation variables
56 // material collection
57 m_materialThicknessInX0 = 0.;
58 m_materialThicknessInL0 = 0.;
59 m_materialThicknessZARho = 0.;
60 m_materialEmulatedIonizationLoss = 0.;
61 m_materialThicknessInX0Sensitive = 0.;
62 m_materialThicknessInX0Passive = 0.;
63 m_materialThicknessInX0Boundary = 0.;
65 m_materialThicknessInX0Bwd = 0.;
66 m_materialThicknessInL0Bwd = 0.;
68 m_materialThicknessInX0Cylinder = 0.;
69 m_materialThicknessInX0Disc = 0.;
70 m_materialThicknessInX0Plane = 0.;
72 if (m_collectMaterial) {
73 m_materialThicknessInX0Accumulated->clear();
74 m_materialThicknessInX0Steps->clear();
75 m_materialThicknessInX0Steps->clear();
76 m_materialPositionX->clear();
77 m_materialPositionY->clear();
78 m_materialPositionZ->clear();
79 m_materialPositionR->clear();
80 m_materialPositionP->clear();
81 m_materialPositionPt->clear();
82 m_materialScaling->clear();
83 m_stepDirection->clear();
86 if (m_collectSensitive) {
87 m_sensitiveLayerIndex->clear();
88 m_sensitiveSurfaceType->clear();
89 m_sensitiveCenterPosX->clear();
90 m_sensitiveCenterPosY->clear();
91 m_sensitiveCenterPosZ->clear();
92 m_sensitiveCenterPosR->clear();
93 m_sensitiveCenterPosPhi->clear();
94 m_sensitiveLocalPosX->clear();
95 m_sensitiveLocalPosY->clear();
96 m_sensitiveLocalPosR->clear();
97 m_sensitiveLocalPosPhi->clear();
98 m_sensitiveDetector->clear();
99 m_sensitiveIsInnermost->clear();
100 m_sensitiveIsNextToInnermost->clear();
101 m_sensitiveBarrelEndcap->clear();
102 m_sensitiveLayerDisc->clear();
103 m_sensitiveEtaModule->clear();
104 m_sensitivePhiModule->clear();
105 m_sensitiveSide->clear();
106 m_sensitiveNearBondGap->clear();
107 m_sensitiveisInside->clear();
108 m_sensitiveisInsideBound->clear();
109 m_materialtInX0AccumulatedUpTo->clear();
112 for (size_t ip = 0; ip < m_parameterNames.size(); ++ip) {
114 m_pPositionX[ip]->clear();
115 m_pPositionY[ip]->clear();
116 m_pPositionZ[ip]->clear();
117 m_pPositionR[ip]->clear();
119 m_pTheta[ip]->clear();
125 Amg::Vector3D momentum(p* sin(theta)* cos(phi), p* sin(theta)* sin(phi), p* cos(theta));
126 // create the start parameters
128 m_smearProductionVertex ? (m_smearFlatOriginT ? (m_d0Min + (m_d0Max - m_d0Min) *
129 Trk::TrkExUnitTestBase::m_flatDist->shoot()) : Trk::
130 TrkExUnitTestBase::m_gaussDist->shoot() * m_sigmaOriginT) : 0.;
132 m_smearProductionVertex ? (m_smearFlatOriginZ ? (m_z0Min + (m_z0Max - m_z0Min) *
133 Trk::TrkExUnitTestBase::m_flatDist->shoot()) : Trk::
134 TrkExUnitTestBase::m_gaussDist->shoot() * m_sigmaOriginZ) : 0.;
135 // create the start parameters
136 P startParameters(d0, z0, phi, theta, q / p, Trk::PerigeeSurface(Amg::Vector3D(0., 0., 0.)));
138 m_startPositionX = startParameters.position().x();
139 m_startPositionY = startParameters.position().y();
140 m_startPositionZ = startParameters.position().z();
141 m_startPositionR = startParameters.position().perp();
144 m_startTheta = theta;
149 // the step surfaces for stepwise extrapolation afterwards
150 std::vector < const Trk::Surface * > stepSurfaces;
152 // position / momentum writer
153 if (!m_posmomWriter.empty() && m_writeTTree) {
155 Trk::PdgToParticleHypothesis pdgFromHypothesis;
157 double mass = Trk::ParticleMasses::mass[m_particleHypothesis];
158 int pdg = pdgFromHypothesis.convert(static_cast<Trk::ParticleHypothesis>(m_particleHypothesis.value()), m_charge);
159 m_posmomWriter->initializeTrack(startParameters.position(),
160 startParameters.momentum(),
164 // setup the extrapolation how you'd like it
165 Trk::ExtrapolationCell < T > ecc(startParameters);
166 ecc.setParticleHypothesis(static_cast<Trk::ParticleHypothesis>(m_particleHypothesis.value()));
167 ecc.addConfigurationMode(Trk::ExtrapolationMode::StopAtBoundary);
168 ecc.addConfigurationMode(Trk::ExtrapolationMode::FATRAS);
169 if (m_collectSensitive) ecc.addConfigurationMode(Trk::ExtrapolationMode::CollectSensitive);
170 if (m_collectPassive) ecc.addConfigurationMode(Trk::ExtrapolationMode::CollectPassive);
171 if (m_collectBoundary) ecc.addConfigurationMode(Trk::ExtrapolationMode::CollectBoundary);
172 if (m_collectMaterial) ecc.addConfigurationMode(Trk::ExtrapolationMode::CollectMaterial);
173 // force a stop in the extrapoaltion mode
174 if (m_pathLimit > 0.) {
175 ecc.pathLimit = m_pathLimit;
176 ecc.addConfigurationMode(Trk::ExtrapolationMode::StopWithPathLimit);
179 ATH_MSG_DEBUG("===> forward extrapolation - collecting information <<===");
180 // call the extrapolation engine
181 Trk::ExtrapolationCode eCode = m_extrapolationEngine->extrapolate(ecc);
182 // save the success code
183 m_endSuccessful = int(eCode.code);
184 // end the parameters if there
185 if (eCode.isSuccess() && ecc.endParameters && m_writeTTree) {
186 m_endPositionX = ecc.endParameters->position().x();
187 m_endPositionY = ecc.endParameters->position().y();
188 m_endPositionZ = ecc.endParameters->position().z();
189 m_endPositionR = ecc.endParameters->position().perp();
190 m_endPhi = ecc.endParameters->momentum().phi();
191 m_endEta = ecc.endParameters->momentum().eta();
192 m_endTheta = ecc.endParameters->momentum().theta();
193 m_endP = ecc.endParameters->momentum().mag();
194 m_endPt = ecc.endParameters->momentum().perp();
195 m_endPathLength = ecc.pathLength;
196 // fill in the step parameters
197 if (!fillStepInformationT < T, P > (ecc, 1, stepSurfaces).isSuccess()) ATH_MSG_VERBOSE(
198 "Somthing went wrong with recording the step information.");
199 // for the memory cleanup
200 const T* feParameters = ecc.endParameters;
201 // do the back extrapolation
202 if (m_backExtrapolation) {
203 // let's restart at the destination
204 ecc.restartAtDestination();
206 ATH_MSG_DEBUG("===> backward extrapolation - collecting information <<===");
207 // call the extrapolation engine
208 Trk::ExtrapolationCode eCodeBwd = m_extrapolationEngine->extrapolate(ecc, &startParameters.associatedSurface());
209 // safe the success code
210 m_backSuccessful = int(eCodeBwd.code);
211 // end the parameters if there
212 if (eCodeBwd.isSuccess() && ecc.endParameters) {
213 m_backPositionX = ecc.endParameters->position().x();
214 m_backPositionY = ecc.endParameters->position().y();
215 m_backPositionZ = ecc.endParameters->position().z();
216 m_backPositionR = ecc.endParameters->position().perp();
217 m_backPhi = ecc.endParameters->momentum().phi();
218 m_backEta = ecc.endParameters->momentum().eta();
219 m_backTheta = ecc.endParameters->momentum().theta();
220 m_backP = ecc.endParameters->momentum().mag();
221 m_backPt = ecc.endParameters->momentum().perp();
222 // fill in the step parameters
223 if (!fillStepInformationT < T, P > (ecc, -1, stepSurfaces).isSuccess()) ATH_MSG_VERBOSE(
224 "Somthing went wrong with recording the step information.");
226 delete ecc.endParameters;
229 // next test - run stepwise extrapolation from one surface to another - this is for Kalman filtering
230 if (m_stepwiseExtrapolation && stepSurfaces.size()) {
231 // the ExtrapolationCell for stepwise
232 Trk::ExtrapolationCell < T > eccSteps(startParameters);
233 eccSteps.setParticleHypothesis(static_cast<Trk::ParticleHypothesis>(m_particleHypothesis.value()));
234 eccSteps.addConfigurationMode(Trk::ExtrapolationMode::StopAtBoundary);
235 if (m_collectMaterial) eccSteps.addConfigurationMode(Trk::ExtrapolationMode::CollectMaterial);
236 std::vector < const T * > stepParameters;
237 // prepare the extrapolation Code
238 Trk::ExtrapolationCode eCodeSteps(Trk::ExtrapolationCode::Unset);
240 double pathLength = 0.;
241 double materialInX0 = 0.;
242 // we step through the surfaces and go from one to the next
243 ATH_MSG_DEBUG("===> stepwise extrapolation - collecting information <<===");
244 for (auto& sf : stepSurfaces) {
245 ATH_MSG_VERBOSE("===> step <===");
246 // do the extrapolation - will create a memory leak, but
247 eCodeSteps = m_extrapolationEngine->extrapolate(eccSteps, sf);
248 if (!eCodeSteps.isSuccess()) {
249 m_endStepSuccessful = int(eCodeSteps.code);
252 // accumulate the path length
253 pathLength += eccSteps.pathLength;
254 // loop over the collected information for material collection
255 for (auto& es : eccSteps.extrapolationSteps) {
256 if (es.material) materialInX0 += es.material->thicknessInX0() * es.materialScaling;
259 eccSteps.restartAtDestination();
260 // for the memory cleanup
261 stepParameters.push_back(eccSteps.endParameters);
263 // the final extrapolation
264 if (eCodeSteps.isSuccess()) {
265 // the final one is still missing
266 ATH_MSG_VERBOSE("===> final step <===");
267 eCodeSteps = m_extrapolationEngine->extrapolate(eccSteps, &feParameters->associatedSurface());
268 // what is the final extrapolation code
269 m_endStepSuccessful = int(eCodeSteps.code);
271 if (eCodeSteps.isSuccess()) {
272 // and now the final step
273 pathLength += eccSteps.pathLength;
275 // loop over the collected information for material collection
276 for (auto& es : eccSteps.extrapolationSteps) {
277 if (es.material) materialInX0 += es.material->thicknessInX0() * es.materialScaling;
279 // stepwise end parameters
280 m_endStepPositionX = eccSteps.endParameters->position().x();
281 m_endStepPositionY = eccSteps.endParameters->position().y();
282 m_endStepPositionZ = eccSteps.endParameters->position().z();
283 m_endStepPositionR = eccSteps.endParameters->position().perp();
284 m_endStepPhi = eccSteps.endParameters->momentum().phi();
285 m_endStepEta = eccSteps.endParameters->momentum().eta();
286 m_endStepTheta = eccSteps.endParameters->momentum().theta();
287 m_endStepP = eccSteps.endParameters->momentum().mag();
288 m_endStepPt = eccSteps.endParameters->momentum().perp();
289 m_endStepPathLength = pathLength;
290 m_endStepThicknessInX0 = materialInX0;
294 for (auto& sp : stepParameters)
298 // now clean up the first end parameters
302 "Extrapolation was not Successful - code : " << eCode.toString() << " -> Printing start parameters.");
303 ATH_MSG_WARNING(" -> start x = " << m_startPositionX);
304 ATH_MSG_WARNING(" -> start y = " << m_startPositionY);
305 ATH_MSG_WARNING(" -> start z = " << m_startPositionZ);
306 ATH_MSG_WARNING(" -> start r = " << m_startPositionR);
307 ATH_MSG_WARNING(" -> start phi = " << m_startPhi);
308 ATH_MSG_WARNING(" -> start eta = " << m_startEta);
309 ATH_MSG_WARNING(" -> start theta = " << m_startTheta);
310 ATH_MSG_WARNING(" -> start pt = " << m_startPt);
311 ATH_MSG_WARNING(" -> start p = " << m_startP);
314 // use the dedicated position / momentum writeter
315 if (!m_posmomWriter.empty()) m_posmomWriter->finalizeTrack();
318 return StatusCode::SUCCESS;
321 template < class T, class P > StatusCode Trk::ExtrapolationEngineTest::fillStepInformationT(
322 Trk::ExtrapolationCell < T >& ecc, int fwdBwd, std::vector < const Trk::Surface* >& stepSurfaces) {
324 // loop over the collected information
325 for (auto& es : ecc.extrapolationSteps) {
326 // continue if we have parameters
327 const T* parameters = es.parameters;
329 // collect for stepwise extrapolation
330 if (fwdBwd && m_stepwiseExtrapolation && parameters->associatedSurface().associatedDetectorElement()) {
331 stepSurfaces.push_back(¶meters->associatedSurface());
333 // there are parameters assigned, so they need to be either
334 // sensitive -> ip = 0
336 // boundary -> ip = 2
338 if (es.stepConfiguration.checkMode(Trk::ExtrapolationMode::CollectPassive)) {
340 } else if (es.stepConfiguration.checkMode(Trk::ExtrapolationMode::CollectBoundary)) {
344 // use the dedicated position / momentum writeter
345 if (!m_posmomWriter.empty()) {
346 m_posmomWriter->recordTrackState(parameters->position(), parameters->momentum());
348 // fill the parameters
349 m_pPositionX[ip]->push_back(parameters->position().x());
350 m_pPositionY[ip]->push_back(parameters->position().y());
351 m_pPositionZ[ip]->push_back(parameters->position().z());
352 m_pPositionR[ip]->push_back(parameters->position().perp());
353 m_pPhi[ip]->push_back(parameters->momentum().phi());
354 m_pTheta[ip]->push_back(parameters->momentum().eta());
355 m_pEta[ip]->push_back(parameters->momentum().theta());
356 m_pP[ip]->push_back(parameters->momentum().mag());
357 m_pPt[ip]->push_back(parameters->momentum().perp());
359 if (ip == 0 && m_collectSensitive) {
360 // fill the layer index and the local variables
361 int idx = parameters->associatedSurface().associatedLayer() ?
362 parameters->associatedSurface().associatedLayer()->layerIndex().value() : -1;
363 m_sensitiveLayerIndex->push_back(idx);
364 int type = int(parameters->associatedSurface().type());
365 m_sensitiveSurfaceType->push_back(type);
366 m_sensitiveLocalPosX->push_back(parameters->localPosition()[Trk::loc1]);
367 m_sensitiveLocalPosY->push_back(parameters->localPosition()[Trk::loc2]);
369 if (parameters->associatedSurface().associatedDetectorElement()) {
370 Identifier id = parameters->associatedSurface().associatedDetectorElement()->identify();
372 m_sensitiveCenterPosX->push_back(parameters->associatedSurface().associatedDetectorElement()->center().x());
373 m_sensitiveCenterPosY->push_back(parameters->associatedSurface().associatedDetectorElement()->center().y());
374 m_sensitiveCenterPosZ->push_back(parameters->associatedSurface().associatedDetectorElement()->center().z());
375 m_sensitiveCenterPosR->push_back(parameters->associatedSurface().associatedDetectorElement()->center().perp());
376 m_sensitiveCenterPosPhi->push_back(parameters->associatedSurface().associatedDetectorElement()->center().phi());
379 if (m_idHelper->is_pixel(id)) {
381 m_sensitiveDetector->push_back(0);
382 m_sensitiveIsInnermost->push_back(m_pixel_ID->layer_disk(id)==0);
383 m_sensitiveIsNextToInnermost->push_back(m_pixel_ID->layer_disk(id)==1);
384 m_sensitiveBarrelEndcap->push_back(m_pixel_ID->barrel_ec(id));
385 m_sensitiveLayerDisc->push_back(m_pixel_ID->layer_disk(id));
386 m_sensitiveEtaModule->push_back(m_pixel_ID->eta_module(id));
387 m_sensitivePhiModule->push_back(m_pixel_ID->phi_module(id));
388 m_sensitiveSide->push_back(0);
389 } else if (m_idHelper->is_sct(id)) {
391 m_sensitiveDetector->push_back(1);
392 m_sensitiveIsInnermost->push_back(0);
393 m_sensitiveIsNextToInnermost->push_back(0);
394 m_sensitiveBarrelEndcap->push_back(m_sct_ID->barrel_ec(id));
395 m_sensitiveLayerDisc->push_back(m_sct_ID->layer_disk(id));
396 m_sensitiveEtaModule->push_back(m_sct_ID->eta_module(id));
397 m_sensitivePhiModule->push_back(m_sct_ID->phi_module(id));
398 m_sensitiveSide->push_back(m_sct_ID->side(id));
399 } else if (m_useHGTD and m_idHelper->is_hgtd(id)) {
401 m_sensitiveDetector->push_back(2);
402 m_sensitiveIsInnermost->push_back(0);
403 m_sensitiveIsNextToInnermost->push_back(0);
404 m_sensitiveBarrelEndcap->push_back(m_hgtd_ID->endcap(id));
405 m_sensitiveLayerDisc->push_back(m_hgtd_ID->layer(id));
406 m_sensitiveEtaModule->push_back(m_hgtd_ID->eta_module(id));
407 m_sensitivePhiModule->push_back(m_hgtd_ID->phi_module(id));
408 m_sensitiveSide->push_back(0);
413 const InDetDD::SiDetectorElement* siElement = dynamic_cast<const InDetDD::SiDetectorElement*> (parameters->associatedSurface().associatedDetectorElement());
414 double phitol = 2.5; double etatol = 5.0;
415 bool nearGap = siElement->nearBondGap(parameters->localPosition(), etatol);
416 InDetDD::SiIntersect siIn = siElement->inDetector(parameters->localPosition(), phitol, etatol);
417 bool isInside = siIn.in();
418 m_sensitiveNearBondGap->push_back(nearGap ? 1 : 0);
419 m_sensitiveisInside->push_back(isInside ? 1 : 0);
420 auto& pBounds = parameters->associatedSurface().bounds();
421 bool isInsideBounds = pBounds.inside(parameters->localPosition(),phitol,etatol);
422 m_sensitiveisInsideBound->push_back(isInsideBounds ? 1 : 0);
424 m_sensitiveNearBondGap->push_back(-1);
425 m_sensitiveisInside->push_back(-1);
426 m_sensitiveisInsideBound->push_back(-1);
429 m_materialtInX0AccumulatedUpTo->push_back(m_materialThicknessInX0Accumulated->empty() ? 0. : m_materialThicknessInX0Accumulated->back());
432 // collect the material if configured to do so
433 if (m_collectMaterial && es.material) {
435 ATH_MSG_VERBOSE("Crossing material with scaling " << es.materialScaling);
436 ATH_MSG_VERBOSE(" material information " << *es.material);
437 // thickness in X0 and L0
438 float tInX0 = es.material->thicknessInX0() * es.materialScaling;
439 float tInL0 = es.material->thicknessInL0() * es.materialScaling;
440 float tzarho = es.material->zOverAtimesRhoTimesD() * es.materialScaling;
442 double path = es.materialScaling * es.material->thickness();
445 float mean = Trk::MaterialInteraction::dE_MPV_ionization(parameters->momentum().mag(),
446 (es.material->material()),
451 // sample the energy loss - m_tRandom.Landau(fabs(mpv),sigma);
452 // fabs(energyLoss)+energyLossSigma*CLHEP::RandLandau::shoot(m_randomEngine);
453 double eloss = mean - sigma * Trk::TrkExUnitTestBase::m_landauDist->shoot();
455 // the accummulated material
457 m_materialThicknessInX0 += tInX0;
458 m_materialThicknessInL0 += tInL0;
459 m_materialThicknessZARho += tzarho;
460 m_materialEmulatedIonizationLoss += eloss;
462 m_materialThicknessInX0Bwd += tInX0;
463 m_materialThicknessInL0Bwd += tInL0;
465 // the stepwise material
466 m_materialThicknessInX0Accumulated->push_back(m_materialThicknessInX0);
467 m_materialThicknessInX0Steps->push_back(tInX0);
468 m_materialThicknessInX0Steps->push_back(tInL0);
469 m_materialPositionX->push_back(es.materialPosition.x());
470 m_materialPositionY->push_back(es.materialPosition.y());
471 m_materialPositionZ->push_back(es.materialPosition.z());
472 m_materialPositionR->push_back(es.materialPosition.perp());
473 m_materialPositionP->push_back(parameters->momentum().mag());
474 m_materialPositionPt->push_back(parameters->momentum().perp());
475 m_materialScaling->push_back(es.materialScaling);
476 m_stepDirection->push_back(fwdBwd);
477 // check what type of material you have
478 if (es.stepConfiguration.checkMode(Trk::ExtrapolationMode::CollectSensitive)) {
479 m_materialThicknessInX0Sensitive += tInX0;
480 } else if (es.stepConfiguration.checkMode(Trk::ExtrapolationMode::CollectPassive)) {
481 m_materialThicknessInX0Passive += tInX0;
482 } else if (es.stepConfiguration.checkMode(Trk::ExtrapolationMode::CollectBoundary)) {
483 m_materialThicknessInX0Boundary += tInX0;
485 // check what type of surface you have
486 if (es.surface && es.surface->type() == Trk::SurfaceType::Cylinder) {
487 m_materialThicknessInX0Cylinder += tInX0;
488 } else if (es.surface && es.surface->type() == Trk::SurfaceType::Disc) {
489 m_materialThicknessInX0Disc += tInX0;
490 } else if (es.surface && es.surface->type() == Trk::SurfaceType::Plane) {
491 m_materialThicknessInX0Plane += tInX0;
494 // delete the parameters if there are any there
499 return StatusCode::SUCCESS;