DD4hep - The AIDA detector description toolkit for high energy physics experiments
DD4hep  Rev:Unversioneddirectory
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Geant4Particle.cpp
Go to the documentation of this file.
1 // $Id: $
2 //==========================================================================
3 // AIDA Detector description implementation for LCD
4 //--------------------------------------------------------------------------
5 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
6 // All rights reserved.
7 //
8 // For the licensing terms see $DD4hepINSTALL/LICENSE.
9 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
10 //
11 // Author : M.Frank
12 //
13 //==========================================================================
14 
15 // Framework include files
16 #include "DD4hep/Printout.h"
17 #include "DD4hep/Primitives.h"
18 #include "DD4hep/InstanceCount.h"
19 #include "DDG4/Geant4Particle.h"
20 #include "TDatabasePDG.h"
21 #include "TParticlePDG.h"
22 #include "G4ParticleTable.hh"
23 #include "G4ParticleDefinition.hh"
24 #include "G4VProcess.hh"
25 #include "G4ChargedGeantino.hh"
26 #include "G4Geantino.hh"
27 
28 #include <iostream>
29 
30 using namespace DD4hep;
31 using namespace DD4hep::Simulation;
33 
36 }
37 
40 : ref(1), id(c.id), g4Parent(c.g4Parent), reason(c.reason), mask(c.mask),
41  steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID),
42  status(c.status), charge(0),
43  vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
44  vex(c.vex), vey(c.vey), vez(c.vez),
45  psx(c.psx), psy(c.psy), psz(c.psz),
46  pex(c.pex), pey(c.pey), pez(c.pez),
47  mass(c.mass), time(c.time), properTime(c.properTime),
48  parents(c.parents), daughters(c.daughters), extension(),
49  process(c.process)//, definition(c.definition)
50 {
52  spin[0] = c.spin[0];
53  spin[1] = c.spin[1];
54  spin[2] = c.spin[2];
55  colorFlow[0] = c.colorFlow[0];
56  colorFlow[1] = c.colorFlow[1];
57 }
58 
61 : ref(1), id(0), g4Parent(0), reason(0), mask(0),
62  steps(0), secondaries(0), pdgID(0),
63  status(0), charge(0),
64  vsx(0.0), vsy(0.0), vsz(0.0),
65  vex(0.0), vey(0.0), vez(0.0),
66  psx(0.0), psy(0.0), psz(0.0),
67  pex(0.0), pey(0.0), pez(0.0),
68  mass(0.0), time(0.0), properTime(0.0),
69  daughters(), extension(), process(0)//, definition(0)
70 {
72  spin[0] = spin[1] = spin[2] = 0;
73  colorFlow[0] = colorFlow[1] = 0;
74 }
75 
78 : ref(1), id(part_id), g4Parent(0), reason(0), mask(0),
79  steps(0), secondaries(0), pdgID(0),
80  status(0), charge(0),
81  vsx(0.0), vsy(0.0), vsz(0.0),
82  vex(0.0), vey(0.0), vez(0.0),
83  psx(0.0), psy(0.0), psz(0.0),
84  pex(0.0), pey(0.0), pez(0.0),
85  mass(0.0), time(0.0), properTime(0.0),
86  daughters(), extension(), process(0)//, definition(0)
87 {
89  spin[0] = spin[1] = spin[2] = 0;
90  colorFlow[0] = colorFlow[1] = 0;
91 }
92 
96  //::printf("************ Delete Geant4Particle[%p]: ID:%d pdgID %d ref:%d\n",(void*)this,id,pdgID,ref);
97 }
98 
100  //::printf("************ Release Geant4Particle[%p]: ID:%d pdgID %d ref:%d\n",(void*)this,id,pdgID,ref-1);
101  if ( --ref <= 0 ) {
102  delete this;
103  }
104 }
105 
108  if ( this != &c ) {
109  id = c.id;
110  g4Parent = c.g4Parent;
111  reason = c.reason;
112  mask = c.mask;
113  status = c.status;
114  charge = c.charge;
115  steps = c.steps;
117  pdgID = c.pdgID;
118  vsx = c.vsx;
119  vsy = c.vsy;
120  vsz = c.vsz;
121  vex = c.vex;
122  vey = c.vey;
123  vez = c.vez;
124  psx = c.psx;
125  psy = c.psy;
126  psz = c.psz;
127  pex = c.pex;
128  pey = c.pey;
129  pez = c.pez;
130  mass = c.mass;
131  time = c.time;
133  process = c.process;
134  //definition = c.definition;
135  daughters = c.daughters;
136  parents = c.parents;
137 #if __cplusplus >= 201103L
138  extension.swap(c.extension);
139 #else
140  extension = c.extension;
141 #endif
142  //dd4hep_ptr<ParticleExtension>(c.extension.release());
143  }
144  return *this;
145 }
146 
148 void Geant4Particle::removeDaughter(int id_daughter) {
149  std::set<int>::iterator j = daughters.find(id_daughter);
150  if ( j != daughters.end() ) daughters.erase(j);
151 }
152 
154 const G4ParticleDefinition* Geant4ParticleHandle::definition() const {
155  G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
156  G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
157  if ( 0 == def && 0 == particle->pdgID ) {
158  if ( 0 == particle->charge )
159  return G4Geantino::Definition();
160  return G4ChargedGeantino::Definition();
161  }
162  return def;
163 }
164 
167  const G4ParticleDefinition* def = definition();
168  if ( def ) {
169  //particle->definition = def;
170  return def->GetParticleName();
171  }
172 #if 0
173  TDatabasePDG* db = TDatabasePDG::Instance();
174  TParticlePDG* pdef = db->GetParticle(particle->pdgID);
175  if ( pdef ) return pdef->GetName();
176 #endif
177  char text[32];
178  ::snprintf(text,sizeof(text),"PDG:%d",particle->pdgID);
179  return text;
180 }
181 
184  const G4ParticleDefinition* def = definition();
185  if ( def ) {
186  //particle->definition = def;
187  return def->GetParticleType();
188  }
189 #if 0
190  TDatabasePDG* db = TDatabasePDG::Instance();
191  TParticlePDG* pdef = db->GetParticle(particle->pdgID);
192  if ( pdef ) return pdef->ParticleClass();
193 #endif
194  char text[32];
195  ::snprintf(text,sizeof(text),"PDG:%d",particle->pdgID);
196  return text;
197 }
198 
201  if ( particle->process ) return particle->process->GetProcessName();
202  else if ( particle->reason&G4PARTICLE_PRIMARY ) return "Primary";
203  else if ( particle->status&G4PARTICLE_GEN_EMPTY ) return "Gen.Empty";
204  else if ( particle->status&G4PARTICLE_GEN_STABLE ) return "Gen.Stable";
205  else if ( particle->status&G4PARTICLE_GEN_DECAYED ) return "Gen.Decay";
206  else if ( particle->status&G4PARTICLE_GEN_DOCUMENTATION ) return "Gen.DOC";
207  return "???";
208 }
209 
212  if ( particle->process ) {
213  return G4VProcess::GetProcessTypeName(particle->process->GetProcessType());
214  }
215  return "";
216 }
217 
219 void Geant4ParticleHandle::offset(int off) const {
220  std::set<int> temp;
222  p->id += off;
223 
224  temp = p->daughters;
225  p->daughters.clear();
226  for(std::set<int>::iterator i=temp.begin(); i != temp.end(); ++i)
227  p->daughters.insert((*i)+off);
228 
229  temp = p->parents;
230  p->parents.clear();
231  for(std::set<int>::iterator i=temp.begin(); i != temp.end(); ++i)
232  p->parents.insert((*i)+off);
233 }
234 
236 void Geant4ParticleHandle::dump1(int level, const std::string& src, const char* tag) const {
237  char text[256];
238  Geant4ParticleHandle p(*this);
239  text[0]=0;
240  if ( p->parents.size() == 1 )
241  ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
242  else if ( p->parents.size() > 1 ) {
243  text[0]='/';text[1]=0;
244  for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
245  ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
246  }
247  printout((DD4hep::PrintLevel)level,src,
248  "+++ %s %4d def [%-11s,%8s] reason:%8d E:%+.2e %3s #Dau:%3d #Par:%3d%-5s",
249  tag, p->id,
250  p.particleName().c_str(),
251  p.particleType().c_str(),
252  p->reason,
253  p.energy(),
254  p->g4Parent>0 ? "Sim" : "Gen",
255  int(p->daughters.size()),
256  int(p->parents.size()),text);
257 }
258 
260 void Geant4ParticleHandle::dump2(int level, const std::string& src, const char* tag, int g4id, bool inrec) const {
261  char text[32];
262  Geant4ParticleHandle p(*this);
263  if ( p->parents.size() == 0 ) text[0]=0;
264  else if ( p->parents.size() == 1 ) ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
265  else if ( p->parents.size() > 1 ) ::snprintf(text,sizeof(text),"/%d..",*(p->parents.begin()));
266  printout((DD4hep::PrintLevel)level,src,
267  "+++ %s %4d G4:%4d [%-12s,%8s] reason:%8d "
268  "E:%+.2e in record:%s #Par:%3d%-5s #Dau:%3d",
269  tag, p->id, g4id,
270  p.particleName().c_str(),
271  p.particleType().c_str(),
272  p->reason,
273  p.energy(),
274  yes_no(inrec),
275  int(p->parents.size()),text,
276  int(p->daughters.size()));
277 }
278 
280 void Geant4ParticleHandle::dumpWithVertex(int level, const std::string& src, const char* tag) const {
281  char text[256];
282  Geant4ParticleHandle p(*this);
283  text[0]=0;
284  if ( p->parents.size() == 1 )
285  ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
286  else if ( p->parents.size() > 1 ) {
287  text[0]='/';text[1]=0;
288  for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
289  ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
290  }
291  printout((DD4hep::PrintLevel)level,src,
292  "+++ %s ID:%3d %-12s status:%08X PDG:%6d Vtx:(%+.2e,%+.2e,%+.2e)[mm] "
293  "time: %+.2e [ns] #Dau:%3d #Par:%1d%-6s",
294  tag,p->id,p.particleName().c_str(),p->status,p->pdgID,
296  p->daughters.size(),
297  p->parents.size(),
298  text);
299 }
300 
301 
303 void Geant4ParticleHandle::dumpWithMomentum(int level, const std::string& src, const char* tag) const {
304  char text[256];
305  Geant4ParticleHandle p(*this);
306  text[0]=0;
307  if ( p->parents.size() == 1 )
308  ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
309  else if ( p->parents.size() > 1 ) {
310  text[0]='/';text[1]=0;
311  for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
312  ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
313  }
314  printout((DD4hep::PrintLevel)level,src,
315  "+++%s ID:%3d %-12s stat:%08X PDG:%6d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
316  "time: %+.2e [ns] #Dau:%3d #Par:%1d%-6s",
317  tag,p->id,p.particleName().c_str(),p->status,p->pdgID,
319  int(p->daughters.size()),
320  int(p->parents.size()),
321  text);
322 }
323 
325 void Geant4ParticleHandle::dumpWithMomentumAndVertex(int level, const std::string& src, const char* tag) const {
326  char text[256];
327  Geant4ParticleHandle p(*this);
328  text[0]=0;
329  if ( p->parents.size() == 1 )
330  ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
331  else if ( p->parents.size() > 1 ) {
332  text[0]='/';text[1]=0;
333  for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
334  ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
335  }
336  printout((DD4hep::PrintLevel)level,src,
337  "+++%s %3d %-12s stat:%08X PDG:%6d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
338  "Vtx:(%+.2e,%+.2e,%+.2e)[mm] #Dau:%3d #Par:%1d%-6s",
339  tag,p->id,p.particleName().c_str(),p->status,p->pdgID,
342  int(p->daughters.size()),
343  int(p->parents.size()),
344  text);
345 }
346 
347 void Geant4ParticleHandle::dump4(int level, const std::string& src, const char* tag) const {
348  Geant4ParticleHandle p(*this);
349  char equiv[32];
350  PropertyMask mask(p->reason);
351  PropertyMask status(p->status);
352  std::string proc_name = p.processName();
353  std::string proc_type = p.processTypeName();
354  int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
355 
356  equiv[0] = 0;
357  if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
358  ::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
359  }
360  printout((DD4hep::PrintLevel)level,src,
361  "+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s] %c%c%c%c -- %c%c%c%c%c%c%c",
362  tag,
363  p->id,
364  p.particleName().c_str(),
365  parent_id,equiv,
368  int(p->daughters.size()),
370  p.energy(),
374  mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
375  p.numParent(),
376  proc_name.c_str(),
377  p->process ? "/" : "",
378  proc_type.c_str(),
379  status.isSet(G4PARTICLE_GEN_EMPTY) ? 'E' : '.',
380  status.isSet(G4PARTICLE_GEN_STABLE) ? 'S' : '.',
381  status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
382  status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.',
383 
384  status.isSet(G4PARTICLE_SIM_CREATED) ? 's' : '.',
385  status.isSet(G4PARTICLE_SIM_BACKSCATTER) ? 'b' : '.',
386  status.isSet(G4PARTICLE_SIM_PARENT_RADIATED) ? 'v' : '.',
387  status.isSet(G4PARTICLE_SIM_DECAY_TRACKER) ? 't' : '.',
388  status.isSet(G4PARTICLE_SIM_DECAY_CALO) ? 'c' : '.',
389  status.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) ? 'l' : '.',
390  status.isSet(G4PARTICLE_SIM_STOPPED) ? 's' : '.'
391  );
392 }
393 
396  clear();
397 }
398 
402  particleMap.clear();
403  equivalentTracks.clear();
404 }
405 
408  int cnt;
409  char text[64];
410  using namespace std;
411  const Geant4ParticleMap* m = this;
412 
413  cnt = 0;
414  cout << "Particle map:" << endl;
415  for(Geant4ParticleMap::ParticleMap::const_iterator i=m->particleMap.begin(); i!=m->particleMap.end();++i) {
416  ::snprintf(text,sizeof(text)," [%-4d:%p]",(*i).second->id,(void*)(*i).second);
417  cout << text;
418  if ( ++cnt == 8 ) {
419  cout << endl;
420  cnt = 0;
421  }
422  }
423  cout << endl;
424 
425  cnt = 0;
426  cout << "Equivalents:" << endl;
427  for(Geant4ParticleMap::TrackEquivalents::const_iterator i=m->equivalentTracks.begin(); i!=m->equivalentTracks.end();++i) {
428  ::snprintf(text,sizeof(text)," [%-5d : %-5d]",(*i).first,(*i).second);
429  cout << text;
430  if ( ++cnt == 8 ) {
431  cout << endl;
432  cnt = 0;
433  }
434  }
435  cout << endl;
436 }
437 
440  clear();
441  particleMap = pm;
442  equivalentTracks = equiv;
443  pm.clear();
444  equiv.clear();
445  //dump();
446 }
447 
450  return !equivalentTracks.empty();
451 }
452 
454 int Geant4ParticleMap::particleID(int g4_id, bool) const {
455  TrackEquivalents::const_iterator iequiv = equivalentTracks.find(g4_id);
456  if ( iequiv != equivalentTracks.end() ) return (*iequiv).second;
457  printout(ERROR,"Geant4ParticleMap","+++ No Equivalent particle for track:%d."
458  " Monte Carlo truth record looks broken!",g4_id);
459  dump();
460  return -1;
461 }
const G4VProcess * process
not persisten. ROOT cannot handle
static const double MeV
Definition: DD4hepUnits.h:146
dd4hep_ptr< ParticleExtension > extension
User data extension if required.
std::map< int, Particle * > ParticleMap
void offset(int off) const
Handlers.
Data structure to store the MC particle information.
void adopt(ParticleMap &pm, TrackEquivalents &equiv)
Adopt particle maps.
void removeDaughter(int id_daughter)
Remove daughter from set.
void dumpWithVertex(int level, const std::string &src, const char *tag) const
Output type 3:+++ "tag" ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] #Par: 0 #Dau: 4.
static void decrement(T *)
Decrement count according to type information.
const char * yes_no(bool value)
Helper function to print booleans in format YES/NO.
Definition: Printout.h:295
int particleID(int track, bool throw_if_not_found=true) const
Access the equivalent track id (shortcut to the usage of TrackEquivalents)
void dump4(int level, const std::string &src, const char *tag) const
void dump2(int level, const std::string &src, const char *tag, int g4id, bool inrec) const
Output type 2:+++ "tag" 20 G4: 7 def:0xde4eaa8 [gamma , gamma] reason: 20 E:+3.304035e+01 in record:Y...
Data structure to manipulate a bitmask held by reference and represented by an integer.
Definition: Primitives.h:526
PrintLevel
Definition: Printout.h:47
Particles parents
The list of daughters of this MC particle.
Data structure to map particles produced during the generation and the simulation.
const G4ParticleDefinition * definition() const
Access the Geant4 particle definition object (expensive!)
virtual ~ParticleExtension()
Default destructor.
std::string particleName() const
Access to the Geant4 particle name.
void clear()
Clear particle maps.
void releaseObjects(M &m)
Definition: Primitives.h:344
std::string particleType() const
Access to the Geant4 particle type.
size_t numParent() const
Accessor to the number of particle parents.
Data structure to access derived MC particle information.
static void increment(T *)
Increment count according to type information.
Definition: InstanceCount.h:98
double energy() const
Scalar particle energy.
void dumpWithMomentumAndVertex(int level, const std::string &src, const char *tag) const
Output type 3:+++ <tag> ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] #Par: 0 #Dau: 4.
virtual ~Geant4ParticleMap()
Default destructor.
Geant4Particle & get_data(Geant4Particle &c)
Assignment operator.
bool isValid() const
Check if the particle map was ever filled (ie. some particle handler was present) ...
virtual ~Geant4Particle()
Default destructor.
void dump1(int level, const std::string &src, const char *tag) const
Various output formats:
std::string processName() const
Access to the creator process name.
ReferenceBitMask< int > PropertyMask
static const double mm
Definition: DD4hepUnits.h:69
int printout(PrintLevel severity, const char *src, const char *fmt,...)
Calls the display action with a given severity level.
Definition: Printout.cpp:111
Geant4Particle * particle
Particle pointer.
ParticleMap particleMap
Mapping of particles of this event.
bool isSet(const T &m) const
Definition: Primitives.h:544
static const double ns
Definition: DD4hepUnits.h:122
TGeoShape TGeoMedium * m
Definition: Volumes.cpp:294
TrackEquivalents equivalentTracks
Map associating the G4Track identifiers with identifiers of existing MCParticles. ...
void dump() const
Dump content.
Geant4Particle()
Default constructor.
void dumpWithMomentum(int level, const std::string &src, const char *tag) const
Output type 3:+++ <tag> ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] #Par: 0 #Dau: 4.
std::string processTypeName() const
Access to the creator process type name.
void release()
Decrease reference count. Deletes object if NULL.