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