10 #include "Acts/Definitions/Units.hpp"
11 #include "Acts/Propagator/ActorList.hpp"
12 #include "Acts/Propagator/Propagator.hpp"
13 #include "Acts/Propagator/MaterialInteractor.hpp"
14 #include "Acts/Propagator/EigenStepper.hpp"
15 #include "Acts/Propagator/AtlasStepper.hpp"
16 #include "Acts/Utilities/VectorHelpers.hpp"
17 #include "Acts/Visualization/ObjVisualization3D.hpp"
18 #include "Acts/Visualization/GeometryView3D.hpp"
19 #include "Acts/EventData/ParticleHypothesis.hpp"
20 #include "Acts/Propagator/detail/SteppingLogger.hpp"
33 using namespace Acts::UnitLiterals;
38 using SegLink_t = std::vector<ElementLink<xAOD::MuonSegmentContainer>>;
41 struct PropagatorRecorder{
57 double atlasPropabsMomentum{0.};
59 double actsPropabsMomentum{0.};
61 double actsStepSize{0.};
81 IdentifierHash MuonDetectorNavTest::layerHash(
const Identifier&
id)
const {
82 return m_r4DetMgr->getReadoutElement(
id)->layerHash(
id);
88 ATH_CHECK(m_fieldCacheCondObjInputKey.initialize());
89 ATH_CHECK(m_truthParticleKey.initialize());
90 ATH_CHECK(m_truthSegLinkKey.initialize());
94 ATH_CHECK(m_truthSegLinkKey.initialize());
99 return StatusCode::SUCCESS;
105 return StatusCode::SUCCESS;
109 const EventContext& ctx{Gaudi::Hive::currentContext()};
110 ATH_MSG_DEBUG(
"Execute in event "<<ctx.eventID().event_number());
122 const Acts::GeometryContext geoContext = gctx->context();
123 const Acts::MagneticFieldContext mfContext = Acts::MagneticFieldContext(fieldCondObj);
126 using Stepper = Acts::EigenStepper<>;
128 using Navigator = Acts::Experimental::DetectorNavigator;
129 using Propagator = Acts::Propagator<Stepper,Navigator>;
130 using ActorList = Acts::ActorList<Acts::detail::SteppingLogger, Acts::MaterialInteractor, Acts::EndOfWorldReached>;
131 using PropagatorOptions = Propagator::Options<ActorList>;
134 PropagatorOptions
options(geoContext, mfContext);
136 options.pathLimit = m_pathLimit;
137 options.stepping.maxStepSize = m_maxStepSize;
138 options.stepping.stepTolerance = m_stepTolerance;
140 options.maxTargetSkipping = m_maxTargetSkipping;
143 auto& materialInteractor =
options.actorList.get<Acts::MaterialInteractor>();
144 materialInteractor.energyLoss =
false;
145 materialInteractor.multipleScattering =
false;
146 materialInteractor.recordInteractions =
false;
150 Acts::Experimental::DetectorNavigator::Config navCfg;
152 navCfg.resolvePassive =
true;
154 auto bField = std::make_shared<ATLASMagneticFieldWrapper>();
155 auto bfieldCache = bField->makeCache(mfContext);
156 auto stepper =
Stepper(bField);
159 Acts::Experimental::DetectorNavigator navigator(navCfg, Acts::getDefaultLogger(
"DetectorNavigator",
ActsTrk::actsLevelVector(msgLevel())));
160 Acts::Propagator<Stepper, Acts::Experimental::DetectorNavigator> propagator(
161 std::move(stepper), std::move(navigator),
167 ATH_MSG_DEBUG(
"Consider truth particle "<<truthParticle->pt()<<
", "<<truthParticle->p4().P()<<
", "<<truthParticle->eta()<<
", "<<truthParticle->phi()
168 <<
", pdgId: "<<truthParticle->pdgId()<<
", status: "<<
HepMC::status(truthParticle)
172 std::vector<PropagatorRecorder> propagatedHits;
173 std::vector<std::pair<Identifier, Amg::Vector3D>> truthHits;
178 Acts::AnyCharge(std::abs(truthParticle->charge())));
181 std::vector<std::pair<const xAOD::MuonSegment*, std::vector<const xAOD::MuonSimHit*>>> muonSegmentWithSimHits;
182 for (
const auto& truthSegLink : segAcc(*truthParticle)){
185 std::vector<const xAOD::MuonSimHit*> muonSimHits{unordedHits.begin(), unordedHits.end()};
187 muonSegmentWithSimHits.emplace_back(seg,muonSimHits);
189 const auto&
particle = truthParticle->p4();
204 Amg::Vector3D magFieldPos{startPropPos.x()+
r, startPropPos.y()+
r, startPropPos.z()};
215 if(m_startFromFirstHit){
216 if (muonSegmentWithSimHits.empty()) {
220 std::ranges::sort(muonSegmentWithSimHits,
221 [](
const auto&
a,
const auto&
b) {
222 return a.first->position().
perp() <
b.first->position().perp();
226 std::vector<const xAOD::MuonSimHit*>& muonSimHits = muonSegmentWithSimHits.front().second;
228 for(
const auto& simHit: muonSimHits){
229 ATH_MSG_VERBOSE(
"Sim Hit of first segment at global position: "<<
Amg::toString(toGlobalTrf(*gctx, simHit->identify())*xAOD::toEigen(simHit->localPosition())));
234 std::ranges::sort(muonSimHits,
238 return globalPos1.norm()<globalPos2.norm();
242 for(
const auto& simHit: muonSimHits){
243 ATH_MSG_VERBOSE(
"Sim Hit of first segment at position: "<<
Amg::toString(toGlobalTrf(*gctx, simHit->identify())*xAOD::toEigen(simHit->localPosition())));
248 startPropPos = toGlobalTrf(*gctx, muonSimHits.front()->identify())*xAOD::toEigen(muonSimHits.front()->localPosition());
249 startPropDir = toGlobalTrf(*gctx, muonSimHits.front()->identify()).linear()*xAOD::toEigen(muonSimHits.front()->localDirection());
250 startPropP = muonSimHits.front()->kineticEnergy();
251 ATH_MSG_VERBOSE(
"Kinetic Energy from the simHit: "<<muonSimHits.front()->kineticEnergy() /
Gaudi::Units::GeV<<
" and mass: "<<muonSimHits.front()->mass()<<
" and energy deposit: "<<muonSimHits.front()->energyDeposit() /
Gaudi::Units::eV);
255 <<
Amg::toString(startPropDir)<<
" and momentum "<<startPropP);
258 Acts::BoundTrackParameters
start = Acts::BoundTrackParameters::createCurvilinear(
259 Acts::VectorHelpers::makeVector4(startPropPos, 0.), startPropDir,
260 truthParticle->charge() / startPropP,
269 const Acts::detail::SteppingLogger::result_type state =
result.value().get<Acts::detail::SteppingLogger::result_type>();
270 const Acts::MaterialInteractor::result_type material =
result.value().get<Acts::MaterialInteractor::result_type>();
271 ATH_MSG_DEBUG(
"Material interactions: " << material.materialInteractions.size());
274 std::vector<Acts::MaterialInteraction> interactions = material.materialInteractions;
275 for(
const Acts::MaterialInteraction& interaction : interactions){
276 ATH_MSG_DEBUG(
"Path Correction: " << interaction.pathCorrection);
281 Acts::ObjVisualization3D
helper{};
282 for (
const auto&
s : state.steps) {
283 if(
s.surface) Acts::GeometryView3D::drawSurface(
helper, *
s.surface, geoContext);
285 helper.write(
std::format(
"event_{:}_{:}_Hits.obj", ctx.eventID().event_number(), truthParticle->index()));
288 ATH_MSG_VERBOSE(
"Number of propagated steps : " << state.steps.size());
289 m_propSteps = state.steps.size();
290 m_propLength =
result.value().pathLength;
291 for(
const auto& state: state.steps){
302 ATH_MSG_VERBOSE(
"Identify propagated hit " << m_idHelperSvc->toString(
ID) <<
" with hit at local position " <<
Amg::toString(toGap*state.position)<<
" and global direction "<<
Amg::toString(state.momentum.unit()));
304 PropagatorRecorder newRecord{};
306 newRecord.actsPropPos = toGap*state.position;
307 newRecord.actsGlobalPos = state.position;
308 newRecord.actsPropDir = toGap.linear()*state.momentum.unit();
309 newRecord.actsPropabsMomentum = state.momentum.norm();
310 newRecord.actsStepSize = state.stepSize.value();
316 if(truthParticle->isMuon()){
321 std::unique_ptr<Trk::TrackParameters> atlasPars{m_extrapolator->extrapolate(ctx, atlasStart,
327 newRecord.atlasPropPos = toGap* atlasPars->position();
328 newRecord.atlasPropDir = toGap.linear()*atlasPars->momentum().unit();
329 newRecord.atlasPropabsMomentum = atlasPars->momentum().mag();
330 newRecord.atlasGlobalPos = atlasPars->position();
332 propagatedHits.emplace_back(newRecord);
335 std::size_t propHitsSize = propagatedHits.size();
338 unsigned int nMatchedTruth = 0, nMatchedProp = 0, nTruth{0};
339 for (
const auto&[
segment, simHits] : muonSegmentWithSimHits) {
345 toGlobalTrf(*gctx,simHit->identify())};
346 const Amg::Vector3D localPos = localTrf*xAOD::toEigen(simHit->localPosition());
347 const Amg::Vector3D globalPos = toGlobalTrf(*gctx, simHit->identify())*xAOD::toEigen(simHit->localPosition());
348 const Amg::Vector3D localDir = localTrf.linear()*xAOD::toEigen(simHit->localDirection());
349 m_detId.push_back(
ID);
350 m_techIdx.push_back(
toInt(m_idHelperSvc->technologyIndex(
ID)));
351 m_gasGapId.push_back(layerHash(
ID));
353 m_truthLoc.push_back(localPos);
354 m_truthDir.push_back(localDir);
355 m_truthGlob.push_back(globalPos);
356 m_startGlob.push_back(startPropPos);
359 auto it_begin = std::ranges::find_if(propagatedHits,
360 [
this,
ID](
const auto& propagatedHit) {
361 return m_idHelperSvc->detElId(
ID) == m_idHelperSvc->detElId(propagatedHit.id) &&
362 layerHash(
ID) == layerHash(propagatedHit.id);
364 m_isPropagated.push_back(it_begin != propagatedHits.end());
366 if (it_begin == propagatedHits.end()) {
367 m_actsPropLoc.push_back(dummyPos);
368 m_actsPropDir.push_back(dummyPos);
369 m_atlasPropLoc.push_back(dummyPos);
370 m_atlasPropDir.push_back(dummyPos);
371 m_actsPropGlob.push_back(dummyPos);
372 m_atlasPropGlob.push_back(dummyPos);
373 m_actsPropMomentum.push_back(0.);
374 m_atlasPropMomentum.push_back(0.);
375 m_actsStepSize.push_back(0.);
378 auto it_end = std::find_if(it_begin, propagatedHits.end(),
379 [
this,
ID](
const auto& propagatedHit) {
380 return m_idHelperSvc->detElId(ID) != m_idHelperSvc->detElId(propagatedHit.id) ||
381 layerHash(ID) != layerHash(propagatedHit.id);
384 auto it = std::min_element(it_begin, it_end,
385 [&localPos](
const PropagatorRecorder&
a,
386 const PropagatorRecorder&
b){
387 return (localPos -
a.actsPropPos).mag() < (localPos -
b.actsPropPos).mag();
391 ATH_MSG_DEBUG(
"Truth hit found in propagated hits: " << m_idHelperSvc->toString(
ID));
392 m_actsPropLoc.push_back(
it->actsPropPos);
393 m_actsPropDir.push_back(
it->actsPropDir);
394 m_actsPropGlob.push_back(
it->actsGlobalPos);
395 m_atlasPropLoc.push_back(
it->atlasPropPos);
396 m_atlasPropDir.push_back(
it->atlasPropDir);
397 m_atlasPropGlob.push_back(
it->atlasGlobalPos);
398 m_actsPropMomentum.push_back(
it->actsPropabsMomentum);
399 m_atlasPropMomentum.push_back(
it->atlasPropabsMomentum);
400 m_actsStepSize.push_back(
it->actsStepSize);
404 m_matchedTruthFraction = 1.f*nMatchedTruth / nTruth;
405 m_matchedPropFraction = 1.f*nMatchedProp / propHitsSize;
406 m_event = ctx.eventID().event_number();
410 return StatusCode::SUCCESS;