ATLAS Offline Software
TruthTestTool.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 "TruthTestTool.h"
6 
7 #include "AtlasHepMC/GenEvent.h"
8 #include "TruthUtils/MagicNumbers.h"
9 
10 #include "GeneratorObjects/McEventCollection.h"
11 #include <TH1.h>
12 #include <TH2D.h>
13 
14 #include <fstream>
15 
16 
17 TruthTestTool::TruthTestTool(const std::string& type, const std::string& name, const IInterface* parent)
18  : SimTestToolBase(type, name, parent),
19  m_n_vert(0), m_n_part(0),
20  m_n_vert_prim(0), m_n_vert_sec(0),
21  m_n_part_prim(0), m_n_part_sec(0),
22  m_x_vert(0), m_y_vert(0), m_z_vert(0),
23  m_n_generations(0),
24  m_vtx_r(0), m_vtx_z(0),
25  m_vtx_zr_indet(0), m_vtx_xy_indet(0),
26  m_meanx_vert(0), m_meany_vert(0), m_meanz_vert(0),
27  m_px_truth(0), m_py_truth(0), m_pz_truth(0), m_pt_truth(0), m_log_pt_truth(0),
28  m_theta_truth(0), m_eta_truth(0), m_phi_truth(0),
29  m_origin(0),
30  m_particle_status(0),m_particle_type(0),
31  m_p_gen(0), m_log_p_gen(0), m_eta_gen(0), m_phi_gen(0),m_pion_mass(0)
32 {
33  declareProperty("DumpMcEvent", m_mcEventDump=10);
34 }
35 
36 
37 StatusCode TruthTestTool::initialize()
38 {
39  ATH_CHECK(m_eventInfoKey.initialize());
40 
41  // initialize standard TruthTestTool
42  std::string histdirname("Truth");
43  if (m_key!="TruthEvent") {
44  // initialize alternative TruthTestTool e.g. for GEN_EVENT
45  histdirname = m_key;
46  }
47  m_path+=histdirname+"/";
48 
49  _TH1D(m_n_vert,"N_vert_all",50,0.,1000.);
50  _SET_TITLE(m_n_vert,"total number of verticies","n","dN/dn");
51  m_n_vert -> StatOverflows();
52  _TH1D(m_n_part,"N_part_all",50,0.,5000.);
53  _SET_TITLE(m_n_part,"total number of particles","n","dN/dn");
54  m_n_part -> StatOverflows();
55 
56  _TH1D(m_n_vert_prim,"N_vert_prim",25,0.,750.);
57  _SET_TITLE(m_n_vert_prim,"number of primary verticies","n","dN/dn");
58  m_n_vert_prim -> StatOverflows();
59  _TH1D(m_n_part_prim,"N_part_prim",25,0.,1500.);
60  _SET_TITLE(m_n_part_prim,"number of primary particles","n","dN/dn");
61  m_n_part_prim -> StatOverflows();
62  _TH1D(m_n_vert_sec,"N_vert_sec",25,0.,750.);
63  _SET_TITLE(m_n_vert_sec,"number of secondary verticies","n","dN/dn");
64  m_n_vert_sec -> StatOverflows();
65  _TH1D(m_n_part_sec,"N_part_sec",25,0.,1500.);
66  _SET_TITLE(m_n_part_sec,"number of secondary particles","n","dN/dn");
67  m_n_part_sec -> StatOverflows();
68 
69  _TH1D(m_x_vert,"X_vert_gen",100,-0.5,0.5);
70  _SET_TITLE(m_x_vert,"x vertex distribution","x [mm]","dN/dx [1/mm]");
71  _TH1D(m_y_vert,"Y_vert_gen",100,0.5,1.5);
72  _SET_TITLE(m_y_vert,"y vertex distribution","y [mm]","dN/dy [1/mm]");
73  _TH1D(m_z_vert,"Z_vert_gen",100,-300.,300.);
74  _SET_TITLE(m_z_vert,"z vertex distribution","z [mm]","dN/dz [1/mm]");
75 
76  _TH1D(m_n_generations,"part_generation",26,-0.5,25.5);
77  _SET_TITLE(m_n_generations,"number of particles per generation (prim,sec,...)","n","dN/dn");
78  m_n_generations -> StatOverflows();
79 
80  _TH1D(m_vtx_r,"vertex_r",150,0.,4500.);
81  _SET_TITLE(m_vtx_r, "GenVertex r positions","r [mm]","dN/dr");
82  _TH1D(m_vtx_z,"vertex_z",150,-25000.,25000.);
83  _SET_TITLE(m_vtx_z, "GenVertex z positions","z [mm]", "dN/dz");
84 
85  _TH2D(m_vtx_zr_indet,"vertex_zr_indet",100,-3200.,3200.,100,0.,1200.);
86  _SET_TITLE(m_vtx_zr_indet, "GenVertex zr positions","z [mm]","r [mm]");
87  _TH2D(m_vtx_xy_indet,"vertex_xy_indet",100,-1200.,1200.,100,-1200.,1200.);
88  _SET_TITLE(m_vtx_xy_indet, "GenVertex xy positions","x [mm]","y [mm]");
89 
90 
91  _TH1D(m_meanx_vert,"MeanX_vert_gen",100,-0.5,0.5);
92  _SET_TITLE(m_meanx_vert,"mean x vertex distribution","x [mm]","dN/dx [1/mm]");
93  _TH1D(m_meany_vert,"MeanY_vert_gen",100,0.5,1.5);
94  _SET_TITLE(m_meany_vert,"mean y vertex distribution","y [mm]","dN/dy [1/mm]");
95  _TH1D(m_meanz_vert,"MeanZ_vert_gen",100,-300.,300.);
96  _SET_TITLE(m_meanz_vert,"mean z vertex distribution","z [mm]","dN/dz [1/mm]");
97 
98  // MC truth (incl. G4 particles created in the inner detector
99  _TH1D(m_pt_truth,"Pt_truth",25,0.,50000.);
100  _SET_TITLE(m_pt_truth,"pt distribution of MC truth","pt [MeV]","dN/dpt [1/MeV]");
101  _TH1D(m_log_pt_truth,"log_Pt_truth",100.,0.,16.);
102  _SET_TITLE(m_log_pt_truth,"log(pt) distribution of MC truth","log(pt [MeV])","dN/dlog(pt) [1/log(MeV)]");
103 
104 
105  _TH1D(m_px_truth,"PX_truth",25,-50000.,50000.);
106  _SET_TITLE(m_px_truth,"px distribution of MC truth","px [MeV]","dN/dp [1/MeV]");
107  _TH1D(m_py_truth,"PY_truth",25,-50000.,50000.);
108  _SET_TITLE(m_py_truth,"py distribution of MC truth","py [MeV]","dN/dp [1/MeV]");
109  _TH1D(m_pz_truth,"PZ_truth",25,-50000.,50000.);
110  _SET_TITLE(m_pz_truth,"pz distribution of MC truth","pz [MeV]","dN/dp [1/MeV]");
111 
112  _TH1D(m_theta_truth,"theta_truth",25,0,M_PI);
113  _SET_TITLE(m_theta_truth,"theta distribution of MC truth","theta","dN/dtheta");
114  _TH1D(m_eta_truth,"eta_truth",25,-5.,5.);
115  _SET_TITLE(m_eta_truth,"eta distribution of MC truth","eta","dN/deta");
116  _TH1D(m_phi_truth,"phi_truth",25,-M_PI,M_PI);
117  _SET_TITLE(m_phi_truth,"phi distribution of MC truth","phi","dN/dphi");
118  _TH1D(m_origin,(histdirname+"_origin").c_str(),2,0.,1.);
119  _TH1D(m_particle_status,"particle_status",5,-0.5,4.5);
120  _SET_TITLE(m_particle_status,"Particle Status (1==undecayed, 2==decayed, 4==beam particle)","GenParticle Status","N GenParticles");
121  _TH1D(m_particle_type,"particle_type",15,-7.5,7.5);
122  _SET_TITLE(m_particle_type,"Type (Type 0 = PDG 22, 1=11, 2=13,15, 3=111,211, 4=130,310,311,321, 5=2122, 6=2212, 7=else)","Particle Type (negative for anti-particles)","N Particles");
123 
124 
125  // generator particles
126  _TH1D(m_eta_gen,"eta_gen",25,-5.,5.);
127  _SET_TITLE(m_eta_gen,"eta distribution of generator particles","eta","dN/deta");
128  _TH1D(m_phi_gen,"phi_gen",25,-M_PI,M_PI);
129  _SET_TITLE(m_phi_gen,"phi distribution of generator particles","phi","dN/dphi");
130  _TH1D(m_p_gen,"p_gen",50,0.,250000.);
131  _SET_TITLE(m_p_gen,"momentum distribution of generator particles","p","dN/dp [1/MeV]");
132  _TH1D(m_log_p_gen,"log_p_gen",50,0.,16.);
133  _SET_TITLE(m_log_p_gen,"momentum distribution of generator particles","log(p)","dN/dlog(p) [1/log(MeV)]");
134 
135 
136  _TH1D(m_pion_mass,"p_pion_mass",50,114.,164.);
137  _SET_TITLE(m_pion_mass,"pion mass distribution (generator particles)","m [MeV]","dN/dm [1/MeV]");
138 
139  return StatusCode::SUCCESS;
140 }
141 
142 
143 StatusCode TruthTestTool::processEvent()
144 {
145  SG::ReadHandle<xAOD::EventInfo> eventInfo (m_eventInfoKey,Gaudi::Hive::currentContext());
146  ATH_CHECK(eventInfo.isValid());
147 
148  const int evtnum(eventInfo->eventNumber());
149 
150  const McEventCollection* mcCollection;
151  if (evtStore()->retrieve(mcCollection,m_key).isSuccess()) {
152  McEventCollection::const_iterator currentGenEventIter = mcCollection->begin();
153 
154  if (currentGenEventIter!=mcCollection->end() ) {
155  {
156  auto vtx = HepMC::signal_process_vertex(*currentGenEventIter);
157 #ifdef HEPMC3
158  if (!vtx && (*currentGenEventIter)->vertices().size()>0) vtx=((*currentGenEventIter)->vertices()).front();
159 #else
160  if (!vtx && (*currentGenEventIter)->vertices_size()>0) vtx=*((*currentGenEventIter)->vertices_begin());
161 #endif
162  if ( vtx )
163  {
164  m_x_vert->Fill( vtx->position().x() );
165  m_y_vert->Fill( vtx->position().y() );
166  m_z_vert->Fill( vtx->position().z() );
167  }
168  }
169 
170  int nvtx=0;
171  int nvtx_sec=0;
172  float mx=0.,my=0.,mz=0.;
173 #ifdef HEPMC3
174  for (const auto& vtx: (*currentGenEventIter)->vertices()) {
175 #else
176  for (HepMC::GenEvent::vertex_const_iterator vtxit=(*currentGenEventIter)->vertices_begin();
177  vtxit!=(*currentGenEventIter)->vertices_end();++vtxit) {
178  auto vtx=*vtxit;
179 #endif
180 
181  double x = vtx->position().x();
182  double y = vtx->position().y();
183  double z = vtx->position().z();
184  double r = std::sqrt(x*x+y*y);
185  m_vtx_r->Fill(r);
186  m_vtx_z->Fill(z);
187  m_vtx_zr_indet->Fill(z, r);
188  m_vtx_xy_indet->Fill(x, y);
189  if ( !HepMC::is_simulation_vertex(vtx) ) {
190  mx+=x;
191  my+=y;
192  mz+=z;
193  ++nvtx;
194  }
195  else {
196  // TODO: add 2D distribution of secondary verticies (and primary ?)
197  ++nvtx_sec;
198  }
199  }
200  if (nvtx>0) {
201  m_meanx_vert->Fill( mx/nvtx );
202  m_meany_vert->Fill( my/nvtx );
203  m_meanz_vert->Fill( mz/nvtx );
204  }
205 
206  m_n_vert->Fill( nvtx + nvtx_sec );
207  m_n_vert_prim->Fill( nvtx);
208  m_n_vert_sec->Fill( nvtx_sec );
209 
210  }
211 
212  for (;currentGenEventIter!=mcCollection->end(); ++currentGenEventIter) {
213  if (evtnum<m_mcEventDump) {
214  char fname[80];
215  sprintf(fname,"%s.event%d.txt",m_key.c_str(),evtnum);
216  std::ofstream of(fname);
217  HepMC::Print::line(of,*(*currentGenEventIter)); // verbose output
218  of.close();
219  }
220 
221  int npart_prim=0, npart_sec=0;
222  for (const auto& currentGenParticle: *(*currentGenEventIter)) {
223  bool is_simulation = HepMC::is_simulation_particle(currentGenParticle);
224  const HepMC::FourVector mom = currentGenParticle->momentum();
225  m_px_truth->Fill( mom.x() );
226  m_py_truth->Fill( mom.y() );
227  m_pz_truth->Fill( mom.z() );
228  m_pt_truth->Fill( mom.perp() );
229  m_log_pt_truth->Fill( mom.perp() > 0 ? log(mom.perp()) : -1 );
230 
231  m_theta_truth->Fill( mom.theta() );
232  m_eta_truth->Fill( mom.eta() );
233  m_phi_truth->Fill( mom.phi() );
234 
235  if(std::abs(currentGenParticle->pdg_id())==211) {
236  m_pion_mass->Fill(mom.m());
237  }
238  m_origin->Fill(is_simulation?1:0);
239 
240  m_particle_status->Fill(currentGenParticle->status());
241 
242  int pdg = currentGenParticle->pdg_id();
243  int particleType = 100;
244  switch ( abs(pdg) ) {
245  case 22:
246  particleType = 0;
247  break;
248  case 11:
249  particleType = 1;
250  break;
251  case 13:
252  case 15:
253  particleType = 2;
254  break;
255  case 111:
256  case 211:
257  particleType = 3;
258  break;
259  case 130:
260  case 310:
261  case 311:
262  case 321:
263  particleType = 4;
264  break;
265  case 2112:
266  particleType = 5;
267  break;
268  case 2212:
269  particleType = 6;
270  break;
271  default:
272  particleType = 7;
273  break;
274  }
275  particleType = (pdg<0) ? -particleType : particleType;
276  m_particle_type->Fill( particleType );
277 
278  if ( !is_simulation ) {
279  double momentum=std::sqrt(mom.x()*mom.x()+mom.y()*mom.y()+mom.z()*mom.z());
280  m_p_gen->Fill( momentum );
281  m_log_p_gen->Fill( std::log(momentum) );
282  m_eta_gen->Fill( mom.eta() );
283  m_phi_gen->Fill( mom.phi() );
284  ++npart_prim;
285  }
286  else {
287  ++npart_sec;
288  }
289  m_n_generations ->Fill(HepMC::generations(currentGenParticle));
290  }
291  m_n_part_prim->Fill(npart_prim);
292  m_n_part_sec->Fill(npart_sec);
293  m_n_part->Fill(npart_prim+npart_sec);
294  }
295  }
296 
297  return StatusCode::SUCCESS;
298 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
TruthTestTool::m_theta_truth
TH1 * m_theta_truth
Definition: TruthTestTool.h:37
beamspotman.r
def r
Definition: beamspotman.py:676
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
GenEvent.h
TruthTestTool::TruthTestTool
TruthTestTool(const std::string &type, const std::string &name, const IInterface *parent)
Definition: TruthTestTool.cxx:17
fitman.my
my
Definition: fitman.py:523
xAOD::EventInfo_v1::eventNumber
uint64_t eventNumber() const
The current event's event number.
TruthTestTool.h
TruthTestTool::m_n_part_sec
TH1 * m_n_part_sec
Definition: TruthTestTool.h:30
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
LB_AnalMapSplitter.of
of
Definition: LB_AnalMapSplitter.py:48
TruthTestTool::initialize
virtual StatusCode initialize() override
Definition: TruthTestTool.cxx:37
TruthTestTool::m_phi_truth
TH1 * m_phi_truth
Definition: TruthTestTool.h:37
M_PI
#define M_PI
Definition: ActiveFraction.h:11
fitman.mx
mx
Definition: fitman.py:520
TruthTestTool::m_pz_truth
TH1 * m_pz_truth
Definition: TruthTestTool.h:36
_SET_TITLE
#define _SET_TITLE(var, title, xaxis, yaxis)
Definition: SimTestHisto.h:93
TruthTestTool::m_meanz_vert
TH1 * m_meanz_vert
Definition: TruthTestTool.h:35
HepMC::Print::line
void line(std::ostream &os, const GenEvent &e)
Definition: GenEvent.h:676
TruthTestTool::m_pt_truth
TH1 * m_pt_truth
Definition: TruthTestTool.h:36
TruthTestTool::m_vtx_zr_indet
TH2 * m_vtx_zr_indet
Definition: TruthTestTool.h:34
x
#define x
particleType
Definition: particleType.h:29
TruthTestTool::m_eventInfoKey
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfoKey
SG key for Event Info.
Definition: TruthTestTool.h:25
_TH1D
#define _TH1D(var, name, nbin, xmin, xmax)
Definition: SimTestHisto.h:47
TruthTestTool::m_n_generations
TH1 * m_n_generations
Definition: TruthTestTool.h:32
TruthTestTool::m_pion_mass
TH1 * m_pion_mass
Definition: TruthTestTool.h:40
TruthTestTool::m_meanx_vert
TH1 * m_meanx_vert
Definition: TruthTestTool.h:35
TruthTestTool::m_particle_type
TH1 * m_particle_type
Definition: TruthTestTool.h:39
TruthTestTool::m_log_p_gen
TH1 * m_log_p_gen
Definition: TruthTestTool.h:40
SimTestToolBase
Definition: SimTestToolBase.h:20
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
TruthTestTool::m_vtx_r
TH1 * m_vtx_r
Definition: TruthTestTool.h:33
TruthTestTool::m_p_gen
TH1 * m_p_gen
Definition: TruthTestTool.h:40
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
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:342
McEventCollection.h
z
#define z
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
fitman.mz
mz
Definition: fitman.py:526
TruthTestTool::m_particle_status
TH1 * m_particle_status
Definition: TruthTestTool.h:39
TruthTestTool::m_n_vert_sec
TH1 * m_n_vert_sec
Definition: TruthTestTool.h:29
test_pyathena.parent
parent
Definition: test_pyathena.py:15
TruthTestTool::m_n_vert
TH1 * m_n_vert
Definition: TruthTestTool.h:28
_TH2D
#define _TH2D(var, name, nbinx, xmin, xmax, nbiny, ymin, ymax)
Definition: SimTestHisto.h:81
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:348
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:33
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
TruthTestTool::m_meany_vert
TH1 * m_meany_vert
Definition: TruthTestTool.h:35
TruthTestTool::m_n_part
TH1 * m_n_part
Definition: TruthTestTool.h:28
TruthTestTool::m_y_vert
TH1 * m_y_vert
Definition: TruthTestTool.h:31
TruthTestTool::m_x_vert
TH1 * m_x_vert
Definition: TruthTestTool.h:31
TruthTestTool::m_vtx_z
TH1 * m_vtx_z
Definition: TruthTestTool.h:33
TruthTestTool::m_z_vert
TH1 * m_z_vert
Definition: TruthTestTool.h:31
TruthTestTool::m_phi_gen
TH1 * m_phi_gen
Definition: TruthTestTool.h:40
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
MagicNumbers.h
TruthTestTool::m_eta_truth
TH1 * m_eta_truth
Definition: TruthTestTool.h:37
TruthTestTool::m_px_truth
TH1 * m_px_truth
Definition: TruthTestTool.h:36
TruthTestTool::m_n_part_prim
TH1 * m_n_part_prim
Definition: TruthTestTool.h:30
TruthTestTool::processEvent
virtual StatusCode processEvent() override
Definition: TruthTestTool.cxx:143
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
python.AthDsoLogger.fname
string fname
Definition: AthDsoLogger.py:67
y
#define y
TruthTestTool::m_vtx_xy_indet
TH2 * m_vtx_xy_indet
Definition: TruthTestTool.h:34
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
TruthTestTool::m_log_pt_truth
TH1 * m_log_pt_truth
Definition: TruthTestTool.h:36
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
TruthTestTool::m_mcEventDump
int m_mcEventDump
Definition: TruthTestTool.h:42
SimTestHisto::m_path
std::string m_path
Definition: SimTestHisto.h:34
TruthTestTool::m_n_vert_prim
TH1 * m_n_vert_prim
Definition: TruthTestTool.h:29
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:345
SimTestToolBase::m_key
std::string m_key
The MC truth key.
Definition: SimTestToolBase.h:34
TruthTestTool::m_py_truth
TH1 * m_py_truth
Definition: TruthTestTool.h:36
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
TruthTestTool::m_eta_gen
TH1 * m_eta_gen
Definition: TruthTestTool.h:40
TruthTestTool::m_origin
TH1 * m_origin
Definition: TruthTestTool.h:38
HepMC::signal_process_vertex
GenVertex * signal_process_vertex(const GenEvent *e)
Definition: GenEvent.h:625
Generated on Thu Nov 7 2024 21:31:10 for ATLAS Offline Software by   doxygen 1.8.18