ATLAS Offline Software
Loading...
Searching...
No Matches
TruthHitAnalysis.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "TruthHitAnalysis.h"
6
7// Section of includes for Truth tests
10
11
12
14 ATH_MSG_DEBUG( "Initializing TruthHitAnalysis" );
15
16 // Grab the Ntuple and histogramming service for the tree
17 ATH_CHECK(m_readKey.initialize());
18
20 m_h_n_vert = new TH1D("h_n_vert","n_vert", 100,200, 1500);
21 m_h_n_vert->StatOverflows();
22 ATH_CHECK(histSvc()->regHist(m_path + m_h_n_vert->GetName(), m_h_n_vert));
23
24 m_h_n_part = new TH1D("h_n_part","n_part", 100,1000, 10000);
25 m_h_n_part->StatOverflows();
26 ATH_CHECK(histSvc()->regHist(m_path + m_h_n_part->GetName(), m_h_n_part));
27
28 m_h_n_vert_prim = new TH1D("h_n_vert_prim","n_vert prim", 100,0, 1000);
29 m_h_n_vert_prim->StatOverflows();
30 ATH_CHECK(histSvc()->regHist(m_path + m_h_n_vert_prim->GetName(), m_h_n_vert_prim));
31
32 m_h_n_part_prim = new TH1D("h_n_part_prim","n_part prim", 100,200, 1500);
33 m_h_n_part_prim->StatOverflows();
34 ATH_CHECK(histSvc()->regHist(m_path + m_h_n_part_prim->GetName(), m_h_n_part_prim));
35
36 m_h_n_vert_sec = new TH1D("h_n_vert_sec","n_vert sec", 100,0, 1000);
37 m_h_n_vert_sec->StatOverflows();
38 ATH_CHECK(histSvc()->regHist(m_path + m_h_n_vert_sec->GetName(), m_h_n_vert_sec));
39
40 m_h_n_part_sec = new TH1D("h_n_part_sec","n_part sec", 100,0, 5000);
41 m_h_n_part_sec->StatOverflows();
42 ATH_CHECK(histSvc()->regHist(m_path + m_h_n_part_sec->GetName(), m_h_n_part_sec));
43
44 m_h_vtx_x = new TH1D("h_vtx_x","vtx_x", 100,-1300, 1300);
45 m_h_vtx_x->StatOverflows();
46 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_x->GetName(), m_h_vtx_x));
47
48 m_h_vtx_y = new TH1D("h_vtx_y","vtx_y", 100,-1200, 1200);
49 m_h_vtx_y->StatOverflows();
50 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_y->GetName(), m_h_vtx_y));
51
52 m_h_vtx_z = new TH1D("h_vtx_z","vtx_z", 100,-5000, 5000);
53 m_h_vtx_z->StatOverflows();
54 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_z->GetName(), m_h_vtx_z));
55
56 m_h_vtx_r = new TH1D("h_vtx_r","vtx_r", 100,0, 1160);
57 m_h_vtx_r->StatOverflows();
58 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_r->GetName(), m_h_vtx_r));
59
60 m_h_vtx_proc_x = new TH1D("h_vtx_proc_x","vtx_proc_x", 100,-1300, 1300);
61 m_h_vtx_proc_x->StatOverflows();
62 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_proc_x->GetName(), m_h_vtx_proc_x));
63
64 m_h_vtx_proc_y = new TH1D("h_vtx_proc_y","vtx_proc_y", 100,-1200, 1200);
65 m_h_vtx_proc_y->StatOverflows();
66 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_proc_y->GetName(), m_h_vtx_proc_y));
67
68 m_h_vtx_proc_z = new TH1D("h_vtx_proc_z","vtx_proc_z", 100,-5000, 5000);
69 m_h_vtx_proc_z->StatOverflows();
70 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_proc_z->GetName(), m_h_vtx_proc_z));
71
72 m_h_vtx_proc_r = new TH1D("h_vtx_proc_r","vtx_proc_r", 100,0, 1160);
73 m_h_vtx_proc_r->StatOverflows();
74 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_proc_r->GetName(), m_h_vtx_proc_r));
75
76 m_h_vtx_prim_x = new TH1D("h_vtx_prim_x","vtx_prim_x", 100,-1300, 1300);
77 m_h_vtx_prim_x->StatOverflows();
78 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_x->GetName(), m_h_vtx_prim_x));
79
80 m_h_vtx_prim_y = new TH1D("h_vtx_prim_y","vtx_prim_y", 100,-1200, 1200);
81 m_h_vtx_prim_y->StatOverflows();
82 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_y->GetName(), m_h_vtx_prim_y));
83
84 m_h_vtx_prim_z = new TH1D("h_vtx_prim_z","vtx_prim_z", 100,-5000, 5000);
85 m_h_vtx_prim_z->StatOverflows();
86 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_z->GetName(), m_h_vtx_prim_z));
87
88 m_h_vtx_prim_r = new TH1D("h_vtx_prim_r","vtx_prim_r", 100,0, 1160);
89 m_h_vtx_prim_r->StatOverflows();
90 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_r->GetName(), m_h_vtx_prim_r));
91
92 m_h_vtx_prim_xy = new TH2D("h_vtx_prim_xy","vtx_prim_xy", 100,-100, 100, 100,-100, 100);
93 m_h_vtx_prim_xy->StatOverflows();
94 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_xy->GetName(), m_h_vtx_prim_xy));
95
96 m_h_vtx_prim_zr = new TH2D("h_vtx_prim_zr","vtx_prim_zr", 100,-1500, 1500, 100,0, 150);
97 m_h_vtx_prim_zr->StatOverflows();
98 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_zr->GetName(), m_h_vtx_prim_zr));
99
100 m_h_vtx_sec_xy = new TH2D("h_vtx_sec_xy","vtx_sec_xy", 100,-1200, 1200, 100,-1200, 1200);
101 m_h_vtx_sec_xy->StatOverflows();
102 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_sec_xy->GetName(), m_h_vtx_sec_xy));
103
104 m_h_vtx_sec_zr = new TH2D("h_vtx_sec_zr","vtx_sec_zr", 100,-6000, 6000, 100,0, 1160);
105 m_h_vtx_sec_zr->StatOverflows();
106 ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_sec_zr->GetName(), m_h_vtx_sec_zr));
107
108 m_h_n_generations = new TH1D("h_n_generations","h_generations", 100,0, 25);
109 m_h_n_generations->StatOverflows();
111
112 m_h_truth_px = new TH1D("h_truth_px","truth_px", 100,0, 4000);
113 m_h_truth_px->StatOverflows();
114 ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_px->GetName(), m_h_truth_px));
115
116 m_h_truth_py = new TH1D("h_truth_py","truth_py", 100,0, 4000);
117 m_h_truth_py->StatOverflows();
118 ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_py->GetName(), m_h_truth_py));
119
120 m_h_truth_pz = new TH1D("h_truth_pz","truth_pz", 100,0, 4000);
121 m_h_truth_pz->StatOverflows();
122 ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_pz->GetName(), m_h_truth_pz));
123
124 m_h_truth_pt = new TH1D("h_truth_pt","truth_pt", 100,0, 4000);
125 m_h_truth_pt->StatOverflows();
126 ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_pt->GetName(), m_h_truth_pt));
127
128 m_h_truth_eta = new TH1D("h_truth_eta","truth_eta", 50,-10, 10);
129 m_h_truth_eta->StatOverflows();
130 ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_eta->GetName(), m_h_truth_eta));
131
132 m_h_truth_phi = new TH1D("h_truth_phi","truth_phi", 25,-3.1416, 3.1416);
133 m_h_truth_phi->StatOverflows();
134 ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_phi->GetName(), m_h_truth_phi));
135
136 m_h_barcode = new TH1D("h_truth_barcode","truth_barcode", 100,0, 300000);
137 m_h_barcode->StatOverflows();
138 ATH_CHECK(histSvc()->regHist(m_path + m_h_barcode->GetName(), m_h_barcode));
139
140 m_h_part_status = new TH1D("h_part_status","part status", 100,0,50);
141 m_h_part_status->StatOverflows();
142 ATH_CHECK(histSvc()->regHist(m_path + m_h_part_status->GetName(), m_h_part_status));
143
144 m_h_part_pdgid = new TH1D("h_part_pdgid","part pdgid", 100,-5000, 5000);
145 m_h_part_pdgid->StatOverflows();
146 ATH_CHECK(histSvc()->regHist(m_path + m_h_part_pdgid->GetName(), m_h_part_pdgid));
147
148 m_h_part_pdgid_sec = new TH1D("h_part_pdgid_sec","part pdgid sec", 100,-5000, 5000);
149 m_h_part_pdgid_sec->StatOverflows();
151
152 m_h_part_eta = new TH1D("h_part_eta","part eta", 100,-10, 10);
153 m_h_part_eta->StatOverflows();
154 ATH_CHECK(histSvc()->regHist(m_path + m_h_part_eta->GetName(), m_h_part_eta));
155
156 m_h_part_phi = new TH1D("h_part_phi","part phi", 100,-3.2, 3.2);
157 m_h_part_phi->StatOverflows();
158 ATH_CHECK(histSvc()->regHist(m_path + m_h_part_phi->GetName(), m_h_part_phi));
159
160 m_h_part_p = new TH1D("h_part_p","part p", 100,0, 5000);
161 m_h_part_p->StatOverflows();
162 ATH_CHECK(histSvc()->regHist(m_path + m_h_part_p->GetName(), m_h_part_p));
163
165 m_tree = new TTree("Truth", "Truth");
166 std::string fullNtupleName = "/" + m_ntupleFileName + "/";
167 ATH_CHECK(histSvc()->regTree(fullNtupleName,m_tree));
168
169 m_tree->Branch("vtx_x", &m_vtx_x);
170 m_tree->Branch("vtx_y", &m_vtx_y);
171 m_tree->Branch("vtx_z", &m_vtx_z);
172 m_tree->Branch("vtx_r", &m_vtx_r);
173 m_tree->Branch("vtx_barcode", &m_vtx_barcode);
174 m_tree->Branch("vtx_proc_x", &m_vtx_proc_x);
175 m_tree->Branch("vtx_proc_y", &m_vtx_proc_y);
176 m_tree->Branch("vtx_proc_z", &m_vtx_proc_z);
177 m_tree->Branch("vtx_proc_r", &m_vtx_proc_r);
178 m_tree->Branch("vtx_proc_barcode", &m_vtx_proc_barcode);
179 m_tree->Branch("vtx_prim_x", &m_vtx_prim_x);
180 m_tree->Branch("vtx_prim_y", &m_vtx_prim_y);
181 m_tree->Branch("vtx_prim_z", &m_vtx_prim_z);
182 m_tree->Branch("vtx_prim_r", &m_vtx_prim_r);
183 m_tree->Branch("vtx_prim_barcode", &m_vtx_prim_barcode);
184 m_tree->Branch("truth_px", &m_truth_px);
185 m_tree->Branch("truth_py", &m_truth_py);
186 m_tree->Branch("truth_pz", &m_truth_pz);
187 m_tree->Branch("truth_pt", &m_truth_pt);
188 m_tree->Branch("truth_eta", &m_truth_eta);
189 m_tree->Branch("truth_phi", &m_truth_phi);
190 m_tree->Branch("barcode", &m_barcode);
191 m_tree->Branch("status", &m_status);
192 m_tree->Branch("pdg_id", &m_pdgid);
193
194 return StatusCode::SUCCESS;
195}
196
197
199 ATH_MSG_DEBUG( "In TruthHitAnalysis::execute()" );
200
201 m_vtx_x->clear();
202 m_vtx_y->clear();
203 m_vtx_z->clear();
204 m_vtx_r->clear();
205 m_vtx_barcode->clear();
206 m_vtx_proc_x->clear();
207 m_vtx_proc_y->clear();
208 m_vtx_proc_z->clear();
209 m_vtx_proc_r->clear();
210 m_vtx_proc_barcode->clear();
211 m_vtx_prim_x->clear();
212 m_vtx_prim_y->clear();
213 m_vtx_prim_z->clear();
214 m_vtx_prim_r->clear();
215 m_vtx_prim_barcode->clear();
216 m_truth_px->clear();
217 m_truth_py->clear();
218 m_truth_pz->clear();
219 m_truth_pt->clear();
220 m_truth_eta->clear();
221 m_truth_phi->clear();
222 m_barcode->clear();
223 m_status->clear();
224 m_pdgid->clear();
225
226 const EventContext& ctx{Gaudi::Hive::currentContext()};
227 const McEventCollection* mcCollection{nullptr};
228 ATH_CHECK(SG::get(mcCollection, m_readKey, ctx));
229
230
231 McEventCollection::const_iterator currentGenEventIter = mcCollection->begin();
232 if (currentGenEventIter != mcCollection->end()) {
233 int nvtx = 0;
234 int nvtx_sec=0;
235
236 const auto &procVtx = HepMC::signal_process_vertex(*currentGenEventIter);
237#ifdef HEPMC3
238 const auto &barcodes = (*currentGenEventIter)->attribute<HepMC::GenEventBarcodes> ("barcodes");
239 std::map<int,int> id_to_barcode_map;
240 if (barcodes) id_to_barcode_map = barcodes->id_to_barcode_map();
241 for (const auto& vtx: (*currentGenEventIter)->vertices()) {
242 int bcode = id_to_barcode_map[vtx->id()];
243#else
244 for (HepMC::GenEvent::vertex_const_iterator vtxit=(*currentGenEventIter)->vertices_begin(); vtxit!=(*currentGenEventIter)->vertices_end(); ++vtxit) {
245 auto vtx=*vtxit;
246 int bcode = HepMC::barcode(*vtx);
247#endif
248 double x = vtx->position().x();
249 double y = vtx->position().y();
250 double z = vtx->position().z();
251 double r = std::sqrt(x*x+y*y);
252 m_h_vtx_x->Fill(x);
253 m_h_vtx_y->Fill(y);
254 m_h_vtx_r->Fill(r);
255 m_h_vtx_z->Fill(z);
256
257 m_vtx_x->push_back(x);
258 m_vtx_y->push_back(y);
259 m_vtx_r->push_back(r);
260 m_vtx_z->push_back(z);
261 m_vtx_barcode->push_back(bcode);
262
263 if (vtx == procVtx) {
264 m_h_vtx_proc_x->Fill(x);
265 m_h_vtx_proc_y->Fill(y);
266 m_h_vtx_proc_r->Fill(r);
267 m_h_vtx_proc_z->Fill(z);
268
269 m_vtx_proc_x->push_back(x);
270 m_vtx_proc_y->push_back(y);
271 m_vtx_proc_r->push_back(r);
272 m_vtx_proc_z->push_back(z);
273 m_vtx_proc_barcode->push_back(bcode);
274 }
275
276 if (!HepMC::is_simulation_vertex(vtx)) {
277 m_h_vtx_prim_x->Fill(x);
278 m_h_vtx_prim_y->Fill(y);
279 m_h_vtx_prim_r->Fill(r);
280 m_h_vtx_prim_z->Fill(z);
281 m_h_vtx_prim_xy->Fill(x,y);
282 m_h_vtx_prim_zr->Fill(z,r);
283
284 m_vtx_prim_x->push_back(x);
285 m_vtx_prim_y->push_back(y);
286 m_vtx_prim_r->push_back(r);
287 m_vtx_prim_z->push_back(z);
288 m_vtx_prim_barcode->push_back(bcode);
289 ++nvtx;
290 }
291 else {
292 m_h_vtx_sec_xy->Fill(x,y);
293 m_h_vtx_sec_zr->Fill(z,r);
294 ++nvtx_sec;
295 }
296 } //End iteration over vertices
297
298 m_h_n_vert->Fill(nvtx+nvtx_sec);
299 m_h_n_vert_prim->Fill(nvtx);
300 m_h_n_vert_sec->Fill(nvtx_sec);
301
302 int npart_prim=0;
303 int npart_sec=0;
304
305 for (auto currentGenParticle: *(*currentGenEventIter)) {
306 const HepMC::FourVector mom = currentGenParticle->momentum();
307#ifdef HEPMC3
308 int currentGenParticlebarcode = id_to_barcode_map[currentGenParticle->id()];
309#else
310 int currentGenParticlebarcode=HepMC::barcode(currentGenParticle);
311#endif
312 m_h_truth_px->Fill(mom.x());
313 m_h_truth_py->Fill(mom.y());
314 m_h_truth_pz->Fill(mom.z());
315 m_h_truth_pt->Fill(mom.perp());
316 m_h_truth_eta->Fill(mom.eta());
317 m_h_truth_phi->Fill(mom.phi());
318 m_h_barcode->Fill(currentGenParticlebarcode);
319 m_h_part_status->Fill(currentGenParticle->status());
320 m_truth_px->push_back(mom.x());
321 m_truth_py->push_back(mom.y());
322 m_truth_pz->push_back(mom.z());
323 m_truth_pt->push_back(mom.perp());
324 m_truth_eta->push_back(mom.eta());
325 m_truth_phi->push_back(mom.phi());
326 m_barcode->push_back(currentGenParticlebarcode);
327 m_status->push_back(currentGenParticle->status());
328
329 int pdg = currentGenParticle->pdg_id();
330 m_pdgid->push_back(pdg);
331
332 if (!HepMC::is_simulation_particle(currentGenParticle)) {
333 m_h_part_pdgid->Fill(pdg);
334 m_h_part_p->Fill(std::sqrt(mom.x()*mom.x()+mom.y()*mom.y()+mom.z()*mom.z()));
335 m_h_part_eta->Fill(mom.eta());
336 m_h_part_phi->Fill(mom.phi());
337 ++npart_prim;
338 m_h_n_generations->Fill(1);
339/* The conditions don't make any sense.
340 if (currentGenParticlebarcode < HepMC::PHOTOSMIN) {
341 m_h_n_generations->Fill(0);
342 }
343 else {
344 m_h_n_generations->Fill(1);
345 }
346*/
347 } //End simulation particle
348 else {
349 m_h_part_pdgid_sec->Fill(pdg);
350 ++npart_sec;
351 const int gen = HepMC::generations(currentGenParticle) + 2;
352 m_h_n_generations->Fill(gen);
353 }
354 } // End iteration over particles
355
356 m_h_n_part_prim->Fill(npart_prim);
357 m_h_n_part_sec->Fill(npart_sec);
358 m_h_n_part->Fill(npart_prim+npart_sec);
359 } // End mcCollection
360
361 m_tree->Fill();
362
363 return StatusCode::SUCCESS;
364}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_DEBUG(x)
#define y
#define x
#define z
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_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
This defines the McEventCollection, which is really just an ObjectVector of McEvent objectsFile: Gene...
std::vector< float > * m_vtx_z
std::vector< float > * m_truth_pt
virtual StatusCode initialize() override final
std::vector< float > * m_truth_py
std::vector< float > * m_vtx_prim_barcode
std::vector< float > * m_vtx_proc_r
std::vector< float > * m_pdgid
std::vector< float > * m_vtx_y
Gaudi::Property< std::string > m_path
std::vector< float > * m_vtx_prim_z
TH1 * m_h_n_vert
Some variables.
virtual StatusCode execute() override final
std::vector< float > * m_vtx_proc_y
std::vector< float > * m_vtx_prim_x
std::vector< float > * m_vtx_x
SG::ReadHandleKey< McEventCollection > m_readKey
std::vector< float > * m_vtx_proc_x
std::vector< float > * m_vtx_barcode
std::vector< float > * m_vtx_proc_z
std::vector< float > * m_vtx_prim_r
std::vector< float > * m_vtx_prim_y
std::vector< float > * m_truth_pz
Gaudi::Property< std::string > m_ntupleFileName
std::vector< float > * m_truth_eta
std::vector< float > * m_barcode
std::vector< float > * m_truth_phi
std::vector< float > * m_status
std::vector< float > * m_vtx_r
std::vector< float > * m_truth_px
std::vector< float > * m_vtx_proc_barcode
int r
Definition globals.cxx:22
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
int barcode(const T *p)
Definition Barcode.h:16
bool is_simulation_particle(const T &p)
Method to establish if a particle (or barcode) was created during the simulation (TODO update to be s...
int generations(const T &p)
Method to return how many interactions a particle has undergone during simulation (TODO migrate to be...
GenVertex * signal_process_vertex(const GenEvent *e)
Definition GenEvent.h:625
const T * get(const ReadCondHandleKey< T > &key, const EventContext &ctx)
Convenience function to retrieve an object given a ReadCondHandleKey.