20 #include "CLHEP/Units/SystemOfUnits.h"
21 #include "CLHEP/Units/PhysicalConstants.h"
24 #include "G4ParticleDefinition.hh"
26 #include "G4PrimaryVertex.hh"
27 #include "G4PrimaryParticle.hh"
34 using namespace DD4hep;
35 using namespace DD4hep::Simulation;
46 v->
time = g4->GetT0();
54 const G4PrimaryParticle* g4p)
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();
78 p->
mass = g4p->GetMass();
79 p->
charge = g4p->GetCharge();
89 G4PrimaryParticle* gp)
91 int pid = int(interaction->
particles.size());
93 G4PrimaryParticle* dau = gp->GetDaughter();
95 int mask = interaction->
mask;
97 interaction->
particles.insert(make_pair(p->id,p));
100 particle_origine->
out.insert(p->id);
106 int vid = int(interaction->
vertices.size());
111 dv->
in.insert(p->id);
112 interaction->
vertices.insert(make_pair(vid,dv));
113 for(; dau; dau = dau->GetNext())
122 const G4PrimaryVertex* gv)
126 int vid = int(interaction->
vertices.size());
127 interaction->
locked =
true;
128 interaction->
mask = mask;
130 interaction->
vertices.insert(make_pair(vid,v));
131 for (G4PrimaryParticle *gp = gv->GetPrimary(); gp; gp = gp->GetNext() )
173 Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
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!");
185 output->
vertices.insert(make_pair((*iv).second->mask,(*iv).second->addRef()));
190 Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
193 for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip ) {
196 mx_id = p->
id+1 > mx_id ? p->
id+1 : mx_id;
202 Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
204 set<int>::iterator i;
206 for(iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv) {
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);
222 typedef vector<Interaction*> Interactions;
225 Interaction* output =
event.
extension<Interaction>();
227 int particle_offset = 0;
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!");
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 )
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);
260 caller->
abortRun(
"Locked interactions may not be boosted!",
261 "Cannot boost interactions with a native G4 primary record!");
263 else if ( alpha != 0.0 ) {
287 double px = betagamma * std::sqrt(e2) + gamma * p->
psx;
307 double dx,
double dy,
double dz,
double dt)
309 Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
310 Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
313 caller->
abortRun(
"Locked interactions may not be smeared!",
314 "Cannot smear interactions with a native G4 primary record!");
340 G4PrimaryParticle* g4 = 0;
341 if ( 0 != p->
pdgID ) {
345 const G4ParticleDefinition* def = p.
definition();
346 g4 =
new G4PrimaryParticle(def, p->
psx, p->
psy, p->
psz);
347 g4->SetCharge(p.
charge());
349 g4->SetMass(p->
mass);
353 static vector< pair<Geant4Particle*,G4PrimaryParticle*> >
355 map<int,G4PrimaryParticle*>& prim,
359 typedef vector< pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
361 visited.insert(p->
id);
364 if ( prim.find(p->
id) == prim.end() ) {
367 res.push_back(make_pair(p,p4));
372 int first_daughter = *(dau.begin());
375 double me = en > std::numeric_limits<double>::epsilon() ? p->
mass / en : 0.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 ) ;
381 if ( !isProperTimeZero ) {
382 map<int,G4PrimaryParticle*>::iterator ip4 = prim.find(p->
id);
383 G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
386 p4->SetProperTime(proper_time);
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());
395 for(Primaries::iterator i=daughters.begin(); i!=daughters.end(); ++i)
396 p4->SetDaughter((*i).second);
398 res.push_back(make_pair(p,p4));
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());
417 typedef vector< pair<Geant4Particle*,G4PrimaryParticle*> > Primaries;
420 Interaction* interaction = context->
event().
extension<Interaction>();
421 Interaction::ParticleMap& pm = interaction->particles;
422 Interaction::VertexMap& vm = interaction->vertices;
423 map<int,G4PrimaryParticle*> prim;
426 if ( interaction->locked ) {
427 caller->
abortRun(
"Locked interactions may not be used to generate primaries!",
428 "Cannot handle a native G4 primary record!");
432 Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
433 for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i) {
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) {
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) {
450 G4PrimaryParticle* p4 = (*j).second;
456 ::snprintf(text,
sizeof(text),
"-> G4Primary[%3d]",num_part);
463 for(map<int,G4PrimaryParticle*>::iterator i=prim.begin(); i!=prim.end(); ++i) {
465 primaries->
insert((*i).second,p);
PrintLevel outputLevel() const
Access the output level.
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.
std::set< int > Particles
Particles out
The list of outgoing particles.
Class modelling a single interaction with multiple primary vertices and particles.
Particles in
The list of incoming particles.
Data structure to manipulate a bitmask held by reference and represented by an integer.
Class modelling a complete primary event with multiple interactions.
Data structure to map primaries to particles.
Data structure to store the MC vertex information.
const std::string & name() const
Access name of the action.
void print(const char *fmt,...) const
Support for messages with variable output level using output level.
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.
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.
int mask
Vertex mask to associate particles from collision.
Data structure to access derived MC particle information.
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.
double energy() const
Scalar particle energy.
int locked
Flag that the event is locked for G4 native generators.
static const double c_light
Generic context to extend user, run and event information.
std::map< int, Vertex * > VertexMap
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:
Geant4Vertex * createPrimary(const G4PrimaryVertex *g4)
Create a vertex object from it's G4 counterpart.
bool isSet(const T &m) const
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.
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.