ATLAS Offline Software
Loading...
Searching...
No Matches
ITkPixelRDOAnalysis.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
5
10
11#include "TTree.h"
12#include "TString.h"
13
14#include <algorithm>
15#include <math.h>
16#include <functional>
17#include <iostream>
18
19namespace ITk
20{
21
23 ATH_MSG_DEBUG( "Initializing ITkPixelRDOAnalysis" );
24
25 // This will check that the properties were initialized
26 // properly by job configuration.
27 ATH_CHECK( m_inputKey.initialize() );
28 ATH_CHECK( m_inputTruthKey.initialize() );
30
31 // Grab PixelID helper
32 ATH_CHECK(detStore()->retrieve(m_pixelID, m_pixelIDName.value()));
33 ATH_CHECK(detStore()->retrieve(m_pixelManager, m_detectorName.value()));
34
35
36 m_tree = new TTree(m_ntupleName.value().c_str(), "ITkPixelRDOAnalysis");
37 ATH_CHECK(histSvc()->regTree(m_ntuplePath.value() + m_ntupleName.value(), m_tree));
38 // PIXEL RDO
39 m_tree->Branch("rdoID", &m_rdoID);
40 m_tree->Branch("rdoWord", &m_rdoWord);
41 m_tree->Branch("barrelEndcap", &m_barrelEndcap);
42 m_tree->Branch("layerDisk", &m_layerDisk);
43 m_tree->Branch("phiModule", &m_phiModule);
44 m_tree->Branch("etaModule", &m_etaModule);
45 m_tree->Branch("phiIndex", &m_phiIndex);
46 m_tree->Branch("etaIndex", &m_etaIndex);
47 m_tree->Branch("isInnermost", &m_isInnermost);
48 m_tree->Branch("isNextToInnermost", &m_isNextToInnermost);
49 m_tree->Branch("ToT", &m_ToT); // time over threshold value (0-255)
50 m_tree->Branch("BCID", &m_BCID); // beam crossing ID
51 m_tree->Branch("LVL1A", &m_LVL1A); // Level1 accept (0-15)
52 m_tree->Branch("LVL1ID", &m_LVL1ID); // ATLAS LVL1 (0-255)
53 // Global coordinates
54 if (m_doPosition) {
55 m_tree->Branch("globalX", &m_globalX);
56 m_tree->Branch("globalY", &m_globalY);
57 m_tree->Branch("globalZ", &m_globalZ);
58 m_tree->Branch("localX", &m_localX);
59 m_tree->Branch("localY", &m_localY);
60 m_tree->Branch("localZ", &m_localZ);
61 }
62
63 // PIXEL SDO DEPOSITS
64 m_tree->Branch("sdoID", &m_sdoID);
65 m_tree->Branch("sdoWord", &m_sdoWord);
66 m_tree->Branch("barrelEndcap_sdo", &m_barrelEndcap_sdo);
67 m_tree->Branch("layerDisk_sdo", &m_layerDisk_sdo);
68 m_tree->Branch("phiModule_sdo", &m_phiModule_sdo);
69 m_tree->Branch("etaModule_sdo", &m_etaModule_sdo);
70 m_tree->Branch("phiIndex_sdo", &m_phiIndex_sdo);
71 m_tree->Branch("etaIndex_sdo", &m_etaIndex_sdo);
72 m_tree->Branch("noise", &m_noise);
73 m_tree->Branch("belowThresh", &m_belowThresh);
74 m_tree->Branch("disabled", &m_disabled);
75 m_tree->Branch("badTOT", &m_badTOT);
76 m_tree->Branch("barcode", &m_barcode);
77 m_tree->Branch("eventIndex", &m_eventIndex);
78 m_tree->Branch("charge", &m_charge);
79 m_tree->Branch("barcode_vec", &m_barcode_vec);
80 m_tree->Branch("eventIndex_vec", &m_eventIndex_vec);
81 m_tree->Branch("charge_vec", &m_charge_vec);
82
83 // HISTOGRAMS
84
86 m_h_rdoID = new TH1F("h_rdoID", "rdoID", 100, 0, 10e17);
87 m_h_rdoWord = new TH1F("h_rdoWord", "rdoWord", 100, 0, 350);
88 m_h_barrelEndcap = new TH1F("h_barrelEndcap", "Barrel or Endcap", 100, -5, 5);
89 m_h_layerDisk = new TH1F("h_layerDisk", "Barrel layer or Endcap disk", 100, 0, 10);
90 m_h_phiModule = new TH1F("h_phiModule", "Phi module", 100, 0, 100);
91 m_h_etaModule = new TH1F("h_etaModule", "Eta module", 100, -50, 50);
92 m_h_phiIndex = new TH1F("h_phiIndex", "Phi index", 820, 0, 820);
93 m_h_etaIndex = new TH1F("h_etaIndex", "Eta index", 820, 0, 820);
94 m_h_ToT = new TH1F("h_ToT", "ToT", 100, 0, 100);
95 m_h_BCID = new TH1F("h_BCID", "BCID", 100, -1.5, 1.5);
96 m_h_LVL1A = new TH1F("h_LVL1A", "LVL1A", 100, -1.5, 1.5);
97 m_h_LVL1ID = new TH1F("h_LVL1ID", "LVL1ID", 100, -1.5, 1.5);
98
100 m_h_brlLayer = new TH1F("h_brlLayer", "Barrel layer", 100, 0, 10);
101 m_h_brlPhiMod = new TH1F("h_brlPhiMod", "Barrel phi module", 100, 0, 100);
102 m_h_brlEtaMod = new TH1F("h_brlEtaMod", "Barrel eta module", 100, -50, 50);
103 m_h_brlPhiIndex = new TH1F("h_brlPhiIndex", "Barrel phi index", 820, 0, 820);
104 m_h_brlEtaIndex = new TH1F("h_brlEtaIndex", "Barrel eta index", 820, 0, 820);
105 m_h_brlToT = new TH1F("h_brlToT", "Barrel ToT", 100, 0, 100);
106 m_h_brlBCID = new TH1F("h_brlBCID", "Barrel BCID", 100, -1.5, 1.5);
107 m_h_brlLVL1A = new TH1F("h_brlLVL1A", "Barrel LVL1A", 100, -1.5, 1.5);
108 m_h_brlLVL1ID = new TH1F("h_brlLVL1ID", "Barrel LVL1ID", 100, -1.5, 1.5);
109
110 // ec histograms
111 m_h_ecDisk = new TH1F("h_ecDisk", "Endcap disk", 100, 0, 10);
112 m_h_ecPhiMod = new TH1F("h_ecPhiMod", "Endcap phi module", 100, 0, 100);
113 m_h_ecEtaMod = new TH1F("h_ecEtaMod", "Endcap eta module", 100, -50, 50);
114 m_h_ecPhiIndex = new TH1F("h_ecPhiIndex", "Endcap phi index", 820, 0, 820);
115 m_h_ecEtaIndex = new TH1F("h_ecEtaIndex", "Endcap eta index", 820, 0, 820);
116 m_h_ecToT = new TH1F("h_ecToT", "EndcapToT", 100, 0, 100);
117 m_h_ecBCID = new TH1F("h_ecBCID", "Endcap BCID", 100, -1.5, 1.5);
118 m_h_ecLVL1A = new TH1F("h_ecLVL1A", "Endcap LVL1A", 100, -1.5, 1.5);
119 m_h_ecLVL1ID = new TH1F("h_ecLVL1ID", "Endcap LVL1ID", 100, -1.5, 1.5);
120
121
123 m_h_sdoID = new TH1F("h_sdoID", "sdoID", 100, 0, 10e17);
124 m_h_sdoWord = new TH1F("h_sdoWord", "sdoWord", 100, 0, 350);
125 m_h_barrelEndcap_sdo = new TH1F("h_barrelEndcap_sdo", "Barrel or Endcap (SDO)", 100, -5, 5);
126 m_h_layerDisk_sdo = new TH1F("h_layerDisk_sdo", "Barrel layer or Endcap disk (SDO)", 100, 0, 10);
127 m_h_phiModule_sdo = new TH1F("h_phiModule_sdo", "Phi module (SDO)", 100, 0, 100);
128 m_h_etaModule_sdo = new TH1F("h_etaModule_sdo", "Eta module (SDO)", 100, -50, 50);
129 m_h_phiIndex_sdo = new TH1F("h_phiIndex_sdo", "Phi index (SDO)", 820, 0, 820);
130 m_h_etaIndex_sdo = new TH1F("h_etaIndex_sdo", "Eta index (SDO)", 820, 0, 820);
131 m_h_barcode = new TH1F("h_barcode", "Barcode (SDO)", 100, 0, 2.2e5);
132 m_h_eventIndex = new TH1F("h_eventIndex", "Event Index (SDO)", 100, 0, 2);
133 m_h_charge = new TH1F("h_charge", "Charge (SDO)", 100, 0, 1e7);
134
135
136 m_h_belowThresh_brl = new TH1F("h_belowThresh_brl", "Below threshold pixels - Barrel; # below threshold pixels; layer", 8, -0.5, 7.5);
137 m_h_disabled_brl = new TH1F("h_disabled_brl", "Disabled pixels - Barrel; # disabled pixels; layer", 8, -0.5, 7.5);
138
139 m_h_belowThresh_ec = new TH1F("h_belowThresh_ec", "Below threshold pixels - Endcap; # below threshold pixels; layer", 8, -0.5, 7.5);
140 m_h_disabled_ec = new TH1F("h_disabled_ec", "Disabled pixels - Endcap; # disabled pixels; layer", 8, -0.5, 7.5);
141
142
143 m_h_rdoID->StatOverflows();
144 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_rdoID->GetName(), m_h_rdoID));
145
146 m_h_rdoWord->StatOverflows();
147 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_rdoWord->GetName(), m_h_rdoWord));
148
149 m_h_barrelEndcap->StatOverflows();
151
152 m_h_layerDisk->StatOverflows();
153 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_layerDisk->GetName(), m_h_layerDisk));
154
155 m_h_phiModule->StatOverflows();
156 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_phiModule->GetName(), m_h_phiModule));
157
158 m_h_etaModule->StatOverflows();
159 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_etaModule->GetName(), m_h_etaModule));
160
161 m_h_phiIndex->StatOverflows();
162 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_phiIndex->GetName(), m_h_phiIndex));
163
164 m_h_etaIndex->StatOverflows();
165 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_etaIndex->GetName(), m_h_etaIndex));
166
167 m_h_ToT->StatOverflows();
168 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ToT->GetName(), m_h_ToT));
169
170 m_h_BCID->StatOverflows();
171 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_BCID->GetName(), m_h_BCID));
172
173 m_h_LVL1A->StatOverflows();
174 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_LVL1A->GetName(), m_h_LVL1A));
175
176 m_h_LVL1ID->StatOverflows();
177 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_LVL1ID->GetName(), m_h_LVL1ID));
178
179 m_h_brlLayer->StatOverflows();
180 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlLayer->GetName(), m_h_brlLayer));
181
182 m_h_brlPhiMod->StatOverflows();
183 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlPhiMod->GetName(), m_h_brlPhiMod));
184
185 m_h_brlEtaMod->StatOverflows();
186 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlEtaMod->GetName(), m_h_brlEtaMod));
187
188 m_h_brlPhiIndex->StatOverflows();
190
191 m_h_brlEtaIndex->StatOverflows();
193
194 m_h_brlToT->StatOverflows();
195 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlToT->GetName(), m_h_brlToT));
196
197 m_h_brlBCID->StatOverflows();
198 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlBCID->GetName(), m_h_brlBCID));
199
200 m_h_brlLVL1A->StatOverflows();
201 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlLVL1A->GetName(), m_h_brlLVL1A));
202
203 m_h_brlLVL1ID->StatOverflows();
204 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_brlLVL1ID->GetName(), m_h_brlLVL1ID));
205
206 m_h_ecDisk->StatOverflows();
207 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecDisk->GetName(), m_h_ecDisk));
208
209 m_h_ecPhiMod->StatOverflows();
210 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecPhiMod->GetName(), m_h_ecPhiMod));
211
212 m_h_ecEtaMod->StatOverflows();
213 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecEtaMod->GetName(), m_h_ecEtaMod));
214
215 m_h_ecPhiIndex->StatOverflows();
216 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecPhiIndex->GetName(), m_h_ecPhiIndex));
217
218 m_h_ecEtaIndex->StatOverflows();
219 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecEtaIndex->GetName(), m_h_ecEtaIndex));
220
221 m_h_ecToT->StatOverflows();
222 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecToT->GetName(), m_h_ecToT));
223
224 m_h_ecBCID->StatOverflows();
225 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecBCID->GetName(), m_h_ecBCID));
226
227 m_h_ecLVL1A->StatOverflows();
228 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecLVL1A->GetName(), m_h_ecLVL1A));
229
230 m_h_ecLVL1ID->StatOverflows();
231 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecLVL1ID->GetName(), m_h_ecLVL1ID));
232
233 m_h_sdoID->StatOverflows();
234 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_sdoID->GetName(), m_h_sdoID));
235
236 m_h_sdoWord->StatOverflows();
237 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_sdoWord->GetName(), m_h_sdoWord));
238
239 m_h_barrelEndcap_sdo->StatOverflows();
241
242 m_h_layerDisk_sdo->StatOverflows();
244
245 m_h_phiModule_sdo->StatOverflows();
247
248 m_h_etaModule_sdo->StatOverflows();
250
251 m_h_phiIndex_sdo->StatOverflows();
253
254 m_h_etaIndex_sdo->StatOverflows();
256
257 m_h_barcode->StatOverflows();
258 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_barcode->GetName(), m_h_barcode));
259
260 m_h_eventIndex->StatOverflows();
261 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_eventIndex->GetName(), m_h_eventIndex));
262
263 m_h_charge->StatOverflows();
264 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_charge->GetName(), m_h_charge));
265
266 m_h_belowThresh_brl->StatOverflows();
268
269 m_h_belowThresh_ec->StatOverflows();
271
272 m_h_disabled_brl->StatOverflows();
274
275 m_h_disabled_ec->StatOverflows();
277
278 for (unsigned int layer=0; layer<33; layer++) {
279 m_h_brlflatPhiIndex_perLayer.emplace_back(new TH1F(("h_brlflatPhiIndex_perLayer"+std::to_string(layer)).c_str(), ("Phi index - Barrel Flat - Layer "+std::to_string(layer)).c_str(), 820, 0, 820));
280 m_h_brlflatPhiIndex_perLayer.back()->StatOverflows();
282
283 m_h_brlflatEtaIndex_perLayer.emplace_back(new TH1F(("h_brlflatEtaIndex_perLayer"+std::to_string(layer)).c_str(), ("Eta index - Barrel Flat - Layer "+std::to_string(layer)).c_str(), 820, 0, 820));
284 m_h_brlflatEtaIndex_perLayer.back()->StatOverflows();
286
287 m_h_brlinclPhiIndex_perLayer.emplace_back(new TH1F(("h_brlinclPhiIndex_perLayer"+std::to_string(layer)).c_str(), ("Phi index - Barrel Inclined - Layer "+std::to_string(layer)).c_str(), 820, 0, 820));
288 m_h_brlinclPhiIndex_perLayer.back()->StatOverflows();
290
291 m_h_brlinclEtaIndex_perLayer.emplace_back(new TH1F(("h_brlinclEtaIndex_perLayer"+std::to_string(layer)).c_str(), ("Eta index - Barrel Inclined - Layer "+std::to_string(layer)).c_str(), 820, 0, 820));
292 m_h_brlinclEtaIndex_perLayer.back()->StatOverflows();
294
295 m_h_ecPhiIndex_perLayer.emplace_back(new TH1F(("h_ecPhiIndex_perLayer"+std::to_string(layer)).c_str(), ("Phi index - Endcap - Layer "+std::to_string(layer)).c_str(), 820, 0, 820));
296 m_h_ecPhiIndex_perLayer.back()->StatOverflows();
297 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecPhiIndex_perLayer.back()->GetName(), m_h_ecPhiIndex_perLayer.back()));
298
299 m_h_ecEtaIndex_perLayer.emplace_back(new TH1F(("h_ecEtaIndex_perLayer"+std::to_string(layer)).c_str(), ("Eta index - Endcap - Layer "+std::to_string(layer)).c_str(), 820, 0, 820));
300 m_h_ecEtaIndex_perLayer.back()->StatOverflows();
301 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_ecEtaIndex_perLayer.back()->GetName(), m_h_ecEtaIndex_perLayer.back()));
302 }
303
304 m_h_phiIndexInnermost = new TH1F("h_PhiIndexInnermost", "Phi index - Innermost Layer ", 820, 0, 820);
305 m_h_phiIndexInnermost->StatOverflows();
307
308 m_h_etaIndexInnermost = new TH1F("h_EtaIndexInnermost", "Eta index - Innermost Layer ", 820, 0, 820);
309 m_h_etaIndexInnermost->StatOverflows();
311
312 m_h_phiIndexNextToInnermost = new TH1F("h_PhiIndexNextToInnermost", "Phi index - Next To Innermost Layer ", 820, 0, 820);
313 m_h_phiIndexNextToInnermost->StatOverflows();
315
316 m_h_etaIndexNextToInnermost = new TH1F("h_EtaIndexNextToInnermost", "Eta index - Next To Innermost Layer ", 820, 0, 820);
317 m_h_etaIndexNextToInnermost->StatOverflows();
319
320 m_h_globalXY = new TH2F("h_globalXY","h_globalXY; x [mm]; y [mm]",700,-350.,350,700,-350.,350);
321 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_globalXY->GetName(), m_h_globalXY));
322 m_h_globalZR = new TH2F("h_globalZR","h_globalZR; z [mm]; r [mm]",6800,-3400.,3400,350,0.,350);
323 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_globalZR->GetName(), m_h_globalZR));
324 m_h_globalX = new TH1F("h_globalX","h_globalX; x [mm]",700,-350.,350.);
325 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_globalX->GetName(), m_h_globalX));
326 m_h_globalY = new TH1F("h_globalY","h_globalY; y [mm]",700,-350.,350.);
327 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_globalY->GetName(), m_h_globalY));
328 m_h_globalZ = new TH1F("h_globalZ","h_globalZ; z [mm]",6800,-3400.,3400.);
329 ATH_CHECK(histSvc()->regHist(m_histPath + m_h_globalZ->GetName(), m_h_globalZ));
330
331 // Special shared ITk histograms
332 std::string xy_name = "h_ITk_xy";
333 auto xy = std::make_unique<TH2D>(xy_name.c_str(), xy_name.c_str(), 2200, -1100, 1100, 2200, -1100, 1100);
334 xy->StatOverflows();
335 ATH_CHECK(histSvc()->regShared(m_sharedHistPath + xy_name, std::move(xy), m_h_globalXY_shared));
336
337 std::string zr_name = "h_ITk_zr";
338 auto zr = std::make_unique<TH2D>(zr_name.c_str(), zr_name.c_str(), 6800, -3400, 3400, 1100, 0, 1100);
339 zr->StatOverflows();
340 ATH_CHECK(histSvc()->regShared(m_sharedHistPath + zr_name, std::move(zr), m_h_globalZR_shared));
341
342 m_h_truthMatchedRDOs = new TH1F("h_TruthMatchedPixelRDOs", "h_TruthMatchedPixelRDOs", 4, 1, 5);
343 TString truthMatchBinLables[4] = { "All RDOs", "Truth Matched", "HS Matched", "Unmatched" };
344 for(unsigned int ibin = 1; ibin < 5; ibin++) {
345 m_h_truthMatchedRDOs->GetXaxis()->SetBinLabel(ibin, truthMatchBinLables[ibin-1]);
346 }
348 return StatusCode::SUCCESS;
349}
350
352 ATH_MSG_DEBUG(" In ITkPixelRDOAnalysis::execute()" );
353
354 m_rdoID->clear();
355 m_rdoWord->clear();
356 m_barrelEndcap->clear();
357 m_layerDisk->clear();
358 m_phiModule->clear();
359 m_etaModule->clear();
360 m_phiIndex->clear();
361 m_etaIndex->clear();
362 m_isInnermost->clear();
363 m_isNextToInnermost->clear();
364 m_ToT->clear();
365 m_BCID->clear();
366 m_LVL1A->clear();
367 m_LVL1ID->clear();
368 m_sdoID->clear();
369 m_sdoWord->clear();
370 m_barrelEndcap_sdo->clear();
371 m_layerDisk_sdo->clear();
372 m_phiModule_sdo->clear();
373 m_etaModule_sdo->clear();
374 m_phiIndex_sdo->clear();
375 m_etaIndex_sdo->clear();
376 m_noise->clear();
377 m_belowThresh->clear();
378 m_disabled->clear();
379 m_badTOT->clear();
380 m_barcode->clear();
381 m_eventIndex->clear();
382 m_charge->clear();
383 m_barcode_vec->clear();
384 m_eventIndex_vec->clear();
385 m_charge_vec->clear();
386 if (m_doPosition) {
387 m_globalX->clear();
388 m_globalY->clear();
389 m_globalZ->clear();
390 m_localX->clear();
391 m_localY->clear();
392 m_localZ->clear();
393 }
394 // Raw Data
395 const EventContext& ctx{Gaudi::Hive::currentContext()};
396
397 const PixelRDO_Container* p_pixelRDO_cont{nullptr};
398 ATH_CHECK(SG::get(p_pixelRDO_cont, m_inputKey, ctx));
399 //Adding SimMap and McEvent here for added truthMatching checks
400 const InDetSimDataCollection* simDataMapPixel{nullptr};
401 ATH_CHECK(SG::get(simDataMapPixel, m_inputTruthKey, ctx));
402 const McEventCollection* mcEventCollection{nullptr};
403 ATH_CHECK(SG::get(mcEventCollection, m_inputMcEventCollectionKey, ctx));
404 bool doTruthMatching = true;
405 const HepMC::GenEvent* hardScatterEvent(nullptr);
406
407 if (mcEventCollection->size()==0){
408 ATH_MSG_WARNING("Failed to retrieve a nonzero sized truth event collection, disabling truthMatching");
409 doTruthMatching = false;
410 }
411 if(doTruthMatching) hardScatterEvent = mcEventCollection->at(0);
412
413 if(p_pixelRDO_cont) {
414 // loop over RDO container
415 PixelRDO_Container::const_iterator rdoCont_itr(p_pixelRDO_cont->begin());
416 const PixelRDO_Container::const_iterator rdoCont_end(p_pixelRDO_cont->end());
417 for ( ; rdoCont_itr != rdoCont_end; ++rdoCont_itr ) {
418 const PixelRDO_Collection* p_pixelRDO_coll(*rdoCont_itr);
419 PixelRDO_Collection::const_iterator rdo_itr(p_pixelRDO_coll->begin());
420 const PixelRDO_Collection::const_iterator rdo_end(p_pixelRDO_coll->end());
421
422 const Identifier rdoIDColl(p_pixelRDO_coll->identify());
423 const int pixBrlEc(m_pixelID->barrel_ec(rdoIDColl));
424 const int pixLayerDisk(m_pixelID->layer_disk(rdoIDColl));
425 const int pixPhiMod(m_pixelID->phi_module(rdoIDColl));
426 const int pixEtaMod(m_pixelID->eta_module(rdoIDColl));
427 const int pixIsInnermost(0);
428 const int pixIsNextToInnermost(0);
429 //These are not yet implemented! Need to
430 //const int pixIsInnermost(m_pixelID->is_innermost(rdoIDColl));
431 //const int pixIsNextToInnermost(m_pixelID->is_nexttoinnermost(rdoIDColl));
432
433 const InDetDD::SiDetectorElement *detEl = m_pixelManager->getDetectorElement(rdoIDColl);
434
435 for ( ; rdo_itr != rdo_end; ++rdo_itr ) {
436 if(doTruthMatching){
437 m_h_truthMatchedRDOs->Fill(1.5);
438 bool findMatch = false;
439 if(simDataMapPixel){
440 InDetSimDataCollection::const_iterator iter = (*simDataMapPixel).find((*rdo_itr)->identify());
441
442 if ( iter != (*simDataMapPixel).end() ) {
443 const InDetSimData& sdo = iter->second;
444 const std::vector< InDetSimData::Deposit >& deposits = sdo.getdeposits();
445 std::vector< InDetSimData::Deposit >::const_iterator nextdeposit = deposits.begin();
446 std::vector< InDetSimData::Deposit >::const_iterator lastdeposit = deposits.end();
447 for( ; nextdeposit!=lastdeposit; ++nextdeposit) {
448 const HepMcParticleLink& particleLink = nextdeposit->first;
449 if(particleLink.isValid() && !findMatch){
450 HepMC::ConstGenParticlePtr genPart(particleLink.cptr());
451 if(genPart->parent_event() == hardScatterEvent) m_h_truthMatchedRDOs->Fill(3.5);
452 m_h_truthMatchedRDOs->Fill(2.5);
453 findMatch = true;
454 }
455 }
456 }
457 }
458 if(!findMatch) m_h_truthMatchedRDOs->Fill(4.5);
459 }
460 const Identifier rdoID((*rdo_itr)->identify());
461 const unsigned int rdoWord((*rdo_itr)->getWord());
462 const int pixPhiIx(m_pixelID->phi_index(rdoID));
463 const int pixEtaIx(m_pixelID->eta_index(rdoID));
464 const int pixToT((*rdo_itr)->getToT());
465 const int pixBCID((*rdo_itr)->getBCID());
466 const int pixLVL1A((*rdo_itr)->getLVL1A());
467 const int pixLVL1ID((*rdo_itr)->getLVL1ID());
468
469 //NB This is the position without lorentz correction
470 //(ATLSWUPGR-103)
471 //NB rationlization as here should probably also be done in master
472 //https://gitlab.cern.ch/atlas/athena/-/merge_requests/33398
473 InDetDD::SiLocalPosition localPos = detEl->rawLocalPositionOfCell(rdoID);
474 Amg::Vector3D globalPos = detEl->globalPosition(localPos);
475 if (m_doPosition) {
476 m_globalX->push_back(globalPos[Amg::x]);
477 m_globalY->push_back(globalPos[Amg::y]);
478 m_globalZ->push_back(globalPos[Amg::z]);
479
480 m_localX->push_back(localPos.xPhi());
481 m_localY->push_back(localPos.xEta());
482 m_localZ->push_back(localPos.xDepth());
483 }
484 float pixelRadius = sqrt(globalPos[Amg::x]*globalPos[Amg::x]+globalPos[Amg::y]*globalPos[Amg::y]);
485 m_h_globalXY->Fill(globalPos[Amg::x],globalPos[Amg::y]);
486 m_h_globalXY_shared->Fill(globalPos[Amg::x],globalPos[Amg::y]);
487 m_h_globalZR->Fill(globalPos[Amg::z],pixelRadius);
488 m_h_globalZR_shared->Fill(globalPos[Amg::z],pixelRadius);
489 m_h_globalX->Fill(globalPos[Amg::x]);
490 m_h_globalY->Fill(globalPos[Amg::y]);
491 m_h_globalZ->Fill(globalPos[Amg::z]);
492 const unsigned long long rdoID_int = rdoID.get_compact();
493 m_rdoID->push_back(rdoID_int);
494 m_rdoWord->push_back(rdoWord);
495 m_barrelEndcap->push_back(pixBrlEc);
496 m_layerDisk->push_back(pixLayerDisk);
497 m_phiModule->push_back(pixPhiMod);
498 m_etaModule->push_back(pixEtaMod);
499 m_isInnermost->push_back(pixIsInnermost);
500 m_isNextToInnermost->push_back(pixIsNextToInnermost);
501 m_phiIndex->push_back(pixPhiIx);
502 m_etaIndex->push_back(pixEtaIx);
503 m_ToT->push_back(pixToT);
504 m_BCID->push_back(pixBCID);
505 m_LVL1A->push_back(pixLVL1A);
506 m_LVL1ID->push_back(pixLVL1ID);
507
508 m_h_rdoID->Fill(rdoID_int);
509 m_h_rdoWord->Fill(rdoWord);
510 m_h_barrelEndcap->Fill(pixBrlEc);
511 m_h_layerDisk->Fill(pixLayerDisk);
512 m_h_phiModule->Fill(pixPhiMod);
513 m_h_etaModule->Fill(pixEtaMod);
514 m_h_phiIndex->Fill(pixPhiIx);
515 m_h_etaIndex->Fill(pixEtaIx);
516 m_h_ToT->Fill(pixToT);
517 m_h_BCID->Fill(pixBCID);
518 m_h_LVL1A->Fill(pixLVL1A);
519 m_h_LVL1ID->Fill(pixLVL1ID);
520
521 if (pixIsInnermost) {
522 m_h_phiIndexInnermost->Fill(pixPhiIx);
523 m_h_etaIndexInnermost->Fill(pixEtaIx);
524 } else if (pixIsNextToInnermost) {
525 m_h_phiIndexNextToInnermost->Fill(pixPhiIx);
526 m_h_etaIndexNextToInnermost->Fill(pixEtaIx);
527 }
528
529 //isInclined not yet implemented
530 //if (detEl->isBarrel() || detEl->isInclined()) {
531 if (detEl->isBarrel()){
532 m_h_brlLayer->Fill(pixLayerDisk);
533 m_h_brlPhiMod->Fill(pixPhiMod);
534 m_h_brlEtaMod->Fill(pixEtaMod);
535 m_h_brlPhiIndex->Fill(pixPhiIx);
536 m_h_brlEtaIndex->Fill(pixEtaIx);
537 m_h_brlToT->Fill(pixToT);
538 m_h_brlBCID->Fill(pixBCID);
539 m_h_brlLVL1A->Fill(pixLVL1A);
540 m_h_brlLVL1ID->Fill(pixLVL1ID);
541
542 //if (detEl->isInclined()) {
543 // m_h_brlinclPhiIndex_perLayer[pixLayerDisk]->Fill(pixPhiIx);
544 // m_h_brlinclEtaIndex_perLayer[pixLayerDisk]->Fill(pixEtaIx);
545 //} else {
546 m_h_brlflatPhiIndex_perLayer[pixLayerDisk]->Fill(pixPhiIx);
547 m_h_brlflatEtaIndex_perLayer[pixLayerDisk]->Fill(pixEtaIx);
548 //}
549
550 }
551 else if (detEl->isEndcap()) {
552 m_h_ecDisk->Fill(pixLayerDisk);
553 m_h_ecPhiMod->Fill(pixPhiMod);
554 m_h_ecEtaMod->Fill(pixEtaMod);
555 m_h_ecPhiIndex->Fill(pixPhiIx);
556 m_h_ecEtaIndex->Fill(pixEtaIx);
557 m_h_ecToT->Fill(pixToT);
558 m_h_ecBCID->Fill(pixBCID);
559 m_h_ecLVL1A->Fill(pixLVL1A);
560 m_h_ecLVL1ID->Fill(pixLVL1ID);
561
562 m_h_ecPhiIndex_perLayer[pixLayerDisk]->Fill(pixPhiIx);
563 m_h_ecEtaIndex_perLayer[pixLayerDisk]->Fill(pixEtaIx);
564
565 }
566 }
567 }
568 }
569
570 // Sim Data
571 if(simDataMapPixel) {
572 // loop over SDO container
573 InDetSimDataCollection::const_iterator sdo_itr(simDataMapPixel->begin());
574 const InDetSimDataCollection::const_iterator sdo_end(simDataMapPixel->end());
575
576 std::vector<int> barcode_vec;
577 std::vector<int> eventIndex_vec;
578 std::vector<float> charge_vec;
579 for ( ; sdo_itr != sdo_end; ++sdo_itr ) {
580 const Identifier sdoID((*sdo_itr).first);
581 const InDetSimData& sdo((*sdo_itr).second);
582 const unsigned long long sdoID_int = sdoID.get_compact();
583 const int sdoWord(sdo.word());
584 const int pixBrlEc_sdo(m_pixelID->barrel_ec(sdoID));
585 const int pixLayerDisk_sdo(m_pixelID->layer_disk(sdoID));
586 const int pixPhiMod_sdo(m_pixelID->phi_module(sdoID));
587 const int pixEtaMod_sdo(m_pixelID->eta_module(sdoID));
588 const int pixPhiIx_sdo(m_pixelID->phi_index(sdoID));
589 const int pixEtaIx_sdo(m_pixelID->eta_index(sdoID));
590 const bool noise(PixelSimHelper::isNoise(sdo));
591 const bool belowThresh(PixelSimHelper::isBelowThreshold(sdo));
592 const bool disabled(PixelSimHelper::isDisabled(sdo));
593 const bool badTOT(PixelSimHelper::hasBadTOT(sdo));
594
595 m_sdoID->push_back(sdoID_int);
596 m_sdoWord->push_back(sdoWord);
597 m_barrelEndcap_sdo->push_back(pixBrlEc_sdo);
598 m_layerDisk_sdo->push_back(pixLayerDisk_sdo);
599 m_phiModule_sdo->push_back(pixPhiMod_sdo);
600 m_etaModule_sdo->push_back(pixEtaMod_sdo);
601 m_phiIndex_sdo->push_back(pixPhiIx_sdo);
602 m_etaIndex_sdo->push_back(pixEtaIx_sdo);
603 m_noise->push_back(noise);
604
605 m_belowThresh->push_back(belowThresh);
606 if (belowThresh) {
607 if (pixBrlEc_sdo==0)
608 m_h_belowThresh_brl->Fill(pixLayerDisk_sdo);
609 else if (abs(pixBrlEc_sdo)==2)
610 m_h_belowThresh_ec->Fill(pixLayerDisk_sdo);
611 }
612
613 m_disabled->push_back(disabled);
614 if (disabled) {
615 if (pixBrlEc_sdo==0)
616 m_h_disabled_brl->Fill(pixLayerDisk_sdo);
617 else if (abs(pixBrlEc_sdo)==2)
618 m_h_disabled_ec->Fill(pixLayerDisk_sdo);
619 }
620
621 m_badTOT->push_back(badTOT);
622
623 m_h_sdoID->Fill(sdoID_int);
624 m_h_sdoWord->Fill(sdoWord);
625 m_h_barrelEndcap_sdo->Fill(pixBrlEc_sdo);
626 m_h_layerDisk_sdo->Fill(pixLayerDisk_sdo);
627 m_h_phiModule_sdo->Fill(pixPhiMod_sdo);
628 m_h_etaModule_sdo->Fill(pixEtaMod_sdo);
629 m_h_phiIndex_sdo->Fill(pixPhiIx_sdo);
630 m_h_etaIndex_sdo->Fill(pixEtaIx_sdo);
631
632 // loop over deposits
633 const std::vector<InDetSimData::Deposit>& deposits = sdo.getdeposits();
634 std::vector<InDetSimData::Deposit>::const_iterator dep_itr(deposits.begin());
635 const std::vector<InDetSimData::Deposit>::const_iterator dep_end(deposits.end());
636 for ( ; dep_itr != dep_end; ++dep_itr ) {
637 const HepMcParticleLink& particleLink = (*dep_itr).first;
638 const int bar(HepMC::barcode(particleLink)); // FIXME barcode-based
639 const int eventIx(particleLink.eventIndex());
640 const int charge((*dep_itr).second);
641
642 m_barcode->push_back(bar);
643 m_eventIndex->push_back(eventIx);
644 m_charge->push_back(charge);
645
646 m_h_barcode->Fill(bar);
647 m_h_eventIndex->Fill(eventIx);
648 m_h_charge->Fill(charge);
649
650 barcode_vec.push_back(bar);
651 eventIndex_vec.push_back(eventIx);
652 charge_vec.push_back(charge);
653 }
654 m_barcode_vec->push_back(barcode_vec);
655 m_eventIndex_vec->push_back(eventIndex_vec);
656 m_charge_vec->push_back(charge_vec);
657 barcode_vec.clear();
658 eventIndex_vec.clear();
659 charge_vec.clear();
660 }
661 }
662
663 m_tree->Fill();
664
665 return StatusCode::SUCCESS;
666}
667
668} // namespace ITk
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
double charge(const T &p)
Definition AtlasPID.h:997
InDetRawDataCollection< PixelRDORawData > PixelRDO_Collection
InDetRawDataContainer< InDetRawDataCollection< PixelRDORawData > > PixelRDO_Container
Handle class for reading from StoreGate.
const ServiceHandle< StoreGateSvc > & detStore() const
const ServiceHandle< ITHistSvc > & histSvc() const
The standard THistSvc (for writing histograms and TTrees and more to a root file) Returns (kind of) a...
DataModel_detail::const_iterator< DataVector > const_iterator
Definition DataVector.h:838
const T * at(size_type n) const
Access an element, as an rvalue.
const_iterator end() const noexcept
const_iterator begin() const noexcept
size_type size() const noexcept
Returns the number of elements in the collection.
SG::ReadHandleKey< InDetSimDataCollection > m_inputTruthKey
std::vector< TH1 * > m_h_ecEtaIndex_perLayer
std::vector< double > * m_globalZ
std::vector< unsigned long long > * m_sdoID
SG::ReadHandleKey< McEventCollection > m_inputMcEventCollectionKey
std::vector< int > * m_etaIndex
std::vector< int > * m_barrelEndcap_sdo
const InDetDD::PixelDetectorManager * m_pixelManager
std::vector< bool > * m_noise
std::vector< int > * m_phiIndex
std::vector< unsigned int > * m_rdoWord
std::vector< double > * m_localX
std::vector< double > * m_localZ
Gaudi::Property< std::string > m_histPath
std::vector< int > * m_LVL1A
std::vector< int > * m_charge
virtual StatusCode execute() override final
std::vector< TH1 * > m_h_ecPhiIndex_perLayer
std::vector< int > * m_eventIndex
std::vector< double > * m_globalX
std::vector< int > * m_barrelEndcap
std::vector< int > * m_barcode
std::vector< int > * m_phiModule_sdo
std::vector< std::vector< int > > * m_barcode_vec
std::vector< std::vector< float > > * m_charge_vec
std::vector< TH1 * > m_h_brlinclEtaIndex_perLayer
LockedHandle< TH2 > m_h_globalXY_shared
std::vector< unsigned long long > * m_rdoID
std::vector< int > * m_layerDisk_sdo
std::vector< int > * m_sdoWord
std::vector< int > * m_ToT
LockedHandle< TH2 > m_h_globalZR_shared
std::vector< double > * m_localY
std::vector< int > * m_BCID
std::vector< std::vector< int > > * m_eventIndex_vec
Gaudi::Property< std::string > m_ntuplePath
std::vector< double > * m_globalY
std::vector< int > * m_phiModule
virtual StatusCode initialize() override final
Gaudi::Property< std::string > m_sharedHistPath
Gaudi::Property< bool > m_doPosition
std::vector< int > * m_layerDisk
std::vector< bool > * m_belowThresh
std::vector< TH1 * > m_h_brlflatEtaIndex_perLayer
std::vector< int > * m_phiIndex_sdo
std::vector< TH1 * > m_h_brlinclPhiIndex_perLayer
std::vector< int > * m_etaIndex_sdo
std::vector< int > * m_LVL1ID
std::vector< int > * m_isInnermost
std::vector< int > * m_etaModule
std::vector< bool > * m_disabled
std::vector< bool > * m_badTOT
std::vector< int > * m_isNextToInnermost
Gaudi::Property< std::string > m_pixelIDName
Gaudi::Property< std::string > m_ntupleName
std::vector< TH1 * > m_h_brlflatPhiIndex_perLayer
Gaudi::Property< std::string > m_detectorName
SG::ReadHandleKey< PixelRDO_Container > m_inputKey
std::vector< int > * m_etaModule_sdo
const_iterator end() const
return const_iterator for end of container
const_iterator begin() const
return const_iterator for first entry
value_type get_compact() const
Get the compact id.
Class to hold geometrical description of a silicon detector element.
Class to represent a position in the natural frame of a silicon sensor, for Pixel and SCT For Pixel: ...
double xDepth() const
position along depth direction:
double xPhi() const
position along phi direction:
double xEta() const
position along eta direction:
HepGeom::Point3D< double > globalPosition(const HepGeom::Point3D< double > &localPos) const
transform a reconstruction local position into a global position (inline):
Amg::Vector2D rawLocalPositionOfCell(const SiCellId &cellId) const
Returns position (center) of cell.
virtual Identifier identify() const override final
int word() const
const std::vector< Deposit > & getdeposits() const
This defines the McEventCollection, which is really just an ObjectVector of McEvent objectsFile: Gene...
static bool isBelowThreshold(const InDetSimData &sdo)
static bool hasBadTOT(const InDetSimData &sdo)
static bool isNoise(const InDetSimData &sdo)
static bool isDisabled(const InDetSimData &sdo)
Eigen::Matrix< double, 3, 1 > Vector3D
int barcode(const T *p)
Definition Barcode.h:16
const GenParticle * ConstGenParticlePtr
Definition GenParticle.h:38
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.