From f96763a7399646d3b72a64a3eded177eaebc3fc7 Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Mon, 15 Sep 2014 09:44:42 +0000
Subject: [PATCH] Improces DDG4 MC thruth and LCIO I/O. Please see
doc/release.notes
---
DDCore/include/DD4hep/Primitives.h | 7 +-
DDEve/DDEve/DDG4IO.C | 28 +-
DDG4/CMakeLists.txt | 2 +-
DDG4/examples/CLICSidSimu.py | 110 ++++--
DDG4/include/DDG4/Geant4Action.h | 100 ++---
DDG4/include/DDG4/Geant4Context.h | 12 +-
DDG4/include/DDG4/Geant4Data.h | 95 ++---
DDG4/include/DDG4/Geant4DataConversion.h | 6 +
DDG4/include/DDG4/Geant4DataDump.h | 4 +-
DDG4/include/DDG4/Geant4GeneratorActionInit.h | 75 ++++
DDG4/include/DDG4/Geant4HitCollection.h | 27 +-
DDG4/include/DDG4/Geant4InteractionMerger.h | 42 +++
.../DDG4/Geant4InteractionVertexBoost.h | 49 +++
.../DDG4/Geant4InteractionVertexSmear.h | 51 +++
DDG4/include/DDG4/Geant4IsotropeGenerator.h | 55 +++
DDG4/include/DDG4/Geant4MonteCarloTruth.h | 3 +-
DDG4/include/DDG4/Geant4Particle.h | 303 +++++++++++++++
DDG4/include/DDG4/Geant4ParticleGun.h | 5 -
DDG4/include/DDG4/Geant4ParticleHandler.h | 30 +-
DDG4/include/DDG4/Geant4ParticlePrint.h | 25 +-
DDG4/include/DDG4/Geant4Primary.h | 150 ++++++++
DDG4/include/DDG4/Geant4PrimaryHandler.h | 42 +++
DDG4/include/DDG4/Geant4Random.h | 69 ++++
DDG4/include/DDG4/Geant4SensDetAction.h | 13 +-
DDG4/include/DDG4/Geant4UserParticleHandler.h | 4 +-
DDG4/include/DDG4/Geant4Vertex.h | 72 ++++
DDG4/lcio/EventReader.cpp | 303 ---------------
DDG4/lcio/EventReader.h | 125 -------
DDG4/lcio/Geant4Output2LCIO.cpp | 28 +-
DDG4/lcio/HepEventReader.cpp | 14 +-
DDG4/lcio/LCIOConversions.cpp | 159 ++++----
DDG4/lcio/LCIOEventReader.cpp | 20 +
DDG4/lcio/LCIOEventReader.h | 82 +++++
...LcioEventReader.cpp => LCIOFileReader.cpp} | 32 +-
DDG4/lcio/LCIOInputAction.cpp | 344 +++++++++++++++++
DDG4/lcio/LCIOInputAction.h | 87 +++++
DDG4/lcio/LCIOMCParticles.cpp | 150 ++++++++
...oStdHepReader.cpp => LCIOStdHepReader.cpp} | 35 +-
DDG4/plugins/Geant4EscapeCounter.cpp | 5 +-
DDG4/plugins/Geant4Factories.cpp | 22 ++
DDG4/python/DDG4Dict.C | 11 +-
DDG4/src/Geant4Action.cpp | 125 +++----
DDG4/src/Geant4Context.cpp | 4 +-
DDG4/src/Geant4Data.cpp | 107 +-----
DDG4/src/Geant4DataConversion.cpp | 14 +-
DDG4/src/Geant4DataDump.cpp | 7 +-
DDG4/src/Geant4Exec.cpp | 8 +-
DDG4/src/Geant4GeneratorActionInit.cpp | 81 ++++
DDG4/src/Geant4HitCollection.cpp | 10 +
DDG4/src/Geant4InteractionMerger.cpp | 93 +++++
DDG4/src/Geant4InteractionVertexBoost.cpp | 99 +++++
DDG4/src/Geant4InteractionVertexSmear.cpp | 84 +++++
DDG4/src/Geant4IsotropeGenerator.cpp | 92 +++++
DDG4/src/Geant4Output2ROOT.cpp | 14 +-
DDG4/src/Geant4Particle.cpp | 315 ++++++++++++++++
DDG4/src/Geant4ParticleGun.cpp | 18 +-
DDG4/src/Geant4ParticleHandler.cpp | 346 ++++++++++++------
DDG4/src/Geant4ParticlePrint.cpp | 121 +++---
DDG4/src/Geant4Primary.cpp | 119 ++++++
DDG4/src/Geant4PrimaryHandler.cpp | 71 ++++
DDG4/src/Geant4ROOTDump.cpp | 7 +-
DDG4/src/Geant4Random.cpp | 82 +++++
DDG4/src/Geant4SensDetAction.cpp | 6 +-
DDG4/src/Geant4TestActions.cpp | 32 +-
DDG4/src/Geant4Vertex.cpp | 61 +++
UtilityApps/src/run_plugin.h | 2 +-
doc/release.notes | 13 +
67 files changed, 3576 insertions(+), 1151 deletions(-)
create mode 100644 DDG4/include/DDG4/Geant4GeneratorActionInit.h
create mode 100644 DDG4/include/DDG4/Geant4InteractionMerger.h
create mode 100644 DDG4/include/DDG4/Geant4InteractionVertexBoost.h
create mode 100644 DDG4/include/DDG4/Geant4InteractionVertexSmear.h
create mode 100644 DDG4/include/DDG4/Geant4IsotropeGenerator.h
create mode 100644 DDG4/include/DDG4/Geant4Particle.h
create mode 100644 DDG4/include/DDG4/Geant4Primary.h
create mode 100644 DDG4/include/DDG4/Geant4PrimaryHandler.h
create mode 100644 DDG4/include/DDG4/Geant4Random.h
create mode 100644 DDG4/include/DDG4/Geant4Vertex.h
delete mode 100644 DDG4/lcio/EventReader.cpp
delete mode 100644 DDG4/lcio/EventReader.h
create mode 100644 DDG4/lcio/LCIOEventReader.cpp
create mode 100644 DDG4/lcio/LCIOEventReader.h
rename DDG4/lcio/{LcioEventReader.cpp => LCIOFileReader.cpp} (66%)
create mode 100644 DDG4/lcio/LCIOInputAction.cpp
create mode 100644 DDG4/lcio/LCIOInputAction.h
create mode 100644 DDG4/lcio/LCIOMCParticles.cpp
rename DDG4/lcio/{LcioStdHepReader.cpp => LCIOStdHepReader.cpp} (62%)
create mode 100644 DDG4/src/Geant4GeneratorActionInit.cpp
create mode 100644 DDG4/src/Geant4InteractionMerger.cpp
create mode 100644 DDG4/src/Geant4InteractionVertexBoost.cpp
create mode 100644 DDG4/src/Geant4InteractionVertexSmear.cpp
create mode 100644 DDG4/src/Geant4IsotropeGenerator.cpp
create mode 100644 DDG4/src/Geant4Particle.cpp
create mode 100644 DDG4/src/Geant4Primary.cpp
create mode 100644 DDG4/src/Geant4PrimaryHandler.cpp
create mode 100644 DDG4/src/Geant4Random.cpp
create mode 100644 DDG4/src/Geant4Vertex.cpp
diff --git a/DDCore/include/DD4hep/Primitives.h b/DDCore/include/DD4hep/Primitives.h
index f68088aec..df05dd66d 100644
--- a/DDCore/include/DD4hep/Primitives.h
+++ b/DDCore/include/DD4hep/Primitives.h
@@ -134,6 +134,11 @@ namespace DD4hep {
delete p;
p = 0;
}
+ /// Helper to delete objects from heap and reset the pointer. Saves many many lines of code
+ template <typename T> inline void deleteObject(T* p) {
+ if (0 != p)
+ delete p;
+ }
/// Helper to delete objects from heap and reset the pointer
template <typename T> inline void destroyObject(T*& p) {
deletePtr(p);
@@ -236,7 +241,7 @@ namespace DD4hep {
void set(const T& m) {
mask |= m;
}
- bool isSet(const T& m) {
+ bool isSet(const T& m) const {
return (mask&m) == m;
}
bool testBit(int bit) const {
diff --git a/DDEve/DDEve/DDG4IO.C b/DDEve/DDEve/DDG4IO.C
index 8c349f570..15dde5880 100644
--- a/DDEve/DDEve/DDG4IO.C
+++ b/DDEve/DDEve/DDG4IO.C
@@ -21,6 +21,7 @@ namespace DD4hep { namespace Geometry {
class G4Step;
class G4StepPoint;
#include "DDG4/Geant4Data.h"
+#include "DDG4/Geant4Particle.h"
#include "DDEve/DDEveEventData.h"
#ifdef __DD4HEP_DDEVE_EXCLUSIVE__
@@ -34,7 +35,7 @@ namespace DD4hep {
* Simulation namespace declaration
*/
namespace Simulation {
- #define NO_CALL { throw "This function shoule never ever be called!"; }
+#define NO_CALL { throw "This function shoule never ever be called!"; }
/// Default constructor
inline SimpleRun::SimpleRun() { }
/// Default destructor
@@ -46,15 +47,22 @@ namespace DD4hep {
/// Default destructor
inline DataExtension::~DataExtension() { }
+ /// Default destructor
+ inline DataExtension::~ParticleExtension() { }
/// Default constructor
- inline Particle::Particle() { }
+ inline Geant4Particle::Geant4Particle() { }
/// Copy constructor
- inline Particle::Particle(const Particle&) { NO_CALL }
+ inline Geant4Particle::Geant4Particle(const Geant4Particle&) { NO_CALL }
/// Default destructor
- inline Particle::~Particle() { }
+ inline Geant4Particle::~Geant4Particle() { }
/// Remove daughter from set
- inline void Particle::removeDaughter(int) { NO_CALL }
-
+ inline void Geant4Particle::removeDaughter(int) { NO_CALL }
+ /// Default destructor
+ inline Geant4PrimaryMap::~Geant4PrimaryMap();
+ /// Default destructor
+ inline Geant4ParticleMap::~Geant4ParticleMap() { }
+ /// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
+ inline int Geant4ParticleMap::particleID(int, bool) const { NO_CALL }
/// Default constructor
inline SimpleHit::SimpleHit() { }
/// Default destructor
@@ -92,8 +100,8 @@ namespace DD4hep {
#pragma link C++ class DD4hep::Simulation::SimpleEvent+;
#pragma link C++ class DD4hep::Simulation::DataExtension+;
-#pragma link C++ class DD4hep::Simulation::Particle+;
-#pragma link C++ class std::vector<DD4hep::Simulation::Particle*>+;
+#pragma link C++ class DD4hep::Simulation::Geant4Particle+;
+#pragma link C++ class std::vector<DD4hep::Simulation::Geant4Particle*>+;
#pragma link C++ class DD4hep::Simulation::SimpleHit+;
#pragma link C++ class std::vector<DD4hep::Simulation::SimpleHit*>+;
@@ -132,7 +140,7 @@ namespace {
if (source ) {
static TClass* cl_calo = gROOT->GetClass(typeid(SimpleCalorimeter::Hit));
static TClass* cl_tracker = gROOT->GetClass(typeid(SimpleTracker::Hit));
- //static TClass* cl_particles = gROOT->GetClass(typeid(Particle));
+ //static TClass* cl_particles = gROOT->GetClass(typeid(Geant4Particle));
void* result = 0;
SimpleHit* hit = (SimpleHit*)source;
const std::type_info& type = typeid(*hit);
@@ -145,7 +153,7 @@ namespace {
void* _convertParticleFunc(void* source, DDEveParticle* p) {
if (source ) {
- Particle* s = (Particle*)source;
+ Geant4Particle* s = (Geant4Particle*)source;
p->id = s->id;
p->vsx = s->vsx;
p->vsy = s->vsy;
diff --git a/DDG4/CMakeLists.txt b/DDG4/CMakeLists.txt
index 80b396893..da5eedb16 100644
--- a/DDG4/CMakeLists.txt
+++ b/DDG4/CMakeLists.txt
@@ -48,7 +48,7 @@ if(DD4HEP_USE_LCIO)
include_directories( ${LCIO_INCLUDE_DIRS} )
file(GLOB lcio_sources lcio/*.cpp)
add_library(DD4hepG4LCIO SHARED ${lcio_sources})
- target_link_libraries(DD4hepG4LCIO DD4hepCore DD4hepG4 ${Geant4_LIBRARIES} ${LCIO_LIBRARIES} ${ROOT_LIBRARIES} Reflex)
+ target_link_libraries(DD4hepG4LCIO DD4hepCore DD4hepG4 ${Geant4_LIBRARIES} ${LCIO_LIBRARIES} ${ROOT_LIBRARIES} Reflex EG)
SET_TARGET_PROPERTIES(DD4hepG4LCIO PROPERTIES VERSION ${DD4hep_VERSION} SOVERSION ${DD4hep_SOVERSION})
dd4hep_generate_rootmap(DD4hepG4LCIO)
endif()
diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index c09e83475..24f690fb4 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -36,42 +36,102 @@ def run():
kernel.runAction().adopt(run1)
# Configure Event actions
- evt2 = DDG4.EventAction(kernel,'Geant4TestEventAction/UserEvent_2')
- evt2.Property_int = 123454321
- evt2.Property_double = 5e15*GeV
- evt2.Property_string = 'Hello_2 from the python setup'
- evt2.enableUI()
- kernel.registerGlobalAction(evt2)
-
- evt1 = DDG4.EventAction(kernel,'Geant4TestEventAction/UserEvent_1')
- evt1.Property_int=01234
- evt1.Property_double=1e11
- evt1.Property_string='Hello_1'
- evt1.enableUI()
-
- kernel.eventAction().adopt(evt1)
- kernel.eventAction().adopt(evt2)
-
prt = DDG4.EventAction(kernel,'Geant4ParticlePrint/ParticlePrint')
- prt.OutputType = 3 # Print both: table and tree
+ prt.OutputLevel = Output.INFO
+ prt.OutputType = 3 # Print both: table and tree
kernel.eventAction().adopt(prt)
# Configure I/O
evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
+ evt_lcio.OutputLevel = Output.ERROR
+ # Setup particle gun
+ #gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=3)
- gen = DDG4.GeneratorAction(kernel,"Geant4TestGeneratorAction/Generate")
+ gen = DDG4.GeneratorAction(kernel,"Geant4GeneratorActionInit/GenerationInit")
kernel.generatorAction().adopt(gen)
- # Setup particle gun
- gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=3)
+ """
+ # First particle generator: pi+
+ gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropPi+");
+ 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)
+ 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
+ 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)
+ kernel.generatorAction().adopt(gen)
+ """
+
+ generator_output_level = Output.INFO
+
+ # 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.Mask = 1
+ gen.MomentumScale = 0.1
+ kernel.generatorAction().adopt(gen)
+ # Install vertex smearing for this interaction
+ gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/Smear1");
+ gen.OutputLevel = 4 #generator_output_level
+ gen.Mask = 1
+ gen.Offset = (-20*mm, -10*mm, -10*mm, 0*ns)
+ gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns)
+ kernel.generatorAction().adopt(gen)
+
+ # 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.Mask = 2
+ gen.MomentumScale = 0.1
+ kernel.generatorAction().adopt(gen)
+ # Install vertex smearing for this interaction
+ gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/Smear2");
+ gen.OutputLevel = generator_output_level
+ gen.Mask = 2
+ 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
+ kernel.generatorAction().adopt(gen)
+
+ # Finally generate Geant4 primaries
+ gen = DDG4.GeneratorAction(kernel,"Geant4PrimaryHandler/PrimaryHandler")
+ gen.OutputLevel = generator_output_level
+ kernel.generatorAction().adopt(gen)
- trk = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler")
- kernel.generatorAction().adopt(trk)
- trk.saveProcesses = ['conv','Decay']
- trk.enableUI()
+ # 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.enableUI()
user = DDG4.Action(kernel,"Geant4UserParticleHandler/UserParticleHandler")
- trk.adopt(user)
+ part.adopt(user)
"""
rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader")
diff --git a/DDG4/include/DDG4/Geant4Action.h b/DDG4/include/DDG4/Geant4Action.h
index 172395aff..1daffe7d0 100644
--- a/DDG4/include/DDG4/Geant4Action.h
+++ b/DDG4/include/DDG4/Geant4Action.h
@@ -67,61 +67,13 @@ namespace DD4hep {
static TypeName split(const std::string& type_name, const std::string& delim);
};
-#if 0
-
- /** @class Geant4UserTrajectory Geant4Action.h DDG4/Geant4Action.h
- *
- * @author M.Frank
- * @version 1.0
- */
- struct Geant4UserTrajectory {
- /// Default constructor
- Geant4UserTrajectory();
- /// Standard destructor
- virtual ~Geant4UserTrajectory();
- /// accessors
- virtual int trackID() const = 0;
- virtual int parentID() const = 0;
- virtual std::string particleName() const = 0;
- virtual double charge() const = 0;
- virtual int pdgID() const = 0;
- virtual G4ThreeVector momentum() const = 0;
- virtual int numPoints() const = 0;
- virtual G4VTrajectoryPoint* point(int i) const = 0;
- };
-
- /** @class Geant4Trajectory Geant4Action.h DDG4/Geant4Action.h
- *
- * @author M.Frank
- * @version 1.0
- */
- struct Geant4Trajectory : public G4VTrajectory {
- std::auto_ptr<Geant4UserTrajectory> trajectory;
- /// Default constructor
- Geant4Trajectory(Geant4UserTrajectory* traj);
- /// Standard destructor
- virtual ~Geant4Trajectory();
- /// Mandatory G4 overloads: Get/Set functions
- virtual G4int GetTrackID() const {return trajectory->trackID;}
- virtual G4int GetParentID() const {return trajectory->parentID();}
- virtual G4String GetParticleName() const {return trajectory->particleName();}
- /// Mandatory G4 overloads: Charge is that of G4DynamicParticle
- virtual G4double GetCharge() const {return trajectory->charge();}
- /// Mandatory G4 overloads: Zero will be returned if the particle does not have PDG code.
- virtual G4int GetPDGEncoding() const {return trajectory->pdgID();}
- /// Mandatory G4 overloads: Momentum at the origin of the track in global coordinate system.
- virtual G4ThreeVector GetInitialMomentum() const {return trajectory->momentum();}
-
- /// Mandatory G4 overloads: Returns the number of trajectory points
- virtual int GetPointEntries() const {return trajectory->numPoints();}
- virtual G4VTrajectoryPoint* GetPoint(G4int i) const {return trajectory->point(i);}
- };
-#endif
-
- /// Default base class for all geant 4 actions and derivates thereof.
+ /// Default base class for all Geant 4 actions and derivates thereof.
/**
- * @author M.Frank
- * @version 1.0
+ * This is a utility class supporting properties, output and access to
+ * event and run objects through the context.
+ *
+ * @author M.Frank
+ * @version 1.0
*/
class Geant4Action {
protected:
@@ -327,39 +279,33 @@ namespace DD4hep {
void installPropertyMessenger();
/// Support for messages with variable output level using output level
- void print(const std::string& fmt, ...) const;
- /// Support for messages with variable output level using output level
- void print(const std::string& tag, const std::string& fmt, ...) const;
+ void print(const char* fmt, ...) const;
+ /// Support for messages with variable output level using output level-1
+ void printM1(const char* fmt, ...) const;
+ /// Support for messages with variable output level using output level-2
+ void printM2(const char* fmt, ...) const;
+ /// Support for messages with variable output level using output level+1
+ void printP1(const char* fmt, ...) const;
+ /// Support for messages with variable output level using output level+2
+ void printP2(const char* fmt, ...) const;
/// Support of debug messages.
- void debug(const std::string& fmt, ...) const;
- /// Support of warning messages.
- void debug(const std::string& tag, const std::string& fmt, ...) const;
+ void debug(const char* fmt, ...) const;
/// Support of info messages.
- void info(const std::string& fmt, ...) const;
- /// Support of warning messages.
- void info(const std::string& tag, const std::string& fmt, ...) const;
+ void info(const char* fmt, ...) const;
/// Support of warning messages.
- void warning(const std::string& fmt, ...) const;
- /// Support of warning messages.
- void warning(const std::string& tag, const std::string& fmt, ...) const;
- /// Support of error messages.
- void error(const std::string& fmt, ...) const;
+ void warning(const char* fmt, ...) const;
/// Support of error messages.
- void error(const std::string& tag, const std::string& fmt, ...) const;
+ void error(const char* fmt, ...) const;
/// Action to support error messages.
- bool error(bool return_value, const std::string& fmt, ...) const;
+ bool error(bool return_value, const char* fmt, ...) const;
/// Support of fatal messages. Throws exception
- void fatal(const std::string& fmt, ...) const;
- /// Support of warning messages.
- void fatal(const std::string& tag, const std::string& fmt, ...) const;
+ void fatal(const char* fmt, ...) const;
/// Support of exceptions: Print fatal message and throw runtime_error.
- void except(const std::string& fmt, ...) const;
- /// Support of warning messages.
- void except(const std::string& tag, const std::string& fmt, ...) const;
+ void except(const char* fmt, ...) const;
/// Abort Geant4 Run by throwing a G4Exception with type RunMustBeAborted
- void abortRun(const std::string& exception, const std::string& fmt, ...) const;
+ void abortRun(const std::string& exception, const char* fmt, ...) const;
/// Access to the main run action sequence from the kernel object
Geant4RunActionSequence& runAction() const;
diff --git a/DDG4/include/DDG4/Geant4Context.h b/DDG4/include/DDG4/Geant4Context.h
index cd59f4aae..50a4c1c2e 100644
--- a/DDG4/include/DDG4/Geant4Context.h
+++ b/DDG4/include/DDG4/Geant4Context.h
@@ -39,6 +39,7 @@ namespace DD4hep {
class Geant4Run;
class Geant4Event;
class Geant4Kernel;
+ class Geant4Random;
class ContextUpdate;
class Geant4RunActionSequence;
class Geant4EventActionSequence;
@@ -110,21 +111,30 @@ namespace DD4hep {
* across different events.
* Hence: They are only useful to extend data of an event.
*
+ * Any random numbers used to process one event should be accessed
+ * from this location. The framework ensures that the same seeded
+ * sequence is used throughout the processing of one single event.
+ *
* @author M.Frank
* @version 1.0
*/
class Geant4Event : public ObjectExtensions {
/// Reference to the original Geant4 event object
const G4Event* m_event;
+ /// Reference to the main random number generator
+ Geant4Random* m_random;
+
public:
/// Intializing constructor
- Geant4Event(const G4Event* run);
+ Geant4Event(const G4Event* run, Geant4Random* rndm);
/// Default destructor
virtual ~Geant4Event();
/// Access the G4Event directly: Automatic type conversion
operator const G4Event&() const { return *m_event; }
/// Access the G4Event directly: Explicit G4Event accessor
const G4Event& event() const { return *m_event; }
+ /// Access the random number generator
+ Geant4Random& random() const { return *m_random; }
/// Add an extension object to the detector element
/** Note:
diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h
index 2009478be..3d51abafe 100644
--- a/DDG4/include/DDG4/Geant4Data.h
+++ b/DDG4/include/DDG4/Geant4Data.h
@@ -13,10 +13,14 @@
#ifndef __DD4HEP_DDEVE_EXCLUSIVE__
#include "DD4hep/Objects.h"
#include "G4Step.hh"
+class G4PrimaryParticle;
#else
-typedef void G4VProcess;
-typedef void G4ParticleDefinition;
+typedef void* G4Step;
+typedef void* G4StepPoint;
+typedef void* G4PrimaryParticle;
#endif
+
+// C/C++ include files
#include <set>
#include <memory>
@@ -30,10 +34,15 @@ namespace DD4hep {
*/
namespace Simulation {
-
+ // Forward type definitions
typedef Geometry::Position Position;
typedef Geometry::Position Direction;
+ // Forward declarations
+ class Geant4Particle;
+ class Geant4Vertex;
+
+
/** @class HitCompare Geant4Data.h DDG4/Geant4Data.h
*
* Base class for hit comparisons.
@@ -112,59 +121,17 @@ namespace DD4hep {
virtual ~DataExtension();
};
- /// Track properties
- enum ParticleProperties {
- G4PARTICLE_CREATED_HIT = 1<<1,
- G4PARTICLE_PRIMARY = 1<<2,
- G4PARTICLE_HAS_SECONDARIES = 1<<3,
- G4PARTICLE_ABOVE_ENERGY_THRESHOLD = 1<<4,
- G4PARTICLE_KEEP_PROCESS = 1<<5,
- G4PARTICLE_KEEP_PARENT = 1<<6,
- G4PARTICLE_CREATED_CALORIMETER_HIT = 1<<7,
- G4PARTICLE_CREATED_TRACKER_HIT = 1<<8,
- G4PARTICLE_KEEP_ALWAYS = 1<<9,
- G4PARTICLE_LAST_NOTHING = 1<<31
- };
-
- /// Data structure to store the MC particle information
- /**
- * @author M.Frank
- * @version 1.0
- */
- class Particle {
- private:
- /// Copy constructor
- Particle(const Particle& c);
- public:
- int id, g4Parent, parent, reason, steps, secondaries, pdgID;
- double vsx, vsy, vsz;
- double vex, vey, vez;
- double psx, psy, psz, pex, pey, pez, energy, time;
- /// The list of daughters of this MC particle
- std::set<int> daughters;
- /// User data extension if required
- std::auto_ptr<DataExtension> extension;
- const G4VProcess *process; //!
- const G4ParticleDefinition *definition; //!
- /// Default constructor
- Particle();
- /// Default destructor
- virtual ~Particle();
- /// Assignment operator
- Particle& get_data(Particle& c);
- /// Remove daughter from set
- void removeDaughter(int id_daughter);
- };
-
- /** @class SimpleHit Geant4Data.h DDG4/Geant4Data.h
+ /// Base class for geant4 hit structures
+ /*
+ * Base class for geant4 hit structures created by the
+ * example sensitive detectors. This is a generic class
+ * only dealing with the cellID. Users may add an extension
+ * object, which normally should not be necessary.
*
- * Base class for geant4 hit structures created by the
- * example sensitive detectors.
- *
- * @author M.Frank
- * @version 1.0
+ * @author M.Frank
+ * @version 1.0
*/
- struct SimpleHit {
+ struct Geant4HitData {
// cellID
long long int cellID;
/// User data extension if required
@@ -217,15 +184,15 @@ namespace DD4hep {
typedef std::vector<MonteCarloContrib> Contributions;
public:
/// Default constructor
- SimpleHit();
+ Geant4HitData();
/// Default destructor
- virtual ~SimpleHit();
+ virtual ~Geant4HitData();
/// Extract the MC contribution for a given hit from the step information
static Contribution extractContribution(G4Step* step);
};
- struct SimpleTracker {
- /** @class SimpleTracker::Hit Geant4Data.h DDG4/Geant4Data.h
+ struct Geant4Tracker {
+ /** @class Geant4Tracker::Hit Geant4Data.h DDG4/Geant4Data.h
*
* Geant4 tracker hit class. Tracker hits contain the momentum
* direction as well as the hit position.
@@ -233,7 +200,7 @@ namespace DD4hep {
* @author M.Frank
* @version 1.0
*/
- struct Hit : public SimpleHit {
+ struct Hit : public Geant4HitData {
/// Hit position
Position position;
/// Hit direction
@@ -259,8 +226,8 @@ namespace DD4hep {
};
};
- struct SimpleCalorimeter {
- /** @class SimpleCalorimeter::Hit Geant4Data.h DDG4/Geant4Data.h
+ struct Geant4Calorimeter {
+ /** @class Geant4Calorimeter::Hit Geant4Data.h DDG4/Geant4Data.h
*
* Geant4 tracker hit class. Calorimeter hits contain the momentum
* direction as well as the hit position.
@@ -268,7 +235,7 @@ namespace DD4hep {
* @author M.Frank
* @version 1.0
*/
- struct Hit : public SimpleHit {
+ struct Hit : public Geant4HitData {
public:
/// Hit position
Position position;
@@ -285,6 +252,10 @@ namespace DD4hep {
virtual ~Hit();
};
};
+ typedef Geant4HitData SimpleHit;
+ typedef Geant4Tracker SimpleTracker;
+ typedef Geant4Calorimeter SimpleCalorimeter;
+
} // End namespace Simulation
} // End namespace DD4hep
#endif // DD4HEP_GEANT4DATA_H
diff --git a/DDG4/include/DDG4/Geant4DataConversion.h b/DDG4/include/DDG4/Geant4DataConversion.h
index 9174c92b8..44a39ff07 100644
--- a/DDG4/include/DDG4/Geant4DataConversion.h
+++ b/DDG4/include/DDG4/Geant4DataConversion.h
@@ -11,6 +11,7 @@
// Framework include files
#include "DD4hep/VolumeManager.h"
+#include "DD4hep/Detector.h"
#include <typeinfo>
/*
@@ -33,7 +34,12 @@ namespace DD4hep {
Geant4ConversionHelper();
/// Default destructor
virtual ~Geant4ConversionHelper();
+ /// Access to the data encoding using the volume manager and a specified volume id
static std::string encoding(Geometry::VolumeManager vm, Geometry::VolumeManager::VolumeID vid);
+ /// Access to the hit encoding in this sensitive detector
+ static std::string encoding(Geometry::Handle<Geometry::SensitiveDetectorObject> sd);
+ /// Access to the hit encoding in this readout object
+ static std::string encoding(Geometry::Readout ro);
};
/** @class Geant4Conversion Geant4DataConversion.h DDG4/Geant4DataConversion.h
diff --git a/DDG4/include/DDG4/Geant4DataDump.h b/DDG4/include/DDG4/Geant4DataDump.h
index e754e65b4..4d7065782 100644
--- a/DDG4/include/DDG4/Geant4DataDump.h
+++ b/DDG4/include/DDG4/Geant4DataDump.h
@@ -12,6 +12,7 @@
// Framework include files
#include "DD4hep/Printout.h"
#include "DDG4/Geant4Data.h"
+#include "DDG4/Geant4Particle.h"
// C/C++ include files
#include <vector>
@@ -34,6 +35,7 @@ namespace DD4hep {
*/
class Geant4DataDump {
public:
+ typedef Geant4Particle Particle;
typedef std::vector<Particle*> Particles;
typedef SimpleTracker::Hit TrackerHit;
@@ -53,7 +55,7 @@ namespace DD4hep {
virtual ~Geant4DataDump();
/// Print a single particle to the output logging using the specified print level
- void print(PrintLevel level, const Particle* p) const;
+ void print(PrintLevel level, Geant4ParticleHandle p) const;
/// Print the particle container to the output logging using the specified print level
void print(PrintLevel level, const std::string& container, const Particles* parts);
diff --git a/DDG4/include/DDG4/Geant4GeneratorActionInit.h b/DDG4/include/DDG4/Geant4GeneratorActionInit.h
new file mode 100644
index 000000000..700b6817f
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4GeneratorActionInit.h
@@ -0,0 +1,75 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4GENERATORACTIONINIT_H
+#define DD4HEP_DDG4_GEANT4GENERATORACTIONINIT_H
+
+// Framework include files
+#include "DDG4/Geant4GeneratorAction.h"
+
+// Forward declarations
+class G4Event;
+class G4Run;
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ /** Geant4Action to collect the MC particle information.
+ *
+ * This action should register all event extension required for the further
+ * processing. We want to avoid that every client has to check if a given
+ * object is present or not and than later install the required data structures.
+ *
+ * These by default are extensions of type:
+ * -- Geant4PrimaryEvent with multiple interaction sections, one for each interaction
+ * This is the MAIN and ONLY information to feed Geant4
+ *
+ * -- Geant4PrimaryInteraction containing the track/vertex information to create
+ * the primary particles for Geant4. This record is build from the Geant4PrimaryEvent
+ * information.
+ * -- Geant4PrimaryMap a map of the GEant4Particles converted to G4PrimaryParticles
+ * to ease particle handling later.
+ * -- Geant4ParticleMap the map of partcles created during the event simulation.
+ * This map has directly the correct particle offsets, so that the merging of
+ * Geant4PrimaryInteraction particles and the simulation particles is easy....
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4GeneratorActionInit : public Geant4GeneratorAction {
+ protected:
+ /// Current run identifier
+ int m_run;
+ /// Counter for total number of events
+ int m_evtTotal;
+ /// Counter for total number of events in current run
+ int m_evtRun;
+ public:
+ /// Standard constructor
+ Geant4GeneratorActionInit(Geant4Context* context, const std::string& nam);
+ /// Default destructor
+ virtual ~Geant4GeneratorActionInit();
+ /// Event generation action callback
+ virtual void operator()(G4Event* event);
+ /// Begin-run action callback
+ void begin(const G4Run* run);
+ /// End-run action callback
+ void end(const G4Run* run);
+
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+
+#endif // DD4HEP_DDG4_GEANT4GENERATORACTIONINIT_H
diff --git a/DDG4/include/DDG4/Geant4HitCollection.h b/DDG4/include/DDG4/Geant4HitCollection.h
index 198677724..1d13183e9 100644
--- a/DDG4/include/DDG4/Geant4HitCollection.h
+++ b/DDG4/include/DDG4/Geant4HitCollection.h
@@ -10,6 +10,7 @@
#define DD4HEP_DDG4_GEANT4HITCOLLECTION_H
// Framework include files
+#include "DD4hep/Handle.h"
#include "DDG4/ComponentUtils.h"
#include "G4VHitsCollection.hh"
#include "G4VHit.hh"
@@ -25,6 +26,9 @@
*/
namespace DD4hep {
+ /// Forward declarations
+ namespace Geometry { class SensitiveDetectorObject; }
+
/*
* Simulation namespace declaration
*/
@@ -167,8 +171,13 @@ namespace DD4hep {
*/
class Geant4HitCollection: public G4VHitsCollection {
public:
- typedef std::vector<Geant4HitWrapper> WrappedHits;
+ /** Local type declarations */
+ /// Hit wrapper
+ typedef std::vector<Geant4HitWrapper> WrappedHits;
+ /// Hit manipulator
typedef Geant4HitWrapper::HitManipulator Manip;
+ /// Sensitive detector
+ typedef Geometry::Handle<Geometry::SensitiveDetectorObject> SensitiveDetector;
/** @class Compare Geant4HitCollection.h DDG4/Geant4HitCollection.h
*
@@ -187,6 +196,8 @@ namespace DD4hep {
protected:
/// The collection of hit pointers in the wrapped format
WrappedHits m_hits;
+ /// Handle to the sensitive detector
+ SensitiveDetector m_detector;
/// The type of the objects in this collection. Set by the constructor
Manip* m_manipulator;
@@ -203,20 +214,26 @@ namespace DD4hep {
public:
/// Initializing constructor (C++ version)
template <typename TYPE>
- Geant4HitCollection(const std::string& det, const std::string& coll)
- : G4VHitsCollection(det, coll), m_manipulator(Geant4HitWrapper::manipulator<TYPE>()) {
+ Geant4HitCollection(const std::string& det, const std::string& coll, SensitiveDetector sd)
+ : G4VHitsCollection(det, coll), m_detector(sd),
+ m_manipulator(Geant4HitWrapper::manipulator<TYPE>()) {
newInstance();
m_hits.reserve(200);
}
/// Initializing constructor
template <typename TYPE>
- Geant4HitCollection(const std::string& det, const std::string& coll, const TYPE*)
- : G4VHitsCollection(det, coll), m_manipulator(Geant4HitWrapper::manipulator<TYPE>()) {
+ Geant4HitCollection(const std::string& det, const std::string& coll, SensitiveDetector sd, const TYPE*)
+ : G4VHitsCollection(det, coll), m_detector(sd),
+ m_manipulator(Geant4HitWrapper::manipulator<TYPE>()) {
newInstance();
m_hits.reserve(200);
}
/// Default destructor
virtual ~Geant4HitCollection();
+ /// Set the sensitive detector
+ void setSensitiveDetector(SensitiveDetector detector);
+ /// Access the sensitive detector
+ SensitiveDetector sensitiveDetector() const;
/// Type information of the object stored
const ComponentCast& type() const;
/// Type information of the vector type for extracting data
diff --git a/DDG4/include/DDG4/Geant4InteractionMerger.h b/DDG4/include/DDG4/Geant4InteractionMerger.h
new file mode 100644
index 000000000..90e2a1690
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4InteractionMerger.h
@@ -0,0 +1,42 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4INTERACTIONMERGER_H
+#define DD4HEP_DDG4_GEANT4INTERACTIONMERGER_H
+
+// Framework include files
+#include "DDG4/Geant4GeneratorAction.h"
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ /** Geant4Action to convert the particle information to Geant4
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4InteractionMerger : public Geant4GeneratorAction {
+ public:
+ /// Standard constructor
+ Geant4InteractionMerger(Geant4Context* context, const std::string& nam);
+ /// Default destructor
+ virtual ~Geant4InteractionMerger();
+ /// Event generation action callback
+ virtual void operator()(G4Event* event);
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+
+#endif // DD4HEP_DDG4_GEANT4INTERACTIONMERGER_H
diff --git a/DDG4/include/DDG4/Geant4InteractionVertexBoost.h b/DDG4/include/DDG4/Geant4InteractionVertexBoost.h
new file mode 100644
index 000000000..fd2017763
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4InteractionVertexBoost.h
@@ -0,0 +1,49 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4INTERACTIONVERTEXBOOST_H
+#define DD4HEP_DDG4_GEANT4INTERACTIONVERTEXBOOST_H
+
+// Framework include files
+#include "DDG4/Geant4GeneratorAction.h"
+
+// ROOT include files
+#include "Math/Vector4D.h"
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+ /// Generate particles isotrop in space around origine (0,0,0)
+ /**
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4InteractionVertexBoost: public Geant4GeneratorAction {
+ protected:
+ /// The constant Lorentz transformation angle
+ double m_angle;
+ /// Property: Unique identifier of the interaction to be modified
+ int m_mask;
+ public:
+ /// Standard constructor
+ Geant4InteractionVertexBoost(Geant4Context* context, const std::string& name);
+ /// Default destructor
+ virtual ~Geant4InteractionVertexBoost();
+ /// Callback to generate primary particles
+ virtual void operator()(G4Event* event);
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif /* DD4HEP_DDG4_GEANT4INTERACTIONVERTEXBOOST_H */
diff --git a/DDG4/include/DDG4/Geant4InteractionVertexSmear.h b/DDG4/include/DDG4/Geant4InteractionVertexSmear.h
new file mode 100644
index 000000000..76cebb4df
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4InteractionVertexSmear.h
@@ -0,0 +1,51 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4INTERACTIONVERTEXSMEAR_H
+#define DD4HEP_DDG4_GEANT4INTERACTIONVERTEXSMEAR_H
+
+// Framework include files
+#include "DDG4/Geant4GeneratorAction.h"
+
+// ROOT include files
+#include "Math/Vector4D.h"
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+ /// Generate particles isotrop in space around origine (0,0,0)
+ /**
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4InteractionVertexSmear: public Geant4GeneratorAction {
+ protected:
+ /// The constant smearing offset
+ ROOT::Math::PxPyPzEVector m_offset;
+ /// The gaussian sigmas to the offset
+ ROOT::Math::PxPyPzEVector m_sigma;
+ /// Property: Unique identifier of the interaction created
+ int m_mask;
+ public:
+ /// Standard constructor
+ Geant4InteractionVertexSmear(Geant4Context* context, const std::string& name);
+ /// Default destructor
+ virtual ~Geant4InteractionVertexSmear();
+ /// Callback to generate primary particles
+ virtual void operator()(G4Event* event);
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif /* DD4HEP_DDG4_GEANT4INTERACTIONVERTEXSMEAR_H */
diff --git a/DDG4/include/DDG4/Geant4IsotropeGenerator.h b/DDG4/include/DDG4/Geant4IsotropeGenerator.h
new file mode 100644
index 000000000..0fa7ba5e0
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4IsotropeGenerator.h
@@ -0,0 +1,55 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4ISOTROPEGENERATOR_H
+#define DD4HEP_DDG4_GEANT4ISOTROPEGENERATOR_H
+
+// Framework include files
+#include "DDG4/Geant4GeneratorAction.h"
+
+// Forward declarations
+class G4ParticleDefinition;
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+ /// Generate particles isotrop in space around origine (0,0,0)
+ /**
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4IsotropeGenerator: public Geant4GeneratorAction {
+ protected:
+ /// Pointer to geant4 particle definition
+ G4ParticleDefinition* m_particle;
+ /// Property: Particle energy
+ double m_energy;
+ /// Property: Particle name
+ std::string m_particleName;
+ /// Property: Desired multiplicity of the particles to be shot
+ int m_multiplicity;
+ /// Property: User mask passed to all particles in the generated interaction
+ int m_mask;
+ public:
+ /// Standard constructor
+ Geant4IsotropeGenerator(Geant4Context* context, const std::string& name);
+ /// Default destructor
+ virtual ~Geant4IsotropeGenerator();
+ /// Callback to generate primary particles
+ virtual void operator()(G4Event* event);
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif /* DD4HEP_DDG4_GEANT4ISOTROPEGENERATOR_H */
diff --git a/DDG4/include/DDG4/Geant4MonteCarloTruth.h b/DDG4/include/DDG4/Geant4MonteCarloTruth.h
index f19f64f54..7c9726b60 100644
--- a/DDG4/include/DDG4/Geant4MonteCarloTruth.h
+++ b/DDG4/include/DDG4/Geant4MonteCarloTruth.h
@@ -31,7 +31,7 @@ namespace DD4hep {
namespace Simulation {
// Forward declarations
- class Particle;
+ class Geant4Particle;
/** @class Geant4MonteCarloTruth Geant4MonteCarloTruth.h DDG4/Geant4MonteCarloTruth.h
*
@@ -42,6 +42,7 @@ namespace DD4hep {
*/
class Geant4MonteCarloTruth {
public:
+ typedef Geant4Particle Particle;
typedef std::map<int,Particle*> ParticleMap;
typedef std::map<int,int> TrackEquivalents;
protected:
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
new file mode 100644
index 000000000..6dd0d3e82
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -0,0 +1,303 @@
+// $Id: Geant4Data.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_GEANT4PARTICLE_H
+#define DD4HEP_GEANT4PARTICLE_H
+
+// Framework include files
+
+// ROOT includes
+#include "Math/Vector4D.h"
+
+// Geant4 forward declarations
+class G4ParticleDefinition;
+class G4VProcess;
+
+// C/C++ include files
+#include <set>
+#include <map>
+#include <memory>
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Particle;
+
+ class ParticleExtension {
+ public:
+ /// Default constructor
+ ParticleExtension() {}
+ /// Default destructor
+ virtual ~ParticleExtension();
+ };
+
+ /// Track properties
+ enum Geant4ParticleProperties {
+ G4PARTICLE_CREATED_HIT = 1<<1,
+ G4PARTICLE_PRIMARY = 1<<2,
+ G4PARTICLE_HAS_SECONDARIES = 1<<3,
+ G4PARTICLE_ABOVE_ENERGY_THRESHOLD = 1<<4,
+ G4PARTICLE_KEEP_PROCESS = 1<<5,
+ G4PARTICLE_KEEP_PARENT = 1<<6,
+ G4PARTICLE_CREATED_CALORIMETER_HIT = 1<<7,
+ G4PARTICLE_CREATED_TRACKER_HIT = 1<<8,
+ G4PARTICLE_KEEP_ALWAYS = 1<<9,
+
+ // Generator status for a given particles: bit 0...3 come from LCIO, rest is internal
+ G4PARTICLE_GEN_EMPTY = 1<<0, // Empty line
+ G4PARTICLE_GEN_STABLE = 1<<1, // undecayed particle, stable in the generator
+ 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_STATUS = 0x3FF, // Mask for generator status (bit 0...9)
+
+ // Simulation status of a given particle
+ G4PARTICLE_SIM_CREATED = 1<<10, // True if the particle has been created by the simulation program (rather than the generator)
+ G4PARTICLE_SIM_BACKSCATTER = 1<<11, // True if the particle is the result of a backscatter from a calorimeter shower.
+ G4PARTICLE_SIM_DECAY_CALO = 1<<12, // True if the particle has interacted in a calorimeter region.
+ G4PARTICLE_SIM_DECAY_TRACKER = 1<<13, // True if the particle has interacted in a tracking region.
+ G4PARTICLE_SIM_STOPPED = 1<<14, // True if the particle has been stopped by the simulation program.
+ G4PARTICLE_SIM_LEFT_DETECTOR = 1<<15, // True if the particle has left the world volume undecayed.
+ G4PARTICLE_SIM_PARENT_RADIATED = 1<<16, // True if the particle's vertex is not the endpoint of the parent particle.
+ G4PARTICLE_SIM_OVERLAY = 1<<17, // True if the particle has been overlayed by the simulation (or digitization) program.
+
+ G4PARTICLE_LAST_NOTHING = 1<<31
+ };
+
+
+
+ /// Data structure to store the MC particle information
+ /**
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4Particle {
+ private:
+ /// Copy constructor
+ Geant4Particle(const Geant4Particle& c);
+ public:
+ typedef std::set<int> Particles;
+ /// Reference counter
+ int ref; //! not persistent
+ int id, g4Parent, reason, usermask;
+ int steps, secondaries, pdgID;
+ int status, colorFlow[2];
+ char charge, _spare[3];
+ float spin[3];
+ // 12 ints + 4 floats should be aligned to 8 bytes....
+ double vsx, vsy, vsz;
+ double vex, vey, vez;
+ double psx, psy, psz;
+ double pex, pey, pez;
+ double mass, time, properTime;
+ /// The list of daughters of this MC particle
+ Particles parents;
+ Particles daughters;
+
+ /// User data extension if required
+ std::auto_ptr<ParticleExtension> extension;
+ const G4VProcess *process; //! not persistent
+ const G4ParticleDefinition *definition; //! not persistent
+ /// Default constructor
+ Geant4Particle();
+ /// Constructor with ID initialization
+ Geant4Particle(int part_id);
+ /// Default destructor
+ virtual ~Geant4Particle();
+ /// Increase reference count
+ Geant4Particle* addRef() {
+ ++ref;
+ return this;
+ }
+ /// Decrease reference count. Deletes object if NULL
+ void release();
+ /// Assignment operator
+ Geant4Particle& get_data(Geant4Particle& c);
+ /// Remove daughter from set
+ void removeDaughter(int id_daughter);
+ };
+
+#ifndef __DD4HEP_DDEVE_EXCLUSIVE__
+ /// Data structure to access MC particle information
+ /**
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4ParticleHandle {
+ protected:
+ /// Particle pointer
+ Geant4Particle* particle;
+ public:
+ /// Default constructor
+ Geant4ParticleHandle(Geant4Particle* part);
+ /// Initializing constructor
+ Geant4ParticleHandle(const Geant4ParticleHandle& h);
+ /// Assignment operator
+ Geant4ParticleHandle& operator=(Geant4Particle* part);
+ /// Overloaded -> operator to access particle details
+ Geant4Particle* operator->() const;
+ /// Conversion operator to pointer
+ operator Geant4Particle*() const;
+ /// Accessor to the number of particle parents
+ size_t numParent() const;
+ /// Accessor to the number of particle daughters
+ size_t numDaughter() const;
+ /// Scalar particle momentum squared
+ double momentum2() const;
+ /// Scalar particle energy
+ double energy() const;
+ /// Scalar particle momentum
+ double momentum() const { return sqrt(momentum2()); }
+ /// Access to the Geant4 particle name
+ std::string particleName() const;
+ /// Access to the Geant4 particle type
+ std::string particleType() const;
+ /// Access to the creator process name
+ std::string processName() const;
+ /// Access to the creator process type name
+ std::string processTypeName() const;
+ /// Access patricle momentum, energy as 4 vector
+ ROOT::Math::PxPyPzM4D<double> pxPyPzM() const;
+ /// Access patricle momentum, energy as 4 vector
+ ROOT::Math::Cartesian3D<double> startVertex() const;
+ /// Access patricle momentum, energy as 4 vector
+ ROOT::Math::Cartesian3D<double> endVertex() const;
+
+ /// Various output formats:
+
+ /// Output type 1:+++ <tag> 10 def:0xde4eaa8 [gamma , gamma] reason: 20 E:+1.017927e+03 #Par: 1/4 #Dau: 2
+ void dump1(int level, const std::string& src, const char* tag) const;
+ /// Output type 2:+++ <tag> 20 G4: 7 def:0xde4eaa8 [gamma , gamma] reason: 20 E:+3.304035e+01 in record:YES #Par: 1/18 #Dau: 0
+ void dump2(int level, const std::string& src, const char* tag, int g4id, bool inrec) 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
+ void dump3(int level, const std::string& src, const char* tag) const;
+
+ /// Handlers
+
+ /// Offset all particle identifiers (id, parents, daughters) by a constant number
+ void offset(int off) const;
+
+ };
+
+ inline Geant4ParticleHandle::Geant4ParticleHandle(const Geant4ParticleHandle& c)
+ : particle(c.particle) {
+ }
+
+ /// Initializing constructor
+ inline Geant4ParticleHandle::Geant4ParticleHandle(Geant4Particle* p)
+ : particle(p) {
+ }
+
+ /// Overloaded -> operator to access particle details
+ inline Geant4Particle* Geant4ParticleHandle::operator->() const {
+ return particle;
+ }
+ /// Assignment operator
+ inline Geant4ParticleHandle& Geant4ParticleHandle::operator=(Geant4Particle* part) {
+ particle = part;
+ return *this;
+ }
+ /// Conversion operator to pointer
+ inline Geant4ParticleHandle::operator Geant4Particle*() const {
+ return particle;
+ }
+ /// Accessor to the number of particle parents
+ inline size_t Geant4ParticleHandle::numParent() const {
+ return particle->parents.size();
+ }
+ /// Accessor to the number of particle daughters
+ inline size_t Geant4ParticleHandle::numDaughter() const {
+ return particle->daughters.size();
+ }
+ /// Access patricle momentum, energy as 4 vector
+ inline ROOT::Math::PxPyPzM4D<double> Geant4ParticleHandle::pxPyPzM() const {
+ const Geant4Particle* p = particle;
+ return ROOT::Math::PxPyPzM4D<double>(p->psx,p->psy,p->psz,p->mass);
+ }
+
+ /// Access patricle momentum, energy as 4 vector
+ inline ROOT::Math::Cartesian3D<double> Geant4ParticleHandle::startVertex() const {
+ const Geant4Particle* p = particle;
+ return ROOT::Math::Cartesian3D<double>(p->vsx,p->vsy,p->vsz);
+ }
+
+ /// Access patricle momentum, energy as 4 vector
+ inline ROOT::Math::Cartesian3D<double> Geant4ParticleHandle::endVertex() const {
+ const Geant4Particle* p = particle;
+ return ROOT::Math::Cartesian3D<double>(p->vex,p->vey,p->vez);
+ }
+
+ /// Scalar particle energy
+ inline double Geant4ParticleHandle::energy() const {
+ const Geant4Particle* p = particle;
+ ROOT::Math::PxPyPzM4D<double> v(p->psx,p->psy,p->psz,p->mass);
+ return v.E();
+ }
+
+ /// Scalar particle momentum squared
+ inline double Geant4ParticleHandle::momentum2() const {
+ const Geant4Particle* p = particle;
+ return (p->psx*p->psx + p->psy*p->psy + p->psz*p->psz);
+ }
+
+ /// Data structure to map particles produced during the generation and the simulation
+ /**
+ * This data structure is added to the Geant4Event data extensions
+ * by the Geant4GenerationInit action.
+ * Particles are added:
+ * - During the generation if the required modules are activated
+ * - At the end of the handling of the MC truth are particles to be kept
+ * are inserted if the required modules are activated such as the
+ * Geant4ParticleHandler.
+ *
+ * Note: This object takes OWNERSHIP of the inserted particles!
+ * beware of double deletion of objects!
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4ParticleMap {
+ public:
+ typedef Geant4Particle Particle;
+ typedef std::map<int,Particle*> ParticleMap;
+ typedef std::map<int,int> TrackEquivalents;
+ /// Mapping of particles of this event
+ ParticleMap particleMap; //! not persistent
+ /// Map associating the G4Track identifiers with identifiers of existing MCParticles
+ TrackEquivalents equivalentTracks;
+
+ /// Default constructor
+ Geant4ParticleMap() {}
+ /// Default destructor
+ virtual ~Geant4ParticleMap();
+ /// Clear particle maps
+ void clear();
+ /// Adopt particle maps
+ void adopt(ParticleMap& pm, TrackEquivalents& equiv);
+ /// Access the particle map
+ const ParticleMap& particles() const { return particleMap; }
+ /// Access the map of track equivalents
+ const TrackEquivalents& equivalents() const { return equivalentTracks; }
+ /// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
+ int particleID(int track, bool throw_if_not_found=true) const;
+ };
+#endif
+
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif // DD4HEP_GEANT4PARTICLE_H
diff --git a/DDG4/include/DDG4/Geant4ParticleGun.h b/DDG4/include/DDG4/Geant4ParticleGun.h
index 69bb5cd2a..52fb2ec75 100644
--- a/DDG4/include/DDG4/Geant4ParticleGun.h
+++ b/DDG4/include/DDG4/Geant4ParticleGun.h
@@ -1,4 +1,3 @@
-// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
//====================================================================
// AIDA Detector description implementation
//--------------------------------------------------------------------
@@ -16,7 +15,6 @@
// Forward declarations
class G4ParticleDefinition;
class G4ParticleGun;
-class TRandom1;
/*
* DD4hep namespace declaration
@@ -42,8 +40,6 @@ namespace DD4hep {
G4ParticleDefinition* m_particle;
/// Pointer to the particle gun itself
G4ParticleGun* m_gun;
- /// Random number generator
- TRandom1* m_rndm;
/// Particle energy
double m_energy;
/// Particle name
@@ -65,4 +61,3 @@ namespace DD4hep {
} // End namespace Simulation
} // End namespace DD4hep
#endif /* DD4HEP_DDG4_GEANT4PARTICLEGUN_H */
-
diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h
index 35f8038af..bafcea190 100644
--- a/DDG4/include/DDG4/Geant4ParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4ParticleHandler.h
@@ -10,7 +10,7 @@
#define DD4HEP_DDG4_GEANT4PARTICLEHANDLER_H
// Framework include files
-#include "DDG4/Geant4Data.h"
+#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4GeneratorAction.h"
#include "DDG4/Geant4MonteCarloTruth.h"
@@ -32,18 +32,26 @@ namespace DD4hep {
// Forward declarations
class Particle;
+ class Geant4PrimaryMap;
class Geant4UserParticleHandler;
/** Geant4Action to collect the MC particle information.
*
+ * Extract the relevant particle information during the simulation step.
*
- * @author M.Frank
- * @version 1.0
+ * @author M.Frank
+ * @version 1.0
*/
- class Geant4ParticleHandler : public Geant4GeneratorAction, public Geant4MonteCarloTruth
- {
+ class Geant4ParticleHandler : public Geant4GeneratorAction, public Geant4MonteCarloTruth {
public:
typedef std::vector<std::string> Processes;
+ struct FindParticleByID {
+ int pid;
+ FindParticleByID(int p) : pid(p) {}
+ inline bool operator()(const std::pair<int,Geant4Particle*>& p) const {
+ return p.second->id == pid;
+ }
+ };
protected:
/// Property: Steer printout at tracking action begin
@@ -52,10 +60,17 @@ namespace DD4hep {
bool m_printEndTracking;
/// Property: Flag to keep all particles generated
bool m_keepAll;
+ /// Property: Flag if the handler is executed in standalone mode and hence must manage particles
+ bool m_ownsParticles;
/// Property: Energy cut below which particles are not collected, but assigned to the parent
double m_kinEnergyCut;
+
+ /// Global particle identifier. Obtained at the begin of the event.
+ int m_globalParticleID, m_initialParticleID;
/// User action pointer
Geant4UserParticleHandler* m_userHandler;
+ /// Primary map
+ Geant4PrimaryMap* m_primaryMap;
/// Local buffer about the 'current' G4Track
Particle m_currTrack;
/// Map with stored MC Particles
@@ -74,6 +89,11 @@ namespace DD4hep {
/// Get proper equivalent track from the particle map according to the given geant4 track ID
int equivalentTrack(int g4_id) const;
+ /// Rebase the simulated tracks, so that they fit to the generator particles
+ void rebaseSimulatedTracks(int base);
+ /// Debugging: Dump Geant4 particle map
+ void dumpMap(const char* tag) const;
+
public:
/// Standard constructor
Geant4ParticleHandler(Geant4Context* context, const std::string& nam);
diff --git a/DDG4/include/DDG4/Geant4ParticlePrint.h b/DDG4/include/DDG4/Geant4ParticlePrint.h
index 1c5e2786e..4429e408b 100644
--- a/DDG4/include/DDG4/Geant4ParticlePrint.h
+++ b/DDG4/include/DDG4/Geant4ParticlePrint.h
@@ -11,7 +11,8 @@
// Framework include files
#include "DDG4/Geant4EventAction.h"
-#include "DDG4/Geant4MonteCarloTruth.h"
+#include "DDG4/Geant4GeneratorAction.h"
+#include "DDG4/Geant4Particle.h"
/*
* DD4hep namespace declaration
@@ -29,21 +30,30 @@ namespace DD4hep {
* @author M.Frank
* @version 1.0
*/
- class Geant4ParticlePrint : public Geant4EventAction {
+ class Geant4ParticlePrint : public Geant4EventAction {
public:
- typedef Geant4MonteCarloTruth::ParticleMap ParticleMap;
- typedef Geant4MonteCarloTruth::TrackEquivalents TrackEquivalents;
+ typedef Geant4ParticleMap::Particle Particle;
+ typedef Geant4ParticleMap::ParticleMap ParticleMap;
+ typedef Geant4ParticleMap::TrackEquivalents TrackEquivalents;
protected:
/// Property: Flag to indicate output type: 1: TABLE, 2:TREE, 3:BOTH (default)
int m_outputType;
+ /// Property: Flag to indicate output type at begin of event
+ bool m_printBegin;
+ /// Property: Flag to indicate output type at end of event
+ bool m_printEnd;
+ /// Property: Flag to indicate output type as part of the generator action
+ bool m_printGeneration;
- void printParticle(const std::string& prefix, const Particle* p) const;
+ void printParticle(const std::string& prefix, Geant4ParticleHandle p) const;
/// Print record of kept particles
void printParticles(const ParticleMap& particles) const;
/// Print tree of kept particles
- void printParticleTree(const ParticleMap& particles, int level, const Particle* p) const;
+ void printParticleTree(const ParticleMap& particles, int level, Geant4ParticleHandle p) const;
/// Print tree of kept particles
void printParticleTree(const ParticleMap& particles) const;
+ /// Print particle table
+ void makePrintout() const;
public:
@@ -55,6 +65,9 @@ namespace DD4hep {
virtual void begin(const G4Event* event);
/// Post-event action callback
virtual void end(const G4Event* event);
+ /// Generation action callback
+ virtual void operator()(G4Event* event);
+
};
} // End namespace Simulation
} // End namespace DD4hep
diff --git a/DDG4/include/DDG4/Geant4Primary.h b/DDG4/include/DDG4/Geant4Primary.h
new file mode 100644
index 000000000..3c3c39c55
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4Primary.h
@@ -0,0 +1,150 @@
+// $Id: Geant4Data.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_GEANT4PRIMARY_H
+#define DD4HEP_GEANT4PRIMARY_H
+
+// Framework include files
+#include "DDG4/Geant4Vertex.h"
+#include "DDG4/Geant4Particle.h"
+
+// C/C++ include files
+#include <set>
+#include <map>
+#include <vector>
+#include <memory>
+
+// Forward declarations
+class G4PrimaryParticle;
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Particle;
+ class Geant4Vertex;
+
+ class PrimaryExtension {
+ public:
+ /// Default constructor
+ PrimaryExtension() {}
+ /// Default destructor
+ virtual ~PrimaryExtension();
+ };
+
+ /// Data structure to map primaries to particles.
+ /**
+ * This data structure is added to the Geant4Event data extensions
+ * by the Geant4GenerationInit action.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4PrimaryMap {
+ public:
+ typedef std::map<const G4PrimaryParticle*,Geant4Particle*> Primaries;
+ /// Mapping of primary particles of this event
+ Primaries primaryMap; //! not persistent
+ /// Default constructor
+ Geant4PrimaryMap() {}
+ /// Default destructor
+ virtual ~Geant4PrimaryMap();
+ };
+
+
+ /// Class modelling a single interaction with multiple primary vertices and particles
+ /*
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4PrimaryInteraction {
+ private:
+ /// Copy constructor
+ Geant4PrimaryInteraction(const Geant4PrimaryInteraction& c);
+ /// Assignment operator
+ Geant4PrimaryInteraction& operator=(const Geant4PrimaryInteraction& c);
+
+ public:
+ typedef Geant4Particle Particle;
+ typedef Geant4Vertex Vertex;
+ typedef std::map<int,Particle*> ParticleMap;
+ typedef std::map<int,Vertex*> VertexMap;
+ typedef std::auto_ptr<PrimaryExtension> ExtensionHandle;
+
+ /// The map of primary vertices for the particles.
+ VertexMap vertices;
+ /// The map of particles participating in this primary interaction
+ ParticleMap particles;
+ /// User data extension if required
+ ExtensionHandle extension;
+ /// User mask to flag the interaction. Also unique identifier
+ int mask;
+ /// Next PID indentifier
+ int next_particle_identifier;
+
+ public:
+ /// Default constructor
+ Geant4PrimaryInteraction();
+ /// Default destructor
+ virtual ~Geant4PrimaryInteraction();
+ /// Access a new particle identifier within the interaction
+ int nextPID();
+ /// Set the next PID value
+ void setNextPID(int value);
+ };
+
+ /// Class modelling a complete primary event with multiple interactions
+ /*
+ * Multiple interactions allow a simple handling of overlay events
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4PrimaryEvent {
+ private:
+ /// Copy constructor
+ Geant4PrimaryEvent(const Geant4PrimaryEvent& c);
+ /// Assignment operator
+ Geant4PrimaryEvent& operator=(const Geant4PrimaryEvent& c);
+
+ public:
+ typedef Geant4PrimaryInteraction Interaction;
+ typedef std::map<int,Interaction*> Interactions;
+ typedef std::auto_ptr<PrimaryExtension> ExtensionHandle;
+
+ protected:
+ /// Set of primary interactions
+ Interactions m_interactions;
+
+ public:
+ /// User data extension if required
+ ExtensionHandle extension;
+
+ public:
+ /// Default constructor
+ Geant4PrimaryEvent();
+ /// Default destructor
+ virtual ~Geant4PrimaryEvent();
+ /// Add a new interaction object to the event
+ void add(int id, Geant4PrimaryInteraction* interaction);
+ /// Retrieve an interaction by it's ID
+ Interaction* get(int id) const;
+ /// Retrieve all intractions
+ std::vector<Interaction*> interactions() const;
+ };
+
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif // DD4HEP_GEANT4PRIMARY_H
diff --git a/DDG4/include/DDG4/Geant4PrimaryHandler.h b/DDG4/include/DDG4/Geant4PrimaryHandler.h
new file mode 100644
index 000000000..0b6873257
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4PrimaryHandler.h
@@ -0,0 +1,42 @@
+// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4PRIMARYHANDLER_H
+#define DD4HEP_DDG4_GEANT4PRIMARYHANDLER_H
+
+// Framework include files
+#include "DDG4/Geant4GeneratorAction.h"
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ /** Geant4Action to convert the particle information to Geant4
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4PrimaryHandler : public Geant4GeneratorAction {
+ public:
+ /// Standard constructor
+ Geant4PrimaryHandler(Geant4Context* context, const std::string& nam);
+ /// Default destructor
+ virtual ~Geant4PrimaryHandler();
+ /// Event generation action callback
+ virtual void operator()(G4Event* event);
+ };
+ } // End namespace Simulation
+} // End namespace DD4hep
+
+#endif // DD4HEP_DDG4_GEANT4PRIMARYHANDLER_H
diff --git a/DDG4/include/DDG4/Geant4Random.h b/DDG4/include/DDG4/Geant4Random.h
new file mode 100644
index 000000000..29b5d3dc8
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4Random.h
@@ -0,0 +1,69 @@
+// $Id: Geant4Data.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_GEANT4RANDOM_H
+#define DD4HEP_GEANT4RANDOM_H
+
+// Framework include files
+
+// C/C++ include files
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Exec;
+
+ /// Mini interface to THE random generator of the application
+ /**
+ * Mini interface to THE random generator of the application.
+ * Necessary, that on every object creates its own instance, but accesses
+ * the main instance avaible throu the Geant4Context.
+ *
+ * This is mandatory to ensure reproducability of the event generation
+ * process. Particular objects may use a dependent generator from
+ * an experiment framework like GAUDI.
+ *
+ * This main interface is supposed to be stable. Unclear however is
+ * if the generation functions will have to become virtual....
+ * Future will tell us.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4Random {
+ friend class Geant4Exec;
+
+ protected:
+ /// Default constructor
+ Geant4Random();
+ /// Default destructor
+ virtual ~Geant4Random();
+ public:
+ void circle(double &x, double &y, double r);
+ double exp(double tau);
+ double gauss(double mean=0, double sigma=1);
+ double landau(double mean=0, double sigma=1);
+ double rndm(int i=0);
+ void rndmArray(int n, float *array);
+ void rndmArray(int n, double *array);
+ void sphere(double &x, double &y, double &z, double r);
+ double uniform(double x1=1);
+ double uniform(double x1, double x2);
+ };
+
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif // DD4HEP_GEANT4RANDOM_H
diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h
index 77674a7c1..9b720a68a 100644
--- a/DDG4/include/DDG4/Geant4SensDetAction.h
+++ b/DDG4/include/DDG4/Geant4SensDetAction.h
@@ -231,7 +231,9 @@ namespace DD4hep {
*/
class Geant4SensDetActionSequence: public Geant4Action {
public:
- typedef G4VHitsCollection* (*create_t)(const std::string&, const std::string&);
+
+ typedef Geometry::SensitiveDetector SensitiveDetector;
+ typedef Geant4HitCollection* (*create_t)(const std::string&, const std::string&, SensitiveDetector);
typedef std::vector<std::pair<std::string, create_t> > HitCollections;
protected:
@@ -248,15 +250,18 @@ namespace DD4hep {
/// The list of sensitive detector objects
Actors<Geant4Sensitive> m_actors;
/// The list of sensitive detector filter objects
- Actors<Geant4Filter> m_filters;
+ Actors<Geant4Filter> m_filters;
/// Hit collection creators
HitCollections m_collections;
+ /// Reference to the sensitive detector element
+ SensitiveDetector m_sensitive;
/// Reference to G4 sensitive detector
Geant4ActionSD* m_detector;
+
/// Create a new typed hit collection
- template <typename TYPE> static G4VHitsCollection* _create(const std::string& det, const std::string& coll) {
- return new Geant4HitCollection(det, coll, (TYPE*) 0);
+ template <typename TYPE> static Geant4HitCollection* _create(const std::string& det, const std::string& coll, SensitiveDetector sd) {
+ return new Geant4HitCollection(det, coll, sd, (TYPE*) 0);
}
public:
diff --git a/DDG4/include/DDG4/Geant4UserParticleHandler.h b/DDG4/include/DDG4/Geant4UserParticleHandler.h
index 5047595a0..745cb3de8 100644
--- a/DDG4/include/DDG4/Geant4UserParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4UserParticleHandler.h
@@ -30,7 +30,7 @@ namespace DD4hep {
namespace Simulation {
// Forward declarations
- class Particle;
+ class Geant4Particle;
class Geant4ParticleHandler;
/** Geant4ParticleHandler user extension action called by the particle handler.
@@ -45,6 +45,8 @@ namespace DD4hep {
* @version 1.0
*/
class Geant4UserParticleHandler : public Geant4Action {
+ public:
+ typedef Geant4Particle Particle;
public:
/// Standard constructor
Geant4UserParticleHandler(Geant4Context* context, const std::string& nam);
diff --git a/DDG4/include/DDG4/Geant4Vertex.h b/DDG4/include/DDG4/Geant4Vertex.h
new file mode 100644
index 000000000..3ad19da25
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4Vertex.h
@@ -0,0 +1,72 @@
+// $Id: Geant4Data.h 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_GEANT4VERTEX_H
+#define DD4HEP_GEANT4VERTEX_H
+
+// C/C++ include files
+#include <set>
+#include <memory>
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Vertex;
+
+ class VertexExtension {
+ public:
+ /// Default constructor
+ VertexExtension() {}
+ /// Default destructor
+ virtual ~VertexExtension();
+ };
+
+ /// Data structure to store the MC vertex information
+ /**
+ * @author M.Frank
+ * @version 1.0
+ */
+ class Geant4Vertex {
+ public:
+ typedef std::set<int> Particles;
+ /// Reference counter
+ int ref; //! not persistent
+ /// Vertex data
+ double x,y,z,time;
+ /// The list of outgoing particles
+ Particles out;
+ /// The list of incoming particles
+ Particles in;
+ /// User data extension if required
+ std::auto_ptr<VertexExtension> extension;
+
+ /// Default constructor
+ Geant4Vertex();
+ /// Copy constructor
+ Geant4Vertex(const Geant4Vertex& c);
+ /// Default destructor
+ virtual ~Geant4Vertex();
+ /// Assignment operator
+ Geant4Vertex& operator=(const Geant4Vertex& c);
+ /// Increase reference count
+ Geant4Vertex* addRef();
+ /// Decrease reference count. Deletes object if NULL
+ void release();
+ };
+
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif // DD4HEP_GEANT4VERTEX_H
diff --git a/DDG4/lcio/EventReader.cpp b/DDG4/lcio/EventReader.cpp
deleted file mode 100644
index bca719711..000000000
--- a/DDG4/lcio/EventReader.cpp
+++ /dev/null
@@ -1,303 +0,0 @@
-// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
-//====================================================================
-// AIDA Detector description implementation for LCD
-//--------------------------------------------------------------------
-//
-// @author P.Kostka (main author)
-// @author M.Frank (code reshuffeling into new DDG4 scheme)
-//
-//====================================================================
-
-// Framework include files
-#include "EventReader.h"
-#include "DD4hep/Plugins.h"
-#include "G4Event.hh"
-#include "G4PrimaryParticle.hh"
-#include "G4ParticleTable.hh"
-#include "IMPL/MCParticleImpl.h"
-#include "IMPL/LCCollectionVec.h"
-#include "Randomize.hh"
-
-using namespace std;
-
-/// Initializing constructor
-DD4hep::lcio::EventReader::EventReader(const string& nam) : m_name(nam) {
-}
-
-/// Default destructor
-DD4hep::lcio::EventReader::~EventReader() {
-}
-
-static void setCharge(IMPL::MCParticleImpl* p) {
- // Initialise the particle charge looking at the actual
- // Particle Table, or set it to -1000 to flag missing
- // information, if PDG unknown in the actual physics list..
- G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
- G4ParticleDefinition* def = tab->FindParticle(p->getPDG());
- double charge = def ? def->GetPDGCharge() : -1000;
- p->setCharge(charge);
-}
-
-/// Standard constructor
-DD4hep::lcio::LcioGeneratorAction::LcioGeneratorAction(Simulation::Geant4Context* context, const string& name)
- : Simulation::Geant4GeneratorAction(context,name)
-{
- declareProperty("Input", m_input);
- declareProperty("Sync", m_SYNCEVT=0);
- declareProperty("zSpread", m_primaryVertexSpreadZ=0.0);
- declareProperty("lorentzAngle",m_lorentzTransformationAngle=0.0);
-}
-
-/// Default destructor
-DD4hep::lcio::LcioGeneratorAction::~LcioGeneratorAction() {
-}
-
-static string issue(int i) {
- stringstream str;
- str << "lcio::LcioGeneratorAction: Event " << i;
- return str.str();
-}
-
-/// Read an event and return a LCCollectionVec of MCParticles.
-DD4hep::lcio::LcioGeneratorAction::Particles *DD4hep::lcio::LcioGeneratorAction::readEvent(int which) {
- int evid = which + m_SYNCEVT;
- if ( 0 == m_reader ) {
- if ( m_input.empty() ) {
- throw runtime_error("LcioGeneratorAction: No inoput file declared!");
- }
- string err;
- Simulation::TypeName tn = Simulation::TypeName::split(m_input,"|");
- try {
- m_reader = PluginService::Create<DD4hep::lcio::EventReader*>(tn.first,tn.second,int(0));
- if ( 0 == m_reader ) {
- abortRun(issue(evid),"Error creating reader plugin.",
- "Failed to create file reader of type %s. Cannot open dataset %s",
- tn.first.c_str(),tn.second.c_str());
- return 0;
- }
- }
- catch(const std::exception& e) {
- err = e.what();
- }
- if ( !err.empty() ) {
- abortRun(issue(evid),err,"Error when reading file %s",m_input.c_str());
- return 0;
- }
- }
- Particles* parts = m_reader->readEvent();
- if ( 0 == parts ) {
- abortRun(issue(evid),"Error when reading file - may be end of file.",
- "Error when reading file %s",m_input.c_str());
- }
- return parts;
-}
-
-/// Callback to generate primary particles
-void DD4hep::lcio::LcioGeneratorAction::operator()(G4Event* evt) {
- //**********************************************************
- // Read in the event
- //**********************************************************
- m_g4ParticleMap.clear();
- m_convertedParticles.clear();
- Particles* primaries = readEvent(evt->GetEventID());
- if ( 0 == primaries ) return;
-
- vector<MCParticle*> col(primaries->size()) ;
- int NHEP = col.size();
-
- // parameters of the Lorentz transformation matrix
- const double zspreadparameter = m_primaryVertexSpreadZ;
- const double alpha = m_lorentzTransformationAngle;
- const double gamma = sqrt(1 + sqr(tan(alpha)));
- const double betagamma = tan(alpha);
-
- double zspread = ( zspreadparameter == 0.0 ? 0.0 : G4RandGauss::shoot(0,zspreadparameter/mm ) );
- //cout << " PrimaryVertexSpreadZ=" << zspreadparameter/mm <<"mm, zspread=" << zspread << "mm" << endl;
- G4ThreeVector particle_position = G4ThreeVector(0,0,zspread);
-
- // build collection of MCParticles
- for(int IHEP=0; IHEP<NHEP; IHEP++ ) {
- Particle* mcp = dynamic_cast<Particle*> ( primaries->getElementAt( IHEP ) );
- setCharge(mcp);
- //Boost to lab frame with crossing angle alpha
- if (alpha != 0) {
-
- double c_light = 299.792;
- const double* p = mcp->getMomentum();
- const double m = mcp->getMass();
- //double e = mcp->getEnergy();
- // after the transformation (boost in x-direction)
- double pPrime[3];
- pPrime[0] = betagamma * sqrt(sqr(p[0])+sqr(p[1])+sqr(p[2])+sqr(m)) + gamma * p[0];
- pPrime[1] = p[1];
- pPrime[2] = p[2];
-
- const double* v = mcp->getVertex();
- const double t = mcp->getTime();
- double vPrime[3];
- double tPrime;
- tPrime = gamma * t + betagamma * v[0] / c_light;
- vPrime[0] = gamma * v[0] + betagamma * c_light * t;
- vPrime[1] = v[1];
- vPrime[2] = v[2];
-
- // py and pz remain the same, E changes implicitly with px
- mcp->setMomentum( pPrime );
- mcp->setVertex( vPrime );
- mcp->setTime( tPrime );
- //cout << IHEP << " alpha=" << alpha <<" gamma=" << gamma << " betagamma=" << betagamma << endl;
- // cout << IHEP << " unboosted momentum: (" << e << "," << p[0] << "," << p[1] << "," << p[2] << ")" << endl;
- // cout << IHEP << " boosted momentum: (" << mcp->getEnergy() << "," << betagamma * sqrt(sqr(p[0])+sqr(p[1])+sqr(p[2])+sqr(m)) + gamma * p[0] << "," << pPrime[1] << "," << pPrime[2] << ")" << endl;
- // cout << IHEP << " (" << mcp->getEnergy()-e << "," << ( betagamma * sqrt(sqr(p[0])+sqr(p[1])+sqr(p[2])+sqr(m)) + gamma * p[0])-p[0] << "," << pPrime[1]-p[1] << "," << pPrime[2]-p[2] << ")" << endl;
- //cout << IHEP << " unboosted vertex: (" << t << "," << v[0] << "," << v[1] << "," << v[2] << ")" << endl;
- // cout << IHEP << " boosted vertex: (" << tPrime << "," << vPrime[0] << "," << vPrime[1] << "," << vPrime[2] << ")" << endl;
- // cout << IHEP << " boosted vertex: (" << tPrime-t << "," << vPrime[0]-v[0] << "," << vPrime[1]-v[1] << "," << vPrime[2]-v[2] << ")" << endl;
- }
- // apply z spread
- if (zspreadparameter != 0){
- const double* v = mcp->getVertex();
- double vPrime[3];
- vPrime[0] = v[0];
- vPrime[1] = v[1];
- vPrime[2] = v[2] + zspread;
- mcp->setVertex(vPrime);
- }
- col[IHEP] = dynamic_cast<MCParticle*>( mcp );
- }
-
- // check if there is at least one particle
- if( NHEP == 0 ) return;
-
- // create G4PrimaryVertex object
- double particle_time = 0.0;
- G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position,particle_time);
- // put initial particles to the vertex
- for( size_t i=0; i< col.size(); i++ ){
- MCParticle* mcp = col[i];
- if ( mcp->getParents().size()==0 ) {
- G4PrimarySet g4set = getRelevantParticles(mcp);
- for (G4PrimarySet::iterator setit=g4set.begin(); setit != g4set.end(); setit++ ){
- vertex->SetPrimary(*setit);
- //cout << "G4PrimaryParticle ("<< (*setit)->GetPDGcode() << ") added to G4PrimaryVertex." << endl;
- }
- }
- }
- // Put the vertex to G4Event object
- evt->AddPrimaryVertex(vertex);
-}
-
-DD4hep::lcio::LcioGeneratorAction::G4PrimarySet
-DD4hep::lcio::LcioGeneratorAction::getRelevantParticles(MCParticle* p) {
- // log each particle which has been called, to avoid double counting and increase efficiency
- m_convertedParticles.insert(p);
- G4ParticleMap::iterator mcpIT;
- G4PrimarySet relevantParticlesSet; //holds all relevant decay particles of p
-
- // logic starts here:
- // 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 (p->getGeneratorStatus() == 1) {
- G4PrimaryParticle* g4p;
- mcpIT = m_g4ParticleMap.find( p );
- if( mcpIT != m_g4ParticleMap.end() ){
- g4p = mcpIT->second;
- }
- else {
- int IDHEP = p->getPDG();
- double PHEP1 = p->getMomentum()[0];
- double PHEP2 = p->getMomentum()[1];
- double PHEP3 = p->getMomentum()[2];
- double PHEP5 = p->getMass();
- // create G4PrimaryParticle object
- g4p = new G4PrimaryParticle(IDHEP, PHEP1*GeV, PHEP2*GeV, PHEP3*GeV);
- g4p->SetMass(PHEP5*GeV);
- m_g4ParticleMap[ p ] = g4p;
- }
- relevantParticlesSet.insert(g4p);
- return relevantParticlesSet;
- }
-
- // 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.
- else if( p->getDaughters().size() > 0 ) { //fg: require that there is at least one daughter - otherwise forget the particle
- cout << "case 2" << endl;
- // calculate proper time
- MCParticle* dp = p->getDaughters()[0];
- double proper_time = fabs((dp->getTime()-p->getTime()) * p->getMass()) / p->getEnergy();
- // fix by S.Morozov for real != 0
- double aPrecision = dp->getTime() * p->getMass() / p->getEnergy();
- double bPrecision = p->getTime() * p->getMass() / p->getEnergy();
- double ten = 10;
- double proper_time_Precision = pow(ten,-DBL_DIG)*fmax(fabs(aPrecision),fabs(bPrecision));
-
- bool isProperTimeZero = ( proper_time <= proper_time_Precision );
-
- // -- remove original --- if (proper_time != 0) {
- if ( isProperTimeZero == false ) {
- G4PrimaryParticle* g4p;
- mcpIT = m_g4ParticleMap.find( p );
- if( mcpIT != m_g4ParticleMap.end() ){
- g4p = mcpIT->second;
- }
- else {
- int IDHEP = p->getPDG();
- double PHEP1 = p->getMomentum()[0];
- double PHEP2 = p->getMomentum()[1];
- double PHEP3 = p->getMomentum()[2];
- double PHEP5 = p->getMass();
- // create G4PrimaryParticle object
- g4p = new G4PrimaryParticle( IDHEP, PHEP1*GeV, PHEP2*GeV, PHEP3*GeV );
- g4p->SetMass( PHEP5*GeV );
- g4p->SetProperTime( proper_time*ns );
- m_g4ParticleMap[ p ] = g4p;
- G4PrimarySet vec3;
- for (size_t i=0; i<p->getDaughters().size(); i++){
- if (m_convertedParticles.count(p->getDaughters()[i])==0) {
- G4PrimarySet vec2 = getRelevantParticles(p->getDaughters()[i]);
- G4PrimarySet::iterator setit;
- for (setit=vec2.begin(); setit != vec2.end(); setit++ ){
- vec3.insert(*setit);
- }
- }
- }
- G4PrimarySet::iterator setit;
- for (setit=vec3.begin(); setit != vec3.end(); setit++ ){
- g4p->SetDaughter(*setit);
- }
- }
- relevantParticlesSet.insert(g4p);
- return relevantParticlesSet;
- }
-
- // 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).
- else {
- for (size_t i=0; i<p->getDaughters().size(); i++){
- if (m_convertedParticles.count(p->getDaughters()[i])==0){
- G4PrimarySet vec2 = getRelevantParticles(p->getDaughters()[i]);
- G4PrimarySet::iterator setit;
- for (setit=vec2.begin(); setit != vec2.end(); setit++ ){
- relevantParticlesSet.insert(*setit);
- }
- }
- }
- return relevantParticlesSet;
- }
- }
- //fg: add a default return
- return relevantParticlesSet;
-}
-
-#include "DDG4/Factories.h"
-DECLARE_GEANT4ACTION_NS(DD4hep::lcio,LcioGeneratorAction)
diff --git a/DDG4/lcio/EventReader.h b/DDG4/lcio/EventReader.h
deleted file mode 100644
index 00219a933..000000000
--- a/DDG4/lcio/EventReader.h
+++ /dev/null
@@ -1,125 +0,0 @@
-// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
-//====================================================================
-// AIDA Detector description implementation for LCD
-//--------------------------------------------------------------------
-//
-//====================================================================
-#ifndef DD4HEP_DDG4_LCIO_EVENTREADER_H
-#define DD4HEP_DDG4_LCIO_EVENTREADER_H
-
-// Framework include files
-#include "DD4hep/Primitives.h"
-#include "DDG4/Geant4GeneratorAction.h"
-
-// C/C++ include files
-#include <set>
-#include <map>
-
-// Forward declarations
-class G4Event;
-class G4PrimaryParticle;
-namespace IO { class LCReader; }
-namespace UTIL { class LCStdHepRdr; }
-namespace EVENT { class MCParticle; }
-namespace IMPL { class MCParticleImpl; }
-namespace IMPL { class LCCollectionVec; }
-
-
-/*
- * DD4hep namespace declaration
- */
-namespace DD4hep {
- /*
- * lcio namespace declaration
- */
- namespace lcio {
-
- /** @class EventReader EventReader.h DDG4/EventReader.h
- *
- * Base class to read lcio files
- *
- * @author P.Kostka (main author)
- * @author M.Frank (code reshuffeling into new DDG4 scheme)
- * @version 1.0
- */
- class EventReader {
- public:
- typedef IMPL::LCCollectionVec Particles;
- protected:
- /// File name
- std::string m_name;
- public:
- /// Initializing constructor
- EventReader(const std::string& nam);
- /// Default destructor
- virtual ~EventReader();
- /// File name
- const std::string& name() const { return m_name; }
- /// Read an event and return a LCCollectionVec of MCParticles.
- virtual Particles *readEvent() = 0;
- };
-
- /** @class LcioGeneratorAction Geant4GeneratorAction.h DDG4/Geant4GeneratorAction.h
- *
- * Concrete implementation of the Geant4 generator action base class
- * populating Geant4 primaries from LCIO and HepStd files.
- *
- * @author P.Kostka (main author)
- * @author M.Frank (code reshuffeling into new DDG4 scheme)
- * @version 1.0
- */
- class LcioGeneratorAction : public Simulation::Geant4GeneratorAction {
- public:
- typedef EVENT::MCParticle MCParticle;
- typedef IMPL::MCParticleImpl Particle;
- typedef IMPL::LCCollectionVec Particles;
- typedef std::set<MCParticle*> LCPrimarySet;
- typedef std::set<G4PrimaryParticle*> G4PrimarySet;
- typedef std::map<MCParticle*,G4PrimaryParticle*> G4ParticleMap;
- protected:
- /// Property: input file
- std::string m_input;
- /// Property: spread of primary vertex in Z
- double m_primaryVertexSpreadZ;
- /// Property: lorentz transformation angle
- double m_lorentzTransformationAngle;
- /// Property: SYNCEVT
- int m_SYNCEVT;
- /// map of partciles already import to Geant4 event
- G4ParticleMap m_g4ParticleMap;
- /// Set of LCIO particles already converted
- LCPrimarySet m_convertedParticles;
- /// Event reader object
- EventReader* m_reader;
-
- /// Read an event and return a LCCollectionVec of MCParticles.
- Particles* readEvent(int which);
- public:
- /// Standard constructor
- LcioGeneratorAction(Simulation::Geant4Context* context, const std::string& name);
- /// Default destructor
- virtual ~LcioGeneratorAction();
- /// Callback to generate primary particles
- virtual void operator()(G4Event*);
- G4PrimarySet getRelevantParticles(MCParticle* p);
- };
- } /* End namespace lcio */
-} /* End namespace DD4hep */
-
-#include "DD4hep/Plugins.h"
-namespace {
- /// Factory to create Geant4 physics constructions
- template <typename P> class Factory<P, DD4hep::lcio::EventReader*(std::string, int)> { public:
- static void Func(void *ret, void*, const std::vector<void*>& a, void*)
- { *(DD4hep::lcio::EventReader**)ret = (DD4hep::lcio::EventReader*)new P(*(std::string*)a[0],*(int*)a[1]);}
- };
-}
-/// Plugin defintion to create event reader objects
-#define DECLARE_LCIO_EVENT_READER(name) \
- PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::lcio::EventReader*(std::string,int))
-
-/// Plugin defintion to create event reader objects
-#define DECLARE_LCIO_EVENT_READER_NS(ns,name) typedef ns::name __##name##__; \
- PLUGINSVC_FACTORY_WITH_ID(__##name##__,std::string(#name),DD4hep::lcio::EventReader*(std::string,int))
-
-#endif /* DD4HEP_DDG4_LCIO_EVENTREADER_H */
diff --git a/DDG4/lcio/Geant4Output2LCIO.cpp b/DDG4/lcio/Geant4Output2LCIO.cpp
index 46019dfa9..ebad8b6a9 100644
--- a/DDG4/lcio/Geant4Output2LCIO.cpp
+++ b/DDG4/lcio/Geant4Output2LCIO.cpp
@@ -39,7 +39,6 @@ namespace DD4hep {
*/
class Geant4Output2LCIO : public Geant4OutputAction {
protected:
- Geometry::VolumeManager m_volMgr;
lcio::LCWriter* m_file;
int m_runNo;
@@ -84,6 +83,8 @@ namespace DD4hep {
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4DataConversion.h"
#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Particle.h"
+#include "DDG4/Geant4Data.h"
//#include "DDG4/Geant4Output2LCIO.h"
#include "G4Event.hh"
@@ -99,14 +100,13 @@ using namespace DD4hep::Simulation;
using namespace DD4hep;
using namespace std;
-
#include "DDG4/Factories.h"
DECLARE_GEANT4ACTION(Geant4Output2LCIO)
/// Standard constructor
Geant4Output2LCIO::Geant4Output2LCIO(Geant4Context* ctxt, const string& nam)
-: Geant4OutputAction(ctxt,nam), m_volMgr(), m_file(0), m_runNo(0)
+: Geant4OutputAction(ctxt,nam), m_file(0), m_runNo(0)
{
InstanceCount::increment(this);
}
@@ -125,9 +125,6 @@ void Geant4Output2LCIO::beginRun(const G4Run* ) {
if ( 0 == m_file && !m_output.empty() ) {
m_file = lcio::LCFactory::getInstance()->createLCWriter();
m_file->open(m_output,lcio::LCIO::WRITE_NEW);
- // Get the volume manager here:
- // in the constructor the geometry may not yet be built!
- m_volMgr = context()->lcdd().volumeManager();
}
}
@@ -159,7 +156,6 @@ void Geant4Output2LCIO::begin(const G4Event* /* event */) {
context()->event().addExtension<lcio::LCEventImpl>( e );
}
-
/// Callback to store the Geant4 event
void Geant4Output2LCIO::saveEvent(OutputContext<G4Event>& ctxt) {
lcio::LCEventImpl* e = context()->event().extension<lcio::LCEventImpl>();
@@ -167,19 +163,31 @@ void Geant4Output2LCIO::saveEvent(OutputContext<G4Event>& ctxt) {
e->setEventNumber(ctxt.context->GetEventID());
e->setDetectorName(context()->lcdd().header().name());
e->setRunNumber(m_runNo);
-
+ lcio::LCEventImpl* evt = context()->event().extension<lcio::LCEventImpl>();
+ Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>(false);
+ if ( part_map ) {
+ print("+++ Saving %d LCIO particles....",int(part_map->particleMap.size()));
+ if ( part_map->particleMap.size() > 0 ) {
+ typedef pair<const Geant4Context*,const Geant4ParticleMap*> _Args;
+ typedef Geant4Conversion<lcio::LCCollectionVec,_Args> _C;
+ const _C& cnv = _C::converter(typeid(Geant4ParticleMap));
+ lcio::LCCollectionVec* col = cnv(_Args(context(),part_map));
+ evt->addCollection(col,"McParticles");
+ }
+ }
}
/// Callback to store each Geant4 hit collection
void Geant4Output2LCIO::saveCollection(OutputContext<G4Event>& /* ctxt */, G4VHitsCollection* collection) {
size_t nhits = collection->GetSize();
std::string hc_nam = collection->GetName();
+ print("+++ Saving LCIO collection %s with %d entries....",hc_nam.c_str(),int(nhits));
if ( nhits > 0 ) {
- typedef pair<Geometry::VolumeManager,G4VHitsCollection*> _Args;
+ typedef pair<const Geant4Context*,G4VHitsCollection*> _Args;
typedef Geant4Conversion<lcio::LCCollectionVec,_Args> _C;
const _C& cnv = _C::converter(typeid(Geant4HitCollection));
lcio::LCEventImpl* evt = context()->event().extension<lcio::LCEventImpl>();
- lcio::LCCollectionVec* col = cnv(_Args(m_volMgr,collection));
+ lcio::LCCollectionVec* col = cnv(_Args(context(),collection));
evt->addCollection(col,hc_nam);
}
}
diff --git a/DDG4/lcio/HepEventReader.cpp b/DDG4/lcio/HepEventReader.cpp
index 3548c10df..60596e130 100644
--- a/DDG4/lcio/HepEventReader.cpp
+++ b/DDG4/lcio/HepEventReader.cpp
@@ -6,7 +6,7 @@
//====================================================================
// Framework include files
-#include "EventReader.h"
+#include "LCIOEventReader.h"
// C/C++ include files
#include <fstream>
@@ -18,7 +18,7 @@ namespace DD4hep {
/*
* lcio namespace declaration
*/
- namespace lcio {
+ namespace Simulation {
/** @class HepEventReader HepEventReader.h DDG4/HepEventReader.h
*
@@ -29,7 +29,7 @@ namespace DD4hep {
* @author M.Frank (code reshuffeling into new DDG4 scheme)
* @version 1.0
*/
- struct HepEventReader : public EventReader {
+ struct HepEventReader : public LCIOEventReader {
protected:
std::ifstream m_input;
int m_format;
@@ -39,7 +39,7 @@ namespace DD4hep {
/// Default destructor
virtual ~HepEventReader();
/// Read an event and return a LCCollectionVec of MCParticles.
- virtual Particles *readEvent();
+ virtual EVENT::LCCollection *readParticles();
};
} /* End namespace lcio */
} /* End namespace DD4hep */
@@ -51,15 +51,15 @@ namespace DD4hep {
#include "lcio.h"
using namespace IMPL;
using namespace EVENT;
+using namespace DD4hep::Simulation;
#define HEPEvt 1
// Factory entry
-typedef DD4hep::lcio::HepEventReader HepEventReader;
DECLARE_LCIO_EVENT_READER(HepEventReader)
/// Initializing constructor
HepEventReader::HepEventReader(const std::string& nam, int fmt)
-: EventReader(nam), m_format(fmt)
+: LCIOEventReader(nam), m_format(fmt)
{
}
@@ -74,7 +74,7 @@ MCParticleImpl* cast_particle(LCObject* obj) {
}
/// Read an event and return a LCCollectionVec of MCParticles.
-HepEventReader::Particles *HepEventReader::readEvent() {
+EVENT::LCCollection *HepEventReader::readParticles() {
static const double c_light = 299.792;// mm/ns
//
// Read the event, check for errors
diff --git a/DDG4/lcio/LCIOConversions.cpp b/DDG4/lcio/LCIOConversions.cpp
index 0cd2b1b1e..0e3a6515e 100644
--- a/DDG4/lcio/LCIOConversions.cpp
+++ b/DDG4/lcio/LCIOConversions.cpp
@@ -13,11 +13,14 @@
#include "DD4hep/Printout.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4DataConversion.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4Data.h"
// LCIO includes
#include "IMPL/LCCollectionVec.h"
//
+#include "IMPL/LCEventImpl.h"
#include "IMPL/ClusterImpl.h"
#include "IMPL/SimTrackerHitImpl.h"
#include "IMPL/SimCalorimeterHitImpl.h"
@@ -38,6 +41,7 @@ namespace DD4hep {
namespace Simulation {
typedef Geometry::VolumeManager VolMgr;
+ typedef Geometry::IDDescriptor IDDescriptor;
/// Data conversion interface calling lower level explicit convetrers
/**
@@ -46,7 +50,7 @@ namespace DD4hep {
*/
template <> lcio::LCCollectionVec*
Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,G4VHitsCollection*>,
+ pair<const Geant4Context*,G4VHitsCollection*>,
Geant4HitCollection>::operator()(const arg_t& args) const {
G4VHitsCollection* c = args.second;
Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(c);
@@ -60,10 +64,10 @@ namespace DD4hep {
"Cannot save the collection entries of:"+c->GetName());
}
- /// Data conversion interface creating lcio::SimTrackerHitImpl from SimpleTracker::Hit structures
+ /// Data conversion interface creating lcio::SimTrackerHitImpl from Geant4Tracker::Hit structures
/**
* This converter is to be used, when the sensitive detectors create fill collections
- * of type Geant4HitCollection with objects of type **SimpleTracker::Hit**.
+ * of type Geant4HitCollection with objects of type **Geant4Tracker::Hit**.
* The original objects are untouched and are automatically when the hosting
* Geant4HitCollection object is released.
*
@@ -72,38 +76,43 @@ namespace DD4hep {
*/
template <> lcio::LCCollectionVec*
Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,Geant4HitCollection*>,
- SimpleTracker::Hit>::operator()(const arg_t& args) const {
-
- Geant4HitCollection* coll = args.second;
- size_t nhits = coll->GetSize();
+ pair<const Geant4Context*,Geant4HitCollection*>,
+ Geant4Tracker::Hit>::operator()(const arg_t& args) const {
+
+ Geant4HitCollection* coll = args.second;
+ size_t nhits = coll->GetSize();
+ string dsc = encoding(coll->sensitiveDetector());
+ Geant4ParticleMap* pm = args.first->event().extension<Geant4ParticleMap>();
+ lcio::LCEventImpl* lc_evt = args.first->event().extension<lcio::LCEventImpl>();
+ EVENT::LCCollection* lc_part = lc_evt->getCollection("McParticles");
lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::SIMTRACKERHIT);
+ UTIL::CellIDEncoder<SimTrackerHit> decoder(dsc,lc_coll);
lc_coll->reserve(nhits);
- if ( nhits > 0 ) {
- SimpleHit* hit = coll->hit(0);
- string dsc = encoding(args.first, hit->cellID);
- UTIL::CellIDEncoder<SimTrackerHit> decoder(dsc,lc_coll);
- for(size_t i=0; i<nhits; ++i) {
- const SimpleTracker::Hit* g4_hit = coll->hit(i);
- double pos[3] = {g4_hit->position.x(), g4_hit->position.y(), g4_hit->position.z()};
- lcio::SimTrackerHitImpl* lc_hit = new lcio::SimTrackerHitImpl;
- lc_hit->setCellID0( (g4_hit->cellID >> 0 ) & 0xFFFFFFFF);
- lc_hit->setCellID1( (g4_hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF);
- lc_hit->setPosition(pos);
- lc_hit->setEDep(g4_hit->energyDeposit);
- lc_hit->setTime(g4_hit->truth.time);
- lc_hit->setMomentum( g4_hit->momentum.x(), g4_hit->momentum.y() , g4_hit->momentum.z() );
- lc_hit->setPathLength( g4_hit->length);
- lc_coll->addElement(lc_hit);
- }
+ for(size_t i=0; i<nhits; ++i) {
+ const Geant4Tracker::Hit* hit = coll->hit(i);
+ const Geant4Tracker::Hit::Contribution& t = hit->truth;
+ int trackID = pm->particleID(t.trackID);
+ EVENT::MCParticle* lc_mcp = (EVENT::MCParticle*)lc_part->getElementAt(trackID);
+ double pos[3] = {hit->position.x(), hit->position.y(), hit->position.z()};
+ float mom[3] = {hit->momentum.x(), hit->momentum.y(), hit->momentum.z()};
+ lcio::SimTrackerHitImpl* lc_hit = new lcio::SimTrackerHitImpl;
+ lc_hit->setCellID0( (hit->cellID >> 0 ) & 0xFFFFFFFF);
+ lc_hit->setCellID1( (hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF);
+ lc_hit->setEDep(hit->energyDeposit);
+ lc_hit->setPathLength(hit->length);
+ lc_hit->setTime(hit->truth.time);
+ lc_hit->setMCParticle(lc_mcp);
+ lc_hit->setPosition(pos);
+ lc_hit->setMomentum(mom);
+ lc_coll->addElement(lc_hit);
}
return lc_coll;
}
- /// Data conversion interface creating lcio::SimCalorimeterHitImpl from SimpleCalorimeter::Hit structures
+ /// Data conversion interface creating lcio::SimCalorimeterHitImpl from Geant4Calorimeter::Hit structures
/**
* This converter is to be used, when the sensitive detectors create fill collections
- * of type Geant4HitCollection with objects of type **SimpleCalorimeter::Hit**.
+ * of type Geant4HitCollection with objects of type **Geant4Calorimeter::Hit**.
* The original objects are untouched and are automatically when the hosting
* Geant4HitCollection object is released.
*
@@ -112,32 +121,43 @@ namespace DD4hep {
*/
template <> lcio::LCCollectionVec*
Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,Geant4HitCollection*>,
- SimpleCalorimeter::Hit>::operator()(const arg_t& args) const {
+ pair<const Geant4Context*,Geant4HitCollection*>,
+ Geant4Calorimeter::Hit>::operator()(const arg_t& args) const {
+ typedef Geant4HitData::Contributions Contributions;
Geant4HitCollection* coll = args.second;
- size_t nhits = coll->GetSize();
+ size_t nhits = coll->GetSize();
+ string dsc = encoding(coll->sensitiveDetector());
+ Geant4ParticleMap* pm = args.first->event().extension<Geant4ParticleMap>();
+ lcio::LCEventImpl* lc_evt = args.first->event().extension<lcio::LCEventImpl>();
+ EVENT::LCCollection* lc_parts = lc_evt->getCollection("McParticles");
lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::SIMCALORIMETERHIT);
+ UTIL::CellIDEncoder<SimCalorimeterHit> decoder(dsc,lc_coll);
lc_coll->setFlag(UTIL::make_bitset32(LCIO::CHBIT_LONG,LCIO::CHBIT_STEP,LCIO::CHBIT_ID1));
lc_coll->reserve(nhits);
- if ( nhits > 0 ) {
- SimpleHit* hit = coll->hit(0);
- string dsc = encoding(args.first, hit->cellID);
- UTIL::CellIDEncoder<SimCalorimeterHit> decoder(dsc,lc_coll);
- for(size_t i=0; i<nhits; ++i) {
- const SimpleCalorimeter::Hit* g4_hit = coll->hit(i);
- float pos[3] = {g4_hit->position.x(), g4_hit->position.y(), g4_hit->position.z()};
- lcio::SimCalorimeterHitImpl* lc_hit = new lcio::SimCalorimeterHitImpl;
- lc_hit->setCellID0( ( g4_hit->cellID >> 0 ) & 0xFFFFFFFF );
- lc_hit->setCellID1( ( g4_hit->cellID >> sizeof( int ) ) & 0xFFFFFFFF );
- lc_hit->setPosition(pos);
- lc_hit->setEnergy( g4_hit->energyDeposit );
- lc_coll->addElement(lc_hit);
+ for(size_t i=0; i<nhits; ++i) {
+ const Geant4Calorimeter::Hit* hit = coll->hit(i);
+ float pos[3] = {hit->position.x(), hit->position.y(), hit->position.z()};
+ lcio::SimCalorimeterHitImpl* lc_hit = new lcio::SimCalorimeterHitImpl;
+ lc_hit->setCellID0(( hit->cellID >> 0 ) & 0xFFFFFFFF);
+ lc_hit->setCellID1(( hit->cellID >> sizeof(int)) & 0xFFFFFFFF); // ????
+ lc_hit->setPosition(pos);
+ lc_hit->setEnergy( hit->energyDeposit );
+ lc_coll->addElement(lc_hit);
+ /// Now add the individual track contributions to the LCIO hit structure
+ for(Contributions::const_iterator j=hit->truth.begin(); j!=hit->truth.end(); ++j) {
+ const Geant4HitData::Contribution& c = *j;
+ int trackID = pm->particleID(c.trackID);
+ EVENT::MCParticle* lc_mcp = (EVENT::MCParticle*)lc_parts->getElementAt(trackID);
+ lc_hit->addMCParticleContribution(lc_mcp,c.deposit,c.time);
}
}
return lc_coll;
}
- template <typename T> lcio::LCCollectionVec* moveEntries(Geant4HitCollection* coll, lcio::LCCollectionVec* lc_coll) {
+ template <typename T>
+ lcio::LCCollectionVec* moveEntries(Geant4HitCollection* coll,
+ lcio::LCCollectionVec* lc_coll)
+ {
size_t nhits = coll->GetSize();
lc_coll->reserve(nhits);
for(size_t i=0; i<nhits; ++i) {
@@ -165,17 +185,13 @@ namespace DD4hep {
*/
template <> lcio::LCCollectionVec*
Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,Geant4HitCollection*>,
- lcio::SimTrackerHitImpl>::operator()(const arg_t& args) const
+ pair<const Geant4Context*,Geant4HitCollection*>,
+ lcio::SimTrackerHitImpl>::operator()(const arg_t& args) const
{
- Geant4HitCollection* coll = args.second;
- lcio::SimTrackerHitImpl* hit = coll->hit(0);
- long long int id1 = hit->getCellID1(), id0=hit->getCellID0();
- long long int cellID = (((id1<<32)&0xFFFFFFFF00000000)|(id0&0xFFFFFFFF));
- string dsc = encoding(args.first, cellID);
- lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::SIMTRACKERHIT);
- UTIL::CellIDEncoder<SimTrackerHit> decoder(dsc,lc_coll);
- return moveEntries<lcio::SimTrackerHitImpl>(args.second,lc_coll);
+ string dsc = encoding(args.second->sensitiveDetector());
+ output_t* lc = new lcio::LCCollectionVec(lcio::LCIO::SIMTRACKERHIT);
+ UTIL::CellIDEncoder<SimTrackerHit> decoder(dsc,lc);
+ return moveEntries<lcio::SimTrackerHitImpl>(args.second,lc);
}
/// Data conversion interface moving lcio::SimCalorimeterHitImpl objects from a Geant4HitCollection to a LCCollectionVec
@@ -193,15 +209,11 @@ namespace DD4hep {
*/
template <> lcio::LCCollectionVec*
Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,Geant4HitCollection*>,
+ pair<const Geant4Context*,Geant4HitCollection*>,
lcio::SimCalorimeterHitImpl>::operator()(const arg_t& args) const
{
- Geant4HitCollection* coll = args.second;
- lcio::SimCalorimeterHitImpl* hit = coll->hit(0);
- output_t* lc = new lcio::LCCollectionVec(lcio::LCIO::SIMCALORIMETERHIT);
- long long int id1 = hit->getCellID1(), id0=hit->getCellID0();
- long long int cellID = (((id1<<32)&0xFFFFFFFF00000000)|(id0&0xFFFFFFFF));
- string dsc = encoding(args.first, cellID);
+ string dsc = encoding(args.second->sensitiveDetector());
+ output_t* lc = new lcio::LCCollectionVec(lcio::LCIO::SIMCALORIMETERHIT);
UTIL::CellIDEncoder<SimCalorimeterHit> decoder(dsc,lc);
lc->setFlag(UTIL::make_bitset32(LCIO::CHBIT_LONG,LCIO::CHBIT_STEP,LCIO::CHBIT_ID1));
return moveEntries<tag_t>(args.second,lc);
@@ -215,41 +227,26 @@ namespace DD4hep {
*/
template <> lcio::LCCollectionVec*
Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,Geant4HitCollection*>,
+ pair<const Geant4Context*,Geant4HitCollection*>,
lcio::ClusterImpl>::operator()(const arg_t& args) const
{
output_t* lc = new lcio::LCCollectionVec(lcio::LCIO::CLUSTER);
return moveEntries<tag_t>(args.second,lc);
}
- /// Data conversion interface moving lcio::MCParticleImpl objects from a Geant4HitCollection to a LCCollectionVec
- /** Same comments apply as for the data mover for lcio::SimCalorimeterHitImpl and lcio::SimTrackerHitImpl
- *
- * @author M.Frank
- * @version 1.0
- */
- template <> lcio::LCCollectionVec*
- Geant4DataConversion<lcio::LCCollectionVec,
- pair<VolMgr,Geant4HitCollection*>,
- lcio::MCParticleImpl>::operator()(const arg_t& args) const
- {
- output_t* lc = new lcio::LCCollectionVec(lcio::LCIO::MCPARTICLE);
- return moveEntries<tag_t>(args.second,lc);
- }
- typedef pair<VolMgr,G4VHitsCollection*> RAW_CONVERSION_ARGS;
- typedef pair<VolMgr,Geant4HitCollection*> CONVERSION_ARGS;
+ typedef pair<const Geant4Context*,G4VHitsCollection*> RAW_CONVERSION_ARGS;
+ typedef pair<const Geant4Context*,Geant4HitCollection*> CONVERSION_ARGS;
template class Geant4Conversion<lcio::LCCollectionVec,RAW_CONVERSION_ARGS>;
DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,RAW_CONVERSION_ARGS,Geant4HitCollection)
template class Geant4Conversion<lcio::LCCollectionVec,CONVERSION_ARGS>;
// Hit converters for simple Geant4Data objects
- DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,SimpleTracker::Hit)
- DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,SimpleCalorimeter::Hit)
+ DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,Geant4Tracker::Hit)
+ DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,Geant4Calorimeter::Hit)
// Hit converters for standard LCIO objects
DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,lcio::SimTrackerHitImpl)
DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,lcio::SimCalorimeterHitImpl)
DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,lcio::ClusterImpl)
- DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,lcio::MCParticleImpl)
} // End namespace Simulation
} // End namespace DD4hep
diff --git a/DDG4/lcio/LCIOEventReader.cpp b/DDG4/lcio/LCIOEventReader.cpp
new file mode 100644
index 000000000..a4ffa7f2e
--- /dev/null
+++ b/DDG4/lcio/LCIOEventReader.cpp
@@ -0,0 +1,20 @@
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// @author P.Kostka (main author)
+// @author M.Frank (code reshuffeling into new DDG4 scheme)
+//
+//====================================================================
+
+// Framework include files
+#include "LCIOEventReader.h"
+
+/// Initializing constructor
+DD4hep::Simulation::LCIOEventReader::LCIOEventReader(const std::string& nam) : m_name(nam) {
+}
+
+/// Default destructor
+DD4hep::Simulation::LCIOEventReader::~LCIOEventReader() {
+}
diff --git a/DDG4/lcio/LCIOEventReader.h b/DDG4/lcio/LCIOEventReader.h
new file mode 100644
index 000000000..b7b2a5eaf
--- /dev/null
+++ b/DDG4/lcio/LCIOEventReader.h
@@ -0,0 +1,82 @@
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_LCIOEVENTREADER_H
+#define DD4HEP_DDG4_LCIOEVENTREADER_H
+
+// C/C++ include files
+#include <string>
+
+// Forward declarations
+class G4Event;
+namespace IO { class LCReader; }
+namespace UTIL { class LCStdHepRdr; }
+namespace EVENT { class MCParticle; }
+namespace EVENT { class LCCollection; }
+namespace IMPL { class MCParticleImpl; }
+namespace IMPL { class LCCollectionVec; }
+
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Particle;
+
+ /** @class EventReader EventReader.h DDG4/EventReader.h
+ *
+ * Base class to read lcio files.
+ *
+ * @author P.Kostka (main author)
+ * @author M.Frank (code reshuffeling into new DDG4 scheme)
+ * @version 1.0
+ */
+ class LCIOEventReader {
+
+ public:
+ typedef Geant4Particle Particle;
+ protected:
+ /// File name
+ std::string m_name;
+
+ public:
+ /// Initializing constructor
+ LCIOEventReader(const std::string& nam);
+ /// Default destructor
+ virtual ~LCIOEventReader();
+ /// File name
+ const std::string& name() const { return m_name; }
+ /// Read an event and return a LCCollectionVec of MCParticles.
+ virtual EVENT::LCCollection *readParticles() = 0;
+ };
+
+ } /* End namespace Simulation */
+} /* End namespace DD4hep */
+
+#include "DD4hep/Plugins.h"
+namespace {
+ /// Factory to create Geant4 physics constructions
+ template <typename P> class Factory<P, DD4hep::Simulation::LCIOEventReader*(std::string, int)> { public:
+ static void Func(void *ret, void*, const std::vector<void*>& a, void*)
+ { *(DD4hep::Simulation::LCIOEventReader**)ret = (DD4hep::Simulation::LCIOEventReader*)new P(*(std::string*)a[0],*(int*)a[1]);}
+ };
+}
+/// Plugin defintion to create event reader objects
+#define DECLARE_LCIO_EVENT_READER(name) \
+ PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::Simulation::LCIOEventReader*(std::string,int))
+
+/// Plugin defintion to create event reader objects
+#define DECLARE_LCIO_EVENT_READER_NS(ns,name) typedef ns::name __##name##__; \
+ PLUGINSVC_FACTORY_WITH_ID(__##name##__,std::string(#name),DD4hep::Simulation::LCIOEventReader*(std::string,int))
+
+#endif /* DD4HEP_DDG4_LCIOEVENTREADER_H */
diff --git a/DDG4/lcio/LcioEventReader.cpp b/DDG4/lcio/LCIOFileReader.cpp
similarity index 66%
rename from DDG4/lcio/LcioEventReader.cpp
rename to DDG4/lcio/LCIOFileReader.cpp
index 25e4bbf89..dfa955cfb 100644
--- a/DDG4/lcio/LcioEventReader.cpp
+++ b/DDG4/lcio/LCIOFileReader.cpp
@@ -6,7 +6,7 @@
//====================================================================
// Framework include files
-#include "EventReader.h"
+#include "LCIOEventReader.h"
/*
* DD4hep namespace declaration
@@ -15,7 +15,7 @@ namespace DD4hep {
/*
* lcio namespace declaration
*/
- namespace lcio {
+ namespace Simulation {
/** @class LcioEventReader LcioEventReader.h DDG4/LcioEventReader.h
*
@@ -25,52 +25,52 @@ namespace DD4hep {
* @author M.Frank (code reshuffeling into new DDG4 scheme)
* @version 1.0
*/
- struct LcioEventReader : public EventReader {
+ struct LCIOFileReader : public LCIOEventReader {
protected:
/// Reference to reader object
IO::LCReader* m_reader;
public:
/// Initializing constructor
- LcioEventReader(const std::string& nam, int);
+ LCIOFileReader(const std::string& nam, int);
/// Default destructor
- virtual ~LcioEventReader();
+ virtual ~LCIOFileReader();
/// Read an event and return a LCCollectionVec of MCParticles.
- virtual Particles *readEvent();
+ virtual EVENT::LCCollection *readParticles();
};
}
}
#include "DD4hep/Printout.h"
+#include "DD4hep/Primitives.h"
#include "lcio.h"
#include "EVENT/LCIO.h"
#include "IO/LCReader.h"
using namespace EVENT;
+using namespace DD4hep::Simulation;
// Factory entry
-typedef DD4hep::lcio::LcioEventReader LcioEventReader;
-DECLARE_LCIO_EVENT_READER(LcioEventReader)
+DECLARE_LCIO_EVENT_READER(LCIOFileReader)
/// Initializing constructor
-LcioEventReader::LcioEventReader(const std::string& nam, int /* arg */)
-: EventReader(nam)
+LCIOFileReader::LCIOFileReader(const std::string& nam, int /* arg */)
+: LCIOEventReader(nam)
{
m_reader = ::lcio::LCFactory::getInstance()->createLCReader();
- printout(INFO,"LcioEventReader","Created file reader. Try to open input %s",nam.c_str());
+ printout(INFO,"LCIOFileReader","Created file reader. Try to open input %s",nam.c_str());
m_reader->open(nam);
}
/// Default destructor
-LcioEventReader::~LcioEventReader() {
- deletePtr(m_reader);
+LCIOFileReader::~LCIOFileReader() {
+ DD4hep::deletePtr(m_reader);
}
/// Read an event and return a LCCollectionVec of MCParticles.
-DD4hep::lcio::EventReader::Particles *LcioEventReader::readEvent() {
+EVENT::LCCollection *LCIOFileReader::readParticles() {
::lcio::LCEvent* evt = m_reader->readNextEvent();
if ( evt ) {
- Particles* mcVec = (Particles*)evt->getCollection(LCIO::MCPARTICLE);
- return mcVec;
+ return evt->getCollection(LCIO::MCPARTICLE);
}
return 0;
}
diff --git a/DDG4/lcio/LCIOInputAction.cpp b/DDG4/lcio/LCIOInputAction.cpp
new file mode 100644
index 000000000..697a6f705
--- /dev/null
+++ b/DDG4/lcio/LCIOInputAction.cpp
@@ -0,0 +1,344 @@
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// @author P.Kostka (main author)
+// @author M.Frank (code reshuffeling into new DDG4 scheme)
+//
+//====================================================================
+
+// Framework include files
+#include "LCIOInputAction.h"
+#include "LCIOEventReader.h"
+#include "DD4hep/Printout.h"
+#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"
+
+using namespace std;
+using namespace DD4hep::Simulation;
+
+typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+
+/// Standard constructor
+LCIOInputAction::LCIOInputAction(Geant4Context* context, const string& name)
+ : Geant4GeneratorAction(context,name), m_reader(0)
+{
+ declareProperty("Input", m_input);
+ declareProperty("Sync", m_SYNCEVT=0);
+ declareProperty("Mask", m_mask = 0);
+ declareProperty("MomentumScale", m_momScale = 1.0);
+ m_needsControl = true;
+}
+
+/// Default destructor
+LCIOInputAction::~LCIOInputAction() {
+}
+
+/// helper to report Geant4 exceptions
+string LCIOInputAction::issue(int i) const {
+ stringstream str;
+ str << "LCIOInputAction[" << name() << "]: Event " << i << " ";
+ return str.str();
+}
+
+/// Read an event and return a LCCollection of MCParticles.
+EVENT::LCCollection* LCIOInputAction::readParticles(int which) {
+ int evid = which + m_SYNCEVT;
+ if ( 0 == m_reader ) {
+ if ( m_input.empty() ) {
+ throw runtime_error("InputAction: No inoput file declared!");
+ }
+ string err;
+ TypeName tn = TypeName::split(m_input,"|");
+ try {
+ m_reader = PluginService::Create<LCIOEventReader*>(tn.first,tn.second,int(0));
+ if ( 0 == m_reader ) {
+ abortRun(issue(evid)+"Error creating reader plugin.",
+ "Failed to create file reader of type %s. Cannot open dataset %s",
+ tn.first.c_str(),tn.second.c_str());
+ return 0;
+ }
+ }
+ catch(const std::exception& e) {
+ err = e.what();
+ }
+ if ( !err.empty() ) {
+ abortRun(issue(evid)+err,"Error when reading file %s",m_input.c_str());
+ return 0;
+ }
+ }
+ EVENT::LCCollection* parts = m_reader->readParticles();
+ if ( 0 == parts ) {
+ abortRun(issue(evid)+"Error when reading file - may be end of file.",
+ "Error when reading file %s",m_input.c_str());
+ }
+ 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());
+ map<EVENT::MCParticle*,int> mcparts;
+ vector<EVENT::MCParticle*> mcpcoll;
+ vector<Particle*> col;
+
+ prim->add(m_mask, inter);
+ // check if there is at least one particle
+ if ( 0 == primaries ) return;
+
+ int NHEP = primaries->getNumberOfElements();
+ // 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);
+ // build collection of MCParticles
+ G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
+ for(int i=0; i<NHEP; ++i ) {
+ EVENT::MCParticle* mcp = mcpcoll[i];
+ Geant4ParticleHandle p(col[i]);
+ const double mom_scale = m_momScale;
+ const double *mom = mcp->getMomentum();
+ const double *vsx = mcp->getVertex();
+ const double *vex = mcp->getEndpoint();
+ const float *spin = mcp->getSpin();
+ 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->psx = mom_scale*mom[0]*GeV;
+ p->psy = mom_scale*mom[1]*GeV;
+ p->psz = mom_scale*mom[2]*GeV;
+ p->time = mcp->getTime()*ns;
+ p->properTime = mcp->getTime()*ns;
+ p->vsx = vsx[0]*mm;
+ p->vsy = vsx[1]*mm;
+ p->vsz = vsx[2]*mm;
+ p->vex = vex[0]*mm;
+ p->vey = vex[1]*mm;
+ p->vez = vex[2]*mm;
+ p->definition = def;
+ p->process = 0;
+ p->spin[0] = spin[0];
+ p->spin[1] = spin[1];
+ p->spin[2] = spin[2];
+ p->colorFlow[0] = color[0];
+ p->colorFlow[0] = color[1];
+ p->mass = mcp->getMass()*GeV;
+ const EVENT::MCParticleVec &par = mcp->getParents(), &dau=mcp->getDaughters();
+ for(int num=dau.size(),k=0; k<num; ++k)
+ p->daughters.insert(GET_ENTRY(mcparts,dau[k]));
+ for(int num=par.size(),k=0; k<num; ++k)
+ 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);
+ if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ if ( mcp->isCreatedInSimulation() ) status.set(G4PARTICLE_SIM_CREATED);
+ if ( mcp->isBackscatter() ) status.set(G4PARTICLE_SIM_BACKSCATTER);
+ if ( mcp->vertexIsNotEndpointOfParent() ) status.set(G4PARTICLE_SIM_PARENT_RADIATED);
+ if ( mcp->isDecayedInTracker() ) status.set(G4PARTICLE_SIM_DECAY_TRACKER);
+ if ( mcp->isDecayedInCalorimeter() ) status.set(G4PARTICLE_SIM_DECAY_CALO);
+ 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);
+ 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"
+DECLARE_GEANT4ACTION(LCIOInputAction)
diff --git a/DDG4/lcio/LCIOInputAction.h b/DDG4/lcio/LCIOInputAction.h
new file mode 100644
index 000000000..9014f266c
--- /dev/null
+++ b/DDG4/lcio/LCIOInputAction.h
@@ -0,0 +1,87 @@
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_LCIOINPUTACTION_H
+#define DD4HEP_DDG4_LCIOINPUTACTION_H
+
+// Framework include files
+#include "DD4hep/Primitives.h"
+#include "DDG4/Geant4Particle.h"
+#include "DDG4/Geant4GeneratorAction.h"
+
+// C/C++ include files
+#include <set>
+#include <map>
+#include <vector>
+
+// Forward declarations
+class G4Event;
+namespace IO { class LCReader; }
+namespace UTIL { class LCStdHepRdr; }
+namespace EVENT { class MCParticle; }
+namespace EVENT { class LCCollection; }
+namespace IMPL { class MCParticleImpl; }
+namespace IMPL { class LCCollectionVec; }
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Particle;
+ class LCIOEventReader;
+
+ /** @class InputAction Geant4GeneratorAction.h DDG4/Geant4GeneratorAction.h
+ *
+ * Concrete implementation of the Geant4 generator action base class
+ * populating Geant4 primaries from LCIO and HepStd files.
+ *
+ * @author P.Kostka (main author)
+ * @author M.Frank (code reshuffeling into new DDG4 scheme)
+ * @version 1.0
+ */
+ class LCIOInputAction : public Simulation::Geant4GeneratorAction {
+
+ public:
+ typedef Geant4Particle Particle;
+
+ protected:
+ /// Property: input file
+ std::string m_input;
+ /// Property: SYNCEVT
+ int m_SYNCEVT;
+ /// Property; interaction mask
+ int m_mask;
+ /// Property: Momentum downscaler for debugging
+ double m_momScale;
+ /// Event reader object
+ LCIOEventReader* m_reader;
+
+ /// Read an event and return a LCCollectionVec of MCParticles.
+ EVENT::LCCollection* readParticles(int which);
+ /// helper to report Geant4 exceptions
+ std::string issue(int i) const;
+
+ public:
+ /// Standard constructor
+ LCIOInputAction(Simulation::Geant4Context* context, const std::string& name);
+ /// Default destructor
+ 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 */
+
+#endif /* DD4HEP_DDG4_LCIOINPUTACTION_H */
diff --git a/DDG4/lcio/LCIOMCParticles.cpp b/DDG4/lcio/LCIOMCParticles.cpp
new file mode 100644
index 000000000..709283155
--- /dev/null
+++ b/DDG4/lcio/LCIOMCParticles.cpp
@@ -0,0 +1,150 @@
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+
+// Framework include files
+#define DDG4_MAKE_INSTANTIATIONS
+#include "DD4hep/LCDD.h"
+#include "DD4hep/Printout.h"
+#include "DDG4/Geant4Particle.h"
+#include "DDG4/Geant4HitCollection.h"
+#include "DDG4/Geant4DataConversion.h"
+#include "DDG4/Geant4MonteCarloTruth.h"
+
+// LCIO includes
+#include "IMPL/LCCollectionVec.h"
+//
+#include "IMPL/ClusterImpl.h"
+#include "IMPL/SimTrackerHitImpl.h"
+#include "IMPL/SimCalorimeterHitImpl.h"
+#include "IMPL/MCParticleImpl.h"
+#include "UTIL/ILDConf.h"
+
+// Geant4 include files
+#include "G4ParticleDefinition.hh"
+#include "G4VProcess.hh"
+
+using namespace std;
+
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ typedef ReferenceBitMask<const int> PropertyMask;
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ typedef Geant4Particle Particle;
+
+ /// Data conversion interface for MC particles to LCIO format
+ /**
+ * @author M.Frank
+ * @version 1.0
+ */
+ template <> lcio::LCCollectionVec*
+ Geant4DataConversion<lcio::LCCollectionVec,
+ pair<const Geant4Context*,const Geant4ParticleMap*>,
+ Geant4ParticleMap>::operator()(const arg_t& args) const {
+ typedef MCParticleImpl MYParticleImpl;
+ typedef Geant4Conversion<output_t,pair<arg_t::first_type,arg_t::second_type> > _C;
+ typedef Geant4ParticleMap::ParticleMap ParticleMap;
+ const ParticleMap& pm = args.second->particleMap;
+ size_t nparts = pm.size();
+ lcio::LCCollectionVec* lc_coll = new lcio::LCCollectionVec(lcio::LCIO::MCPARTICLE);
+ lc_coll->reserve(nparts);
+ if ( nparts > 0 ) {
+ size_t cnt = 0;
+ map<int,int> p_ids;
+ vector<const Particle*> p_part(pm.size(),0);
+ vector<MYParticleImpl*> p_lcio(pm.size(),0);
+ // First create the particles
+ for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end();++i, ++cnt) {
+ int id = (*i).first;
+ const Particle* p = (*i).second;
+ PropertyMask mask(p->status);
+ const G4ParticleDefinition* def = p->definition;
+ MYParticleImpl* q = (MYParticleImpl*)new lcio::MCParticleImpl();
+ q->setPDG(p->pdgID);
+ q->setMomentum(&p->psx);
+ q->setVertex(&p->vsx);
+ q->setEndpoint(&p->vex);
+ q->setTime(p->time);
+ q->setMass(p->mass);
+ q->setCharge(def ? def->GetPDGCharge()/3.0 : 0); // Charge(e+) = 1 !
+
+ // Set generator status
+ if ( mask.isSet(G4PARTICLE_GEN_EMPTY) ) q->setGeneratorStatus(0);
+ else if ( mask.isSet(G4PARTICLE_GEN_STABLE) ) q->setGeneratorStatus(1);
+ else if ( mask.isSet(G4PARTICLE_GEN_DECAYED) ) q->setGeneratorStatus(2);
+ else if ( mask.isSet(G4PARTICLE_GEN_DOCUMENTATION) ) q->setGeneratorStatus(3);
+
+ // Set simulation status
+ q->setCreatedInSimulation( mask.isSet(G4PARTICLE_SIM_CREATED) );
+ q->setBackscatter( mask.isSet(G4PARTICLE_SIM_BACKSCATTER) );
+ q->setVertexIsNotEndpointOfParent( mask.isSet(G4PARTICLE_SIM_PARENT_RADIATED) );
+ q->setDecayedInTracker( mask.isSet(G4PARTICLE_SIM_DECAY_TRACKER) );
+ q->setDecayedInCalorimeter( mask.isSet(G4PARTICLE_SIM_DECAY_CALO) );
+ q->setHasLeftDetector( mask.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) );
+ q->setStopped( mask.isSet(G4PARTICLE_SIM_STOPPED) );
+ q->setOverlay( false );
+
+ q->setSpin(p->spin);
+ q->setColorFlow(p->colorFlow);
+
+ lc_coll->addElement(q);
+ p_ids[id] = cnt;
+ p_part[cnt] = p;
+ p_lcio[cnt] = q;
+ }
+
+ // Now establish parent-daughter relationships
+ for(size_t i=0, n=p_ids.size(); i<n; ++i) {
+ map<int,int>::iterator k;
+ const Particle* p = p_part[i];
+ MYParticleImpl* q = p_lcio[i];
+ const Particle::Particles& dau = p->daughters;
+ for(Particle::Particles::const_iterator j=dau.begin(); j!=dau.end(); ++j) {
+ int idau = *j;
+ if ( (k=p_ids.find(idau)) == p_ids.end() ) { // Error!!!
+ printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find daughter with ID:%d",p->id,idau);
+ continue;
+ }
+ int iqdau = (*k).second;
+ MYParticleImpl* qdau = p_lcio[iqdau];
+ qdau->addParent(q);
+ }
+ const Particle::Particles& par = p->parents;
+ for(Particle::Particles::const_iterator j=par.begin(); j!=par.end(); ++j) {
+ int ipar = *j; // A parent ID iof -1 means NO parent, because a base of 0 is perfectly leagal!
+ if ( ipar>=0 && (k=p_ids.find(ipar)) == p_ids.end() ) { // Error!!!
+ printout(FATAL,"Geant4Conversion","+++ Particle %d: FAILED to find parent with ID:%d",p->id,ipar);
+ continue;
+ }
+ int iqpar = (*k).second;
+ MYParticleImpl* qpar = p_lcio[iqpar];
+ q->addParent(qpar);
+ }
+ }
+ }
+ return lc_coll;
+ }
+
+ typedef pair<const Geant4Context*,const Geant4ParticleMap*> CONVERSION_ARGS;
+ template class Geant4Conversion<lcio::LCCollectionVec,CONVERSION_ARGS>;
+ DECLARE_GEANT4_HITCONVERTER(lcio::LCCollectionVec,CONVERSION_ARGS,Geant4ParticleMap)
+ } // End namespace Simulation
+} // End namespace DD4hep
+
+
+
diff --git a/DDG4/lcio/LcioStdHepReader.cpp b/DDG4/lcio/LCIOStdHepReader.cpp
similarity index 62%
rename from DDG4/lcio/LcioStdHepReader.cpp
rename to DDG4/lcio/LCIOStdHepReader.cpp
index 0bf0b6b31..8da479885 100644
--- a/DDG4/lcio/LcioStdHepReader.cpp
+++ b/DDG4/lcio/LCIOStdHepReader.cpp
@@ -6,18 +6,18 @@
//====================================================================
// Framework include files
-#include "EventReader.h"
+#include "LCIOEventReader.h"
/*
* DD4hep namespace declaration
*/
namespace DD4hep {
/*
- * lcio namespace declaration
+ * Simulation namespace declaration
*/
- namespace lcio {
+ namespace Simulation {
- /** @class StdHepReader StdHepReader.h DDG4/StdHepReader.h
+ /** @class LCIOStdHepReader LCIOStdHepReader.h DDG4/LCIOStdHepReader.h
*
* Base class to read StdHep files with lcio
*
@@ -25,17 +25,17 @@ namespace DD4hep {
* @author M.Frank (code reshuffeling into new DDG4 scheme)
* @version 1.0
*/
- struct StdHepReader : public EventReader {
+ struct LCIOStdHepReader : public LCIOEventReader {
protected:
- /// Reference to reader object
+ /// Reference to Reader object
UTIL::LCStdHepRdr* m_reader;
public:
/// Initializing constructor
- StdHepReader(const std::string& nam, int);
+ LCIOStdHepReader(const std::string& nam, int);
/// Default destructor
- virtual ~StdHepReader();
+ virtual ~LCIOStdHepReader();
/// Read an event and return a LCCollectionVec of MCParticles.
- virtual Particles *readEvent();
+ virtual EVENT::LCCollection *readParticles();
};
} /* End namespace lcio */
} /* End namespace DD4hep */
@@ -44,25 +44,26 @@ namespace DD4hep {
#include "lcio.h"
#include "EVENT/LCIO.h"
#include "UTIL/LCStdHepRdr.h"
+#include "DD4hep/Primitives.h"
+
+using namespace DD4hep::Simulation;
// Factory entry
-typedef DD4hep::lcio::StdHepReader LcioStdHepReader;
-DECLARE_LCIO_EVENT_READER(LcioStdHepReader)
+DECLARE_LCIO_EVENT_READER(LCIOStdHepReader)
/// Initializing constructor
-DD4hep::lcio::StdHepReader::StdHepReader(const std::string& nam, int)
- : EventReader(nam)
+LCIOStdHepReader::LCIOStdHepReader(const std::string& nam, int)
+ : LCIOEventReader(nam)
{
m_reader = new ::lcio::LCStdHepRdr(m_name.c_str());
}
/// Default destructor
-DD4hep::lcio::StdHepReader::~StdHepReader() {
- deletePtr(m_reader);
+LCIOStdHepReader::~LCIOStdHepReader() {
+ DD4hep::deletePtr(m_reader);
}
/// Read an event and return a LCCollectionVec of MCParticles.
-DD4hep::lcio::EventReader::Particles *DD4hep::lcio::StdHepReader::readEvent() {
+EVENT::LCCollection *LCIOStdHepReader::readParticles() {
return m_reader->readEvent();
}
-
diff --git a/DDG4/plugins/Geant4EscapeCounter.cpp b/DDG4/plugins/Geant4EscapeCounter.cpp
index 4b9a46ee3..a520ebd59 100644
--- a/DDG4/plugins/Geant4EscapeCounter.cpp
+++ b/DDG4/plugins/Geant4EscapeCounter.cpp
@@ -111,8 +111,9 @@ bool Geant4EscapeCounter::process(G4Step* step, G4TouchableHistory* /* history *
coll->add(hit);
mark(h.track);
- print(name(),"+++ Track:%4d %8.2f MeV [%s] %s Geant4 path:%s",
- h.trkID(),h.trkEnergy()/MeV,th.name().c_str(),th.creatorName().c_str(),path.c_str());
+ print("+++ Track:%4d %8.2f MeV [%s] %s Geant4 path:%s",
+ h.trkID(),h.trkEnergy()/MeV,th.name().c_str(),
+ th.creatorName().c_str(),path.c_str());
// Kill track, so that it does no longer participate in the propagation
h.track->SetTrackStatus(fStopAndKill);
return true;
diff --git a/DDG4/plugins/Geant4Factories.cpp b/DDG4/plugins/Geant4Factories.cpp
index f517b8ea5..e4cfc7fa9 100644
--- a/DDG4/plugins/Geant4Factories.cpp
+++ b/DDG4/plugins/Geant4Factories.cpp
@@ -68,6 +68,28 @@ DECLARE_GEANT4ACTION(Geant4Output2ROOT)
DECLARE_GEANT4ACTION(Geant4ParticleGun)
//=============================
+#include "DDG4/Geant4GeneratorActionInit.h"
+DECLARE_GEANT4ACTION(Geant4GeneratorActionInit)
+
+//=============================
+#include "DDG4/Geant4IsotropeGenerator.h"
+DECLARE_GEANT4ACTION(Geant4IsotropeGenerator)
+
+//=============================
+#include "DDG4/Geant4InteractionVertexSmear.h"
+DECLARE_GEANT4ACTION(Geant4InteractionVertexSmear)
+
+//=============================
+#include "DDG4/Geant4InteractionVertexBoost.h"
+DECLARE_GEANT4ACTION(Geant4InteractionVertexBoost)
+
+//=============================
+#include "DDG4/Geant4InteractionMerger.h"
+DECLARE_GEANT4ACTION(Geant4InteractionMerger)
+
+//=============================
+#include "DDG4/Geant4PrimaryHandler.h"
+DECLARE_GEANT4ACTION(Geant4PrimaryHandler)
//=============================
#include "DDG4/Geant4TestActions.h"
diff --git a/DDG4/python/DDG4Dict.C b/DDG4/python/DDG4Dict.C
index 70ede136d..245d85783 100644
--- a/DDG4/python/DDG4Dict.C
+++ b/DDG4/python/DDG4Dict.C
@@ -10,6 +10,7 @@
//
//====================================================================
// FRamework include files
+#include "DDG4/Geant4Particle.h"
#include "DDG4/Geant4Data.h"
#include <vector>
@@ -26,8 +27,12 @@ using namespace DD4hep::Simulation;
#pragma link C++ class SimpleEvent+;
//#pragma link C++ class SimpleEvent::Seeds+;
#pragma link C++ class DataExtension+;
+#pragma link C++ class ParticleExtension+;
+#pragma link C++ class auto_ptr<DataExtension>+;
+#pragma link C++ class auto_ptr<ParticleExtension>+;
+
#pragma link C++ class SimpleHit+;
-#pragma link C++ class Particle+;
+#pragma link C++ class Geant4Particle+;
#pragma link C++ class std::vector<SimpleHit*>+;
#pragma link C++ class SimpleHit::Contribution+;
#pragma link C++ class SimpleHit::Contributions+;
@@ -37,7 +42,7 @@ using namespace DD4hep::Simulation;
#pragma link C++ class SimpleCalorimeter+;
#pragma link C++ class SimpleCalorimeter::Hit+;
#pragma link C++ class std::vector<SimpleCalorimeter::Hit*>+;
-#pragma link C++ class std::vector<Particle*>+;
+#pragma link C++ class std::vector<Geant4Particle*>+;
//#pragma link C++ class ;
#endif
@@ -218,6 +223,7 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation;
#pragma link C++ class Geant4ActionCreation;
#pragma link C++ class Geant4Action;
#pragma link C++ class Geant4Kernel;
+#pragma link C++ class Geant4Kernel::PhaseSelector;
#pragma link C++ class Geant4Context;
#pragma link C++ class KernelHandle;
@@ -249,6 +255,7 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation;
#pragma link C++ class Geant4Sensitive;
#pragma link C++ class Geant4SensDetActionSequence;
#pragma link C++ class Geant4ActionPhase;
+
#endif
int Geant4Dict() {
diff --git a/DDG4/src/Geant4Action.cpp b/DDG4/src/Geant4Action.cpp
index 5897c58fb..f98389b19 100644
--- a/DDG4/src/Geant4Action.cpp
+++ b/DDG4/src/Geant4Action.cpp
@@ -21,7 +21,7 @@ using namespace std;
using namespace DD4hep;
using namespace DD4hep::Simulation;
-TypeName TypeName::split(const string& type_name, const std::string& delim) {
+TypeName TypeName::split(const string& type_name, const string& delim) {
size_t idx = type_name.find(delim);
string typ = type_name, nam = type_name;
if (idx != string::npos) {
@@ -75,10 +75,8 @@ long Geant4Action::addRef() {
long Geant4Action::release() {
long count = --m_refCount;
if (m_refCount <= 0) {
- cout << "Geant4Action: Deleting object " << name()
- << " of type " << typeName(typeid(*this))
- << " Ptr:" << (void*)this
- << endl;
+ print("Geant4Action: Deleting object %s of type %s Pointer:%p",
+ m_name.c_str(),typeName(typeid(*this)).c_str(),(void*)this);
delete this;
}
return count;
@@ -98,12 +96,12 @@ Geant4Action& Geant4Action::setProperties(PropertyConfigurator& setup) {
}
/// Check property for existence
-bool Geant4Action::hasProperty(const std::string& name) const {
+bool Geant4Action::hasProperty(const string& name) const {
return m_properties.exists(name);
}
/// Access single property
-Property& Geant4Action::property(const std::string& name) {
+Property& Geant4Action::property(const string& name) {
return properties()[name];
}
@@ -144,91 +142,94 @@ void Geant4Action::enableUI() {
}
/// Support for messages with variable output level using output level
-void Geant4Action::print(const std::string& fmt, ...) const {
- if ( outputLevel() >= printLevel() ) {
+void Geant4Action::print(const char* fmt, ...) const {
+ int level = outputLevel();
+ if ( level >= printLevel() ) {
va_list args;
va_start(args, fmt);
- DD4hep::printout((PrintLevel)outputLevel(),m_name, fmt, args);
+ DD4hep::printout((PrintLevel)level, m_name.c_str(), fmt, args);
va_end(args);
}
}
-/// Support for messages with variable output level using output level
-void Geant4Action::print(const std::string& tag, const std::string& fmt, ...) const {
- if ( outputLevel() >= printLevel() ) {
+/// Support for messages with variable output level using output level-1
+void Geant4Action::printM1(const char* fmt, ...) const {
+ int level = max(outputLevel()-1,(int)VERBOSE);
+ if ( level >= printLevel() ) {
va_list args;
va_start(args, fmt);
- DD4hep::printout((PrintLevel)outputLevel(),tag, fmt, args);
+ DD4hep::printout((PrintLevel)level, m_name.c_str(), fmt, args);
va_end(args);
}
}
-/// Support of debug messages.
-void Geant4Action::debug(const string& fmt, ...) const {
- va_list args;
- va_start(args, fmt);
- DD4hep::printout(DD4hep::DEBUG, "Geant4Action", fmt, args);
- va_end(args);
+/// Support for messages with variable output level using output level-2
+void Geant4Action::printM2(const char* fmt, ...) const {
+ int level = max(outputLevel()-2,(int)VERBOSE);
+ if ( level >= printLevel() ) {
+ va_list args;
+ va_start(args, fmt);
+ DD4hep::printout((PrintLevel)level, m_name.c_str(), fmt, args);
+ va_end(args);
+ }
}
-/// Support for messages with variable output level using output level
-void Geant4Action::debug(const std::string& tag, const std::string& fmt, ...) const {
- va_list args;
- va_start(args, fmt);
- DD4hep::printout(DD4hep::DEBUG, tag, fmt, args);
- va_end(args);
+/// Support for messages with variable output level using output level-1
+void Geant4Action::printP1(const char* fmt, ...) const {
+ int level = min(outputLevel()+1,(int)FATAL);
+ if ( level >= printLevel() ) {
+ va_list args;
+ va_start(args, fmt);
+ DD4hep::printout((PrintLevel)level, m_name.c_str(), fmt, args);
+ va_end(args);
+ }
}
-/// Support of info messages.
-void Geant4Action::info(const string& fmt, ...) const {
- va_list args;
- va_start(args, fmt);
- DD4hep::printout(DD4hep::INFO, "Geant4Action", fmt, args);
- va_end(args);
+/// Support for messages with variable output level using output level-2
+void Geant4Action::printP2(const char* fmt, ...) const {
+ int level = min(outputLevel()+2,(int)FATAL);
+ if ( level >= printLevel() ) {
+ va_list args;
+ va_start(args, fmt);
+ DD4hep::printout((PrintLevel)level, m_name.c_str(), fmt, args);
+ va_end(args);
+ }
}
-/// Support for messages with variable output level using output level
-void Geant4Action::info(const std::string& tag, const std::string& fmt, ...) const {
+/// Support of debug messages.
+void Geant4Action::debug(const char* fmt, ...) const {
va_list args;
va_start(args, fmt);
- DD4hep::printout(DD4hep::INFO, tag, fmt, args);
+ DD4hep::printout(DD4hep::DEBUG, "Geant4Action", fmt, args);
va_end(args);
}
-/// Support of warning messages.
-void Geant4Action::warning(const string& fmt, ...) const {
+/// Support of info messages.
+void Geant4Action::info(const char* fmt, ...) const {
va_list args;
va_start(args, fmt);
- DD4hep::printout(DD4hep::WARNING, "Geant4Action", fmt, args);
+ DD4hep::printout(DD4hep::INFO, "Geant4Action", fmt, args);
va_end(args);
}
-/// Support for messages with variable output level using output level
-void Geant4Action::warning(const std::string& tag, const std::string& fmt, ...) const {
+/// Support of warning messages.
+void Geant4Action::warning(const char* fmt, ...) const {
va_list args;
va_start(args, fmt);
- DD4hep::printout(DD4hep::WARNING, tag, fmt, args);
+ DD4hep::printout(DD4hep::WARNING, "Geant4Action", fmt, args);
va_end(args);
}
/// Action to support error messages.
-void Geant4Action::error(const string& fmt, ...) const {
+void Geant4Action::error(const char* fmt, ...) const {
va_list args;
va_start(args, fmt);
DD4hep::printout(DD4hep::ERROR, "Geant4Action", fmt, args);
va_end(args);
}
-/// Support for messages with variable output level using output level
-void Geant4Action::error(const std::string& tag, const std::string& fmt, ...) const {
- va_list args;
- va_start(args, fmt);
- DD4hep::printout(DD4hep::ERROR, tag, fmt, args);
- va_end(args);
-}
-
/// Action to support error messages.
-bool Geant4Action::error(bool return_value, const string& fmt, ...) const {
+bool Geant4Action::error(bool return_value, const char* fmt, ...) const {
va_list args;
va_start(args, fmt);
DD4hep::printout(DD4hep::ERROR, "Geant4Action", fmt, args);
@@ -237,7 +238,7 @@ bool Geant4Action::error(bool return_value, const string& fmt, ...) const {
}
/// Support of fatal messages. Throws exception if required.
-void Geant4Action::fatal(const string& fmt, ...) const {
+void Geant4Action::fatal(const char* fmt, ...) const {
string err;
va_list args;
va_start(args, fmt);
@@ -245,16 +246,8 @@ void Geant4Action::fatal(const string& fmt, ...) const {
va_end(args);
}
-/// Support for messages with variable output level using output level
-void Geant4Action::fatal(const std::string& tag, const std::string& fmt, ...) const {
- va_list args;
- va_start(args, fmt);
- DD4hep::printout(DD4hep::FATAL, tag, fmt, args);
- va_end(args);
-}
-
/// Support of exceptions: Print fatal message and throw runtime_error.
-void Geant4Action::except(const string& fmt, ...) const {
+void Geant4Action::except(const char* fmt, ...) const {
string err;
va_list args;
va_start(args, fmt);
@@ -264,16 +257,8 @@ void Geant4Action::except(const string& fmt, ...) const {
throw runtime_error(err);
}
-/// Support for messages with variable output level using output level
-void Geant4Action::except(const std::string& tag, const std::string& fmt, ...) const {
- va_list args;
- va_start(args, fmt);
- DD4hep::printout(DD4hep::FATAL, tag, fmt, args);
- va_end(args);
-}
-
/// Abort Geant4 Run by throwing a G4Exception with type RunMustBeAborted
-void Geant4Action::abortRun(const string& exception, const string& fmt, ...) const {
+void Geant4Action::abortRun(const string& exception, const char* fmt, ...) const {
string desc, typ = typeName(typeid(*this));
string issuer = name()+" ["+typ+"]";
va_list args;
diff --git a/DDG4/src/Geant4Context.cpp b/DDG4/src/Geant4Context.cpp
index a31d4b3ab..a3dc1a652 100644
--- a/DDG4/src/Geant4Context.cpp
+++ b/DDG4/src/Geant4Context.cpp
@@ -29,8 +29,8 @@ Geant4Run::~Geant4Run() {
}
/// Intializing constructor
-Geant4Event::Geant4Event(const G4Event* evt)
-: ObjectExtensions(typeid(Geant4Event)), m_event(evt)
+Geant4Event::Geant4Event(const G4Event* evt, Geant4Random* rnd)
+ : ObjectExtensions(typeid(Geant4Event)), m_event(evt), m_random(rnd)
{
InstanceCount::increment(this);
}
diff --git a/DDG4/src/Geant4Data.cpp b/DDG4/src/Geant4Data.cpp
index 0b1c79506..79f54cd37 100644
--- a/DDG4/src/Geant4Data.cpp
+++ b/DDG4/src/Geant4Data.cpp
@@ -8,10 +8,12 @@
//====================================================================
// Framework include files
+#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4Data.h"
// Geant4 include files
+#include "G4Step.hh"
#include "G4Allocator.hh"
#include "G4OpticalPhoton.hh"
@@ -45,92 +47,19 @@ SimpleEvent::~SimpleEvent() {
DataExtension::~DataExtension() {
}
-/// Copy constructor
-Particle::Particle(const Particle& c)
- : id(c.id), g4Parent(c.g4Parent), parent(c.parent), reason(c.reason),
- steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID),
- vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
- vex(c.vex), vey(c.vey), vez(c.vez),
- psx(c.psx), psy(c.psy), psz(c.psz),
- pex(c.pex), pey(c.pey), pez(c.pez),
- energy(c.energy), time(c.time),
- daughters(c.daughters), extension(),
- process(c.process), definition(c.definition)
-{
- InstanceCount::increment(this);
-}
-
-/// Default constructor
-Particle::Particle()
- : id(0), g4Parent(0), parent(0), reason(0),
- steps(0), secondaries(0), pdgID(0),
- vsx(0.0), vsy(0.0), vsz(0.0),
- vex(0.0), vey(0.0), vez(0.0),
- psx(0.0), psy(0.0), psz(0.0),
- pex(0.0), pey(0.0), pez(0.0),
- energy(0.0), time(0.0),
- daughters(), extension(), process(0), definition(0)
-{
- InstanceCount::increment(this);
-}
-
-/// Default destructor
-Particle::~Particle() {
- InstanceCount::decrement(this);
-}
-
-/// Assignment operator
-Particle& Particle::get_data(Particle& c) {
- if ( this != &c ) {
- id = c.id;
- g4Parent = c.g4Parent;
- parent = c.parent;
- reason = c.reason;
- steps = c.steps;
- secondaries = c.secondaries;
- pdgID = c.pdgID;
- vsx = c.vsx;
- vsy = c.vsy;
- vsz = c.vsz;
- vex = c.vex;
- vey = c.vey;
- vez = c.vez;
- psx = c.psx;
- psy = c.psy;
- psz = c.psz;
- pex = c.pex;
- pey = c.pey;
- pez = c.pez;
- energy = c.energy;
- time = c.time;
- process = c.process;
- definition = c.definition;
- daughters = c.daughters;
- extension = c.extension;
- //auto_ptr<DataExtension>(c.extension.release());
- }
- return *this;
-}
-
-/// Remove daughter from set
-void Particle::removeDaughter(int id_daughter) {
- set<int>::iterator j = daughters.find(id_daughter);
- if ( j != daughters.end() ) daughters.erase(j);
-}
-
/// Default constructor
-SimpleHit::SimpleHit()
+Geant4HitData::Geant4HitData()
: cellID(0), extension() {
InstanceCount::increment(this);
}
/// Default destructor
-SimpleHit::~SimpleHit() {
+Geant4HitData::~Geant4HitData() {
InstanceCount::decrement(this);
}
/// Extract the MC contribution for a given hit from the step information
-SimpleHit::Contribution SimpleHit::extractContribution(G4Step* step) {
+Geant4HitData::Contribution Geant4HitData::extractContribution(G4Step* step) {
G4Track* trk = step->GetTrack();
double energy_deposit =
(trk->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition()) ?
@@ -140,24 +69,24 @@ SimpleHit::Contribution SimpleHit::extractContribution(G4Step* step) {
}
/// Default constructor
-SimpleTracker::Hit::Hit()
- : SimpleHit(), position(), momentum(), length(0.0), truth(), energyDeposit(0.0) {
+Geant4Tracker::Hit::Hit()
+ : Geant4HitData(), position(), momentum(), length(0.0), truth(), energyDeposit(0.0) {
InstanceCount::increment(this);
}
/// Standard initializing constructor
-SimpleTracker::Hit::Hit(int track_id, int pdg_id, double deposit, double time_stamp)
- : SimpleHit(), position(), momentum(), length(0.0), truth(track_id, pdg_id, deposit, time_stamp), energyDeposit(deposit) {
+Geant4Tracker::Hit::Hit(int track_id, int pdg_id, double deposit, double time_stamp)
+ : Geant4HitData(), position(), momentum(), length(0.0), truth(track_id, pdg_id, deposit, time_stamp), energyDeposit(deposit) {
InstanceCount::increment(this);
}
/// Default destructor
-SimpleTracker::Hit::~Hit() {
+Geant4Tracker::Hit::~Hit() {
InstanceCount::decrement(this);
}
/// Assignment operator
-SimpleTracker::Hit& SimpleTracker::Hit::operator=(const Hit& c) {
+Geant4Tracker::Hit& Geant4Tracker::Hit::operator=(const Hit& c) {
if ( &c != this ) {
position = c.position;
momentum = c.momentum;
@@ -168,7 +97,7 @@ SimpleTracker::Hit& SimpleTracker::Hit::operator=(const Hit& c) {
}
/// Clear hit content
-SimpleTracker::Hit& SimpleTracker::Hit::clear() {
+Geant4Tracker::Hit& Geant4Tracker::Hit::clear() {
position.SetXYZ(0, 0, 0);
momentum.SetXYZ(0, 0, 0);
length = 0.0;
@@ -177,7 +106,7 @@ SimpleTracker::Hit& SimpleTracker::Hit::clear() {
}
/// Store Geant4 point and step information into tracker hit structure.
-SimpleTracker::Hit& SimpleTracker::Hit::storePoint(G4Step* step, G4StepPoint* pnt) {
+Geant4Tracker::Hit& Geant4Tracker::Hit::storePoint(G4Step* step, G4StepPoint* pnt) {
G4Track* trk = step->GetTrack();
G4ThreeVector pos = pnt->GetPosition();
G4ThreeVector mom = pnt->GetMomentum();
@@ -193,18 +122,18 @@ SimpleTracker::Hit& SimpleTracker::Hit::storePoint(G4Step* step, G4StepPoint* pn
}
/// Default constructor (for ROOT)
-SimpleCalorimeter::Hit::Hit()
- : SimpleHit(), position(), truth(), energyDeposit(0) {
+Geant4Calorimeter::Hit::Hit()
+ : Geant4HitData(), position(), truth(), energyDeposit(0) {
InstanceCount::increment(this);
}
/// Standard constructor
-SimpleCalorimeter::Hit::Hit(const Position& pos)
- : SimpleHit(), position(pos), truth(), energyDeposit(0) {
+Geant4Calorimeter::Hit::Hit(const Position& pos)
+ : Geant4HitData(), position(pos), truth(), energyDeposit(0) {
InstanceCount::increment(this);
}
/// Default destructor
-SimpleCalorimeter::Hit::~Hit() {
+Geant4Calorimeter::Hit::~Hit() {
InstanceCount::decrement(this);
}
diff --git a/DDG4/src/Geant4DataConversion.cpp b/DDG4/src/Geant4DataConversion.cpp
index 99e68301c..0ca9c6a04 100644
--- a/DDG4/src/Geant4DataConversion.cpp
+++ b/DDG4/src/Geant4DataConversion.cpp
@@ -20,9 +20,21 @@ Geant4ConversionHelper::Geant4ConversionHelper() {
Geant4ConversionHelper::~Geant4ConversionHelper() {
}
+/// Access to the data encoding using the volume manager and a specified volume id
std::string Geant4ConversionHelper::encoding(Geometry::VolumeManager vm, Geometry::VolumeManager::VolumeID vid) {
Geometry::PlacedVolume pv = vm.lookupPlacement(vid);
Geometry::SensitiveDetector sd = pv.volume().sensitiveDetector();
- Geometry::IDDescriptor id = sd.readout().idSpec();
+ return encoding(sd);
+}
+
+/// Access to the hit encoding in this sensitive detector
+std::string Geant4ConversionHelper::encoding(Geometry::Handle<Geometry::SensitiveDetectorObject> sd) {
+ Geometry::IDDescriptor id = Geometry::SensitiveDetector(sd).readout().idSpec();
+ return id.fieldDescription();
+}
+
+/// Access to the hit encoding in this readout object
+std::string Geant4ConversionHelper::encoding(Geometry::Readout ro) {
+ Geometry::IDDescriptor id = ro.idSpec();
return id.fieldDescription();
}
diff --git a/DDG4/src/Geant4DataDump.cpp b/DDG4/src/Geant4DataDump.cpp
index 176dcdae3..ad55dfd1c 100644
--- a/DDG4/src/Geant4DataDump.cpp
+++ b/DDG4/src/Geant4DataDump.cpp
@@ -26,17 +26,18 @@ Geant4DataDump::~Geant4DataDump() {
}
/// Print a single particle to the output logging using the specified print level
-void Geant4DataDump::print(PrintLevel level, const Particle* p) const {
+void Geant4DataDump::print(PrintLevel level, Geant4ParticleHandle p) const {
PropertyMask mask(p->reason);
+ int parent = p->parents.empty() ? -1 : *p->parents.begin();
printout(level, m_tag, " +++ TrackID: %6d %12d %6d %-7s %3s %5d %6s %8.3g %-4s %-7s %-7s %-3s",
p->id,
p->pdgID,
- p->parent,
+ parent,
yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
int(p->daughters.size()),
yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
- p->energy,
+ p.energy(),
yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
diff --git a/DDG4/src/Geant4Exec.cpp b/DDG4/src/Geant4Exec.cpp
index a37dbfa2b..d3d160a48 100644
--- a/DDG4/src/Geant4Exec.cpp
+++ b/DDG4/src/Geant4Exec.cpp
@@ -28,8 +28,9 @@
#include "DDG4/Geant4StackingAction.h"
#include "DDG4/Geant4GeneratorAction.h"
#include "DDG4/Geant4PhysicsList.h"
-#include "DDG4/Geant4Kernel.h"
#include "DDG4/Geant4UIManager.h"
+#include "DDG4/Geant4Kernel.h"
+#include "DDG4/Geant4Random.h"
#include <memory>
#include <stdexcept>
@@ -46,6 +47,7 @@ namespace DD4hep {
namespace {
Geant4Context* s_globalContext = 0;
+ Geant4Random* s_globalRandom = 0;
}
Geant4Context* ddg4_globalContext() {
@@ -98,7 +100,7 @@ namespace DD4hep {
}
}
void createClientContext(const G4Event* evt) {
- Geant4Event* e = new Geant4Event(evt);
+ Geant4Event* e = new Geant4Event(evt,s_globalRandom);
m_activeContext->setEvent(e);
setContextToClients();
}
@@ -291,6 +293,8 @@ int Geant4Exec::configure(Geant4Kernel& kernel) {
CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine);
Geometry::LCDD& lcdd = kernel.lcdd();
Geant4Context* ctx = s_globalContext = new Geant4Context(&kernel);
+ // For now do this:
+ /* Geant4Random* rnd = */ s_globalRandom = new Geant4Random();
// Construct the default run manager
G4RunManager& runManager = kernel.runManager();
diff --git a/DDG4/src/Geant4GeneratorActionInit.cpp b/DDG4/src/Geant4GeneratorActionInit.cpp
new file mode 100644
index 000000000..cb3560d83
--- /dev/null
+++ b/DDG4/src/Geant4GeneratorActionInit.cpp
@@ -0,0 +1,81 @@
+// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+// Framework include files
+#include "DD4hep/Printout.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Kernel.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4RunAction.h"
+#include "DDG4/Geant4GeneratorActionInit.h"
+
+#include "G4Run.hh"
+
+using namespace DD4hep::Simulation;
+
+/// Standard constructor
+Geant4GeneratorActionInit::Geant4GeneratorActionInit(Geant4Context* context, const std::string& nam)
+ : Geant4GeneratorAction(context,nam), m_run(0), m_evtTotal(0), m_evtRun(0)
+{
+ InstanceCount::increment(this);
+ context->kernel().runAction().callAtEnd(this,&Geant4GeneratorActionInit::end);
+ context->kernel().runAction().callAtBegin(this,&Geant4GeneratorActionInit::begin);
+}
+
+/// Default destructor
+Geant4GeneratorActionInit::~Geant4GeneratorActionInit() {
+ InstanceCount::decrement(this);
+}
+
+/// Begin-run action callback
+void Geant4GeneratorActionInit::begin(const G4Run* run) {
+ m_evtRun = 0;
+ m_run = run->GetRunID();
+}
+
+/// End-run action callback
+void Geant4GeneratorActionInit::end(const G4Run* /* run */) {
+ printP1("+++ Funished run %d after %d events (% events in total)",m_run,m_evtRun,m_evtTotal);
+ m_evtRun = 0;
+ m_run = 0;
+}
+
+/// Event generation action callback
+void Geant4GeneratorActionInit::operator()(G4Event* /* event */) {
+ /// Update event counters
+ ++m_evtTotal;
+ ++m_evtRun;
+ /// + Printout
+ print("+++ Initializing event %d. Within run:%d event %d.",m_evtTotal,m_run,m_evtRun);
+
+ /**
+ * This action should register all event extension required for the further
+ * processing. We want to avoid that every client has to check if a given
+ * object is present or not and than later install the required data structures.
+ */
+ context()->event().addExtension(new Geant4PrimaryMap());
+
+ // The final set of created particles in the simulation.
+ context()->event().addExtension(new Geant4ParticleMap());
+
+ //
+ // The Geant4PrimaryEvent extension contains a whole set of
+ // Geant4PrimaryInteraction objects each may represent a complete
+ // interaction. Particles and vertices may be unbiased.
+ // This is the input to the translator forming the final
+ // Geant4PrimaryInteraction (see below) containing rebiased particle
+ // and vertex maps.
+ Geant4PrimaryEvent* evt = new Geant4PrimaryEvent();
+ context()->event().addExtension(evt);
+ //
+ // The Geant4PrimaryInteraction extension contains the final
+ // vertices and particles ready to be injected to Geant4.
+ Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
+ inter->setNextPID(-1);
+ context()->event().addExtension(inter);
+}
diff --git a/DDG4/src/Geant4HitCollection.cpp b/DDG4/src/Geant4HitCollection.cpp
index b1a1b4604..b8e110c07 100644
--- a/DDG4/src/Geant4HitCollection.cpp
+++ b/DDG4/src/Geant4HitCollection.cpp
@@ -75,6 +75,16 @@ Geant4HitCollection::~Geant4HitCollection() {
InstanceCount::decrement(this);
}
+/// Set the sensitive detector
+void Geant4HitCollection::setSensitiveDetector(SensitiveDetector detector) {
+ m_detector = detector;
+}
+
+/// Access the sensitive detector
+Geant4HitCollection::SensitiveDetector Geant4HitCollection::sensitiveDetector() const {
+ return m_detector;
+}
+
/// Type information of the object stored
const ComponentCast& Geant4HitCollection::type() const {
return m_manipulator->cast;
diff --git a/DDG4/src/Geant4InteractionMerger.cpp b/DDG4/src/Geant4InteractionMerger.cpp
new file mode 100644
index 000000000..e8b4489ae
--- /dev/null
+++ b/DDG4/src/Geant4InteractionMerger.cpp
@@ -0,0 +1,93 @@
+// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $
+//====================================================================
+// 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/Geant4InteractionMerger.h"
+
+#include <stdexcept>
+
+using namespace std;
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+/// Standard constructor
+Geant4InteractionMerger::Geant4InteractionMerger(Geant4Context* context, const std::string& nam)
+ : Geant4GeneratorAction(context,nam)
+{
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4InteractionMerger::~Geant4InteractionMerger() {
+ InstanceCount::decrement(this);
+}
+
+void rebaseParticles(Geant4PrimaryInteraction::ParticleMap& particles, int &offset) {
+ Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
+ int mx_id = offset;
+ // Now move begin and end-vertex of all primary and generator particles accordingly
+ for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip ) {
+ Geant4ParticleHandle p((*ip).second);
+ p.offset(offset);
+ mx_id = p->id+1 > mx_id ? p->id+1 : mx_id;
+ }
+ offset = mx_id;
+}
+
+void rebaseVertices(Geant4PrimaryInteraction::VertexMap& vertices, int part_offset) {
+ Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
+ set<int> in, out;
+ set<int>::iterator i;
+ // Now move begin and end-vertex of all primary vertices accordingly
+ for(iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv) {
+ Geant4Vertex* v = (*iv).second;
+ in = v->in;
+ out = v->out;
+ for(in=v->in, i=in.begin(), v->in.clear(); i != in.end(); ++i)
+ v->in.insert((*i)+part_offset);
+ for(out=v->out, i=out.begin(), v->out.clear(); i != out.end(); ++i)
+ v->out.insert((*i)+part_offset);
+ }
+}
+
+void 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()));
+}
+
+/// Event generation action callback
+void Geant4InteractionMerger::operator()(G4Event* /* event */) {
+ typedef Geant4PrimaryEvent::Interaction Interaction;
+ typedef vector<Interaction*> Interactions;
+ Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
+ Interaction* output = context()->event().extension<Interaction>();
+ Interactions inter = evt->interactions();
+ int particle_offset = 0, vertex_offset = 0;
+
+ for(Interactions::const_iterator i=inter.begin(); i != inter.end(); ++i) {
+ Interaction* interaction = *i;
+ vertex_offset = particle_offset;
+ rebaseParticles(interaction->particles,particle_offset);
+ rebaseVertices(interaction->vertices,vertex_offset);
+ appendInteraction(output,interaction);
+ }
+ output->setNextPID(particle_offset);
+ Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
+ debug("+++ Merged MC input record from %d interactions:",(int)inter.size());
+ for( ip=output->particles.begin(), ipend=output->particles.end(); ip != ipend; ++ip )
+ Geant4ParticleHandle((*ip).second).dump1(DEBUG,name(),"Merged particles");
+}
diff --git a/DDG4/src/Geant4InteractionVertexBoost.cpp b/DDG4/src/Geant4InteractionVertexBoost.cpp
new file mode 100644
index 000000000..d9b889cd1
--- /dev/null
+++ b/DDG4/src/Geant4InteractionVertexBoost.cpp
@@ -0,0 +1,99 @@
+// $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/Geant4Random.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4InteractionVertexBoost.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+
+#include "G4ParticleDefinition.hh"
+
+// C/C++ include files
+#include <stdexcept>
+#include <cmath>
+
+using namespace std;
+using namespace DD4hep::Simulation;
+#define SQR(x) (x*x)
+
+/// Standard constructor
+Geant4InteractionVertexBoost::Geant4InteractionVertexBoost(Geant4Context* context, const string& name)
+ : Geant4GeneratorAction(context, name)
+{
+ InstanceCount::increment(this);
+ declareProperty("Angle", m_angle = 0);
+ declareProperty("Mask", m_mask = 0);
+ m_needsControl = true;
+}
+
+/// Default destructor
+Geant4InteractionVertexBoost::~Geant4InteractionVertexBoost() {
+ InstanceCount::decrement(this);
+}
+
+/// Callback to generate primary particles
+void Geant4InteractionVertexBoost::operator()(G4Event*) {
+ typedef Geant4PrimaryEvent::Interaction Interaction;
+ typedef Geant4Particle Particle;
+ typedef Geant4Particle Vertex;
+ Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
+ Interaction* inter = evt->get(m_mask);
+
+ if ( inter ) {
+ Interaction::VertexMap::iterator iv;
+ Interaction::ParticleMap::iterator ip;
+ double alpha = m_angle;
+ double gamma = std::sqrt(1 + SQR(tan(alpha)));
+ double betagamma = std::tan(alpha);
+
+ if ( alpha != 0.0 ) {
+ // Now move begin and end-vertex of all primary vertices accordingly
+ for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
+ Geant4Vertex* v = (*iv).second;
+ double t = gamma * t + betagamma * v->x / c_light;
+ double x = gamma * v->x + betagamma * c_light * v->time;
+ double y = v->y;
+ double z = v->z;
+ v->x = x;
+ v->y = y;
+ v->z = z;
+ v->time = t;
+ }
+
+ // Now move begin and end-vertex of all primary and generator particles accordingly
+ for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
+ Particle* p = (*ip).second;
+ double t = gamma * p->time + betagamma * p->vsx / c_light;
+ double x = gamma * p->vsx + betagamma * c_light * p->time;
+ double y = p->vsx;
+ double z = p->vsz;
+
+ double m = p->definition->GetPDGMass();
+ double e2 = SQR(p->psx)+SQR(p->psy)+SQR(p->psz)+SQR(m);
+ double px = betagamma * std::sqrt(e2) + gamma * p->psx;
+ double py = p->psy;
+ double pz = p->psz;
+
+ p->vsx = x;
+ p->vsy = y;
+ p->vsz = z;
+ p->time = t;
+
+ p->psx += px;
+ p->psy += py;
+ p->psz += pz;
+ }
+ }
+ }
+}
diff --git a/DDG4/src/Geant4InteractionVertexSmear.cpp b/DDG4/src/Geant4InteractionVertexSmear.cpp
new file mode 100644
index 000000000..6a2c562f0
--- /dev/null
+++ b/DDG4/src/Geant4InteractionVertexSmear.cpp
@@ -0,0 +1,84 @@
+// $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/Geant4Random.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4InteractionVertexSmear.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+// C/C++ include files
+#include <stdexcept>
+#include <cmath>
+
+using namespace std;
+using namespace DD4hep::Simulation;
+
+/// Standard constructor
+Geant4InteractionVertexSmear::Geant4InteractionVertexSmear(Geant4Context* context, const string& name)
+ : Geant4GeneratorAction(context, name)
+{
+ InstanceCount::increment(this);
+ declareProperty("Offset", m_offset);
+ declareProperty("Sigma", m_sigma);
+ declareProperty("Mask", m_mask = 0);
+ m_needsControl = true;
+}
+
+/// Default destructor
+Geant4InteractionVertexSmear::~Geant4InteractionVertexSmear() {
+ InstanceCount::decrement(this);
+}
+
+/// Callback to generate primary particles
+void Geant4InteractionVertexSmear::operator()(G4Event*) {
+ typedef Geant4PrimaryEvent::Interaction Interaction;
+ typedef Geant4Particle Particle;
+ typedef Geant4Particle Vertex;
+ Geant4Random& rndm = context()->event().random();
+ Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
+ Interaction* inter = evt->get(m_mask);
+
+ if ( inter ) {
+ Interaction::VertexMap::iterator iv;
+ Interaction::ParticleMap::iterator ip;
+ double dx = rndm.gauss(m_offset.x(),m_sigma.x());
+ double dy = rndm.gauss(m_offset.y(),m_sigma.y());
+ 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);
+
+ // Now move begin and end-vertex of all primary vertices accordingly
+ for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
+ Geant4Vertex* v = (*iv).second;
+ v->x += dx;
+ v->y += dy;
+ v->z += dz;
+ v->time += dt;
+ }
+ // Now move begin and end-vertex of all primary and generator particles accordingly
+ for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
+ Particle* p = (*ip).second;
+ p->vsx += dx;
+ p->vsy += dy;
+ p->vsz += dz;
+ p->vex += dx;
+ p->vey += dy;
+ p->vez += dz;
+ p->time += dt;
+ }
+ }
+ else {
+ print("+++ No interaction of type %d present.",m_mask);
+ }
+}
diff --git a/DDG4/src/Geant4IsotropeGenerator.cpp b/DDG4/src/Geant4IsotropeGenerator.cpp
new file mode 100644
index 000000000..bd6d4ef64
--- /dev/null
+++ b/DDG4/src/Geant4IsotropeGenerator.cpp
@@ -0,0 +1,92 @@
+// $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/Geant4Random.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4IsotropeGenerator.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+
+#include "G4ParticleTable.hh"
+#include "G4ParticleDefinition.hh"
+
+// C/C++ include files
+#include <stdexcept>
+#include <cmath>
+
+using namespace std;
+using namespace DD4hep::Simulation;
+
+/// Standard constructor
+Geant4IsotropeGenerator::Geant4IsotropeGenerator(Geant4Context* context, const string& name)
+ : Geant4GeneratorAction(context, name), m_particle(0)
+{
+ 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);
+}
+
+/// Default destructor
+Geant4IsotropeGenerator::~Geant4IsotropeGenerator() {
+ InstanceCount::decrement(this);
+}
+
+/// Callback to generate primary particles
+void Geant4IsotropeGenerator::operator()(G4Event*) {
+ typedef Geant4Particle Particle;
+
+ if (0 == m_particle || m_particle->GetParticleName() != m_particleName.c_str()) {
+ G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
+ m_particle = particleTable->FindParticle(m_particleName);
+ if (0 == m_particle) {
+ throw runtime_error("Geant4IsotropeGenerator: Bad particle type:"+m_particleName+"!");
+ }
+ }
+ Geant4Event& evt = context()->event();
+ Geant4Random& rnd = evt.random();
+ Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
+ Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
+ prim->add(m_mask, inter);
+
+ Geant4Vertex* vtx = new Geant4Vertex();
+ inter->vertices.insert(make_pair(inter->vertices.size(),vtx));
+ for(int i=0; i<m_multiplicity; ++i) {
+ double phi = 2*M_PI*rnd.rndm();
+ double theta = M_PI*rnd.rndm();
+ double x1 = std::sin(theta)*std::cos(phi);
+ double x2 = std::sin(theta)*std::sin(phi);
+ double x3 = std::cos(theta);
+ double momentum = rnd.rndm()*m_energy;
+
+ 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->pdgID = m_particle->GetPDGEncoding();
+ p->definition = m_particle;
+ p->psx = x1*momentum;
+ p->psy = x2*momentum;
+ p->psz = x3*momentum;
+ p->mass = m_particle->GetPDGMass();
+ inter->particles.insert(make_pair(p->id,p));
+ vtx->out.insert(p->id);
+ printout(INFO,name(),"Particle [%d] %s %.3f GeV direction:(%6.3f %6.3f %6.3f)",
+ p->id, m_particleName.c_str(), momentum/GeV, x1, x2, x3);
+
+ }
+}
diff --git a/DDG4/src/Geant4Output2ROOT.cpp b/DDG4/src/Geant4Output2ROOT.cpp
index 118465c28..85a933ad7 100644
--- a/DDG4/src/Geant4Output2ROOT.cpp
+++ b/DDG4/src/Geant4Output2ROOT.cpp
@@ -13,7 +13,7 @@
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4HitCollection.h"
#include "DDG4/Geant4Output2ROOT.h"
-#include "DDG4/Geant4MonteCarloTruth.h"
+#include "DDG4/Geant4Particle.h"
#include "DDG4/Geant4Data.h"
// Geant4 include files
#include "G4HCofThisEvent.hh"
@@ -135,12 +135,12 @@ void Geant4Output2ROOT::commit(OutputContext<G4Event>& ctxt) {
/// Callback to store the Geant4 event
void Geant4Output2ROOT::saveEvent(OutputContext<G4Event>& /* ctxt */) {
- Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false);
- if ( truth ) {
+ Geant4ParticleMap* parts = context()->event().extension<Geant4ParticleMap>();
+ if ( parts ) {
typedef Geant4HitWrapper::HitManipulator Manip;
- typedef Geant4MonteCarloTruth::ParticleMap ParticleMap;
- Manip* manipulator = Geant4HitWrapper::manipulator<Particle>();
- const ParticleMap& pm = truth->particles();
+ typedef Geant4ParticleMap::ParticleMap ParticleMap;
+ Manip* manipulator = Geant4HitWrapper::manipulator<Geant4Particle>();
+ const ParticleMap& pm = parts->particles();
vector<void*> particles;
for(ParticleMap::const_iterator i=pm.begin(); i!=pm.end(); ++i) {
particles.push_back((ParticleMap::mapped_type*)(*i).second);
@@ -157,7 +157,7 @@ void Geant4Output2ROOT::saveCollection(OutputContext<G4Event>& /* ctxt */, G4VHi
if (coll) {
coll->getHitsUnchecked(hits);
size_t nhits = coll->GetSize();
- Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false);
+ Geant4ParticleMap* truth = context()->event().extension<Geant4ParticleMap>();
if ( truth && nhits > 0 ) {
try {
for(size_t i=0; i<nhits; ++i) {
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
new file mode 100644
index 000000000..bd2953641
--- /dev/null
+++ b/DDG4/src/Geant4Particle.cpp
@@ -0,0 +1,315 @@
+// $Id: Geant4Hits.cpp 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+
+// Framework include files
+#include "DD4hep/Printout.h"
+#include "DD4hep/Primitives.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Particle.h"
+#include "TDatabasePDG.h"
+#include "TParticlePDG.h"
+#include "G4ParticleTable.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4VProcess.hh"
+
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+/// Default destructor
+ParticleExtension::~ParticleExtension() {
+}
+
+/// Copy constructor
+Geant4Particle::Geant4Particle(const Geant4Particle& c)
+ : ref(1), id(c.id), g4Parent(c.g4Parent), reason(c.reason), usermask(c.usermask),
+ steps(c.steps), secondaries(c.secondaries), pdgID(c.pdgID),
+ status(c.status), charge(0),
+ vsx(c.vsx), vsy(c.vsy), vsz(c.vsz),
+ vex(c.vex), vey(c.vey), vez(c.vez),
+ psx(c.psx), psy(c.psy), psz(c.psz),
+ pex(c.pex), pey(c.pey), pez(c.pez),
+ mass(c.mass), time(c.time), properTime(c.properTime),
+ parents(c.parents), daughters(c.daughters), extension(),
+ process(c.process), definition(c.definition)
+{
+ InstanceCount::increment(this);
+ spin[0] = c.spin[0];
+ spin[1] = c.spin[1];
+ spin[2] = c.spin[2];
+ colorFlow[0] = c.colorFlow[0];
+ colorFlow[1] = c.colorFlow[1];
+}
+
+/// Default constructor
+Geant4Particle::Geant4Particle()
+ : ref(1), id(0), g4Parent(0), reason(0), usermask(0),
+ steps(0), secondaries(0), pdgID(0),
+ status(0), charge(0),
+ vsx(0.0), vsy(0.0), vsz(0.0),
+ vex(0.0), vey(0.0), vez(0.0),
+ psx(0.0), psy(0.0), psz(0.0),
+ pex(0.0), pey(0.0), pez(0.0),
+ mass(0.0), time(0.0), properTime(0.0),
+ daughters(), extension(), process(0), definition(0)
+{
+ InstanceCount::increment(this);
+ spin[0] = spin[1] = spin[2] = 0;
+ colorFlow[0] = colorFlow[1] = 0;
+}
+
+/// Constructor with ID initialization
+Geant4Particle::Geant4Particle(int part_id)
+ : ref(1), id(part_id), g4Parent(0), reason(0), usermask(0),
+ steps(0), secondaries(0), pdgID(0),
+ status(0), charge(0),
+ vsx(0.0), vsy(0.0), vsz(0.0),
+ vex(0.0), vey(0.0), vez(0.0),
+ psx(0.0), psy(0.0), psz(0.0),
+ pex(0.0), pey(0.0), pez(0.0),
+ mass(0.0), time(0.0), properTime(0.0),
+ daughters(), extension(), process(0), definition(0)
+{
+ InstanceCount::increment(this);
+ spin[0] = spin[1] = spin[2] = 0;
+ colorFlow[0] = colorFlow[1] = 0;
+}
+
+/// Default destructor
+Geant4Particle::~Geant4Particle() {
+ InstanceCount::decrement(this);
+ //::printf("************ Delete Geant4Particle[%p]: ID:%d pdgID %d ref:%d\n",(void*)this,id,pdgID,ref);
+}
+
+void Geant4Particle::release() {
+ //::printf("************ Release Geant4Particle[%p]: ID:%d pdgID %d ref:%d\n",(void*)this,id,pdgID,ref-1);
+ if ( --ref <= 0 ) {
+ delete this;
+ }
+}
+
+/// Assignment operator
+Geant4Particle& Geant4Particle::get_data(Geant4Particle& c) {
+ if ( this != &c ) {
+ id = c.id;
+ g4Parent = c.g4Parent;
+ reason = c.reason;
+ usermask = c.usermask;
+ status = c.status;
+ charge = c.charge;
+ steps = c.steps;
+ secondaries = c.secondaries;
+ pdgID = c.pdgID;
+ vsx = c.vsx;
+ vsy = c.vsy;
+ vsz = c.vsz;
+ vex = c.vex;
+ vey = c.vey;
+ vez = c.vez;
+ psx = c.psx;
+ psy = c.psy;
+ psz = c.psz;
+ pex = c.pex;
+ pey = c.pey;
+ pez = c.pez;
+ mass = c.mass;
+ time = c.time;
+ properTime = c.properTime;
+ process = c.process;
+ definition = c.definition;
+ daughters = c.daughters;
+ parents = c.parents;
+ extension = c.extension;
+ //auto_ptr<ParticleExtension>(c.extension.release());
+ }
+ return *this;
+}
+
+/// Remove daughter from set
+void Geant4Particle::removeDaughter(int id_daughter) {
+ std::set<int>::iterator j = daughters.find(id_daughter);
+ if ( j != daughters.end() ) daughters.erase(j);
+}
+
+/// Access to the Geant4 particle name
+std::string Geant4ParticleHandle::particleName() const {
+ if ( particle->definition ) return particle->definition->GetParticleName();
+ G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
+ G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
+ if ( def ) {
+ particle->definition = def;
+ return def->GetParticleName();
+ }
+#if 0
+ TDatabasePDG* db = TDatabasePDG::Instance();
+ TParticlePDG* pdef = db->GetParticle(particle->pdgID);
+ if ( pdef ) return pdef->GetName();
+#endif
+ char text[32];
+ ::snprintf(text,sizeof(text),"PDG:%d",particle->pdgID);
+ return text;
+}
+
+/// Access to the Geant4 particle type
+std::string Geant4ParticleHandle::particleType() const {
+ if ( particle->definition ) return particle->definition->GetParticleType();
+ G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
+ G4ParticleDefinition* def = tab->FindParticle(particle->pdgID);
+ if ( def ) {
+ particle->definition = def;
+ return def->GetParticleType();
+ }
+#if 0
+ TDatabasePDG* db = TDatabasePDG::Instance();
+ TParticlePDG* pdef = db->GetParticle(particle->pdgID);
+ if ( pdef ) return pdef->ParticleClass();
+#endif
+ char text[32];
+ ::snprintf(text,sizeof(text),"PDG:%d",particle->pdgID);
+ return text;
+}
+
+/// Access to the creator process name
+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";
+ return "???";
+}
+
+/// Access to the creator process type name
+std::string Geant4ParticleHandle::processTypeName() const {
+ if ( particle->process ) {
+ return G4VProcess::GetProcessTypeName(particle->process->GetProcessType());
+ }
+ return "";
+}
+
+/// Offset all particle identifiers (id, parents, daughters) by a constant number
+void Geant4ParticleHandle::offset(int off) const {
+ std::set<int> temp;
+ Geant4Particle* p = particle;
+ p->id += off;
+
+ temp = p->daughters;
+ p->daughters.clear();
+ for(std::set<int>::iterator i=temp.begin(); i != temp.end(); ++i)
+ p->daughters.insert((*i)+off);
+
+ temp = p->parents;
+ p->parents.clear();
+ for(std::set<int>::iterator i=temp.begin(); i != temp.end(); ++i)
+ p->parents.insert((*i)+off);
+}
+
+/// Output type 1:+++ <tag> 10 def:0xde4eaa8 [gamma , gamma] reason: 20 E:+1.017927e+03 #Par: 1/4 #Dau: 2
+void Geant4ParticleHandle::dump1(int level, const std::string& src, const char* tag) const {
+ char text[256];
+ Geant4ParticleHandle p(*this);
+ text[0]=0;
+ if ( p->parents.size() == 1 )
+ ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
+ else if ( p->parents.size() > 1 ) {
+ text[0]='/';text[1]=0;
+ for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
+ ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
+ }
+ printout((DD4hep::PrintLevel)level,src,
+ "+++ %s %4d def [%-11s,%8s] reason:%8d E:%+.2e %3s #Dau:%3d #Par:%3d%-5s",
+ tag, p->id,
+ p.particleName().c_str(),
+ p.particleType().c_str(),
+ p->reason,
+ p.energy(),
+ p->g4Parent>0 ? "Sim" : "Gen",
+ int(p->daughters.size()),
+ int(p->parents.size()),text);
+}
+
+/// Output type 2:+++ <tag> 20 G4: 7 def:0xde4eaa8 [gamma , gamma] reason: 20 E:+3.304035e+01 in record:YES #Par: 1/18 #Dau: 0
+void Geant4ParticleHandle::dump2(int level, const std::string& src, const char* tag, int g4id, bool inrec) const {
+ char text[32];
+ Geant4ParticleHandle p(*this);
+ if ( p->parents.size() == 0 ) text[0]=0;
+ 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 "
+ "E:%+.2e in record:%s #Par:%3d%-5s #Dau:%3d",
+ tag, p->id, g4id,
+ p.particleName().c_str(),
+ p.particleType().c_str(),
+ p->reason,
+ p.energy(),
+ yes_no(inrec),
+ int(p->parents.size()),text,
+ int(p->daughters.size()));
+}
+
+/// 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
+void Geant4ParticleHandle::dump3(int level, const std::string& src, const char* tag) const {
+ char text[256];
+ Geant4ParticleHandle p(*this);
+ text[0]=0;
+ if ( p->parents.size() == 1 )
+ ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
+ else if ( p->parents.size() > 1 ) {
+ text[0]='/';text[1]=0;
+ for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
+ ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
+ }
+ printout((DD4hep::PrintLevel)level,src,
+ "+++ %s ID:%3d %-12s status:%08X typ:%9d Vtx:(%+.2e,%+.2e,%+.2e)[mm] "
+ "time: %+.2e [ns] #Dau:%3d #Par:%1d%-6s",
+ tag,p->id,p.particleName().c_str(),p->status,p->pdgID,
+ p->vsx/mm,p->vsy/mm,p->vsz/mm,p->time/ns,
+ p->daughters.size(),
+ p->parents.size(),
+ text);
+}
+
+/// Default destructor
+Geant4ParticleMap::~Geant4ParticleMap() {
+ clear();
+}
+
+/// Clear particle maps
+void Geant4ParticleMap::clear() {
+ releaseObjects(particleMap)();
+ particleMap.clear();
+ equivalentTracks.clear();
+}
+
+/// Adopt particle maps
+void Geant4ParticleMap::adopt(ParticleMap& pm, TrackEquivalents& equiv) {
+ clear();
+ particleMap = pm;
+ equivalentTracks = equiv;
+ pm.clear();
+ equiv.clear();
+}
+
+/// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
+int Geant4ParticleMap::particleID(int g4_id, bool) const {
+ int equiv_id = g4_id;
+ if ( g4_id != 0 ) {
+ ParticleMap::const_iterator ipar;
+ while( (ipar=particleMap.find(equiv_id)) == particleMap.end() ) {
+ TrackEquivalents::const_iterator iequiv = equivalentTracks.find(equiv_id);
+ equiv_id = (iequiv == equivalentTracks.end()) ? -1 : (*iequiv).second;
+ if ( equiv_id < 0 ) {
+ printout(INFO,"Geant4ParticleMap",
+ "+++ No Equivalent particle for track:%d last known is:%d",
+ g4_id,equiv_id);
+ break;
+ }
+ }
+ }
+ return equiv_id;
+}
diff --git a/DDG4/src/Geant4ParticleGun.cpp b/DDG4/src/Geant4ParticleGun.cpp
index bfefea6fa..86af0110c 100644
--- a/DDG4/src/Geant4ParticleGun.cpp
+++ b/DDG4/src/Geant4ParticleGun.cpp
@@ -10,6 +10,8 @@
// Framework include files
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Random.h"
#include "DDG4/Geant4ParticleGun.h"
#include "CLHEP/Units/SystemOfUnits.h"
@@ -17,8 +19,6 @@
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
-#include "TRandom1.h"
-
// C/C++ include files
#include <stdexcept>
#include <limits>
@@ -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_rndm(0), m_shotNo(0)
+ m_particle(0), m_gun(0), m_shotNo(0)
{
InstanceCount::increment(this);
m_needsControl = true;
@@ -46,8 +46,6 @@ Geant4ParticleGun::Geant4ParticleGun(Geant4Context* context, const string& name)
Geant4ParticleGun::~Geant4ParticleGun() {
if (m_gun)
delete m_gun;
- if ( m_rndm )
- delete m_rndm;
InstanceCount::decrement(this);
}
@@ -60,13 +58,13 @@ void Geant4ParticleGun::operator()(G4Event* event) {
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
m_particle = particleTable->FindParticle(m_particleName);
if (0 == m_particle) {
- throw runtime_error("Bad particle type:"+m_particleName+"!");
+ throw runtime_error("Geant4ParticleGun: Bad particle type:"+m_particleName+"!");
}
}
if ( m_isotrop ) {
- if ( 0 == m_rndm ) m_rndm = new TRandom1();
- double phi = 2*M_PI*m_rndm->Rndm();
- double theta = M_PI*m_rndm->Rndm();
+ Geant4Random& rnd = context()->event().random();
+ double phi = 2*M_PI*rnd.rndm();
+ double theta = M_PI*rnd.rndm();
double x1 = std::sin(theta)*std::cos(phi);
double x2 = std::sin(theta)*std::sin(phi);
double x3 = std::cos(theta);
@@ -78,7 +76,7 @@ void Geant4ParticleGun::operator()(G4Event* event) {
m_direction.SetXYZ(m_direction.X()/r, m_direction.Y()/r, m_direction.Z()/r);
}
}
- print("Geant4ParticleGun","Shoot [%d] %.3f GeV %s pos:(%.3f %.3f %.3f)[mm] dir:(%6.3f %6.3f %6.3f)",
+ 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_direction.X(),m_direction.Y(), m_direction.Z());
m_gun->SetParticleDefinition(m_particle);
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index eb9e70c63..c844cfac6 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -7,7 +7,7 @@
//
//====================================================================
// Framework include files
-#include "DD4hep/Printout.h"
+//#include "DD4hep/Printout.h"
#include "DD4hep/Primitives.h"
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4StepHandler.h"
@@ -30,6 +30,7 @@
#include <set>
#include <stdexcept>
+#include <algorithm>
using namespace std;
using namespace DD4hep;
@@ -53,8 +54,8 @@ namespace {
}
/// Standard constructor
-Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* context, const std::string& nam)
- : Geant4GeneratorAction(context,nam), Geant4MonteCarloTruth(), m_userHandler(0)
+Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* context, const string& nam)
+ : Geant4GeneratorAction(context,nam), Geant4MonteCarloTruth(), m_userHandler(0), m_primaryMap(0)
{
//generatorAction().adopt(this);
eventAction().callAtBegin(this,&Geant4ParticleHandler::beginEvent);
@@ -62,12 +63,12 @@ Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* context, const std::
trackingAction().callAtFinal(this,&Geant4ParticleHandler::end,CallbackSequence::FRONT);
trackingAction().callUpFront(this,&Geant4ParticleHandler::begin,CallbackSequence::FRONT);
steppingAction().call(this,&Geant4ParticleHandler::step);
-
- declareProperty("printEndTracking", m_printEndTracking = false);
- declareProperty("printStartTracking", m_printStartTracking = false);
- declareProperty("keepAllParticles", m_keepAll = false);
- declareProperty("saveProcesses", m_processNames);
- declareProperty("minimalKineticEnergy",m_kinEnergyCut = 100e0*MeV);
+ m_globalParticleID = 0;
+ declareProperty("PrintEndTracking", m_printEndTracking = false);
+ declareProperty("PrintStartTracking", m_printStartTracking = false);
+ declareProperty("KeepAllParticles", m_keepAll = false);
+ declareProperty("SaveProcesses", m_processNames);
+ declareProperty("MinimalKineticEnergy",m_kinEnergyCut = 100e0*MeV);
InstanceCount::increment(this);
}
@@ -98,8 +99,7 @@ bool Geant4ParticleHandler::adopt(Geant4UserParticleHandler* action) {
/// Clear particle maps
void Geant4ParticleHandler::clear() {
- for(ParticleMap::iterator i=m_particleMap.begin(); i!=m_particleMap.end(); ++i)
- delete (*i).second;
+ releaseObjects(m_particleMap)();
m_particleMap.clear();
m_equivalentTracks.clear();
}
@@ -144,7 +144,6 @@ void Geant4ParticleHandler::mark(const G4Track* track) {
mask.set(G4PARTICLE_CREATED_CALORIMETER_HIT);
}
else if ( !mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT) ) {
- //printout(INFO,name(),"+++ Create TRACKER HIT: id=%d",m_currTrack.id);
mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
}
}
@@ -152,7 +151,7 @@ void Geant4ParticleHandler::mark(const G4Track* track) {
/// Event generation action callback
void Geant4ParticleHandler::operator()(G4Event* event) {
typedef Geant4MonteCarloTruth _MC;
- printout(INFO,name(),"+++ Add EVENT extension of type Geant4ParticleHandler.....");
+ info("+++ Event:%d Add EVENT extension of type Geant4ParticleHandler.....",event->GetEventID());
context()->event().addExtension((_MC*)this, typeid(_MC), 0);
clear();
/// Call the user particle handler
@@ -163,14 +162,14 @@ void Geant4ParticleHandler::operator()(G4Event* event) {
/// User stepping callback
void Geant4ParticleHandler::step(const G4Step* step, G4SteppingManager* mgr) {
- typedef std::vector<const G4Track*> _Sec;
+ typedef vector<const G4Track*> _Sec;
Geant4StepHandler h(step);
++m_currTrack.steps;
const G4ThreeVector& v = h.postPosG4();
m_currTrack.vex = v.x();
m_currTrack.vey = v.y();
m_currTrack.vez = v.z();
- if ( m_currTrack.energy > m_kinEnergyCut ) {
+ if ( h.trkKineEnergy() > m_kinEnergyCut ) {
//
// Tracks below the energy threshold are NOT stored.
// If these tracks produce hits or are selected due to another signature,
@@ -193,18 +192,63 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
double kine = h.kineticEnergy();
G4ThreeVector m = track->GetMomentum();
const G4ThreeVector& v = h.vertex();
+ int reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0;
+ G4PrimaryParticle* prim = primary(h.id(),context()->event().event());
+ Particle* prim_part = 0;
+
+
+ if ( prim ) {
+ Geant4PrimaryMap::Primaries::const_iterator iprim = m_primaryMap->primaryMap.find(prim);
+ if ( iprim == m_primaryMap->primaryMap.end() ) {
+ except("+++ Tracking preaction: Primary particle without generator particle!");
+ }
+ prim_part = (*iprim).second;
+ reason |= (G4PARTICLE_PRIMARY|G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
+ m_particleMap[h.id()] = prim_part->addRef();
+ }
- // Extract here all the information from the start of the tracking action
- // which we will need later to create a proper MC particle.
- m_currTrack.id = h.id();
- m_currTrack.reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0;
+ if ( prim_part ) {
+ m_currTrack.id = prim_part->id;
+ m_currTrack.reason = prim_part->reason|reason;
+ m_currTrack.usermask = prim_part->usermask;
+ m_currTrack.status = prim_part->status;
+ m_currTrack.spin[0] = prim_part->spin[0];
+ m_currTrack.spin[1] = prim_part->spin[1];
+ m_currTrack.spin[2] = prim_part->spin[2];
+ m_currTrack.colorFlow[0] = prim_part->colorFlow[0];
+ m_currTrack.colorFlow[1] = prim_part->colorFlow[1];
+ m_currTrack.parents = prim_part->parents;
+ m_currTrack.daughters = prim_part->daughters;
+ m_currTrack.definition = prim_part->definition;
+ m_currTrack.pdgID = prim_part->pdgID;
+ m_currTrack.mass = prim_part->mass;
+ }
+ else {
+ m_currTrack.id = m_globalParticleID;
+ m_currTrack.reason = reason;
+ m_currTrack.usermask = 0;
+ m_currTrack.status |= G4PARTICLE_SIM_CREATED;
+ m_currTrack.spin[0] = 0;
+ m_currTrack.spin[1] = 0;
+ m_currTrack.spin[2] = 0;
+ m_currTrack.colorFlow[0] = 0;
+ m_currTrack.colorFlow[1] = 0;
+ m_currTrack.parents.clear();
+ m_currTrack.daughters.clear();
+ m_currTrack.definition = h.trackDef();
+ m_currTrack.pdgID = h.trackDef()->GetPDGEncoding();
+ m_currTrack.mass = h.trackDef()->GetPDGMass();
+ // Once the daughter gets simulated, the parent is already done....
+ ParticleMap::iterator ipar = m_particleMap.find(h.parent());
+ if ( ipar != m_particleMap.end() ) {
+ Particle* p_par = (*ipar).second;
+ m_currTrack.parents.insert(p_par->id);
+ }
+ ++m_globalParticleID;
+ }
m_currTrack.steps = 0;
m_currTrack.secondaries = 0;
- m_currTrack.pdgID = h.trackDef()->GetPDGEncoding();
- m_currTrack.energy = kine;
m_currTrack.g4Parent = h.parent();
- m_currTrack.parent = h.parent();
- m_currTrack.definition = h.trackDef();
m_currTrack.process = h.creatorProcess();
m_currTrack.time = h.globalTime();
m_currTrack.vsx = v.x();
@@ -219,7 +263,6 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
m_currTrack.pex = 0.0;
m_currTrack.pey = 0.0;
m_currTrack.pez = 0.0;
- m_currTrack.daughters.clear();
// If the creator process of the track is in the list of process products to be kept, set the proper flag
if ( m_currTrack.process ) {
Processes::iterator i=find(m_processNames.begin(),m_processNames.end(),m_currTrack.process->GetProcessName());
@@ -230,6 +273,7 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
if ( m_keepAll ) {
PropertyMask(m_currTrack.reason).set(G4PARTICLE_KEEP_ALWAYS);
}
+
/// Initial update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->begin(track, m_currTrack);
@@ -239,18 +283,8 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
/// Post-track action callback
void Geant4ParticleHandler::end(const G4Track* track) {
Geant4TrackHandler h(track);
- int id = h.id();
+ int g4_id = h.id();
PropertyMask mask(m_currTrack.reason);
- G4PrimaryParticle* prim = primary(id,context()->event().event());
-
- if ( prim ) {
- m_currTrack.reason |= (G4PARTICLE_PRIMARY|G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
- const G4ParticleDefinition* def = m_currTrack.definition;
- printout(INFO,name(),"+++ Add Primary particle %d def:%p [%s, %s] reason:%d",
- id, def,
- def ? def->GetParticleName().c_str() : "UNKNOWN",
- def ? def->GetParticleType().c_str() : "UNKNOWN", m_currTrack.reason);
- }
if ( !mask.isNull() ) {
// These are candate tracks with a probability to be stored due to their properties:
@@ -261,85 +295,162 @@ void Geant4ParticleHandler::end(const G4Track* track) {
// - to be kept due to creator process
//
// Update vertex end point and final momentum
+ Geant4ParticleHandle ph(&m_currTrack);
G4ThreeVector m = track->GetMomentum();
- m_currTrack.pex = m.x();
- m_currTrack.pey = m.y();
- m_currTrack.pez = m.z();
-
+ ph->pex = m.x();
+ ph->pey = m.y();
+ ph->pez = m.z();
+ m_equivalentTracks[g4_id] = g4_id;
+ // Add this track to it's parents list of daughters
+ ParticleMap::iterator ipar = m_particleMap.find(ph->g4Parent);
+ if ( ipar != m_particleMap.end() ) {
+ Particle* p_par = (*ipar).second;
+ p_par->daughters.insert(ph->id);
+ ph->usermask |= p_par->usermask;
+ 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....
+ }
+ }
/// Final update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->begin(track, m_currTrack);
}
-
+ ParticleMap::iterator ip = m_particleMap.find(g4_id);
+ if ( mask.isSet(G4PARTICLE_PRIMARY) ) {
+ ph.dump2(outputLevel()-1,name(),"Add Primary",h.id(),ip!=m_particleMap.end());
+ }
// Create a new MC particle from the current track information saved in the pre-tracking action
- Particle* p = m_particleMap[id] = new Particle();
+ Particle* p = ip==m_particleMap.end() ? (m_particleMap[g4_id] = new Particle()) : (*ip).second;
p->get_data(m_currTrack);
- // Add this track to it's parents list of daughters
- ParticleMap::iterator ipar = m_particleMap.find(p->g4Parent);
- if ( ipar != m_particleMap.end() ) {
- Particle* p_par = (*ipar).second;
- p_par->daughters.insert(id);
- }
- m_equivalentTracks[id] = id;
-#if 0
- /// Some printout
- const char* hit = mask.isSet(G4PARTICLE_CREATED_HIT) ? "CREATED_HIT" : "";
- const char* sec = mask.isSet(G4PARTICLE_HAS_SECONDARIES) ? "SECONDARIES" : "";
- const char* prim = mask.isSet(G4PARTICLE_PRIMARY) ? "PRIMARY" : "";
- printout(INFO,name(),"+++ Tracking postaction: %d/%d par:%d [%s, %s] created by:%s E:%e state:%s %s %s %s",
- id,int(m_particleMap.size()),
- h.parent(),h.name().c_str(),h.type().c_str(),h.creatorName().c_str(),
- p->energy, h.statusName(), prim, hit, sec);
-#endif
}
else {
// These are tracks without any special properties.
//
// We will not store them on the record, but have to memorise the
// track identifier in order to restore the history for the created hits.
- m_equivalentTracks[id] = h.parent();
+ m_equivalentTracks[g4_id] = m_currTrack.g4Parent;
}
}
/// Pre-event action callback
void Geant4ParticleHandler::beginEvent(const G4Event* event) {
+ Geant4PrimaryInteraction* interaction = context()->event().extension<Geant4PrimaryInteraction>();
+ info("+++ Event %d Begin event action. Access event related information.",event->GetEventID());
+ m_primaryMap = context()->event().extension<Geant4PrimaryMap>();
+ m_globalParticleID = interaction->nextPID();
+ m_particleMap.clear();
+ m_equivalentTracks.clear();
/// Call the user particle handler
if ( m_userHandler ) {
m_userHandler->begin(event);
}
}
+/// Debugging: Dump Geant4 particle map
+void Geant4ParticleHandler::dumpMap(const char* tag) const {
+ for(ParticleMap::const_iterator iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i)
+ Geant4ParticleHandle((*i).second).dump1(INFO,name(),tag);
+}
+
/// Post-event action callback
void Geant4ParticleHandler::endEvent(const G4Event* event) {
int count = 0;
+ int level = outputLevel();
do {
- printout(INFO,name(),"+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
+ if ( level <= VERBOSE ) dumpMap("Particle");
+ print("+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
} while( recombineParents() > 0 );
- // Update the particle's parents and replace the original Geant4 track with the
- // equivalent particle still present in the record.
- for(ParticleMap::const_iterator ipar, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
- Particle* p = (*i).second;
- p->parent = equivalentTrack(p->g4Parent);
- if ( p->parent > 0 ) {
- ipar = m_particleMap.find(p->parent);
- if ( ipar != iend ) {
- set<int>& daughters = (*ipar).second->daughters;
- if ( daughters.find(p->id) == daughters.end() ) daughters.insert(p->id);
- }
- }
- }
- /// No we have to updte the map of equivalent tracks and assign the 'equivalentTrack' entry
- for(TrackEquivalents::iterator i=m_equivalentTracks.begin(); i!=m_equivalentTracks.end(); ++i) {
- int& track_id = (*i).second;
- int equiv_id = equivalentTrack(track_id);
- track_id = equiv_id;
- }
+
+ if ( level <= VERBOSE ) dumpMap("Recombined");
+ // Rebase the simulated tracks, so that they fit to the generator particles
+ rebaseSimulatedTracks(0);
+ if ( level <= DEBUG ) dumpMap("Rebased");
// Consistency check....
checkConsistency();
/// Call the user particle handler
if ( m_userHandler ) {
m_userHandler->end(event);
}
+ // Now export the data to the final record.
+ Geant4ParticleMap* part_map = context()->event().extension<Geant4ParticleMap>();
+ part_map->adopt(m_particleMap, m_equivalentTracks);
+ m_primaryMap = 0;
+ clear();
+}
+
+/// Rebase the simulated tracks, so that they fit to the generator particles
+void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
+ /// No we have to update the map of equivalent tracks and assign the 'equivalentTrack' entry
+ TrackEquivalents equivalents, orgParticles;
+ ParticleMap finalParticles;
+ ParticleMap::const_iterator ipar, iend, i;
+ int count;
+
+ Geant4PrimaryInteraction* interaction = context()->event().extension<Geant4PrimaryInteraction>();
+ ParticleMap& pm = interaction->particles;
+
+ // (1.0) Copy the pre-defined particle mapping for the simulated tracks
+ // It is assumed the mapping is ZERO based without holes.
+ for(count = 0, iend=pm.end(), i=pm.begin(); i!=iend; ++i) {
+ Particle* p = (*i).second;
+ orgParticles[p->id] = p->id;
+ finalParticles[p->id] = p;
+ if ( p->id > count ) count = p->id;
+ if ( (p->reason&G4PARTICLE_PRIMARY) != G4PARTICLE_PRIMARY ) {
+ p->addRef();
+ }
+ }
+ // (1.1) Define the new particle mapping for the simulated tracks
+ for(++count, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
+ Particle* p = (*i).second;
+ if ( (p->reason&G4PARTICLE_PRIMARY) != G4PARTICLE_PRIMARY ) {
+ //if ( orgParticles.find(p->id) == orgParticles.end() ) {
+ orgParticles[p->id] = count;
+ finalParticles[count] = p;
+ p->id = count;
+ ++count;
+ }
+ }
+ // (2) Re-evaluate the corresponding geant4 track equivalents using the new mapping
+ for(TrackEquivalents::iterator i=m_equivalentTracks.begin(),iend=m_equivalentTracks.end(); i!=iend; ++i) {
+ int g4_equiv = (*i).first;
+ ParticleMap::const_iterator ipar;
+ while( (ipar=m_particleMap.find(g4_equiv)) == m_particleMap.end() ) {
+ TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(g4_equiv);
+ if ( iequiv == iend )
+ break; // ERROR !! Will be handled by printout below because ipar==end()
+ g4_equiv = (*iequiv).second;
+ }
+ if ( ipar != m_particleMap.end() )
+ equivalents[(*i).first] = (*ipar).second->id; // requires (1) !
+ else
+ error("+++ No Equivalent particle for track:%d last known is:%d",(*i).second,g4_equiv);
+ }
+ m_equivalentTracks = equivalents;
+ equivalents.clear();
+
+ // (3) Update the particle's parents and replace the original Geant4 track with the
+ // equivalent particle still present in the record.
+ for(ParticleMap::const_iterator ipar, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
+ Particle* p = (*i).second;
+ set<int> daughters = p->daughters;
+ 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 ( p->g4Parent > 0 ) {
+ p->parents.clear();
+ p->parents.insert(equivalentTrack(p->g4Parent)); // Requires (2)
+ ipar = m_particleMap.find(p->g4Parent);
+ if ( ipar != iend ) {
+ set<int>& dau = (*ipar).second->daughters;
+ if ( dau.find(p->id) == dau.end() ) dau.insert(p->id);
+ }
+ }
+ }
+ }
+ m_particleMap = finalParticles;
}
/// Clean the monte carlo record. Remove all unwanted stuff.
@@ -350,11 +461,12 @@ int Geant4ParticleHandler::recombineParents() {
/// Need to start from BACK, to clean first the latest produced stuff.
/// Otherwise the daughter list of the earlier produced tracks would not be empty!
for(ParticleMap::reverse_iterator i=m_particleMap.rbegin(); i!=m_particleMap.rend(); ++i) {
- Particle* par = (*i).second;
- set<int>& daughters = par->daughters;
+ int g4_id = (*i).first;
+ Particle* p = (*i).second;
+ set<int>& daughters = p->daughters;
if ( daughters.size() == 0 ) {
- PropertyMask mask(par->reason);
- int id = par->id;
+ PropertyMask mask(p->reason);
+ int id = p->id;
bool secondaries = mask.isSet(G4PARTICLE_HAS_SECONDARIES);
bool tracker_track = mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT);
bool calo_track = mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT);
@@ -362,7 +474,6 @@ int Geant4ParticleHandler::recombineParents() {
bool low_energy = !mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD);
bool keep_process = mask.isSet(G4PARTICLE_KEEP_PROCESS);
bool keep_parent = mask.isSet(G4PARTICLE_KEEP_PARENT);
- int parent_id = par->g4Parent;
bool remove_me = false;
if ( id == break_trackID ) { // Used for debugging to set break point
@@ -377,7 +488,7 @@ int Geant4ParticleHandler::recombineParents() {
//continue;
}
else if ( keep_process ) {
- ParticleMap::iterator ip = m_particleMap.find(parent_id);
+ ParticleMap::iterator ip = m_particleMap.find(p->g4Parent);
if ( ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
PropertyMask parent_mask(parent_part->reason);
@@ -408,19 +519,19 @@ int Geant4ParticleHandler::recombineParents() {
/// Remove this track from the list and also do the cleanup in the parent's children list
if ( remove_me ) {
- ParticleMap::iterator ip = m_particleMap.find(parent_id);
- remove.insert(id);
- m_equivalentTracks[id] = parent_id;
+ ParticleMap::iterator ip = m_particleMap.find(p->g4Parent);
+ remove.insert(g4_id);
+ m_equivalentTracks[g4_id] = p->g4Parent;
if ( ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
PropertyMask(parent_part->reason).set(mask.value());
// Remove track from parent's list of daughters
parent_part->removeDaughter(id);
- parent_part->steps += par->steps;
- parent_part->secondaries += par->secondaries;
+ parent_part->steps += p->steps;
+ parent_part->secondaries += p->secondaries;
/// Update of the particle using the user handler
if ( m_userHandler ) {
- m_userHandler->combine(*par, *parent_part);
+ m_userHandler->combine(*p, *parent_part);
}
}
}
@@ -429,7 +540,7 @@ int Geant4ParticleHandler::recombineParents() {
for(set<int>::const_iterator r=remove.begin(); r!=remove.end();++r) {
ParticleMap::iterator ir = m_particleMap.find(*r);
if ( ir != m_particleMap.end() ) {
- delete (*ir).second;
+ (*ir).second->release();
m_particleMap.erase(ir);
}
}
@@ -444,42 +555,46 @@ void Geant4ParticleHandler::checkConsistency() const {
for(ParticleMap::const_iterator j, i=m_particleMap.begin(); i!=m_particleMap.end(); ++i) {
Particle* p = (*i).second;
PropertyMask mask(p->reason);
+ PropertyMask status(p->status);
set<int>& daughters = p->daughters;
// For all particles, the set of daughters must be contained in the record.
for(set<int>::const_iterator id=daughters.begin(); id!=daughters.end(); ++id) {
int id_dau = *id;
- j = m_particleMap.find(id_dau);
- if ( j == m_particleMap.end() ) {
+ if ( (j=m_particleMap.find(id_dau)) == m_particleMap.end() ) {
++num_errors;
- printout(INFO,name(),"+++ Particle:%d Daughter %d is not in particle map!",p->id,id_dau);
+ error("+++ Particle:%d Daughter %d is not in particle map!",p->id,id_dau);
}
}
- // For all particles except the primaries, the parent must be contained in the record.
- if ( !mask.isSet(G4PARTICLE_PRIMARY) ) {
- int id_par = p->parent;
- j = m_particleMap.find(id_par);
- if ( j == m_particleMap.end() ) {
+ // 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) ) {
+ 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();
+ char parent_list[1024];
+ parent_list[0] = 0;
+ if ( !in_map || !in_parent_list ) {
++num_errors;
- printout(INFO,name(),"+++ Particle:%d Parent %d is not in particle map!",p->id,id_par);
+ for(set<int>::const_iterator ip=p->parents.begin(); ip!=p->parents.end();++ip)
+ ::snprintf(parent_list+strlen(parent_list),sizeof(parent_list)-strlen(parent_list),"%d ",*ip);
+ error("+++ Particle:%d Parent %d (G4id:%d) In record:%s In parent list:%s [%s]",
+ p->id,parent_id,p->g4Parent,yes_no(in_map),yes_no(in_parent_list),parent_list);
}
}
}
/// No we have to check the consistency of the map of equivalent tracks used to assign the
/// proper MC particle to the created hits
- for(TrackEquivalents::const_iterator i=m_equivalentTracks.begin(); i!=m_equivalentTracks.end(); ++i) {
- int track_id = (*i).second;
- int equiv_id = equivalentTrack(track_id);
- ParticleMap::const_iterator j = m_particleMap.find(equiv_id);
- if ( j == m_particleMap.end() ) {
- int g4_id = (*i).first;
+ for(TrackEquivalents::const_iterator i=m_equivalentTracks.begin(), iend=m_equivalentTracks.end(); i!=iend; ++i) {
+ int g4_id = (*i).first;
+ int equiv_id = equivalentTrack(g4_id);
+ if ( equiv_id < 0 ) {
++num_errors;
- printout(INFO,name(),"+++ Geant4 Track:%d Track:%d Equivalent track %d is not in particle map!",g4_id,track_id,equiv_id);
}
}
if ( num_errors > 0 ) {
- printout(ERROR,name(),"+++ Consistency check failed. Found %d problems.",num_errors);
+ except("+++ Consistency check failed. Found %d problems.",num_errors);
}
}
@@ -487,17 +602,10 @@ void Geant4ParticleHandler::checkConsistency() const {
int Geant4ParticleHandler::equivalentTrack(int g4_id) const {
int equiv_id = g4_id;
if ( g4_id != 0 ) {
- ParticleMap::const_iterator ipar;
- while( (ipar=m_particleMap.find(equiv_id)) == m_particleMap.end() ) {
- TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(equiv_id);
- equiv_id = (iequiv == m_equivalentTracks.end()) ? -1 : (*iequiv).second;
- if ( equiv_id < 0 ) {
- printout(INFO,name(),"+++ No Equivalent particle for track:%d last known is:%d",g4_id,equiv_id);
- break;
- }
- }
+ TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(equiv_id);
+ return (iequiv == m_equivalentTracks.end()) ? -1 : (*iequiv).second;
}
- return equiv_id;
+ return -1;
}
/// Access the equivalent track id (shortcut to the usage of TrackEquivalents)
diff --git a/DDG4/src/Geant4ParticlePrint.cpp b/DDG4/src/Geant4ParticlePrint.cpp
index 8a7768f07..ad6aa9395 100644
--- a/DDG4/src/Geant4ParticlePrint.cpp
+++ b/DDG4/src/Geant4ParticlePrint.cpp
@@ -13,10 +13,6 @@
#include "DDG4/Geant4ParticlePrint.h"
#include "DDG4/Geant4Data.h"
-#include "G4ParticleDefinition.hh"
-#include "G4VProcess.hh"
-#include "G4Event.hh"
-
using namespace std;
using namespace DD4hep;
using namespace DD4hep::Simulation;
@@ -28,6 +24,9 @@ Geant4ParticlePrint::Geant4ParticlePrint(Geant4Context* context, const std::stri
: Geant4EventAction(context,nam)
{
declareProperty("OutputType",m_outputType=3);
+ declareProperty("PrintBegin",m_printBegin=false);
+ declareProperty("PrintEnd", m_printEnd=true);
+ declareProperty("PrintGen", m_printGeneration=true);
InstanceCount::increment(this);
}
@@ -36,53 +35,63 @@ Geant4ParticlePrint::~Geant4ParticlePrint() {
InstanceCount::decrement(this);
}
+/// Print particle table
+void Geant4ParticlePrint::makePrintout() const {
+ Geant4ParticleMap* parts = context()->event().extension<Geant4ParticleMap>();
+ if ( parts ) {
+ const ParticleMap& particles = parts->particles();
+ if ( (m_outputType&2) != 0 ) printParticleTree(particles); // Tree printout....
+ if ( (m_outputType&1) != 0 ) printParticles(particles); // Table printout....
+ return;
+ }
+ except("No Geant4MonteCarloTruth object present in the event structure to access the particle information!", c_name());
+}
+
+/// Generation action callback
+void Geant4ParticlePrint::operator()(G4Event* ) {
+ if ( m_printGeneration ) makePrintout();
+}
+
/// Pre-event action callback
void Geant4ParticlePrint::begin(const G4Event* ) {
+ if ( m_printBegin ) makePrintout();
}
/// Post-event action callback
void Geant4ParticlePrint::end(const G4Event* ) {
- Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>();
- if ( truth ) {
- const ParticleMap& particles = truth->particles();
- if ( (m_outputType&1) != 0 ) printParticles(particles); // Table printout....
- if ( (m_outputType&2) != 0 ) printParticleTree(particles); // Tree printout....
- return;
- }
- except("No Geant4MonteCarloTruth object present in the event structure to access the particle information!", c_name());
+ if ( m_printEnd ) makePrintout();
}
-void Geant4ParticlePrint::printParticle(const std::string& prefix, const Particle* p) const {
+void Geant4ParticlePrint::printParticle(const std::string& prefix, Geant4ParticleHandle p) const {
char equiv[32];
PropertyMask mask(p->reason);
- const char* proc_name = "???";
- const char* proc_type = p->process ? G4VProcess::GetProcessTypeName(p->process->GetProcessType()).c_str() : "";
-
- if ( p->process ) proc_name = p->process->GetProcessName().c_str();
- else if ( mask.isSet(G4PARTICLE_PRIMARY) ) proc_name = "Primary";
+ string proc_name = p.processName();
+ string proc_type = p.processTypeName();
+ int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
equiv[0] = 0;
- if ( p->g4Parent != p->parent ) {
+ if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
}
- printout(INFO,name(),"+++ %sID: %6d %12s %6d%-7s %7s %3s %5d %6s %8.3g %-4s %-7s %-7s %-3s [%s%s%s]",
- prefix.c_str(),
- p->id,
- p->definition->GetParticleName().c_str(),
- p->parent,equiv,
- yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
- yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
- int(p->daughters.size()),
- yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
- p->energy,
- yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
- yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
- yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
- mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
- proc_name,
- p->process ? "/" : "",
- proc_type
- );
+ print("+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s]",
+ prefix.c_str(),
+ p->id,
+ p.particleName().c_str(),
+ parent_id,equiv,
+ yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
+ yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
+ int(p->daughters.size()),
+ yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
+ p.energy(),
+ yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
+ mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
+ p.numParent(),
+ proc_name.c_str(),
+ p->process ? "/" : "",
+ proc_type.c_str()
+ );
}
/// Print record of kept particles
@@ -95,11 +104,11 @@ void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const {
int num_calo_hits = 0;
int num_tracker_hits = 0;
- printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
- "Primary Secondary Energy %-8s Calo Tracker Process/Par Details",
- int(particles.size()),"ParticleType","","in [MeV]");
+ print("+++ MC Particles #Tracks:%7d ParticleType Parent/Geant4 "
+ "Primary Secondary Energy in [MeV] Calo Tracker Process/Par Details",
+ int(particles.size()));
for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
- const Particle* p = (*i).second;
+ Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->reason);
printParticle("MC Particle Track",p);
num_secondaries += int(p->daughters.size());
@@ -110,17 +119,15 @@ void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const {
if ( mask.isSet(G4PARTICLE_KEEP_PARENT) ) ++num_parent;
else if ( mask.isSet(G4PARTICLE_KEEP_PROCESS) ) ++num_process;
}
- printout(INFO,name(),"+++ MC Particles #Tracks:%6d %-12s Parent%-7s "
- "Primary Secondary Energy %-8s Calo Tracker Process/Parent",
- int(particles.size()),"","","");
- printout(INFO,name(),"+++ MC Particle Summary: %6s %-12s %6s%-7s "
- "%7d %9d %6d %-8s %4d %7d %7d %6d",
- "","","","",
- num_primary, num_secondaries, num_energy,"",
- num_calo_hits,num_tracker_hits,num_process,num_parent);
+ print("+++ MC Particles #Tracks:%7d ParticleType Parent/Geant4 "
+ "Primary Secondary Energy Calo Tracker Process/Par",
+ int(particles.size()));
+ print("+++ MC Particle Summary: %7d %10d %7d %7d %9d %5d %6d",
+ num_primary, num_secondaries, num_energy,
+ num_calo_hits,num_tracker_hits,num_process,num_parent);
}
-void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int level, const Particle* p) const {
+void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int level, Geant4ParticleHandle p) const {
char txt[32];
size_t len = sizeof(txt)-1;
// Ensure we do not overwrite the array
@@ -128,7 +135,7 @@ void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int le
::snprintf(txt,sizeof(txt),"%5d ",level);
::memset(txt+6,' ',len-6);
- txt[len] = 0;
+ txt[len-1] = 0;
txt[len-2] = '>';
txt[level+6]='+';
::memset(txt+level+6+1,'-',len-level-3-6);
@@ -138,19 +145,19 @@ void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int le
// For all particles, the set of daughters must be contained in the record.
for(set<int>::const_iterator id=daughters.begin(); id!=daughters.end(); ++id) {
int id_dau = *id;
- const Particle* dau = (*particles.find(id_dau)).second;
+ Geant4ParticleHandle dau = (*particles.find(id_dau)).second;
printParticleTree(particles,level+1,dau);
}
}
/// Print tree of kept particles
void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles) const {
- printout(INFO,name(),"+++ MC Particle Parent daughter relationships. [%d particles]",int(particles.size()));
- printout(INFO,name(),"+++ MC Particles %12s #Tracks:%7d %-12s Parent%-7s "
- "Primary Secondary Energy %-8s Calo Tracker Process/Par Details",
- "",int(particles.size()),"ParticleType","","in [MeV]");
+ print("+++ MC Particle Parent daughter relationships. [%d particles]",int(particles.size()));
+ print("+++ MC Particles %12s #Tracks:%7d %-12s Parent%-7s "
+ "Primary Secondary Energy %-8s Calo Tracker Process/Par Details",
+ "",int(particles.size()),"ParticleType","","in [MeV]");
for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
- const Particle* p = (*i).second;
+ Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->reason);
if ( mask.isSet(G4PARTICLE_PRIMARY) ) printParticleTree(particles,0,p);
}
diff --git a/DDG4/src/Geant4Primary.cpp b/DDG4/src/Geant4Primary.cpp
new file mode 100644
index 000000000..bfd24d1b1
--- /dev/null
+++ b/DDG4/src/Geant4Primary.cpp
@@ -0,0 +1,119 @@
+// $Id: Geant4Hits.cpp 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+
+// Framework include files
+#include "DD4hep/Primitives.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Primary.h"
+
+// C/C++ include files
+#include <stdexcept>
+#include <cstdio>
+
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+/// Default destructor
+PrimaryExtension::~PrimaryExtension() {
+}
+
+/// Default destructor
+Geant4PrimaryMap::~Geant4PrimaryMap() {
+ releaseObjects(primaryMap)();
+}
+
+/// Default constructor
+Geant4PrimaryInteraction::Geant4PrimaryInteraction()
+ : mask(0), next_particle_identifier(-1)
+{
+}
+
+/// Copy constructor
+Geant4PrimaryInteraction::Geant4PrimaryInteraction(const Geant4PrimaryInteraction&)
+ : mask(0), next_particle_identifier(-1)
+{
+}
+
+/// Assignment operator
+Geant4PrimaryInteraction& Geant4PrimaryInteraction::operator=(const Geant4PrimaryInteraction& c) {
+ if ( &c == this ) {}
+ return *this;
+}
+
+/// Default destructor
+Geant4PrimaryInteraction::~Geant4PrimaryInteraction() {
+ releaseObjects(vertices)();
+ releaseObjects(particles)();
+}
+
+/// Access a new particle identifier within the interaction
+int Geant4PrimaryInteraction::nextPID() {
+ return ++next_particle_identifier;
+}
+
+/// Access a new particle identifier within the interaction
+void Geant4PrimaryInteraction::setNextPID(int new_value) {
+ next_particle_identifier = new_value-1;
+}
+
+/// Default constructor
+Geant4PrimaryEvent::Geant4PrimaryEvent()
+{
+}
+
+/// Copy constructor
+Geant4PrimaryEvent::Geant4PrimaryEvent(const Geant4PrimaryEvent&)
+{
+}
+
+/// Assignment operator
+Geant4PrimaryEvent& Geant4PrimaryEvent::operator=(const Geant4PrimaryEvent& c) {
+ if ( &c == this ) {}
+ return *this;
+}
+
+/// Default destructor
+Geant4PrimaryEvent::~Geant4PrimaryEvent() {
+ destroyObjects(m_interactions)();
+}
+
+/// Add a new interaction object to the event
+void Geant4PrimaryEvent::add(int id, Geant4PrimaryInteraction* interaction) {
+ if ( interaction ) {
+ Interactions::iterator i = m_interactions.find(id);
+ if ( i == m_interactions.end() ) {
+ interaction->mask = id;
+ m_interactions.insert(std::make_pair(id,interaction));
+ return;
+ }
+ char text[132];
+ ::snprintf(text,sizeof(text),"Geant4PrimaryEvent: Interaction with ID '%d' "
+ "exists and cannot be added twice!",id);
+ throw std::runtime_error(text);
+ }
+ throw std::runtime_error("Geant4PrimaryEvent: CANNOT add invalid Interaction!");
+}
+
+/// Retrieve an interaction by it's ID
+Geant4PrimaryEvent::Interaction* Geant4PrimaryEvent::get(int mask) const {
+ Interactions::const_iterator i = m_interactions.find(mask);
+ if ( i != m_interactions.end() ) {
+ return (*i).second;
+ }
+ return 0;
+}
+
+/// Retrieve all intractions
+std::vector<Geant4PrimaryEvent::Interaction*> Geant4PrimaryEvent::interactions() const {
+ std::vector<Interaction*> v;
+ for(Interactions::const_iterator i=m_interactions.begin(); i!=m_interactions.end(); ++i)
+ v.push_back((*i).second);
+ return v;
+}
+
diff --git a/DDG4/src/Geant4PrimaryHandler.cpp b/DDG4/src/Geant4PrimaryHandler.cpp
new file mode 100644
index 000000000..f45ca9ed7
--- /dev/null
+++ b/DDG4/src/Geant4PrimaryHandler.cpp
@@ -0,0 +1,71 @@
+// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $
+//====================================================================
+// 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/Geant4PrimaryHandler.h"
+
+#include "G4Event.hh"
+#include "G4PrimaryVertex.hh"
+#include "G4PrimaryParticle.hh"
+
+#include <stdexcept>
+
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+/// Standard constructor
+Geant4PrimaryHandler::Geant4PrimaryHandler(Geant4Context* context, const std::string& nam)
+ : Geant4GeneratorAction(context,nam)
+{
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4PrimaryHandler::~Geant4PrimaryHandler() {
+ InstanceCount::decrement(this);
+}
+
+/// 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;
+
+ 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);
+ }
+ 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);
+ }
+}
diff --git a/DDG4/src/Geant4ROOTDump.cpp b/DDG4/src/Geant4ROOTDump.cpp
index 3d48b9358..457b0c7db 100644
--- a/DDG4/src/Geant4ROOTDump.cpp
+++ b/DDG4/src/Geant4ROOTDump.cpp
@@ -12,6 +12,7 @@
#include "DD4hep/Factories.h"
#include "DD4hep/Primitives.h"
#include "DDG4/Geant4DataDump.h"
+#include "DDG4/Geant4Particle.h"
#include "TInterpreter.h"
#include "TSystem.h"
@@ -127,7 +128,7 @@ static long dump_root(DD4hep::Geometry::LCDD&, int argc, char** argv) {
if ( data.first == cl_particles ) {
Particles* parts = (Particles*)data.second;
dump.print(INFO, (*i).first, parts);
- for_each(parts->begin(), parts->end(), DestroyObject<Particle*>());
+ for_each(parts->begin(), parts->end(), DestroyObject<Geant4Particle*>());
}
}
for(ENTRIES::const_iterator i=event.begin(); i!=event.end(); ++i) {
@@ -137,12 +138,12 @@ static long dump_root(DD4hep::Geometry::LCDD&, int argc, char** argv) {
else if ( data.first == cl_tracker ) {
TrackerHits* hits = (TrackerHits*)data.second;
dump.print(INFO, (*i).first, hits);
- for_each(hits->begin(), hits->end(), DestroyObject<SimpleTracker::Hit*>());
+ for_each(hits->begin(), hits->end(), DestroyObject<Geant4Tracker::Hit*>());
}
else if ( data.first == cl_calo ) {
CalorimeterHits* hits = (CalorimeterHits*)data.second;
dump.print(INFO, (*i).first, hits);
- for_each(hits->begin(), hits->end(), DestroyObject<SimpleCalorimeter::Hit*>());
+ for_each(hits->begin(), hits->end(), DestroyObject<Geant4Calorimeter::Hit*>());
}
if ( data.first ) data.first->Destructor(data.second);
}
diff --git a/DDG4/src/Geant4Random.cpp b/DDG4/src/Geant4Random.cpp
new file mode 100644
index 000000000..2e729ff7c
--- /dev/null
+++ b/DDG4/src/Geant4Random.cpp
@@ -0,0 +1,82 @@
+// $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/Geant4Random.h"
+
+#include "TRandom1.h"
+
+// C/C++ include files
+#include <cmath>
+
+using namespace std;
+using namespace DD4hep::Simulation;
+
+/// Default constructor
+Geant4Random::Geant4Random() {
+ if ( 0 == gRandom ) gRandom = new TRandom1();
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4Random::~Geant4Random() {
+ InstanceCount::decrement(this);
+}
+
+void Geant4Random::circle(double &x, double &y, double r) {
+ gRandom->Circle(x,y,r);
+}
+
+//
+double Geant4Random::exp(double tau) {
+ return gRandom->Exp(tau);
+}
+
+//
+double Geant4Random::gauss(double mean, double sigma) {
+ return gRandom->Gaus(mean,sigma);
+}
+
+//
+double Geant4Random::landau(double mean, double sigma) {
+ return gRandom->Landau(mean,sigma);
+}
+
+//
+double Geant4Random::rndm(int i) {
+ return gRandom->Rndm(i);
+}
+
+//
+void Geant4Random::rndmArray(int n, float *array) {
+ gRandom->RndmArray(n,array);
+}
+
+//
+void Geant4Random::rndmArray(int n, double *array) {
+ gRandom->RndmArray(n,array);
+}
+
+//
+void Geant4Random::sphere(double &x, double &y, double &z, double r) {
+ gRandom->Sphere(x,y,z,r);
+}
+
+//
+double Geant4Random::uniform(double x1) {
+ return gRandom->Uniform(x1);
+}
+
+//
+double Geant4Random::uniform(double x1, double x2) {
+ return gRandom->Uniform(x1,x2);
+}
+
diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp
index 7f190989d..287d7978a 100644
--- a/DDG4/src/Geant4SensDetAction.cpp
+++ b/DDG4/src/Geant4SensDetAction.cpp
@@ -217,8 +217,8 @@ Geant4SensDetActionSequence::Geant4SensDetActionSequence(Geant4Context* context,
m_needsControl = true;
context->sensitiveActions().insert(name(), this);
/// Update the sensitive detector type, so that the proper instance is created
- Geometry::SensitiveDetector sd = context->lcdd().sensitiveDetector(nam);
- sd.setType("Geant4SensDet");
+ m_sensitive = context->lcdd().sensitiveDetector(nam);
+ m_sensitive.setType("Geant4SensDet");
InstanceCount::increment(this);
}
@@ -334,7 +334,7 @@ void Geant4SensDetActionSequence::begin(G4HCofThisEvent* hce) {
m_actors(ContextUpdate(context()));
for (size_t count = 0; count < m_collections.size(); ++count) {
const std::pair<string, create_t>& cr = m_collections[count];
- G4VHitsCollection* c = (*cr.second)(name(), cr.first);
+ Geant4HitCollection* c = (*cr.second)(name(), cr.first, m_sensitive);
int id = m_detector->GetCollectionID(count);
m_hce->AddHitsCollection(id, c);
}
diff --git a/DDG4/src/Geant4TestActions.cpp b/DDG4/src/Geant4TestActions.cpp
index dae7a52e5..ef0eaf0b0 100644
--- a/DDG4/src/Geant4TestActions.cpp
+++ b/DDG4/src/Geant4TestActions.cpp
@@ -61,7 +61,7 @@ Geant4TestGeneratorAction::~Geant4TestGeneratorAction() {
/// Callback to generate primary particles
void Geant4TestGeneratorAction::operator()(G4Event* evt) {
- PRINT(name(), "%s> calling Geant4TestGeneratorAction(event_id=%d Context: run=%p evt=%p)",
+ PRINT("%s> calling Geant4TestGeneratorAction(event_id=%d Context: run=%p evt=%p)",
m_type.c_str(), evt->GetEventID(), &context()->run(), &context()->event());
}
@@ -78,25 +78,25 @@ Geant4TestRunAction::~Geant4TestRunAction() {
/// begin-of-run callback
void Geant4TestRunAction::begin(const G4Run* run) {
- PRINT(name(), "%s> calling begin(run_id=%d,num_event=%d Context:%p)", m_type.c_str(), run->GetRunID(),
+ PRINT("%s> calling begin(run_id=%d,num_event=%d Context:%p)", m_type.c_str(), run->GetRunID(),
run->GetNumberOfEventToBeProcessed(), &context()->run());
}
/// End-of-run callback
void Geant4TestRunAction::end(const G4Run* run) {
- PRINT(name(), "%s> calling end(run_id=%d, num_event=%d Context:%p)",
+ PRINT("%s> calling end(run_id=%d, num_event=%d Context:%p)",
m_type.c_str(), run->GetRunID(), run->GetNumberOfEvent(), &context()->run());
}
/// begin-of-event callback
void Geant4TestRunAction::beginEvent(const G4Event* evt) {
- PRINT(name(), "%s> calling beginEvent(event_id=%d Context: run=%p evt=%p)",
+ PRINT("%s> calling beginEvent(event_id=%d Context: run=%p evt=%p)",
m_type.c_str(), evt->GetEventID(), &context()->run(), &context()->event());
}
/// End-of-event callback
void Geant4TestRunAction::endEvent(const G4Event* evt) {
- PRINT(name(), "%s> calling endEvent(event_id=%d Context: run=%p evt=%p)",
+ PRINT("%s> calling endEvent(event_id=%d Context: run=%p evt=%p)",
m_type.c_str(), evt->GetEventID(), &context()->run(), &context()->event());
}
@@ -112,7 +112,7 @@ Geant4TestEventAction::~Geant4TestEventAction() {
}
/// begin-of-event callback
void Geant4TestEventAction::begin(const G4Event* evt) {
- PRINT(name(), "%s> calling begin(event_id=%d Context: run=%p (%d) evt=%p (%d))",
+ PRINT("%s> calling begin(event_id=%d Context: run=%p (%d) evt=%p (%d))",
m_type.c_str(), evt->GetEventID(),
&context()->run(), context()->run().run().GetRunID(),
&context()->event(), context()->event().event().GetEventID());
@@ -120,7 +120,7 @@ void Geant4TestEventAction::begin(const G4Event* evt) {
/// End-of-event callback
void Geant4TestEventAction::end(const G4Event* evt) {
- PRINT(name(), "%s> calling end(event_id=%d Context: run=%p (%d) evt=%p (%d))",
+ PRINT("%s> calling end(event_id=%d Context: run=%p (%d) evt=%p (%d))",
m_type.c_str(), evt->GetEventID(), &context()->run(), &context()->event(),
&context()->run(), context()->run().run().GetRunID(),
&context()->event(), context()->event().event().GetEventID());
@@ -128,14 +128,14 @@ void Geant4TestEventAction::end(const G4Event* evt) {
/// begin-of-run callback
void Geant4TestEventAction::beginRun(const G4Run* run) {
- PRINT(name(), "%s> calling beginRun(run_id=%d,num_event=%d Context:%p)",
+ PRINT("%s> calling beginRun(run_id=%d,num_event=%d Context:%p)",
m_type.c_str(), run->GetRunID(),
run->GetNumberOfEventToBeProcessed(), &context()->run());
}
/// End-of-run callback
void Geant4TestEventAction::endRun(const G4Run* run) {
- PRINT(name(), "%s> calling endRun(run_id=%d, num_event=%d Context:%p)",
+ PRINT("%s> calling endRun(run_id=%d, num_event=%d Context:%p)",
m_type.c_str(), run->GetRunID(),
run->GetNumberOfEvent(), &context()->run());
}
@@ -152,7 +152,7 @@ Geant4TestTrackAction::~Geant4TestTrackAction() {
}
/// Begin-of-tracking callback
void Geant4TestTrackAction::begin(const G4Track* trk) {
- PRINT(name(), "%s> calling begin(track=%d, parent=%d, position=(%f,%f,%f) Context: run=%p evt=%p)",
+ PRINT("%s> calling begin(track=%d, parent=%d, position=(%f,%f,%f) Context: run=%p evt=%p)",
m_type.c_str(), trk->GetTrackID(),
trk->GetParentID(), trk->GetPosition().x(), trk->GetPosition().y(), trk->GetPosition().z(),
&context()->run(), &context()->event());
@@ -160,7 +160,7 @@ void Geant4TestTrackAction::begin(const G4Track* trk) {
/// End-of-tracking callback
void Geant4TestTrackAction::end(const G4Track* trk) {
- PRINT(name(), "%s> calling end(track=%d, parent=%d, position=(%f,%f,%f) Context: run=%p evt=%p)",
+ PRINT("%s> calling end(track=%d, parent=%d, position=(%f,%f,%f) Context: run=%p evt=%p)",
m_type.c_str(), trk->GetTrackID(),
trk->GetParentID(), trk->GetPosition().x(), trk->GetPosition().y(), trk->GetPosition().z(),
&context()->run(), &context()->event());
@@ -178,7 +178,7 @@ Geant4TestStepAction::~Geant4TestStepAction() {
}
/// User stepping callback
void Geant4TestStepAction::operator()(const G4Step*, G4SteppingManager*) {
- PRINT(name(), "%s> calling operator()", m_type.c_str());
+ PRINT("%s> calling operator()", m_type.c_str());
}
/// Standard constructor with initializing arguments
@@ -186,7 +186,7 @@ Geant4TestSensitive::Geant4TestSensitive(Geant4Context* c, const std::string& n,
: Geant4Sensitive(c, n, det, lcdd), Geant4TestBase(this, "Geant4TestSensitive") {
InstanceCount::increment(this);
m_collectionID = defineCollection < TestHit > (n);
- PRINT(name(), "%s> Collection ID is %d", m_type.c_str(), int(m_collectionID));
+ PRINT("%s> Collection ID is %d", m_type.c_str(), int(m_collectionID));
}
/// Default destructor
@@ -197,7 +197,7 @@ Geant4TestSensitive::~Geant4TestSensitive() {
/// Begin-of-tracking callback
void Geant4TestSensitive::begin(G4HCofThisEvent* hce) {
Geant4HitCollection* c = collectionByID(m_collectionID);
- PRINT(name(), "%s> calling begin(num_coll=%d, coll=%s Context: run=%p evt=%p)",
+ PRINT("%s> calling begin(num_coll=%d, coll=%s Context: run=%p evt=%p)",
m_type.c_str(), hce->GetNumberOfCollections(),
c ? c->GetName().c_str() : "None", &context()->run(), &context()->event());
}
@@ -205,7 +205,7 @@ void Geant4TestSensitive::begin(G4HCofThisEvent* hce) {
/// End-of-tracking callback
void Geant4TestSensitive::end(G4HCofThisEvent* hce) {
Geant4HitCollection* c = collection(m_collectionID);
- PRINT(name(), "%s> calling end(num_coll=%d, coll=%s Context: run=%p evt=%p)",
+ PRINT("%s> calling end(num_coll=%d, coll=%s Context: run=%p evt=%p)",
m_type.c_str(), hce->GetNumberOfCollections(),
c ? c->GetName().c_str() : "None", &context()->run(), &context()->event());
}
@@ -213,7 +213,7 @@ void Geant4TestSensitive::end(G4HCofThisEvent* hce) {
/// Method for generating hit(s) using the information of G4Step object.
bool Geant4TestSensitive::process(G4Step* step, G4TouchableHistory*) {
Geant4HitCollection* c = collection(m_collectionID);
- PRINT(name(), "%s> calling process(track=%d, dE=%f, dT=%f len=%f, First,last in Vol=(%c,%c), coll=%s Context: run=%p evt=%p)",
+ PRINT("%s> calling process(track=%d, dE=%f, dT=%f len=%f, First,last in Vol=(%c,%c), coll=%s Context: run=%p evt=%p)",
m_type.c_str(), step->GetTrack()->GetTrackID(),
step->GetTotalEnergyDeposit(), step->GetDeltaTime(),
step->GetStepLength(), step->IsFirstStepInVolume() ? 'Y' : 'N',
diff --git a/DDG4/src/Geant4Vertex.cpp b/DDG4/src/Geant4Vertex.cpp
new file mode 100644
index 000000000..ee5658ccb
--- /dev/null
+++ b/DDG4/src/Geant4Vertex.cpp
@@ -0,0 +1,61 @@
+// $Id: Geant4Hits.cpp 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+
+// Framework include files
+#include "DD4hep/Printout.h"
+#include "DD4hep/InstanceCount.h"
+#include "DDG4/Geant4Vertex.h"
+
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+/// Default destructor
+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)
+{
+ InstanceCount::increment(this);
+}
+
+/// Default constructor
+Geant4Vertex::Geant4Vertex()
+ : ref(1), x(0), y(0), z(0), time(0)
+{
+ InstanceCount::increment(this);
+}
+
+/// Default destructor
+Geant4Vertex::~Geant4Vertex() {
+ InstanceCount::decrement(this);
+}
+
+/// Assignment operator
+Geant4Vertex& Geant4Vertex::operator=(const Geant4Vertex& c) {
+ if ( this != &c ) {
+ x = c.x;
+ y = c.y;
+ z = c.z;
+ time = c.time;
+ in = c.in;
+ out = c.out;
+ }
+ return *this;
+}
+
+Geant4Vertex* Geant4Vertex::addRef() {
+ ++ref;
+ return this;
+}
+
+void Geant4Vertex::release() {
+ if ( --ref <= 0 ) delete this;
+}
diff --git a/UtilityApps/src/run_plugin.h b/UtilityApps/src/run_plugin.h
index 48a40eeec..4dc37b4f3 100644
--- a/UtilityApps/src/run_plugin.h
+++ b/UtilityApps/src/run_plugin.h
@@ -28,7 +28,7 @@ using namespace DD4hep::Geometry;
//______________________________________________________________________________
namespace {
- LCDD& dd4hep_instance(const char* name="") {
+ LCDD& dd4hep_instance(const char* /* name */ ="") {
#if 0
#include "Reflex/PluginService.h"
try {
diff --git a/doc/release.notes b/doc/release.notes
index 05774538e..bfa80c226 100644
--- a/doc/release.notes
+++ b/doc/release.notes
@@ -4,6 +4,19 @@ DD4hep ---- Release Notes
=================================
+2014/08/15 Markus Frank
+-----------------------
+ - DDG4: Impreoved LCIO handling for DDG4
+ - DDG4: Input handling:
+ - LCIO input file reading with multiple *independent* inputs
+ to support overlay, multiple interactions, etc.
+ - Primary verex smearing independent for each input
+ - Primary vertex boosts independent for each input
+ - DDG4: Output handling
+ - MC Particle handling
+ - Improved MC truth handling for produced Hits.
+ Still needs revisiting.
+
--------
| v00-10 |
--------
--
GitLab