ATLAS Offline Software
SiHitAnalysis.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 "SiHitAnalysis.h"
6 
7 
8 #include "StoreGate/ReadHandle.h"
9 #include "GeoAdaptors/GeoSiHit.h"
11 #include "AtlasHepMC/GenVertex.h"
12 #include "AtlasHepMC/GenParticle.h"
13 
14 #include "TH1.h"
15 #include "TH2.h"
16 #include "TTree.h"
17 
18 
20 {
21  ATH_MSG_DEBUG("Initializing SiHitAnalysis");
22 
23  std::string detName("Pixel");
24  std::string ntupName("SiPixel");
25 
26  if (m_hitsContainerKey.key()=="PixelHits") {
27  detName = "Pixel";
28  ntupName = "SiPixel";
29  }
30  else if (m_hitsContainerKey.key()=="ITkPixelHits") {
31  detName = "ITkPixel";
32  ntupName = "SiITkPixel";
33  }
34  else if (m_hitsContainerKey.key()=="PLR_Hits") {
35  detName = "PLR";
36  ntupName = "SiPLR";
37  }
38  else if (m_hitsContainerKey.key()=="SCT_Hits") {
39  detName = "SCT";
40  ntupName = "SiSCT";
41  }
42  else if (m_hitsContainerKey.key()=="ITkStripHits") {
43  detName = "ITkStrip";
44  ntupName = "SiITkStrip";
45  }
46  else if (m_hitsContainerKey.key()=="BCMHits") {
47  detName = "BCM";
48  ntupName = "SiBCM";
49  }
50  else if (m_hitsContainerKey.key()=="BLMHits") {
51  detName = "BLM";
52  ntupName = "SiBLM";
53  }
54  else if (m_hitsContainerKey.key()=="HGTD_Hits") {
55  detName = "HGTD";
56  ntupName = "SiHGTD";
57  }
58  else {
59  ATH_MSG_ERROR("SiHitsAnalysis for " << m_hitsContainerKey.key() << " not supported!!!");
60  return StatusCode::FAILURE;
61  }
62 
64  float bin_down = -600;
65  float bin_up = 600;
66  float radius_up = 600;
67  float radius_down = 200;
68  float z_max = 1200;
69  if (detName == "Pixel") {
70  bin_down = -170;
71  bin_up = 170;
72  radius_up = 170;
73  radius_down = 0;
74  } else if (detName == "ITkPixel") {
75  bin_down = -350;
76  bin_up = 350;
77  radius_up = 350;
78  radius_down = 0;
79  z_max = 3400;
80  } else if (detName == "ITkStrip") {
81  bin_down = -1100;
82  bin_up = 1100;
83  radius_up = 1100;
84  radius_down = 0;
85  z_max = 3400;
86  } else if (detName == "HGTD") {
87  bin_down = -1000;
88  bin_up = 1000;
89  radius_up = 700;
90  radius_down = 0;
91  z_max = 4000;
92  } else if (detName == "PLR") {
93  bin_down = -125;
94  bin_up = 125;
95  radius_up = 125;
96  radius_down = 0;
97  z_max = 3000;
98  }
99  m_h_hits_x = new TH1D(("h_"+detName+"_x").c_str(),("h_"+detName+"_x").c_str(), 100,bin_down, bin_up);
100  m_h_hits_x->StatOverflows();
101 
102  m_h_hits_y = new TH1D(("h_"+detName+"_y").c_str(), ("h_"+detName+"_y").c_str(), 100,bin_down,bin_up);
103  m_h_hits_y->StatOverflows();
104 
105  m_h_hits_z = new TH1D(("h_"+detName+"_z").c_str(), ("h_"+detName+"_z").c_str(), 200,-z_max,z_max);
106  m_h_hits_z->StatOverflows();
107 
108  m_h_hits_r = new TH1D(("h_"+detName+"_r").c_str(), ("h_"+detName+"_r").c_str(), 100,radius_down,radius_up);
109  m_h_hits_r->StatOverflows();
110 
111  m_h_xy = new TH2D(("h_"+detName+"_xy").c_str(), ("h_"+detName+"_xy").c_str(), 100,bin_down,bin_up,100, bin_down, bin_up);
112  m_h_xy->StatOverflows();
113 
114  m_h_zr = new TH2D(("h_"+detName+"_zr").c_str(), ("h_"+detName+"_zr").c_str(), 100,-z_max, z_max, 100, radius_down, radius_up);
115  m_h_zr->StatOverflows();
116 
117  m_h_hits_time = new TH1D(("h_"+detName+"_time").c_str(), ("h_"+detName+"_time").c_str(), 100,0,500);
118  m_h_hits_time->StatOverflows();
119 
120  m_h_hits_eloss = new TH1D(("h_"+detName+"_eloss").c_str(), ("h_"+detName+"_eloss").c_str(), 100,0,50);
121  m_h_hits_eloss->StatOverflows();
122 
123  m_h_hits_step = new TH1D(("h_"+detName+"_step").c_str(), ("h_"+detName+"_step").c_str(), 100,0,50);
124  m_h_hits_step->StatOverflows();
125 
126  m_h_hits_barcode = new TH1D(("h_"+detName+"_barcode").c_str(), ("h_"+detName+"_barcode").c_str(), 200,0,250000);
127  m_h_hits_barcode->StatOverflows();
128 
129  m_h_time_eloss = new TH2D(("h_"+detName+"_time_eloss").c_str(), ("h_"+detName+" Eloss vs. time").c_str(),100, 0,500,100,0,50);
130  m_h_time_eloss->StatOverflows();
131 
132  m_h_z_eloss = new TH2D(("h_"+detName+"_z_eloss").c_str(), ("h_"+detName+" Eloss vs. z").c_str(),100, -z_max,z_max,100,0,50);
133  m_h_z_eloss->StatOverflows();
134 
135  m_h_r_eloss = new TH2D(("h_"+detName+"_r_eloss").c_str(), ("h_"+detName+ " Eloss vs. r").c_str(),100, radius_down,radius_down,100,0,50);
136  m_h_r_eloss->StatOverflows();
137 
138  m_h_barrel_endcap = new TH1D(("h_"+detName+"_barrel_endcap").c_str(), ("h_"+detName+ " barrel/endcap").c_str(), 10, -5, 5);
139  m_h_barrel_endcap->StatOverflows();
140 
141  m_h_layer_disk = new TH1D(("h_"+detName+"_layer_disk").c_str(), ("h_"+detName+ " layer/disk").c_str(), 10, 0, 10);
142  m_h_layer_disk->StatOverflows();
143 
144  m_h_module_eta = new TH1D(("h_"+detName+"_module_eta").c_str(), ("h_"+detName+ " module in #eta").c_str(), 100, 0, 100);
145  m_h_module_eta->StatOverflows();
146 
147  m_h_module_phi = new TH1D(("h_"+detName+"_module_phi").c_str(), ("h_"+detName+ " module in #phi").c_str(), 100, 0, 100);
148  m_h_module_phi->StatOverflows();
149 
150  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_x->GetName(), m_h_hits_x));
151  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_y->GetName(), m_h_hits_y));
152  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_z->GetName(), m_h_hits_z));
153  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_r->GetName(), m_h_hits_r));
154  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_xy->GetName(), m_h_xy));
155  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_zr->GetName(), m_h_zr));
156  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_time->GetName(), m_h_hits_time));
157  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_eloss->GetName(), m_h_hits_eloss));
158  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_step->GetName(), m_h_hits_step));
159  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_hits_barcode->GetName(), m_h_hits_barcode));
160 
161  //To be filled only when the expert mode is on.
162  if (m_expert.value()) {
163  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_time_eloss->GetName(), m_h_time_eloss));
164  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_z_eloss->GetName(), m_h_z_eloss));
165  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_r_eloss->GetName(), m_h_r_eloss));
166  }
167 
169  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_layer_disk->GetName(), m_h_layer_disk));
170  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_module_eta->GetName(), m_h_module_eta));
171  ATH_CHECK(histSvc()->regHist(m_histPath + m_h_module_phi->GetName(), m_h_module_phi));
172 
173  // Special shared ITk histograms
174  if (detName.find("ITk") != std::string::npos) {
175  std::string xy_name = "h_ITk_xy";
176  auto xy = std::make_unique<TH2D>(xy_name.c_str(), xy_name.c_str(), 2200, -1100, 1100, 2200, -1100, 1100);
177  xy->StatOverflows();
178  ATH_CHECK(histSvc()->regShared(m_histPath + xy_name, std::move(xy), m_h_xy_shared));
179 
180  std::string zr_name = "h_ITk_zr";
181  auto zr = std::make_unique<TH2D>(zr_name.c_str(), zr_name.c_str(), 6800, -3400, 3400, 1100, 0, 1100);
182  zr->StatOverflows();
183  ATH_CHECK(histSvc()->regShared(m_histPath + zr_name, std::move(zr), m_h_zr_shared));
184  }
185 
187  m_tree = new TTree(ntupName.c_str(), ntupName.c_str());
188  std::string fullNtupleName = "/" + m_ntuplePath + "/" + detName;
189  ATH_CHECK(histSvc()->regTree(fullNtupleName,m_tree));
190 
191  if (m_tree){
192  m_tree->Branch((detName+"_x").c_str(), &m_hits_x);
193  m_tree->Branch((detName+"_y").c_str(), &m_hits_y);
194  m_tree->Branch((detName+"_z").c_str(), &m_hits_z);
195  m_tree->Branch((detName+"_r").c_str(), &m_hits_r);
196  m_tree->Branch((detName+"_time").c_str(), &m_hits_time);
197  m_tree->Branch((detName+"_eloss").c_str(), &m_hits_eloss);
198  m_tree->Branch((detName+"_step").c_str(), &m_hits_step);
199  m_tree->Branch((detName+"_barcode").c_str(), &m_hits_barcode);
200  if (m_extraTruthBranches) {
201  m_tree->Branch((detName+"_pdgId").c_str(), &m_hits_pdgId);
202  m_tree->Branch((detName+"_pT").c_str(), &m_hits_pT);
203  m_tree->Branch((detName+"_eta").c_str(), &m_hits_eta);
204  m_tree->Branch((detName+"_phi").c_str(), &m_hits_phi);
205  m_tree->Branch((detName+"_prodVtx_x").c_str(), &m_hits_prodVtx_x);
206  m_tree->Branch((detName+"_prodVtx_y").c_str(), &m_hits_prodVtx_y);
207  m_tree->Branch((detName+"_prodVtx_z").c_str(), &m_hits_prodVtx_z);
208  }
209 
210  m_tree->Branch((detName+"_barrel_endcap").c_str(), &m_barrel_endcap);
211  m_tree->Branch((detName+"_layer_disk").c_str(), &m_layer_disk);
212  m_tree->Branch((detName+"_module_eta").c_str(), &m_module_eta);
213  m_tree->Branch((detName+"_module_phi").c_str(), &m_module_phi);
214  }
215  else {
216  ATH_MSG_ERROR("No tree found!");
217  }
218 
219  // Initialize keys
220  ATH_CHECK(m_hitsContainerKey.initialize());
221 
222  return StatusCode::SUCCESS;
223 }
224 
225 
227 {
228  ATH_MSG_DEBUG("In SiHitAnalysis::execute()");
229 
230  m_hits_x->clear();
231  m_hits_y->clear();
232  m_hits_z->clear();
233  m_hits_r->clear();
234  m_hits_time->clear();
235  m_hits_eloss->clear();
236  m_hits_step->clear();
237  m_hits_barcode->clear();
238  if (m_extraTruthBranches) {
239  m_hits_pdgId->clear();
240  m_hits_pT->clear();
241  m_hits_eta->clear();
242  m_hits_phi->clear();
243  m_hits_prodVtx_x->clear();
244  m_hits_prodVtx_y->clear();
245  m_hits_prodVtx_z->clear();
246  }
247 
248  m_barrel_endcap->clear();
249  m_layer_disk->clear();
250  m_module_eta->clear();
251  m_module_phi->clear();
252 
253  const EventContext&ctx {Gaudi::Hive::currentContext()};
254  const SiHitCollection* hitCollection{nullptr};
255  ATH_CHECK(SG::get(hitCollection, m_hitsContainerKey, ctx));
256  ATH_MSG_INFO("Event contains " << hitCollection->size() << " entries in " << m_hitsContainerKey.key());
257  for (const SiHit &hit : *hitCollection) {
258  GeoSiHit ghit(hit);
259  HepGeom::Point3D<double> p = ghit.getGlobalPosition();
260  m_h_hits_x->Fill(p.x());
261  m_h_hits_y->Fill(p.y());
262  m_h_hits_z->Fill(p.z());
263  m_h_hits_r->Fill(p.perp());
264  m_h_xy->Fill(p.x(), p.y());
265  m_h_zr->Fill(p.z(), p.perp());
266  if (m_h_xy_shared.get() != nullptr) {
267  m_h_xy_shared->Fill(p.x(), p.y());
268  }
269  if (m_h_zr_shared.get() != nullptr) {
270  m_h_zr_shared->Fill(p.z(), p.perp());
271  }
272  m_h_hits_eloss->Fill(hit.energyLoss());
273  m_h_hits_time->Fill(hit.meanTime());
274  double step_length = (hit.localStartPosition() - hit.localEndPosition()).mag();
275  m_h_hits_step->Fill(step_length);
276  m_h_hits_barcode->Fill(HepMC::barcode(hit.particleLink()));
277 
278  if (m_expert.value()) {
279  m_h_time_eloss->Fill(hit.meanTime(), hit.energyLoss());
280  if (hit.getBarrelEndcap() == 0) {
281  m_h_z_eloss->Fill(p.z(), hit.energyLoss());
282  }
283  else {
284  m_h_r_eloss->Fill(p.perp(), hit.energyLoss());
285  }
286  }
287 
288  m_h_barrel_endcap->Fill(hit.getBarrelEndcap());
289  m_h_layer_disk->Fill(hit.getLayerDisk());
290  m_h_module_eta->Fill(hit.getEtaModule());
291  m_h_module_phi->Fill(hit.getPhiModule());
292 
293  m_hits_x->push_back(p.x());
294  m_hits_y->push_back(p.y());
295  m_hits_z->push_back(p.z());
296  m_hits_r->push_back(p.perp());
297  m_hits_eloss->push_back(hit.energyLoss());
298  m_hits_time->push_back(hit.meanTime());
299  m_hits_step->push_back(step_length);
300  m_hits_barcode->push_back(HepMC::barcode(hit.particleLink()));
301  if (m_extraTruthBranches) {
302  const auto& HMPL = hit.particleLink();
303  if (HMPL.isValid()) {
304  auto part = HMPL.cptr();
305  m_hits_pdgId->push_back(part->pdg_id());
306  m_hits_pT->push_back(part->momentum().perp());
307  m_hits_eta->push_back(part->momentum().eta());
308  m_hits_phi->push_back(part->momentum().phi());
309  m_hits_prodVtx_x->push_back(part->production_vertex()->position().x());
310  m_hits_prodVtx_y->push_back(part->production_vertex()->position().y());
311  m_hits_prodVtx_z->push_back(part->production_vertex()->position().z());
312  }
313  else {
314  m_hits_pdgId->push_back(-9999);
315  m_hits_pT->push_back(-9999);
316  m_hits_eta->push_back(-9999);
317  m_hits_phi->push_back(-9999);
318  m_hits_prodVtx_x->push_back(-9999);
319  m_hits_prodVtx_y->push_back(-9999);
320  m_hits_prodVtx_z->push_back(-9999);
321  }
322  }
323 
324  m_barrel_endcap->push_back(hit.getBarrelEndcap());
325  m_layer_disk->push_back(hit.getLayerDisk());
326  m_module_eta->push_back(hit.getEtaModule());
327  m_module_phi->push_back(hit.getPhiModule());
328  } // End while hits
329 
330  m_tree->Fill();
331 
332  return StatusCode::SUCCESS;
333 }
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
SiHitAnalysis::m_expert
Gaudi::Property< bool > m_expert
Definition: SiHitAnalysis.h:78
AthHistogramAlgorithm::histSvc
const ServiceHandle< ITHistSvc > & histSvc() const
The standard THistSvc (for writing histograms and TTrees and more to a root file) Returns (kind of) a...
Definition: AthHistogramAlgorithm.h:113
SiHitAnalysis::m_h_hits_y
TH1 * m_h_hits_y
Definition: SiHitAnalysis.h:34
SiHitAnalysis::m_hits_x
std::vector< float > * m_hits_x
Definition: SiHitAnalysis.h:54
SiHitAnalysis::m_module_phi
std::vector< int > * m_module_phi
Definition: SiHitAnalysis.h:71
SiHitAnalysis::m_hits_y
std::vector< float > * m_hits_y
Definition: SiHitAnalysis.h:55
SiHitAnalysis::m_hits_eta
std::vector< float > * m_hits_eta
Definition: SiHitAnalysis.h:64
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
SiHitAnalysis::m_hits_pT
std::vector< float > * m_hits_pT
Definition: SiHitAnalysis.h:63
SiHitAnalysis::m_barrel_endcap
std::vector< int > * m_barrel_endcap
Definition: SiHitAnalysis.h:69
GenVertex.h
SiHitAnalysis::m_hits_prodVtx_y
std::vector< float > * m_hits_prodVtx_y
Definition: SiHitAnalysis.h:67
SiHitAnalysis::m_h_hits_eloss
TH1 * m_h_hits_eloss
Definition: SiHitAnalysis.h:42
AtlasHitsVector< SiHit >
SiHitAnalysis::m_module_eta
std::vector< int > * m_module_eta
Definition: SiHitAnalysis.h:72
SiHitAnalysis::m_hits_z
std::vector< float > * m_hits_z
Definition: SiHitAnalysis.h:56
SiHitAnalysis::m_extraTruthBranches
Gaudi::Property< bool > m_extraTruthBranches
Definition: SiHitAnalysis.h:79
SiHitAnalysis::m_hits_barcode
std::vector< float > * m_hits_barcode
Definition: SiHitAnalysis.h:61
GenParticle.h
GeoSiHit
Definition: GeoSiHit.h:30
SiHit
Definition: SiHit.h:19
SiHitAnalysis::m_hits_pdgId
std::vector< float > * m_hits_pdgId
Definition: SiHitAnalysis.h:62
SiHitAnalysis::m_h_barrel_endcap
TH1 * m_h_barrel_endcap
Definition: SiHitAnalysis.h:48
SiHitAnalysis::m_hitsContainerKey
SG::ReadHandleKey< SiHitCollection > m_hitsContainerKey
Definition: SiHitAnalysis.h:74
SiHitAnalysis::m_h_r_eloss
TH2 * m_h_r_eloss
Definition: SiHitAnalysis.h:47
SiHitAnalysis::m_hits_phi
std::vector< float > * m_hits_phi
Definition: SiHitAnalysis.h:65
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
GeoSiHit.h
SiHitAnalysis::m_h_hits_x
TH1 * m_h_hits_x
Some variables.
Definition: SiHitAnalysis.h:33
SiHitAnalysis.h
SiHitAnalysis::m_hits_r
std::vector< float > * m_hits_r
Definition: SiHitAnalysis.h:57
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
SG::get
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.
Definition: ReadCondHandle.h:287
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
SiHitAnalysis::execute
virtual StatusCode execute() override
Definition: SiHitAnalysis.cxx:226
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
SiHitAnalysis::m_h_hits_time
TH1 * m_h_hits_time
Definition: SiHitAnalysis.h:41
SiHitAnalysis::m_h_module_phi
TH1 * m_h_module_phi
Definition: SiHitAnalysis.h:50
SiHitAnalysis::m_h_hits_z
TH1 * m_h_hits_z
Definition: SiHitAnalysis.h:35
SiHitAnalysis::initialize
virtual StatusCode initialize() override
Definition: SiHitAnalysis.cxx:19
SiHitAnalysis::m_hits_prodVtx_z
std::vector< float > * m_hits_prodVtx_z
Definition: SiHitAnalysis.h:68
SiHitAnalysis::m_hits_eloss
std::vector< float > * m_hits_eloss
Definition: SiHitAnalysis.h:59
SiHitAnalysis::m_tree
TTree * m_tree
Definition: SiHitAnalysis.h:53
SiHitAnalysis::m_h_hits_barcode
TH1 * m_h_hits_barcode
Definition: SiHitAnalysis.h:44
SiHitAnalysis::m_h_hits_r
TH1 * m_h_hits_r
Definition: SiHitAnalysis.h:36
SiHitAnalysis::m_hits_time
std::vector< float > * m_hits_time
Definition: SiHitAnalysis.h:58
SiHitAnalysis::m_layer_disk
std::vector< int > * m_layer_disk
Definition: SiHitAnalysis.h:70
SiHitAnalysis::m_h_zr_shared
LockedHandle< TH2 > m_h_zr_shared
Definition: SiHitAnalysis.h:40
SiHitAnalysis::m_h_xy_shared
LockedHandle< TH2 > m_h_xy_shared
Definition: SiHitAnalysis.h:39
SiHitAnalysis::m_h_zr
TH2 * m_h_zr
Definition: SiHitAnalysis.h:38
GeoSiHit::getGlobalPosition
HepGeom::Point3D< double > getGlobalPosition() const
SiHitAnalysis::m_hits_prodVtx_x
std::vector< float > * m_hits_prodVtx_x
Definition: SiHitAnalysis.h:66
SiHitAnalysis::m_h_xy
TH2 * m_h_xy
Definition: SiHitAnalysis.h:37
ReadHandle.h
Handle class for reading from StoreGate.
python.LArCondContChannels.detName
detName
Definition: LArCondContChannels.py:665
SiHitAnalysis::m_histPath
Gaudi::Property< std::string > m_histPath
Definition: SiHitAnalysis.h:76
SiHitAnalysis::m_hits_step
std::vector< float > * m_hits_step
Definition: SiHitAnalysis.h:60
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
SiHitAnalysis::m_h_time_eloss
TH2 * m_h_time_eloss
Definition: SiHitAnalysis.h:45
SiHitAnalysis::m_h_layer_disk
TH1 * m_h_layer_disk
Definition: SiHitAnalysis.h:49
SiHitAnalysis::m_ntuplePath
Gaudi::Property< std::string > m_ntuplePath
Definition: SiHitAnalysis.h:77
SiHitAnalysis::m_h_z_eloss
TH2 * m_h_z_eloss
Definition: SiHitAnalysis.h:46
SiHitAnalysis::m_h_hits_step
TH1 * m_h_hits_step
Definition: SiHitAnalysis.h:43
SiHitAnalysis::m_h_module_eta
TH1 * m_h_module_eta
Definition: SiHitAnalysis.h:51