ATLAS Offline Software
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
8 #include "AtlasHepMC/GenEvent.h"
9 #include "TruthUtils/MagicNumbers.h"
10 
11 
12 
13 StatusCode TruthHitAnalysis::initialize() {
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_prim_xy = new TH2D("h_vtx_prim_xy","vtx_prim_xy", 100,-100, 100, 100,-100, 100);
61  m_h_vtx_prim_xy->StatOverflows();
62  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_xy->GetName(), m_h_vtx_prim_xy));
63 
64  m_h_vtx_prim_zr = new TH2D("h_vtx_prim_zr","vtx_prim_zr", 100,-1500, 1500, 100,0, 150);
65  m_h_vtx_prim_zr->StatOverflows();
66  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_prim_zr->GetName(), m_h_vtx_prim_zr));
67 
68  m_h_vtx_sec_xy = new TH2D("h_vtx_sec_xy","vtx_sec_xy", 100,-1200, 1200, 100,-1200, 1200);
69  m_h_vtx_sec_xy->StatOverflows();
70  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_sec_xy->GetName(), m_h_vtx_sec_xy));
71 
72  m_h_vtx_sec_zr = new TH2D("h_vtx_sec_zr","vtx_sec_zr", 100,-6000, 6000, 100,0, 1160);
73  m_h_vtx_sec_zr->StatOverflows();
74  ATH_CHECK(histSvc()->regHist(m_path + m_h_vtx_sec_zr->GetName(), m_h_vtx_sec_zr));
75 
76  m_h_n_generations = new TH1D("h_n_generations","h_generations", 100,0, 25);
77  m_h_n_generations->StatOverflows();
78  ATH_CHECK(histSvc()->regHist(m_path + m_h_n_generations->GetName(), m_h_n_generations));
79 
80  m_h_truth_px = new TH1D("h_truth_px","truth_px", 100,0, 4000);
81  m_h_truth_px->StatOverflows();
82  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_px->GetName(), m_h_truth_px));
83 
84  m_h_truth_py = new TH1D("h_truth_py","truth_py", 100,0, 4000);
85  m_h_truth_py->StatOverflows();
86  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_py->GetName(), m_h_truth_py));
87 
88  m_h_truth_pz = new TH1D("h_truth_pz","truth_pz", 100,0, 4000);
89  m_h_truth_pz->StatOverflows();
90  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_pz->GetName(), m_h_truth_pz));
91 
92  m_h_truth_pt = new TH1D("h_truth_pt","truth_pt", 100,0, 4000);
93  m_h_truth_pt->StatOverflows();
94  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_pt->GetName(), m_h_truth_pt));
95 
96  m_h_truth_eta = new TH1D("h_truth_eta","truth_eta", 50,-10, 10);
97  m_h_truth_eta->StatOverflows();
98  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_eta->GetName(), m_h_truth_eta));
99 
100  m_h_truth_phi = new TH1D("h_truth_phi","truth_phi", 25,-3.1416, 3.1416);
101  m_h_truth_phi->StatOverflows();
102  ATH_CHECK(histSvc()->regHist(m_path + m_h_truth_phi->GetName(), m_h_truth_phi));
103 
104  m_h_barcode = new TH1D("h_truth_barcode","truth_barcode", 100,0, 300000);
105  m_h_barcode->StatOverflows();
106  ATH_CHECK(histSvc()->regHist(m_path + m_h_barcode->GetName(), m_h_barcode));
107 
108  m_h_part_status = new TH1D("h_part_status","part status", 100,0,50);
109  m_h_part_status->StatOverflows();
110  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_status->GetName(), m_h_part_status));
111 
112  m_h_part_pdgid = new TH1D("h_part_pdgid","part pdgid", 100,-5000, 5000);
113  m_h_part_pdgid->StatOverflows();
114  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_pdgid->GetName(), m_h_part_pdgid));
115 
116  m_h_part_pdgid_sec = new TH1D("h_part_pdgid_sec","part pdgid sec", 100,-5000, 5000);
117  m_h_part_pdgid_sec->StatOverflows();
118  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_pdgid_sec->GetName(), m_h_part_pdgid_sec));
119 
120  m_h_part_eta = new TH1D("h_part_eta","part eta", 100,-10, 10);
121  m_h_part_eta->StatOverflows();
122  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_eta->GetName(), m_h_part_eta));
123 
124  m_h_part_phi = new TH1D("h_part_phi","part phi", 100,-3.2, 3.2);
125  m_h_part_phi->StatOverflows();
126  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_phi->GetName(), m_h_part_phi));
127 
128  m_h_part_p = new TH1D("h_part_p","part p", 100,0, 5000);
129  m_h_part_p->StatOverflows();
130  ATH_CHECK(histSvc()->regHist(m_path + m_h_part_p->GetName(), m_h_part_p));
131 
133  m_tree = new TTree("Truth", "Truth");
134  std::string fullNtupleName = "/" + m_ntupleFileName + "/";
135  ATH_CHECK(histSvc()->regTree(fullNtupleName,m_tree));
136 
137  m_tree->Branch("vtx_x", &m_vtx_x);
138  m_tree->Branch("vtx_y", &m_vtx_y);
139  m_tree->Branch("vtx_z", &m_vtx_z);
140  m_tree->Branch("vtx_r", &m_vtx_r);
141  m_tree->Branch("vtx_barcode", &m_vtx_barcode);
142  m_tree->Branch("truth_px", &m_truth_px);
143  m_tree->Branch("truth_py", &m_truth_py);
144  m_tree->Branch("truth_pz", &m_truth_pz);
145  m_tree->Branch("truth_pt", &m_truth_pt);
146  m_tree->Branch("truth_eta", &m_truth_eta);
147  m_tree->Branch("truth_phi", &m_truth_phi);
148  m_tree->Branch("barcode", &m_barcode);
149  m_tree->Branch("status", &m_status);
150  m_tree->Branch("pdg_id", &m_pdgid);
151 
152  return StatusCode::SUCCESS;
153 }
154 
155 
156 StatusCode TruthHitAnalysis::execute() {
157  ATH_MSG_DEBUG( "In TruthHitAnalysis::execute()" );
158 
159  m_vtx_x->clear();
160  m_vtx_y->clear();
161  m_vtx_z->clear();
162  m_vtx_r->clear();
163  m_vtx_barcode->clear();
164  m_truth_px->clear();
165  m_truth_py->clear();
166  m_truth_pz->clear();
167  m_truth_pt->clear();
168  m_truth_eta->clear();
169  m_truth_phi->clear();
170  m_barcode->clear();
171  m_status->clear();
172  m_pdgid->clear();
173 
174  const EventContext& ctx{Gaudi::Hive::currentContext()};
175  const McEventCollection* mcCollection{nullptr};
176  ATH_CHECK(SG::get(mcCollection, m_readKey, ctx));
177 
178 
179  McEventCollection::const_iterator currentGenEventIter = mcCollection->begin();
180  if (currentGenEventIter != mcCollection->end()) {
181  int nvtx = 0;
182  int nvtx_sec=0;
183 #ifdef HEPMC3
184  const auto &barcodes = (*currentGenEventIter)->attribute<HepMC::GenEventBarcodes> ("barcodes");
185  std::map<int,int> id_to_barcode_map;
186  if (barcodes) id_to_barcode_map = barcodes->id_to_barcode_map();
187  for (const auto& vtx: (*currentGenEventIter)->vertices()) {
188  int bcode = id_to_barcode_map[vtx->id()];
189 #else
190  for (HepMC::GenEvent::vertex_const_iterator vtxit=(*currentGenEventIter)->vertices_begin(); vtxit!=(*currentGenEventIter)->vertices_end(); ++vtxit) {
191  auto vtx=*vtxit;
192  int bcode = HepMC::barcode(*vtx);
193 #endif
194  double x = vtx->position().x();
195  double y = vtx->position().y();
196  double z = vtx->position().z();
197  double r = std::sqrt(x*x+y*y);
198  m_h_vtx_x->Fill(x);
199  m_h_vtx_y->Fill(y);
200  m_h_vtx_r->Fill(r);
201  m_h_vtx_z->Fill(z);
202 
203  m_vtx_x->push_back(x);
204  m_vtx_y->push_back(y);
205  m_vtx_r->push_back(r);
206  m_vtx_z->push_back(z);
207  m_vtx_barcode->push_back(bcode);
208 
209  if (!HepMC::is_simulation_vertex(vtx)) {
210  m_h_vtx_prim_xy->Fill(x,y);
211  m_h_vtx_prim_zr->Fill(z,r);
212  ++nvtx;
213  }
214  else {
215  m_h_vtx_sec_xy->Fill(x,y);
216  m_h_vtx_sec_zr->Fill(z,r);
217  ++nvtx_sec;
218  }
219  } //End iteration over vertices
220 
221  m_h_n_vert->Fill(nvtx+nvtx_sec);
222  m_h_n_vert_prim->Fill(nvtx);
223  m_h_n_vert_sec->Fill(nvtx_sec);
224 
225  int npart_prim=0;
226  int npart_sec=0;
227 
228  for (auto currentGenParticle: *(*currentGenEventIter)) {
229  const HepMC::FourVector mom = currentGenParticle->momentum();
230 #ifdef HEPMC3
231  int currentGenParticlebarcode = id_to_barcode_map[currentGenParticle->id()];
232 #else
233  int currentGenParticlebarcode=HepMC::barcode(currentGenParticle);
234 #endif
235  m_h_truth_px->Fill(mom.x());
236  m_h_truth_py->Fill(mom.y());
237  m_h_truth_pz->Fill(mom.z());
238  m_h_truth_pt->Fill(mom.perp());
239  m_h_truth_eta->Fill(mom.eta());
240  m_h_truth_phi->Fill(mom.phi());
241  m_h_barcode->Fill(currentGenParticlebarcode);
242  m_h_part_status->Fill(currentGenParticle->status());
243  m_truth_px->push_back(mom.x());
244  m_truth_py->push_back(mom.y());
245  m_truth_pz->push_back(mom.z());
246  m_truth_pt->push_back(mom.perp());
247  m_truth_eta->push_back(mom.eta());
248  m_truth_phi->push_back(mom.phi());
249  m_barcode->push_back(currentGenParticlebarcode);
250  m_status->push_back(currentGenParticle->status());
251 
252  int pdg = currentGenParticle->pdg_id();
253  m_pdgid->push_back(pdg);
254 
255  if (!HepMC::is_simulation_particle(currentGenParticle)) {
256  m_h_part_pdgid->Fill(pdg);
257  m_h_part_p->Fill(std::sqrt(mom.x()*mom.x()+mom.y()*mom.y()+mom.z()*mom.z()));
258  m_h_part_eta->Fill(mom.eta());
259  m_h_part_phi->Fill(mom.phi());
260  ++npart_prim;
261  m_h_n_generations->Fill(1);
262 /* The conditions don't make any sense.
263  if (currentGenParticlebarcode < HepMC::PHOTOSMIN) {
264  m_h_n_generations->Fill(0);
265  }
266  else {
267  m_h_n_generations->Fill(1);
268  }
269 */
270  } //End simulation particle
271  else {
272  m_h_part_pdgid_sec->Fill(pdg);
273  ++npart_sec;
274  const int gen = HepMC::generations(currentGenParticle) + 2;
275  m_h_n_generations->Fill(gen);
276  }
277  } // End iteration over particles
278 
279  m_h_n_part_prim->Fill(npart_prim);
280  m_h_n_part_sec->Fill(npart_sec);
281  m_h_n_part->Fill(npart_prim+npart_sec);
282  } // End mcCollection
283 
284  m_tree->Fill();
285 
286  return StatusCode::SUCCESS;
287 }
TruthHitAnalysis::m_h_truth_eta
TH1 * m_h_truth_eta
Definition: TruthHitAnalysis.h:45
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
beamspotman.r
def r
Definition: beamspotman.py:672
TruthHitAnalysis::m_h_n_vert_prim
TH1 * m_h_n_vert_prim
Definition: TruthHitAnalysis.h:28
TruthHitAnalysis::m_h_vtx_prim_zr
TH2 * m_h_vtx_prim_zr
Definition: TruthHitAnalysis.h:37
TruthHitAnalysis::m_h_truth_py
TH1 * m_h_truth_py
Definition: TruthHitAnalysis.h:42
TruthHitAnalysis::m_h_truth_pz
TH1 * m_h_truth_pz
Definition: TruthHitAnalysis.h:43
TruthHitAnalysis::m_pdgid
std::vector< float > * m_pdgid
Definition: TruthHitAnalysis.h:68
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
GenEvent.h
TruthHitAnalysis::m_h_part_status
TH1 * m_h_part_status
Definition: TruthHitAnalysis.h:48
TruthHitAnalysis::m_h_part_p
TH1 * m_h_part_p
Definition: TruthHitAnalysis.h:53
TruthHitAnalysis::m_h_barcode
TH1 * m_h_barcode
Definition: TruthHitAnalysis.h:47
TruthHitAnalysis::m_h_truth_pt
TH1 * m_h_truth_pt
Definition: TruthHitAnalysis.h:44
TruthHitAnalysis::m_truth_py
std::vector< float > * m_truth_py
Definition: TruthHitAnalysis.h:61
TruthHitAnalysis::m_h_truth_px
TH1 * m_h_truth_px
Definition: TruthHitAnalysis.h:41
TruthHitAnalysis::m_barcode
std::vector< float > * m_barcode
Definition: TruthHitAnalysis.h:66
TruthHitAnalysis::m_h_part_pdgid
TH1 * m_h_part_pdgid
Definition: TruthHitAnalysis.h:49
TruthHitAnalysis::m_tree
TTree * m_tree
Definition: TruthHitAnalysis.h:70
TruthHitAnalysis::m_h_vtx_sec_zr
TH2 * m_h_vtx_sec_zr
Definition: TruthHitAnalysis.h:39
TruthHitAnalysis::m_h_vtx_prim_xy
TH2 * m_h_vtx_prim_xy
Definition: TruthHitAnalysis.h:36
TruthHitAnalysis::m_h_n_part_prim
TH1 * m_h_n_part_prim
Definition: TruthHitAnalysis.h:29
TruthHitAnalysis::m_h_n_vert
TH1 * m_h_n_vert
Some variables.
Definition: TruthHitAnalysis.h:26
TruthHitAnalysis::m_h_n_part_sec
TH1 * m_h_n_part_sec
Definition: TruthHitAnalysis.h:31
TruthHitAnalysis::m_path
Gaudi::Property< std::string > m_path
Definition: TruthHitAnalysis.h:72
TruthHitAnalysis::m_truth_phi
std::vector< float > * m_truth_phi
Definition: TruthHitAnalysis.h:65
TruthHitAnalysis::execute
virtual StatusCode execute() override final
Definition: TruthHitAnalysis.cxx:156
x
#define x
TruthHitAnalysis::m_h_vtx_x
TH1 * m_h_vtx_x
Definition: TruthHitAnalysis.h:32
TruthHitAnalysis::m_h_vtx_r
TH1 * m_h_vtx_r
Definition: TruthHitAnalysis.h:35
TruthHitAnalysis::m_vtx_r
std::vector< float > * m_vtx_r
Definition: TruthHitAnalysis.h:58
master.gen
gen
Definition: master.py:32
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
TruthHitAnalysis.h
TruthHitAnalysis::m_vtx_barcode
std::vector< float > * m_vtx_barcode
Definition: TruthHitAnalysis.h:59
HepMC::is_simulation_particle
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...
Definition: MagicNumbers.h:354
python.Dumpers.barcodes
def barcodes(beg, end, sz)
Definition: Dumpers.py:2831
z
#define z
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
TruthHitAnalysis::m_h_part_pdgid_sec
TH1 * m_h_part_pdgid_sec
Definition: TruthHitAnalysis.h:50
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
HepMC::barcode
int barcode(const T *p)
Definition: Barcode.h:16
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
HepMC::is_simulation_vertex
bool is_simulation_vertex(const T &v)
Method to establish if the vertex was created during simulation (TODO migrate to be based on status).
Definition: MagicNumbers.h:360
TruthHitAnalysis::m_readKey
SG::ReadHandleKey< McEventCollection > m_readKey
Definition: TruthHitAnalysis.h:74
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
McEventCollection
This defines the McEventCollection, which is really just an ObjectVector of McEvent objects.
Definition: McEventCollection.h:32
TruthHitAnalysis::m_h_n_vert_sec
TH1 * m_h_n_vert_sec
Definition: TruthHitAnalysis.h:30
TruthHitAnalysis::m_h_n_part
TH1 * m_h_n_part
Definition: TruthHitAnalysis.h:27
MagicNumbers.h
TruthHitAnalysis::m_h_vtx_z
TH1 * m_h_vtx_z
Definition: TruthHitAnalysis.h:34
TruthHitAnalysis::m_status
std::vector< float > * m_status
Definition: TruthHitAnalysis.h:67
TruthHitAnalysis::m_truth_px
std::vector< float > * m_truth_px
Definition: TruthHitAnalysis.h:60
TruthHitAnalysis::m_h_vtx_y
TH1 * m_h_vtx_y
Definition: TruthHitAnalysis.h:33
TruthHitAnalysis::m_h_vtx_sec_xy
TH2 * m_h_vtx_sec_xy
Definition: TruthHitAnalysis.h:38
y
#define y
TruthHitAnalysis::m_vtx_y
std::vector< float > * m_vtx_y
Definition: TruthHitAnalysis.h:56
TruthHitAnalysis::m_h_truth_phi
TH1 * m_h_truth_phi
Definition: TruthHitAnalysis.h:46
TruthHitAnalysis::m_vtx_z
std::vector< float > * m_vtx_z
Definition: TruthHitAnalysis.h:57
TruthHitAnalysis::m_h_part_eta
TH1 * m_h_part_eta
Definition: TruthHitAnalysis.h:51
HepMC::generations
int generations(const T &p)
Method to return how many interactions a particle has undergone during simulation (TODO migrate to be...
Definition: MagicNumbers.h:357
TruthHitAnalysis::m_h_n_generations
TH1 * m_h_n_generations
Definition: TruthHitAnalysis.h:40
TruthHitAnalysis::m_h_part_phi
TH1 * m_h_part_phi
Definition: TruthHitAnalysis.h:52
TruthHitAnalysis::m_truth_pt
std::vector< float > * m_truth_pt
Definition: TruthHitAnalysis.h:63
TruthHitAnalysis::m_truth_eta
std::vector< float > * m_truth_eta
Definition: TruthHitAnalysis.h:64
TruthHitAnalysis::m_ntupleFileName
Gaudi::Property< std::string > m_ntupleFileName
Definition: TruthHitAnalysis.h:73
TruthHitAnalysis::initialize
virtual StatusCode initialize() override final
Definition: TruthHitAnalysis.cxx:13
TruthHitAnalysis::m_vtx_x
std::vector< float > * m_vtx_x
Definition: TruthHitAnalysis.h:55
TruthHitAnalysis::m_truth_pz
std::vector< float > * m_truth_pz
Definition: TruthHitAnalysis.h:62
Generated on Mon Sep 29 2025 21:23:00 for ATLAS Offline Software by   doxygen 1.8.18