diff --git a/DDCore/include/DD4hep/Primitives.h b/DDCore/include/DD4hep/Primitives.h
index df05dd66d2bae20782cdc645fad14616210ce9da..c91bf84fae8e1e6f210eaca9af511651dff9288d 100644
--- a/DDCore/include/DD4hep/Primitives.h
+++ b/DDCore/include/DD4hep/Primitives.h
@@ -241,9 +241,15 @@ namespace DD4hep {
void set(const T& m) {
mask |= m;
}
+ void clear(const T& m) {
+ mask &= ~m;
+ }
bool isSet(const T& m) const {
return (mask&m) == m;
}
+ bool anySet(const T& m) const {
+ return (mask&m) != 0;
+ }
bool testBit(int bit) const {
T m = T(1)<<bit;
return isSet(m);
diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index 267958c7f4deb0117497b3e5be9ba044c81b2621..2be162aebb71568408a9b57f9436e1ea7d451981 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -52,56 +52,63 @@ def run():
gen = DDG4.GeneratorAction(kernel,"Geant4GeneratorActionInit/GenerationInit")
kernel.generatorAction().adopt(gen)
+ #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
"""
Generation of isotrope tracks using the DDG4 partcle gun:
- """
+
# Setup particle gun
gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=3)
gun.OutputLevel = 5 # generator_output_level
+ gun.Mask = 1
+ #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+ """
- gen = DDG4.GeneratorAction(kernel,"Geant4PrimaryConverter/GunConverter");
- gen.OutputLevel = 5 # generator_output_level
- kernel.generatorAction().adopt(gen)
-
+ #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
"""
Generation of isotrope tracks of a given multiplicity with overlay:
# First particle generator: pi+
gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropPi+");
- gen.particle = 'pi+'
- gen.energy = 100 * GeV
- gen.multiplicity = 2
- gen.mask = 1
+ gen.Particle = 'pi+'
+ gen.Energy = 100 * GeV
+ gen.Multiplicity = 2
+ gen.Mask = 1
kernel.generatorAction().adopt(gen)
# Install vertex smearing for this interaction
gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearPi+");
- gen.mask = 1
- gen.offset = (20*mm, 10*mm, 10*mm, 0*ns)
- gen.sigma = (4*mm, 1*mm, 1*mm, 0*ns)
+ gen.Mask = 1
+ gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns)
+ gen.Sigma = (4*mm, 1*mm, 1*mm, 0*ns)
kernel.generatorAction().adopt(gen)
# Second particle generator: e-
gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropE-");
- gen.particle = 'e-'
- gen.energy = 25 * GeV
- gen.multiplicity = 3
- gen.mask = 2
+ gen.Particle = 'e-'
+ gen.Energy = 25 * GeV
+ gen.Multiplicity = 3
+ gen.Mask = 2
kernel.generatorAction().adopt(gen)
# Install vertex smearing for this interaction
gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearE-");
- gen.mask = 2
- gen.offset = (-20*mm, -10*mm, -10*mm, 0*ns)
- gen.sigma = (12*mm, 8*mm, 8*mm, 0*ns)
+ gen.Mask = 2
+ gen.Offset = (-20*mm, -10*mm, -10*mm, 0*ns)
+ gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns)
kernel.generatorAction().adopt(gen)
+ #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
"""
+ #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV
"""
Generation of primary particles from LCIO input files
+ """
# First particle file reader
gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO1");
gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_1.stdhep"
- gen.OutputLevel = generator_output_level
+ #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/qq_gen_128_999.stdhep"
+ gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/smuonLR_PointK_3TeV_BS_noBkg_run0001.stdhep"
+ gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbbb_3TeV.stdhep"
+ gen.OutputLevel = 4 # generator_output_level
gen.Mask = 1
gen.MomentumScale = 0.1
kernel.generatorAction().adopt(gen)
@@ -116,7 +123,8 @@ def run():
# Second particle file reader
gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO2");
gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_2.stdhep"
- gen.OutputLevel = generator_output_level
+ #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbnn_3TeV_01.stdhep"
+ gen.OutputLevel = 4 # generator_output_level
gen.Mask = 2
gen.MomentumScale = 0.1
kernel.generatorAction().adopt(gen)
@@ -127,24 +135,24 @@ def run():
gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns)
gen.Sigma = (2*mm, 1*mm, 1*mm, 0*ns)
kernel.generatorAction().adopt(gen)
+ #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
# Merge all existing interaction records
gen = DDG4.GeneratorAction(kernel,"Geant4InteractionMerger/InteractionMerger")
- gen.OutputLevel = generator_output_level
+ gen.OutputLevel = 4 #generator_output_level
kernel.generatorAction().adopt(gen)
# Finally generate Geant4 primaries
gen = DDG4.GeneratorAction(kernel,"Geant4PrimaryHandler/PrimaryHandler")
- gen.OutputLevel = generator_output_level
+ gen.OutputLevel = 4 #generator_output_level
kernel.generatorAction().adopt(gen)
- """
# And handle the simulation particles.
part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler")
kernel.generatorAction().adopt(part)
#part.SaveProcesses = ['conv','Decay']
part.SaveProcesses = ['Decay']
- part.OutputLevel = generator_output_level
+ part.OutputLevel = 4 # generator_output_level
part.enableUI()
user = DDG4.Action(kernel,"Geant4UserParticleHandler/UserParticleHandler")
part.adopt(user)
diff --git a/DDG4/include/DDG4/Geant4InteractionMerger.h b/DDG4/include/DDG4/Geant4InteractionMerger.h
index 90e2a16909976f23276ea6738f9dae674b19309f..a71fa1a743016d300ddfbc77df349599fdde691a 100644
--- a/DDG4/include/DDG4/Geant4InteractionMerger.h
+++ b/DDG4/include/DDG4/Geant4InteractionMerger.h
@@ -22,12 +22,17 @@ namespace DD4hep {
*/
namespace Simulation {
+ // Forward declarations
+ class Geant4PrimaryInteraction;
+
/** Geant4Action to convert the particle information to Geant4
*
* @author M.Frank
* @version 1.0
*/
class Geant4InteractionMerger : public Geant4GeneratorAction {
+ /// Append input interaction to global output
+ void appendInteraction(Geant4PrimaryInteraction* output, Geant4PrimaryInteraction* input);
public:
/// Standard constructor
Geant4InteractionMerger(Geant4Context* context, const std::string& nam);
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index 0b023766c58f475729b6c7096430a3a4a09a8e02..5382f8ceff1a4bd470fd9a8df40d01270d2cea02 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -62,8 +62,9 @@ namespace DD4hep {
G4PARTICLE_GEN_DECAYED = 1<<2, // particle decayed in the generator
G4PARTICLE_GEN_DOCUMENTATION = 1<<3, // documentation line
- G4PARTICLE_GEN_CREATED = 1<<4, // Particle generated by generator
- G4PARTICLE_GEN_HISTORY = 1<<5, // Particle is not a primary, but generation history
+ G4PARTICLE_GEN_GENERATOR = // Particle comes from generator
+ ( G4PARTICLE_GEN_EMPTY+G4PARTICLE_GEN_STABLE+
+ G4PARTICLE_GEN_DECAYED+G4PARTICLE_GEN_DOCUMENTATION ),
G4PARTICLE_GEN_STATUS = 0x3FF, // Mask for generator status (bit 0...9)
// Simulation status of a given particle
@@ -94,7 +95,7 @@ namespace DD4hep {
typedef std::set<int> Particles;
/// Reference counter
int ref; //! not persistent
- int id, g4Parent, reason, usermask;
+ int id, g4Parent, reason, mask;
int steps, secondaries, pdgID;
int status, colorFlow[2];
char charge, _spare[3];
diff --git a/DDG4/include/DDG4/Geant4ParticleGun.h b/DDG4/include/DDG4/Geant4ParticleGun.h
index 52fb2ec7509cd7066a9c973e89b6d10ad133db8e..253ccca46c1db45fd160b0ff20e4c254a71995fb 100644
--- a/DDG4/include/DDG4/Geant4ParticleGun.h
+++ b/DDG4/include/DDG4/Geant4ParticleGun.h
@@ -14,7 +14,6 @@
// Forward declarations
class G4ParticleDefinition;
-class G4ParticleGun;
/*
* DD4hep namespace declaration
@@ -34,22 +33,23 @@ namespace DD4hep {
*/
class Geant4ParticleGun: public Geant4GeneratorAction {
protected:
- /// Position and shooting direction of the gun
+ /// Property: Position and shooting direction of the gun
ROOT::Math::XYZVector m_position, m_direction;
- /// Pointer to geant4 particle definition
- G4ParticleDefinition* m_particle;
- /// Pointer to the particle gun itself
- G4ParticleGun* m_gun;
- /// Particle energy
+ /// Property: Particle energy
double m_energy;
- /// Particle name
+ /// Property: Particle name
std::string m_particleName;
- /// Desired multiplicity of the particles to be shot
+ /// Property: Desired multiplicity of the particles to be shot
int m_multiplicity;
+ /// Property: Interaction mask indentifier
+ int m_mask;
+ /// Property: Isotrope particles?
+ bool m_isotrop;
+
+ /// Pointer to geant4 particle definition
+ G4ParticleDefinition* m_particle;
/// Shot number in sequence
int m_shotNo;
- /// Isotrope particles?
- bool m_isotrop;
public:
/// Standard constructor
Geant4ParticleGun(Geant4Context* context, const std::string& name);
diff --git a/DDG4/include/DDG4/Geant4Primary.h b/DDG4/include/DDG4/Geant4Primary.h
index 3bbf70131b28a0c5a16c228d4a011e4b37bd28e6..98628ebb4606357da7f3e84b94385771b73e844d 100644
--- a/DDG4/include/DDG4/Geant4Primary.h
+++ b/DDG4/include/DDG4/Geant4Primary.h
@@ -103,6 +103,8 @@ namespace DD4hep {
int nextPID();
/// Set the next PID value
void setNextPID(int value);
+ /// Apply mask to all contained vertices (max. 1) and particles
+ bool applyMask();
};
/// Class modelling a complete primary event with multiple interactions
diff --git a/DDG4/include/DDG4/Geant4PrimaryConverter.h b/DDG4/include/DDG4/Geant4PrimaryConverter.h
deleted file mode 100644
index 9fea3351202b49a27523ce39749ed8816aaaeb69..0000000000000000000000000000000000000000
--- a/DDG4/include/DDG4/Geant4PrimaryConverter.h
+++ /dev/null
@@ -1,47 +0,0 @@
-//====================================================================
-// AIDA Detector description implementation
-//--------------------------------------------------------------------
-//
-// Author : M.Frank
-//
-//====================================================================
-#ifndef DD4HEP_DDG4_GEANT4PRIMARYCONVERTER_H
-#define DD4HEP_DDG4_GEANT4PRIMARYCONVERTER_H
-
-// Framework include files
-#include "DDG4/Geant4GeneratorAction.h"
-#include "Math/Vector3D.h"
-
-/*
- * DD4hep namespace declaration
- */
-namespace DD4hep {
-
- /*
- * Simulation namespace declaration
- */
- namespace Simulation {
-
- /** @class Geant4Primaryconverter Geant4Primaryconverter.h DDG4/Geant4Primaryconverter.h
- *
- * Convert G4PrimaryVertix and G4PrimaryParticle constructs into
- * Geant4Primary records. Required for native G4 particle guns.
- *
- * @author M.Frank
- * @version 1.0
- */
- class Geant4PrimaryConverter: public Geant4GeneratorAction {
- protected:
- /// Interaction mask
- int m_mask;
- public:
- /// Standard constructor
- Geant4PrimaryConverter(Geant4Context* context, const std::string& name);
- /// Default destructor
- virtual ~Geant4PrimaryConverter();
- /// Callback to generate primary particles
- virtual void operator()(G4Event* event);
- };
- } // End namespace Simulation
-} // End namespace DD4hep
-#endif /* DD4HEP_DDG4_GEANT4PRIMARYCONVERTER_H */
diff --git a/DDG4/include/DDG4/Geant4Vertex.h b/DDG4/include/DDG4/Geant4Vertex.h
index 3ad19da25b892c4cb5b67001fef6711580e7feaa..222ea5752fea441eb1d76d36af110f071ec868af 100644
--- a/DDG4/include/DDG4/Geant4Vertex.h
+++ b/DDG4/include/DDG4/Geant4Vertex.h
@@ -44,6 +44,8 @@ namespace DD4hep {
typedef std::set<int> Particles;
/// Reference counter
int ref; //! not persistent
+ /// Vertex mask to associate particles from collision
+ int mask;
/// Vertex data
double x,y,z,time;
/// The list of outgoing particles
diff --git a/DDG4/lcio/LCIOInputAction.cpp b/DDG4/lcio/LCIOInputAction.cpp
index 45fb86588709bd071ac56229731d760e518646b0..b952b17b155f54ab058c7a5f4b9c998283635dae 100644
--- a/DDG4/lcio/LCIOInputAction.cpp
+++ b/DDG4/lcio/LCIOInputAction.cpp
@@ -15,16 +15,25 @@
#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4Context.h"
#include "G4ParticleTable.hh"
-#include "G4Event.hh"
#include "EVENT/MCParticle.h"
#include "EVENT/LCCollection.h"
-#include "Randomize.hh"
+
+#include "G4Event.hh"
using namespace std;
using namespace DD4hep::Simulation;
-
typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+namespace {
+ inline int GET_ENTRY(const map<EVENT::MCParticle*,int>& mcparts, EVENT::MCParticle* part) {
+ map<EVENT::MCParticle*,int>::const_iterator ip=mcparts.find(part);
+ if ( ip == mcparts.end() ) {
+ throw runtime_error("Unknown particle identifier look-up!");
+ }
+ return (*ip).second;
+ }
+}
+
/// Standard constructor
LCIOInputAction::LCIOInputAction(Geant4Context* context, const string& name)
: Geant4GeneratorAction(context,name), m_reader(0)
@@ -81,36 +90,14 @@ EVENT::LCCollection* LCIOInputAction::readParticles(int which) {
return parts;
}
-namespace {
- inline int GET_ENTRY(const map<EVENT::MCParticle*,int>& mcparts, EVENT::MCParticle* part) {
- map<EVENT::MCParticle*,int>::const_iterator ip=mcparts.find(part);
- if ( ip == mcparts.end() ) {
- throw runtime_error("Unknown particle identifier look-up!");
- }
- return (*ip).second;
- }
- struct FindVertex {
- double x,y,z,t;
- FindVertex(double xx, double yy, double zz, double tt) { x=xx; y=yy; z=zz; t=tt; }
- bool operator()(const Geant4Vertex* o) const {
- if ( fabs(o->x - x) > numeric_limits<double>::epsilon() ) return false;
- if ( fabs(o->y - y) > numeric_limits<double>::epsilon() ) return false;
- if ( fabs(o->z - z) > numeric_limits<double>::epsilon() ) return false;
- if ( fabs(o->time - t) > numeric_limits<double>::epsilon() ) return false;
- return true;
- }
- };
-}
-
/// Callback to generate primary particles
void LCIOInputAction::operator()(G4Event* event) {
Geant4Event& evt = context()->event();
- Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
- Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
- EVENT::LCCollection* primaries = readParticles(event->GetEventID());
+ Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
+ Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
+ EVENT::LCCollection* primaries = readParticles(event->GetEventID());
map<EVENT::MCParticle*,int> mcparts;
- vector<EVENT::MCParticle*> mcpcoll;
- vector<Particle*> col;
+ vector<EVENT::MCParticle*> mcpcoll;
prim->add(m_mask, inter);
// check if there is at least one particle
@@ -120,21 +107,25 @@ void LCIOInputAction::operator()(G4Event* event) {
// check if there is at least one particle
if ( NHEP == 0 ) return;
- col.resize(NHEP,0);
mcpcoll.resize(NHEP,0);
for(int i=0; i<NHEP; ++i ) {
EVENT::MCParticle* p=dynamic_cast<EVENT::MCParticle*>(primaries->getElementAt(i));
mcparts[p] = i;
mcpcoll[i] = p;
- col[i] = new Particle(i);
}
print("+++ Particle interaction with %d generator particles ++++++++++++++++++++++++",NHEP);
+ Geant4Vertex* vtx = new Geant4Vertex();
+ vtx->x = 0;
+ vtx->y = 0;
+ vtx->z = 0;
+ vtx->time = 0;
+ inter->vertices.insert(make_pair(m_mask,vtx));
// build collection of MCParticles
G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
for(int i=0; i<NHEP; ++i ) {
EVENT::MCParticle* mcp = mcpcoll[i];
- Geant4ParticleHandle p(col[i]);
+ Geant4ParticleHandle p(new Particle(i));
const double mom_scale = m_momScale;
const double *mom = mcp->getMomentum();
const double *vsx = mcp->getVertex();
@@ -143,9 +134,6 @@ void LCIOInputAction::operator()(G4Event* event) {
const int *color = mcp->getColorFlow();
G4ParticleDefinition* def = tab->FindParticle(mcp->getPDG());
PropertyMask status(p->status);
-
- p->reason = 0;
- p->usermask = 0;
p->pdgID = mcp->getPDG();
p->charge = int(mcp->getCharge()*3.0);
p->psx = mom_scale*mom[0]*GeV;
@@ -174,7 +162,6 @@ void LCIOInputAction::operator()(G4Event* event) {
p->parents.insert(GET_ENTRY(mcparts,par[k]));
int genStatus = mcp->getGeneratorStatus();
- status.set(G4PARTICLE_GEN_CREATED);
if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE);
if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
@@ -187,158 +174,15 @@ void LCIOInputAction::operator()(G4Event* event) {
if ( mcp->hasLeftDetector() ) status.set(G4PARTICLE_SIM_LEFT_DETECTOR);
if ( mcp->isStopped() ) status.set(G4PARTICLE_SIM_STOPPED);
if ( mcp->isOverlay() ) status.set(G4PARTICLE_SIM_OVERLAY);
+ if ( par.size() == 0 ) {
+ if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
+ vtx->in.insert(p->id); // Beam particles and primary quarks etc.
+ else
+ vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
+ }
inter->particles.insert(make_pair(p->id,p));
p.dump3(outputLevel()-1,name(),"+->");
}
-
- // put initial particles to the vertex
- vector<Geant4Vertex*> vertices;
- double time_end = 0.0;
-
- for( size_t i=0; i< col.size(); i++ ) {
- Particle* p = col[i];
- if ( p->time > time_end ) time_end = p->time;
- }
-
- // Create Primary Vertex object
- set<int> missing;
- for( size_t i=0; i< col.size(); i++ ) {
- Geant4ParticleHandle p(col[i]);
- PropertyMask reason(p->reason);
- PropertyMask status(p->status);
- bool empty = status.isSet(G4PARTICLE_GEN_EMPTY);
- bool stable = status.isSet(G4PARTICLE_GEN_STABLE);
- bool decayed = status.isSet(G4PARTICLE_GEN_DECAYED);
- bool docu = status.isSet(G4PARTICLE_GEN_DOCUMENTATION);
-
- if ( docu || empty || decayed ) {
- status.set(G4PARTICLE_GEN_HISTORY);
- }
- else if ( stable ) {
- missing.insert(i);
- reason.set(G4PARTICLE_PRIMARY);
- }
- else if ( !stable ) {
- double ene = p.energy();
- double proper_time = std::fabs((time_end - p->time) * p->mass) / ene;
- double precision = std::pow(10E0,-DBL_DIG)*std::max(time_end,p->time)*p->mass/ene;
- if ( proper_time >= precision ) {
- p->properTime = proper_time;
- }
- reason.set(G4PARTICLE_PRIMARY);
- missing.insert(i);
- }
- else { // Catchall: will not make it to the primary record.....
- status.set(G4PARTICLE_GEN_HISTORY);
- }
- }
- set<Geant4Particle*> missing_def;
- for(set<int>::iterator i=missing.begin(); i != missing.end(); ++i) {
- Geant4ParticleHandle p(col[*i]);
- Geant4Vertex* vtx = 0;
- PropertyMask reason(p->reason);
-
- if ( 0 == p->definition ) {
- missing_def.insert(p);
- continue;
- }
- vector<Geant4Vertex*>::iterator iv=find_if(vertices.begin(),vertices.end(),
- FindVertex(p->vsx,p->vsy,p->vsz,p->time));
- if ( iv == vertices.end() ) {
- vtx = new Geant4Vertex();
- vtx->x = p->vsx;
- vtx->y = p->vsy;
- vtx->z = p->vsz;
- vtx->time = p->time;
- vertices.push_back(vtx);
- inter->vertices.insert(make_pair(inter->vertices.size(),vtx));
- vtx->out.insert(*i);
- }
- else {
- (*iv)->out.insert(*i);
- }
- }
- if ( missing_def.size() > 0 ) {
- for(set<Geant4Particle*>::const_iterator i=missing_def.begin(); i!= missing_def.end(); ++i)
- Geant4ParticleHandle(*i).dump1(WARNING,name(),"!!!! Missing definition");
- //abortRun(issue(event->GetEventID()),"Error creating primary particles.",
- // "+++ Missing particle definitions for primary particles. Geant4 cannot handle this!");
- }
-
-#if 0
- // put initial particles to the vertex
- for( size_t i=0; i< col.size(); i++ ) {
- Particle* p = col[i];
- if ( p->parents.size()==0 ) {
- set<int> outgoing = relevantParticles(col, p);
- prim_vtx->out.insert(outgoing.begin(),outgoing.end());
- }
- }
-#endif
-}
-
-set<int> LCIOInputAction::relevantParticles(const vector<Particle*>& particles, Geant4ParticleHandle p) const {
- set<int> result;
- PropertyMask mask(p->status);
- // CASE1: if p is a stable particle: search for it in G4ParticleMap
- // if it is already there: get G4PrimaryParticle version of p from G4ParticleMap
- // else: create G4PrimaryParticle version of p and write it to G4ParticleMap
- // in the end: insert this single G4PrimaryParticle verion of p to the
- // relevant particle list and return this "list".
- if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) {
- result.insert(p->id);
- return result;
- }
-
- if ( p->daughters.size() > 0 ) {
- // calculate proper time
- Particle::Particles::iterator id;
- int idau = *p->daughters.begin();
- Particle* dp = particles[idau];
- double ten = 10;
- double mass = p->definition->GetPDGMass();
- double proper_time = std::fabs((dp->time - p->time) * mass) / p.energy();
- double aPrecision = dp->time * mass / p.energy();
- double bPrecision = p->time * mass / p.energy();
- double proper_time_precision = std::pow(ten,-DBL_DIG)
- * std::max(std::fabs(aPrecision),std::fabs(bPrecision));
- bool isProperTimeZero = ( proper_time <= proper_time_precision );
-
- // CASE2: if p is not stable: get first decay daughter and calculate the proper time of p
- // if proper time is not zero: particle is "relevant", since it carries vertex information
- // if p is already in G4ParticleMap: take it
- // otherwise: create G4PrimaryParticle version of p and write it to G4ParticleMap
- // now collect all relevant particles of all daughters and setup "relevant mother-
- // daughter-relations" between relevant decay particles and G4PrimaryParticle version of p
- // in the end: insert only the G4PrimaryParticle version of p to the
- // relevant particle list and return this "list". The added particle has now the correct pre-assigned
- // decay products and (hopefully) the right lifetime.
- if ( isProperTimeZero == false ) {
- result.insert(p->id);
- p->properTime = proper_time;
- for (id=p->daughters.begin(); id!=p->daughters.end(); ++id) {
- int pid = *id;
- dp = particles[pid];
- set<int> dau = relevantParticles(particles,dp);
- result.insert(pid);
- result.insert(dau.begin(),dau.end());
- }
- return result;
- }
- // CASE3: if p is not stable AND has zero lifetime: forget about p and retrieve all relevant
- // decay particles of all daughters of p. In this case this step of recursion is just there for
- // collecting all relevant decay products of daughters (and return them).
- for (id=p->daughters.begin(); id!=p->daughters.end(); ++id) {
- int pid = *id;
- dp = particles[pid];
- set<int> dau = relevantParticles(particles,dp);
- result.insert(dau.begin(),dau.end());
- }
- }
- // This logic sucks!
- // - What happens to decayed beauty particles with a displaced vertex, which are stable?
- // - What happend to anything stable comping NOT from the very primary vertex.....
- return result;
}
#include "DDG4/Factories.h"
diff --git a/DDG4/lcio/LCIOInputAction.h b/DDG4/lcio/LCIOInputAction.h
index 9014f266cf605774bdc196913113e77250a1cc61..36e9c3d246e552b41bc95733b2a58ea4fcfe7088 100644
--- a/DDG4/lcio/LCIOInputAction.h
+++ b/DDG4/lcio/LCIOInputAction.h
@@ -78,8 +78,6 @@ namespace DD4hep {
virtual ~LCIOInputAction();
/// Callback to generate primary particles
virtual void operator()(G4Event*);
-
- std::set<int> relevantParticles(const std::vector<Particle*>& particles, Geant4ParticleHandle p) const;
};
} /* End namespace lcio */
} /* End namespace DD4hep */
diff --git a/DDG4/plugins/Geant4Factories.cpp b/DDG4/plugins/Geant4Factories.cpp
index 899f3b1c37d57635641f200759f3ab213c03da91..e4cfc7fa9cb232f6537dbdd3ddc7012d914b2f44 100644
--- a/DDG4/plugins/Geant4Factories.cpp
+++ b/DDG4/plugins/Geant4Factories.cpp
@@ -91,10 +91,6 @@ DECLARE_GEANT4ACTION(Geant4InteractionMerger)
#include "DDG4/Geant4PrimaryHandler.h"
DECLARE_GEANT4ACTION(Geant4PrimaryHandler)
-//=============================
-#include "DDG4/Geant4PrimaryConverter.h"
-DECLARE_GEANT4ACTION(Geant4PrimaryConverter)
-
//=============================
#include "DDG4/Geant4TestActions.h"
namespace DD4hep { namespace Simulation {
diff --git a/DDG4/src/Geant4InteractionMerger.cpp b/DDG4/src/Geant4InteractionMerger.cpp
index e8b4489ae11ea0efef49115cabadf5c9c476981b..840c7efcc157c1486317cdb9dd562fd39e2a8831 100644
--- a/DDG4/src/Geant4InteractionMerger.cpp
+++ b/DDG4/src/Geant4InteractionMerger.cpp
@@ -58,15 +58,22 @@ void rebaseVertices(Geant4PrimaryInteraction::VertexMap& vertices, int part_offs
}
}
-void appendInteraction(Geant4PrimaryInteraction* output, Geant4PrimaryInteraction* input) {
+/// Append input interaction to global output
+void Geant4InteractionMerger::appendInteraction(Geant4PrimaryInteraction* output, Geant4PrimaryInteraction* input) {
Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
for( ip=input->particles.begin(), ipend=input->particles.end(); ip != ipend; ++ip ) {
Geant4Particle* p = (*ip).second;
output->particles.insert(make_pair(p->id,p->addRef()));
}
- Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
- for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv )
- output->vertices.insert(make_pair(output->vertices.size(),(*iv).second->addRef()));
+ Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
+ for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv ) {
+ ivfnd = output->vertices.find((*iv).second->mask);
+ if ( ivfnd != output->vertices.end() ) {
+ abortRun("Duplicate primary interaction identifier!",
+ "Cannot handle 2 interactions with identical identifiers!");
+ }
+ output->vertices.insert(make_pair((*iv).second->mask,(*iv).second->addRef()));
+ }
}
/// Event generation action callback
@@ -76,11 +83,15 @@ void Geant4InteractionMerger::operator()(G4Event* /* event */) {
Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
Interaction* output = context()->event().extension<Interaction>();
Interactions inter = evt->interactions();
- int particle_offset = 0, vertex_offset = 0;
+ int particle_offset = 0;
for(Interactions::const_iterator i=inter.begin(); i != inter.end(); ++i) {
Interaction* interaction = *i;
- vertex_offset = particle_offset;
+ int vertex_offset = particle_offset;
+ if ( !interaction->applyMask() ) {
+ abortRun("Found single interaction with multiple primary vertices!",
+ "Cannot merge individual interactions with more than one primary!");
+ }
rebaseParticles(interaction->particles,particle_offset);
rebaseVertices(interaction->vertices,vertex_offset);
appendInteraction(output,interaction);
diff --git a/DDG4/src/Geant4InteractionVertexSmear.cpp b/DDG4/src/Geant4InteractionVertexSmear.cpp
index 6a2c562f0572200a13cac1652d4397cb667a0ce7..68a0f824fd121fd21c8cd3f5aa1a7b3d245a22b6 100644
--- a/DDG4/src/Geant4InteractionVertexSmear.cpp
+++ b/DDG4/src/Geant4InteractionVertexSmear.cpp
@@ -56,7 +56,9 @@ void Geant4InteractionVertexSmear::operator()(G4Event*) {
double dz = rndm.gauss(m_offset.z(),m_sigma.z());
double dt = rndm.gauss(m_offset.t(),m_sigma.t());
- print("+++ Smearing primary vertex for interaction type %d by (%+.2e mm, %+.2e mm, %+.2e mm, %+.2e ns)",m_mask,dx,dy,dz,dt);
+ print("+++ Smearing primary vertex for interaction type %d (%d Vertices, %d particles) "
+ "by (%+.2e mm, %+.2e mm, %+.2e mm, %+.2e ns)",
+ m_mask,int(inter->vertices.size()),int(inter->particles.size()),dx,dy,dz,dt);
// Now move begin and end-vertex of all primary vertices accordingly
for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
diff --git a/DDG4/src/Geant4IsotropeGenerator.cpp b/DDG4/src/Geant4IsotropeGenerator.cpp
index bd6d4ef6412ee3a3eec0ff5242ce8cf81612fa10..dc4ef90f711d3ae82ad666a5dc8f8ee3f22e239a 100644
--- a/DDG4/src/Geant4IsotropeGenerator.cpp
+++ b/DDG4/src/Geant4IsotropeGenerator.cpp
@@ -32,10 +32,10 @@ Geant4IsotropeGenerator::Geant4IsotropeGenerator(Geant4Context* context, const s
{
InstanceCount::increment(this);
m_needsControl = true;
- declareProperty("particle", m_particleName = "e-");
- declareProperty("energy", m_energy = 50 * MeV);
- declareProperty("multiplicity", m_multiplicity = 1);
- declareProperty("mask", m_mask = 0);
+ declareProperty("Particle", m_particleName = "e-");
+ declareProperty("Energy", m_energy = 50 * MeV);
+ declareProperty("Multiplicity", m_multiplicity = 1);
+ declareProperty("Mask", m_mask = 0);
}
/// Default destructor
@@ -61,6 +61,7 @@ void Geant4IsotropeGenerator::operator()(G4Event*) {
prim->add(m_mask, inter);
Geant4Vertex* vtx = new Geant4Vertex();
+ vtx->mask = m_mask;
inter->vertices.insert(make_pair(inter->vertices.size(),vtx));
for(int i=0; i<m_multiplicity; ++i) {
double phi = 2*M_PI*rnd.rndm();
@@ -72,11 +73,8 @@ void Geant4IsotropeGenerator::operator()(G4Event*) {
Particle* p = new Particle();
p->id = inter->nextPID();
- p->reason |= G4PARTICLE_PRIMARY;
- p->status |= G4PARTICLE_GEN_CREATED;
p->status |= G4PARTICLE_GEN_STABLE;
- p->usermask = m_mask;
- p->reason = G4PARTICLE_PRIMARY;
+ p->mask = m_mask;
p->pdgID = m_particle->GetPDGEncoding();
p->definition = m_particle;
p->psx = x1*momentum;
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
index bd2953641906fa007ad96e61b0c0857136e46f8f..cb70c4434cad0a5a2fda90cf5036da0b3daffebf 100644
--- a/DDG4/src/Geant4Particle.cpp
+++ b/DDG4/src/Geant4Particle.cpp
@@ -27,7 +27,7 @@ ParticleExtension::~ParticleExtension() {
/// Copy constructor
Geant4Particle::Geant4Particle(const Geant4Particle& c)
- : ref(1), id(c.id), g4Parent(c.g4Parent), reason(c.reason), usermask(c.usermask),
+ : ref(1), id(c.id), g4Parent(c.g4Parent), reason(c.reason), mask(c.mask),
steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID),
status(c.status), charge(0),
vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
@@ -48,7 +48,7 @@ Geant4Particle::Geant4Particle(const Geant4Particle& c)
/// Default constructor
Geant4Particle::Geant4Particle()
- : ref(1), id(0), g4Parent(0), reason(0), usermask(0),
+ : ref(1), id(0), g4Parent(0), reason(0), mask(0),
steps(0), secondaries(0), pdgID(0),
status(0), charge(0),
vsx(0.0), vsy(0.0), vsz(0.0),
@@ -65,7 +65,7 @@ Geant4Particle::Geant4Particle()
/// Constructor with ID initialization
Geant4Particle::Geant4Particle(int part_id)
- : ref(1), id(part_id), g4Parent(0), reason(0), usermask(0),
+ : ref(1), id(part_id), g4Parent(0), reason(0), mask(0),
steps(0), secondaries(0), pdgID(0),
status(0), charge(0),
vsx(0.0), vsy(0.0), vsz(0.0),
@@ -99,7 +99,7 @@ Geant4Particle& Geant4Particle::get_data(Geant4Particle& c) {
id = c.id;
g4Parent = c.g4Parent;
reason = c.reason;
- usermask = c.usermask;
+ mask = c.mask;
status = c.status;
charge = c.charge;
steps = c.steps;
@@ -178,8 +178,10 @@ std::string Geant4ParticleHandle::particleType() const {
std::string Geant4ParticleHandle::processName() const {
if ( particle->process ) return particle->process->GetProcessName();
else if ( particle->reason&G4PARTICLE_PRIMARY ) return "Primary";
- else if ( particle->status&G4PARTICLE_GEN_HISTORY ) return "Gen.History";
- else if ( particle->status&G4PARTICLE_GEN_CREATED ) return "Generator";
+ else if ( particle->status&G4PARTICLE_GEN_EMPTY ) return "Gen.Empty";
+ else if ( particle->status&G4PARTICLE_GEN_STABLE ) return "Gen.Stable";
+ else if ( particle->status&G4PARTICLE_GEN_DECAYED ) return "Gen.Decay";
+ else if ( particle->status&G4PARTICLE_GEN_DOCUMENTATION ) return "Gen.DOC";
return "???";
}
@@ -240,7 +242,7 @@ void Geant4ParticleHandle::dump2(int level, const std::string& src, const char*
else if ( p->parents.size() == 1 ) ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
else if ( p->parents.size() > 1 ) ::snprintf(text,sizeof(text),"/%d..",*(p->parents.begin()));
printout((DD4hep::PrintLevel)level,src,
- "+++ %s %4d G4:%4d def [%-11s,%8s] reason:%8d "
+ "+++ %s %4d G4:%4d [%-12s,%8s] reason:%8d "
"E:%+.2e in record:%s #Par:%3d%-5s #Dau:%3d",
tag, p->id, g4id,
p.particleName().c_str(),
diff --git a/DDG4/src/Geant4ParticleGun.cpp b/DDG4/src/Geant4ParticleGun.cpp
index 86af0110c390ed9039554d97d12a9b5e7d85f2cf..290545cf7f6a2842a0161f82976e578f10fde7c3 100644
--- a/DDG4/src/Geant4ParticleGun.cpp
+++ b/DDG4/src/Geant4ParticleGun.cpp
@@ -12,10 +12,10 @@
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4Context.h"
#include "DDG4/Geant4Random.h"
+#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4ParticleGun.h"
#include "CLHEP/Units/SystemOfUnits.h"
-#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
@@ -30,7 +30,7 @@ using namespace DD4hep::Simulation;
/// Standard constructor
Geant4ParticleGun::Geant4ParticleGun(Geant4Context* context, const string& name)
: Geant4GeneratorAction(context, name), m_position(0,0,0), m_direction(1,1,0.3),
- m_particle(0), m_gun(0), m_shotNo(0)
+ m_particle(0), m_shotNo(0)
{
InstanceCount::increment(this);
m_needsControl = true;
@@ -40,20 +40,19 @@ Geant4ParticleGun::Geant4ParticleGun(Geant4Context* context, const string& name)
declareProperty("position", m_position);
declareProperty("direction", m_direction);
declareProperty("isotrop", m_isotrop = false);
+ declareProperty("Mask", m_mask = 1);
}
/// Default destructor
Geant4ParticleGun::~Geant4ParticleGun() {
- if (m_gun)
- delete m_gun;
InstanceCount::decrement(this);
}
/// Callback to generate primary particles
-void Geant4ParticleGun::operator()(G4Event* event) {
- if (0 == m_gun) {
- m_gun = new G4ParticleGun(m_multiplicity);
- }
+void Geant4ParticleGun::operator()(G4Event* ) {
+ typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+ Geant4Event& evt = context()->event();
+ Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
if (0 == m_particle || m_particle->GetParticleName() != m_particleName.c_str()) {
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
m_particle = particleTable->FindParticle(m_particleName);
@@ -77,12 +76,44 @@ void Geant4ParticleGun::operator()(G4Event* event) {
}
}
print("Shoot [%d] %.3f GeV %s pos:(%.3f %.3f %.3f)[mm] dir:(%6.3f %6.3f %6.3f)",
- m_shotNo, m_energy/GeV, m_particleName.c_str(), m_position.X(), m_position.Y(), m_position.Z(),
+ m_shotNo, m_energy/GeV, m_particleName.c_str(),
+ m_position.X(), m_position.Y(), m_position.Z(),
m_direction.X(),m_direction.Y(), m_direction.Z());
- m_gun->SetParticleDefinition(m_particle);
- m_gun->SetParticleEnergy(m_energy);
- m_gun->SetParticleMomentumDirection(G4ThreeVector(m_direction.X(), m_direction.Y(), m_direction.Z()));
- m_gun->SetParticlePosition(G4ThreeVector(m_position.X(), m_position.Y(), m_position.Z()));
- m_gun->GeneratePrimaryVertex(event);
+
+ Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
+ Geant4Vertex* vtx = new Geant4Vertex();
+ inter->vertices.insert(make_pair(m_mask,vtx));
+ prim->add(m_mask, inter);
+
+ for(int i=0; i<m_multiplicity; ++i) {
+ Geant4ParticleHandle p = new Geant4Particle(i);
+ p->reason = 0;
+ p->pdgID = m_particle->GetPDGEncoding();
+ p->psx = m_direction.X()*m_energy;
+ p->psy = m_direction.Y()*m_energy;
+ p->psz = m_direction.Z()*m_energy;
+ p->time = 0;
+ p->properTime = 0;
+ p->vsx = m_position.X();
+ p->vsy = m_position.Y();
+ p->vsz = m_position.Z();
+ p->vex = m_position.X();
+ p->vey = m_position.Y();
+ p->vez = m_position.Z();
+ p->definition = m_particle;
+ p->process = 0;
+ p->spin[0] = 0;
+ p->spin[1] = 0;
+ p->spin[2] = 0;
+ p->colorFlow[0] = 0;
+ p->colorFlow[0] = 0;
+ p->mass = m_particle->GetPDGMass();
+ p->charge = m_particle->GetPDGCharge();
+ PropertyMask status(p->status);
+ status.set(G4PARTICLE_GEN_STABLE);
+ vtx->out.insert(p->id);
+ inter->particles.insert(make_pair(p->id,p));
+ p.dump3(outputLevel()-1,name(),"+->");
+ }
++m_shotNo;
}
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index c647d2ba3d85c05de67b4b561cde229c49d22c5f..2b559ff0984355c7fabaa474f8173a868097a586 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -219,7 +219,7 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
if ( prim_part ) {
m_currTrack.id = prim_part->id;
m_currTrack.reason = prim_part->reason|reason;
- m_currTrack.usermask = prim_part->usermask;
+ m_currTrack.mask = prim_part->mask;
m_currTrack.status = prim_part->status;
m_currTrack.spin[0] = prim_part->spin[0];
m_currTrack.spin[1] = prim_part->spin[1];
@@ -235,7 +235,7 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
else {
m_currTrack.id = m_globalParticleID;
m_currTrack.reason = reason;
- m_currTrack.usermask = 0;
+ m_currTrack.mask = 0;
m_currTrack.status |= G4PARTICLE_SIM_CREATED;
m_currTrack.spin[0] = 0;
m_currTrack.spin[1] = 0;
@@ -315,7 +315,7 @@ void Geant4ParticleHandler::end(const G4Track* track) {
if ( ipar != m_particleMap.end() ) {
Particle* p_par = (*ipar).second;
p_par->daughters.insert(ph->id);
- ph->usermask |= p_par->usermask;
+ ph->mask |= p_par->mask;
G4ThreeVector dist(ph->vsx-p_par->vex,ph->vsy-p_par->vey,ph->vsz-p_par->vez);
if ( dist.r() > 1e-15 ) {
// Set flag that the end vertex of the parent is not the start vertex of this track....
@@ -447,7 +447,8 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
p->daughters.clear();
for(set<int>::iterator id=daughters.begin(); id!=daughters.end(); ++id)
p->daughters.insert(orgParticles[*id]); // Requires (1)
- if ( 0 == (p->status&G4PARTICLE_GEN_HISTORY) ) {
+ //if ( 0 == (p->status&G4PARTICLE_GEN_GENERATOR) )
+ {
if ( p->g4Parent > 0 ) {
p->parents.clear();
p->parents.insert(equivalentTrack(p->g4Parent)); // Requires (2)
@@ -576,7 +577,7 @@ void Geant4ParticleHandler::checkConsistency() const {
}
// We assume that particles from the generator have consistent parents
// For all other particles except the primaries, the parent must be contained in the record.
- if ( !mask.isSet(G4PARTICLE_PRIMARY) && !status.isSet(G4PARTICLE_GEN_HISTORY) ) {
+ if ( !mask.isSet(G4PARTICLE_PRIMARY) && !status.anySet(G4PARTICLE_GEN_GENERATOR) ) {
int parent_id = equivalentTrack(p->g4Parent);
bool in_map = (j=m_particleMap.find(parent_id)) != m_particleMap.end();
bool in_parent_list = p->parents.find(parent_id) != p->parents.end();
diff --git a/DDG4/src/Geant4ParticlePrint.cpp b/DDG4/src/Geant4ParticlePrint.cpp
index ad6aa93956cc12edfb5695bd1eac80f0ac19a991..c81c0dcbddde3af3f8057f032e15836f2d3f186b 100644
--- a/DDG4/src/Geant4ParticlePrint.cpp
+++ b/DDG4/src/Geant4ParticlePrint.cpp
@@ -65,6 +65,7 @@ void Geant4ParticlePrint::end(const G4Event* ) {
void Geant4ParticlePrint::printParticle(const std::string& prefix, Geant4ParticleHandle p) const {
char equiv[32];
PropertyMask mask(p->reason);
+ PropertyMask status(p->status);
string proc_name = p.processName();
string proc_type = p.processTypeName();
int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
@@ -73,7 +74,7 @@ void Geant4ParticlePrint::printParticle(const std::string& prefix, Geant4Particl
if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
}
- print("+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s]",
+ print("+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s] %c%c%c%c",
prefix.c_str(),
p->id,
p.particleName().c_str(),
@@ -90,7 +91,11 @@ void Geant4ParticlePrint::printParticle(const std::string& prefix, Geant4Particl
p.numParent(),
proc_name.c_str(),
p->process ? "/" : "",
- proc_type.c_str()
+ proc_type.c_str(),
+ status.isSet(G4PARTICLE_GEN_EMPTY) ? 'E' : '.',
+ status.isSet(G4PARTICLE_GEN_STABLE) ? 'S' : '.',
+ status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
+ status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.'
);
}
diff --git a/DDG4/src/Geant4Primary.cpp b/DDG4/src/Geant4Primary.cpp
index bfd24d1b1604779470204c426b9a3f7be64ae229..38e8dbf0fca38fe23f30a3cd70302083ae416d19 100644
--- a/DDG4/src/Geant4Primary.cpp
+++ b/DDG4/src/Geant4Primary.cpp
@@ -62,6 +62,21 @@ void Geant4PrimaryInteraction::setNextPID(int new_value) {
next_particle_identifier = new_value-1;
}
+/// Apply mask to all contained vertices and particles
+bool Geant4PrimaryInteraction::applyMask() {
+ if ( vertices.size() <= 1 ) {
+ Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
+ for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip )
+ (*ip).second->mask = mask;
+
+ Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
+ for( iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv )
+ (*iv).second->mask = mask;
+ return true;
+ }
+ return false;
+}
+
/// Default constructor
Geant4PrimaryEvent::Geant4PrimaryEvent()
{
diff --git a/DDG4/src/Geant4PrimaryConverter.cpp b/DDG4/src/Geant4PrimaryConverter.cpp
deleted file mode 100644
index 4554a8fdf8e86db95b21ec53fdce51a9db35649a..0000000000000000000000000000000000000000
--- a/DDG4/src/Geant4PrimaryConverter.cpp
+++ /dev/null
@@ -1,101 +0,0 @@
-// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
-//====================================================================
-// AIDA Detector description implementation for LCD
-//--------------------------------------------------------------------
-//
-// Author : M.Frank
-//
-//====================================================================
-
-// Framework include files
-#include "DD4hep/Printout.h"
-#include "DD4hep/InstanceCount.h"
-#include "DDG4/Geant4Primary.h"
-#include "DDG4/Geant4Context.h"
-#include "DDG4/Geant4PrimaryConverter.h"
-
-#include "G4Event.hh"
-#include "G4PrimaryVertex.hh"
-#include "G4PrimaryParticle.hh"
-
-// C/C++ include files
-#include <limits>
-
-using namespace std;
-using namespace DD4hep::Simulation;
-
-/// Standard constructor
-Geant4PrimaryConverter::Geant4PrimaryConverter(Geant4Context* context, const string& name)
- : Geant4GeneratorAction(context, name)
-{
- InstanceCount::increment(this);
- declareProperty("Mask", m_mask = 1);
- m_needsControl = true;
-}
-
-/// Default destructor
-Geant4PrimaryConverter::~Geant4PrimaryConverter() {
- InstanceCount::decrement(this);
-}
-
-/// Callback to generate primary particles
-void Geant4PrimaryConverter::operator()(G4Event* event) {
- typedef DD4hep::ReferenceBitMask<int> PropertyMask;
- Geant4Event& evt = context()->event();
- Geant4PrimaryMap* primaries = evt.extension<Geant4PrimaryMap>();
- Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
- Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
- Geant4PrimaryInteraction* output = evt.extension<Geant4PrimaryInteraction>();
- const G4Event& e = *event;
-
- prim->add(m_mask, inter);
- print("+++ Particle interaction converted from Geant4 primaries:");
-
- for(int cnt=0, i=0, ni=e.GetNumberOfPrimaryVertex(); i<ni; ++i) {
- G4PrimaryVertex* v4 = e.GetPrimaryVertex(i);
- Geant4Vertex* vtx = new Geant4Vertex();
- vtx->x = v4->GetX0()*mm;
- vtx->y = v4->GetY0()*mm;
- vtx->z = v4->GetZ0()*mm;
- vtx->time = v4->GetT0()*ns;
- inter->vertices.insert(make_pair(inter->vertices.size(),vtx));
- output->vertices.insert(make_pair(output->vertices.size(),vtx->addRef()));
- for(int j=0, nj=v4->GetNumberOfParticle(); j<nj; ++j) {
- const G4PrimaryParticle* p4 = v4->GetPrimary(j);
- Geant4ParticleHandle p = new Geant4Particle(cnt);
- p->reason = 0;
- p->usermask = 0;
- p->pdgID = p4->GetPDGcode();
- p->psx = p4->GetPx()*GeV;
- p->psy = p4->GetPy()*GeV;
- p->psz = p4->GetPz()*GeV;
- p->time = v4->GetT0()*ns;
- p->properTime = v4->GetT0()*ns;
- p->vsx = v4->GetX0()*mm;
- p->vsy = v4->GetY0()*mm;
- p->vsz = v4->GetZ0()*mm;
- p->vex = v4->GetX0()*mm;
- p->vey = v4->GetY0()*mm;
- p->vez = v4->GetZ0()*mm;
- p->definition = p4->GetParticleDefinition();
- p->process = 0;
- p->spin[0] = 0;
- p->spin[1] = 0;
- p->spin[2] = 0;
- p->colorFlow[0] = 0;
- p->colorFlow[0] = 0;
- p->mass = p4->GetMass();
- p->charge = p4->GetCharge();
- PropertyMask status(p->status);
- status.set(G4PARTICLE_GEN_CREATED);
- status.set(G4PARTICLE_GEN_STABLE);
- p.dump3(outputLevel()-1,name(),"+->");
-
- vtx->out.insert(p->id);
- inter->particles.insert(make_pair(p->id,p));
- output->particles.insert(make_pair(p->id,p->addRef()));
- primaries->primaryMap.insert(make_pair(p4,p->addRef()));
- cnt++;
- }
- }
-}
diff --git a/DDG4/src/Geant4PrimaryHandler.cpp b/DDG4/src/Geant4PrimaryHandler.cpp
index f45ca9ed79b96929a075c65112b7b342fdfceb6f..4e677b3bef6403b7c207b28793659479d069cb0c 100644
--- a/DDG4/src/Geant4PrimaryHandler.cpp
+++ b/DDG4/src/Geant4PrimaryHandler.cpp
@@ -18,9 +18,13 @@
#include <stdexcept>
+using namespace std;
using namespace DD4hep;
using namespace DD4hep::Simulation;
+typedef Geant4PrimaryInteraction Interaction;
+typedef ReferenceBitMask<int> PropertyMask;
+
/// Standard constructor
Geant4PrimaryHandler::Geant4PrimaryHandler(Geant4Context* context, const std::string& nam)
: Geant4GeneratorAction(context,nam)
@@ -33,39 +37,108 @@ Geant4PrimaryHandler::~Geant4PrimaryHandler() {
InstanceCount::decrement(this);
}
+static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p) {
+ G4PrimaryParticle* g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz);
+ g4->SetMass(p->mass);
+ return g4;
+}
+
+typedef map<Geant4Particle*,G4PrimaryParticle*> Primaries;
+Primaries getRelevant(set<int>& visited,
+ map<int,G4PrimaryParticle*>& prim,
+ Interaction::ParticleMap& pm,
+ const Geant4ParticleHandle p)
+{
+ Primaries res;
+ visited.insert(p->id);
+ PropertyMask status(p->status);
+ if ( status.isSet(G4PARTICLE_GEN_STABLE) ) {
+ if ( prim.find(p->id) == prim.end() ) {
+ G4PrimaryParticle* p4 = createG4Primary(p);
+ prim[p->id] = p4;
+ res.insert(make_pair(p,p4));
+ }
+ }
+ else if ( p->daughters.size() > 0 ) {
+ const Geant4Particle::Particles& dau = p->daughters;
+ int first_daughter = *(dau.begin());
+ Geant4ParticleHandle dp = pm[first_daughter];
+ double me = p->mass / p.energy();
+ // fix by S.Morozov for real != 0
+ double proper_time = fabs(dp->time-p->time) * me;
+ double proper_time_Precision = pow(10,-DBL_DIG)*me*fmax(fabs(p->time),fabs(dp->time));
+ bool isProperTimeZero = ( proper_time <= proper_time_Precision ) ;
+ // -- remove original --- if (proper_time != 0) {
+ if ( !isProperTimeZero ) {
+ map<int,G4PrimaryParticle*>::iterator ip4 = prim.find(p->id);
+ G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
+ if ( !p4 ) {
+ p4 = createG4Primary(p);
+ p4->SetProperTime(proper_time);
+ prim[p->id] = p4;
+ Primaries daughters;
+ for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
+ if ( visited.find(*i) == visited.end() ) {
+ Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
+ daughters.insert(tmp.begin(),tmp.end());
+ }
+ }
+ for(Primaries::iterator i=daughters.begin(); i!=daughters.end(); ++i)
+ p4->SetDaughter((*i).second);
+ }
+ res.insert(make_pair(p,p4));
+ }
+ else {
+ for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
+ if ( visited.find(*i) == visited.end() ) {
+ Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
+ res.insert(tmp.begin(),tmp.end());
+ }
+ }
+ }
+ }
+ return res;
+}
+
/// Event generation action callback
void Geant4PrimaryHandler::operator()(G4Event* event) {
- typedef Geant4PrimaryInteraction Interaction;
Geant4PrimaryMap* primaries = context()->event().extension<Geant4PrimaryMap>();
Interaction* interaction = context()->event().extension<Interaction>();
Interaction::ParticleMap& pm = interaction->particles;
Interaction::VertexMap& vm = interaction->vertices;
- Interaction::VertexMap::const_iterator iv, ivend;
- int num_vtx = 0, num_part = 0;
+ map<int,G4PrimaryParticle*> prim;
+ set<int> visited;
- for(iv=vm.begin(), ivend=vm.end(); iv != ivend; ++iv, ++num_vtx, num_part=0) {
- Geant4Vertex* v = (*iv).second;
- G4PrimaryVertex* g4 = new G4PrimaryVertex(v->x,v->y,v->z,v->time);
- print("+++++ G4PrimaryVertex %3d at (%+.2e,%+.2e,%+.2e) [mm] %+.2e [ns] %d particles",
- num_vtx,v->x/mm,v->y/mm,v->z/mm,v->time/ns,int(v->out.size()));
- // Generate Geant4 primaries coming from this vertex
- for(Geant4Vertex::Particles::const_iterator j=v->out.begin(); j!=v->out.end(); ++j, ++num_part) {
- // Same particle cannot come from 2 vertices! Hence it must ALWAYS be recreated
- Interaction::ParticleMap::const_iterator ip = pm.find(*j);
- if ( ip == pm.end() ) { // ERROR. may not happen. Something went wrong in the gathering.
- const char* text = "+++ Fatal inconsistency in the Geant4PrimaryInteraction record.";
- printout(ERROR,name(),text);
- throw std::runtime_error(name()+std::string(" ")+text);
+ Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
+ for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i) {
+ int num_part = 0;
+ Geant4Vertex* v = (*i).second;
+ G4PrimaryVertex* v4 = new G4PrimaryVertex(v->x,v->y,v->z,v->time);
+ event->AddPrimaryVertex(v4);
+ print("+++++ G4PrimaryVertex at (%+.2e,%+.2e,%+.2e) [mm] %+.2e [ns]",
+ v->x/mm,v->y/mm,v->z/mm,v->time/ns);
+ for(Geant4Vertex::Particles::const_iterator ip=v->out.begin(); ip!=v->out.end(); ++ip) {
+ Geant4ParticleHandle p = pm[*ip];
+ if ( p->parents.size() == 0 ) {
+ Primaries relevant = getRelevant(visited,prim,pm,p);
+ for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j) {
+ Geant4ParticleHandle r = (*j).first;
+ G4PrimaryParticle* p4 = (*j).second;
+ PropertyMask reason(r->reason);
+ reason.set(G4PARTICLE_PRIMARY);
+ v4->SetPrimary(p4);
+ printM1("+++ +-> G4Primary[%3d] ID:%3d type:%9d/%-12s "
+ "Momentum:(%+.2e,%+.2e,%+.2e) [GeV] time:%+.2e [ns] #Par:%3d #Dau:%3d",
+ num_part,r->id,r->pdgID,r.particleName().c_str(),
+ r->psx/GeV,r->psy/GeV,r->psz/GeV,r->time/ns,
+ int(r->parents.size()),int(r->daughters.size()));
+ ++num_part;
+ }
}
- Geant4Particle* p = (*ip).second;
- G4PrimaryParticle* g4part = new G4PrimaryParticle(p->pdgID,p->psx,p->psy,p->psz);
- g4part->SetMass(p->mass);
- g4->SetPrimary(g4part);
- printM1("+++ +-> G4Primary[%3d] ID:%3d type:%9d Momentum:(%+.2e,%+.2e,%+.2e) [GeV] time:%+.2e [ns] #Par:%3d #Dau:%3d",
- num_part,p->id,p->pdgID,p->psx/GeV,p->psy/GeV,p->psz/GeV,p->time/ns,
- int(p->parents.size()),int(p->daughters.size()));
- primaries->primaryMap[g4part] = p->addRef();
}
- event->AddPrimaryVertex(g4);
+ }
+ for(map<int,G4PrimaryParticle*>::iterator i=prim.begin(); i!=prim.end(); ++i) {
+ Geant4ParticleHandle p = pm[(*i).first];
+ primaries->primaryMap[(*i).second] = p->addRef();
}
}
diff --git a/DDG4/src/Geant4UIManager.cpp b/DDG4/src/Geant4UIManager.cpp
index 97ebd1d9b0900587de40a854cac3cef4754f87a4..459626226c8322aae76b43588a25e0eb120edc11 100644
--- a/DDG4/src/Geant4UIManager.cpp
+++ b/DDG4/src/Geant4UIManager.cpp
@@ -108,6 +108,6 @@ void Geant4UIManager::start() {
/// Stop and release resources
void Geant4UIManager::stop() {
- deletePtr(m_vis);
- deletePtr(m_ui);
+ // deletePtr(m_vis);
+ //deletePtr(m_ui);
}
diff --git a/DDG4/src/Geant4Vertex.cpp b/DDG4/src/Geant4Vertex.cpp
index ee5658ccb3c55610b1ad55bb8ff047e07c1d53d9..d1dc57ef90358760fc813929a8fefb13e6bd06c8 100644
--- a/DDG4/src/Geant4Vertex.cpp
+++ b/DDG4/src/Geant4Vertex.cpp
@@ -21,14 +21,14 @@ VertexExtension::~VertexExtension() {
/// Copy constructor
Geant4Vertex::Geant4Vertex(const Geant4Vertex& c)
- : ref(1), x(c.x), y(c.y), z(c.z), time(c.time), out(c.out), in(c.in)
+ : ref(1), mask(c.mask), x(c.x), y(c.y), z(c.z), time(c.time), out(c.out), in(c.in)
{
InstanceCount::increment(this);
}
/// Default constructor
Geant4Vertex::Geant4Vertex()
- : ref(1), x(0), y(0), z(0), time(0)
+ : ref(1), mask(0), x(0), y(0), z(0), time(0)
{
InstanceCount::increment(this);
}
@@ -41,6 +41,7 @@ Geant4Vertex::~Geant4Vertex() {
/// Assignment operator
Geant4Vertex& Geant4Vertex::operator=(const Geant4Vertex& c) {
if ( this != &c ) {
+ mask = c.mask;
x = c.x;
y = c.y;
z = c.z;