diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index 37b259bc86dcf506a6a5780b7f2fd625df1cb526..ef58298f4d7b04967d47ccb0586a9e7018a0da38 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -43,8 +43,8 @@ def run():
# Configure I/O
evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
- evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
- evt_lcio.OutputLevel = Output.ERROR
+ #evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M'))
+ #evt_lcio.OutputLevel = Output.ERROR
generator_output_level = Output.INFO
@@ -57,7 +57,7 @@ def run():
Generation of isotrope tracks using the DDG4 partcle gun:
# Setup particle gun
- gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=3)
+ gun = simple.setupGun('Gun','pi-',energy=10*GeV,isotrop=True,multiplicity=1)
gun.OutputLevel = 5 # generator_output_level
gun.Mask = 1
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -101,16 +101,15 @@ def run():
"""
Generation of primary particles from LCIO input files
"""
-
# First particle file reader
gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO1");
- gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_1.stdhep"
+ #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_1.stdhep"
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/qq_gen_128_999.stdhep"
- gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/smuonLR_PointK_3TeV_BS_noBkg_run0001.stdhep"
+ #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/smuonLR_PointK_3TeV_BS_noBkg_run0001.stdhep"
gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbbb_3TeV.stdhep"
+ #gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_pi-_5GeV.slcio"
gen.OutputLevel = 4 # generator_output_level
gen.Mask = 1
- gen.MomentumScale = 0.1
kernel.generatorAction().adopt(gen)
# Install vertex smearing for this interaction
gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/Smear1");
@@ -135,6 +134,7 @@ def run():
gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns)
gen.Sigma = (2*mm, 1*mm, 1*mm, 0*ns)
kernel.generatorAction().adopt(gen)
+
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
# Merge all existing interaction records
@@ -152,7 +152,10 @@ def run():
kernel.generatorAction().adopt(part)
#part.SaveProcesses = ['conv','Decay']
part.SaveProcesses = ['Decay']
- part.OutputLevel = 4 # generator_output_level
+ part.MinimalKineticEnergy = 100*MeV
+ part.Z_trackers = 165.70*cm # EcalEndcap_zmin
+ part.R_trackers = 126.50*cm # EcalBarrel_rmin
+ part.OutputLevel = 5 # generator_output_level
part.enableUI()
#user = DDG4.Action(kernel,"Geant4UserParticleHandler/UserParticleHandler")
#part.adopt(user)
diff --git a/DDG4/include/DDG4/Geant4HitCollection.h b/DDG4/include/DDG4/Geant4HitCollection.h
index 3978dbd59db56880c656275282e633ba61a9155c..0a642b69cea04ae1a85627e6f6cb23bab6eb8349 100644
--- a/DDG4/include/DDG4/Geant4HitCollection.h
+++ b/DDG4/include/DDG4/Geant4HitCollection.h
@@ -288,6 +288,51 @@ namespace DD4hep {
void getHitsUnchecked(std::vector<void*>& result);
};
+
+ /** @class PositionCompare Geant4HitCollection.h DDG4/Geant4HitCollection.h
+ *
+ * Class for hit matching using the hit position.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ template<typename TYPE, typename POS> struct PositionCompare : public Geant4HitCollection::Compare {
+ const POS& pos;
+ /// Constructor
+ PositionCompare(const POS& p) : pos(p) { }
+ /// Comparison function
+ virtual void* operator()(const Geant4HitWrapper& w) const;
+ };
+
+ template <typename TYPE, typename POS>
+ void* PositionCompare<TYPE,POS>::operator()(const Geant4HitWrapper& w) const {
+ TYPE* h = w;
+ return pos == h->position ? h : 0;
+ }
+
+ /** @class PositionCompare Geant4HitCollection.h DDG4/Geant4HitCollection.h
+ *
+ * Class for hit matching using the hit's cell identifier.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ template<typename TYPE> struct CellIDCompare : public Geant4HitCollection::Compare {
+ long long int id;
+ /// Constructor
+ CellIDCompare(long long int i) : id(i) { }
+ /// Comparison function.
+ virtual void* operator()(const Geant4HitWrapper& w) const;
+ };
+
+ template <typename TYPE>
+ void* CellIDCompare<TYPE>::operator()(const Geant4HitWrapper& w) const {
+ TYPE* h = w;
+ if ( id == h->cellID )
+ return h;
+ return 0;
+ }
+
} // End namespace Simulation
} // End namespace DD4hep
diff --git a/DDG4/include/DDG4/Geant4InputHandling.h b/DDG4/include/DDG4/Geant4InputHandling.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6779249f7d4cf10c3307cd8c121676cc13f3289
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4InputHandling.h
@@ -0,0 +1,55 @@
+//====================================================================
+// AIDA Detector description implementation
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+#ifndef DD4HEP_DDG4_GEANT4INPUTHANDLING_H
+#define DD4HEP_DDG4_GEANT4INPUTHANDLING_H
+
+// Framework include files
+#include "DDG4/Geant4Primary.h"
+
+// Forward declarations
+class G4Event;
+
+/*
+ * DD4hep namespace declaration
+ */
+namespace DD4hep {
+
+ /*
+ * Simulation namespace declaration
+ */
+ namespace Simulation {
+
+ // Forward declarations
+ class Geant4Action;
+ class Geant4Context;
+
+ /// Initialize the generation of one event
+ int generationInitialization(const Geant4Action* caller,const Geant4Context* context);
+
+ /// Merge all interactions present in the context
+ int mergeInteractions(const Geant4Action* caller,const Geant4Context* context);
+
+ /// Boost particles of one interaction identified by its mask
+ int boostInteraction(const Geant4Action* caller,
+ Geant4PrimaryEvent::Interaction* inter,
+ double alpha);
+
+ /// Smear the primary vertex of an interaction
+ int smearInteraction(const Geant4Action* caller,
+ Geant4PrimaryEvent::Interaction* inter,
+ double dx, double dy, double dz, double dt);
+
+
+ /// Generate all primary vertices corresponding to the merged interaction
+ int generatePrimaries(const Geant4Action* caller,
+ const Geant4Context* context,
+ G4Event* event);
+
+ } // End namespace Simulation
+} // End namespace DD4hep
+#endif /* DD4HEP_DDG4_GEANT4INPUTHANDLING_H */
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index f79faa0e283c2d9c7f448410048b67e94733900b..0dc6a874d41a03f3b5523e087d838dd4edf79581 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -191,6 +191,7 @@ namespace DD4hep {
void dump2(int level, const std::string& src, const char* tag, int g4id, bool inrec) const;
/// Output type 3:+++ <tag> ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] #Par: 0 #Dau: 4
void dump3(int level, const std::string& src, const char* tag) const;
+ void dump4(int level, const std::string& src, const char* tag) const;
/// Handlers
diff --git a/DDG4/include/DDG4/Geant4ParticleGun.h b/DDG4/include/DDG4/Geant4ParticleGun.h
index 253ccca46c1db45fd160b0ff20e4c254a71995fb..17750f03bc875323ac2338054e807249279c1111 100644
--- a/DDG4/include/DDG4/Geant4ParticleGun.h
+++ b/DDG4/include/DDG4/Geant4ParticleGun.h
@@ -24,6 +24,7 @@ namespace DD4hep {
* Simulation namespace declaration
*/
namespace Simulation {
+
/** @class Geant4ParticleGun Geant4ParticleGun.h DDG4/Geant4ParticleGun.h
*
* Implementation wrapper of the Geant4 particle gun
@@ -45,6 +46,8 @@ namespace DD4hep {
int m_mask;
/// Property: Isotrope particles?
bool m_isotrop;
+ /// Property: Standalone mode: includes interaction merging and primary generation
+ bool m_standalone;
/// Pointer to geant4 particle definition
G4ParticleDefinition* m_particle;
diff --git a/DDG4/include/DDG4/Geant4ParticleHandler.h b/DDG4/include/DDG4/Geant4ParticleHandler.h
index e3f390f224445436f89067efcbeb706ee3e15e27..ff291107d6c393222167003f101980886468cfe2 100644
--- a/DDG4/include/DDG4/Geant4ParticleHandler.h
+++ b/DDG4/include/DDG4/Geant4ParticleHandler.h
@@ -75,6 +75,9 @@ namespace DD4hep {
/// Property: Energy cut below which particles are not collected, but assigned to the parent
double m_kinEnergyCut;
+ double m_zTracker, m_rTracker;
+
+
/// Global particle identifier. Obtained at the begin of the event.
int m_globalParticleID, m_initialParticleID;
/// User action pointer
diff --git a/DDG4/include/DDG4/Geant4ParticlePrint.h b/DDG4/include/DDG4/Geant4ParticlePrint.h
index a9a44a697bf5f5fe84d6a5d83e577a55545f3543..808ef2e5d1f69a3e2a34d0d6d4a6899a6de2c80a 100644
--- a/DDG4/include/DDG4/Geant4ParticlePrint.h
+++ b/DDG4/include/DDG4/Geant4ParticlePrint.h
@@ -14,6 +14,9 @@
#include "DDG4/Geant4GeneratorAction.h"
#include "DDG4/Geant4Particle.h"
+// Forward declarations
+class G4Event;
+
/*
* DD4hep namespace declaration
*/
@@ -44,16 +47,18 @@ namespace DD4hep {
bool m_printEnd;
/// Property: Flag to indicate output type as part of the generator action
bool m_printGeneration;
+ /// Property: Flag to indicate output of hit data in tree
+ bool m_printHits;
- void printParticle(const std::string& prefix, Geant4ParticleHandle p) const;
+ void printParticle(const std::string& prefix, const G4Event* e, Geant4ParticleHandle p) const;
/// Print record of kept particles
- void printParticles(const ParticleMap& particles) const;
+ void printParticles(const G4Event* e, const ParticleMap& particles) const;
/// Print tree of kept particles
- void printParticleTree(const ParticleMap& particles, int level, Geant4ParticleHandle p) const;
+ void printParticleTree(const G4Event* e, const ParticleMap& particles, int level, Geant4ParticleHandle p) const;
/// Print tree of kept particles
- void printParticleTree(const ParticleMap& particles) const;
+ void printParticleTree(const G4Event* e, const ParticleMap& particles) const;
/// Print particle table
- void makePrintout(int event_id) const;
+ void makePrintout(const G4Event* e) const;
public:
diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h
index 2ec461117ac90d5e18eebbdcbdcb6f1ac6509a62..86107bbc0312c1069f32540dc190f8f42d9f4628 100644
--- a/DDG4/include/DDG4/Geant4SensDetAction.h
+++ b/DDG4/include/DDG4/Geant4SensDetAction.h
@@ -30,6 +30,12 @@ class G4VReadOutGeometry;
*/
namespace DD4hep {
+ // Forward declarations
+ namespace Geometry {
+ class LCDD;
+ class DetElement;
+ }
+
/*
* Simulation namespace declaration
*/
@@ -405,6 +411,44 @@ namespace DD4hep {
return sequence().defineCollection<TYPE>(this, coll_name);
}
+ /** @class Geant4SensitiveAction Geant4SensDetAction.h DDG4/Geant4SensDetAction.h
+ *
+ * Templated implementation to realize sensitive detectors.
+ * Default implementations for all functions are provided in the file
+ * DDG4/Geant4SensDetAction.inl.
+ *
+ * Users may override any of the templated callbacks or the of the virtual functions
+ * of the base class using explicit template specialization.
+ * An example may be found in DDG4/plugins/eant4SDActions.
+ *
+ * @author M.Frank
+ * @version 1.0
+ */
+ template <typename T> class Geant4SensitiveAction : public Geant4Sensitive {
+ protected:
+ /// Collection identifier
+ size_t m_collectionID;
+ public:
+ /// Standard , initializing constructor
+ Geant4SensitiveAction(Geant4Context* context,
+ const std::string& name,
+ Geometry::DetElement det,
+ Geometry::LCDD& lcdd);
+ /// Default destructor
+ virtual ~Geant4SensitiveAction();
+ /// Define collections created by this sensitivie action object
+ virtual void defineCollections() {}
+ /// G4VSensitiveDetector interface: Method invoked at the begining of each event.
+ virtual void begin(G4HCofThisEvent* hce);
+ /// G4VSensitiveDetector interface: Method invoked at the end of each event.
+ virtual void end(G4HCofThisEvent* hce);
+ /// G4VSensitiveDetector interface: Method for generating hit(s) using the G4Step object.
+ virtual bool process(G4Step* step,G4TouchableHistory* history);
+ /// G4VSensitiveDetector interface: Method invoked if the event was aborted.
+ virtual void clear(G4HCofThisEvent* hce);
+ };
+
+
} // End namespace Simulation
} // End namespace DD4hep
diff --git a/DDG4/include/DDG4/Geant4SensDetAction.inl b/DDG4/include/DDG4/Geant4SensDetAction.inl
new file mode 100644
index 0000000000000000000000000000000000000000..cd63b2f0bfa0584654f2b7582abe8766336ee227
--- /dev/null
+++ b/DDG4/include/DDG4/Geant4SensDetAction.inl
@@ -0,0 +1,57 @@
+// $Id: Geant4Field.cpp 513 2013-04-05 14:31:53Z gaede $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+
+// Framework include files
+#ifndef DD4HEP_DDG4_GEANT4SENSDETACTION_H
+#include "DDG4/Geant4SensDetAction.h"
+#endif
+#include "DD4hep/InstanceCount.h"
+
+namespace DD4hep {
+ namespace Simulation {
+
+ /// Standard , initializing constructor
+ template <typename T>
+ Geant4SensitiveAction<T>::Geant4SensitiveAction(Geant4Context* context,
+ const std::string& name,
+ Geometry::DetElement det,
+ Geometry::LCDD& lcdd)
+ : Geant4Sensitive(context,name,det,lcdd), m_collectionID(0)
+ {
+ defineCollections();
+ InstanceCount::increment(this);
+ }
+
+ /// Default destructor
+ template <typename T> Geant4SensitiveAction<T>::~Geant4SensitiveAction() {
+ InstanceCount::decrement(this);
+ }
+
+ /// G4VSensitiveDetector interface: Method invoked at the begining of each event.
+ template <typename T> void Geant4SensitiveAction<T>::begin(G4HCofThisEvent* hce) {
+ Geant4Sensitive::begin(hce);
+ }
+
+ /// G4VSensitiveDetector interface: Method invoked at the end of each event.
+ template <typename T> void Geant4SensitiveAction<T>::end(G4HCofThisEvent* hce) {
+ Geant4Sensitive::end(hce);
+ }
+
+ /// G4VSensitiveDetector interface: Method for generating hit(s) using the G4Step object.
+ template <typename T> bool Geant4SensitiveAction<T>::process(G4Step* step,G4TouchableHistory* history) {
+ return Geant4Sensitive::process(step,history);
+ }
+
+ /// G4VSensitiveDetector interface: Method invoked if the event was aborted.
+ template <typename T> void Geant4SensitiveAction<T>::clear(G4HCofThisEvent* hce) {
+ Geant4Sensitive::clear(hce);
+ }
+
+ }
+}
diff --git a/DDG4/include/DDG4/Geant4TrackHandler.h b/DDG4/include/DDG4/Geant4TrackHandler.h
index c5553631c7b095931a60112c30592627233631b5..35e93a04b059831ea0ce3ee94e185a102bdb0b13 100644
--- a/DDG4/include/DDG4/Geant4TrackHandler.h
+++ b/DDG4/include/DDG4/Geant4TrackHandler.h
@@ -137,6 +137,9 @@ namespace DD4hep {
G4VPhysicalVolume* vol() const {
return track->GetVolume();
}
+ G4ThreeVector momentum() const {
+ return track->GetMomentum();
+ }
/// Next physical volume of the track
G4VPhysicalVolume* nextVol() const {
return track->GetNextVolume();
diff --git a/DDG4/plugins/Geant4SDActions.cpp b/DDG4/plugins/Geant4SDActions.cpp
index 608c6ca48c27ee7cae059a3e2ad0cbb68a985d35..511ceca1fc8bf741790ec1c93d665fd219626e35 100644
--- a/DDG4/plugins/Geant4SDActions.cpp
+++ b/DDG4/plugins/Geant4SDActions.cpp
@@ -7,91 +7,13 @@
//
//====================================================================
// Framework include files
-#include "DDG4/Geant4SensDetAction.h"
+#include "DDG4/Geant4SensDetAction.inl"
#include "DDG4/Geant4Data.h"
#include "DD4hep/Printout.h"
using namespace std;
-/*
- * DD4hep namespace declaration
- */
-namespace DD4hep {
-
- /*
- * Simulation namespace declaration
- */
- namespace Simulation {
-
- template <typename T> class Geant4SensitiveAction : public Geant4Sensitive {
- protected:
- /// Collection identifiers
- size_t m_collectionID;
- public:
- typedef SimpleHit::Contribution HitContribution;
- /// Standard , initializing constructor
- Geant4SensitiveAction(Geant4Context* context, const string& name, DetElement det, LCDD& lcdd);
- /// Default destructor
- virtual ~Geant4SensitiveAction();
- /// Define collections created by this sensitivie action object
- virtual void defineCollections() {}
- /// G4VSensitiveDetector interface: Method invoked at the begining of each event.
- virtual void begin(G4HCofThisEvent* hce);
- /// G4VSensitiveDetector interface: Method invoked at the end of each event.
- virtual void end(G4HCofThisEvent* hce);
- /// G4VSensitiveDetector interface: Method for generating hit(s) using the G4Step object.
- virtual bool process(G4Step* step,G4TouchableHistory* history);
- /// G4VSensitiveDetector interface: Method invoked if the event was aborted.
- virtual void clear(G4HCofThisEvent* hce);
- };
- }
-}
-#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4TouchableHandler.h"
-
-namespace DD4hep {
- namespace Simulation {
-
- /// Standard , initializing constructor
- template <typename T>
- Geant4SensitiveAction<T>::Geant4SensitiveAction(Geant4Context* context,
- const string& name,
- DetElement det,
- LCDD& lcdd)
- : Geant4Sensitive(context,name,det,lcdd), m_collectionID(0)
- {
- defineCollections();
- InstanceCount::increment(this);
- }
-
- /// Default destructor
- template <typename T> Geant4SensitiveAction<T>::~Geant4SensitiveAction() {
- InstanceCount::decrement(this);
- }
-
- /// G4VSensitiveDetector interface: Method invoked at the begining of each event.
- template <typename T> void Geant4SensitiveAction<T>::begin(G4HCofThisEvent* hce) {
- Geant4Sensitive::begin(hce);
- }
-
- /// G4VSensitiveDetector interface: Method invoked at the end of each event.
- template <typename T> void Geant4SensitiveAction<T>::end(G4HCofThisEvent* hce) {
- Geant4Sensitive::end(hce);
- }
-
- /// G4VSensitiveDetector interface: Method for generating hit(s) using the G4Step object.
- template <typename T> bool Geant4SensitiveAction<T>::process(G4Step* step,G4TouchableHistory* history) {
- return Geant4Sensitive::process(step,history);
- }
-
- /// G4VSensitiveDetector interface: Method invoked if the event was aborted.
- template <typename T> void Geant4SensitiveAction<T>::clear(G4HCofThisEvent* hce) {
- Geant4Sensitive::clear(hce);
- }
-
- }
-}
-
#include "DDG4/Geant4StepHandler.h"
#include "DDG4/Geant4VolumeManager.h"
#include "DDG4/Geant4Mapping.h"
@@ -102,30 +24,7 @@ namespace DD4hep {
namespace DD4hep {
namespace Simulation {
- template<typename TYPE> struct PositionCompare : public Geant4HitCollection::Compare {
- const Position& pos;
- Geant4HitWrapper::HitManipulator* manip;
- /// Constructor
- PositionCompare(const Position& p) : pos(p), manip(Geant4HitWrapper::HitManipulator::instance<TYPE>()) { }
- /// Comparison function
- virtual void* operator()(const Geant4HitWrapper& w) const {
- TYPE* h = w;//(const TYPE*)manip->castHit<TYPE>(w);
- return pos == h->position ? h : 0;
- }
- };
- template<typename TYPE> struct CellIDCompare : public Geant4HitCollection::Compare {
- long long int id;
- Geant4HitWrapper::HitManipulator* manip;
- /// Constructor
- CellIDCompare(long long int i) : id(i), manip(Geant4HitWrapper::HitManipulator::instance<TYPE>()) { }
- /// Comparison function
- virtual void* operator()(const Geant4HitWrapper& w) const {
- TYPE* h = w;//(const TYPE*)manip->castHit<TYPE>(w);
- if ( id == h->cellID )
- return h;
- return 0;
- }
- };
+ typedef Geant4HitData::Contribution HitContribution;
/// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
/// Geant4SensitiveAction<Geant4Tracker>
@@ -243,7 +142,7 @@ namespace DD4hep {
HitCollection* coll = collection(m_collectionID);
HitContribution contrib = Hit::extractContribution(step);
Position pos = h.prePos();
- Hit* hit = coll->find<Hit>(PositionCompare<Hit>(pos));
+ Hit* hit = coll->find<Hit>(PositionCompare<Hit,Position>(pos));
if ( !hit ) {
hit = new Hit(pos);
hit->cellID = volumeID(step);
diff --git a/DDG4/plugins/Geant4XMLSetup.cpp b/DDG4/plugins/Geant4XMLSetup.cpp
index 1f34bd040c9c49fa563a1661f5d478df186456d0..79b3635848478e3a7ca3f7ec5655f7abc9b732ec 100644
--- a/DDG4/plugins/Geant4XMLSetup.cpp
+++ b/DDG4/plugins/Geant4XMLSetup.cpp
@@ -426,7 +426,15 @@ namespace DD4hep {
template <> void Converter<XMLSetup>::operator()(xml_h seq) const {
xml_elt_t compact(seq);
// First execute the basic setup from the plugins module
- ROOT::Reflex::PluginService::Create<NamedObject*>("geant4_xml_setup",&lcdd,&seq);
+ long result = ROOT::Reflex::PluginService::Create<long>("geant4_XML_reader",&lcdd,&seq);
+ if ( 0 == result ) {
+ throw runtime_error("DD4hep: Failed to locate plugin to interprete files of type"
+ " \"" + seq.tag() + "\" - no factory of type geant4_XML_reader.");
+ }
+ result = *(long*) result;
+ if (result != 1) {
+ throw runtime_error("DD4hep: Failed to parse the XML tag " + seq.tag() + " with the plugin geant4_XML_reader");
+ }
xml_coll_t(compact,_Unicode(kernel)).for_each(Converter<Kernel>(lcdd,param));
// Now deal with the new stuff.....
xml_coll_t(compact,_Unicode(actions) ).for_each(_Unicode(action),Converter<Action>(lcdd,param));
diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py
index be43da09fb577a5fa2aff09789f2bbfdc4a60778..f863a9222ba17689946f71367bdeccee9c39c828 100644
--- a/DDG4/python/DDG4.py
+++ b/DDG4/python/DDG4.py
@@ -225,13 +225,17 @@ class Simple:
self.kernel.generatorAction().add(gun)
return gun
- def setupCshUI(self):
+ def setupUI(typ='csh',vis=False,ui=True):
# Configure UI
ui = Action(self.kernel,"Geant4UIManager/UI")
- ui.HaveVIS = True
- ui.HaveUI = True
- ui.SessionType = 'csh'
+ ui.HaveVIS = vis
+ ui.HaveUI = ui
+ ui.SessionType = typ
self.kernel.registerGlobalAction(ui)
+
+ def setupCshUI(self,typ='csh',vis=False,ui=True):
+ self.setupUI('csh',vis,ui)
+
def setupROOTOutput(self,name,output):
evt_root = EventAction(self.kernel,'Geant4Output2ROOT/'+name)
evt_root.Control = True
diff --git a/DDG4/src/Geant4GeneratorActionInit.cpp b/DDG4/src/Geant4GeneratorActionInit.cpp
index 5194dfb366f1581616535a6a17fd4665cfb880f7..09ee12e1d1c837a34eba79ee74f1775686994f5e 100644
--- a/DDG4/src/Geant4GeneratorActionInit.cpp
+++ b/DDG4/src/Geant4GeneratorActionInit.cpp
@@ -10,9 +10,9 @@
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4Kernel.h"
-#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4RunAction.h"
#include "DDG4/Geant4GeneratorActionInit.h"
+#include "DDG4/Geant4InputHandling.h"
#include "G4Run.hh"
@@ -52,30 +52,5 @@ void Geant4GeneratorActionInit::operator()(G4Event* /* event */) {
++m_evtRun;
/// + Printout
print("+++ Initializing event %d. Within run:%d event %d.",m_evtTotal,m_run,m_evtRun);
-
- /**
- * This action should register all event extension required for the further
- * processing. We want to avoid that every client has to check if a given
- * object is present or not and than later install the required data structures.
- */
- context()->event().addExtension(new Geant4PrimaryMap());
-
- // The final set of created particles in the simulation.
- context()->event().addExtension(new Geant4ParticleMap());
-
- //
- // The Geant4PrimaryEvent extension contains a whole set of
- // Geant4PrimaryInteraction objects each may represent a complete
- // interaction. Particles and vertices may be unbiased.
- // This is the input to the translator forming the final
- // Geant4PrimaryInteraction (see below) containing rebiased particle
- // and vertex maps.
- Geant4PrimaryEvent* evt = new Geant4PrimaryEvent();
- context()->event().addExtension(evt);
- //
- // The Geant4PrimaryInteraction extension contains the final
- // vertices and particles ready to be injected to Geant4.
- Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
- inter->setNextPID(-1);
- context()->event().addExtension(inter);
+ generationInitialization(this,context());
}
diff --git a/DDG4/src/Geant4HitCollection.cpp b/DDG4/src/Geant4HitCollection.cpp
index b1a1b4604b598f61e8d61ac139c7404a18a39e97..f6b9c94c5019eb9fca50e7b8ce5430469a207475 100644
--- a/DDG4/src/Geant4HitCollection.cpp
+++ b/DDG4/src/Geant4HitCollection.cpp
@@ -10,6 +10,7 @@
// Framework include files
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4HitCollection.h"
+#include "DDG4/Geant4Data.h"
#include "G4Allocator.hh"
using namespace DD4hep;
diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a2180100431e02634e6ad8dfb0105e20b838216
--- /dev/null
+++ b/DDG4/src/Geant4InputHandling.cpp
@@ -0,0 +1,338 @@
+// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $
+//====================================================================
+// AIDA Detector description implementation for LCD
+//--------------------------------------------------------------------
+//
+// Author : M.Frank
+//
+//====================================================================
+// Framework include files
+#include "DDG4/Geant4InputHandling.h"
+#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4Context.h"
+#include "DDG4/Geant4Action.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include "CLHEP/Units/PhysicalConstants.h"
+
+#include "G4ParticleDefinition.hh"
+#include "G4Event.hh"
+#include "G4PrimaryVertex.hh"
+#include "G4PrimaryParticle.hh"
+
+// C/C++ include files
+#include <stdexcept>
+#include <cmath>
+
+using namespace std;
+using namespace DD4hep;
+using namespace DD4hep::Simulation;
+
+typedef ReferenceBitMask<int> PropertyMask;
+
+/// Initialize the generation of one event
+int DD4hep::Simulation::generationInitialization(const Geant4Action* /* caller */,
+ const Geant4Context* context)
+{
+ /**
+ * This action should register all event extension required for the further
+ * processing. We want to avoid that every client has to check if a given
+ * object is present or not and than later install the required data structures.
+ */
+ context->event().addExtension(new Geant4PrimaryMap());
+
+ // The final set of created particles in the simulation.
+ context->event().addExtension(new Geant4ParticleMap());
+
+ //
+ // The Geant4PrimaryEvent extension contains a whole set of
+ // Geant4PrimaryInteraction objects each may represent a complete
+ // interaction. Particles and vertices may be unbiased.
+ // This is the input to the translator forming the final
+ // Geant4PrimaryInteraction (see below) containing rebiased particle
+ // and vertex maps.
+ Geant4PrimaryEvent* evt = new Geant4PrimaryEvent();
+ context->event().addExtension(evt);
+ //
+ // The Geant4PrimaryInteraction extension contains the final
+ // vertices and particles ready to be injected to Geant4.
+ Geant4PrimaryInteraction* inter = new Geant4PrimaryInteraction();
+ inter->setNextPID(-1);
+ context->event().addExtension(inter);
+ return 1;
+}
+
+/// Append input interaction to global output
+static void appendInteraction(const Geant4Action* caller,
+ Geant4PrimaryInteraction* output,
+ Geant4PrimaryInteraction* input)
+{
+ Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
+ for( ip=input->particles.begin(), ipend=input->particles.end(); ip != ipend; ++ip ) {
+ Geant4Particle* p = (*ip).second;
+ output->particles.insert(make_pair(p->id,p->addRef()));
+ }
+ Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
+ for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv ) {
+ ivfnd = output->vertices.find((*iv).second->mask);
+ if ( ivfnd != output->vertices.end() ) {
+ caller->abortRun("Duplicate primary interaction identifier!",
+ "Cannot handle 2 interactions with identical identifiers!");
+ }
+ output->vertices.insert(make_pair((*iv).second->mask,(*iv).second->addRef()));
+ }
+}
+
+static void rebaseParticles(Geant4PrimaryInteraction::ParticleMap& particles, int &offset) {
+ Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
+ int mx_id = offset;
+ // Now move begin and end-vertex of all primary and generator particles accordingly
+ for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip ) {
+ Geant4ParticleHandle p((*ip).second);
+ p.offset(offset);
+ mx_id = p->id+1 > mx_id ? p->id+1 : mx_id;
+ }
+ offset = mx_id;
+}
+
+static void rebaseVertices(Geant4PrimaryInteraction::VertexMap& vertices, int part_offset) {
+ Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
+ set<int> in, out;
+ set<int>::iterator i;
+ // Now move begin and end-vertex of all primary vertices accordingly
+ for(iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv) {
+ Geant4Vertex* v = (*iv).second;
+ in = v->in;
+ out = v->out;
+ for(in=v->in, i=in.begin(), v->in.clear(); i != in.end(); ++i)
+ v->in.insert((*i)+part_offset);
+ for(out=v->out, i=out.begin(), v->out.clear(); i != out.end(); ++i)
+ v->out.insert((*i)+part_offset);
+ }
+}
+
+/// Merge all interactions present in the context
+int DD4hep::Simulation::mergeInteractions(const Geant4Action* caller,
+ const Geant4Context* context)
+{
+ typedef Geant4PrimaryEvent::Interaction Interaction;
+ typedef vector<Interaction*> Interactions;
+ Geant4Event& event = context->event();
+ Geant4PrimaryEvent* evt = event.extension<Geant4PrimaryEvent>();
+ Interaction* output = event.extension<Interaction>();
+ Interactions inter = evt->interactions();
+ int particle_offset = 0;
+
+ for(Interactions::const_iterator i=inter.begin(); i != inter.end(); ++i) {
+ Interaction* interaction = *i;
+ int vertex_offset = particle_offset;
+ if ( !interaction->applyMask() ) {
+ caller->abortRun("Found single interaction with multiple primary vertices!",
+ "Cannot merge individual interactions with more than one primary!");
+ }
+ rebaseParticles(interaction->particles,particle_offset);
+ rebaseVertices(interaction->vertices,vertex_offset);
+ appendInteraction(caller,output,interaction);
+ }
+ output->setNextPID(particle_offset);
+ Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
+ caller->debug("+++ Merging MC input record from %d interactions:",(int)inter.size());
+ for( ip=output->particles.begin(), ipend=output->particles.end(); ip != ipend; ++ip )
+ Geant4ParticleHandle((*ip).second).dump1(DEBUG,caller->name(),"Merged particles");
+ return 1;
+}
+
+/// Boost particles of one interaction identified by its mask
+int DD4hep::Simulation::boostInteraction(const Geant4Action* /* caller */,
+ Geant4PrimaryEvent::Interaction* inter,
+ double alpha)
+{
+#define SQR(x) (x*x)
+ Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
+ Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
+ double gamma = std::sqrt(1 + SQR(tan(alpha)));
+ double betagamma = std::tan(alpha);
+
+ if ( alpha != 0.0 ) {
+ // Now move begin and end-vertex of all primary vertices accordingly
+ for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
+ Geant4Vertex* v = (*iv).second;
+ double t = gamma * v->time + betagamma * v->x / c_light;
+ double x = gamma * v->x + betagamma * c_light * v->time;
+ double y = v->y;
+ double z = v->z;
+ v->x = x;
+ v->y = y;
+ v->z = z;
+ v->time = t;
+ }
+
+ // Now move begin and end-vertex of all primary and generator particles accordingly
+ for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
+ Geant4Particle* p = (*ip).second;
+ double t = gamma * p->time + betagamma * p->vsx / c_light;
+ double x = gamma * p->vsx + betagamma * c_light * p->time;
+ double y = p->vsx;
+ double z = p->vsz;
+
+ double m = p->definition->GetPDGMass();
+ double e2 = SQR(p->psx)+SQR(p->psy)+SQR(p->psz)+SQR(m);
+ double px = betagamma * std::sqrt(e2) + gamma * p->psx;
+ double py = p->psy;
+ double pz = p->psz;
+
+ p->vsx = x;
+ p->vsy = y;
+ p->vsz = z;
+ p->time = t;
+
+ p->psx += px;
+ p->psy += py;
+ p->psz += pz;
+ }
+ }
+ return 1;
+}
+
+/// Smear the primary vertex of an interaction
+int DD4hep::Simulation::smearInteraction(const Geant4Action* /* caller */,
+ Geant4PrimaryEvent::Interaction* inter,
+ double dx, double dy, double dz, double dt)
+{
+ Geant4PrimaryEvent::Interaction::VertexMap::iterator iv;
+ Geant4PrimaryEvent::Interaction::ParticleMap::iterator ip;
+
+ // Now move begin and end-vertex of all primary vertices accordingly
+ for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
+ Geant4Vertex* v = (*iv).second;
+ v->x += dx;
+ v->y += dy;
+ v->z += dz;
+ v->time += dt;
+ }
+ // Now move begin and end-vertex of all primary and generator particles accordingly
+ for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
+ Geant4Particle* p = (*ip).second;
+ p->vsx += dx;
+ p->vsy += dy;
+ p->vsz += dz;
+ p->vex += dx;
+ p->vey += dy;
+ p->vez += dz;
+ p->time += dt;
+ }
+ return 1;
+}
+
+static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p) {
+ G4PrimaryParticle* g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz);
+ g4->SetMass(p->mass);
+ return g4;
+}
+
+static map<Geant4Particle*,G4PrimaryParticle*>
+getRelevant(set<int>& visited,
+ map<int,G4PrimaryParticle*>& prim,
+ Geant4PrimaryInteraction::ParticleMap& pm,
+ const Geant4ParticleHandle p)
+{
+ typedef map<Geant4Particle*,G4PrimaryParticle*> Primaries;
+ Primaries res;
+ visited.insert(p->id);
+ PropertyMask status(p->status);
+ if ( status.isSet(G4PARTICLE_GEN_STABLE) ) {
+ if ( prim.find(p->id) == prim.end() ) {
+ G4PrimaryParticle* p4 = createG4Primary(p);
+ prim[p->id] = p4;
+ res.insert(make_pair(p,p4));
+ }
+ }
+ else if ( p->daughters.size() > 0 ) {
+ const Geant4Particle::Particles& dau = p->daughters;
+ int first_daughter = *(dau.begin());
+ Geant4ParticleHandle dp = pm[first_daughter];
+ double en = p.energy();
+ double me = en > std::numeric_limits<double>::epsilon() ? p->mass / en : 0.0;
+ // fix by S.Morozov for real != 0
+ double proper_time = fabs(dp->time-p->time) * me;
+ double proper_time_Precision = pow(10.,-DBL_DIG)*me*fmax(fabs(p->time),fabs(dp->time));
+ bool isProperTimeZero = ( proper_time <= proper_time_Precision ) ;
+ // -- remove original --- if (proper_time != 0) {
+ if ( !isProperTimeZero ) {
+ map<int,G4PrimaryParticle*>::iterator ip4 = prim.find(p->id);
+ G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
+ if ( !p4 ) {
+ p4 = createG4Primary(p);
+ p4->SetProperTime(proper_time);
+ prim[p->id] = p4;
+ Primaries daughters;
+ for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
+ if ( visited.find(*i) == visited.end() ) {
+ Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
+ daughters.insert(tmp.begin(),tmp.end());
+ }
+ }
+ for(Primaries::iterator i=daughters.begin(); i!=daughters.end(); ++i)
+ p4->SetDaughter((*i).second);
+ }
+ res.insert(make_pair(p,p4));
+ }
+ else {
+ for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
+ if ( visited.find(*i) == visited.end() ) {
+ Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
+ res.insert(tmp.begin(),tmp.end());
+ }
+ }
+ }
+ }
+ return res;
+}
+
+/// Generate all primary vertices corresponding to the merged interaction
+int DD4hep::Simulation::generatePrimaries(const Geant4Action* caller,
+ const Geant4Context* context,
+ G4Event* event)
+{
+ typedef map<Geant4Particle*,G4PrimaryParticle*> Primaries;
+ typedef Geant4PrimaryInteraction Interaction;
+ Geant4PrimaryMap* primaries = context->event().extension<Geant4PrimaryMap>();
+ Interaction* interaction = context->event().extension<Interaction>();
+ Interaction::ParticleMap& pm = interaction->particles;
+ Interaction::VertexMap& vm = interaction->vertices;
+ map<int,G4PrimaryParticle*> prim;
+ set<int> visited;
+
+ Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
+ for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i) {
+ int num_part = 0;
+ Geant4Vertex* v = (*i).second;
+ G4PrimaryVertex* v4 = new G4PrimaryVertex(v->x,v->y,v->z,v->time);
+ event->AddPrimaryVertex(v4);
+ caller->print("+++++ G4PrimaryVertex at (%+.2e,%+.2e,%+.2e) [mm] %+.2e [ns]",
+ v->x/mm,v->y/mm,v->z/mm,v->time/ns);
+ for(Geant4Vertex::Particles::const_iterator ip=v->out.begin(); ip!=v->out.end(); ++ip) {
+ Geant4ParticleHandle p = pm[*ip];
+ if ( p->parents.size() == 0 ) {
+ Primaries relevant = getRelevant(visited,prim,pm,p);
+ for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j) {
+ Geant4ParticleHandle r = (*j).first;
+ G4PrimaryParticle* p4 = (*j).second;
+ PropertyMask reason(r->reason);
+ reason.set(G4PARTICLE_PRIMARY);
+ v4->SetPrimary(p4);
+ caller->printM1("+++ +-> G4Primary[%3d] ID:%3d type:%9d/%-12s "
+ "Momentum:(%+.2e,%+.2e,%+.2e) [GeV] time:%+.2e [ns] #Par:%3d #Dau:%3d",
+ num_part,r->id,r->pdgID,r.particleName().c_str(),
+ r->psx/GeV,r->psy/GeV,r->psz/GeV,r->time/ns,
+ int(r->parents.size()),int(r->daughters.size()));
+ ++num_part;
+ }
+ }
+ }
+ }
+ for(map<int,G4PrimaryParticle*>::iterator i=prim.begin(); i!=prim.end(); ++i) {
+ Geant4ParticleHandle p = pm[(*i).first];
+ primaries->primaryMap[(*i).second] = p->addRef();
+ }
+ return 1;
+}
diff --git a/DDG4/src/Geant4InteractionMerger.cpp b/DDG4/src/Geant4InteractionMerger.cpp
index 840c7efcc157c1486317cdb9dd562fd39e2a8831..7423bca01b612e88d85ce1abc0eb7914d656aa3a 100644
--- a/DDG4/src/Geant4InteractionMerger.cpp
+++ b/DDG4/src/Geant4InteractionMerger.cpp
@@ -9,13 +9,9 @@
// Framework include files
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
-#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4InputHandling.h"
#include "DDG4/Geant4InteractionMerger.h"
-#include <stdexcept>
-
-using namespace std;
-using namespace DD4hep;
using namespace DD4hep::Simulation;
/// Standard constructor
@@ -30,75 +26,10 @@ Geant4InteractionMerger::~Geant4InteractionMerger() {
InstanceCount::decrement(this);
}
-void rebaseParticles(Geant4PrimaryInteraction::ParticleMap& particles, int &offset) {
- Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
- int mx_id = offset;
- // Now move begin and end-vertex of all primary and generator particles accordingly
- for( ip=particles.begin(), ipend=particles.end(); ip != ipend; ++ip ) {
- Geant4ParticleHandle p((*ip).second);
- p.offset(offset);
- mx_id = p->id+1 > mx_id ? p->id+1 : mx_id;
- }
- offset = mx_id;
-}
-
-void rebaseVertices(Geant4PrimaryInteraction::VertexMap& vertices, int part_offset) {
- Geant4PrimaryInteraction::VertexMap::iterator iv, ivend;
- set<int> in, out;
- set<int>::iterator i;
- // Now move begin and end-vertex of all primary vertices accordingly
- for(iv=vertices.begin(), ivend=vertices.end(); iv != ivend; ++iv) {
- Geant4Vertex* v = (*iv).second;
- in = v->in;
- out = v->out;
- for(in=v->in, i=in.begin(), v->in.clear(); i != in.end(); ++i)
- v->in.insert((*i)+part_offset);
- for(out=v->out, i=out.begin(), v->out.clear(); i != out.end(); ++i)
- v->out.insert((*i)+part_offset);
- }
-}
-
-/// Append input interaction to global output
-void Geant4InteractionMerger::appendInteraction(Geant4PrimaryInteraction* output, Geant4PrimaryInteraction* input) {
- Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
- for( ip=input->particles.begin(), ipend=input->particles.end(); ip != ipend; ++ip ) {
- Geant4Particle* p = (*ip).second;
- output->particles.insert(make_pair(p->id,p->addRef()));
- }
- Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
- for( iv=input->vertices.begin(), ivend=input->vertices.end(); iv != ivend; ++iv ) {
- ivfnd = output->vertices.find((*iv).second->mask);
- if ( ivfnd != output->vertices.end() ) {
- abortRun("Duplicate primary interaction identifier!",
- "Cannot handle 2 interactions with identical identifiers!");
- }
- output->vertices.insert(make_pair((*iv).second->mask,(*iv).second->addRef()));
- }
-}
-
/// Event generation action callback
void Geant4InteractionMerger::operator()(G4Event* /* event */) {
- typedef Geant4PrimaryEvent::Interaction Interaction;
- typedef vector<Interaction*> Interactions;
Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
- Interaction* output = context()->event().extension<Interaction>();
- Interactions inter = evt->interactions();
- int particle_offset = 0;
-
- for(Interactions::const_iterator i=inter.begin(); i != inter.end(); ++i) {
- Interaction* interaction = *i;
- int vertex_offset = particle_offset;
- if ( !interaction->applyMask() ) {
- abortRun("Found single interaction with multiple primary vertices!",
- "Cannot merge individual interactions with more than one primary!");
- }
- rebaseParticles(interaction->particles,particle_offset);
- rebaseVertices(interaction->vertices,vertex_offset);
- appendInteraction(output,interaction);
- }
- output->setNextPID(particle_offset);
- Geant4PrimaryInteraction::ParticleMap::iterator ip, ipend;
- debug("+++ Merged MC input record from %d interactions:",(int)inter.size());
- for( ip=output->particles.begin(), ipend=output->particles.end(); ip != ipend; ++ip )
- Geant4ParticleHandle((*ip).second).dump1(DEBUG,name(),"Merged particles");
+ const std::vector<Geant4PrimaryEvent::Interaction*>& inter = evt->interactions();
+ debug("+++ Merging MC input record from %d interactions:",(int)inter.size());
+ mergeInteractions(this,context());
}
diff --git a/DDG4/src/Geant4InteractionVertexBoost.cpp b/DDG4/src/Geant4InteractionVertexBoost.cpp
index d9b889cd1f856ddea522cc6f062e86d9d3bd1aee..1960fd9d1822ecbaa225c323504ef1e0d37d887b 100644
--- a/DDG4/src/Geant4InteractionVertexBoost.cpp
+++ b/DDG4/src/Geant4InteractionVertexBoost.cpp
@@ -10,25 +10,13 @@
// Framework include files
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
-#include "DDG4/Geant4Random.h"
-#include "DDG4/Geant4Context.h"
-#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4InputHandling.h"
#include "DDG4/Geant4InteractionVertexBoost.h"
-#include "CLHEP/Units/SystemOfUnits.h"
-#include "CLHEP/Units/PhysicalConstants.h"
-#include "G4ParticleDefinition.hh"
-
-// C/C++ include files
-#include <stdexcept>
-#include <cmath>
-
-using namespace std;
using namespace DD4hep::Simulation;
-#define SQR(x) (x*x)
/// Standard constructor
-Geant4InteractionVertexBoost::Geant4InteractionVertexBoost(Geant4Context* context, const string& name)
+Geant4InteractionVertexBoost::Geant4InteractionVertexBoost(Geant4Context* context, const std::string& name)
: Geant4GeneratorAction(context, name)
{
InstanceCount::increment(this);
@@ -44,56 +32,9 @@ Geant4InteractionVertexBoost::~Geant4InteractionVertexBoost() {
/// Callback to generate primary particles
void Geant4InteractionVertexBoost::operator()(G4Event*) {
- typedef Geant4PrimaryEvent::Interaction Interaction;
- typedef Geant4Particle Particle;
- typedef Geant4Particle Vertex;
- Geant4PrimaryEvent* evt = context()->event().extension<Geant4PrimaryEvent>();
- Interaction* inter = evt->get(m_mask);
-
+ Geant4PrimaryEvent::Interaction* inter =
+ context()->event().extension<Geant4PrimaryEvent>()->get(m_mask);
if ( inter ) {
- Interaction::VertexMap::iterator iv;
- Interaction::ParticleMap::iterator ip;
- double alpha = m_angle;
- double gamma = std::sqrt(1 + SQR(tan(alpha)));
- double betagamma = std::tan(alpha);
-
- if ( alpha != 0.0 ) {
- // Now move begin and end-vertex of all primary vertices accordingly
- for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
- Geant4Vertex* v = (*iv).second;
- double t = gamma * t + betagamma * v->x / c_light;
- double x = gamma * v->x + betagamma * c_light * v->time;
- double y = v->y;
- double z = v->z;
- v->x = x;
- v->y = y;
- v->z = z;
- v->time = t;
- }
-
- // Now move begin and end-vertex of all primary and generator particles accordingly
- for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
- Particle* p = (*ip).second;
- double t = gamma * p->time + betagamma * p->vsx / c_light;
- double x = gamma * p->vsx + betagamma * c_light * p->time;
- double y = p->vsx;
- double z = p->vsz;
-
- double m = p->definition->GetPDGMass();
- double e2 = SQR(p->psx)+SQR(p->psy)+SQR(p->psz)+SQR(m);
- double px = betagamma * std::sqrt(e2) + gamma * p->psx;
- double py = p->psy;
- double pz = p->psz;
-
- p->vsx = x;
- p->vsy = y;
- p->vsz = z;
- p->time = t;
-
- p->psx += px;
- p->psy += py;
- p->psz += pz;
- }
- }
+ boostInteraction(this, inter, m_angle);
}
}
diff --git a/DDG4/src/Geant4InteractionVertexSmear.cpp b/DDG4/src/Geant4InteractionVertexSmear.cpp
index 68a0f824fd121fd21c8cd3f5aa1a7b3d245a22b6..a531b6490d69e395ccc8664f96af17e04885e919 100644
--- a/DDG4/src/Geant4InteractionVertexSmear.cpp
+++ b/DDG4/src/Geant4InteractionVertexSmear.cpp
@@ -12,19 +12,16 @@
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4Random.h"
#include "DDG4/Geant4Context.h"
-#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4InputHandling.h"
#include "DDG4/Geant4InteractionVertexSmear.h"
-#include "CLHEP/Units/SystemOfUnits.h"
// C/C++ include files
-#include <stdexcept>
#include <cmath>
-using namespace std;
using namespace DD4hep::Simulation;
/// Standard constructor
-Geant4InteractionVertexSmear::Geant4InteractionVertexSmear(Geant4Context* context, const string& name)
+Geant4InteractionVertexSmear::Geant4InteractionVertexSmear(Geant4Context* context, const std::string& name)
: Geant4GeneratorAction(context, name)
{
InstanceCount::increment(this);
@@ -49,36 +46,14 @@ void Geant4InteractionVertexSmear::operator()(G4Event*) {
Interaction* inter = evt->get(m_mask);
if ( inter ) {
- Interaction::VertexMap::iterator iv;
- Interaction::ParticleMap::iterator ip;
double dx = rndm.gauss(m_offset.x(),m_sigma.x());
double dy = rndm.gauss(m_offset.y(),m_sigma.y());
double dz = rndm.gauss(m_offset.z(),m_sigma.z());
double dt = rndm.gauss(m_offset.t(),m_sigma.t());
-
print("+++ Smearing primary vertex for interaction type %d (%d Vertices, %d particles) "
"by (%+.2e mm, %+.2e mm, %+.2e mm, %+.2e ns)",
m_mask,int(inter->vertices.size()),int(inter->particles.size()),dx,dy,dz,dt);
-
- // Now move begin and end-vertex of all primary vertices accordingly
- for(iv=inter->vertices.begin(); iv != inter->vertices.end(); ++iv) {
- Geant4Vertex* v = (*iv).second;
- v->x += dx;
- v->y += dy;
- v->z += dz;
- v->time += dt;
- }
- // Now move begin and end-vertex of all primary and generator particles accordingly
- for(ip=inter->particles.begin(); ip != inter->particles.end(); ++ip) {
- Particle* p = (*ip).second;
- p->vsx += dx;
- p->vsy += dy;
- p->vsz += dz;
- p->vex += dx;
- p->vey += dy;
- p->vez += dz;
- p->time += dt;
- }
+ smearInteraction(this,inter,dx,dy,dz,dt);
}
else {
print("+++ No interaction of type %d present.",m_mask);
diff --git a/DDG4/src/Geant4Output2ROOT.cpp b/DDG4/src/Geant4Output2ROOT.cpp
index 85a933ad7136b74d0c7ff158e9936baced9222ec..dc44761191df2ae57d1d098a4eed4324a3b80f9d 100644
--- a/DDG4/src/Geant4Output2ROOT.cpp
+++ b/DDG4/src/Geant4Output2ROOT.cpp
@@ -161,18 +161,18 @@ void Geant4Output2ROOT::saveCollection(OutputContext<G4Event>& /* ctxt */, G4VHi
if ( truth && nhits > 0 ) {
try {
for(size_t i=0; i<nhits; ++i) {
- SimpleHit* h = coll->hit(i);
- SimpleTracker::Hit* trk_hit = dynamic_cast<SimpleTracker::Hit*>(h);
+ Geant4HitData* h = coll->hit(i);
+ Geant4Tracker::Hit* trk_hit = dynamic_cast<Geant4Tracker::Hit*>(h);
if ( 0 != trk_hit ) {
- SimpleHit::Contribution& t = trk_hit->truth;
+ Geant4HitData::Contribution& t = trk_hit->truth;
int trackID = t.trackID;
t.trackID = truth->particleID(trackID);
}
- SimpleCalorimeter::Hit* cal_hit = dynamic_cast<SimpleCalorimeter::Hit*>(h);
+ Geant4Calorimeter::Hit* cal_hit = dynamic_cast<Geant4Calorimeter::Hit*>(h);
if ( 0 != cal_hit ) {
- SimpleHit::Contributions& c = cal_hit->truth;
- for(SimpleHit::Contributions::iterator j=c.begin(); j!=c.end(); ++j) {
- SimpleHit::Contribution& t = *j;
+ Geant4HitData::Contributions& c = cal_hit->truth;
+ for(Geant4HitData::Contributions::iterator j=c.begin(); j!=c.end(); ++j) {
+ Geant4HitData::Contribution& t = *j;
int trackID = t.trackID;
t.trackID = truth->particleID(trackID);
}
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
index cb70c4434cad0a5a2fda90cf5036da0b3daffebf..69b10132ae13f543344e10d8ac5542e1d11f0f3a 100644
--- a/DDG4/src/Geant4Particle.cpp
+++ b/DDG4/src/Geant4Particle.cpp
@@ -20,6 +20,7 @@
using namespace DD4hep;
using namespace DD4hep::Simulation;
+typedef ReferenceBitMask<int> PropertyMask;
/// Default destructor
ParticleExtension::~ParticleExtension() {
@@ -276,6 +277,45 @@ void Geant4ParticleHandle::dump3(int level, const std::string& src, const char*
text);
}
+void Geant4ParticleHandle::dump4(int level, const std::string& src, const char* tag) const {
+ Geant4ParticleHandle p(*this);
+ char equiv[32];
+ PropertyMask mask(p->reason);
+ PropertyMask status(p->status);
+ std::string proc_name = p.processName();
+ std::string proc_type = p.processTypeName();
+ int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
+
+ equiv[0] = 0;
+ if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
+ ::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
+ }
+ printout((DD4hep::PrintLevel)level,src,
+ "+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s] %c%c%c%c",
+ tag,
+ p->id,
+ p.particleName().c_str(),
+ parent_id,equiv,
+ yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
+ yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
+ int(p->daughters.size()),
+ yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
+ p.energy(),
+ yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
+ mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
+ p.numParent(),
+ proc_name.c_str(),
+ p->process ? "/" : "",
+ proc_type.c_str(),
+ status.isSet(G4PARTICLE_GEN_EMPTY) ? 'E' : '.',
+ status.isSet(G4PARTICLE_GEN_STABLE) ? 'S' : '.',
+ status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
+ status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.'
+ );
+}
+
/// Default destructor
Geant4ParticleMap::~Geant4ParticleMap() {
clear();
diff --git a/DDG4/src/Geant4ParticleGun.cpp b/DDG4/src/Geant4ParticleGun.cpp
index 290545cf7f6a2842a0161f82976e578f10fde7c3..216ec2d4943bfed05a14bbcc5dc49becee30588e 100644
--- a/DDG4/src/Geant4ParticleGun.cpp
+++ b/DDG4/src/Geant4ParticleGun.cpp
@@ -14,6 +14,7 @@
#include "DDG4/Geant4Random.h"
#include "DDG4/Geant4Primary.h"
#include "DDG4/Geant4ParticleGun.h"
+#include "DDG4/Geant4InputHandling.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "G4ParticleTable.hh"
@@ -41,6 +42,7 @@ Geant4ParticleGun::Geant4ParticleGun(Geant4Context* context, const string& name)
declareProperty("direction", m_direction);
declareProperty("isotrop", m_isotrop = false);
declareProperty("Mask", m_mask = 1);
+ declareProperty("Standalone", m_standalone = true);
}
/// Default destructor
@@ -49,8 +51,12 @@ Geant4ParticleGun::~Geant4ParticleGun() {
}
/// Callback to generate primary particles
-void Geant4ParticleGun::operator()(G4Event* ) {
+void Geant4ParticleGun::operator()(G4Event* event) {
typedef DD4hep::ReferenceBitMask<int> PropertyMask;
+ if ( m_standalone ) {
+ generationInitialization(this,context());
+ }
+
Geant4Event& evt = context()->event();
Geant4PrimaryEvent* prim = evt.extension<Geant4PrimaryEvent>();
if (0 == m_particle || m_particle->GetParticleName() != m_particleName.c_str()) {
@@ -116,4 +122,9 @@ void Geant4ParticleGun::operator()(G4Event* ) {
p.dump3(outputLevel()-1,name(),"+->");
}
++m_shotNo;
+
+ if ( m_standalone ) {
+ mergeInteractions(this,context());
+ generatePrimaries(this,context(),event);
+ }
}
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index 67d3eb05a82da7e912a9869490e6f947cc4306c9..1b7f1c3b085e281d65152863844bf5d2f2fe3299 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -69,6 +69,8 @@ Geant4ParticleHandler::Geant4ParticleHandler(Geant4Context* context, const strin
declareProperty("KeepAllParticles", m_keepAll = false);
declareProperty("SaveProcesses", m_processNames);
declareProperty("MinimalKineticEnergy",m_kinEnergyCut = 100e0*MeV);
+ declareProperty("Z_trackers",m_zTracker=1e100);
+ declareProperty("R_trackers",m_rTracker=1e100);
InstanceCount::increment(this);
}
@@ -155,6 +157,7 @@ void Geant4ParticleHandler::mark(const G4Track* track) {
else if ( !mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT) ) {
mask.set(G4PARTICLE_CREATED_TRACKER_HIT);
}
+ //Geant4ParticleHandle(&m_currTrack).dump4(outputLevel(),vol->GetName(),"hit created by particle");
}
/// Event generation action callback
@@ -172,13 +175,8 @@ void Geant4ParticleHandler::operator()(G4Event* event) {
/// User stepping callback
void Geant4ParticleHandler::step(const G4Step* step, G4SteppingManager* mgr) {
typedef vector<const G4Track*> _Sec;
- Geant4StepHandler h(step);
++m_currTrack.steps;
- const G4ThreeVector& v = h.postPosG4();
- m_currTrack.vex = v.x();
- m_currTrack.vey = v.y();
- m_currTrack.vez = v.z();
- if ( h.trkKineEnergy() > m_kinEnergyCut ) {
+ if ( (m_currTrack.reason&G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) {
//
// Tracks below the energy threshold are NOT stored.
// If these tracks produce hits or are selected due to another signature,
@@ -199,13 +197,14 @@ void Geant4ParticleHandler::step(const G4Step* step, G4SteppingManager* mgr) {
void Geant4ParticleHandler::begin(const G4Track* track) {
Geant4TrackHandler h(track);
double kine = h.kineticEnergy();
- G4ThreeVector m = track->GetMomentum();
+ G4ThreeVector m = h.momentum();
const G4ThreeVector& v = h.vertex();
int reason = (kine > m_kinEnergyCut) ? G4PARTICLE_ABOVE_ENERGY_THRESHOLD : 0;
G4PrimaryParticle* prim = primary(h.id(),context()->event().event());
Particle* prim_part = 0;
+
if ( prim ) {
Geant4PrimaryMap::Primaries::const_iterator iprim = m_primaryMap->primaryMap.find(prim);
if ( iprim == m_primaryMap->primaryMap.end() ) {
@@ -247,12 +246,6 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
m_currTrack.definition = h.trackDef();
m_currTrack.pdgID = h.trackDef()->GetPDGEncoding();
m_currTrack.mass = h.trackDef()->GetPDGMass();
- // Once the daughter gets simulated, the parent is already done....
- ParticleMap::iterator ipar = m_particleMap.find(h.parent());
- if ( ipar != m_particleMap.end() ) {
- Particle* p_par = (*ipar).second;
- m_currTrack.parents.insert(p_par->id);
- }
++m_globalParticleID;
}
m_currTrack.steps = 0;
@@ -292,35 +285,35 @@ void Geant4ParticleHandler::begin(const G4Track* track) {
/// Post-track action callback
void Geant4ParticleHandler::end(const G4Track* track) {
Geant4TrackHandler h(track);
+ Geant4ParticleHandle ph(&m_currTrack);
int g4_id = h.id();
+ int track_reason = m_currTrack.reason;
PropertyMask mask(m_currTrack.reason);
+ double r_prod = sqrt(ph->vsx*ph->vsx + ph->vsy*ph->vsy);
+ double z_prod = fabs(ph->vsz);
+ if ( r_prod > m_rTracker || z_prod > m_zTracker ) {
+ m_currTrack.reason = 0;
+ }
+
+ // These are candate tracks with a probability to be stored due to their properties:
+ // - primary particle
+ // - hits created
+ // - secondaries
+ // - above energy threshold
+ // - to be kept due to creator process
+ //
if ( !mask.isNull() ) {
- // These are candate tracks with a probability to be stored due to their properties:
- // - primary particle
- // - hits created
- // - secondaries
- // - above energy threshold
- // - to be kept due to creator process
- //
// Update vertex end point and final momentum
- Geant4ParticleHandle ph(&m_currTrack);
G4ThreeVector m = track->GetMomentum();
+ const G4ThreeVector& p = track->GetPosition();
ph->pex = m.x();
ph->pey = m.y();
ph->pez = m.z();
+ ph->vex = p.x();
+ ph->vey = p.y();
+ ph->vez = p.z();
m_equivalentTracks[g4_id] = g4_id;
- // Add this track to it's parents list of daughters
- ParticleMap::iterator ipar = m_particleMap.find(ph->g4Parent);
- if ( ipar != m_particleMap.end() ) {
- Particle* p_par = (*ipar).second;
- p_par->daughters.insert(ph->id);
- ph->mask |= p_par->mask;
- G4ThreeVector dist(ph->vsx-p_par->vex,ph->vsy-p_par->vey,ph->vsz-p_par->vez);
- if ( dist.r() > 1e-15 ) {
- // Set flag that the end vertex of the parent is not the start vertex of this track....
- }
- }
/// Final update of the particle using the user handler
if ( m_userHandler ) {
m_userHandler->begin(track, m_currTrack);
@@ -330,15 +323,30 @@ void Geant4ParticleHandler::end(const G4Track* track) {
ph.dump2(outputLevel()-1,name(),"Add Primary",h.id(),ip!=m_particleMap.end());
}
// Create a new MC particle from the current track information saved in the pre-tracking action
- Particle* p = ip==m_particleMap.end() ? (m_particleMap[g4_id] = new Particle()) : (*ip).second;
- p->get_data(m_currTrack);
+ Particle* part = 0;
+ if ( ip==m_particleMap.end() ) part = m_particleMap[g4_id] = new Particle();
+ else part = (*ip).second;
+ part->get_data(m_currTrack);
}
else {
// These are tracks without any special properties.
//
// We will not store them on the record, but have to memorise the
// track identifier in order to restore the history for the created hits.
- m_equivalentTracks[g4_id] = m_currTrack.g4Parent;
+ int pid = m_currTrack.g4Parent;
+ m_equivalentTracks[g4_id] = pid;
+ // Need to find the last stored particle and OR this particle's mask
+ // with the mask of the last stored particle
+ TrackEquivalents::const_iterator iequiv, iend = m_equivalentTracks.end();
+ ParticleMap::iterator ip;
+ for(ip=m_particleMap.find(pid); ip == m_particleMap.end(); ip=m_particleMap.find(pid)) {
+ if ((iequiv=m_equivalentTracks.find(pid)) == iend) break; // ERROR
+ pid = (*iequiv).second;
+ }
+ if ( ip != m_particleMap.end() )
+ (*ip).second->reason |= track_reason;
+ else
+ ph.dump3(outputLevel()+3,name(),"No real parent present");
}
}
@@ -358,8 +366,9 @@ void Geant4ParticleHandler::beginEvent(const G4Event* event) {
/// Debugging: Dump Geant4 particle map
void Geant4ParticleHandler::dumpMap(const char* tag) const {
- for(ParticleMap::const_iterator iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i)
- Geant4ParticleHandle((*i).second).dump1(INFO,name(),tag);
+ for(ParticleMap::const_iterator iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
+ Geant4ParticleHandle((*i).second).dump4(INFO,name(),tag);
+ }
}
/// Post-event action callback
@@ -438,30 +447,27 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
else
error("+++ No Equivalent particle for track:%d last known is:%d",(*i).second,g4_equiv);
}
- m_equivalentTracks = equivalents;
- equivalents.clear();
- // (3) Update the particle's parents and replace the original Geant4 track with the
+ // (3) Compute the particle's parents and daughters.
+ // Replace the original Geant4 track with the
// equivalent particle still present in the record.
for(ParticleMap::const_iterator ipar, iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
Particle* p = (*i).second;
- set<int> daughters = p->daughters;
- p->daughters.clear();
- for(set<int>::iterator id=daughters.begin(); id!=daughters.end(); ++id)
- p->daughters.insert(orgParticles[*id]); // Requires (1)
- //if ( 0 == (p->status&G4PARTICLE_GEN_GENERATOR) )
- {
- if ( p->g4Parent > 0 ) {
- p->parents.clear();
- p->parents.insert(equivalentTrack(p->g4Parent)); // Requires (2)
- ipar = m_particleMap.find(p->g4Parent);
- if ( ipar != iend ) {
- set<int>& dau = (*ipar).second->daughters;
- if ( dau.find(p->id) == dau.end() ) dau.insert(p->id);
- }
+ if ( p->g4Parent > 0 ) {
+ int equiv_id = equivalents[p->g4Parent];
+ if ( (ipar=finalParticles.find(equiv_id)) != finalParticles.end() ) {
+ Particle* q = (*ipar).second;
+ q->daughters.insert(p->id);
+ p->parents.insert(q->id);
+ }
+ else {
+ error("+++ Inconsistency in particle record: Geant4 parent %d "
+ "of particle %d (equiv:%d) not in record!",
+ p->g4Parent,p->id,equiv_id);
}
}
}
+ m_equivalentTracks = equivalents;
m_particleMap = finalParticles;
}
@@ -537,8 +543,6 @@ int Geant4ParticleHandler::recombineParents() {
if ( ip != m_particleMap.end() ) {
Particle* parent_part = (*ip).second;
PropertyMask(parent_part->reason).set(mask.value());
- // Remove track from parent's list of daughters
- parent_part->removeDaughter(id);
parent_part->steps += p->steps;
parent_part->secondaries += p->secondaries;
/// Update of the particle using the user handler
@@ -615,7 +619,7 @@ int Geant4ParticleHandler::equivalentTrack(int g4_id) const {
int equiv_id = g4_id;
if ( g4_id != 0 ) {
TrackEquivalents::const_iterator iequiv = m_equivalentTracks.find(equiv_id);
- return (iequiv == m_equivalentTracks.end()) ? -1 : (*iequiv).second;
+ if ( iequiv != m_equivalentTracks.end() ) return (*iequiv).second;
}
return -1;
}
diff --git a/DDG4/src/Geant4ParticlePrint.cpp b/DDG4/src/Geant4ParticlePrint.cpp
index af5b4e3efd07befb698e2b25fa462a5ef3e02799..27e082177162c85bcdf306b7de8d580e69aa58a2 100644
--- a/DDG4/src/Geant4ParticlePrint.cpp
+++ b/DDG4/src/Geant4ParticlePrint.cpp
@@ -12,6 +12,7 @@
#include "DD4hep/InstanceCount.h"
#include "DDG4/Geant4ParticlePrint.h"
#include "DDG4/Geant4Data.h"
+#include "DDG4/Geant4HitCollection.h"
#include <cstring>
#include "G4Event.hh"
@@ -30,6 +31,7 @@ Geant4ParticlePrint::Geant4ParticlePrint(Geant4Context* context, const std::stri
declareProperty("PrintBegin",m_printBegin=false);
declareProperty("PrintEnd", m_printEnd=true);
declareProperty("PrintGen", m_printGeneration=true);
+ declareProperty("PrintHits", m_printHits=false);
InstanceCount::increment(this);
}
@@ -39,13 +41,13 @@ Geant4ParticlePrint::~Geant4ParticlePrint() {
}
/// Print particle table
-void Geant4ParticlePrint::makePrintout(int event_id) const {
+void Geant4ParticlePrint::makePrintout(const G4Event* e) const {
Geant4ParticleMap* parts = context()->event().extension<Geant4ParticleMap>();
if ( parts ) {
const ParticleMap& particles = parts->particles();
- print("+++ ******** MC Particle Printout for event ID:%d ********",event_id);
- if ( (m_outputType&2) != 0 ) printParticleTree(particles); // Tree printout....
- if ( (m_outputType&1) != 0 ) printParticles(particles); // Table printout....
+ print("+++ ******** MC Particle Printout for event ID:%d ********",e->GetEventID());
+ if ( (m_outputType&2) != 0 ) printParticleTree(e,particles); // Tree printout....
+ if ( (m_outputType&1) != 0 ) printParticles(0,particles); // Table printout....
return;
}
except("No Geant4MonteCarloTruth object present in the event structure to access the particle information!", c_name());
@@ -53,20 +55,20 @@ void Geant4ParticlePrint::makePrintout(int event_id) const {
/// Generation action callback
void Geant4ParticlePrint::operator()(G4Event* e) {
- if ( m_printGeneration ) makePrintout(e->GetEventID());
+ if ( m_printGeneration ) makePrintout(e);
}
/// Pre-event action callback
void Geant4ParticlePrint::begin(const G4Event* e) {
- if ( m_printBegin ) makePrintout(e->GetEventID());
+ if ( m_printBegin ) makePrintout(e);
}
/// Post-event action callback
void Geant4ParticlePrint::end(const G4Event* e) {
- if ( m_printEnd ) makePrintout(e->GetEventID());
+ if ( m_printEnd ) makePrintout(e);
}
-void Geant4ParticlePrint::printParticle(const std::string& prefix, Geant4ParticleHandle p) const {
+void Geant4ParticlePrint::printParticle(const std::string& prefix, const G4Event* e, Geant4ParticleHandle p) const {
char equiv[32];
PropertyMask mask(p->reason);
PropertyMask status(p->status);
@@ -101,10 +103,45 @@ void Geant4ParticlePrint::printParticle(const std::string& prefix, Geant4Particl
status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.'
);
+ if ( e && m_printHits ) {
+ Geant4ParticleMap* truth = context()->event().extension<Geant4ParticleMap>();
+ G4HCofThisEvent* hc = e->GetHCofThisEvent();
+ for (int ihc=0, nhc=hc->GetNumberOfCollections(); ihc<nhc; ++ihc) {
+ G4VHitsCollection* c = hc->GetHC(ihc);
+ Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(c);
+ size_t nhits = coll->GetSize();
+ for(size_t i=0; i<nhits; ++i) {
+ Geant4HitData* h = coll->hit(i);
+ Geant4Tracker::Hit* trk_hit = dynamic_cast<Geant4Tracker::Hit*>(h);
+ if ( 0 != trk_hit ) {
+ Geant4HitData::Contribution& t = trk_hit->truth;
+ int trackID = t.trackID;
+ int trueID = truth->particleID(trackID);
+ if ( trueID == p->id ) {
+ print("+++ %20s %s[%d] (%+.2e,%+.2e,%+.2e)[mm]","",c->GetName().c_str(),i,
+ trk_hit->position.x(),trk_hit->position.y(),trk_hit->position.z());
+ }
+ }
+ Geant4Calorimeter::Hit* cal_hit = dynamic_cast<Geant4Calorimeter::Hit*>(h);
+ if ( 0 != cal_hit ) {
+ Geant4HitData::Contributions& contrib = cal_hit->truth;
+ for(Geant4HitData::Contributions::iterator j=contrib.begin(); j!=contrib.end(); ++j) {
+ Geant4HitData::Contribution& t = *j;
+ int trackID = t.trackID;
+ int trueID = truth->particleID(trackID);
+ if ( trueID == p->id ) {
+ print("+++ %20s %s[%d] (%+.2e,%+.2e,%+.2e)[mm]","",c->GetName().c_str(),i,
+ cal_hit->position.x(),cal_hit->position.y(),cal_hit->position.z());
+ }
+ }
+ }
+ }
+ }
+ }
}
/// Print record of kept particles
-void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const {
+void Geant4ParticlePrint::printParticles(const G4Event* e, const ParticleMap& particles) const {
int num_energy = 0;
int num_parent = 0;
int num_process = 0;
@@ -119,7 +156,7 @@ void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const {
for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->reason);
- printParticle("MC Particle Track",p);
+ printParticle("MC Particle Track",e, p);
num_secondaries += int(p->daughters.size());
if ( mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD) ) ++num_energy;
if ( mask.isSet(G4PARTICLE_PRIMARY) ) ++num_primary;
@@ -136,7 +173,7 @@ void Geant4ParticlePrint::printParticles(const ParticleMap& particles) const {
num_calo_hits,num_tracker_hits,num_process,num_parent);
}
-void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int level, Geant4ParticleHandle p) const {
+void Geant4ParticlePrint::printParticleTree(const G4Event* e, const ParticleMap& particles, int level, Geant4ParticleHandle p) const {
char txt[32];
size_t len = sizeof(txt)-1;
// Ensure we do not overwrite the array
@@ -149,18 +186,18 @@ void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles, int le
txt[level+6]='+';
::memset(txt+level+6+1,'-',len-level-3-6);
- printParticle(txt, p);
+ printParticle(txt, e, p);
const set<int>& daughters = p->daughters;
// For all particles, the set of daughters must be contained in the record.
for(set<int>::const_iterator id=daughters.begin(); id!=daughters.end(); ++id) {
int id_dau = *id;
Geant4ParticleHandle dau = (*particles.find(id_dau)).second;
- printParticleTree(particles,level+1,dau);
+ printParticleTree(e, particles, level+1, dau);
}
}
/// Print tree of kept particles
-void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles) const {
+void Geant4ParticlePrint::printParticleTree(const G4Event* e, const ParticleMap& particles) const {
print("+++ MC Particle Parent daughter relationships. [%d particles]",int(particles.size()));
print("+++ MC Particles %12s #Tracks:%7d %-12s Parent%-7s "
"Primary Secondary Energy %-8s Calo Tracker Process/Par Details",
@@ -168,6 +205,6 @@ void Geant4ParticlePrint::printParticleTree(const ParticleMap& particles) const
for(ParticleMap::const_iterator i=particles.begin(); i!=particles.end(); ++i) {
Geant4ParticleHandle p = (*i).second;
PropertyMask mask(p->reason);
- if ( mask.isSet(G4PARTICLE_PRIMARY) ) printParticleTree(particles,0,p);
+ if ( mask.isSet(G4PARTICLE_PRIMARY) ) printParticleTree(e, particles, 0, p);
}
}
diff --git a/DDG4/src/Geant4PrimaryHandler.cpp b/DDG4/src/Geant4PrimaryHandler.cpp
index bdbd15065825e1ed58b493cac6c485a6b861ee3e..687aed5d7e58635d6d715d2c0e105c457dc9b9d8 100644
--- a/DDG4/src/Geant4PrimaryHandler.cpp
+++ b/DDG4/src/Geant4PrimaryHandler.cpp
@@ -9,22 +9,11 @@
// Framework include files
#include "DD4hep/Printout.h"
#include "DD4hep/InstanceCount.h"
-#include "DDG4/Geant4Primary.h"
+#include "DDG4/Geant4InputHandling.h"
#include "DDG4/Geant4PrimaryHandler.h"
-#include "G4Event.hh"
-#include "G4PrimaryVertex.hh"
-#include "G4PrimaryParticle.hh"
-
-#include <stdexcept>
-
-using namespace std;
-using namespace DD4hep;
using namespace DD4hep::Simulation;
-typedef Geant4PrimaryInteraction Interaction;
-typedef ReferenceBitMask<int> PropertyMask;
-
/// Standard constructor
Geant4PrimaryHandler::Geant4PrimaryHandler(Geant4Context* context, const std::string& nam)
: Geant4GeneratorAction(context,nam)
@@ -37,109 +26,7 @@ Geant4PrimaryHandler::~Geant4PrimaryHandler() {
InstanceCount::decrement(this);
}
-static G4PrimaryParticle* createG4Primary(const Geant4ParticleHandle p) {
- G4PrimaryParticle* g4 = new G4PrimaryParticle(p->pdgID, p->psx, p->psy, p->psz);
- g4->SetMass(p->mass);
- return g4;
-}
-
-typedef map<Geant4Particle*,G4PrimaryParticle*> Primaries;
-Primaries getRelevant(set<int>& visited,
- map<int,G4PrimaryParticle*>& prim,
- Interaction::ParticleMap& pm,
- const Geant4ParticleHandle p)
-{
- Primaries res;
- visited.insert(p->id);
- PropertyMask status(p->status);
- if ( status.isSet(G4PARTICLE_GEN_STABLE) ) {
- if ( prim.find(p->id) == prim.end() ) {
- G4PrimaryParticle* p4 = createG4Primary(p);
- prim[p->id] = p4;
- res.insert(make_pair(p,p4));
- }
- }
- else if ( p->daughters.size() > 0 ) {
- const Geant4Particle::Particles& dau = p->daughters;
- int first_daughter = *(dau.begin());
- Geant4ParticleHandle dp = pm[first_daughter];
- double en = p.energy();
- double me = en > std::numeric_limits<double>::epsilon() ? p->mass / en : 0.0;
- // fix by S.Morozov for real != 0
- double proper_time = fabs(dp->time-p->time) * me;
- double proper_time_Precision = pow(10.,-DBL_DIG)*me*fmax(fabs(p->time),fabs(dp->time));
- bool isProperTimeZero = ( proper_time <= proper_time_Precision ) ;
- // -- remove original --- if (proper_time != 0) {
- if ( !isProperTimeZero ) {
- map<int,G4PrimaryParticle*>::iterator ip4 = prim.find(p->id);
- G4PrimaryParticle* p4 = (ip4 == prim.end()) ? 0 : (*ip4).second;
- if ( !p4 ) {
- p4 = createG4Primary(p);
- p4->SetProperTime(proper_time);
- prim[p->id] = p4;
- Primaries daughters;
- for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
- if ( visited.find(*i) == visited.end() ) {
- Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
- daughters.insert(tmp.begin(),tmp.end());
- }
- }
- for(Primaries::iterator i=daughters.begin(); i!=daughters.end(); ++i)
- p4->SetDaughter((*i).second);
- }
- res.insert(make_pair(p,p4));
- }
- else {
- for(Geant4Particle::Particles::const_iterator i=dau.begin(); i!=dau.end(); ++i) {
- if ( visited.find(*i) == visited.end() ) {
- Primaries tmp = getRelevant(visited,prim,pm,pm[*i]);
- res.insert(tmp.begin(),tmp.end());
- }
- }
- }
- }
- return res;
-}
-
/// Event generation action callback
void Geant4PrimaryHandler::operator()(G4Event* event) {
- Geant4PrimaryMap* primaries = context()->event().extension<Geant4PrimaryMap>();
- Interaction* interaction = context()->event().extension<Interaction>();
- Interaction::ParticleMap& pm = interaction->particles;
- Interaction::VertexMap& vm = interaction->vertices;
- map<int,G4PrimaryParticle*> prim;
- set<int> visited;
-
- Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
- for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i) {
- int num_part = 0;
- Geant4Vertex* v = (*i).second;
- G4PrimaryVertex* v4 = new G4PrimaryVertex(v->x,v->y,v->z,v->time);
- event->AddPrimaryVertex(v4);
- print("+++++ G4PrimaryVertex at (%+.2e,%+.2e,%+.2e) [mm] %+.2e [ns]",
- v->x/mm,v->y/mm,v->z/mm,v->time/ns);
- for(Geant4Vertex::Particles::const_iterator ip=v->out.begin(); ip!=v->out.end(); ++ip) {
- Geant4ParticleHandle p = pm[*ip];
- if ( p->parents.size() == 0 ) {
- Primaries relevant = getRelevant(visited,prim,pm,p);
- for(Primaries::const_iterator j=relevant.begin(); j!= relevant.end(); ++j) {
- Geant4ParticleHandle r = (*j).first;
- G4PrimaryParticle* p4 = (*j).second;
- PropertyMask reason(r->reason);
- reason.set(G4PARTICLE_PRIMARY);
- v4->SetPrimary(p4);
- printM1("+++ +-> G4Primary[%3d] ID:%3d type:%9d/%-12s "
- "Momentum:(%+.2e,%+.2e,%+.2e) [GeV] time:%+.2e [ns] #Par:%3d #Dau:%3d",
- num_part,r->id,r->pdgID,r.particleName().c_str(),
- r->psx/GeV,r->psy/GeV,r->psz/GeV,r->time/ns,
- int(r->parents.size()),int(r->daughters.size()));
- ++num_part;
- }
- }
- }
- }
- for(map<int,G4PrimaryParticle*>::iterator i=prim.begin(); i!=prim.end(); ++i) {
- Geant4ParticleHandle p = pm[(*i).first];
- primaries->primaryMap[(*i).second] = p->addRef();
- }
+ generatePrimaries(this, context(), event);
}
diff --git a/DDG4/src/Geant4UIManager.cpp b/DDG4/src/Geant4UIManager.cpp
index 459626226c8322aae76b43588a25e0eb120edc11..97ebd1d9b0900587de40a854cac3cef4754f87a4 100644
--- a/DDG4/src/Geant4UIManager.cpp
+++ b/DDG4/src/Geant4UIManager.cpp
@@ -108,6 +108,6 @@ void Geant4UIManager::start() {
/// Stop and release resources
void Geant4UIManager::stop() {
- // deletePtr(m_vis);
- //deletePtr(m_ui);
+ deletePtr(m_vis);
+ deletePtr(m_ui);
}