ATLAS Offline Software
Loading...
Searching...
No Matches
MuonTruthAlgs/src/TruthHitSummaryAlg.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
6
11
12namespace Muon {
13 using namespace MuonStationIndex;
14 // Initialize method:
16 ATH_CHECK(m_muonTruth.initialize());
17 ATH_CHECK(m_PRD_TruthNames.initialize());
18 ATH_CHECK(m_idHelperSvc.retrieve());
19
20 ATH_CHECK(m_nprecLayersKey.initialize());
21 ATH_CHECK(m_nphiLayersKey.initialize());
22 ATH_CHECK(m_ntrigEtaLayersKey.initialize());
23 ATH_CHECK(m_innerSmallHitsKey.initialize());
24 ATH_CHECK(m_innerLargeHitsKey.initialize());
25 ATH_CHECK(m_middleSmallHitsKey.initialize());
26 ATH_CHECK(m_middleLargeHitsKey.initialize());
27 ATH_CHECK(m_outerSmallHitsKey.initialize());
28 ATH_CHECK(m_outerLargeHitsKey.initialize());
31 ATH_CHECK(m_phiLayer1HitsKey.initialize());
32 ATH_CHECK(m_phiLayer2HitsKey.initialize());
33 ATH_CHECK(m_phiLayer3HitsKey.initialize());
34 ATH_CHECK(m_phiLayer4HitsKey.initialize());
35 ATH_CHECK(m_etaLayer1HitsKey.initialize());
36 ATH_CHECK(m_etaLayer2HitsKey.initialize());
37 ATH_CHECK(m_etaLayer3HitsKey.initialize());
38 ATH_CHECK(m_etaLayer4HitsKey.initialize());
39 ATH_CHECK(m_truthMdtHitsKey.initialize(m_idHelperSvc->hasMDT()));
40 ATH_CHECK(m_truthTgcHitsKey.initialize(m_idHelperSvc->hasTGC()));
41 ATH_CHECK(m_truthRpcHitsKey.initialize(m_idHelperSvc->hasRPC()));
42 ATH_CHECK(m_truthCscHitsKey.initialize(m_idHelperSvc->hasCSC()));
43 ATH_CHECK(m_truthStgcHitsKey.initialize(m_idHelperSvc->hasSTGC()));
44 ATH_CHECK(m_truthMMHitsKey.initialize(m_idHelperSvc->hasMM()));
45
46 return StatusCode::SUCCESS;
47 }
48
49 // Execute method:
50 StatusCode TruthHitSummaryAlg::execute(const EventContext& ctx) const {
51 // skip if no input data found
52 SG::ReadHandle muonTruthContainer(m_muonTruth, ctx);
53 ATH_CHECK(muonTruthContainer.isPresent());
54
55 summaryDecors myDecors{this, ctx};
56
57 // loop over truth coll (muon only)
58 for (const xAOD::TruthParticle* truthParticle : *muonTruthContainer) {
59 ChamberIdMap ids{};
60 ATH_CHECK(addHitCounts(ctx, *truthParticle, ids, myDecors));
61 ATH_CHECK(addHitIDVectors(*truthParticle, ids, myDecors));
62 }
63 return StatusCode::SUCCESS;
64 }
65
66 StatusCode TruthHitSummaryAlg::addHitCounts(const EventContext& ctx,
67 const xAOD::TruthParticle& truthParticle,
68 ChamberIdMap& ids, summaryDecors& myDecors) const {
69
70 std::vector<unsigned int> nprecHitsPerChamberLayer;
71 std::vector<unsigned int> nphiHitsPerChamberLayer;
72 std::vector<unsigned int> ntrigEtaHitsPerChamberLayer;
73
74 ntrigEtaHitsPerChamberLayer.resize(toInt(PhiIndex::PhiIndexMax));
75 nprecHitsPerChamberLayer.resize(toInt(ChIndex::ChIndexMax));
76 nphiHitsPerChamberLayer.resize(toInt(PhiIndex::PhiIndexMax));
77
78 ATH_MSG_DEBUG("addHitCounts: unique ID " << HepMC::uniqueID(truthParticle));
79 auto truthParticleHistory = HepMC::simulation_history(&truthParticle, -1); // Returns a list of unique IDs
80 // loop over detector technologies
82 ATH_CHECK(col.isPresent());
83
84 // loop over trajectories
85 for (const std::pair<Identifier, HepMcParticleLink> trajectory : *col) {
86 // check if gen particle same as input
87 if (std::ranges::find(truthParticleHistory, HepMC::uniqueID(trajectory.second)) == truthParticleHistory.end()) {
88 continue;
89 }
90 const Identifier& id = trajectory.first;
91 bool measPhi = m_idHelperSvc->measuresPhi(id);
92 bool isTgc = m_idHelperSvc->isTgc(id);
93 ChIndex chIndex = !isTgc ? m_idHelperSvc->chamberIndex(id) : ChIndex::ChUnknown;
94
95 // add identifier to map
96 if (m_idHelperSvc->isTgc(id)) { // TGCS should be added to both EIL and EIS
97 PhiIndex index = m_idHelperSvc->phiIndex(id);
98 if (index == PhiIndex::T4) {
99 ids[ChIndex::EIS].push_back(id);
100 ids[ChIndex::EIL].push_back(id);
101 } else {
102 ids[ChIndex::EMS].push_back(id);
103 ids[ChIndex::EML].push_back(id);
104 }
105 } else {
106 ids[m_idHelperSvc->chamberIndex(id)].push_back(id);
107 }
108 if (m_idHelperSvc->issTgc(id)) {
109 if (measPhi) {
110 PhiIndex index = m_idHelperSvc->phiIndex(id);
111 ++nphiHitsPerChamberLayer.at(toInt(index));
112 } else {
113 ++nprecHitsPerChamberLayer.at(toInt(chIndex));
114 }
115 } else if (m_idHelperSvc->isMM(id)) {
116 ++nprecHitsPerChamberLayer.at(toInt(chIndex));
117 } else if (m_idHelperSvc->isTrigger(id)) {
118 PhiIndex index = m_idHelperSvc->phiIndex(id);
119 if (index != PhiIndex::PhiUnknown) {
120 if (measPhi)
121 ++nphiHitsPerChamberLayer.at(toInt(index));
122 else
123 ++ntrigEtaHitsPerChamberLayer.at(toInt(index));
124 }
125 } else {
126 if (measPhi) {
127 PhiIndex index = m_idHelperSvc->phiIndex(id);
128 ++nphiHitsPerChamberLayer.at(toInt(index));
129 } else {
130 ++nprecHitsPerChamberLayer.at(toInt(chIndex));
131 }
132 }
133 }
134 }
135
136 uint8_t innerSmallHits = nprecHitsPerChamberLayer[toInt(ChIndex::BIS)] +
137 nprecHitsPerChamberLayer[toInt(ChIndex::EIS)] +
138 nprecHitsPerChamberLayer[toInt(ChIndex::CSS)];
139
140 uint8_t innerLargeHits = nprecHitsPerChamberLayer[toInt(ChIndex::BIL)] +
141 nprecHitsPerChamberLayer[toInt(ChIndex::EIL)] +
142 nprecHitsPerChamberLayer[toInt(ChIndex::CSL)];
143
144 uint8_t middleSmallHits = nprecHitsPerChamberLayer[toInt(ChIndex::BMS)] +
145 nprecHitsPerChamberLayer[toInt(ChIndex::EMS)];
146
147 uint8_t middleLargeHits = nprecHitsPerChamberLayer[toInt(ChIndex::BML)] +
148 nprecHitsPerChamberLayer[toInt(ChIndex::EML)];
149
150 uint8_t outerSmallHits = nprecHitsPerChamberLayer[toInt(ChIndex::BOS)] +
151 nprecHitsPerChamberLayer[toInt(ChIndex::EOS)];
152
153 uint8_t outerLargeHits = nprecHitsPerChamberLayer[toInt(ChIndex::BML)] +
154 nprecHitsPerChamberLayer[toInt(ChIndex::EOL)];
155
156 uint8_t extendedSmallHits = nprecHitsPerChamberLayer[toInt(ChIndex::EES)] +
157 nprecHitsPerChamberLayer[toInt(ChIndex::BEE)];
158
159 uint8_t extendedLargeHits = nprecHitsPerChamberLayer[toInt(ChIndex::EEL)];
160
161 uint8_t phiLayer1Hits = nphiHitsPerChamberLayer[toInt(PhiIndex::BM1)] +
162 nphiHitsPerChamberLayer[toInt(PhiIndex::T4)] +
163 nphiHitsPerChamberLayer[toInt(PhiIndex::CSC)] +
164 nphiHitsPerChamberLayer[toInt(PhiIndex::STGC1)] +
165 nphiHitsPerChamberLayer[toInt(PhiIndex::STGC2)];
166
167 uint8_t phiLayer2Hits = nphiHitsPerChamberLayer[toInt(PhiIndex::BM2)] +
168 nphiHitsPerChamberLayer[toInt(PhiIndex::T1)];
169
170 uint8_t phiLayer3Hits = nphiHitsPerChamberLayer[toInt(PhiIndex::BO1)] +
171 nphiHitsPerChamberLayer[toInt(PhiIndex::T2)];
172
173 uint8_t phiLayer4Hits = nphiHitsPerChamberLayer[toInt(PhiIndex::BO2)] +
174 nphiHitsPerChamberLayer[toInt(PhiIndex::T3)];
175
176 uint8_t etaLayer1Hits = ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BM1)] +
177 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T4)]+
178 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::CSC)] +
179 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::STGC1)] +
180 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::STGC2)];
181
182 uint8_t etaLayer2Hits = ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BM2)] +
183 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T1)];
184
185 uint8_t etaLayer3Hits = ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BO1)] +
186 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T2)];
187
188 uint8_t etaLayer4Hits = ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BO2)] +
189 ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T3)];
190
191 uint8_t nprecLayers = 0;
192 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::BIS)] + nprecHitsPerChamberLayer[toInt(ChIndex::BIL)] > 3);
193 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::BMS)] + nprecHitsPerChamberLayer[toInt(ChIndex::BML)] > 2);
194 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::BOS)] + nprecHitsPerChamberLayer[toInt(ChIndex::BOL)] > 2);
195 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::EIS)] + nprecHitsPerChamberLayer[toInt(ChIndex::EIL)] > 3);
196 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::EMS)] + nprecHitsPerChamberLayer[toInt(ChIndex::EML)] > 2);
197 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::EOS)] + nprecHitsPerChamberLayer[toInt(ChIndex::EOL)] > 2);
198 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::EES)] + nprecHitsPerChamberLayer[toInt(ChIndex::EEL)] > 3);
199 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::CSS)] + nprecHitsPerChamberLayer[toInt(ChIndex::CSL)] > 2);
200 nprecLayers += (nprecHitsPerChamberLayer[toInt(ChIndex::BEE)] > 3);
201
202 uint8_t nphiLayers = 0;
203 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::BM1)] > 0);
204 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::BM2)] > 0);
205 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::BO1)] > 0);
206 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::BO2)] > 0);
207 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::T1)] > 0);
208 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::T2)] > 0);
209 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::T3)] > 0);
210 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::T4)] > 0);
211 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::CSC)] > 2);
212 nphiLayers += (nphiHitsPerChamberLayer[toInt(PhiIndex::STGC1)] + nphiHitsPerChamberLayer[toInt(PhiIndex::STGC2)] > 3);
213
214 uint8_t ntrigEtaLayers = 0;
215 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BM1)] > 0);
216 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BM2)] > 0);
217 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BO1)] > 0);
218 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::BO2)] > 0);
219 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T1)] > 0);
220 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T2)] > 0);
221 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T3)] > 0);
222 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::T4)] > 0);
223 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::CSC)] > 2);
224 ntrigEtaLayers += (ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::STGC1)] + ntrigEtaHitsPerChamberLayer[toInt(PhiIndex::STGC2)] > 3);
225 // copy hit counts onto TruthParticle
226 (*myDecors.nprecLayersDecor)(truthParticle) = nprecLayers;
227 (*myDecors.nphiLayersDecor)(truthParticle) = nphiLayers;
228 (*myDecors.ntrigEtaLayersDecor)(truthParticle) = ntrigEtaLayers;
229 (*myDecors.innerSmallHitsDecor)(truthParticle) = innerSmallHits;
230 (*myDecors.innerLargeHitsDecor)(truthParticle) = innerLargeHits;
231 (*myDecors.middleSmallHitsDecor)(truthParticle) = middleSmallHits;
232 (*myDecors.middleLargeHitsDecor)(truthParticle) = middleLargeHits;
233 (*myDecors.outerSmallHitsDecor)(truthParticle) = outerSmallHits;
234 (*myDecors.outerLargeHitsDecor)(truthParticle) = outerLargeHits;
235 (*myDecors.extendedSmallHitsDecor)(truthParticle) = extendedSmallHits;
236 (*myDecors.extendedLargeHitsDecor)(truthParticle) = extendedLargeHits;
237
238 (*myDecors.phiLayer1HitsDecor)(truthParticle) = phiLayer1Hits;
239 (*myDecors.phiLayer2HitsDecor)(truthParticle) = phiLayer2Hits;
240 (*myDecors.phiLayer3HitsDecor)(truthParticle) = phiLayer3Hits;
241 (*myDecors.phiLayer4HitsDecor)(truthParticle) = phiLayer4Hits;
242
243 (*myDecors.etaLayer1HitsDecor)(truthParticle) = etaLayer1Hits;
244 (*myDecors.etaLayer2HitsDecor)(truthParticle) = etaLayer2Hits;
245 (*myDecors.etaLayer3HitsDecor)(truthParticle) = etaLayer3Hits;
246 (*myDecors.etaLayer4HitsDecor)(truthParticle) = etaLayer4Hits;
247
248
249 if (msgLvl(MSG::DEBUG)) {
250 ATH_MSG_DEBUG("Precision layers " << static_cast<int>(nprecLayers) << " phi layers " << static_cast<int>(nphiLayers)
251 << " triggerEta layers " << static_cast<int>(ntrigEtaLayers));
252
253 if (nprecLayers > 0) {
254 msg(MSG::VERBOSE) << " Precision chambers ";
255
256 for (int index = 0; index < static_cast<int>(nprecHitsPerChamberLayer.size()); ++index) {
257 if (nprecHitsPerChamberLayer[index] > 0)
258 msg(MSG::VERBOSE) << " " << chName(static_cast<ChIndex>(index))
259 << " hits " << nprecHitsPerChamberLayer[index];
260 }
261 }
262 if (nphiLayers > 0) {
263 msg(MSG::VERBOSE) << endmsg << " Phi chambers ";
264 for (int index = 0; index < static_cast<int>(nphiHitsPerChamberLayer.size()); ++index) {
265 if (nphiHitsPerChamberLayer[index] > 0)
266 msg(MSG::VERBOSE) << " " << phiName(static_cast<PhiIndex>(index))
267 << " hits " << nphiHitsPerChamberLayer[index];
268 }
269 }
270
271 if (ntrigEtaLayers > 0) {
272 msg(MSG::VERBOSE) << endmsg << " Trigger Eta ";
273 for (int index = 0; index < static_cast<int>(ntrigEtaHitsPerChamberLayer.size()); ++index) {
274 if (ntrigEtaHitsPerChamberLayer[index] > 0)
275 msg(MSG::VERBOSE) << " " << phiName(static_cast<PhiIndex>(index))
276 << " hits " << ntrigEtaHitsPerChamberLayer[index];
277 }
278 }
279 msg(MSG::VERBOSE) << endmsg;
280 }
281 return StatusCode::SUCCESS;
282 }
283
285 const ChamberIdMap& ids,
286 summaryDecors& myDecors) const {
287 std::vector<unsigned long long> mdtTruthHits{};
288 std::vector<unsigned long long> tgcTruthHits{};
289 std::vector<unsigned long long> rpcTruthHits{};
290 std::vector<unsigned long long> stgcTruthHits{};
291 std::vector<unsigned long long> cscTruthHits{};
292 std::vector<unsigned long long> mmTruthHits{};
293
294 // loop over chamber layers
295 int nEI = 0, nEM = 0;
296 for (const auto& lay : ids) {
297 // loop over hits
298 if (lay.first == ChIndex::EIS || lay.first == ChIndex::EIL) nEI++;
299 if (lay.first == ChIndex::EMS || lay.first == ChIndex::EML) nEM++;
300 for (const Identifier& id : lay.second) {
301 if (m_idHelperSvc->isMdt(id))
302 mdtTruthHits.push_back(id.get_compact());
303 else if (m_idHelperSvc->isCsc(id))
304 cscTruthHits.push_back(id.get_compact());
305 else if (m_idHelperSvc->isTgc(id)) {
306 if ((lay.first == ChIndex::EIS || lay.first == ChIndex::EIL) && nEI > 1)
307 continue; // otherwise we double-count
308 if ((lay.first == ChIndex::EMS || lay.first == ChIndex::EML) && nEM > 1)
309 continue; // otherwise we double-count
310 tgcTruthHits.push_back(id.get_compact());
311 } else if (m_idHelperSvc->issTgc(id))
312 stgcTruthHits.push_back(id.get_compact());
313 else if (m_idHelperSvc->isRpc(id))
314 rpcTruthHits.push_back(id.get_compact());
315 else if (m_idHelperSvc->isMM(id))
316 mmTruthHits.push_back(id.get_compact());
317 }
318 }
319 auto attatchHits = [&truthParticle](const WriteDecor_llvec& dec,
320 std::vector<unsigned long long>& hits) {
321 if (dec) {
322 (*dec)(truthParticle) = std::move(hits);
323 }
324 };
325 attatchHits(myDecors.truthMdtHitsDecor, mdtTruthHits);
326 attatchHits(myDecors.truthTgcHitsDecor, tgcTruthHits);
327 attatchHits(myDecors.truthRpcHitsDecor, rpcTruthHits);
328 attatchHits(myDecors.truthCscHitsDecor, cscTruthHits);
329 attatchHits(myDecors.truthStgcHitsDecor, stgcTruthHits);
330 attatchHits(myDecors.truthMMHitsDecor, mmTruthHits);
331 ATH_MSG_VERBOSE("Added " << mdtTruthHits.size() << " mdt truth hits, " << cscTruthHits.size() << " csc truth hits, "
332 << rpcTruthHits.size() << " rpc truth hits, and " << tgcTruthHits.size() << " tgc truth hits");
333 return StatusCode::SUCCESS;
334 }
335} // namespace Muon
#define endmsg
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_DEBUG(x)
bool msgLvl(const MSG::Level lvl) const
std::map< Muon::MuonStationIndex::ChIndex, std::vector< Identifier > > ChamberIdMap
This map contains all the hits corresponding to truth muons classified by chamber layer that recorded...
virtual StatusCode execute(const EventContext &ctx) const override
StatusCode addHitIDVectors(const xAOD::TruthParticle &truthParticle, const ChamberIdMap &ids, summaryDecors &myDecors) const
This function collapses the information given by the ChamberIdMap, which contains all hits from a tru...
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_muonTruth
Key of the container of truth muons.
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Handle for muonIdHelper service.
std::unique_ptr< SG::WriteDecorHandle< xAOD::TruthParticleContainer, llvec > > WriteDecor_llvec
SG::ReadHandleKeyArray< PRD_MultiTruthCollection > m_PRD_TruthNames
Keys for the containers of truth hits, grouped by detector technology.
StatusCode addHitCounts(const EventContext &ctx, const xAOD::TruthParticle &truthParticle, ChamberIdMap &ids, summaryDecors &myDecors) const
This is the actual function that decorates truth muons with hit counts.
bool isPresent() const
Is the referenced object present in SG?
int uniqueID(const T &p)
std::deque< int > simulation_history(const T &p, const int direction)
Function to calculate all the descendants(direction=1)/ancestors(direction=-1) of the particle.
ChIndex chIndex(const std::string &index)
convert ChIndex name string to enum
PhiIndex
enum to classify the different phi layers in the muon spectrometer
constexpr int toInt(const EnumType enumVal)
const std::string & chName(ChIndex index)
convert ChIndex into a string
const std::string & phiName(PhiIndex index)
convert PhiIndex into a string
ChIndex
enum to classify the different chamber layers in the muon spectrometer
NRpcCablingAlg reads raw condition data and writes derived condition data to the condition store.
Definition index.py:1
TruthParticle_v1 TruthParticle
Typedef to implementation.