ATLAS Offline Software
Loading...
Searching...
No Matches
MeasurementToTruthAssociationAlg.icc
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "MeasurementToTruthAssociationAlg.h"
6
7
8#include <iomanip>
9#include <cmath>
10#include <type_traits>
11#include <typeinfo>
12
13namespace ActsTrk {
14
15
16
17 template <class T_MeasurementCollection, class T_SimDataCollection, class T_TruthEventCollection, bool IsDebug>
18 StatusCode MeasurementToTruthAssociationAlg<T_MeasurementCollection,
19 T_SimDataCollection,
20 T_TruthEventCollection,
21 IsDebug>::initialize()
22 {
23 ATH_CHECK( m_measurementKey.initialize() );
24 ATH_CHECK( m_simDataKey.initialize() );
25 for (unsigned int i=0; i<m_statRDO.size(); ++i) {
26 m_statRDO[i]=0;
27 }
28 ATH_CHECK( m_associationOutKey.initialize() );
29 if constexpr(std::is_same<T_TruthEventCollection, void>::value) {
30 ATH_CHECK( m_truthEventCollectionKey.initialize(false) );
31 }
32 else {
33 ATH_MSG_INFO( " Truth event " << typeid(m_truthEventCollectionKey).name()
34 << " " << m_truthEventCollectionKey.key());
35 ATH_CHECK( m_truthEventCollectionKey.initialize() );
36 }
37 return StatusCode::SUCCESS;
38 }
39
40 template <class T_MeasurementCollection, class T_SimDataCollection, class T_TruthEventCollection, bool IsDebug>
41 StatusCode MeasurementToTruthAssociationAlg<T_MeasurementCollection,
42 T_SimDataCollection,
43 T_TruthEventCollection,
44 IsDebug>::finalize()
45 {
46 std::array<std::string, kNCategories> names {
47 "Measurements without SimData", // kNoTruth
48 "Measurements with SimData", // kHasSimHit
49 "Deposits without xAOD TruthParticle", // kInvalidTruthLink
50 "SimData without deposits above threshold", // kHasSimHitNoParticle
51 "Associated truth exceeds small vector size" // kBeyondSmallVectorSize
52 };
53 auto max_name_iter = std::max_element(names.begin(),names.end(),[](std::string &a,std::string &b) { return a.size() < b.size(); } );
54 for (unsigned int i=0; i<names.size(); ++i) {
55 ATH_MSG_INFO( "RDO truth stat " << std::left << std::setw(max_name_iter->size()) << names[i] << std::right << " " << m_statRDO[i]);
56 }
57 if constexpr(IsDebug) {
58 ATH_MSG_INFO("Truth particles per RDO " << m_stat.m_particlesPerMeasurement << std::endl
59 << m_stat.m_particlesPerMeasurement.histogramToString() );
60 ATH_MSG_INFO("Measurements per particle " << m_stat.m_measurementsPerParticle << std::endl
61 << m_stat.m_measurementsPerParticle.histogramToString() );
62 ATH_MSG_INFO("Log10 of deposited energy per RDO " << m_stat.m_depositedEnergy << std::endl
63 << m_stat.m_depositedEnergy.histogramToString() );
64 }
65 ATH_MSG_INFO("Deposits HS: " << m_depositCounts[1]
66 << " PU: " << m_depositCounts[0]
67 << " ; without truth particle HS: "
68 << " " << m_depositCounts[1+2]
69 << " PU: " << m_depositCounts[0+2] );
70
71 return StatusCode::SUCCESS;
72 }
73
74 template <class T_MeasurementCollection, class T_SimDataCollection, class T_TruthEventCollection, bool IsDebug>
75 StatusCode MeasurementToTruthAssociationAlg<T_MeasurementCollection,
76 T_SimDataCollection,
77 T_TruthEventCollection,
78 IsDebug>::execute(const EventContext &ctx) const
79 {
80 std::unique_ptr<MeasurementToTruthParticleAssociation>
81 association( std::make_unique<MeasurementToTruthParticleAssociation>() );
82
83 const T_SimDataCollection *simData{};
84 ATH_CHECK(SG::get(simData, m_simDataKey, ctx));
85
86 ATH_MSG_DEBUG("Retrieving measurement collection with key: " << m_measurementKey.key());
87 const T_MeasurementCollection *measurements{nullptr};
88 ATH_CHECK(SG::get(measurements, m_measurementKey, ctx));
89 ATH_MSG_DEBUG("Retrieved " << measurements->size() << " input measurements" );
90 association->resize( measurements->size() );
91 association->setSourceContainer(measurements, ctx);
92
93 const T_TruthEventCollection *truth_event_collection=nullptr;
94 if constexpr(!std::is_same_v<T_TruthEventCollection, void>) {
95 ATH_CHECK(SG::get(truth_event_collection, m_truthEventCollectionKey, ctx));
96 }
97
98 auto deposit_to_truth_map = makeDepositToTruthParticleMap(truth_event_collection);
99
100 typename std::conditional<IsDebug,Dbg::HistTemp,Dbg::Empty>::type stat(m_stat);
101 typename std::conditional<IsDebug,
102 std::unordered_map<const xAOD::TruthParticle *,unsigned int>,
103 Dbg::Empty>::type
104 hits_per_truthparticle;
105 std::array<std::size_t,4> depositCounts {0u, 0u, 0u, 0u};
106 std::array<unsigned int, kNCategories> truth_stat{};
107
108 std::vector<float> summed_contribution;
109 for ( const auto *measurement : *measurements ) {
110 summed_contribution.clear();
111 for (const auto& a_rdo : getRDOList(*measurement) ) {
112 typename T_SimDataCollection::const_iterator sim_data_iter(simData->find(Identifier(a_rdo)));
113
114 if(sim_data_iter != simData->end() ) {
115 auto deposits_for_measurement = getSimDataDeposits(*simData, sim_data_iter);
116 unsigned int n_particles=0;
117 for (const auto &deposit : deposits_for_measurement ) {
118 auto *truth_particle = deposit_to_truth_map.getTruthParticle(deposit);
119 ++(depositCounts.at(deposit_to_truth_map.isHardScatter(deposit)+(truth_particle ? 0u : 2u)));
120 if (!truth_particle) {
121 ++truth_stat[kInvalidTruthLink];
122 continue;
123 }
124 if constexpr(IsDebug) { stat.m_depositedEnergy.add( std::log10(getDepositedEnergy(deposit)) ); }
125
126 if (getDepositedEnergy(deposit) >= m_depositedEnergyMin.value()) {
127 if (measurement->index() >= association->size()) {
128 throw std::range_error("Measurement index out of range");
129 }
130 ParticleVector::iterator
131 particle_iter = std::find( (*association)[ measurement->index() ].begin(),
132 (*association)[ measurement->index() ].end(),
133 truth_particle );
134 if ( particle_iter == (*association)[ measurement->index() ].end()
135 || *particle_iter != truth_particle) {
136 particle_iter = (*association)[ measurement->index() ].insert( particle_iter, truth_particle);
137 ++n_particles;
138 }
139 unsigned int idx = (/*dest_iter*/ particle_iter - (*association)[ measurement->index() ].begin());
140 if (idx >= summed_contribution.size() ) {
141 summed_contribution.resize(idx+1, 0.);
142 }
143 summed_contribution[idx] += getDepositedEnergy(deposit);
144 if constexpr(IsDebug) {
145 std::unordered_map<const xAOD::TruthParticle *,unsigned int>::iterator
146 part_hits_iter = hits_per_truthparticle.find(truth_particle);
147 if (part_hits_iter == hits_per_truthparticle.end()) {
148 hits_per_truthparticle.insert(std::make_pair(truth_particle,1));
149 }
150 else {
151 ++(part_hits_iter->second);
152 }
153 }
154 }
155 }
156 if constexpr(IsDebug) { stat.m_particlesPerMeasurement.add(n_particles); }
157 ++truth_stat[ n_particles == 0 ? kHasSimHitNoParticle : kHasSimHit];
158 if (n_particles>ActsTrk::NTruthParticlesPerMeasurement) {
159 ++truth_stat[kBeyondSmallVectorSize];
160 }
161 }
162 else {
163 ++truth_stat[kNoTruth];
164 }
165 }
166 if ( (*association)[ measurement->index() ].size()>1) {
167 // move truth particle with largest contribution to the measurement to the beginning of the list
168 std::vector<float>::const_iterator iter_max=std::max_element(summed_contribution.begin(),summed_contribution.end());
169 unsigned int idx = iter_max - summed_contribution.begin();
170 if (idx>0) {
171 std::swap((*association)[ measurement->index() ][0],
172 (*association)[ measurement->index() ][idx]);
173 }
174 }
175 }
176 if constexpr(IsDebug) {
177 std::lock_guard<std::mutex> lock(m_stat.m_mutex);
178 m_stat.m_depositedEnergy += stat.m_depositedEnergy;
179 m_stat.m_particlesPerMeasurement += stat.m_particlesPerMeasurement;
180 for (const auto &particle_hit_counts : hits_per_truthparticle) {
181 m_stat.m_measurementsPerParticle.add( particle_hit_counts.second);
182 }
183 }
184 for (unsigned int i=0; i<m_depositCounts.size(); ++i) {
185 m_depositCounts[i] += depositCounts[i];
186 }
187 for (unsigned int i=0; i<m_statRDO.size(); ++i) {
188 m_statRDO[i]+=truth_stat[i];
189 }
190
191 SG::WriteHandle<MeasurementToTruthParticleAssociation> associationOutHandle(m_associationOutKey, ctx);
192 if (associationOutHandle.record( std::move(association)).isFailure()) {
193 ATH_MSG_ERROR("Failed to record measurement to truth assocition with key " << m_associationOutKey.key() );
194 return StatusCode::FAILURE;
195 }
196
197 return StatusCode::SUCCESS;
198 }
199
200
201
202
203} // namespace
Generated on for ATLAS Offline Software by doxygen 1.17.0