From 3e7e0095e07d08bea94e855c3dfc9c61d05ac7f6 Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Fri, 17 Oct 2014 08:26:38 +0000
Subject: [PATCH] Improve DDG4 input handling. Add HepMC ascii reader.
---
DDCore/include/DD4hep/Primitives.h | 61 +-
DDCore/src/XML/XMLElements.cpp | 2 +-
DDEve/lcio/LCIOEventHandler.cpp | 2 -
DDEve/src/HitActors.cpp | 32 +-
DDG4/examples/readHEPMC.py | 32 +
DDG4/include/DDG4/Factories.h | 19 +-
.../include/DDG4/Geant4DetectorConstruction.h | 39 +-
DDG4/include/DDG4/Geant4EventAction.h | 15 +-
DDG4/include/DDG4/Geant4GeneratorActionInit.h | 4 +-
DDG4/include/DDG4/Geant4InputAction.h | 58 +-
DDG4/include/DDG4/Geant4InteractionMerger.h | 5 +-
DDG4/include/DDG4/Geant4Mapping.h | 5 +
DDG4/include/DDG4/Geant4ParticleGun.h | 15 +
DDG4/include/DDG4/Geant4ParticleHandler.h | 17 +
DDG4/include/DDG4/Geant4PrimaryHandler.h | 5 +-
DDG4/include/DDG4/Geant4SensDetAction.h | 18 +-
DDG4/include/DDG4/Geant4StepHandler.h | 2 +-
DDG4/include/DDG4/Geant4UIManager.h | 11 +-
DDG4/lcio/Geant4Output2LCIO.cpp | 2 +-
DDG4/lcio/HepEventReader.cpp | 232 -------
DDG4/lcio/LCIOConversions.cpp | 3 +-
DDG4/lcio/LCIOEventReader.cpp | 127 +++-
DDG4/lcio/LCIOEventReader.h | 41 +-
DDG4/lcio/LCIOFileReader.cpp | 29 +-
DDG4/lcio/LCIOInputAction.cpp | 191 ------
DDG4/lcio/LCIOStdHepReader.cpp | 35 +-
DDG4/plugins/Geant4EscapeCounter.cpp | 5 +-
...Reader.cpp => Geant4EventReaderHepEvt.cpp} | 77 ++-
DDG4/plugins/Geant4EventReaderHepMC.cpp | 632 ++++++++++++++++++
DDG4/plugins/Geant4SensDet.cpp | 10 +
DDG4/python/DDG4Dict.C | 5 +
DDG4/src/Geant4Data.cpp | 2 +-
DDG4/src/Geant4DetectorConstruction.cpp | 98 ++-
DDG4/src/Geant4Exec.cpp | 4 +-
DDG4/src/Geant4InputAction.cpp | 61 +-
DDG4/src/Geant4InteractionVertexSmear.cpp | 2 -
DDG4/src/Geant4ParticleHandler.cpp | 17 +-
DDG4/src/Geant4SensDetAction.cpp | 2 +-
UtilityApps/src/teve_display.cpp | 2 +-
cmake/DD4hepConfig.cmake.in | 4 +-
40 files changed, 1307 insertions(+), 616 deletions(-)
create mode 100644 DDG4/examples/readHEPMC.py
delete mode 100644 DDG4/lcio/HepEventReader.cpp
delete mode 100644 DDG4/lcio/LCIOInputAction.cpp
rename DDG4/plugins/{Geant4HepEventReader.cpp => Geant4EventReaderHepEvt.cpp} (75%)
create mode 100644 DDG4/plugins/Geant4EventReaderHepMC.cpp
diff --git a/DDCore/include/DD4hep/Primitives.h b/DDCore/include/DD4hep/Primitives.h
index 18559678d..719de23fa 100644
--- a/DDCore/include/DD4hep/Primitives.h
+++ b/DDCore/include/DD4hep/Primitives.h
@@ -149,6 +149,29 @@ namespace DD4hep {
}
};
+ /// Operator to select second element of a pair
+ template <typename T> struct Select2nd {
+ typedef T arg_t;
+ typedef typename T::second_type result_t;
+ /// Operator function
+ const result_t& operator()(const arg_t &p) const { return p.second; }
+ };
+ /// Generator to create Operator to select value elements of a map
+ template <typename T> Select2nd<typename T::value_type> select2nd(const T&)
+ { return Select2nd<typename T::value_type>(); }
+
+ /// Operator to select the first element of a pair
+ template <typename T> struct Select1st {
+ typedef T arg_t;
+ typedef typename T::first_type result_t;
+ /// Operator function
+ const result_t& operator()(const arg_t &p) const { return p.first; }
+ };
+ /// Generator to create Operator to select key values of a map
+ template <typename T> Select1st<typename T::value_type> select1st(const T&)
+ { return Select1st<typename T::value_type>(); }
+
+
/// map Functor to delete objects from heap
template <typename M> struct DestroyObjects {
M& object;
@@ -168,6 +191,9 @@ namespace DD4hep {
template <typename M> DestroyObjects<M> destroyObjects(M& m) {
return DestroyObjects<M>(m);
}
+ template <typename M> DestroyObjects<M> destroy2nd(M& m) {
+ return DestroyObjects<M>(m);
+ }
/// map Functor to delete objects from heap
template <typename M> struct DestroyFirst {
@@ -188,6 +214,9 @@ namespace DD4hep {
template <typename M> DestroyFirst<M> destroyFirst(M& m) {
return DestroyFirst<M>(m);
}
+ template <typename M> DestroyFirst<M> destroy1st(M& m) {
+ return DestroyFirst<M>(m);
+ }
/// Helper to delete objects from heap and reset the pointer. Saves many many lines of code
template <typename T> inline void releasePtr(T*& p) {
@@ -196,13 +225,13 @@ namespace DD4hep {
p = 0;
}
- /// Functor to delete objects from heap and reset the pointer
+ /// Functor to release objects from heap and reset the pointer
template <typename T> struct ReleaseObject {
void operator()(T& p) const {
releasePtr(p);
}
};
- /// map Functor to delete objects from heap
+ /// Map Functor to release objects from heap
template <typename M> struct ReleaseObjects {
M& object;
ReleaseObjects(M& m)
@@ -218,9 +247,37 @@ namespace DD4hep {
for_each(object.begin(),object.end(),(*this));
}
};
+ template <typename M> ReleaseObject<typename M::value_type> releaseObject(M&) {
+ return ReleaseObject<typename M::value_type>();
+ }
template <typename M> ReleaseObjects<M> releaseObjects(M& m) {
return ReleaseObjects<M>(m);
}
+ template <typename M> ReleaseObjects<M> release2nd(M& m) {
+ return ReleaseObjects<M>(m);
+ }
+
+ /// Functor to delete objects from heap and reset the pointer
+ template <typename T> struct ReferenceObject {
+ typedef T arg_t;
+ T operator()(T p) const {
+ if ( p ) p->addRef();
+ return p;
+ }
+ };
+ /// Functor to delete objects from heap and reset the pointer
+ template <typename M> struct ReferenceObjects {
+ typedef typename M::second_type result_t;
+ result_t operator()(const M& p) const {
+ return ReferenceObject<result_t>()(p.second);
+ }
+ };
+ template <typename M> ReferenceObject<M> referenceObject(M&) {
+ return ReferenceObject<typename M::value_type>();
+ }
+ template <typename M> ReferenceObjects<typename M::value_type> reference2nd(M&) {
+ return ReferenceObjects<typename M::value_type>();
+ }
/// Data structure to manipulate a bitmask held by reference and represented by an integer
diff --git a/DDCore/src/XML/XMLElements.cpp b/DDCore/src/XML/XMLElements.cpp
index 246c3a236..441d14853 100644
--- a/DDCore/src/XML/XMLElements.cpp
+++ b/DDCore/src/XML/XMLElements.cpp
@@ -889,8 +889,8 @@ static unsigned int adler32(unsigned int adler, const XmlChar* xml_buff, size_t
/// Checksum (sub-)tree of a xml document/tree
typedef unsigned int (fcn_t)(unsigned int, const XmlChar*, size_t);
unsigned int Handle_t::checksum(unsigned int param, fcn_t fcn) const {
- typedef std::map<std::string, std::string> StringMap;
#ifdef DD4HEP_USE_TINYXML
+ typedef std::map<std::string, std::string> StringMap;
TiXmlNode* n = Xml(m_node).n;
if ( n ) {
if ( 0 == fcn ) fcn = adler32;
diff --git a/DDEve/lcio/LCIOEventHandler.cpp b/DDEve/lcio/LCIOEventHandler.cpp
index d5a2c7e9c..bfee6b19f 100644
--- a/DDEve/lcio/LCIOEventHandler.cpp
+++ b/DDEve/lcio/LCIOEventHandler.cpp
@@ -99,7 +99,6 @@ EventHandler::CollectionType LCIOEventHandler::collectionType(const std::string&
/// Call functor on hit collection
size_t LCIOEventHandler::collectionLoop(const std::string& collection, DDEveHitActor& actor) {
- typedef std::vector<void*> _P;
Branches::const_iterator i = m_branches.find(collection);
if ( i != m_branches.end() ) {
LCCollection* c = (*i).second;
@@ -121,7 +120,6 @@ size_t LCIOEventHandler::collectionLoop(const std::string& collection, DDEveHitA
/// Loop over collection and extract particle data
size_t LCIOEventHandler::collectionLoop(const std::string& collection, DDEveParticleActor& actor) {
- typedef std::vector<void*> _P;
Branches::const_iterator i = m_branches.find(collection);
if ( i != m_branches.end() ) {
LCCollection* c = (*i).second;
diff --git a/DDEve/src/HitActors.cpp b/DDEve/src/HitActors.cpp
index 3b630850b..f76a48d95 100644
--- a/DDEve/src/HitActors.cpp
+++ b/DDEve/src/HitActors.cpp
@@ -122,14 +122,14 @@ void BoxsetCreator::operator()(const DDEveHit& hit) {
float s2Y = 0.5*(scale(0)*std::cos(phi)+scale(2)*std::sin(phi));
float s1Z = scale(1)/2.0;
float s2Z = s1Z;
- float coords[24]= { p.X()+s1X, p.Y()+s1Y, p.Z()-s1Z,
- p.X()+s1X, p.Y()+s1Y, p.Z()+s1Z,
- p.X()-s2X, p.Y()-s2Y, p.Z()+s2Z,
- p.X()-s2X, p.Y()-s2Y, p.Z()-s2Z,
- p.X()+s2X, p.Y()+s2Y, p.Z()-s2Z,
- p.X()+s2X, p.Y()+s2Y, p.Z()+s2Z,
- p.X()-s1X, p.Y()-s1Y, p.Z()+s1Z,
- p.X()-s1X, p.Y()-s1Y, p.Z()-s1Z};
+ float coords[24]= { float(p.X()+s1X), float(p.Y()+s1Y), float(p.Z()-s1Z),
+ float(p.X()+s1X), float(p.Y()+s1Y), float(p.Z()+s1Z),
+ float(p.X()-s2X), float(p.Y()-s2Y), float(p.Z()+s2Z),
+ float(p.X()-s2X), float(p.Y()-s2Y), float(p.Z()-s2Z),
+ float(p.X()+s2X), float(p.Y()+s2Y), float(p.Z()-s2Z),
+ float(p.X()+s2X), float(p.Y()+s2Y), float(p.Z()+s2Z),
+ float(p.X()-s1X), float(p.Y()-s1Y), float(p.Z()+s1Z),
+ float(p.X()-s1X), float(p.Y()-s1Y), float(p.Z()-s1Z) };
++count;
deposit += hit.deposit*MEV_2_GEV;
boxset->AddBox(coords);
@@ -149,14 +149,14 @@ void TowersetCreator::operator()(const DDEveHit& hit) {
float s1Z = scale(1)/2.0;
float s2Z = s1Z;
p = TVector3(hit.x*MM_2_CM-s1X, hit.y*MM_2_CM-s1Y, hit.z*MM_2_CM-s1Z);
- float coords[24]= { p.X()+s1X, p.Y()+s1Y, p.Z()-s1Z,
- p.X()+s1X, p.Y()+s1Y, p.Z()+s1Z,
- p.X()-s2X, p.Y()-s2Y, p.Z()+s2Z,
- p.X()-s2X, p.Y()-s2Y, p.Z()-s2Z,
- p.X()+s2X, p.Y()+s2Y, p.Z()-s2Z,
- p.X()+s2X, p.Y()+s2Y, p.Z()+s2Z,
- p.X()-s1X, p.Y()-s1Y, p.Z()+s1Z,
- p.X()-s1X, p.Y()-s1Y, p.Z()-s1Z};
+ float coords[24]= { float(p.X()+s1X), float(p.Y()+s1Y), float(p.Z()-s1Z),
+ float(p.X()+s1X), float(p.Y()+s1Y), float(p.Z()+s1Z),
+ float(p.X()-s2X), float(p.Y()-s2Y), float(p.Z()+s2Z),
+ float(p.X()-s2X), float(p.Y()-s2Y), float(p.Z()-s2Z),
+ float(p.X()+s2X), float(p.Y()+s2Y), float(p.Z()-s2Z),
+ float(p.X()+s2X), float(p.Y()+s2Y), float(p.Z()+s2Z),
+ float(p.X()-s1X), float(p.Y()-s1Y), float(p.Z()+s1Z),
+ float(p.X()-s1X), float(p.Y()-s1Y), float(p.Z()-s1Z) };
++count;
deposit += hit.deposit*MEV_2_GEV;
boxset->AddBox(coords);
diff --git a/DDG4/examples/readHEPMC.py b/DDG4/examples/readHEPMC.py
new file mode 100644
index 000000000..1098db629
--- /dev/null
+++ b/DDG4/examples/readHEPMC.py
@@ -0,0 +1,32 @@
+#
+#
+import os, time, DDG4
+from DDG4 import OutputLevel as Output
+from SystemOfUnits import *
+#
+#
+"""
+
+DD4hep simulation example setup using the python configuration
+
+@author M.Frank
+@version 1.0
+
+"""
+def run():
+ kernel = DDG4.Kernel()
+ lcdd = kernel.lcdd()
+ simple = DDG4.Simple(kernel)
+
+ gen = DDG4.GeneratorAction(kernel,"Geant4InputAction/Input")
+ kernel.generatorAction().adopt(gen)
+ gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt";
+ parts = gen.new_particles()
+ gen.readParticles(0,parts);
+ parts.clear()
+ print 132*'*'
+ print 132*'*'
+ gen.readParticles(1,parts);
+
+if __name__ == "__main__":
+ run()
diff --git a/DDG4/include/DDG4/Factories.h b/DDG4/include/DDG4/Factories.h
index 99fa0fdc7..200d7e501 100644
--- a/DDG4/include/DDG4/Factories.h
+++ b/DDG4/include/DDG4/Factories.h
@@ -9,9 +9,7 @@
#ifndef DDG4_FACTORIES_H
#define DDG4_FACTORIES_H
-#ifndef __CINT__
-#include "Reflex/PluginService.h"
-#endif
+#include "DD4hep/Plugins.h"
#include "RVersion.h"
// Framework include files
@@ -45,6 +43,7 @@ namespace DD4hep {
class Geant4Converter;
class Geant4Sensitive;
class Geant4UserPhysics;
+ class Geant4EventReader;
class Geant4PhysicsListActionSequence;
/// Templated factory method to invoke setup action
@@ -162,6 +161,12 @@ namespace {
}
};
+ /// Factory template to create Geant4 event reader objects
+ template <typename P> class Factory<P, DD4hep::Simulation::Geant4EventReader*(std::string)> { public:
+ static void Func(void *ret, void*, const std::vector<void*>& a, void*)
+ { *(DD4hep::Simulation::Geant4EventReader**)ret = (DD4hep::Simulation::Geant4EventReader*)new P(*(std::string*)a[0]);}
+ };
+
}
#define DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(name,func) \
@@ -213,4 +218,12 @@ namespace {
template <> long Geant4SetupAction<DD4hep::Simulation::xml_g4_setup_##name>::create(LCDD& l,const DD4hep::Simulation::Geant4Converter& e, const std::map<std::string,std::string>& a) {return func(l,e,a);} }} \
PLUGINSVC_FACTORY_WITH_ID(xml_g4_setup_##name,std::string(#name "_Geant4_action"),long(DD4hep::Geometry::LCDD*,const DD4hep::Simulation::Geant4Converter*,const std::map<std::string,std::string>*))
+/// Plugin defintion to create event reader objects
+#define DECLARE_GEANT4_EVENT_READER(name) \
+ PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::Simulation::Geant4EventReader*(std::string))
+
+/// Plugin defintion to create event reader objects
+#define DECLARE_GEANT4_EVENT_READER_NS(ns,name) typedef ns::name __##name##__; \
+ PLUGINSVC_FACTORY_WITH_ID(__##name##__,std::string(#name),DD4hep::Simulation::Geant4EventReader*(std::string))
+
#endif // DDG4_FACTORIES_H
diff --git a/DDG4/include/DDG4/Geant4DetectorConstruction.h b/DDG4/include/DDG4/Geant4DetectorConstruction.h
index 0d127d732..61ada9413 100644
--- a/DDG4/include/DDG4/Geant4DetectorConstruction.h
+++ b/DDG4/include/DDG4/Geant4DetectorConstruction.h
@@ -8,7 +8,12 @@
#ifndef DD4HEP_GEANT4DETECTORCONSTRUCTION_H
#define DD4HEP_GEANT4DETECTORCONSTRUCTION_H
+// Framework include files
#include "DD4hep/Printout.h"
+#include "DDG4/Geant4Action.h"
+#include "DDG4/Geant4GeometryInfo.h"
+
+// Geant4 include files
#include "G4VUserDetectorConstruction.hh"
/// Namespace for the AIDA detector description toolkit
@@ -27,21 +32,51 @@ namespace DD4hep {
/// Class to create Geant4 detector geometry from TGeo representation in memory
/**
+ * On demand (ie. when calling "Construct") the DD4hep geometry is converted
+ * to Geant4 with all volumes, assemblies, shapes, materials etc.
+ * The actuak work is performed by the Geant4Converter class called by this method.
+ *
* \author M.Frank
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
- class Geant4DetectorConstruction : public G4VUserDetectorConstruction {
+ class Geant4DetectorConstruction : public Geant4Action, public G4VUserDetectorConstruction {
public:
+ /// Instance accessor
+ static Geant4DetectorConstruction* instance(Geant4Kernel& kernel);
/// Initializing constructor for DDG4
Geant4DetectorConstruction(Geant4Kernel& kernel);
/// Initializing constructor for other clients
Geant4DetectorConstruction(Geometry::LCDD& lcdd);
/// Default destructor
virtual ~Geant4DetectorConstruction();
- /// Geometry construction callback
+ /// Geometry construction callback: Invoke the conversion to Geant4
G4VPhysicalVolume* Construct();
+ //@{ Accessor to the various geant4 maps after construction
+
+ /// Access to the converted materials
+ const Geant4GeometryMaps::MaterialMap& materials() const;
+ /// Access to the converted elements
+ const Geant4GeometryMaps::ElementMap& elements() const;
+ /// Access to the converted shapes
+ const Geant4GeometryMaps::SolidMap& shapes() const;
+ /// Access to the converted volumes
+ const Geant4GeometryMaps::VolumeMap& volumes() const;
+ /// Access to the converted placements
+ const Geant4GeometryMaps::PlacementMap& placements() const;
+ /// Access to the converted assemblys
+ const Geant4GeometryMaps::AssemblyMap& assemblies() const;
+
+ /// Access to the converted limit sets
+ const Geant4GeometryMaps::LimitMap& limits() const;
+ /// Access to the converted regions
+ const Geant4GeometryMaps::RegionMap& regions() const;
+ /// Access to the converted sensitive detectors
+ const Geant4GeometryMaps::SensDetMap& sensitives() const;
+
+ //@}
+
private:
/// Printlevel used for the geometry conversion
PrintLevel m_outputLevel;
diff --git a/DDG4/include/DDG4/Geant4EventAction.h b/DDG4/include/DDG4/Geant4EventAction.h
index 071f8ba5b..f9866c3d7 100644
--- a/DDG4/include/DDG4/Geant4EventAction.h
+++ b/DDG4/include/DDG4/Geant4EventAction.h
@@ -24,10 +24,17 @@ namespace DD4hep {
/// Concrete basic implementation of the Geant4 event action
/**
* The EventAction is called for every event.
- * During the callback all particles are created which form the
- * microscopic kinematic action of the particle collision.
- * This input may either origin directly from an event generator
- * program or come from file.
+ *
+ * This class is the base class for all user actions, which have
+ * to hook into the begin- and end-of-event actions.
+ * Typical use cases are the collection/computation of event
+ * related properties.
+ *
+ * Examples of this functionality may include for example:
+ * - Reset variables summing event related information in the
+ * begin-event callback.
+ * - Monitoring activities such as filling histograms
+ * from hits collected during the end-event action.
*
* \author M.Frank
* \version 1.0
diff --git a/DDG4/include/DDG4/Geant4GeneratorActionInit.h b/DDG4/include/DDG4/Geant4GeneratorActionInit.h
index dda4dc771..706c6b42b 100644
--- a/DDG4/include/DDG4/Geant4GeneratorActionInit.h
+++ b/DDG4/include/DDG4/Geant4GeneratorActionInit.h
@@ -36,9 +36,9 @@ namespace DD4hep {
* -- 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
+ * -- Geant4PrimaryMap a map of the Geant4Particles converted to G4PrimaryParticles
* to ease particle handling later.
- * -- Geant4ParticleMap the map of partcles created during the event simulation.
+ * -- Geant4ParticleMap the map of particles 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....
*
diff --git a/DDG4/include/DDG4/Geant4InputAction.h b/DDG4/include/DDG4/Geant4InputAction.h
index faa7ef808..ba4912828 100644
--- a/DDG4/include/DDG4/Geant4InputAction.h
+++ b/DDG4/include/DDG4/Geant4InputAction.h
@@ -13,6 +13,7 @@
// C/C++ include files
#include <vector>
+#include <memory>
// Forward declarations
class G4Event;
@@ -36,12 +37,23 @@ namespace DD4hep {
public:
typedef Geant4Particle Particle;
+ typedef std::vector<Particle*> Particles;
+ /// Status codes of the event reader object. Anything with NOT low-bit set is an error.
+ enum EventReaderStatus {
+ EVENT_READER_ERROR=0,
+ EVENT_READER_OK=1,
+ EVENT_READER_NO_DIRECT=2,
+ EVENT_READER_NO_PRIMARIES=4,
+ EVENT_READER_NO_FACTORY=6,
+ EVENT_READER_IO_ERROR=8
+ };
protected:
- /// File name
+ /// File name to be opened and read
std::string m_name;
- /// Flag if direct event access is supported
+ /// Flag if direct event access is supported. To be explicitly set by subclass constructors
bool m_directAccess;
-
+ /// Current event number
+ int m_currEvent;
public:
/// Initializing constructor
Geant4EventReader(const std::string& nam);
@@ -51,8 +63,21 @@ namespace DD4hep {
const std::string& name() const { return m_name; }
/// Flag if direct event access (by event sequence number) is supported (Default: false)
bool hasDirectAccess() const { return m_directAccess; }
+ /// Move to the indicated event number.
+ /** For pure sequential access, the default implementation
+ * will skip events one by one.
+ * For technologies supporting direct event access the default
+ * implementation will be empty.
+ *
+ * @return
+ */
+ virtual EventReaderStatus moveToEvent(int event_number);
+ /// Skip event. To be implemented for sequential sources
+ virtual EventReaderStatus skipEvent();
/// Read an event and fill a vector of MCParticles.
- virtual int readParticles(int event_number, std::vector<Particle*>& particles) = 0;
+ /** The additional argument
+ */
+ virtual EventReaderStatus readParticles(int event_number, Particles& particles) = 0;
};
/// Generic input action capable of using the Geant4EventReader class.
@@ -69,7 +94,7 @@ namespace DD4hep {
public:
typedef Geant4Particle Particle;
-
+ typedef std::vector<Particle*> Particles;
protected:
/// Property: input file
std::string m_input;
@@ -82,36 +107,21 @@ namespace DD4hep {
/// Event reader object
Geant4EventReader* m_reader;
+ public:
/// Read an event and return a LCCollectionVec of MCParticles.
- int readParticles(int event_number, std::vector<Particle*>& particles);
+ int readParticles(int event_number, Particles& particles);
/// helper to report Geant4 exceptions
std::string issue(int i) const;
- public:
/// Standard constructor
Geant4InputAction(Geant4Context* context, const std::string& name);
/// Default destructor
virtual ~Geant4InputAction();
+ /// Create particle vector
+ Particles* new_particles() const { return new Particles; }
/// Callback to generate primary particles
virtual void operator()(G4Event* event);
};
} /* End namespace Simulation */
} /* End namespace DD4hep */
-
-#include "DD4hep/Plugins.h"
-namespace {
- /// Factory template to create Geant4 event reader objects
- template <typename P> class Factory<P, DD4hep::Simulation::Geant4EventReader*(std::string)> { public:
- static void Func(void *ret, void*, const std::vector<void*>& a, void*)
- { *(DD4hep::Simulation::Geant4EventReader**)ret = (DD4hep::Simulation::Geant4EventReader*)new P(*(std::string*)a[0]);}
- };
-}
-/// Plugin defintion to create event reader objects
-#define DECLARE_GEANT4_EVENT_READER(name) \
- PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::Simulation::Geant4EventReader*(std::string))
-
-/// Plugin defintion to create event reader objects
-#define DECLARE_GEANT4_EVENT_READER_NS(ns,name) typedef ns::name __##name##__; \
- PLUGINSVC_FACTORY_WITH_ID(__##name##__,std::string(#name),DD4hep::Simulation::Geant4EventReader*(std::string))
-
#endif /* DD4HEP_DDG4_GEANT4INPUTACTION_H */
diff --git a/DDG4/include/DDG4/Geant4InteractionMerger.h b/DDG4/include/DDG4/Geant4InteractionMerger.h
index 32c97877a..6e2009c3c 100644
--- a/DDG4/include/DDG4/Geant4InteractionMerger.h
+++ b/DDG4/include/DDG4/Geant4InteractionMerger.h
@@ -22,7 +22,10 @@ namespace DD4hep {
class Geant4PrimaryInteraction;
/// Geant4Action to merge several independent interaction to one
- /** Geant4Action to convert the particle information to Geant4
+ /** Merge all interactions created by each \tt{Geant4InputAction} into one single
+ * record. The input records are taken from the item \tt{Geant4PrimaryEvent}
+ * and are merged into the \tt{Geant4PrimaryInteraction} object attached to the
+ * \tt{Geant4Event} event context.
*
* \author M.Frank
* \version 1.0
diff --git a/DDG4/include/DDG4/Geant4Mapping.h b/DDG4/include/DDG4/Geant4Mapping.h
index c5c54f2c5..060efd4a8 100644
--- a/DDG4/include/DDG4/Geant4Mapping.h
+++ b/DDG4/include/DDG4/Geant4Mapping.h
@@ -67,6 +67,11 @@ namespace DD4hep {
return *m_dataPtr;
}
+ /// Access to the data pointer
+ Geant4GeometryInfo* ptr() const {
+ return m_dataPtr;
+ }
+
/// Create and attach new data block. Delete old data block if present.
Geant4GeometryInfo& init();
diff --git a/DDG4/include/DDG4/Geant4ParticleGun.h b/DDG4/include/DDG4/Geant4ParticleGun.h
index efc9734f0..0e620e86c 100644
--- a/DDG4/include/DDG4/Geant4ParticleGun.h
+++ b/DDG4/include/DDG4/Geant4ParticleGun.h
@@ -23,6 +23,21 @@ namespace DD4hep {
/// Implementation of a particle gun using Geant4Particles.
/**
+ * The {\tt{Geant4ParticleGun}} is a tool to shoot a number of
+ * particles with identical properties into a given region of the
+ * detector to be simulated.
+ *
+ * The particle gun is a input source like any other and participates
+ * in the general input stage merging process like any other input
+ * e.g. from file. Hence, there may be several particle guns present
+ * each generating it's own primary vertex. Use the mask property to
+ * ensure each gun generates it's own, well identified primary vertex.
+ *
+ * There is one 'user lazyness' support though:
+ * If there is only one particle gun in use, the property 'Standalone',
+ * which by default is set to true invokes the interaction merging and he
+ * Geant4 primary generation directly.
+ *
* \author M.Frank
* \version 1.0
* \ingroup DD4HEP_SIMULATION
diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h
index 4c64b7a20..9844eb836 100644
--- a/DDG4/include/DDG4/Geant4ParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4ParticleHandler.h
@@ -33,6 +33,23 @@ namespace DD4hep {
/**
* Extract the relevant particle information during the simulation step.
*
+ * This procedure works as follows:
+ * -- At the beginning of the event generation the object registers itself as
+ * Monte-Carlo truth handler to the event context.
+ * -- At the begin of each track action a particle candidate is created and filled
+ * with all properties known at this time.
+ * -- At each stepping action a flag is set if the step produced secondaries.
+ * -- Sensitive detectors call the MC truth handler if a hit was created.
+ * This fact is remembered.
+ * -- At the end of the tracking action a first decision is taken if the candidate is to be
+ * kept for the final record.
+ * -- At the end of the event action finally all particles are reduced to the
+ * final record. This logic can be overridden by a user handler to be attached.
+ * .
+ * Any of these actions may be intercepted by a {\tt{Geant4UserParticleHandler}}
+ * attached to the particle handler.
+ * See class {\tt{Geant4UserParticleHandler}} for details.
+ *
* \author M.Frank
* \version 1.0
* \ingroup DD4HEP_SIMULATION
diff --git a/DDG4/include/DDG4/Geant4PrimaryHandler.h b/DDG4/include/DDG4/Geant4PrimaryHandler.h
index 0ea34e2d2..202e12327 100644
--- a/DDG4/include/DDG4/Geant4PrimaryHandler.h
+++ b/DDG4/include/DDG4/Geant4PrimaryHandler.h
@@ -19,7 +19,10 @@ namespace DD4hep {
namespace Simulation {
/// Geant4Action to convert the particle information to Geant4
- /**
+ /** Convert the primary interaction (object \tt{Geant4PrimaryInteraction} object
+ * attached to the \tt{Geant4Event} event context) and pass the result
+ * to Geant4 for simulation.
+ *
* \author M.Frank
* \version 1.0
* \ingroup DD4HEP_SIMULATION
diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h
index 70f8a96fd..c440b5dee 100644
--- a/DDG4/include/DDG4/Geant4SensDetAction.h
+++ b/DDG4/include/DDG4/Geant4SensDetAction.h
@@ -32,6 +32,7 @@ namespace DD4hep {
namespace Geometry {
class LCDD;
class DetElement;
+ class SensitiveDetector;
}
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
@@ -57,6 +58,9 @@ namespace DD4hep {
/// Default destructor
virtual ~Geant4ActionSD();
public:
+ /// Definition: Sensitive detector type
+ typedef Geometry::SensitiveDetector SensitiveDetector;
+
/// Initialize the usage of a hit collection. Returns the collection identifier
virtual size_t defineCollection(const std::string& name) = 0;
/// Access to the readout geometry of the sensitive detector
@@ -65,10 +69,14 @@ namespace DD4hep {
virtual G4int GetCollectionID(G4int i) = 0;
/// Is the detector active?
virtual bool isActive() const = 0;
+ /// Access to the LCDD sensitive detector handle
+ virtual SensitiveDetector sensitiveDetector() const = 0;
/// G4VSensitiveDetector internals: Access to the detector path name
virtual std::string path() const = 0;
/// G4VSensitiveDetector internals: Access to the detector path name
virtual std::string fullPath() const = 0;
+ /// Access to the sensitive type of the detector
+ virtual const std::string& sensitiveType() const = 0;
};
/// Base class to construct filters for Geant4 sensitive detectors
@@ -99,7 +107,9 @@ namespace DD4hep {
typedef Geometry::Readout Readout;
typedef Geometry::DetElement DetElement;
typedef Geometry::Segmentation Segmentation;
+ /// Definition: Sensitive detector type
typedef Geometry::SensitiveDetector SensitiveDetector;
+
typedef Geant4StepHandler StepHandler;
typedef Geant4HitCollection HitCollection;
@@ -284,7 +294,8 @@ namespace DD4hep {
SensitiveDetector m_sensitive;
/// Reference to G4 sensitive detector
Geant4ActionSD* m_detector;
-
+ /// The true sensitive type of the detector
+ std::string m_sensitiveType;
/// Create a new typed hit collection
template <typename TYPE> static Geant4HitCollection* _create(const std::string& det, const std::string& coll, Geant4Sensitive* sd) {
return new Geant4HitCollection(det, coll, sd, (TYPE*) 0);
@@ -297,6 +308,11 @@ namespace DD4hep {
/// Default destructor
virtual ~Geant4SensDetActionSequence();
+ /// Access to the sensitive type of the detector
+ virtual const std::string& sensitiveType() const {
+ return m_sensitiveType;
+ }
+
/// Called at construction time of the sensitive detector to declare all hit collections
size_t defineCollections(Geant4ActionSD* sens_det);
diff --git a/DDG4/include/DDG4/Geant4StepHandler.h b/DDG4/include/DDG4/Geant4StepHandler.h
index 834914bba..b17d69c79 100644
--- a/DDG4/include/DDG4/Geant4StepHandler.h
+++ b/DDG4/include/DDG4/Geant4StepHandler.h
@@ -31,7 +31,7 @@ namespace DD4hep {
/// Helper class to ease the extraction of information from a G4Step object.
/**
*
- * Tiny helper/utility class to easily access Geant4 step information.
+ * Helper/utility class to easily access Geant4 step information.
* Born by lazyness: Avoid typing millions of statements!
*
* \author M.Frank
diff --git a/DDG4/include/DDG4/Geant4UIManager.h b/DDG4/include/DDG4/Geant4UIManager.h
index 26bd64c48..e9e417a13 100644
--- a/DDG4/include/DDG4/Geant4UIManager.h
+++ b/DDG4/include/DDG4/Geant4UIManager.h
@@ -26,10 +26,15 @@ namespace DD4hep {
namespace Simulation {
/// Standard UI interface implementation with configuration using property options
- /** @class Geant4UIManager Geant4UIManager.h DDG4/Geant4UIManager.h
+ /** The {\tt{Geant4UIManager}} is a component attached to the {\tt{Geant4Kernel}} object.
*
- * @author M.Frank
- * @version 1.0
+ * All properties of all {\tt{Geant4Action}} instances may be exported to
+ * Geant4 messengers and {\em{may}} hence be accessible directly from the Geant4
+ * prompt. To export properties from any action, call the {\tt{enableUI()}}
+ * method of the action.
+ *
+ * \author M.Frank
+ * \version 1.0
*/
class Geant4UIManager : public Geant4Action, virtual public Geant4Call {
protected:
diff --git a/DDG4/lcio/Geant4Output2LCIO.cpp b/DDG4/lcio/Geant4Output2LCIO.cpp
index aff0ec9ad..4f3544639 100644
--- a/DDG4/lcio/Geant4Output2LCIO.cpp
+++ b/DDG4/lcio/Geant4Output2LCIO.cpp
@@ -183,7 +183,7 @@ lcio::LCCollectionVec* Geant4Output2LCIO::saveParticles(Geant4ParticleMap* parti
MCParticleImpl* q = new lcio::MCParticleImpl();
q->setPDG(p->pdgID);
- float ps_fa[3] = { p->psx/GeV, p->psy/GeV, p->psz/GeV } ;
+ float ps_fa[3] = {float(p->psx/GeV),float(p->psy/GeV),float(p->psz/GeV)};
q->setMomentum( ps_fa );
double vs_fa[3] = { p->vsx/mm, p->vsy/mm, p->vsz/mm } ;
diff --git a/DDG4/lcio/HepEventReader.cpp b/DDG4/lcio/HepEventReader.cpp
deleted file mode 100644
index a1124beba..000000000
--- a/DDG4/lcio/HepEventReader.cpp
+++ /dev/null
@@ -1,232 +0,0 @@
-// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
-//====================================================================
-// AIDA Detector description implementation for LCD
-//--------------------------------------------------------------------
-//
-//====================================================================
-
-// Framework include files
-#include "LCIOEventReader.h"
-
-// C/C++ include files
-#include <fstream>
-
-/// Namespace for the AIDA detector description toolkit
-namespace DD4hep {
-
- /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
- namespace Simulation {
-
- /// File reader to load StdHep format (ASCII)
- /**
- *
- * Class to populate Geant4 primary particles and vertices from a
- * file in HepEvent format (ASCII)
- *
- * \author P.Kostka (main author)
- * \author M.Frank (code reshuffeling into new DDG4 scheme)
- * \version 1.0
- * \ingroup DD4HEP_SIMULATION
- */
- class HepEventReader : public LCIOEventReader {
- protected:
- std::ifstream m_input;
- int m_format;
- public:
- /// Initializing constructor
- HepEventReader(const std::string& nam, int fmt);
- /// Default destructor
- virtual ~HepEventReader();
- /// Read an event and return a LCCollectionVec of MCParticles.
- virtual EVENT::LCCollection *readParticles();
- };
- } /* End namespace lcio */
-} /* End namespace DD4hep */
-
-
-#include "IMPL/LCCollectionVec.h"
-#include "IMPL/MCParticleImpl.h"
-#include "EVENT/LCIO.h"
-#include "lcio.h"
-using namespace IMPL;
-using namespace EVENT;
-using namespace DD4hep::Simulation;
-#define HEPEvt 1
-
-// Factory entry
-DECLARE_LCIO_EVENT_READER(HepEventReader)
-
-/// Initializing constructor
-HepEventReader::HepEventReader(const std::string& nam, int fmt)
-: LCIOEventReader(nam), m_format(fmt)
-{
-}
-
-/// Default destructor
-HepEventReader::~HepEventReader() {
- m_input.close();
-}
-
-/// Helper function
-MCParticleImpl* cast_particle(LCObject* obj) {
- return dynamic_cast<MCParticleImpl*>(obj);
-}
-
-/// Read an event and return a LCCollectionVec of MCParticles.
-EVENT::LCCollection *HepEventReader::readParticles() {
- static const double c_light = 299.792;// mm/ns
- //
- // Read the event, check for errors
- //
- int NHEP; // number of entries
- m_input >> NHEP;
- if( m_input.eof() ) {
- //
- // End of File :: ??? Exception ???
- // -> FG: EOF is not an exception as it happens for every file at the end !
- return 0;
- }
-
- //
- // Create a Collection Vector
- IMPL::LCCollectionVec* mcVec = new IMPL::LCCollectionVec(LCIO::MCPARTICLE);
- //
- // Loop over particles
- int ISTHEP; // status code
- int IDHEP; // PDG code
- int JMOHEP1; // first mother
- int JMOHEP2; // last mother
- int JDAHEP1; // first daughter
- int JDAHEP2; // last daughter
- double PHEP1; // px in GeV/c
- double PHEP2; // py in GeV/c
- double PHEP3; // pz in GeV/c
- double PHEP4; // energy in GeV
- double PHEP5; // mass in GeV/c**2
- double VHEP1; // x vertex position in mm
- double VHEP2; // y vertex position in mm
- double VHEP3; // z vertex position in mm
- double VHEP4; // production time in mm/c
-
- //???? Memory leak ?
- std::vector<int> *daughter1 = new std::vector<int> ();
- std::vector<int> *daughter2 = new std::vector<int> ();
-
- for( int IHEP=0; IHEP<NHEP; IHEP++ ) {
- if ( m_format == HEPEvt)
- m_input >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
- >> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5;
- else
- m_input >> ISTHEP >> IDHEP
- >> JMOHEP1 >> JMOHEP2
- >> JDAHEP1 >> JDAHEP2
- >> PHEP1 >> PHEP2 >> PHEP3
- >> PHEP4 >> PHEP5
- >> VHEP1 >> VHEP2 >> VHEP3
- >> VHEP4;
-
- if(m_input.eof())
- return 0;
- //
- // Create a MCParticle and fill it from stdhep info
- MCParticleImpl* mcp = new MCParticleImpl();
- //
- // PDGID
- mcp->setPDG(IDHEP);
- //
- // Momentum vector
- float p0[3] = {PHEP1,PHEP2,PHEP3};
- mcp->setMomentum(p0);
- //
- // Mass
- mcp->setMass(PHEP5);
- //
- // Vertex
- // (missing information in HEPEvt files)
- double v0[3] = {0.,0.,0.};
- mcp->setVertex(v0);
- //
- // Generator status
- mcp->setGeneratorStatus(ISTHEP);
- //
- // Simulator status 0 until simulator acts on it
- mcp->setSimulatorStatus(0);
- //
- // Creation time (note the units)
- // (No information in HEPEvt files)
- mcp->setTime(0./c_light);
- //
- // Add the particle to the collection vector
- mcVec->push_back(mcp);
- //
- // Keep daughters information for later
- daughter1->push_back(JDAHEP1);
- daughter2->push_back(JDAHEP2);
- }// End loop over particles
- //
- // Now make a second loop over the particles, checking the daughter
- // information. This is not always consistent with parent
- // information, and this utility assumes all parents listed are
- // parents and all daughters listed are daughters
- //
- for( int IHEP=0; IHEP<NHEP; IHEP++ ) {
- struct Particle {
- MCParticleImpl* m_part;
- Particle(MCParticleImpl* p) : m_part(p) {}
- void addParent(MCParticleImpl* p) { m_part->addParent(p); }
- MCParticleImpl* findParent(MCParticleImpl* p) {
- int np = m_part->getParents().size();
- for(int ip=0;ip < np;ip++) {
- MCParticleImpl* q = cast_particle(p->getParents()[ip]); // cast necessary ? only if virt inheritance
- if ( q == p ) return p;
- }
- return 0;
- }
- };
-
- //
- // Get the MCParticle
- //
- MCParticleImpl* mcp = cast_particle(mcVec->getElementAt(IHEP));
- //
- // Get the daughter information, discarding extra information
- // sometimes stored in daughter variables.
- //
- int fd = daughter1->operator[](IHEP) - 1;
- int ld = daughter2->operator[](IHEP) - 1;
-
- //
- // As with the parents, look for range, 2 discreet or 1 discreet
- // daughter.
- if( (fd > -1) && (ld > -1) ) {
- if(ld >= fd) {
- for(int id=fd;id<ld+1;id++) {
- //
- // Get the daughter, and see if it already lists this particle as
- // a parent. If not, add this particle as a parent
- //
- Particle part(cast_particle(mcVec->getElementAt(id)));
- if ( !part.findParent(mcp) ) part.addParent(mcp);
- }
- }
- else { // Same logic, discrete cases
- Particle part_fd(cast_particle(mcVec->getElementAt(fd)));
- if ( !part_fd.findParent(mcp) ) part_fd.addParent(mcp);
-
- Particle part_ld(cast_particle(mcVec->getElementAt(ld)));
- if ( !part_ld.findParent(mcp) ) part_ld.addParent(mcp);
- }
- }
- else if(fd > -1) {
- Particle part(cast_particle(mcVec->getElementAt(fd)));
- if ( !part.findParent(mcp) ) part.addParent(mcp);
- }
- else if(ld > -1) {
- Particle part(cast_particle(mcVec->getElementAt(ld)));
- if ( !part.findParent(mcp) ) part.addParent(mcp);
- }
- } // End second loop over particles
-
- return mcVec;
-}
-
diff --git a/DDG4/lcio/LCIOConversions.cpp b/DDG4/lcio/LCIOConversions.cpp
index 03d2611af..d873ae155 100644
--- a/DDG4/lcio/LCIOConversions.cpp
+++ b/DDG4/lcio/LCIOConversions.cpp
@@ -29,6 +29,7 @@
//
#include "UTIL/Operators.h"
#include "UTIL/ILDConf.h"
+
using namespace std;
//==================================================================================
@@ -170,7 +171,7 @@ namespace DD4hep {
}
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()};
+ float pos[3] = {float(hit->position.x()), float(hit->position.y()), float(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); // ????
diff --git a/DDG4/lcio/LCIOEventReader.cpp b/DDG4/lcio/LCIOEventReader.cpp
index a4ffa7f2e..f51da1a04 100644
--- a/DDG4/lcio/LCIOEventReader.cpp
+++ b/DDG4/lcio/LCIOEventReader.cpp
@@ -10,11 +10,134 @@
// Framework include files
#include "LCIOEventReader.h"
+#include "DD4hep/Printout.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Factories.h"
+
+#include "G4ParticleTable.hh"
+#include "EVENT/MCParticle.h"
+#include "EVENT/LCCollection.h"
+
+#include "G4Event.hh"
+
+using namespace std;
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+
+// Neede for backwards compatibility:
+namespace DD4hep{namespace Simulation{typedef Geant4InputAction LCIOInputAction;}}
+DECLARE_GEANT4ACTION(LCIOInputAction)
+
+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;
+ }
+}
+
/// Initializing constructor
-DD4hep::Simulation::LCIOEventReader::LCIOEventReader(const std::string& nam) : m_name(nam) {
+LCIOEventReader::LCIOEventReader(const string& nam)
+ : Geant4EventReader(nam), m_numEvent(0)
+{
}
/// Default destructor
-DD4hep::Simulation::LCIOEventReader::~LCIOEventReader() {
+LCIOEventReader::~LCIOEventReader() {
+}
+
+/// Read an event and fill a vector of MCParticles.
+LCIOEventReader::EventReaderStatus
+LCIOEventReader::readParticles(int event_number, vector<Particle*>& particles) {
+ EVENT::LCCollection* primaries = 0;
+ map<EVENT::MCParticle*,int> mcparts;
+ vector<EVENT::MCParticle*> mcpcoll;
+ EventReaderStatus ret = EVENT_READER_OK;
+
+ if ( hasDirectAccess() )
+ ret = readParticles(event_number,&primaries);
+ else if ( m_numEvent == event_number )
+ ret = readParticles(event_number,&primaries);
+ else
+ ret = EVENT_READER_NO_DIRECT;
+
+ ++m_numEvent;
+ if ( ret != EVENT_READER_OK ) return ret;
+
+ int NHEP = primaries->getNumberOfElements();
+ // check if there is at least one particle
+ if ( NHEP == 0 ) return EVENT_READER_NO_PRIMARIES;
+
+ 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;
+ }
+
+ // build collection of MCParticles
+ G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
+ for(int i=0; i<NHEP; ++i ) {
+ EVENT::MCParticle* mcp = mcpcoll[i];
+ Geant4ParticleHandle p(new Particle(i));
+ 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->pdgID = mcp->getPDG();
+ p->charge = int(mcp->getCharge()*3.0);
+ p->psx = mom[0]*GeV;
+ p->psy = mom[1]*GeV;
+ p->psz = 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();
+ if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
+ else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE);
+ else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
+ else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ else {
+ cout << " #### WARNING - LCIOInputAction : unknown generator status : "
+ << genStatus << " -> ignored ! " << endl;
+ }
+ 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);
+ particles.push_back(p);
+ }
+ return EVENT_READER_OK;
}
+
diff --git a/DDG4/lcio/LCIOEventReader.h b/DDG4/lcio/LCIOEventReader.h
index 8cb923d73..5d2e67626 100644
--- a/DDG4/lcio/LCIOEventReader.h
+++ b/DDG4/lcio/LCIOEventReader.h
@@ -7,8 +7,8 @@
#ifndef DD4HEP_DDG4_LCIOEVENTREADER_H
#define DD4HEP_DDG4_LCIOEVENTREADER_H
-// C/C++ include files
-#include <string>
+// Framework include files
+#include "DDG4/Geant4InputAction.h"
// Forward declarations
namespace EVENT { class LCCollection; }
@@ -19,9 +19,6 @@ namespace DD4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace Simulation {
- // Forward declarations
- class Geant4Particle;
-
/// Base class to read lcio files.
/**
* \author P.Kostka (main author)
@@ -29,43 +26,21 @@ namespace DD4hep {
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
- class LCIOEventReader {
-
- public:
- typedef Geant4Particle Particle;
+ class LCIOEventReader : public Geant4EventReader {
protected:
- /// File name
- std::string m_name;
-
+ /// Event counter
+ int m_numEvent;
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 fill a vector of MCParticles.
+ virtual EventReaderStatus readParticles(int event_number, std::vector<Particle*>& particles);
/// Read an event and return a LCCollectionVec of MCParticles.
- virtual EVENT::LCCollection *readParticles() = 0;
+ virtual EventReaderStatus readParticles(int event_number, EVENT::LCCollection** particles) = 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/LCIOFileReader.cpp b/DDG4/lcio/LCIOFileReader.cpp
index 0b3db53bd..307b877fe 100644
--- a/DDG4/lcio/LCIOFileReader.cpp
+++ b/DDG4/lcio/LCIOFileReader.cpp
@@ -32,30 +32,27 @@ namespace DD4hep {
IO::LCReader* m_reader;
public:
/// Initializing constructor
- LCIOFileReader(const std::string& nam, int);
+ LCIOFileReader(const std::string& nam);
/// Default destructor
virtual ~LCIOFileReader();
- /// Read an event and return a LCCollectionVec of MCParticles.
- virtual EVENT::LCCollection *readParticles();
+ /// Read an event and fill a vector of MCParticles.
+ virtual EventReaderStatus readParticles(int event_number, EVENT::LCCollection** particles);
};
}
}
#endif // DD4HEP_DDG4_LCIOFILEREADER_H
#include "DD4hep/Printout.h"
-#include "DD4hep/Primitives.h"
-#include "lcio.h"
-#include "EVENT/LCIO.h"
-#include "IO/LCReader.h"
+#include "DDG4/Factories.h"
+#include "UTIL/ILDConf.h"
-using namespace EVENT;
using namespace DD4hep::Simulation;
// Factory entry
-DECLARE_LCIO_EVENT_READER(LCIOFileReader)
+DECLARE_GEANT4_EVENT_READER_NS(DD4hep::Simulation,LCIOFileReader)
/// Initializing constructor
-LCIOFileReader::LCIOFileReader(const std::string& nam, int /* arg */)
+DD4hep::Simulation::LCIOFileReader::LCIOFileReader(const std::string& nam)
: LCIOEventReader(nam)
{
m_reader = ::lcio::LCFactory::getInstance()->createLCReader();
@@ -64,16 +61,18 @@ LCIOFileReader::LCIOFileReader(const std::string& nam, int /* arg */)
}
/// Default destructor
-LCIOFileReader::~LCIOFileReader() {
+DD4hep::Simulation::LCIOFileReader::~LCIOFileReader() {
DD4hep::deletePtr(m_reader);
}
-/// Read an event and return a LCCollectionVec of MCParticles.
-EVENT::LCCollection *LCIOFileReader::readParticles() {
+/// Read an event and fill a vector of MCParticles.
+Geant4EventReader::EventReaderStatus
+DD4hep::Simulation::LCIOFileReader::readParticles(int /*event_number */, EVENT::LCCollection** particles) {
::lcio::LCEvent* evt = m_reader->readNextEvent();
if ( evt ) {
- return evt->getCollection(LCIO::MCPARTICLE);
+ *particles = evt->getCollection(LCIO::MCPARTICLE);
+ if ( *particles ) return EVENT_READER_OK;
}
- return 0;
+ return EVENT_READER_ERROR;
}
diff --git a/DDG4/lcio/LCIOInputAction.cpp b/DDG4/lcio/LCIOInputAction.cpp
deleted file mode 100644
index cd95c7e65..000000000
--- a/DDG4/lcio/LCIOInputAction.cpp
+++ /dev/null
@@ -1,191 +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 "LCIOInputAction.h"
-#include "LCIOEventReader.h"
-#include "DD4hep/Printout.h"
-#include "DDG4/Geant4Primary.h"
-#include "DDG4/Geant4Context.h"
-#include "G4ParticleTable.hh"
-#include "EVENT/MCParticle.h"
-#include "EVENT/LCCollection.h"
-
-#include "G4Event.hh"
-
-using namespace std;
-using namespace DD4hep::Simulation;
-typedef DD4hep::ReferenceBitMask<int> PropertyMask;
-
-namespace {
- inline int GET_ENTRY(const map<EVENT::MCParticle*,int>& mcparts, EVENT::MCParticle* part) {
- map<EVENT::MCParticle*,int>::const_iterator ip=mcparts.find(part);
- if ( ip == mcparts.end() ) {
- throw runtime_error("Unknown particle identifier look-up!");
- }
- return (*ip).second;
- }
-}
-
-/// Standard constructor
-LCIOInputAction::LCIOInputAction(Geant4Context* context, const string& name)
- : Geant4GeneratorAction(context,name), m_reader(0)
-{
- 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;
-}
-
-/// 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;
-
- 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;
-
- 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;
- }
-
- print("+++ Particle interaction with %d generator particles ++++++++++++++++++++++++",NHEP);
- Geant4Vertex* vtx = new Geant4Vertex();
- vtx->x = 0;
- vtx->y = 0;
- vtx->z = 0;
- vtx->time = 0;
- inter->vertices.insert(make_pair(m_mask,vtx));
- // build collection of MCParticles
- G4ParticleTable* tab = G4ParticleTable::GetParticleTable();
- for(int i=0; i<NHEP; ++i ) {
- EVENT::MCParticle* mcp = mcpcoll[i];
- Geant4ParticleHandle p(new Particle(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->pdgID = mcp->getPDG();
- p->charge = int(mcp->getCharge()*3.0);
- 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();
- if ( genStatus == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
- else if ( genStatus == 1 ) status.set(G4PARTICLE_GEN_STABLE);
- else if ( genStatus == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
- else if ( genStatus == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
- else { std::cout << " #### WARNING - LCIOInputAction : unknown generator status : " << genStatus << " -> ignored ! " << std::endl ; }
-
- 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);
- if ( par.size() == 0 ) {
- if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
- vtx->in.insert(p->id); // Beam particles and primary quarks etc.
- else
- vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
- }
- inter->particles.insert(make_pair(p->id,p));
- p.dump3(outputLevel()-1,name(),"+->");
- }
-}
-
-#include "DDG4/Factories.h"
-DECLARE_GEANT4ACTION(LCIOInputAction)
diff --git a/DDG4/lcio/LCIOStdHepReader.cpp b/DDG4/lcio/LCIOStdHepReader.cpp
index 2458ba353..626b8838b 100644
--- a/DDG4/lcio/LCIOStdHepReader.cpp
+++ b/DDG4/lcio/LCIOStdHepReader.cpp
@@ -6,11 +6,11 @@
//====================================================================
#ifndef DD4HEP_DDG4_LCIOSTDHEPREADER_H
#define DD4HEP_DDG4_LCIOSTDHEPREADER_H
+
// Framework include files
#include "LCIOEventReader.h"
-
-// Forward declarations
-namespace UTIL { class LCStdHepRdr; }
+// LCIO include files
+#include "UTIL/LCStdHepRdr.h"
/// Namespace for the AIDA detector description toolkit
namespace DD4hep {
@@ -31,31 +31,27 @@ namespace DD4hep {
UTIL::LCStdHepRdr* m_reader;
public:
/// Initializing constructor
- LCIOStdHepReader(const std::string& nam, int);
+ LCIOStdHepReader(const std::string& nam);
/// Default destructor
virtual ~LCIOStdHepReader();
- /// Read an event and return a LCCollectionVec of MCParticles.
- virtual EVENT::LCCollection *readParticles();
+ /// Read an event and fill a vector of MCParticles.
+ virtual EventReaderStatus readParticles(int event_number, EVENT::LCCollection** particles);
};
} /* End namespace lcio */
} /* End namespace DD4hep */
#endif /* DD4HEP_DDG4_LCIOSTDHEPREADER_H */
-#include "lcio.h"
-#include "EVENT/LCIO.h"
-#include "UTIL/LCStdHepRdr.h"
-#include "DD4hep/Primitives.h"
+#include "DDG4/Factories.h"
+// Factory entry
+DECLARE_GEANT4_EVENT_READER_NS(DD4hep::Simulation,LCIOStdHepReader)
using namespace DD4hep::Simulation;
-// Factory entry
-DECLARE_LCIO_EVENT_READER(LCIOStdHepReader)
-
/// Initializing constructor
-LCIOStdHepReader::LCIOStdHepReader(const std::string& nam, int)
+LCIOStdHepReader::LCIOStdHepReader(const std::string& nam)
: LCIOEventReader(nam)
{
- m_reader = new ::lcio::LCStdHepRdr(m_name.c_str());
+ m_reader = new UTIL::LCStdHepRdr(m_name.c_str());
}
/// Default destructor
@@ -63,7 +59,10 @@ LCIOStdHepReader::~LCIOStdHepReader() {
DD4hep::deletePtr(m_reader);
}
-/// Read an event and return a LCCollectionVec of MCParticles.
-EVENT::LCCollection *LCIOStdHepReader::readParticles() {
- return m_reader->readEvent();
+/// Read an event and fill a vector of MCParticles.
+Geant4EventReader::EventReaderStatus
+LCIOStdHepReader::readParticles(int /*event_number */, EVENT::LCCollection** particles) {
+ *particles = m_reader->readEvent();
+ if ( 0 == *particles ) return EVENT_READER_ERROR;
+ return EVENT_READER_OK;
}
diff --git a/DDG4/plugins/Geant4EscapeCounter.cpp b/DDG4/plugins/Geant4EscapeCounter.cpp
index bbab87fbf..8bc974034 100644
--- a/DDG4/plugins/Geant4EscapeCounter.cpp
+++ b/DDG4/plugins/Geant4EscapeCounter.cpp
@@ -74,8 +74,7 @@ using namespace DD4hep::Simulation;
/// Standard constructor
Geant4EscapeCounter::Geant4EscapeCounter(Geant4Context* context, const string& nam, DetElement det, LCDD& lcdd)
- // : Geant4SteppingAction(context, nam)
- : Geant4Sensitive(context, nam, det, lcdd)
+: Geant4Sensitive(context, nam, det, lcdd)
{
string coll_name = name()+"Hits";
m_needsControl = true;
@@ -91,8 +90,6 @@ Geant4EscapeCounter::~Geant4EscapeCounter() {
/// G4VSensitiveDetector interface: Method for generating hit(s) using the information of G4Step object.
bool Geant4EscapeCounter::process(G4Step* step, G4TouchableHistory* /* history */) {
- typedef SimpleHit::Contribution HitContribution;
- typedef vector<string> _V;
Geant4StepHandler h(step);
Geant4TrackHandler th(h.track);
Geant4TouchableHandler handler(step);
diff --git a/DDG4/plugins/Geant4HepEventReader.cpp b/DDG4/plugins/Geant4EventReaderHepEvt.cpp
similarity index 75%
rename from DDG4/plugins/Geant4HepEventReader.cpp
rename to DDG4/plugins/Geant4EventReaderHepEvt.cpp
index 666a91af5..ff4c3e053 100644
--- a/DDG4/plugins/Geant4HepEventReader.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepEvt.cpp
@@ -20,25 +20,26 @@ namespace DD4hep {
/// Class to populate Geant4 primaries from StdHep files.
/**
* Class to populate Geant4 primary particles and vertices from a
- * file in StdHep format (ASCII)
+ * file in HEPEvt format (ASCII)
*
* \author P.Kostka (main author)
* \author M.Frank (code reshuffeling into new DDG4 scheme)
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
- class Geant4HepEventReader : public Geant4EventReader {
+ class Geant4EventReaderHepEvt : public Geant4EventReader {
+
protected:
std::ifstream m_input;
- int m_format;
+ int m_format;
public:
/// Initializing constructor
- Geant4HepEventReader(const std::string& nam);
+ Geant4EventReaderHepEvt(const std::string& nam, int format);
/// Default destructor
- virtual ~Geant4HepEventReader();
+ virtual ~Geant4EventReaderHepEvt();
/// Read an event and fill a vector of MCParticles.
- virtual int readParticles(int event_number, std::vector<Particle*>& particles);
+ virtual EventReaderStatus readParticles(int event_number, std::vector<Particle*>& particles);
};
} /* End namespace Simulation */
} /* End namespace DD4hep */
@@ -49,38 +50,72 @@ namespace DD4hep {
//--------------------------------------------------------------------
//
//====================================================================
-// #include "DDG4/Geant4HepEventReader"
+// #include "DDG4/Geant4EventReaderHepEvt"
+#include "DDG4/Factories.h"
+#include "DD4hep/Printout.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include <cerrno>
+using namespace std;
using namespace CLHEP;
using namespace DD4hep::Simulation;
typedef DD4hep::ReferenceBitMask<int> PropertyMask;
-#define HEPEvt 1
+#define HEPEvtShort 1
+#define HEPEvtLong 2
+namespace {
+ class Geant4EventReaderHepEvtShort : public Geant4EventReaderHepEvt {
+ public:
+ /// Initializing constructor
+ Geant4EventReaderHepEvtShort(const string& nam) : Geant4EventReaderHepEvt(nam,HEPEvtShort) {}
+ /// Default destructor
+ virtual ~Geant4EventReaderHepEvtShort() {}
+ };
+ class Geant4EventReaderHepEvtLong : public Geant4EventReaderHepEvt {
+ public:
+ /// Initializing constructor
+ Geant4EventReaderHepEvtLong(const string& nam) : Geant4EventReaderHepEvt(nam,HEPEvtLong) {}
+ /// Default destructor
+ virtual ~Geant4EventReaderHepEvtLong() {}
+ };
+}
// Factory entry
-DECLARE_GEANT4_EVENT_READER(Geant4HepEventReader)
+DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepEvtShort)
+// Factory entry
+DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepEvtLong)
+
/// Initializing constructor
-Geant4HepEventReader::Geant4HepEventReader(const std::string& nam)
-: Geant4EventReader(nam), m_input(), m_format(HEPEvt)
+Geant4EventReaderHepEvt::Geant4EventReaderHepEvt(const string& nam, int format)
+: Geant4EventReader(nam), m_input(), m_format(format)
{
- // Need to set format from input specification!!!
-
// Now open the input file:
+ m_input.open(nam.c_str(),ifstream::in);
+ if ( !m_input.good() ) {
+ string err = "+++ Geant4EventReaderHepEvt: Failed to open input stream:"+nam+
+ " Error:"+string(strerror(errno));
+ throw runtime_error(err);
+ }
}
/// Default destructor
-Geant4HepEventReader::~Geant4HepEventReader() {
+Geant4EventReaderHepEvt::~Geant4EventReaderHepEvt() {
m_input.close();
}
/// Read an event and fill a vector of MCParticles.
-int Geant4HepEventReader::readParticles(int /* event_number */, std::vector<Particle*>& particles) {
+Geant4EventReader::EventReaderStatus
+Geant4EventReaderHepEvt::readParticles(int /* event_number */, vector<Particle*>& particles) {
+
+ // First check the input file status
+ if ( !m_input.good() || m_input.eof() ) {
+ return EVENT_READER_IO_ERROR;
+ }
+
//static const double c_light = 299.792;// mm/ns
//
// Read the event, check for errors
@@ -91,7 +126,7 @@ int Geant4HepEventReader::readParticles(int /* event_number */, std::vector<Part
//
// End of File :: ??? Exception ???
// -> FG: EOF is not an exception as it happens for every file at the end !
- return EIO;
+ return EVENT_READER_IO_ERROR;
}
//
// Loop over particles
@@ -111,11 +146,11 @@ int Geant4HepEventReader::readParticles(int /* event_number */, std::vector<Part
double VHEP3; // z vertex position in mm
double VHEP4; // production time in mm/c
- std::vector<int> daughter1;
- std::vector<int> daughter2;
+ vector<int> daughter1;
+ vector<int> daughter2;
for( int IHEP=0; IHEP<NHEP; IHEP++ ) {
- if ( m_format == HEPEvt)
+ if ( m_format == HEPEvtShort )
m_input >> ISTHEP >> IDHEP >> JDAHEP1 >> JDAHEP2
>> PHEP1 >> PHEP2 >> PHEP3 >> PHEP5;
else
@@ -128,7 +163,7 @@ int Geant4HepEventReader::readParticles(int /* event_number */, std::vector<Part
>> VHEP4;
if(m_input.eof())
- return EIO;
+ return EVENT_READER_IO_ERROR;
//
// Create a MCParticle and fill it from stdhep info
Particle* p = new Particle(IHEP);
@@ -236,6 +271,6 @@ int Geant4HepEventReader::readParticles(int /* event_number */, std::vector<Part
if ( !part.findParent(mcp) ) part.addParent(mcp);
}
} // End second loop over particles
- return 0;
+ return EVENT_READER_OK;
}
diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp
new file mode 100644
index 000000000..07b35190c
--- /dev/null
+++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp
@@ -0,0 +1,632 @@
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+//====================================================================
+
+// Framework include files
+#include "DD4hep/IoStreams.h"
+#include "DDG4/Geant4InputAction.h"
+
+// C/C++ include files
+#include <fstream>
+#include <boost/iostreams/stream.hpp>
+
+/// Namespace for the AIDA detector description toolkit
+namespace DD4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace Simulation {
+
+ namespace HepMC { class EventStream; }
+ /// Class to populate Geant4 primaries from StdHep files.
+ /**
+ * Class to populate Geant4 primary particles and vertices from a
+ * file in HEPEvt format (ASCII)
+ *
+ * \author P.Kostka (main author)
+ * \author M.Frank (code reshuffeling into new DDG4 scheme)
+ * \version 1.0
+ * \ingroup DD4HEP_SIMULATION
+ */
+ class Geant4EventReaderHepMC : public Geant4EventReader {
+ typedef boost::iostreams::stream<dd4hep_file_source<int> > in_stream;
+ //typedef boost::iostreams::stream<dd4hep_file_source<TFile*> > in_stream;
+ typedef HepMC::EventStream EventStream;
+ protected:
+ in_stream m_input;
+ EventStream* m_events;
+ public:
+ /// Initializing constructor
+ Geant4EventReaderHepMC(const std::string& nam);
+ /// Default destructor
+ virtual ~Geant4EventReaderHepMC();
+ /// Read an event and fill a vector of MCParticles.
+ virtual EventReaderStatus readParticles(int event_number, std::vector<Particle*>& particles);
+ };
+ } /* End namespace Simulation */
+} /* End namespace DD4hep */
+
+// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+//====================================================================
+// #include "DDG4/Geant4EventReaderHepMC"
+
+#include "DDG4/Factories.h"
+#include "DD4hep/Printout.h"
+#include "DDG4/Geant4Primary.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+
+#include <cerrno>
+#include <algorithm>
+
+using namespace std;
+using namespace CLHEP;
+using namespace DD4hep::Simulation;
+typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+
+// Factory entry
+DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepMC)
+
+/// Namespace for the AIDA detector description toolkit
+namespace DD4hep {
+
+ /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
+ namespace Simulation {
+
+ namespace HepMC {
+
+ class EventHeader {
+ public:
+ int id;
+ int num_vertices;
+ int bp1, bp2;
+ int signal_process_id;
+ int signal_process_vertex;
+ float scale;
+ float alpha_qcd;
+ float alpha_qed;
+ vector<float> weights;
+ vector<long> random;
+ EventHeader() {}
+ };
+
+ /// The known_io enum is used to track which type of input is being read
+ enum known_io { gen=1, ascii, extascii, ascii_pdt, extascii_pdt };
+
+ class EventStream {
+ public:
+ istream& instream;
+ // io information
+ string key;
+ double mom_unit, pos_unit;
+ int io_type;
+ float xsection, xsection_err;
+ //WeightContainer weights;
+ typedef std::map<int,Geant4Vertex*> Vertices;
+ Vertices m_vertices;
+ typedef std::map<int,Geant4Particle*> Particles;
+ Particles m_particles;
+
+ EventHeader header;
+
+ EventStream(istream& in) : instream(in), io_type(0)
+ { use_default_units(); }
+ Particles& particles() { return m_particles; }
+ Vertices& vertices() { return m_vertices; }
+ void add_vertex(int id, Geant4Vertex* v)
+ { m_vertices.insert(make_pair(id,v)); }
+ void add_particle(Geant4Particle* p);
+ void set_io(int typ, const string& k)
+ { io_type = typ; key = k; }
+ void use_default_units()
+ { mom_unit = MeV; pos_unit = mm; }
+ Geant4Vertex* vertex(int i);
+ void fix_particles();
+ bool read();
+ void clear();
+ };
+
+ char get_input(istream& is, istringstream& iline);
+ int find_event_end(istream & is);
+ int read_weight_names(EventStream &, istringstream& iline);
+ int read_particle(EventStream &info, istringstream& iline, Geant4Particle * p);
+ int read_vertex(EventStream &info, istream& is, istringstream & iline);
+ int read_event_header(EventStream &info, istringstream & input, EventHeader& header);
+ int read_cross_section(EventStream &info, istringstream & input);
+ int read_units(EventStream &info, istringstream & input);
+ int read_heavy_ion(EventStream &, istringstream & input);
+ int read_pdf(EventStream &, istringstream & input);
+ }
+ }
+}
+
+/// Initializing constructor
+Geant4EventReaderHepMC::Geant4EventReaderHepMC(const string& nam)
+: Geant4EventReader(nam), m_input(), m_events(0)
+{
+ // Now open the input file:
+ m_input.open(nam.c_str(),BOOST_IOS::in|BOOST_IOS::binary);
+ if ( m_input < 0 ) {
+ string err = "+++ Geant4EventReaderHepMC: Failed to open input stream:"+nam+
+ " Error:"+string(strerror(errno));
+ throw runtime_error(err);
+ }
+ m_events = new HepMC::EventStream(m_input);
+}
+
+/// Default destructor
+Geant4EventReaderHepMC::~Geant4EventReaderHepMC() {
+ delete m_events;
+ m_events = 0;
+ m_input.close();
+}
+namespace DD4hep{
+ template <typename M> ReferenceObjects<typename M::value_type> __reference2nd(M&) {
+ return ReferenceObjects<typename M::value_type>();
+ }
+}
+/// Read an event and fill a vector of MCParticles.
+Geant4EventReaderHepMC::EventReaderStatus
+Geant4EventReaderHepMC::readParticles(int /* ev_id */, Particles& output) {
+ // make sure the stream is good
+ if ( m_input.eof() || m_input.fail() ) {
+ cerr << "streaming input: end of stream found setting badbit." << endl;
+ m_input.clear(ios::badbit);
+ return EVENT_READER_IO_ERROR;
+ }
+
+ if ( m_events->read() ) {
+ EventStream::Particles& p = m_events->particles();
+ output.reserve(p.size());
+ transform(p.begin(),p.end(),back_inserter(output),reference2nd(p));
+ m_events->clear();
+ for(Particles::const_iterator k=output.begin(); k != output.end(); ++k) {
+ Geant4ParticleHandle p(*k);
+ printout(INFO,m_name,
+ "+++ %s ID:%3d status:%08X typ:%9d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
+ "time: %+.2e [ns] #Dau:%3d #Par:%1d",
+ "",p->id,p->status,p->pdgID,
+ p->psx/MeV,p->psy/MeV,p->psz/MeV,p->time/ns,
+ p->daughters.size(),
+ p->parents.size());
+ //output.push_back(p);
+ }
+ return EVENT_READER_OK;
+ }
+ return EVENT_READER_IO_ERROR;
+}
+
+
+void HepMC::EventStream::fix_particles() {
+ EventStream::Particles::iterator i;
+ std::set<int>::const_iterator id, ip;
+ for(i=m_particles.begin(); i != m_particles.end(); ++i) {
+ Geant4ParticleHandle p((*i).second);
+ int end_vtx_id = p->secondaries;
+ p->secondaries = 0;
+ Geant4Vertex* v = vertex(end_vtx_id);
+ if ( v ) {
+ p->vex = v->x;
+ p->vey = v->y;
+ p->vez = v->z;
+ v->in.insert(p->id);
+ for(id=v->out.begin(); id!=v->out.end();++id)
+ p->daughters.insert(*id);
+ }
+ }
+ EventStream::Vertices::iterator j;
+ for(j=m_vertices.begin(); j != m_vertices.end(); ++j) {
+ Geant4Vertex* v = (*j).second;
+ for (id=v->in.begin(); id!=v->in.end();++id) {
+ for (ip=v->out.begin(); ip!=v->out.end();++ip) {
+ EventStream::Particles::iterator ipp = m_particles.find(*ip);
+ (*ipp).second->parents.insert(*id);
+ }
+ }
+ }
+}
+
+char HepMC::get_input(istream& is, istringstream& iline) {
+ char value = is.peek();
+ if ( !is ) { // make sure the stream is valid
+ cerr << "StreamHelpers: setting badbit." << endl;
+ is.clear(ios::badbit);
+ return -1;
+ }
+ string line, firstc;
+ getline(is,line);
+ if ( !is ) { // make sure the stream is valid
+ cerr << "StreamHelpers: setting badbit." << endl;
+ is.clear(ios::badbit);
+ return -1;
+ }
+ iline.str(line);
+ if ( line.length()>0 && line[1] == ' ' ) iline >> firstc;
+ return iline ? value : -1;
+}
+
+int HepMC::find_event_end(istream & is) {
+ string line;
+ while ( is ) {
+ char val = is.peek();
+ if( val == 'E' ) { // next event
+ return 1;
+ }
+ getline(is,line);
+ if ( !is.good() ) return 0;
+ }
+ return 0;
+}
+
+int HepMC::read_weight_names(EventStream&, istringstream&) {
+#if 0
+ int HepMC::read_weight_names(EventStream& info, istringstream& iline)
+ size_t name_size = 0;
+ iline >> name_size;
+ info.weights.names.clear();
+ info.weights.weights.clear();
+
+ string name;
+ WeightContainer namedWeight;
+ string::size_type i1 = line.find("\""), i2, len = line.size();
+ for(size_t ii = 0; ii < name_size; ++ii) {
+ // weight names may contain blanks
+ if(i1 >= len) {
+ cout << "debug: attempting to read past the end of the named weight line " << endl;
+ cout << "debug: We should never get here" << endl;
+ cout << "debug: Looking for the end of this event" << endl;
+ find_event_end(is);
+ }
+ i2 = line.find("\"",i1+1);
+ name = line.substr(i1+1,i2-i1-1);
+ if ( namedWeight.names.find(name) == namedWeight.names.end() )
+ namedWeight.names[name] = namedWeight.names.size();
+ namedWeight.weights[ii] = info.weights.weights[ii];
+ i1 = line.find("\"",i2+1);
+ }
+
+ info.weights.weights = namedWeight.weights;
+ info.weights.names = namedWeight.names;
+#endif
+ return 1;
+}
+
+int HepMC::read_particle(EventStream &info, istringstream& iline, Geant4Particle * p) {
+ float ene = 0., theta = 0., phi = 0;
+ int size = 0;
+
+ // check that the input stream is still OK after reading item
+ iline >> p->id >> p->pdgID >> p->psx >> p->psy >> p->psz >> ene;
+ if ( !iline )
+ return 0;
+ else if( info.io_type != ascii )
+ iline >> p->mass;
+ else
+ p->mass = sqrt(fabs(ene*ene - p->psx*p->psx + p->psy*p->psy + p->psz*p->psz));
+
+ // Reuse here the secondaries to store the end-vertex ID
+ iline >> p->status >> theta >> phi >> p->secondaries >> size;
+ if(!iline)
+ return 0;
+
+ // read flow patterns if any exist
+ for (int i = 0; i < size; ++i ) {
+ iline >> p->colorFlow[0] >> p->colorFlow[1];
+ if(!iline) return 0;
+ }
+ return 1;
+}
+
+int HepMC::read_vertex(EventStream &info, istream& is, istringstream & iline) {
+ int id=0, dummy=0, num_orphans_in=0, num_particles_out=0, weights_size=0;
+ vector<float> weights;
+ Geant4Vertex* v = new Geant4Vertex();
+ Geant4Particle* p;
+
+ iline >> id >> dummy >> v->x >> v->y >> v->z >> v->time
+ >> num_orphans_in >> num_particles_out >> weights_size;
+ if(!iline)
+ return 0;
+ weights.resize(weights_size);
+ for (int i1 = 0; i1 < weights_size; ++i1) {
+ iline >> weights[i1];
+ if(!iline)
+ return 0;
+ }
+ info.add_vertex(id,v);
+ //cout << "Add Vertex:" << id << endl;
+
+ for(char value = is.peek(); value=='P'; value=is.peek()) {
+ value = get_input(is,iline);
+ if( !iline || value < 0 )
+ return 0;
+
+ read_particle(info, iline,p = new Geant4Particle());
+ if(!iline) {
+ cerr << "Failed to read particle!" << endl;
+ delete p;
+ return 0;
+ }
+ info.add_particle(p);
+ p->pex = p->psx;
+ p->pey = p->psy;
+ p->pez = p->psz;
+ if ( --num_orphans_in >= 0 ) {
+ p->vex = v->x;
+ p->vey = v->y;
+ p->vez = v->z;
+ v->in.insert(p->id);
+ //cout << "Add INGOING Particle:" << p->id << endl;
+ }
+ else if ( num_particles_out >= 0 ) {
+ v->out.insert(p->id);
+ p->vsx = v->x;
+ p->vsy = v->y;
+ p->vsz = v->z;
+ //cout << "Add OUTGOING Particle:" << p->id << endl;
+ }
+ else {
+ delete p;
+ throw runtime_error("Invalid number of particles....");
+ }
+ }
+ return 1;
+}
+
+int HepMC::read_event_header(EventStream &info, istringstream & input, EventHeader& header) {
+ // read values into temp variables, then fill GenEvent
+ int random_states_size = 0, nmpi = -1;
+ input >> header.id;
+ if( info.io_type == gen || info.io_type == extascii ) {
+ input >> nmpi;
+ if( input.fail() ) return 0;
+ //MSF set_mpi( nmpi );
+ }
+ input >> header.scale;
+ input >> header.alpha_qcd;
+ input >> header.alpha_qed;
+ input >> header.signal_process_id;
+ input >> header.signal_process_vertex;
+ input >> header.num_vertices;
+ if( info.io_type == gen || info.io_type == extascii )
+ input >> header.bp1 >> header.bp2;
+
+ input >> random_states_size;
+ cout << input.str() << endl;
+ input.clear();
+ if( input.fail() ) return 0;
+
+ header.random.resize(random_states_size);
+ for(int i = 0; i < random_states_size; ++i )
+ input >> header.random[i];
+
+ size_t weights_size = 0;
+ input >> weights_size;
+ if( input.fail() ) return 0;
+
+ vector<float> wgt(weights_size);
+ for(size_t ii = 0; ii < weights_size; ++ii )
+ input >> wgt[ii];
+ if( input.fail() ) return 0;
+
+ // weight names will be added later if they exist
+ if( weights_size > 0 ) header.weights = wgt;
+ return 1;
+}
+
+int HepMC::read_cross_section(EventStream &info, istringstream & input) {
+ input >> info.xsection >> info.xsection_err;
+ return input.fail() ? 0 : 1;
+}
+
+int HepMC::read_units(EventStream &info, istringstream & input) {
+ if( info.io_type == gen ) {
+ string mom, pos;
+ input >> mom >> pos;
+ if ( !input.fail() ) {
+ if ( mom == "KEV" ) info.mom_unit = keV;
+ else if ( mom == "MEV" ) info.mom_unit = MeV;
+ else if ( mom == "GEV" ) info.mom_unit = GeV;
+ else if ( mom == "TEV" ) info.mom_unit = TeV;
+
+ if ( pos == "MM" ) info.pos_unit = mm;
+ else if ( pos == "CM" ) info.pos_unit = cm;
+ else if ( pos == "M" ) info.pos_unit = m;
+ }
+ }
+ return input.fail() ? 0 : 1;
+}
+
+int HepMC::read_heavy_ion(EventStream &, istringstream & input) {
+ // read values into temp variables, then create a new HeavyIon object
+ int nh =0, np =0, nt =0, nc =0,
+ neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
+ float impact = 0., plane = 0., xcen = 0., inel = 0.;
+ input >> nh >> np >> nt >> nc >> neut >> prot >> nw >> nwn >> nwnw;
+ input >> impact >> plane >> xcen >> inel;
+ cerr << "Reading heavy ion, but igoring data!" << endl;
+ /*
+ ion->set_Ncoll_hard(nh);
+ ion->set_Npart_proj(np);
+ ion->set_Npart_targ(nt);
+ ion->set_Ncoll(nc);
+ ion->set_spectator_neutrons(neut);
+ ion->set_spectator_protons(prot);
+ ion->set_N_Nwounded_collisions(nw);
+ ion->set_Nwounded_N_collisions(nwn);
+ ion->set_Nwounded_Nwounded_collisions(nwnw);
+ ion->set_impact_parameter(impact);
+ ion->set_event_plane_angle(plane);
+ ion->set_eccentricity(xcen);
+ ion->set_sigma_inel_NN(inel);
+ */
+ return input.fail() ? 0 : 1;
+}
+
+int HepMC::read_pdf(EventStream &, istringstream & input) {
+ // read values into temp variables, then create a new PdfInfo object
+ int id1 =0, id2 =0, pdf_id1=0, pdf_id2=0;
+ double x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.;
+ input >> id1 ;
+ if ( input.fail() )
+ return 0;
+ // check now for empty PdfInfo line
+ if( id1 == 0 )
+ return 0;
+ // continue reading
+ input >> id2 >> x1 >> x2 >> scale >> pdf1 >> pdf2;
+ if ( input.fail() )
+ return 0;
+ // check to see if we are at the end of the line
+ if( !input.eof() ) {
+ input >> pdf_id1 >> pdf_id2;
+ }
+ cerr << "Reading pdf, but igoring data!" << endl;
+ /*
+ pdf->set_id1( id1 );
+ pdf->set_id2( id2 );
+ pdf->set_pdf_id1( pdf_id1 );
+ pdf->set_pdf_id2( pdf_id2 );
+ pdf->set_x1( x1 );
+ pdf->set_x2( x2 );
+ pdf->set_scalePDF( scale );
+ pdf->set_pdf1( pdf1 );
+ pdf->set_pdf2( pdf2 );
+ */
+ return input.fail() ? 0 : 1;
+}
+
+void HepMC::EventStream::add_particle(Geant4Particle* p) {
+ p->id = m_particles.size();
+ m_particles.insert(make_pair(p->id,p));
+}
+
+Geant4Vertex* HepMC::EventStream::vertex(int i) {
+ Vertices::iterator it=m_vertices.find(i);
+ return (it==m_vertices.end()) ? 0 : (*it).second;
+}
+
+void HepMC::EventStream::clear() {
+ releaseObjects(m_vertices)();
+ releaseObjects(m_particles)();
+}
+
+bool HepMC::EventStream::read() {
+ EventStream& info = *this;
+
+ // OK - now ready to start reading the event, so set the header flag
+ // The flag will be set to false when we reach the end of the header
+ m_particles.clear();
+ m_vertices.clear();
+
+ bool event_read = false;
+ while( instream.good() ) {
+ char value = instream.peek();
+ if ( value == 'E' && event_read ) break;
+ else if ( event_read ) break;
+ else if ( instream.eof() ) return false;
+
+ istringstream input_line;
+ value = get_input(instream,input_line);
+
+ // On failure switch to end
+ if( !input_line || value < 0 )
+ goto Skip;
+
+ switch( value ) {
+ case 'H': {
+ int iotype = 0;
+ string key;
+ input_line >> key;
+ // search for event listing key before first event only.
+ key = key.substr(0,key.find('\r'));
+ if ( key == "H" && (this->io_type == gen || this->io_type == extascii) ) {
+ read_heavy_ion(info, input_line);
+ break;
+ }
+ else if( key == "HepMC::IO_GenEvent-START_EVENT_LISTING" )
+ this->set_io(gen,key);
+ else if( key == "HepMC::IO_Ascii-START_EVENT_LISTING" )
+ this->set_io(ascii,key);
+ else if( key == "HepMC::IO_ExtendedAscii-START_EVENT_LISTING" )
+ this->set_io(extascii,key);
+ else if( key == "HepMC::IO_Ascii-START_PARTICLE_DATA" )
+ this->set_io(ascii_pdt,key);
+ else if( key == "HepMC::IO_ExtendedAscii-START_PARTICLE_DATA" )
+ this->set_io(extascii_pdt,key);
+ else if( key == "HepMC::IO_GenEvent-END_EVENT_LISTING" )
+ iotype = gen;
+ else if( key == "HepMC::IO_Ascii-END_EVENT_LISTING" )
+ iotype = ascii;
+ else if( key == "HepMC::IO_ExtendedAscii-END_EVENT_LISTING" )
+ iotype = extascii;
+ else if( key == "HepMC::IO_Ascii-END_PARTICLE_DATA" )
+ iotype = ascii_pdt;
+ else if( key == "HepMC::IO_ExtendedAscii-END_PARTICLE_DATA" )
+ iotype = extascii_pdt;
+
+ if( iotype != 0 && this->io_type != iotype ) {
+ cerr << "GenEvent::find_end_key: iotype keys have changed. "
+ << "MALFORMED INPUT" << endl;
+ instream.clear(ios::badbit);
+ return false;
+ }
+ continue;
+ }
+ case 'E': // deal with the event line
+ if ( !read_event_header(info, input_line, this->header) )
+ goto Skip;
+ event_read = true;
+ continue;
+
+ case 'N': // get weight names
+ if ( !read_weight_names(info, input_line) )
+ goto Skip;
+ continue;
+
+ case 'U': // get unit information if it exists
+ if ( !read_units(info, input_line) )
+ goto Skip;
+ continue;
+
+ case 'C': // we have a GenCrossSection line
+ if ( !read_cross_section(info, input_line) )
+ goto Skip;
+ continue;
+
+ case 'V': // Read vertex with particles
+ if ( !read_vertex(info, instream, input_line) )
+ goto Skip;
+ continue;
+
+ case 'F': // Read PDF
+ if ( !read_pdf(info, input_line) )
+ goto Skip;
+ continue;
+
+ case 'P': // we should not find this line
+ cerr << "streaming input: found unexpected Particle line." << endl;
+ continue;
+
+ default: // ignore everything else
+ continue;
+ }
+ continue;
+ Skip:
+ printout(WARNING,"HepMC::EventStream","+++ Skip event with ID: %d",this->header.id);
+ releaseObjects(vertices())();
+ releaseObjects(particles())();
+ find_event_end(instream);
+ event_read = false;
+ if ( instream.eof() ) return false;
+ }
+ this->fix_particles();
+ releaseObjects(vertices())();
+ return true;
+}
diff --git a/DDG4/plugins/Geant4SensDet.cpp b/DDG4/plugins/Geant4SensDet.cpp
index 6e1456cd3..b91f0de30 100644
--- a/DDG4/plugins/Geant4SensDet.cpp
+++ b/DDG4/plugins/Geant4SensDet.cpp
@@ -64,6 +64,9 @@ namespace DD4hep {
virtual public Geant4ActionSD,
virtual public RefCountedSequence<Geant4SensDetActionSequence>
{
+ protected:
+ /// Access to the geant4 sensitive detector handle
+ SensitiveDetector m_sensitive;
public:
/// Constructor. The detector element is identified by the name
Geant4SensDet(const std::string& nam, Geometry::LCDD& lcdd)
@@ -71,6 +74,7 @@ namespace DD4hep {
Geant4Action(0,nam), Geant4ActionSD(nam), Base()
{
Geant4Kernel& kernel = Geant4Kernel::access(lcdd);
+ m_sensitive = lcdd.sensitiveDetector(nam);
m_context = Geant4Context(&kernel);
m_outputLevel = kernel.getOutputLevel(nam);
_aquire(kernel.sensitiveAction(nam));
@@ -98,6 +102,12 @@ namespace DD4hep {
/// Access to the readout geometry of the sensitive detector
virtual G4VReadOutGeometry* readoutGeometry() const
{ return this->G4VSensitiveDetector::GetROgeometry(); }
+ /// Access to the LCDD sensitive detector handle
+ virtual SensitiveDetector sensitiveDetector() const
+ { return m_sensitive; }
+ /// Access to the sensitive type of the detector
+ virtual const std::string& sensitiveType() const
+ { return m_sequence->sensitiveType(); }
/// Callback if the sequence should be accepted or filtered off
virtual G4bool Accept(const G4Step* step) const
{ return m_sequence->accept(step); }
diff --git a/DDG4/python/DDG4Dict.C b/DDG4/python/DDG4Dict.C
index 0221e930f..22c665d51 100644
--- a/DDG4/python/DDG4Dict.C
+++ b/DDG4/python/DDG4Dict.C
@@ -39,6 +39,7 @@ using namespace DD4hep::Simulation;
#include "DDG4/Geant4Config.h"
#include "DDG4/Geant4DataDump.h"
+#include "DDG4/Geant4InputAction.h"
#include <iostream>
namespace DD4hep {
@@ -188,6 +189,8 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation;
// CINT configuration for DDG4
#if defined(__MAKECINT__)
#pragma link C++ class PropertyResult;
+#pragma link C++ class Geant4InputAction::Particles;
+#pragma link C++ class auto_ptr<Geant4InputAction::Particles>;
#pragma link C++ class ActionHandle;
#pragma link C++ class FilterHandle;
@@ -234,6 +237,8 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation;
#pragma link C++ class Geant4GeneratorActionSequence;
#pragma link C++ class Geant4GeneratorAction;
+#pragma link C++ class Geant4InputAction;
+#pragma link C++ class Geant4EventReader;
#pragma link C++ class Geant4PhysicsListActionSequence;
#pragma link C++ class Geant4PhysicsList;
diff --git a/DDG4/src/Geant4Data.cpp b/DDG4/src/Geant4Data.cpp
index c255e06dc..240f39951 100644
--- a/DDG4/src/Geant4Data.cpp
+++ b/DDG4/src/Geant4Data.cpp
@@ -66,7 +66,7 @@ Geant4HitData::Contribution Geant4HitData::extractContribution(G4Step* step) {
(h.trackDef() == G4OpticalPhoton::OpticalPhotonDefinition()) ? h.trkEnergy() : h.totalEnergy();
const G4ThreeVector& pre = h.prePosG4();
const G4ThreeVector& post = h.postPosG4();
- float pos[] = { (pre.x()+post.x())/2.0, (pre.y()+post.y())/2.0, (pre.z()+post.z())/2.0 };
+ float pos[] = {float((pre.x()+post.x())/2.0),float((pre.y()+post.y())/2.0),float((pre.z()+post.z())/2.0) };
Contribution contrib(h.trkID(),h.trkPdgID(),deposit,h.trkTime(),pos);
return contrib;
}
diff --git a/DDG4/src/Geant4DetectorConstruction.cpp b/DDG4/src/Geant4DetectorConstruction.cpp
index 80d847ea6..15ee8bfdb 100644
--- a/DDG4/src/Geant4DetectorConstruction.cpp
+++ b/DDG4/src/Geant4DetectorConstruction.cpp
@@ -21,25 +21,42 @@
using namespace std;
using namespace DD4hep;
-using namespace DD4hep::Geometry;
+using namespace DD4hep::Simulation;
-DD4hep::Simulation::Geant4DetectorConstruction::Geant4DetectorConstruction(Geometry::LCDD& lcdd)
- : m_outputLevel(PrintLevel(printLevel()-1)), m_lcdd(lcdd), m_world(0) {
+namespace {
+ static Geant4DetectorConstruction* s_instance = 0;
}
-DD4hep::Simulation::Geant4DetectorConstruction::Geant4DetectorConstruction(DD4hep::Simulation::Geant4Kernel& kernel)
- : m_lcdd(kernel.lcdd()), m_world(0) {
+/// Instance accessor
+Geant4DetectorConstruction* Geant4DetectorConstruction::instance(Geant4Kernel& kernel) {
+ if ( 0 == s_instance ) s_instance = new Geant4DetectorConstruction(kernel);
+ return s_instance;
+}
+
+/// Initializing constructor for other clients
+Geant4DetectorConstruction::Geant4DetectorConstruction(Geometry::LCDD& lcdd)
+ : Geant4Action(0,"DetectorConstruction"), m_lcdd(lcdd), m_world(0)
+{
+ m_outputLevel = PrintLevel(printLevel()-1);
+ s_instance = this;
+}
+
+/// Initializing constructor for DDG4
+Geant4DetectorConstruction::Geant4DetectorConstruction(Geant4Kernel& kernel)
+ : Geant4Action(0,"DetectorConstruction"), m_lcdd(kernel.lcdd()), m_world(0)
+{
m_outputLevel = kernel.getOutputLevel("Geant4Converter");
+ s_instance = this;
}
/// Default destructor
-DD4hep::Simulation::Geant4DetectorConstruction::~Geant4DetectorConstruction() {
+Geant4DetectorConstruction::~Geant4DetectorConstruction() {
+ s_instance = 0;
}
-G4VPhysicalVolume* DD4hep::Simulation::Geant4DetectorConstruction::Construct() {
- typedef Simulation::Geant4Converter Geant4Converter;
+G4VPhysicalVolume* Geant4DetectorConstruction::Construct() {
Geant4Mapping& g4map = Geant4Mapping::instance();
- DetElement world = m_lcdd.world();
+ Geometry::DetElement world = m_lcdd.world();
Geant4Converter conv(m_lcdd, m_outputLevel);
Geant4GeometryInfo* info = conv.create(world).detach();
g4map.attach(info);
@@ -57,3 +74,66 @@ G4VPhysicalVolume* DD4hep::Simulation::Geant4DetectorConstruction::Construct() {
#endif
return m_world;
}
+
+/// Access to the converted regions
+const Geant4GeometryMaps::RegionMap& Geant4DetectorConstruction::regions() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Regions;
+ throw runtime_error("+++ Geant4DetectorConstruction::regions: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted sensitive detectors
+const Geant4GeometryMaps::SensDetMap& Geant4DetectorConstruction::sensitives() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4SensDets;
+ throw runtime_error("+++ Geant4DetectorConstruction::sensitives: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted volumes
+const Geant4GeometryMaps::VolumeMap& Geant4DetectorConstruction::volumes() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Volumes;
+ throw runtime_error("+++ Geant4DetectorConstruction::volumes: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted shapes
+const Geant4GeometryMaps::SolidMap& Geant4DetectorConstruction::shapes() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Solids;
+ throw runtime_error("+++ Geant4DetectorConstruction::shapes: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted limit sets
+const Geant4GeometryMaps::LimitMap& Geant4DetectorConstruction::limits() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Limits;
+ throw runtime_error("+++ Geant4DetectorConstruction::limits: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted assemblies
+const Geant4GeometryMaps::AssemblyMap& Geant4DetectorConstruction::assemblies() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4AssemblyVolumes;
+ throw runtime_error("+++ Geant4DetectorConstruction::assemblies: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted placements
+const Geant4GeometryMaps::PlacementMap& Geant4DetectorConstruction::placements() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Placements;
+ throw runtime_error("+++ Geant4DetectorConstruction::placements: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted materials
+const Geant4GeometryMaps::MaterialMap& Geant4DetectorConstruction::materials() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Materials;
+ throw runtime_error("+++ Geant4DetectorConstruction::materials: Access not possible. Geometry is not yet converted!");
+}
+
+/// Access to the converted elements
+const Geant4GeometryMaps::ElementMap& Geant4DetectorConstruction::elements() const {
+ Geant4GeometryInfo* p = Geant4Mapping::instance().ptr();
+ if ( p ) return p->g4Elements;
+ throw runtime_error("+++ Geant4DetectorConstruction::elements: Access not possible. Geometry is not yet converted!");
+}
diff --git a/DDG4/src/Geant4Exec.cpp b/DDG4/src/Geant4Exec.cpp
index 78a9198a3..60bc26658 100644
--- a/DDG4/src/Geant4Exec.cpp
+++ b/DDG4/src/Geant4Exec.cpp
@@ -280,8 +280,6 @@ namespace DD4hep {
runAction->releaseContextFromClients();
destroyClientContext(evt);
}
-
-
}
}
@@ -315,7 +313,7 @@ int Geant4Exec::configure(Geant4Kernel& kernel) {
"You sure you loaded the geometry properly?",int(lcdd.detectors().size()));
}
// Get the detector constructed
- Geant4DetectorConstruction* detector = new Geant4DetectorConstruction(kernel);
+ Geant4DetectorConstruction* detector = Geant4DetectorConstruction::instance(kernel);
runManager.SetUserInitialization(detector);
G4VUserPhysicsList* physics = 0;
diff --git a/DDG4/src/Geant4InputAction.cpp b/DDG4/src/Geant4InputAction.cpp
index f833caa6c..161a4e9de 100644
--- a/DDG4/src/Geant4InputAction.cpp
+++ b/DDG4/src/Geant4InputAction.cpp
@@ -9,11 +9,10 @@
//====================================================================
// Framework include files
-#include "DD4hep/Printout.h"
-#include "DD4hep/Primitives.h"
-#include "DDG4/Geant4InputAction.h"
+#include "DD4hep/Plugins.h"
#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4InputAction.h"
#include "G4Event.hh"
@@ -23,7 +22,7 @@ typedef DD4hep::ReferenceBitMask<int> PropertyMask;
/// Initializing constructor
Geant4EventReader::Geant4EventReader(const std::string& nam)
- : m_name(nam), m_directAccess(false)
+: m_name(nam), m_directAccess(false), m_currEvent(0)
{
}
@@ -31,6 +30,41 @@ Geant4EventReader::Geant4EventReader(const std::string& nam)
Geant4EventReader::~Geant4EventReader() {
}
+/// Skip event. To be implemented for sequential sources
+Geant4EventReader::EventReaderStatus Geant4EventReader::skipEvent() {
+ if ( hasDirectAccess() ) {
+ ++m_currEvent;
+ return EVENT_READER_OK;
+ }
+ std::vector<Particle*> particles;
+ ++m_currEvent;
+ EventReaderStatus sc = readParticles(m_currEvent,particles);
+ for_each(particles.begin(),particles.end(),deleteObject<Particle>);
+ return sc;
+}
+
+/// Move to the indicated event number.
+Geant4EventReader::EventReaderStatus
+Geant4EventReader::moveToEvent(int event_number) {
+ if ( m_currEvent == event_number-1 ) {
+ return EVENT_READER_OK;
+ }
+ else if ( hasDirectAccess() ) {
+ m_currEvent = event_number-1;
+ return EVENT_READER_OK;
+ }
+ else if ( event_number<m_currEvent ) {
+ return EVENT_READER_ERROR;
+ }
+ else {
+ for(int i=m_currEvent; i<event_number;++i)
+ skipEvent();
+ m_currEvent = event_number-1;
+ return EVENT_READER_OK;
+ }
+ return EVENT_READER_ERROR;
+}
+
/// Standard constructor
Geant4InputAction::Geant4InputAction(Geant4Context* context, const string& name)
: Geant4GeneratorAction(context,name), m_reader(0)
@@ -65,22 +99,29 @@ int Geant4InputAction::readParticles(int evt_number, std::vector<Particle*>& par
try {
m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
if ( 0 == m_reader ) {
+ PluginDebug dbg();
+ m_reader = PluginService::Create<Geant4EventReader*>(tn.first,tn.second);
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;
+ return Geant4EventReader::EVENT_READER_NO_FACTORY;
}
}
catch(const exception& e) {
err = e.what();
}
if ( !err.empty() ) {
- abortRun(issue(evid)+err,"Error when reading file %s",m_input.c_str());
- return 0;
+ abortRun(issue(evid)+err,"Error when creating reader for file %s",m_input.c_str());
+ return Geant4EventReader::EVENT_READER_NO_FACTORY;
}
}
- int status = m_reader->readParticles(evid,particles);
- if ( 0 != status ) {
+ int status = m_reader->moveToEvent(evid);
+ if ( Geant4EventReader::EVENT_READER_OK != status ) {
+ abortRun(issue(evid)+"Error when moving to event - may be end of file.",
+ "Error when reading file %s",m_input.c_str());
+ }
+ status = m_reader->readParticles(evid,particles);
+ if ( Geant4EventReader::EVENT_READER_OK != status ) {
abortRun(issue(evid)+"Error when reading file - may be end of file.",
"Error when reading file %s",m_input.c_str());
}
@@ -95,7 +136,7 @@ void Geant4InputAction::operator()(G4Event* event) {
Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
int result = readParticles(event->GetEventID(),primaries);
- if ( result != 0 ) { // handle I/O error, but how?
+ if ( result != Geant4EventReader::EVENT_READER_OK ) { // handle I/O error, but how?
return;
}
prim->add(m_mask, inter);
diff --git a/DDG4/src/Geant4InteractionVertexSmear.cpp b/DDG4/src/Geant4InteractionVertexSmear.cpp
index a531b6490..5e33448d1 100644
--- a/DDG4/src/Geant4InteractionVertexSmear.cpp
+++ b/DDG4/src/Geant4InteractionVertexSmear.cpp
@@ -39,8 +39,6 @@ Geant4InteractionVertexSmear::~Geant4InteractionVertexSmear() {
/// 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);
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index 90a700d9c..c2b3155d7 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -13,6 +13,7 @@
#include "DDG4/Geant4StepHandler.h"
#include "DDG4/Geant4TrackHandler.h"
#include "DDG4/Geant4EventAction.h"
+#include "DDG4/Geant4SensDetAction.h"
#include "DDG4/Geant4TrackingAction.h"
#include "DDG4/Geant4SteppingAction.h"
#include "DDG4/Geant4ParticleHandler.h"
@@ -148,14 +149,18 @@ void Geant4ParticleHandler::mark(const G4Track* track) {
PropertyMask mask(m_currTrack.reason);
mask.set(G4PARTICLE_CREATED_HIT);
/// Check if the track origines from the calorimeter.
- // If yes, flag it, because it is a candidate fro removal.
- G4VPhysicalVolume* vol = track->GetVolume();
- if ( strstr(vol->GetName().c_str(),"cal") ) { // just for test!
+ // If yes, flag it, because it is a candidate for removal.
+ G4LogicalVolume* vol = track->GetVolume()->GetLogicalVolume();
+ G4VSensitiveDetector* g4 = vol->GetSensitiveDetector();
+ Geant4ActionSD* sd = dynamic_cast<Geant4ActionSD*>(g4);
+ string typ = sd ? sd->sensitiveType() : string();
+ if ( typ == "calorimeter" )
mask.set(G4PARTICLE_CREATED_CALORIMETER_HIT);
- }
- else if ( !mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT) ) {
+ else if ( typ == "tracker" )
mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
- }
+ else // Assume by default "tracker"
+ mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
+
//Geant4ParticleHandle(&m_currTrack).dump4(outputLevel(),vol->GetName(),"hit created by particle");
}
diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp
index 46a3918ae..9ebcd47ff 100644
--- a/DDG4/src/Geant4SensDetAction.cpp
+++ b/DDG4/src/Geant4SensDetAction.cpp
@@ -197,7 +197,6 @@ long long int Geant4Sensitive::volumeID(G4Step* s) {
/// Returns the cellID(volumeID+local coordinate encoding) of the sensitive volume corresponding to the step
long long int Geant4Sensitive::cellID(G4Step* s) {
StepHandler h(s);
- typedef DDSegmentation::Vector3D _V;
Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager();
VolumeID volID = volMgr.volumeID(h.preTouchable());
if ( m_segmentation.isValid() ) {
@@ -219,6 +218,7 @@ Geant4SensDetActionSequence::Geant4SensDetActionSequence(Geant4Context* context,
context->sensitiveActions().insert(name(), this);
/// Update the sensitive detector type, so that the proper instance is created
m_sensitive = context->lcdd().sensitiveDetector(nam);
+ m_sensitiveType = m_sensitive.type();
m_sensitive.setType("Geant4SensDet");
InstanceCount::increment(this);
}
diff --git a/UtilityApps/src/teve_display.cpp b/UtilityApps/src/teve_display.cpp
index 4433192ca..0460b6e09 100644
--- a/UtilityApps/src/teve_display.cpp
+++ b/UtilityApps/src/teve_display.cpp
@@ -75,7 +75,7 @@ static long teve_display(LCDD& lcdd, int /* argc */, char** /* argv */) {
TEveGeoTopNode* tn = new TEveGeoTopNode(mgr, mgr->GetTopNode());
tn->SetVisLevel(4);
- EvNavHandler *fh = new EvNavHandler;
+ /* EvNavHandler *fh = */ new EvNavHandler;
// // ---- try to set transparency - does not seem to work ...
// TGeoNode* node1 = gGeoManager->GetTopNode();
diff --git a/cmake/DD4hepConfig.cmake.in b/cmake/DD4hepConfig.cmake.in
index 99ae5065c..3ea3dde09 100644
--- a/cmake/DD4hepConfig.cmake.in
+++ b/cmake/DD4hepConfig.cmake.in
@@ -50,7 +50,7 @@ INCLUDE( ${DD4hep_ROOT}/cmake/DD4hepMacros.cmake )
# additional components are set by cmake in variable PKG_FIND_COMPONENTS
# first argument should be the package name
if(@DD4HEP_USE_GEANT4@)
- CHECK_PACKAGE_LIBS( DD4hep DD4hepCore DDSegmentation DD4hepRec DD4hepPlugins DD4hepG4 )
+ CHECK_PACKAGE_LIBS( DD4hep DD4hepCore DDSegmentation DD4hepRec DD4hepG4 )
#--- if geant 4 was built with CLHEP we need to export this to client packages
if( @GEANT4_USE_CLHEP@)
@@ -60,7 +60,7 @@ if(@DD4HEP_USE_GEANT4@)
set( Geant4_DIR @Geant4_DIR@)
else()
- CHECK_PACKAGE_LIBS( DD4hep DD4hepCore DDSegmentation DD4hepRec DD4hepPlugins )
+ CHECK_PACKAGE_LIBS( DD4hep DD4hepCore DDSegmentation DD4hepRec )
endif()
--
GitLab