ATLAS Offline Software
Loading...
Searching...
No Matches
StripRDOAnalysis.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "StripRDOAnalysis.h"
6#include "SCT_ReadoutGeometry/SCT_DetectorManager.h"
7#include "InDetReadoutGeometry/SiDetectorElement.h"
8#include "ReadoutGeometryBase/SiLocalPosition.h"
9#include "StoreGate/ReadHandle.h"
10#include "TTree.h"
11#include "TString.h"
12
13#include <algorithm>
14#include <math.h>
15#include <functional>
16#include <iostream>
17
18namespace ITk
19{
20
21StripRDOAnalysis::StripRDOAnalysis(const std::string& name, ISvcLocator *pSvcLocator)
22 : AthAlgorithm(name, pSvcLocator)
23{
24}
25
26StatusCode StripRDOAnalysis::initialize() {
27 ATH_MSG_DEBUG( "Initializing SCT_RDOAnalysis" );
28
29 // This will check that the properties were initialized
30 // properly by job configuration.
31 ATH_CHECK( m_inputKey.initialize() );
32 ATH_CHECK( m_inputTruthKey.initialize() );
33 ATH_CHECK( m_inputMcEventCollectionKey.initialize() );
34
35 // Grab SCT_ID helper
36 ATH_CHECK(detStore()->retrieve(m_sctID, "SCT_ID"));
37
38 ATH_CHECK(detStore()->retrieve(m_SCT_Manager, "ITkStrip"));
39 // Grab Ntuple and histogramming service for tree
40 ATH_CHECK(m_thistSvc.retrieve());
41
42 m_tree = new TTree(m_ntupleName.value().c_str(), "ITkStripRDOAnalysis");
43 ATH_CHECK(m_thistSvc->regTree(m_ntuplePath.value() + m_ntupleName.value(), m_tree));
44 if (m_tree) {
45 // ITk Strip RDO
46 m_tree->Branch("rdoID", &m_rdoID);
47 m_tree->Branch("rdoWord", &m_rdoWord);
48 m_tree->Branch("barrelEndcap", &m_barrelEndcap);
49 m_tree->Branch("layerDisk", &m_layerDisk);
50 m_tree->Branch("phiModule", &m_phiModule);
51 m_tree->Branch("etaModule", &m_etaModule);
52 m_tree->Branch("side", &m_side);
53 m_tree->Branch("strip", &m_strip);
54 m_tree->Branch("row", &m_row);
55 m_tree->Branch("groupSize", &m_groupSize);
56 // Global coordinates
57 if (m_doPosition) {
58 m_tree->Branch("globalX0", &m_globalX0);
59 m_tree->Branch("globalY0", &m_globalY0);
60 m_tree->Branch("globalZ0", &m_globalZ0);
61 m_tree->Branch("globalX1", &m_globalX1);
62 m_tree->Branch("globalY1", &m_globalY1);
63 m_tree->Branch("globalZ1", &m_globalZ1);
64 m_tree->Branch("localX", &m_localX);
65 m_tree->Branch("localY", &m_localY);
66 m_tree->Branch("localZ", &m_localZ);
67 }
68 // ITk Strip SDO deposits
69 m_tree->Branch("sdoID", &m_sdoID);
70 m_tree->Branch("sdoWord", &m_sdoWord);
71 m_tree->Branch("barrelEndcap_sdo", &m_barrelEndcap_sdo);
72 m_tree->Branch("layerDisk_sdo", &m_layerDisk_sdo);
73 m_tree->Branch("phiModule_sdo", &m_phiModule_sdo);
74 m_tree->Branch("etaModule_sdo", &m_etaModule_sdo);
75 m_tree->Branch("side_sdo", &m_side_sdo);
76 m_tree->Branch("strip_sdo", &m_strip_sdo);
77 m_tree->Branch("row_sdo", &m_row_sdo);
78 m_tree->Branch("noise", &m_noise);
79 m_tree->Branch("belowThresh", &m_belowThresh);
80 m_tree->Branch("disabled", &m_disabled);
81 m_tree->Branch("barcode", &m_barcode);
82 m_tree->Branch("eventIndex", &m_eventIndex);
83 m_tree->Branch("charge", &m_charge);
84 m_tree->Branch("barcode_vec", &m_barcode_vec);
85 m_tree->Branch("eventIndex_vec", &m_eventIndex_vec);
86 m_tree->Branch("charge_vec", &m_charge_vec);
87 }
88 else {
89 ATH_MSG_ERROR("No tree found!");
90 }
91
92 // HISTOGRAMS
93 m_h_rdoID = new TH1F("h_rdoID", "rdoID", 100, 0, 25e17);
94 m_h_rdoID->StatOverflows();
95 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_rdoID->GetName(), m_h_rdoID));
96
97
98 m_h_rdoWord = new TH1F("h_rdoWord", "rdoWord", 100, 0, 17e6);
99 m_h_rdoWord->StatOverflows();
100 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_rdoWord->GetName(), m_h_rdoWord));
101
102 m_h_barrelEndcap = new TH1F("h_barrelEndcap", "Barrel or Endcap", 100, -3, 3);
103 m_h_barrelEndcap->StatOverflows();
104 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_barrelEndcap->GetName(), m_h_barrelEndcap));
105
106 m_h_layerDisk = new TH1F("h_layerDisk", "Barrel layer or Endcap disk", 100, 0, 10);
107 m_h_layerDisk->StatOverflows();
108 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_layerDisk->GetName(), m_h_layerDisk));
109
110 m_h_phiModule = new TH1F("h_phiModule", "Phi module", 100, 0, 60);
111 m_h_phiModule->StatOverflows();
112 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_phiModule->GetName(), m_h_phiModule));
113
114 m_h_etaModule = new TH1F("h_etaModule", "Eta module", 121, -60, 60);
115 m_h_etaModule->StatOverflows();
116 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_etaModule->GetName(), m_h_etaModule));
117
118 m_h_side = new TH1F("h_side", "Side", 100, 0, 1.5);
119 m_h_side->StatOverflows();
120 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_side->GetName(), m_h_side));
121
122 m_h_strip = new TH1F("h_strip", "Strip", 100, 0, 800);
123 m_h_strip->StatOverflows();
124 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_strip->GetName(), m_h_strip));
125
126 m_h_row = new TH1F("h_row", "Row", 100, 0, 4.5);
127 m_h_row->StatOverflows();
128 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_row->GetName(), m_h_row));
129
130 m_h_groupSize = new TH1F("h_groupSize", "Group size", 100, 0, 150);
131 m_h_groupSize->StatOverflows();
132 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_groupSize->GetName(), m_h_groupSize));
133
134 m_h_phi_v_eta = new TH2F("h_phi_v_eta", "Phi module vs eta module", 100, -7, 7, 100, 0, 60);
135 m_h_phi_v_eta->StatOverflows();
136 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_phi_v_eta->GetName(), m_h_phi_v_eta));
137
138 m_h_brlLayer = new TH1F("h_brlLayer", "Barrel layer", 100, 0, 10);
139 m_h_brlLayer->StatOverflows();
140 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brlLayer->GetName(), m_h_brlLayer));
141
142 m_h_brlPhiMod = new TH1F("h_brlPhiMod", "Barrel phi module", 100, 0, 60);
143 m_h_brlPhiMod->StatOverflows();
144 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brlPhiMod->GetName(), m_h_brlPhiMod));
145
146 m_h_brlEtaMod = new TH1F("h_brlEtaMod", "Barrel eta module", 121, -60, 60);
147 m_h_brlEtaMod->StatOverflows();
148 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brlEtaMod->GetName(), m_h_brlEtaMod));
149
150 m_h_brlSide = new TH1F("h_brlSide", "Barrel side", 100, 0, 1.5);
151 m_h_brlSide->StatOverflows();
152 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brlSide->GetName(), m_h_brlSide));
153
154 m_h_brlStrip = new TH1F("h_brlStrip", "Barrel strip", 100, 0, 800);
155 m_h_brlStrip->StatOverflows();
156 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brlStrip->GetName(), m_h_brlStrip));
157
158 m_h_brlGroupSize = new TH1F("h_brlGroupSize", "Barrel group size", 100, 0, 150);
159 m_h_brlGroupSize->StatOverflows();
160 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brlGroupSize->GetName(), m_h_brlGroupSize));
161
162 m_h_brl_phi_v_eta = new TH2F("h_brl_phi_v_eta", "Barrel phi module vs eta module", 100, -7, 7, 100, 0, 60);
163 m_h_brl_phi_v_eta->StatOverflows();
164 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brl_phi_v_eta->GetName(), m_h_brl_phi_v_eta));
165
166 m_h_ecDisk = new TH1F("h_ecDisk", "Endcap disk", 100, 0, 10);
167 m_h_ecDisk->StatOverflows();
168 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ecDisk->GetName(), m_h_ecDisk));
169
170 m_h_ecPhiMod = new TH1F("h_ecPhiMod", "Endcap phi module", 100, 0, 60);
171 m_h_ecPhiMod->StatOverflows();
172 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ecPhiMod->GetName(), m_h_ecPhiMod));
173
174 m_h_ecEtaMod = new TH1F("h_ecEtaMod", "Endcap eta module", 21, 0, 20);
175 m_h_ecEtaMod->StatOverflows();
176 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ecEtaMod->GetName(), m_h_ecEtaMod));
177
178 m_h_ecSide = new TH1F("h_ecSide", "Endcap side", 100, 0, 1.5);
179 m_h_ecSide->StatOverflows();
180 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ecSide->GetName(), m_h_ecSide));
181
182 m_h_ecStrip = new TH1F("h_ecStrip", "Endcap strip", 100, 0, 800);
183 m_h_ecStrip->StatOverflows();
184 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ecStrip->GetName(), m_h_ecStrip));
185
186 m_h_ecGroupSize = new TH1F("h_ecGroupSize", "Endcap group size", 100, 0, 150);
187 m_h_ecGroupSize->StatOverflows();
188 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ecGroupSize->GetName(), m_h_ecGroupSize));
189
190 m_h_ec_phi_v_eta = new TH2F("h_ec_phi_v_eta", "Endcap phi module vs eta module", 100, -7.5, 7.5, 100, 0, 60);
191 m_h_ec_phi_v_eta->StatOverflows();
192 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ec_phi_v_eta->GetName(), m_h_ec_phi_v_eta));
193
194 m_h_sdoID = new TH1F("h_sdoID", "sdoID", 100, 0, 1e18);
195 m_h_sdoID->StatOverflows();
196 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_sdoID->GetName(), m_h_sdoID));
197
198 m_h_sdoWord = new TH1F("h_sdoWord", "sdoWord", 100, 0, 1e7);
199 m_h_sdoWord->StatOverflows();
200 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_sdoWord->GetName(), m_h_sdoWord));
201
202 m_h_barrelEndcap_sdo = new TH1F("h_barrelEndcap_sdo", "Barrel or Endcap (SDO)", 100, -3, 3);
203 m_h_barrelEndcap_sdo->StatOverflows();
204 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_barrelEndcap_sdo->GetName(), m_h_barrelEndcap_sdo));
205
206 m_h_layerDisk_sdo = new TH1F("h_layerDisk_sdo", "Barrel layer or Endcap disk (SDO)", 100, 0, 10);
207 m_h_layerDisk_sdo->StatOverflows();
208 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_layerDisk_sdo->GetName(), m_h_layerDisk_sdo));
209
210 m_h_phiModule_sdo = new TH1F("h_phiModule_sdo", "Phi module (SDO)", 100, 0, 60);
211 m_h_phiModule_sdo->StatOverflows();
212 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_phiModule_sdo->GetName(), m_h_phiModule_sdo));
213
214 m_h_etaModule_sdo = new TH1F("h_etaModule_sdo", "Eta module (SDO)", 121, -60, 60);
215 m_h_etaModule_sdo->StatOverflows();
216 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_etaModule_sdo->GetName(), m_h_etaModule_sdo));
217
218 m_h_side_sdo = new TH1F("h_side_sdo", "Side (SDO)", 100, 0, 1.5);
219 m_h_side_sdo->StatOverflows();
220 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_side_sdo->GetName(), m_h_side_sdo));
221
222 m_h_strip_sdo = new TH1F("h_strip_sdo", "Strip (SDO)", 100, 0, 800);
223 m_h_strip_sdo->StatOverflows();
224 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_strip_sdo->GetName(), m_h_strip_sdo));
225
226 m_h_row_sdo = new TH1F("h_row_sdo", "Row (SDO)", 100, 0, 4.5);
227 m_h_row_sdo->StatOverflows();
228 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_row_sdo->GetName(), m_h_row_sdo));
229
230 m_h_barcode = new TH1F("h_barcode", "Barcode (SDO)", 100, 0, 2.2e5);
231 m_h_barcode->StatOverflows();
232 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_barcode->GetName(), m_h_barcode));
233
234 m_h_eventIndex = new TH1F("h_eventIndex", "Event index (SDO)", 100, 0, 10);
235 m_h_eventIndex->StatOverflows();
236 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_eventIndex->GetName(), m_h_eventIndex));
237
238 m_h_charge = new TH1F("h_charge", "Charge (SDO)", 100, 0, 6e6);
239 m_h_charge->StatOverflows();
240 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_charge->GetName(), m_h_charge));
241
242 m_h_phi_v_eta_sdo = new TH2F("h_phi_v_eta_sdo", "Phi module vs eta module (SDO)", 100, -7, 7, 100, 0, 60);
243 m_h_phi_v_eta_sdo->StatOverflows();
244 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_phi_v_eta_sdo->GetName(), m_h_phi_v_eta_sdo));
245
246 m_h_belowThresh_brl = new TH1F("h_belowThresh_brl", "Below threshold strips - Barrel; # below threshold strips; layer", 8, -0.5, 7.5);
247 m_h_belowThresh_brl->StatOverflows();
248 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_belowThresh_brl->GetName(), m_h_belowThresh_brl));
249
250 m_h_belowThresh_ec = new TH1F("h_belowThresh_ec", "Below threshold strips - Endcap; # below threshold strips; layer", 8, -0.5, 7.5);
251 m_h_belowThresh_ec->StatOverflows();
252 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_belowThresh_ec->GetName(), m_h_belowThresh_ec));
253
254 m_h_disabled_brl = new TH1F("h_disabled_brl", "Disabled strips - Barrel; # disabled strips; layer", 8, -0.5, 7.5);
255 m_h_disabled_brl->StatOverflows();
256 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_disabled_brl->GetName(), m_h_disabled_brl));
257
258 m_h_disabled_ec = new TH1F("h_disabled_ec", "Disabled strips - Endcap; # disabled strips; layer", 8, -0.5, 7.5);
259 m_h_disabled_ec->StatOverflows();
260 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_disabled_ec->GetName(), m_h_disabled_ec));
261
262 for (unsigned int layer=0; layer<4; layer++) {
263 m_h_brl_strip_perLayer.emplace_back(new TH1F(("h_brl_strip_perLayer"+std::to_string(layer)).c_str(), ("Strip index - Barrel - Layer "+std::to_string(layer)).c_str(), 1300, 0, 1300));
264 m_h_brl_strip_perLayer.back()->StatOverflows();
265 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_brl_strip_perLayer.back()->GetName(), m_h_brl_strip_perLayer.back()));
266 }
267
268 for (unsigned int layer=0; layer<9; layer++) {
269 m_h_ec_strip_perLayer.emplace_back(new TH1F(("h_ec_strip_perLayer"+std::to_string(layer)).c_str(), ("Strip index - Barrel - Layer "+std::to_string(layer)).c_str(), 1300, 0, 1300));
270 m_h_ec_strip_perLayer.back()->StatOverflows();
271 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_ec_strip_perLayer.back()->GetName(), m_h_ec_strip_perLayer.back()));
272 }
273
274 m_h_globalXY = new TH2F("h_globalXY","h_globalXY; x [mm]; y [mm]",2200,-1100.,1100.,2200,1100.,1100.);
275 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_globalXY->GetName(), m_h_globalXY));
276 m_h_globalZR = new TH2F("h_globalZR","h_globalZR; z [mm]; r [mm]",6800,-3400.,3400.,1100,0.,1100.);
277 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_globalZR->GetName(), m_h_globalZR));
278 m_h_globalX = new TH1F("h_globalX","h_globalX; x [mm]",2200,-1100.,1100.);
279 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_globalX->GetName(), m_h_globalX));
280 m_h_globalY = new TH1F("h_globalY","h_globalY; y [mm]",2200,-1100.,1100.);
281 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_globalY->GetName(), m_h_globalY));
282 m_h_globalZ = new TH1F("h_globalZ","h_globalZ; z [mm]",6800,-3400.,3400.);
283 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_globalZ->GetName(), m_h_globalZ));
284
285 // Special shared ITk histograms
286 std::string xy_name = "h_ITk_xy";
287 auto xy = std::make_unique<TH2D>(xy_name.c_str(), xy_name.c_str(), 2200, -1100, 1100, 2200, -1100, 1100);
288 xy->StatOverflows();
289 ATH_CHECK(m_thistSvc->regShared(m_sharedHistPath + xy_name, std::move(xy), m_h_globalXY_shared));
290
291 std::string zr_name = "h_ITk_zr";
292 auto zr = std::make_unique<TH2D>(zr_name.c_str(), zr_name.c_str(), 6800, -3400, 3400, 1100, 0, 1100);
293 zr->StatOverflows();
294 ATH_CHECK(m_thistSvc->regShared(m_sharedHistPath + zr_name, std::move(zr), m_h_globalZR_shared));
295
296 m_h_truthMatchedRDOs = new TH1F("h_TruthMatchedITkStripRDOs", "h_TruthMatchedITkStripRDOs", 4, 1, 5);
297 TString truthMatchBinLables[4] = { "All RDOs", "Truth Matched", "HS Matched", "Unmatched" };
298 for(unsigned int ibin = 1; ibin < 5; ibin++) {
299 m_h_truthMatchedRDOs->GetXaxis()->SetBinLabel(ibin, truthMatchBinLables[ibin-1]);
300 }
301 ATH_CHECK(m_thistSvc->regHist(m_histPath + m_h_truthMatchedRDOs->GetName(), m_h_truthMatchedRDOs));
302
303 return StatusCode::SUCCESS;
304}
305
306StatusCode StripRDOAnalysis::execute() {
307 ATH_MSG_DEBUG( "In StripRDOAnalysis::execute()" );
308
309 m_rdoID->clear();
310 m_rdoWord->clear();
311 m_barrelEndcap->clear();
312 m_layerDisk->clear();
313 m_phiModule->clear();
314 m_etaModule->clear();
315 m_side->clear();
316 m_strip->clear();
317 m_row->clear();
318 m_groupSize->clear();
319 if (m_doPosition) {
320 m_globalX0->clear();
321 m_globalY0->clear();
322 m_globalZ0->clear();
323 m_globalX1->clear();
324 m_globalY1->clear();
325 m_globalZ1->clear();
326 m_localX->clear();
327 m_localY->clear();
328 m_localZ->clear();
329 }
330 m_sdoID->clear();
331 m_sdoWord->clear();
332 m_barrelEndcap_sdo->clear();
333 m_layerDisk_sdo->clear();
334 m_phiModule_sdo->clear();
335 m_etaModule_sdo->clear();
336 m_side_sdo->clear();
337 m_strip_sdo->clear();
338 m_row_sdo->clear();
339 m_noise->clear();
340 m_belowThresh->clear();
341 m_disabled->clear();
342 m_barcode->clear();
343 m_eventIndex->clear();
344 m_charge->clear();
345 m_barcode_vec->clear();
346 m_eventIndex_vec->clear();
347 m_charge_vec->clear();
348
349 // RawData
350 SG::ReadHandle<SCT_RDO_Container> p_SCT_RDO_cont (m_inputKey);
351 //Adding SimMap and McEvent here for added truthMatching checks
352 SG::ReadHandle<InDetSimDataCollection> simDataMapSCT (m_inputTruthKey);
353 SG::ReadHandle<McEventCollection> mcEventCollection (m_inputMcEventCollectionKey);
354
355 const HepMC::GenEvent* hardScatterEvent(nullptr);
356 bool doTruthMatching = true;
357 if (mcEventCollection->size()==0){
358 ATH_MSG_WARNING("Failed to retrieve a nonzero sized truth event collection, disabling truthMatching");
359 doTruthMatching = false;
360 }
361 if(doTruthMatching) hardScatterEvent = mcEventCollection->at(0);
362
363 if(p_SCT_RDO_cont.isValid()) {
364 // loop over RDO container
365 SCT_RDO_Container::const_iterator rdoCont_itr(p_SCT_RDO_cont->begin());
366 const SCT_RDO_Container::const_iterator rdoCont_end(p_SCT_RDO_cont->end());
367
368 for ( ; rdoCont_itr != rdoCont_end; ++rdoCont_itr ) {
369 const SCT_RDO_Collection* p_SCT_RDO_coll(*rdoCont_itr);
370 SCT_RDO_Collection::const_iterator rdo_itr(p_SCT_RDO_coll->begin());
371 const SCT_RDO_Collection::const_iterator rdo_end(p_SCT_RDO_coll->end());
372
373 for ( ; rdo_itr != rdo_end; ++rdo_itr ) {
374 if(doTruthMatching){
375 m_h_truthMatchedRDOs->Fill(1.5);
376 bool findMatch = false;
377 if(simDataMapSCT.isValid()){
378 InDetSimDataCollection::const_iterator iter = (*simDataMapSCT).find((*rdo_itr)->identify());
379
380 if ( iter != (*simDataMapSCT).end() ) {
381 const InDetSimData& sdo = iter->second;
382 const std::vector< InDetSimData::Deposit >& deposits = sdo.getdeposits();
383 std::vector< InDetSimData::Deposit >::const_iterator nextdeposit = deposits.begin();
384 std::vector< InDetSimData::Deposit >::const_iterator lastdeposit = deposits.end();
385 for( ; nextdeposit!=lastdeposit; ++nextdeposit) {
386 const HepMcParticleLink& particleLink = nextdeposit->first;
387 if(particleLink.isValid() && !findMatch){
388 HepMC::ConstGenParticlePtr genPart(particleLink.cptr());
389 if(genPart->parent_event() == hardScatterEvent) m_h_truthMatchedRDOs->Fill(3.5);
390 m_h_truthMatchedRDOs->Fill(2.5);
391 findMatch = true;
392 }
393 }
394 }
395 }
396 if(!findMatch) m_h_truthMatchedRDOs->Fill(4.5);
397 }
398 const Identifier rdoID((*rdo_itr)->identify());
399 const unsigned int rdoWord((*rdo_itr)->getWord());
400 const int sctBrlEc(m_sctID->barrel_ec(rdoID));
401 const int sctLayerDisk(m_sctID->layer_disk(rdoID));
402 const int sctPhiMod(m_sctID->phi_module(rdoID));
403 const int sctEtaMod(m_sctID->eta_module(rdoID));
404 const int sctSide(m_sctID->side(rdoID));
405 const int sctStrip(m_sctID->strip(rdoID));
406 const int sctRow(m_sctID->row(rdoID));
407 const int sctGroupSize((*rdo_itr)->getGroupSize());
408
409 const unsigned long long rdoID_int = rdoID.get_compact();
410
411
412 const InDetDD::SiDetectorElement *detEl = m_SCT_Manager->getDetectorElement(rdoID);
413
414 if(!detEl) {
415 ATH_MSG_WARNING("No Element found for ID "<<m_sctID->show_to_string(rdoID)<<" - skipping!");
416 continue;
417 }
418
419 //Using lorentz-angle corrected version here will result in inconsistencies (ATLSWUPGR-103)
420 //NB rationlization as here should probably also be done in master
421 //https://gitlab.cern.ch/atlas/athena/-/merge_requests/33398
422
423 Amg::Vector2D localPos = detEl->rawLocalPositionOfCell(rdoID);
424
425 std::pair<Amg::Vector3D, Amg::Vector3D> endsOfStrip = detEl->endsOfStrip(localPos);
426
427 if (m_doPosition) {
428
429 m_globalX0->push_back(endsOfStrip.first.x());
430 m_globalY0->push_back(endsOfStrip.first.y());
431 m_globalZ0->push_back(endsOfStrip.first.z());
432
433 m_globalX1->push_back(endsOfStrip.second.x());
434 m_globalY1->push_back(endsOfStrip.second.y());
435 m_globalZ1->push_back(endsOfStrip.second.z());
436
437 m_localX->push_back(localPos.x());
438 m_localY->push_back(localPos.y());
439 m_localZ->push_back(0.0);
440
441 }
442 float stripradius0 = sqrt(endsOfStrip.first.x()*endsOfStrip.first.x()+endsOfStrip.first.y()*endsOfStrip.first.y());
443 float stripradius1 = sqrt(endsOfStrip.second.x()*endsOfStrip.second.x()+endsOfStrip.second.y()*endsOfStrip.second.y());
444
445 m_h_globalXY->Fill(endsOfStrip.first.x(),endsOfStrip.first.y());
446 m_h_globalXY->Fill(endsOfStrip.second.x(),endsOfStrip.second.y());
447 m_h_globalXY_shared->Fill(endsOfStrip.first.x(),endsOfStrip.first.y());
448 m_h_globalXY_shared->Fill(endsOfStrip.second.x(),endsOfStrip.second.y());
449 m_h_globalZR->Fill(endsOfStrip.first.z(),stripradius0);
450 m_h_globalZR->Fill(endsOfStrip.second.z(),stripradius1);
451 m_h_globalZR_shared->Fill(endsOfStrip.first.z(),stripradius0);
452 m_h_globalZR_shared->Fill(endsOfStrip.second.z(),stripradius1);
453 m_h_globalX->Fill(endsOfStrip.first.x());
454 m_h_globalY->Fill(endsOfStrip.first.y());
455 m_h_globalZ->Fill(endsOfStrip.first.z());
456 m_h_globalX->Fill(endsOfStrip.second.x());
457 m_h_globalY->Fill(endsOfStrip.second.y());
458 m_h_globalZ->Fill(endsOfStrip.second.z());
459
460 m_rdoID->push_back(rdoID_int);
461 m_rdoWord->push_back(rdoWord);
462 m_barrelEndcap->push_back(sctBrlEc);
463 m_layerDisk->push_back(sctLayerDisk);
464 m_phiModule->push_back(sctPhiMod);
465 m_etaModule->push_back(sctEtaMod);
466 m_side->push_back(sctSide);
467 m_row->push_back(sctRow);
468 m_strip->push_back(sctStrip);
469 m_groupSize->push_back(sctGroupSize);
470
471 m_h_rdoID->Fill(rdoID_int);
472 m_h_rdoWord->Fill(rdoWord);
473 m_h_barrelEndcap->Fill(sctBrlEc);
474 m_h_layerDisk->Fill(sctLayerDisk);
475 m_h_phiModule->Fill(sctPhiMod);
476 m_h_etaModule->Fill(sctEtaMod);
477 m_h_side->Fill(sctSide);
478 m_h_strip->Fill(sctStrip);
479 m_h_row->Fill(sctRow);
480 m_h_groupSize->Fill(sctGroupSize);
481 m_h_phi_v_eta->Fill(sctEtaMod, sctPhiMod);
482
483 if (sctBrlEc == 0) {
484 m_h_brlLayer->Fill(sctLayerDisk);
485 m_h_brlPhiMod->Fill(sctPhiMod);
486 m_h_brlEtaMod->Fill(sctEtaMod);
487 m_h_brlSide->Fill(sctSide);
488 m_h_brlStrip->Fill(sctStrip);
489 m_h_brlGroupSize->Fill(sctGroupSize);
490 m_h_brl_phi_v_eta->Fill(sctEtaMod, sctPhiMod);
491 m_h_brl_strip_perLayer[sctLayerDisk]->Fill(sctStrip);
492 m_h_brl_strip_perLayer[sctLayerDisk]->Fill(sctStrip);
493 }
494 else if (abs(sctBrlEc) == 2) {
495 m_h_ecDisk->Fill(sctLayerDisk);
496 m_h_ecPhiMod->Fill(sctPhiMod);
497 m_h_ecEtaMod->Fill(sctEtaMod);
498 m_h_ecSide->Fill(sctSide);
499 m_h_ecStrip->Fill(sctStrip);
500 m_h_ecGroupSize->Fill(sctGroupSize);
501 m_h_ec_phi_v_eta->Fill(sctEtaMod, sctPhiMod);
502 m_h_ec_strip_perLayer[sctLayerDisk]->Fill(sctStrip);
503 m_h_ec_strip_perLayer[sctLayerDisk]->Fill(sctStrip);
504 }
505 }
506 }
507 }
508
509 // SimData
510 if(simDataMapSCT.isValid()) {
511 // loop over SDO container
512 InDetSimDataCollection::const_iterator sdo_itr(simDataMapSCT->begin());
513 const InDetSimDataCollection::const_iterator sdo_end(simDataMapSCT->end());
514
515 std::vector<int> barcode_vec;
516 std::vector<int> eventIndex_vec;
517 std::vector<float> charge_vec;
518 for ( ; sdo_itr != sdo_end; ++sdo_itr ) {
519 const Identifier sdoID((*sdo_itr).first);
520 const InDetSimData& sdo((*sdo_itr).second);
521 const unsigned long long sdoID_int = sdoID.get_compact();
522 const int sdoWord(sdo.word());
523 const int sctBrlEc_sdo(m_sctID->barrel_ec(sdoID));
524 const int sctLayerDisk_sdo(m_sctID->layer_disk(sdoID));
525 const int sctPhiMod_sdo(m_sctID->phi_module(sdoID));
526 const int sctEtaMod_sdo(m_sctID->eta_module(sdoID));
527 const int sctSide_sdo(m_sctID->side(sdoID));
528 const int sctStrip_sdo(m_sctID->strip(sdoID));
529 const int sctRow_sdo(m_sctID->row(sdoID));
530 const bool noise(SCT_SimHelper::isNoise(sdo));
531 const bool belowThresh(SCT_SimHelper::isBelowThreshold(sdo));
532 const bool disabled(SCT_SimHelper::isDisabled(sdo));
533
534 m_sdoID->push_back(sdoID_int);
535 m_sdoWord->push_back(sdoWord);
536 m_barrelEndcap_sdo->push_back(sctBrlEc_sdo);
537 m_layerDisk_sdo->push_back(sctLayerDisk_sdo);
538 m_phiModule_sdo->push_back(sctPhiMod_sdo);
539 m_etaModule_sdo->push_back(sctEtaMod_sdo);
540 m_side_sdo->push_back(sctSide_sdo);
541 m_strip_sdo->push_back(sctStrip_sdo);
542 m_row_sdo->push_back(sctRow_sdo);
543 m_noise->push_back(noise);
544 m_belowThresh->push_back(belowThresh);
545 m_disabled->push_back(disabled);
546
547 if (belowThresh) {
548 if (sctBrlEc_sdo==0)
549 m_h_belowThresh_brl->Fill(sctLayerDisk_sdo);
550 else if (abs(sctBrlEc_sdo)==2)
551 m_h_belowThresh_ec->Fill(sctLayerDisk_sdo);
552 }
553
554 if (disabled) {
555 if (sctBrlEc_sdo==0)
556 m_h_disabled_brl->Fill(sctLayerDisk_sdo);
557 else if (abs(sctBrlEc_sdo)==2)
558 m_h_disabled_ec->Fill(sctLayerDisk_sdo);
559 }
560
561 m_h_sdoID->Fill(sdoID_int);
562 m_h_sdoWord->Fill(sdoWord);
563 m_h_barrelEndcap_sdo->Fill(sctBrlEc_sdo);
564 m_h_layerDisk_sdo->Fill(sctLayerDisk_sdo);
565 m_h_phiModule_sdo->Fill(sctPhiMod_sdo);
566 m_h_etaModule_sdo->Fill(sctEtaMod_sdo);
567 m_h_side_sdo->Fill(sctSide_sdo);
568 m_h_strip_sdo->Fill(sctStrip_sdo);
569 m_h_row_sdo->Fill(sctRow_sdo);
570 m_h_phi_v_eta_sdo->Fill(sctEtaMod_sdo, sctPhiMod_sdo);
571
572 // loop over deposits
573 // InDetSimData::Deposit typedef for std::pair<HepMCParticleLink,float>
574 const std::vector<InDetSimData::Deposit>& deposits = sdo.getdeposits();
575 std::vector<InDetSimData::Deposit>::const_iterator dep_itr(deposits.begin());
576 const std::vector<InDetSimData::Deposit>::const_iterator dep_end(deposits.end());
577
578 for ( ; dep_itr != dep_end; ++dep_itr ) {
579 const HepMcParticleLink& particleLink = (*dep_itr).first;
580 const int bar(HepMC::barcode(particleLink)); // FIXME barcode-based
581 const int eventIx(particleLink.eventIndex());
582 const float charge((*dep_itr).second);
583
584 m_barcode->push_back(bar);
585 m_eventIndex->push_back(eventIx);
586 m_charge->push_back(charge);
587
588 m_h_barcode->Fill(bar);
589 m_h_eventIndex->Fill(eventIx);
590 m_h_charge->Fill(charge);
591
592 barcode_vec.push_back(bar);
593 eventIndex_vec.push_back(eventIx);
594 charge_vec.push_back(charge);
595 }
596 m_barcode_vec->push_back(barcode_vec);
597 m_eventIndex_vec->push_back(eventIndex_vec);
598 m_charge_vec->push_back(charge_vec);
599 barcode_vec.clear();
600 eventIndex_vec.clear();
601 charge_vec.clear();
602 }
603 }
604
605 if (m_tree) {
606 m_tree->Fill();
607 }
608
609 return StatusCode::SUCCESS;
610}
611
612} // namespace ITk
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
charge
double charge(const T &p)
Definition AtlasPID.h:997
SCT_DetectorManager.h
SCT_RDO_Collection
InDetRawDataCollection< SCT_RDORawData > SCT_RDO_Collection
Definition SCT_RDO_Collection.h:25
ReadHandle.h
Handle class for reading from StoreGate.
AthAlgorithm::AthAlgorithm
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters:
Definition AthAlgorithm.cxx:24
AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
Definition AthCommonDataStore.h:95
InDetRawDataCollection< SCT_RDORawData >::const_iterator
DataModel_detail::const_iterator< DataVector > const_iterator
Definition DataVector.h:838
DataVector< RawDataT >::end
const_iterator end() const noexcept
DataVector< RawDataT >::begin
const_iterator begin() const noexcept
ITk::StripRDOAnalysis::m_charge
std::vector< float > * m_charge
Definition StripRDOAnalysis.h:104
ITk::StripRDOAnalysis::m_phiModule_sdo
std::vector< int > * m_phiModule_sdo
Definition StripRDOAnalysis.h:92
ITk::StripRDOAnalysis::m_globalX1
std::vector< double > * m_globalX1
Definition StripRDOAnalysis.h:78
ITk::StripRDOAnalysis::m_rdoID
std::vector< unsigned long long > * m_rdoID
Definition StripRDOAnalysis.h:62
ITk::StripRDOAnalysis::m_h_globalZR_shared
LockedHandle< TH2 > m_h_globalZR_shared
Definition StripRDOAnalysis.h:164
ITk::StripRDOAnalysis::m_barrelEndcap
std::vector< int > * m_barrelEndcap
Definition StripRDOAnalysis.h:65
ITk::StripRDOAnalysis::StripRDOAnalysis
StripRDOAnalysis(const std::string &name, ISvcLocator *pSvcLocator)
Definition StripRDOAnalysis.cxx:21
ITk::StripRDOAnalysis::m_sdoWord
std::vector< int > * m_sdoWord
Definition StripRDOAnalysis.h:88
ITk::StripRDOAnalysis::m_eventIndex_vec
std::vector< std::vector< int > > * m_eventIndex_vec
Definition StripRDOAnalysis.h:106
ITk::StripRDOAnalysis::m_barcode
std::vector< int > * m_barcode
Definition StripRDOAnalysis.h:102
ITk::StripRDOAnalysis::m_h_brl_strip_perLayer
std::vector< TH1 * > m_h_brl_strip_perLayer
Definition StripRDOAnalysis.h:158
ITk::StripRDOAnalysis::m_inputTruthKey
SG::ReadHandleKey< InDetSimDataCollection > m_inputTruthKey
Definition StripRDOAnalysis.h:48
ITk::StripRDOAnalysis::m_belowThresh
std::vector< bool > * m_belowThresh
Definition StripRDOAnalysis.h:99
ITk::StripRDOAnalysis::m_phiModule
std::vector< int > * m_phiModule
Definition StripRDOAnalysis.h:67
ITk::StripRDOAnalysis::m_charge_vec
std::vector< std::vector< float > > * m_charge_vec
Definition StripRDOAnalysis.h:107
ITk::StripRDOAnalysis::m_side_sdo
std::vector< int > * m_side_sdo
Definition StripRDOAnalysis.h:94
ITk::StripRDOAnalysis::m_groupSize
std::vector< int > * m_groupSize
Definition StripRDOAnalysis.h:73
ITk::StripRDOAnalysis::m_barrelEndcap_sdo
std::vector< int > * m_barrelEndcap_sdo
Definition StripRDOAnalysis.h:90
ITk::StripRDOAnalysis::m_localZ
std::vector< double > * m_localZ
Definition StripRDOAnalysis.h:83
ITk::StripRDOAnalysis::m_etaModule
std::vector< int > * m_etaModule
Definition StripRDOAnalysis.h:68
ITk::StripRDOAnalysis::m_h_globalXY_shared
LockedHandle< TH2 > m_h_globalXY_shared
Definition StripRDOAnalysis.h:163
ITk::StripRDOAnalysis::m_doPosition
Gaudi::Property< bool > m_doPosition
Definition StripRDOAnalysis.h:57
ITk::StripRDOAnalysis::m_localX
std::vector< double > * m_localX
Definition StripRDOAnalysis.h:81
ITk::StripRDOAnalysis::m_globalZ0
std::vector< double > * m_globalZ0
Definition StripRDOAnalysis.h:77
ITk::StripRDOAnalysis::m_sdoID
std::vector< unsigned long long > * m_sdoID
Definition StripRDOAnalysis.h:87
ITk::StripRDOAnalysis::m_layerDisk
std::vector< int > * m_layerDisk
Definition StripRDOAnalysis.h:66
ITk::StripRDOAnalysis::m_row
std::vector< int > * m_row
Definition StripRDOAnalysis.h:71
ITk::StripRDOAnalysis::m_side
std::vector< int > * m_side
Definition StripRDOAnalysis.h:69
ITk::StripRDOAnalysis::m_inputKey
SG::ReadHandleKey< SCT_RDO_Container > m_inputKey
Definition StripRDOAnalysis.h:47
ITk::StripRDOAnalysis::m_ntuplePath
Gaudi::Property< std::string > m_ntuplePath
Definition StripRDOAnalysis.h:55
ITk::StripRDOAnalysis::m_globalX0
std::vector< double > * m_globalX0
Definition StripRDOAnalysis.h:75
ITk::StripRDOAnalysis::m_row_sdo
std::vector< int > * m_row_sdo
Definition StripRDOAnalysis.h:96
ITk::StripRDOAnalysis::m_barcode_vec
std::vector< std::vector< int > > * m_barcode_vec
Definition StripRDOAnalysis.h:105
ITk::StripRDOAnalysis::m_histPath
Gaudi::Property< std::string > m_histPath
Definition StripRDOAnalysis.h:53
ITk::StripRDOAnalysis::m_ntupleName
Gaudi::Property< std::string > m_ntupleName
Definition StripRDOAnalysis.h:56
ITk::StripRDOAnalysis::m_globalY1
std::vector< double > * m_globalY1
Definition StripRDOAnalysis.h:79
ITk::StripRDOAnalysis::m_inputMcEventCollectionKey
SG::ReadHandleKey< McEventCollection > m_inputMcEventCollectionKey
Definition StripRDOAnalysis.h:49
ITk::StripRDOAnalysis::m_layerDisk_sdo
std::vector< int > * m_layerDisk_sdo
Definition StripRDOAnalysis.h:91
ITk::StripRDOAnalysis::m_rdoWord
std::vector< unsigned int > * m_rdoWord
Definition StripRDOAnalysis.h:63
ITk::StripRDOAnalysis::initialize
virtual StatusCode initialize() override final
Definition StripRDOAnalysis.cxx:26
ITk::StripRDOAnalysis::m_noise
std::vector< bool > * m_noise
Definition StripRDOAnalysis.h:98
ITk::StripRDOAnalysis::m_globalZ1
std::vector< double > * m_globalZ1
Definition StripRDOAnalysis.h:80
ITk::StripRDOAnalysis::m_strip_sdo
std::vector< int > * m_strip_sdo
Definition StripRDOAnalysis.h:95
ITk::StripRDOAnalysis::m_localY
std::vector< double > * m_localY
Definition StripRDOAnalysis.h:82
ITk::StripRDOAnalysis::m_strip
std::vector< int > * m_strip
Definition StripRDOAnalysis.h:70
ITk::StripRDOAnalysis::m_sharedHistPath
Gaudi::Property< std::string > m_sharedHistPath
Definition StripRDOAnalysis.h:54
ITk::StripRDOAnalysis::m_SCT_Manager
const InDetDD::SCT_DetectorManager * m_SCT_Manager
Definition StripRDOAnalysis.h:51
ITk::StripRDOAnalysis::m_etaModule_sdo
std::vector< int > * m_etaModule_sdo
Definition StripRDOAnalysis.h:93
ITk::StripRDOAnalysis::m_disabled
std::vector< bool > * m_disabled
Definition StripRDOAnalysis.h:100
ITk::StripRDOAnalysis::execute
virtual StatusCode execute() override final
Definition StripRDOAnalysis.cxx:306
ITk::StripRDOAnalysis::m_h_ec_strip_perLayer
std::vector< TH1 * > m_h_ec_strip_perLayer
Definition StripRDOAnalysis.h:159
ITk::StripRDOAnalysis::m_globalY0
std::vector< double > * m_globalY0
Definition StripRDOAnalysis.h:76
ITk::StripRDOAnalysis::m_thistSvc
ServiceHandle< ITHistSvc > m_thistSvc
Definition StripRDOAnalysis.h:59
ITk::StripRDOAnalysis::m_eventIndex
std::vector< int > * m_eventIndex
Definition StripRDOAnalysis.h:103
IdentifiableContainerMT::const_iterator
Definition IdentifiableContainerMT.h:80
Identifier::get_compact
value_type get_compact() const
Get the compact id.
InDetDD::SiDetectorElement
Class to hold geometrical description of a silicon detector element.
Definition SiDetectorElement.h:109
InDetDD::SiDetectorElement::endsOfStrip
std::pair< Amg::Vector3D, Amg::Vector3D > endsOfStrip(const Amg::Vector2D &position) const
Special method for SCT to retrieve the two ends of a "strip" Returned coordinates are in global frame...
Definition SiDetectorElement.cxx:339
InDetDD::SolidStateDetectorElementBase::rawLocalPositionOfCell
Amg::Vector2D rawLocalPositionOfCell(const SiCellId &cellId) const
Returns position (center) of cell.
Definition SolidStateDetectorElementBase.cxx:230
InDetSimData::word
int word() const
Definition InDetSimData.h:69
InDetSimData::getdeposits
const std::vector< Deposit > & getdeposits() const
Definition InDetSimData.h:74
SCT_SimHelper::isNoise
static bool isNoise(const InDetSimData &sdo)
Definition SCT_SimHelper.h:93
SCT_SimHelper::isBelowThreshold
static bool isBelowThreshold(const InDetSimData &sdo)
Definition SCT_SimHelper.h:98
SCT_SimHelper::isDisabled
static bool isDisabled(const InDetSimData &sdo)
Definition SCT_SimHelper.h:103
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition GeoPrimitives.h:48
HepMC::barcode
int barcode(const T *p)
Definition Barcode.h:16
HepMC::ConstGenParticlePtr
const GenParticle * ConstGenParticlePtr
Definition GenParticle.h:38
Generated on for ATLAS Offline Software by doxygen 1.14.0