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
Geant4InputHandling.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
17 #include "DDG4/Geant4Primary.h"
18 #include "DDG4/Geant4Context.h"
19 #include "DDG4/Geant4Action.h"
20 #include "CLHEP/Units/SystemOfUnits.h"
21 #include "CLHEP/Units/PhysicalConstants.h"
22 
23 // Geant4 include files
24 #include "G4ParticleDefinition.hh"
25 #include "G4Event.hh"
26 #include "G4PrimaryVertex.hh"
27 #include "G4PrimaryParticle.hh"
28 
29 // C/C++ include files
30 #include <stdexcept>
31 #include <cmath>
32 
33 using namespace std;
34 using namespace DD4hep;
35 using namespace DD4hep::Simulation;
36 
38 
39 
41 Geant4Vertex* DD4hep::Simulation::createPrimary(const G4PrimaryVertex* g4) {
42  Geant4Vertex* v = new Geant4Vertex();
43  v->x = g4->GetX0();
44  v->y = g4->GetY0();
45  v->z = g4->GetZ0();
46  v->time = g4->GetT0();
47  return v;
48 }
49 
53  const Geant4Vertex* v,
54  const G4PrimaryParticle* g4p)
55 {
57  p->id = particle_id;
58  p->reason = 0;
59  p->pdgID = g4p->GetPDGcode();
60  p->psx = g4p->GetPx();
61  p->psy = g4p->GetPy();
62  p->psz = g4p->GetPz();
63  p->time = g4p->GetProperTime();
64  p->properTime = g4p->GetProperTime();
65  p->vsx = v->x;
66  p->vsy = v->y;
67  p->vsz = v->z;
68  p->vex = v->x;
69  p->vey = v->y;
70  p->vez = v->z;
71  //p->definition = g4p->GetG4code();
72  //p->process = 0;
73  p->spin[0] = 0;
74  p->spin[1] = 0;
75  p->spin[2] = 0;
76  p->colorFlow[0] = 0;
77  p->colorFlow[0] = 0;
78  p->mass = g4p->GetMass();
79  p->charge = g4p->GetCharge();
80  PropertyMask status(p->status);
81  status.set(G4PARTICLE_GEN_STABLE);
82  return p;
83 }
84 
87  Geant4PrimaryInteraction* interaction,
88  Geant4Vertex* particle_origine,
89  G4PrimaryParticle* gp)
90 {
91  int pid = int(interaction->particles.size());
92  Geant4Particle* p = createPrimary(pid,particle_origine,gp);
93  G4PrimaryParticle* dau = gp->GetDaughter();
94  PropertyMask status(p->status);
95  int mask = interaction->mask;
96 
97  interaction->particles.insert(make_pair(p->id,p));
98  status.set(G4PARTICLE_PRIMARY);
99  p->mask = mask;
100  particle_origine->out.insert(p->id);
101  // Insert pair in map. Do not forget to increase reference count!
102  pm->insert(gp,p);
103 
104  if ( dau ) {
105  Geant4Vertex* dv = new Geant4Vertex(*particle_origine);
106  int vid = int(interaction->vertices.size());
107  PropertyMask reason(p->reason);
108  reason.set(G4PARTICLE_HAS_SECONDARIES);
109 
110  dv->mask = mask;
111  dv->in.insert(p->id);
112  interaction->vertices.insert(make_pair(vid,dv));
113  for(; dau; dau = dau->GetNext())
114  collectPrimaries(pm, interaction, dv, dau);
115  }
116 }
117 
121  Geant4PrimaryMap* pm,
122  const G4PrimaryVertex* gv)
123 {
126  int vid = int(interaction->vertices.size());
127  interaction->locked = true;
128  interaction->mask = mask;
129  v->mask = mask;
130  interaction->vertices.insert(make_pair(vid,v));
131  for (G4PrimaryParticle *gp = gv->GetPrimary(); gp; gp = gp->GetNext() )
132  collectPrimaries(pm, interaction, v, gp);
133  return interaction;
134 }
135 
138  const Geant4Context* context)
139 {
145  context->event().addExtension(new Geant4PrimaryMap());
146 
147  // The final set of created particles in the simulation.
148  context->event().addExtension(new Geant4ParticleMap());
149 
150  //
151  // The Geant4PrimaryEvent extension contains a whole set of
152  // Geant4PrimaryInteraction objects each may represent a complete
153  // interaction. Particles and vertices may be unbiased.
154  // This is the input to the translator forming the final
155  // Geant4PrimaryInteraction (see below) containing rebiased particle
156  // and vertex maps.
158  context->event().addExtension(evt);
159  //
160  // The Geant4PrimaryInteraction extension contains the final
161  // vertices and particles ready to be injected to Geant4.
163  inter->setNextPID(-1);
164  context->event().addExtension(inter);
165  return 1;
166 }
167 
169 static void appendInteraction(const Geant4Action* caller,
170  Geant4PrimaryInteraction* output,
172 {
173  Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
174  for( ip=input->particles.begin(), ipend=input->particles.end(); ip != ipend; ++ip ) {
175  Geant4Particle* p = (*ip).second;
176  output->particles.insert(make_pair(p->id,p->addRef()));
177  }
178  Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
179  for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv ) {
180  ivfnd = output->vertices.find((*iv).second->mask);
181  if ( ivfnd != output->vertices.end() ) {
182  caller->abortRun("Duplicate primary interaction identifier!",
183  "Cannot handle 2 interactions with identical identifiers!");
184  }
185  output->vertices.insert(make_pair((*iv).second->mask,(*iv).second->addRef()));
186  }
187 }
188 
189 static void rebaseParticles(Geant4PrimaryInteraction::ParticleMap& particles, int &offset) {
190  Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
191  int mx_id = offset;
192  // Now move begin and end-vertex of all primary and generator particles accordingly
193  for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip ) {
194  Geant4ParticleHandle p((*ip).second);
195  p.offset(offset);
196  mx_id = p->id+1 > mx_id ? p->id+1 : mx_id;
197  }
198  offset = mx_id;
199 }
200 
201 static void rebaseVertices(Geant4PrimaryInteraction::VertexMap& vertices, int part_offset) {
202  Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
203  set<int> in, out;
204  set<int>::iterator i;
205  // Now move begin and end-vertex of all primary vertices accordingly
206  for(iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv) {
207  Geant4Vertex* v = (*iv).second;
208  in = v->in;
209  out = v->out;
210  for(in=v->in, i=in.begin(), v->in.clear(); i != in.end(); ++i)
211  v->in.insert((*i)+part_offset);
212  for(out=v->out, i=out.begin(), v->out.clear(); i != out.end(); ++i)
213  v->out.insert((*i)+part_offset);
214  }
215 }
216 
219  const Geant4Context* context)
220 {
221  typedef Geant4PrimaryEvent::Interaction Interaction;
222  typedef vector<Interaction*> Interactions;
223  Geant4Event& event = context->event();
225  Interaction* output = event.extension<Interaction>();
226  Interactions inter = evt->interactions();
227  int particle_offset = 0;
228 
229  for(Interactions::const_iterator i=inter.begin(); i != inter.end(); ++i) {
230  Interaction* interaction = *i;
231  int vertex_offset = particle_offset;
232  if ( !interaction->applyMask() ) {
233  caller->abortRun("Found single interaction with multiple primary vertices!",
234  "Cannot merge individual interactions with more than one primary!");
235  }
236  rebaseParticles(interaction->particles,particle_offset);
237  rebaseVertices(interaction->vertices,vertex_offset);
238  appendInteraction(caller,output,interaction);
239  }
240  output->setNextPID(particle_offset);
241  Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
242  caller->debug("+++ Merging MC input record from %d interactions:",(int)inter.size());
243  for( ip=output->particles.begin(), ipend=output->particles.end(); ip != ipend; ++ip )
244  Geant4ParticleHandle((*ip).second).dump1(DEBUG,caller->name(),"Merged particles");
245  return 1;
246 }
247 
251  double alpha)
252 {
253 #define SQR(x) (x*x)
254  Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
255  Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
256  double gamma = std::sqrt(1 + SQR(tan(alpha)));
257  double betagamma = std::tan(alpha);
258 
259  if ( inter->locked ) {
260  caller->abortRun("Locked interactions may not be boosted!",
261  "Cannot boost interactions with a native G4 primary record!");
262  }
263  else if ( alpha != 0.0 ) {
264  // Now move begin and end-vertex of all primary vertices accordingly
265  for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
266  Geant4Vertex* v = (*iv).second;
267  double t = gamma * v->time + betagamma * v->x / CLHEP::c_light;
268  double x = gamma * v->x + betagamma * CLHEP::c_light * v->time;
269  double y = v->y;
270  double z = v->z;
271  v->x = x;
272  v->y = y;
273  v->z = z;
274  v->time = t;
275  }
276 
277  // Now move begin and end-vertex of all primary and generator particles accordingly
278  for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
279  Geant4ParticleHandle p = (*ip).second;
280  double t = gamma * p->time + betagamma * p->vsx / CLHEP::c_light;
281  double x = gamma * p->vsx + betagamma * CLHEP::c_light * p->time;
282  double y = p->vsx;
283  double z = p->vsz;
284 
285  double m = p->mass;
286  double e2 = SQR(p->psx)+SQR(p->psy)+SQR(p->psz)+SQR(m);
287  double px = betagamma * std::sqrt(e2) + gamma * p->psx;
288  double py = p->psy;
289  double pz = p->psz;
290 
291  p->vsx = x;
292  p->vsy = y;
293  p->vsz = z;
294  p->time = t;
295 
296  p->psx = px;
297  p->psy = py;
298  p->psz = pz;
299  }
300  }
301  return 1;
302 }
303 
307  double dx, double dy, double dz, double dt)
308 {
309  Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
310  Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
311 
312  if ( inter->locked ) {
313  caller->abortRun("Locked interactions may not be smeared!",
314  "Cannot smear interactions with a native G4 primary record!");
315  }
316 
317  // Now move begin and end-vertex of all primary vertices accordingly
318  for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
319  Geant4Vertex* v = (*iv).second;
320  v->x += dx;
321  v->y += dy;
322  v->z += dz;
323  v->time += dt;
324  }
325  // Now move begin and end-vertex of all primary and generator particles accordingly
326  for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
327  Geant4Particle* p = (*ip).second;
328  p->vsx += dx;
329  p->vsy += dy;
330  p->vsz += dz;
331  p->vex += dx;
332  p->vey += dy;
333  p->vez += dz;
334  p->time += dt;
335  }
336  return 1;
337 }
338 
339 static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p) {
340  G4PrimaryParticle* g4 = 0;
341  if ( 0 != p->pdgID ) {
342  g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz);
343  }
344  else {
345  const G4ParticleDefinition* def = p.definition();
346  g4 = new G4PrimaryParticle(def, p->psx, p->psy, p->psz);
347  g4->SetCharge(p.charge());
348  }
349  g4->SetMass(p->mass);
350  return g4;
351 }
352 
353 static vector< pair<Geant4Particle*,G4PrimaryParticle*> >
354 getRelevant(set<int>& visited,
355  map<int,G4PrimaryParticle*>& prim,
357  const Geant4ParticleHandle p)
358 {
359  typedef vector< pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
360  Primaries res;
361  visited.insert(p->id);
362  PropertyMask status(p->status);
363  if ( status.isSet(G4PARTICLE_GEN_STABLE) ) {
364  if ( prim.find(p->id) == prim.end() ) {
365  G4PrimaryParticle* p4 = createG4Primary(p);
366  prim[p->id] = p4;
367  res.push_back(make_pair(p,p4));
368  }
369  }
370  else if ( p->daughters.size() > 0 ) {
371  const Geant4Particle::Particles& dau = p->daughters;
372  int first_daughter = *(dau.begin());
373  Geant4ParticleHandle dp = pm[first_daughter];
374  double en = p.energy();
375  double me = en > std::numeric_limits<double>::epsilon() ? p->mass / en : 0.0;
376  // fix by S.Morozov for real != 0
377  double proper_time = fabs(dp->time-p->time) * me;
378  double proper_time_Precision = pow(10.,-DBL_DIG)*me*fmax(fabs(p->time),fabs(dp->time));
379  bool isProperTimeZero = ( proper_time <= proper_time_Precision ) ;
380  // -- remove original --- if (proper_time != 0) {
381  if ( !isProperTimeZero ) {
382  map<int,G4PrimaryParticle*>::iterator ip4 = prim.find(p->id);
383  G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
384  if ( !p4 ) {
385  p4 = createG4Primary(p);
386  p4->SetProperTime(proper_time);
387  prim[p->id] = p4;
388  Primaries daughters;
389  for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
390  if ( visited.find(*i) == visited.end() ) {
391  Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
392  daughters.insert(daughters.end(), tmp.begin(),tmp.end());
393  }
394  }
395  for(Primaries::iterator i=daughters.begin(); i!=daughters.end(); ++i)
396  p4->SetDaughter((*i).second);
397  }
398  res.push_back(make_pair(p,p4));
399  }
400  else {
401  for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
402  if ( visited.find(*i) == visited.end() ) {
403  Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
404  res.insert(res.end(), tmp.begin(),tmp.end());
405  }
406  }
407  }
408  }
409  return res;
410 }
411 
414  const Geant4Context* context,
415  G4Event* event)
416 {
417  typedef vector< pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
418  typedef Geant4PrimaryInteraction Interaction;
419  Geant4PrimaryMap* primaries = context->event().extension<Geant4PrimaryMap>();
420  Interaction* interaction = context->event().extension<Interaction>();
421  Interaction::ParticleMap& pm = interaction->particles;
422  Interaction::VertexMap& vm = interaction->vertices;
423  map<int,G4PrimaryParticle*> prim;
424  set<int> visited;
425 
426  if ( interaction->locked ) {
427  caller->abortRun("Locked interactions may not be used to generate primaries!",
428  "Cannot handle a native G4 primary record!");
429  return 0;
430  }
431  else {
432  Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
433  for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i) {
434  int num_part = 0;
435  Geant4Vertex* v = (*i).second;
436  G4PrimaryVertex* v4 = new G4PrimaryVertex(v->x,v->y,v->z,v->time);
437  event->AddPrimaryVertex(v4);
438  caller->print("+++++ G4PrimaryVertex at (%+.2e,%+.2e,%+.2e) [mm] %+.2e [ns]",
440  for(Geant4Vertex::Particles::const_iterator ip=v->out.begin(); ip!=v->out.end(); ++ip) {
441  Geant4ParticleHandle p = pm[*ip];
442  if ( p->daughters.size() > 0 ) {
443  PropertyMask mask(p->reason);
445  }
446  if ( p->parents.size() == 0 ) {
447  Primaries relevant = getRelevant(visited,prim,pm,p);
448  for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j) {
449  Geant4ParticleHandle r = (*j).first;
450  G4PrimaryParticle* p4 = (*j).second;
451  PropertyMask reason(r->reason);
452  char text[64];
453 
454  reason.set(G4PARTICLE_PRIMARY);
455  v4->SetPrimary(p4);
456  ::snprintf(text,sizeof(text),"-> G4Primary[%3d]",num_part);
457  r.dumpWithMomentum(caller->outputLevel()-1,caller->name(),text);
458  ++num_part;
459  }
460  }
461  }
462  }
463  for(map<int,G4PrimaryParticle*>::iterator i=prim.begin(); i!=prim.end(); ++i) {
464  Geant4ParticleHandle p = pm[(*i).first];
465  primaries->insert((*i).second,p);
466  }
467  }
468  return 1;
469 }
PrintLevel outputLevel() const
Access the output level.
Definition: Geant4Action.h:287
VertexMap vertices
The map of primary vertices for the particles.
Geant4Event & event() const
Access the geant4 event – valid only between BeginEvent() and EndEvent()!
Geant4Particle * addRef()
Increase reference count.
User event context for DDG4.
void offset(int off) const
Handlers.
void setNextPID(int value)
Set the next PID value.
Data structure to store the MC particle information.
void set(const T &m)
Definition: Primitives.h:535
Particles out
The list of outgoing particles.
Definition: Geant4Vertex.h:56
ReferenceBitMask< int > PropertyMask
Class modelling a single interaction with multiple primary vertices and particles.
Definition: Geant4Primary.h:96
Particles in
The list of incoming particles.
Definition: Geant4Vertex.h:58
Data structure to manipulate a bitmask held by reference and represented by an integer.
Definition: Primitives.h:526
#define SQR(x)
Class modelling a complete primary event with multiple interactions.
Data structure to map primaries to particles.
Definition: Geant4Primary.h:65
Data structure to store the MC vertex information.
Definition: Geant4Vertex.h:46
const std::string & name() const
Access name of the action.
Definition: Geant4Action.h:271
void print(const char *fmt,...) const
Support for messages with variable output level using output level.
static G4PrimaryParticle * createG4Primary(const Geant4ParticleHandle p)
int smearInteraction(const Geant4Action *caller, Geant4PrimaryEvent::Interaction *inter, double dx, double dy, double dz, double dt)
Smear the primary vertex of an interaction.
Particles parents
The list of daughters of this MC particle.
void * addExtension(void *ptr, const std::type_info &info, destruct_t dtor)
Add an extension object to the detector element.
ParticleMap particles
The map of particles participating in this primary interaction.
Data structure to map particles produced during the generation and the simulation.
const G4ParticleDefinition * definition() const
Access the Geant4 particle definition object (expensive!)
int boostInteraction(const Geant4Action *caller, Geant4PrimaryEvent::Interaction *inter, double alpha)
Boost particles of one interaction identified by its mask.
static vector< pair< Geant4Particle *, G4PrimaryParticle * > > getRelevant(set< int > &visited, map< int, G4PrimaryParticle * > &prim, Geant4PrimaryInteraction::ParticleMap &pm, const Geant4ParticleHandle p)
static void rebaseVertices(Geant4PrimaryInteraction::VertexMap &vertices, int part_offset)
T * extension(bool alert=true)
Access to type safe extension object. Exception is thrown if the object is invalid.
int mergeInteractions(const Geant4Action *caller, const Geant4Context *context)
Merge all interactions present in the context.
int generatePrimaries(const Geant4Action *caller, const Geant4Context *context, G4Event *event)
Generate all primary vertices corresponding to the merged interaction.
std::vector< Geant4PrimaryInteraction * > interactions() const
Retrieve all interactions.
static void collectPrimaries(Geant4PrimaryMap *pm, Geant4PrimaryInteraction *interaction, Geant4Vertex *particle_origine, G4PrimaryParticle *gp)
Helper to recursively build a DDG4 interaction from an existing G4 interaction (primary vertex) ...
int mask
Vertex mask to associate particles from collision.
Definition: Geant4Vertex.h:52
Data structure to access derived MC particle information.
static void rebaseParticles(Geant4PrimaryInteraction::ParticleMap &particles, int &offset)
void insert(G4PrimaryParticle *g4_particle, Geant4Particle *particle)
Add a new object pair (G4 primary particle, DDG4 particle) into the maps.
int mask
User mask to flag the interaction. Also unique identifier.
static void appendInteraction(const Geant4Action *caller, Geant4PrimaryInteraction *output, Geant4PrimaryInteraction *input)
Append input interaction to global output.
double energy() const
Scalar particle energy.
int locked
Flag that the event is locked for G4 native generators.
static const double c_light
Definition: DD4hepUnits.h:328
#define DEBUG
Generic context to extend user, run and event information.
View * v
Definition: MultiView.cpp:30
std::map< int, Particle * > ParticleMap
void abortRun(const std::string &exception, const char *fmt,...) const
Abort Geant4 Run by throwing a G4Exception with type RunMustBeAborted.
ExtensionHandle extension
User data extension if required.
void dump1(int level, const std::string &src, const char *tag) const
Various output formats:
static const double mm
Definition: DD4hepUnits.h:69
Geant4Vertex * createPrimary(const G4PrimaryVertex *g4)
Create a vertex object from it's G4 counterpart.
bool isSet(const T &m) const
Definition: Primitives.h:544
int generationInitialization(const Geant4Action *caller, const Geant4Context *context)
Initialize the generation of one event.
Default base class for all Geant 4 actions and derivates thereof.
Definition: Geant4Action.h:91
static const double ns
Definition: DD4hepUnits.h:122
TGeoShape TGeoMedium * m
Definition: Volumes.cpp:294
double charge() const
Geant4 charge of the particle.
void debug(const char *fmt,...) const
Support of debug messages.
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.