diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py index 43fb731097025b90c516f90cb3e55b67575f5b43..9cfb759830ca869f53e3e0eec25b668c9fa6a78a 100644 --- a/DDG4/examples/CLICSidSimu.py +++ b/DDG4/examples/CLICSidSimu.py @@ -14,8 +14,9 @@ from SystemOfUnits import * """ def run(): kernel = DDG4.Kernel() + kernel.UI = "UI" kernel.loadGeometry("file:../DD4hep.trunk/DDExamples/CLICSiD/compact/compact.xml") - kernel.loadXML("DDG4_field.xml") + kernel.loadXML("file:../DD4hep.trunk/DDG4/examples/DDG4_field.xml") lcdd = kernel.lcdd() print '+++ List of sensitive detectors:' @@ -25,6 +26,14 @@ def run(): if sd.isValid(): print '+++ %-32s type:%s'%(o.name(), sd.type(), ) + # Configure UI + ui = DDG4.Action(kernel,"Geant4UIManager/UI") + ui.HaveVIS = True + ui.HaveUI = True + ui.SessionType = 'csh' + kernel.registerGlobalAction(ui) + + # Configure Run actions run1 = DDG4.RunAction(kernel,'Geant4TestRunAction/RunInit') run1.Property_int = 12345 @@ -52,6 +61,13 @@ def run(): kernel.eventAction().add(evt1) kernel.eventAction().add(evt2) + trk = DDG4.Action(kernel,"Geant4TrackPersistency/MonteCarloTruthHandler") + mc = DDG4.Action(kernel,"Geant4MonteCarloRecordManager/MonteCarloRecordManager") + kernel.registerGlobalAction(trk) + kernel.registerGlobalAction(mc) + trk.release() + mc.release() + # Configure I/O evt_root = DDG4.EventAction(kernel,'Geant4Output2ROOT/RootOutput') evt_root.Control = True @@ -70,8 +86,18 @@ def run(): gun.energy = 0.5*GeV gun.particle = 'e-' gun.multiplicity = 1 + gun.position = (0.15*mm,0.12*mm,0.1*cm) gun.enableUI() kernel.generatorAction().add(gun) + """ + rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader") + rdr.zSpread = 0.0 + rdr.lorentzAngle = 0.0 + rdr.OutputLevel = DDG4.OutputLevel.INFO + rdr.Input = "LcioEventReader|test.data" + rdr.enableUI() + kernel.generatorAction().add(rdr) + """ # Setup global filters fur use in sensntive detectors f1 = DDG4.Filter(kernel,'GeantinoRejectFilter/GeantinoRejector') diff --git a/DDG4/examples/initAClick.C b/DDG4/examples/initAClick.C index d5c5691cdac7a0b0d62460fa7d99e244484a4a7d..f1b5f674d6f7c47ceafaa42b8793b21609593d17 100644 --- a/DDG4/examples/initAClick.C +++ b/DDG4/examples/initAClick.C @@ -9,10 +9,12 @@ string make_str(const char* data) { } void initAClick() { + string rootsys = make_str(gSystem->Getenv("ROOTSYS")); string g4_base = make_str(gSystem->Getenv("Geant4_DIR")); string dd4hep = make_str(gSystem->Getenv("DD4hep_DIR")); string inc = " -I"+dd4hep+"/include -I"+g4_base+"/include/Geant4 -Wno-shadow -g -O0"; - string libs = " -L"+dd4hep+"/lib -lDD4hepCore -lDD4hepG4 -lDDSegmentation"; + string libs = (" -L"+rootsys+"/lib"); + libs += " -lCore -lCint -lMathCore -L"+dd4hep+"/lib -lDD4hepCore -lDD4hepG4 -lDDSegmentation"; libs += (" -L"+g4_base+"/lib -L"+g4_base+"/lib64 -lG4event -lG4tracking -lG4particles"); gSystem->AddIncludePath(inc.c_str()); gSystem->AddLinkedLibs(libs.c_str()); diff --git a/DDG4/examples/sequences.xml b/DDG4/examples/sequences.xml index 12977fe0f0da1e642ee4eb7ca8da29787ca77f39..909624ab8899c755d3444784578a542a0145127f 100644 --- a/DDG4/examples/sequences.xml +++ b/DDG4/examples/sequences.xml @@ -65,7 +65,7 @@ <properties particle="opticalphoton"/> </filter> <filter name="EnergyDepositMinimumCut"> - <properties Cut="10"/> + <properties Cut="10*keV"/> </filter> </filters> @@ -92,7 +92,7 @@ </sequence> <sequence name="Geant4GeneratorActionSequence/GeneratorAction"> <action name="Geant4ParticleGun/Gun"> - <properties energy="500" + <properties energy="500*GeV" particle="'e-'" multiplicity="1"/> </action> diff --git a/DDG4/include/DDG4/ComponentProperties_inl.h b/DDG4/include/DDG4/ComponentProperties_inl.h index e15fbfa0e9536918f60595807d3f8216e9cdfc8d..a31d8afcbd400b73ccfdb7e0861c3b293167dce7 100644 --- a/DDG4/include/DDG4/ComponentProperties_inl.h +++ b/DDG4/include/DDG4/ComponentProperties_inl.h @@ -49,6 +49,10 @@ namespace DD4hep { virtual std::string str(const void* ptr) const; /// PropertyGrammar overload: Retrieve value from string virtual bool fromString(void* ptr, const std::string& value) const; + /// Evaluate string value if possible before calling boost::spirit + virtual int evaluate(void* ptr, const std::string& value) const; + /// Pre-parse string + virtual std::string pre_parse(const std::string& in) const; }; /// Standarsd constructor @@ -69,24 +73,40 @@ namespace DD4hep { template <typename TYPE> const std::type_info& Grammar<TYPE>::type() const { return typeid(TYPE); } + + /// Evaluate string value if possible before calling boost::spirit + template <typename TYPE> int Grammar<TYPE>::evaluate(void*, const std::string&) const { + return 0; + } + /// Pre-parse string + template <typename TYPE> std::string Grammar<TYPE>::pre_parse(const std::string& in) const{ + return in; + } + /// PropertyGrammar overload: Retrieve value from string template <typename TYPE> bool Grammar<TYPE>::fromString(void* ptr, const std::string& str) const { -#ifdef DD4HEP_USE_BOOST + int sc = 0; TYPE temp; - int sc = Parsers::parse(temp,str); - //std::cout << "Converting value: " << str << " to type " << typeid(TYPE).name() << std::endl; - if ( sc ) { +#ifdef DD4HEP_USE_BOOST + sc = Parsers::parse(temp,str); +#endif + if ( !sc ) sc = evaluate(&temp,str); +#if 0 + std::cout << "Sc=" << sc << " Converting value: " << str + << " to type " << typeid(TYPE).name() + << std::endl; +#endif + if ( sc ) { *(TYPE*)ptr = temp; return true; } +#ifndef DD4HEP_USE_BOOST + throw std::runtime_error("This version of DD4HEP is not compiled to use boost::spirit.\n" + "To enable elaborated property handling set DD4HEP_USE_BOOST=ON\n" + "and BOOST_INCLUDE_DIR=<boost include path>"); +#else PropertyGrammar::invalidConversion(str, typeid(TYPE)); return false; -#else - if (!ptr || str.length() == 0) { - } - throw std::runtime_error("This version of DD4HEP is not compiled to use boost::spirit.\n" - "To enable elaborated property handling set DD4HEP_USE_BOOST=ON\n" - "and BOOST_INCLUDE_DIR=<boost include path>"); #endif } diff --git a/DDG4/include/DDG4/Factories.h b/DDG4/include/DDG4/Factories.h index fc70eea44e05bec23ecdb1c0015fe3534f02524b..e17abb912f857b5991f7866c71aba9ae894f85c1 100644 --- a/DDG4/include/DDG4/Factories.h +++ b/DDG4/include/DDG4/Factories.h @@ -173,12 +173,16 @@ namespace { DECLARE_EXTERNAL_GEANT4SENSITIVEDETECTOR(name,DD4hep::Simulation::__sd_create__##name) #endif -#define DECLARE_GEANT4SENSITIVE(name) namespace DD4hep { namespace Simulation { }} using DD4hep::Simulation::name; \ +#define DECLARE_GEANT4SENSITIVE_NS(ns,name) using ns::name; \ PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::Simulation::Geant4Sensitive*(DD4hep::Simulation::Geant4Context*,std::string,DD4hep::Geometry::DetElement*,DD4hep::Geometry::LCDD*)) +#define DECLARE_GEANT4SENSITIVE(name) DECLARE_GEANT4SENSITIVE_NS(DD4hep::Simulation,name) /// Plugin defintion to create Geant4Action objects -#define DECLARE_GEANT4ACTION(name) namespace DD4hep { namespace Simulation { }} using DD4hep::Simulation::name; \ +#define DECLARE_GEANT4ACTION_NS(ns,name) using ns::name; \ PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::Simulation::Geant4Action*(DD4hep::Simulation::Geant4Context*,std::string)) +/// Plugin defintion to create Geant4Action objects +#define DECLARE_GEANT4ACTION(name) DECLARE_GEANT4ACTION_NS(DD4hep::Simulation,name) + /// Plugin definition to create Geant4 stepper objects #define DECLARE_GEANT4_STEPPER(name) PLUGINSVC_FACTORY_WITH_ID(G4##name,std::string(#name),G4MagIntegratorStepper*(G4EquationOfMotion*)) #define DECLARE_GEANT4_MAGSTEPPER(name) PLUGINSVC_FACTORY_WITH_ID(G4##name,std::string(#name),G4MagIntegratorStepper*(G4Mag_EqRhs*)) diff --git a/DDG4/include/DDG4/Geant4Action.h b/DDG4/include/DDG4/Geant4Action.h index dd7db86d6344a07bc571142b754a630dda16a3a3..1f358dafb945cc5792a5e6b81db13f9741e25dd4 100644 --- a/DDG4/include/DDG4/Geant4Action.h +++ b/DDG4/include/DDG4/Geant4Action.h @@ -62,7 +62,10 @@ namespace DD4hep { TypeName(const std::string& typ, const std::string& nam) : std::pair<std::string, std::string>(typ, nam) { } + /// Split string pair according to default delimiter ('/') static TypeName split(const std::string& type_name); + /// Split string pair according to custom delimiter + static TypeName split(const std::string& type_name, const std::string& delim); }; /** @class Geant4TrackInformation Geant4Action.h DDG4/Geant4Action.h @@ -288,7 +291,7 @@ namespace DD4hep { /// Access to the UI messenger Geant4UIMessenger* control() const; /// Enable and install UI messenger - void enableUI(); + virtual void enableUI(); /// Declare property template <typename T> Geant4Action& declareProperty(const std::string& nam, T& val); /// Declare property @@ -322,6 +325,9 @@ namespace DD4hep { /// Support of exceptions: Print fatal message and throw runtime_error. void except(const std::string& fmt, ...) const; + /// Abort Geant4 Run by throwing a G4Exception with type RunMustBeAborted + void abortRun(const std::string& exception, const std::string& fmt, ...) const; + /// Access to the main run action sequence from the kernel object Geant4RunActionSequence& runAction() const; /// Access to the main event action sequence from the kernel object @@ -334,6 +340,10 @@ namespace DD4hep { Geant4StackingActionSequence& stackingAction() const; /// Access to the main generator action sequence from the kernel object Geant4GeneratorActionSequence& generatorAction() const; + /// Access to the Track Persistency Manager from the kernel object + Geant4MonteCarloTruth& mcTruthMgr() const; + /// Access to the MC record manager from the kernel object + Geant4MonteCarloRecordManager& mcRecordMgr() const; }; /// Declare property diff --git a/DDG4/include/DDG4/Geant4Call.h b/DDG4/include/DDG4/Geant4Call.h new file mode 100644 index 0000000000000000000000000000000000000000..1f517ad846918696116416a0d389836828373277 --- /dev/null +++ b/DDG4/include/DDG4/Geant4Call.h @@ -0,0 +1,39 @@ +// $Id: Geant4Field.cpp 875 2013-11-04 16:15:14Z markus.frank@cern.ch $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4CALL_H +#define DD4HEP_DDG4_GEANT4CALL_H + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + /** @class Geant4Call Geant4Call.h DDG4/Geant4Call.h + * + * Callback interface class with argument + * + * @author M.Frank + * @version 1.0 + */ + class Geant4Call { + public: + /// Default destructor + virtual ~Geant4Call(); + /// Default callback with argument + virtual void operator()(void* param) = 0; + }; + + } // End namespace Simulation +} // End namespace DD4hep +#endif // DD4HEP_DDG4_GEANT4CALL_H diff --git a/DDG4/include/DDG4/Geant4Callback.h b/DDG4/include/DDG4/Geant4Callback.h index d8d0af242fb35b3942e429a1ff4f78cee6ccc2f0..218268c405801cfc77f5d03aeb52a7380aee49b8 100644 --- a/DDG4/include/DDG4/Geant4Callback.h +++ b/DDG4/include/DDG4/Geant4Callback.h @@ -282,12 +282,21 @@ namespace DD4hep { struct CallbackSequence { typedef std::vector<Callback> Callbacks; + enum Location { FRONT, END }; Callbacks callbacks; + /// Default constructor CallbackSequence() { } + /// Copy constructor CallbackSequence(const CallbackSequence& c) : callbacks(c.callbacks) { } + /// Assignment operator + CallbackSequence& operator=(const CallbackSequence& c) { + if ( this != & c ) callbacks = c.callbacks; + return *this; + } + //template <typename TYPE, typename R, typename OBJECT> // CallbackSequence(const std::vector<TYPE*>& objects, R (TYPE::value_type::*pmf)()) { //} @@ -297,8 +306,11 @@ namespace DD4hep { void clear() { callbacks.clear(); } - void add(const Callback& cb) { - callbacks.push_back(cb); + void add(const Callback& cb,Location where) { + if ( where == CallbackSequence::FRONT ) + callbacks.insert(callbacks.begin(),cb); + else + callbacks.insert(callbacks.end(),cb); } void operator()() const; template <typename A0> void operator()(A0 a0) const; @@ -308,19 +320,19 @@ namespace DD4hep { static void checkTypes(const std::type_info& typ1, const std::type_info& typ2, void* test); template <typename TYPE, typename R, typename OBJECT> - void add(TYPE* pointer, R (OBJECT::*pmf)()) { + void add(TYPE* pointer, R (OBJECT::*pmf)(),Location where=CallbackSequence::END) { checkTypes(typeid(TYPE), typeid(OBJECT), dynamic_cast<OBJECT*>(pointer)); - add(Callback(pointer).make(pmf)); + add(Callback(pointer).make(pmf),where); } template <typename TYPE, typename R, typename OBJECT, typename A> - void add(TYPE* pointer, R (OBJECT::*pmf)(A)) { + void add(TYPE* pointer, R (OBJECT::*pmf)(A),Location where=CallbackSequence::END) { checkTypes(typeid(TYPE), typeid(OBJECT), dynamic_cast<OBJECT*>(pointer)); - add(Callback(pointer).make(pmf)); + add(Callback(pointer).make(pmf),where); } template <typename TYPE, typename R, typename OBJECT, typename A1, typename A2> - void add(TYPE* pointer, R (OBJECT::*pmf)(A1, A2)) { + void add(TYPE* pointer, R (OBJECT::*pmf)(A1, A2),Location where=CallbackSequence::END) { checkTypes(typeid(TYPE), typeid(OBJECT), dynamic_cast<OBJECT*>(pointer)); - add(Callback(pointer).make(pmf)); + add(Callback(pointer).make(pmf),where); } }; diff --git a/DDG4/include/DDG4/Geant4Context.h b/DDG4/include/DDG4/Geant4Context.h index 47f3291149df166e3c0ee2287e48f118aa9817d6..320e587b5b793a4bd7feee58a00ea5c9b99212fc 100644 --- a/DDG4/include/DDG4/Geant4Context.h +++ b/DDG4/include/DDG4/Geant4Context.h @@ -38,6 +38,8 @@ namespace DD4hep { class Geant4StackingActionSequence; class Geant4GeneratorActionSequence; class Geant4SensDetSequences; + class Geant4MonteCarloTruth; + class Geant4MonteCarloRecordManager; /** @class Geant4Context Geant4Context.h DDG4/Geant4Context.h * @@ -48,22 +50,18 @@ namespace DD4hep { typedef Geometry::LCDD LCDD; Geant4Kernel* m_kernel; - /// Reference to Geant4 track manager - G4TrackingManager* m_trackMgr; /// Default constructor Geant4Context(Geant4Kernel* kernel); /// Default destructor virtual ~Geant4Context(); - /// Access the tracking manager - G4TrackingManager* trackMgr() const { - return m_trackMgr; - } /// Access to the kernel object Geant4Kernel& kernel() { return *m_kernel; } /// Access to detector description LCDD& lcdd() const; + /// Access the tracking manager + G4TrackingManager* trackMgr() const; /// Create a user trajectory virtual G4VTrajectory* createTrajectory(const G4Track* track) const; /// Access to the main run action sequence from the kernel object @@ -80,6 +78,12 @@ namespace DD4hep { Geant4GeneratorActionSequence& generatorAction() const; /// Access to the sensitive detector sequences from the kernel object Geant4SensDetSequences& sensitiveActions() const; + + /// Access to the Track Manager from the kernel object + Geant4MonteCarloTruth& mcTruthMgr() const; + /// Access to the MC record manager from the kernel object (if instantiated!) + Geant4MonteCarloRecordManager& mcRecordMgr() const; + }; } // End namespace Simulation diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h index 893b4cc553fe43c3bbc38cbe4d92f1856cdbc867..fe43cf7760fac1589479dee51116ac0205fd65aa 100644 --- a/DDG4/include/DDG4/Geant4Data.h +++ b/DDG4/include/DDG4/Geant4Data.h @@ -134,7 +134,17 @@ namespace DD4hep { MonteCarloContrib(const MonteCarloContrib& c) : trackID(c.trackID), pdgID(c.pdgID), deposit(c.deposit), time(c.time) { } - // Clear data content + /// Assignment operator + MonteCarloContrib& operator=(const MonteCarloContrib& c) { + if ( this != &c ) { + trackID = c.trackID; + pdgID = c.pdgID; + deposit = c.deposit; + time = c.time; + } + return *this; + } + /// Clear data content void clear() { time = deposit = 0.0; pdgID = trackID = -1; diff --git a/DDG4/include/DDG4/Geant4EventAction.h b/DDG4/include/DDG4/Geant4EventAction.h index c379c62b66435837afb981a0dd2ca1be7a4bf3fa..80544246f639a4a20fb117d3e0eb1a5e8985a913 100644 --- a/DDG4/include/DDG4/Geant4EventAction.h +++ b/DDG4/include/DDG4/Geant4EventAction.h @@ -65,6 +65,8 @@ namespace DD4hep { CallbackSequence m_begin; /// Callback sequence for event finalization action CallbackSequence m_end; + /// Callback sequence for event finalization action + CallbackSequence m_final; /// The list of action objects to be called Actors<Geant4EventAction> m_actors; public: @@ -82,6 +84,11 @@ namespace DD4hep { void callAtEnd(Q* p, void (T::*f)(const G4Event*)) { m_end.add(p, f); } + /// Register event-cleanup callback (after end-of-event callback -- unordered) + template <typename Q, typename T> + void callAtFinal(Q* p, void (T::*f)(const G4Event*)) { + m_final.add(p, f); + } /// Add an actor responding to all callbacks. Sequence takes ownership. void adopt(Geant4EventAction* action); /// Begin-of-event callback diff --git a/DDG4/include/DDG4/Geant4HitCollection.h b/DDG4/include/DDG4/Geant4HitCollection.h index b64d570466850292152c6ad229b33cf3e17ec839..f48db4dd24b51ee095991399e7f1bda2d30e9f1d 100644 --- a/DDG4/include/DDG4/Geant4HitCollection.h +++ b/DDG4/include/DDG4/Geant4HitCollection.h @@ -98,18 +98,18 @@ namespace DD4hep { public: /// Default constructor - Geant4HitWrapper() { + Geant4HitWrapper() : G4VHit() { m_data.second = manipulator<InvalidHit>(); m_data.first = 0; } /// Copy constructor - Geant4HitWrapper(const Geant4HitWrapper& v) { + Geant4HitWrapper(const Geant4HitWrapper& v) : G4VHit() { m_data = v.m_data; v.m_data.first = 0; //v.release(); } /// Copy constructor - Geant4HitWrapper(const Wrapper& v) { + Geant4HitWrapper(const Wrapper& v) : G4VHit() { m_data = v; } /// Default destructor @@ -140,9 +140,11 @@ namespace DD4hep { } /// Assignment transfers the pointer ownership Geant4HitWrapper& operator=(const Geant4HitWrapper& v) { - m_data = v.m_data; - v.m_data.first = 0; - //v.release(); + if ( this != & v ) { + m_data = v.m_data; + v.m_data.first = 0; + //v.release(); + } return *this; } /// Automatic conversion to the desired type diff --git a/DDG4/include/DDG4/Geant4Hits.h b/DDG4/include/DDG4/Geant4Hits.h index 7a0c9f8833daf60204421182e36f4178dbaeec29..b3512c6cfd3262f3710bccc71189561e2014238a 100644 --- a/DDG4/include/DDG4/Geant4Hits.h +++ b/DDG4/include/DDG4/Geant4Hits.h @@ -102,6 +102,16 @@ namespace DD4hep { MonteCarloContrib(const MonteCarloContrib& c) : trackID(c.trackID), pdgID(c.pdgID), deposit(c.deposit), time(c.time) { } + /// Assignment operator + MonteCarloContrib& operator=(const MonteCarloContrib& c) { + if ( this != &c ) { + trackID = c.trackID; + pdgID = c.pdgID; + deposit = c.deposit; + time = c.time; + } + return *this; + } void clear() { time = deposit = 0.0; pdgID = trackID = -1; @@ -112,7 +122,7 @@ namespace DD4hep { public: /// Standard constructor - Geant4Hit() { + Geant4Hit() : cellID(0) { } /// Default destructor virtual ~Geant4Hit() { diff --git a/DDG4/include/DDG4/Geant4Kernel.h b/DDG4/include/DDG4/Geant4Kernel.h index 0c7c4ba6d6c745d88f13136e6d9d22324ae25b18..3fe94f2a94d24bee764eda54bdc295652f79af23 100644 --- a/DDG4/include/DDG4/Geant4Kernel.h +++ b/DDG4/include/DDG4/Geant4Kernel.h @@ -52,6 +52,8 @@ namespace DD4hep { class Geant4PhysicsListActionSequence; class Geant4SensDetActionSequence; class Geant4SensDetSequences; + class Geant4MonteCarloTruth; + class Geant4MonteCarloRecordManager; /** @class Invoke Geant4Kernel.h DDG4/Geant4Kernel.h * @@ -91,6 +93,13 @@ namespace DD4hep { Geant4SensDetSequences* m_sensDetActions; /// Reference to the geant4 physics list Geant4PhysicsListActionSequence* m_physicsList; + /// Reference to track persistency manager + Geant4MonteCarloTruth* m_mcTruthMgr; + /// Reference to MC record manager + Geant4MonteCarloRecordManager* m_mcRecordMgr; + /// Reference to Geant4 track manager + G4TrackingManager* m_trackMgr; + /// Action phases Phases m_phases; /// Globally registered actions @@ -99,8 +108,12 @@ namespace DD4hep { GlobalActions m_globalFilters; /// Detector description object LCDD& m_lcdd; - /// Name of the G4UI command tree + /// Property: Name of the G4UI command tree std::string m_controlName; + /// Property: Name of the UI action. Must be member of the global actions + std::string m_uiName; + /// Property: Number of events to be executed in batch mode + long m_numEvent; /// Helper to register an action sequence template <typename C> bool registerSequence(C*& seq, const std::string& name); @@ -122,6 +135,8 @@ namespace DD4hep { PhaseSelector(Geant4Kernel* kernel); /// Copy constructor PhaseSelector(const PhaseSelector& c); + /// Assignment operator + PhaseSelector& operator=(const PhaseSelector& c); /// Phase access to the map Geant4ActionPhase& operator[](const std::string& name) const; } phase; @@ -157,12 +172,28 @@ namespace DD4hep { LCDD& lcdd() const { return m_lcdd; } + /// Access the tracking manager + G4TrackingManager* trackMgr() const { + return m_trackMgr; + } + /// Access the tracking manager + void setTrackMgr(G4TrackingManager* mgr) { + m_trackMgr = mgr; + } /// Access to the Geant4 run manager G4RunManager& runManager(); /// Access the command directory const std::string& directoryName() const { return m_controlName; } + /// Declare property + template <typename T> Geant4Kernel& declareProperty(const std::string& nam, T& val); + /// Declare property + template <typename T> Geant4Kernel& declareProperty(const char* nam, T& val); + /// Check property for existence + bool hasProperty(const std::string& name) const; + /// Access single property + Property& property(const std::string& name); /// Register action by name to be retrieved when setting up and connecting action objects /** Note: registered actions MUST be unique. @@ -258,6 +289,10 @@ namespace DD4hep { Geant4PhysicsListActionSequence& physicsList() { return *physicsList(true); } + /// Access to the Track Manager from the kernel object + Geant4MonteCarloTruth& mcTruthMgr(); + /// Access to the MC record manager from the kernel object (if instantiated!) + Geant4MonteCarloRecordManager& mcRecordMgr(); /// Construct detector geometry using lcdd plugin void loadGeometry(const std::string& compact_file); @@ -268,7 +303,18 @@ namespace DD4hep { void terminate(); void loadXML(const char* fname); }; + /// Declare property + template <typename T> Geant4Kernel& Geant4Kernel::declareProperty(const std::string& nam, T& val) { + m_properties.add(nam, val); + return *this; + } + /// Declare property + template <typename T> Geant4Kernel& Geant4Kernel::declareProperty(const char* nam, T& val) { + m_properties.add(nam, val); + return *this; + } + struct Geant4Exec { static int configure(Geant4Kernel& kernel); static int initialize(Geant4Kernel& kernel); diff --git a/DDG4/include/DDG4/Geant4MonteCarloRecordManager.h b/DDG4/include/DDG4/Geant4MonteCarloRecordManager.h new file mode 100644 index 0000000000000000000000000000000000000000..b328e9b49ceae989af5bcd2e9c10de0b081efa19 --- /dev/null +++ b/DDG4/include/DDG4/Geant4MonteCarloRecordManager.h @@ -0,0 +1,49 @@ +// $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 +#ifndef DD4HEP_DDG4_GEANT4MONTECARLORECORDMANAGER_H +#define DD4HEP_DDG4_GEANT4MONTECARLORECORDMANAGER_H + +// Framework include files +#include "DDG4/Geant4TrackPersistency.h" + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + // Forward declarations + + /** @class Invoke Geant4TrackPersistency.h DDG4/Geant4TrackPersistency.h + * + * Default base class for all geant 4 actions and derivates thereof. + * + * @author M.Frank + * @version 1.0 + */ + class Geant4MonteCarloRecordManager : public Geant4Action { + public: + /// Flag to indicate if the track information should be collected + bool m_collectInfo; + public: + /// Standard constructor + Geant4MonteCarloRecordManager(Geant4Context* context, const std::string& nam); + /// Default destructor + virtual ~Geant4MonteCarloRecordManager(); + /// Save G4Track data + virtual void save(const Geant4TrackPersistency::TrackInfo& track); + }; + } // End namespace Simulation +} // End namespace DD4hep +#endif // DD4HEP_DDG4_GEANT4MONTECARLORECORDMANAGER_H diff --git a/DDG4/include/DDG4/Geant4MonteCarloTruth.h b/DDG4/include/DDG4/Geant4MonteCarloTruth.h new file mode 100644 index 0000000000000000000000000000000000000000..cce5b1296496e994cf73504c5e264dc87afc2e10 --- /dev/null +++ b/DDG4/include/DDG4/Geant4MonteCarloTruth.h @@ -0,0 +1,57 @@ +// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $ +//==================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4MONTECARLOTRUTH_H +#define DD4HEP_DDG4_GEANT4MONTECARLOTRUTH_H + +// Framework include files + +// C/C++ include files + +// Forward declarations +class G4Step; +class G4Track; +class G4Event; + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + /** @class Geant4MonteCarloTruth Geant4MonteCarloTruth.h DDG4/Geant4MonteCarloTruth.h + * + * Default Interface class to handle monte carlo truth records + * + * @author M.Frank + * @version 1.0 + */ + class Geant4MonteCarloTruth { + protected: + /// Standard constructor + Geant4MonteCarloTruth(); + public: + /// Default destructor + virtual ~Geant4MonteCarloTruth(); + /// Mark a Geant4 track to be kept for later MC truth analysis + virtual void mark(const G4Track* track) = 0; + /// Store a track, with a flag + virtual void mark(const G4Track* track, bool created_hit) = 0; + /// Mark a Geant4 track of the step to be kept for later MC truth analysis + virtual void mark(const G4Step* step) = 0; + /// Store a track produced in a step to be kept for later MC truth analysis + virtual void mark(const G4Step* step, bool created_hit) = 0; + }; + } // End namespace Simulation +} // End namespace DD4hep + +#endif // DD4HEP_DDG4_GEANT4MONTECARLOTRUTH_H diff --git a/DDG4/include/DDG4/Geant4ParticleGun.h b/DDG4/include/DDG4/Geant4ParticleGun.h index 7d6fa51fd991c745c820492f3eeaa68fc5e06b1f..4a06591852632be030290b4ecafed6d66139a3d1 100644 --- a/DDG4/include/DDG4/Geant4ParticleGun.h +++ b/DDG4/include/DDG4/Geant4ParticleGun.h @@ -11,6 +11,7 @@ // Framework include files #include "DDG4/Geant4GeneratorAction.h" +#include "Math/Vector3D.h" // Forward declarations class G4ParticleDefinition; @@ -35,9 +36,7 @@ namespace DD4hep { class Geant4ParticleGun: public Geant4GeneratorAction { protected: /// Position and shooting direction of the gun - struct { - double x, y, z; - } m_position, m_direction; + ROOT::Math::XYZVector m_position, m_direction; /// Particle energy double m_energy; /// Desired multiplicity of the particles to be shot diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h index 3c9a640b76268e72339b46351d3731387daaae9b..06c5a2e518fe0672c6c9349d940b86aee0211e71 100644 --- a/DDG4/include/DDG4/Geant4SensDetAction.h +++ b/DDG4/include/DDG4/Geant4SensDetAction.h @@ -152,6 +152,12 @@ namespace DD4hep { return detector().readoutGeometry(); } + /// Mark the track to be kept for MC truth propagation during hit processing + void mark(const G4Track* track) const; + + /// Mark the track of this step to be kept for MC truth propagation during hit processing + void mark(const G4Step* step) const; + /// Access to the hosting sequence Geant4SensDetActionSequence& sequence() const; diff --git a/DDG4/include/DDG4/Geant4TrackHandler.h b/DDG4/include/DDG4/Geant4TrackHandler.h index 6a19d5b944c2c9bb9cd2fa681ac44a0db54a3734..de0dafd72493664957c48adf2703c566af5c6937 100644 --- a/DDG4/include/DDG4/Geant4TrackHandler.h +++ b/DDG4/include/DDG4/Geant4TrackHandler.h @@ -6,8 +6,8 @@ // Author : M.Frank // //==================================================================== -#ifndef DD4HEP_GEANT4TRACKHANDLER_H -#define DD4HEP_GEANT4TRACKHANDLER_H +#ifndef DD4HEP_DDG4_GEANT4TRACKHANDLER_H +#define DD4HEP_DDG4_GEANT4TRACKHANDLER_H // Framework include files #include "DDG4/Defs.h" @@ -17,6 +17,10 @@ #include "G4TrajectoryPoint.hh" #include "G4VTouchable.hh" #include "G4VSensitiveDetector.hh" +#include "G4ParticleDefinition.hh" +#include "G4VProcess.hh" + +#include <stdexcept> // Forward declarations class G4VTouchableHandle; @@ -52,23 +56,67 @@ namespace DD4hep { /// Initializing constructor Geant4TrackHandler(const G4Track* t) : track(t) { + /// Should test here if the track pointer is valid to avoind any later trouble + if ( 0 == t ) { + throw std::runtime_error("Geant4TrackHandler: NULL pointer passed to constructor!"); + } + } + const char* statusName() const { + switch( track->GetTrackStatus() ) { + case fAlive: return "Alive"; + case fStopButAlive: return "StopButAlive"; + case fStopAndKill: return "StopAndKill"; + case fKillTrackAndSecondaries: return "KillTrackAndSecondaries"; + case fSuspend: return "Suspend"; + case fPostponeToNextEvent: return "PostponeToNextEvent"; + default: return "UNKNOWN"; + } } + /// Conversion to G4Track operator const G4Track*() const { return track; } + /// Track's identifier + int id() const { + return track->GetTrackID(); + } + /// Track's parent identifier + int parent() const { + return track->GetParentID(); + } /// Track's particle definition G4ParticleDefinition* trackDef() const { return track->GetDefinition(); } - /// Track position + /// Track's particle name + const std::string& name() const { + return trackDef()->GetParticleName(); + } + /// Track's position const G4ThreeVector& position() const { return track->GetPosition(); } - /// Track energy + /// Track's vertex position, where the track was created + const G4ThreeVector& vertex() const { + return track->GetVertexPosition(); + } + /// Track global time + double globalTime() const { + return track->GetGlobalTime(); + } + /// Track proper time + double properTime() const { + return track->GetProperTime(); + } + /// Track's energy double energy() const { return track->GetTotalEnergy(); } + /// Track's kinetic energy + double kineticEnergy() const { + return track->GetKineticEnergy(); + } /// Track velocity double velocity() const { return track->GetVelocity(); @@ -93,7 +141,6 @@ namespace DD4hep { const G4LogicalVolume* vertexVol() const { return track->GetLogicalVolumeAtVertex(); } - /// Touchable of the track const Touchable& touchable() const { return track->GetTouchableHandle(); @@ -102,17 +149,20 @@ namespace DD4hep { const Touchable& nextTouchable() const { return track->GetNextTouchableHandle(); } - /// Physical process of the track generation const G4VProcess* creatorProcess() const { return track->GetCreatorProcess(); } + /// Physical process of the track generation + const std::string creatorName() const { + const G4VProcess* p = creatorProcess(); + if ( p ) return p->GetProcessName(); + return ""; + } /// User information block Info* userInfo() const { return track->GetUserInformation(); } - /// Set user information block (const mis-match) - //void setUserInfo(Info* info) { track->SetUserInformation(info); } /// Specific user information block template <typename T> T* info() const { return (T*) userInfo(); @@ -125,7 +175,7 @@ namespace DD4hep { G4int stepNumber() const { return track->GetCurrentStepNumber(); } - + /// Access the PDG particle identification int pdgID() const { G4ParticleDefinition* def = trackDef(); return def ? def->GetPDGEncoding() : 0; @@ -134,5 +184,4 @@ namespace DD4hep { } // End namespace Simulation } // End namespace DD4hep - -#endif // DD4HEP_GEANT4HITS_H +#endif // DD4HEP_DDG4_GEANT4TRACKHANDLER_H diff --git a/DDG4/include/DDG4/Geant4TrackManager.h b/DDG4/include/DDG4/Geant4TrackManager.h new file mode 100644 index 0000000000000000000000000000000000000000..872136ceffdef21776cf1f76ea21474a4bb25327 --- /dev/null +++ b/DDG4/include/DDG4/Geant4TrackManager.h @@ -0,0 +1,99 @@ +// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $ +//==================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4TRACKMANAGER_H +#define DD4HEP_DDG4_GEANT4TRACKMANAGER_H + +// Framework include files +#include "DDG4/Geant4Action.h" +#include "DDG4/Geant4TrackingAction.h" +#include "Math/PxPyPzE4D.h" +#include "G4VUserTrackInformation.hh" + +// C/C++ include files +#include <map> +#include <string> +#include <typeinfo> + +// Forward declarations +class G4Track; + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + // Forward declarations + + /** @class Geant4TrackManager Geant4TrackManager.h DDG4/Geant4TrackManager.h + * + * Default base class for all geant 4 actions and derivates thereof. + * + * @author M.Frank + * @version 1.0 + */ + class Geant4TrackManager : public Geant4Action { + public: + typedef ROOT::Math::PxPyPzE4D<double> FourMomentum; + + /** @class Geant4TrackManager::TrackInfo Geant4TrackManager.h DDG4/Geant4TrackManager.h + * + * + * + * @author M.Frank + * @version 1.0 + */ + class TrackInfo : public G4VUserTrackInformation { + public: + /// Pointer to self + Geant4TrackManager* manager; + /// Pointer to the track + const G4Track* track; + /// Flag to store the track + bool store; + /// Initial 4-momentum at the beginning of the tracking action + FourMomentum initialMomentum; + + public: + /// Standard constructor + TrackInfo(); + /// Default destructor + virtual ~TrackInfo(); + /// Clear alkl stored information for this track + void set(const G4Track* trk); + }; + protected: + + /// Information block for current track + Geant4TrackManager::TrackInfo m_current; + + public: + /// Standard constructor + Geant4TrackManager(Geant4Context* context, const std::string& nam); + /// Default destructor + virtual ~Geant4TrackManager(); + /// Access the Geant4 tracking manager. Only use between tracking pre- and post action + G4TrackingManager* trackMgr() const { + return m_context->trackMgr(); + } + /// Store a track + void mark(const G4Track* track); + /// Pre-track action callback + virtual void begin(const G4Track* track); + /// Post-track action callback + virtual void end(const G4Track* track); + }; + } // End namespace Simulation +} // End namespace DD4hep + +#endif // DD4HEP_DDG4_GEANT4TRACKMANAGER_H diff --git a/DDG4/include/DDG4/Geant4TrackPersistency.h b/DDG4/include/DDG4/Geant4TrackPersistency.h new file mode 100644 index 0000000000000000000000000000000000000000..c097905b3b99c5e04bb69f3766c23a14583ca779 --- /dev/null +++ b/DDG4/include/DDG4/Geant4TrackPersistency.h @@ -0,0 +1,127 @@ +// $Id: Geant4Hits.h 513 2013-04-05 14:31:53Z gaede $ +//==================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4TRACKPERSISTENCY_H +#define DD4HEP_DDG4_GEANT4TRACKPERSISTENCY_H + +// Framework include files +#include "DDG4/Geant4Action.h" +#include "DDG4/Geant4MonteCarloTruth.h" +#include "Math/PxPyPzE4D.h" +#include "G4VUserTrackInformation.hh" + +// C/C++ include files +#include <map> +#include <string> +#include <typeinfo> + +// Forward declarations +class G4Step; +class G4Track; +class G4Event; +class G4SteppingManager; + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + // Forward declarations + + /** @class Geant4TrackPersistency Geant4TrackPersistency.h DDG4/Geant4TrackPersistency.h + * + * Default base class for all geant 4 actions and derivates thereof. + * + * @author M.Frank + * @version 1.0 + */ + class Geant4TrackPersistency : public Geant4Action, public Geant4MonteCarloTruth { + public: + typedef ROOT::Math::PxPyPzE4D<double> FourMomentum; + typedef std::map<int,void*> TrackMap; + + /** @class Geant4TrackPersistency::TrackInfo Geant4TrackPersistency.h DDG4/Geant4TrackPersistency.h + * + * + * + * @author M.Frank + * @version 1.0 + */ + class TrackInfo : public G4VUserTrackInformation { + public: + /// Pointer to self + Geant4TrackPersistency* manager; + /// Pointer to the track + const G4Track* track; + /// Pointer to MC tracking history + void* info; + /// Flag to store the track + bool store; + /// Flag to know if this track created a Geant4 hit + bool createdHit; + /// Initial 4-momentum at the beginning of the tracking action + FourMomentum initialMomentum; + double initialEnergy; + + public: + /// Standard constructor + TrackInfo(); + /// Default destructor + virtual ~TrackInfo(); + /// Clear all stored information for this track + void set(const G4Track* trk, void* opt); + }; + protected: + + /// Information block for current track + Geant4TrackPersistency::TrackInfo m_current; + /// Property: Steer printout at tracking action begin + bool m_printStartTracking; + /// Property: Steer printout at tracking action end + bool m_printEndTracking; + + public: + TrackMap m_trackMap; + + public: + /// Standard constructor + Geant4TrackPersistency(Geant4Context* context, const std::string& nam); + /// Default destructor + virtual ~Geant4TrackPersistency(); + /// Access the Geant4 tracking manager. Only use between tracking pre- and post action + G4TrackingManager* trackMgr() const { return m_context->trackMgr(); } + /// User stepping callback + virtual void step(const G4Step* step, G4SteppingManager* mgr); + /// Pre-event action callback + virtual void beginEvent(const G4Event* event); + /// Post-event action callback + virtual void endEvent(const G4Event* event); + /// Pre-track action callback + virtual void begin(const G4Track* track); + /// Post-track action callback + virtual void end(const G4Track* track); + + /// Mark a Geant4 track to be kept for later MC truth analysis + virtual void mark(const G4Track* track); + /// Store a track + virtual void mark(const G4Track* track, bool created_hit); + /// Mark a Geant4 track of the step to be kept for later MC truth analysis + virtual void mark(const G4Step* step); + /// Store a track produced in a step to be kept for later MC truth analysis + virtual void mark(const G4Step* step, bool created_hit); + + }; + } // End namespace Simulation +} // End namespace DD4hep + +#endif // DD4HEP_DDG4_GEANT4TRACKPERSISTENCY_H diff --git a/DDG4/include/DDG4/Geant4TrackingAction.h b/DDG4/include/DDG4/Geant4TrackingAction.h index a896a1e63c037e567b1aa793fca50de8c41dcaf8..e8cc401ffa56c682670c8703db298da8f29cc923 100644 --- a/DDG4/include/DDG4/Geant4TrackingAction.h +++ b/DDG4/include/DDG4/Geant4TrackingAction.h @@ -42,10 +42,12 @@ namespace DD4hep { Geant4TrackingAction(Geant4Context* context, const std::string& name = ""); /// Default destructor virtual ~Geant4TrackingAction(); - /// Access the tracking manager + /// Access the Geant4 tracking manager. Only use between tracking pre- and post action G4TrackingManager* trackMgr() const { return m_context->trackMgr(); } + /// Mark the track to be kept for MC truth propagation + void mark(const G4Track* track) const; /// Get the valid Geant4 tarck information Geant4TrackInformation* trackInfo(G4Track* track) const; /// Mark all children of the track to be stored @@ -71,10 +73,14 @@ namespace DD4hep { */ class Geant4TrackingActionSequence: public Geant4Action { protected: + /// Callback sequence for pre tracking action + CallbackSequence m_front; /// Callback sequence for pre tracking action CallbackSequence m_begin; /// Callback sequence for post tracking action CallbackSequence m_end; + /// Callback sequence for pre tracking action + CallbackSequence m_final; /// The list of action objects to be called Actors<Geant4TrackingAction> m_actors; public: @@ -82,13 +88,29 @@ namespace DD4hep { Geant4TrackingActionSequence(Geant4Context* context, const std::string& name); /// Default destructor virtual ~Geant4TrackingActionSequence(); + /// Register Pre-track action callback before anything else + template <typename Q, typename T> + void callUpFront(Q* p, void (T::*f)(const G4Track*), + CallbackSequence::Location where=CallbackSequence::END) { + m_front.add(p, f, where); + } /// Register Pre-track action callback - template <typename Q, typename T> void callAtBegin(Q* p, void (T::*f)(const G4Track*)) { - m_begin.add(p, f); + template <typename Q, typename T> + void callAtBegin(Q* p, void (T::*f)(const G4Track*), + CallbackSequence::Location where=CallbackSequence::END) { + m_begin.add(p, f, where); + } + /// Register Post-track action callback + template <typename Q, typename T> + void callAtEnd(Q* p, void (T::*f)(const G4Track*), + CallbackSequence::Location where=CallbackSequence::END) { + m_end.add(p, f, where); } /// Register Post-track action callback - template <typename Q, typename T> void callAtEnd(Q* p, void (T::*f)(const G4Track*)) { - m_end.add(p, f); + template <typename Q, typename T> + void callAtFinal(Q* p, void (T::*f)(const G4Track*), + CallbackSequence::Location where=CallbackSequence::END) { + m_final.add(p, f, where); } /// Add an actor responding to all callbacks. Sequence takes ownership. void adopt(Geant4TrackingAction* action); diff --git a/DDG4/include/DDG4/Geant4UIManager.h b/DDG4/include/DDG4/Geant4UIManager.h new file mode 100644 index 0000000000000000000000000000000000000000..7c4985e4f51596ff31566550b3d2bf691a42639c --- /dev/null +++ b/DDG4/include/DDG4/Geant4UIManager.h @@ -0,0 +1,73 @@ +// $Id: Geant4Field.cpp 875 2013-11-04 16:15:14Z markus.frank@cern.ch $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4UIMANAGER_H +#define DD4HEP_DDG4_GEANT4UIMANAGER_H + +// Framework include files +#include "DDG4/Geant4Call.h" +#include "DDG4/Geant4Action.h" + +/// Forward declarations +class G4VisManager; +class G4UImanager; +class G4UIExecutive; + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + + /* + * Simulation namespace declaration + */ + namespace Simulation { + + /** @class Geant4UIManager Geant4UIManager.h DDG4/Geant4UIManager.h + * + * Standard UI interface implementation + * + * @author M.Frank + * @version 1.0 + */ + class Geant4UIManager : public Geant4Action, virtual public Geant4Call { + protected: + /// Name of the default session type (="cmd") + std::string m_sessionType; + /// Property: Name of the UI macro file + std::string m_uiSetup; + /// Property: Name of the visualization macro file + std::string m_visSetup; + /// Property: Flag to instantiate Vis manager (default=false, unless m_visSetup set) + bool m_haveVis; + /// Property: Flag to instantiate UI (default=true) + bool m_haveUI; + /// Reference to Geant4 visualtion manager + G4VisManager* m_vis; + /// Reference to Geant4 UI manager + G4UIExecutive* m_ui; + public: + /// Initializing constructor + Geant4UIManager(Geant4Context* context, const std::string& name); + /// Default destructor + virtual ~Geant4UIManager(); + /// Start visualization + G4VisManager* startVis(); + /// Start UI + G4UIExecutive* startUI(); + /// Start manager & session + void start(); + /// Stop and release resources + void stop(); + /// Run UI + virtual void operator()(void* param); + }; + + } // End namespace Simulation +} // End namespace DD4hep +#endif // DD4HEP_DDG4_GEANT4UIMANAGER_H diff --git a/DDG4/include/DDG4/Geant4VolumeManager.h b/DDG4/include/DDG4/Geant4VolumeManager.h index ed6a89d0335537dd37d0ba33a4e7f13f15af35b2..4107fcf35a91e4e20630b0d258a1f8fe23140dda 100644 --- a/DDG4/include/DDG4/Geant4VolumeManager.h +++ b/DDG4/include/DDG4/Geant4VolumeManager.h @@ -78,7 +78,15 @@ namespace DD4hep { template <typename Q> Geant4VolumeManager(const Geometry::Handle<Q>& e) : Base(e), m_isValid(false) { } - + /// Assignment operator + Geant4VolumeManager& operator=(const Geant4VolumeManager& c) { + if ( this != &c ) { + m_element = c.m_element; + m_isValid = c.m_isValid; + } + return *this; + } + /// Helper: Generate placement path from touchable object PlacementPath placementPath(const G4VTouchable* touchable, bool exception = true) const; diff --git a/DDG4/include/DDG4/Parsers.h b/DDG4/include/DDG4/Parsers.h index ec5c5b459ff723e9b105c04677c82393843290f3..eacdb2061e0df50c9752e4802aae76c0a218c654 100644 --- a/DDG4/include/DDG4/Parsers.h +++ b/DDG4/include/DDG4/Parsers.h @@ -11,6 +11,9 @@ #include <list> #include <set> #include <map> +#include "Math/Point3D.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" // ============================================================================ #define PARSERS_DECL_FOR_SINGLE(Type) \ @@ -503,6 +506,178 @@ namespace DD4hep { // return 1; // RETURN } + + // ======================================================================== + /** parse 3D-point + * + * Valid representations of 3D-point: + * + * - a'la python tuple with 3 elements ("canonical") + * - a'la python list with 3 elements + * - tuple or list with named ordered fields + * + * @code + * + * " (1,2,3) " + * " [1,2,3] " + * " [ x : 1, 2, Z:3 ] " + * " [ pX : 1 , PY : 2, 3] " + * + * @endcode + * + * Valid keys for names fields: + * + * @code + * + * "x", "X" , "pX" , "Px" , "PX " + * "y", "Y" , "pY" , "Py" , "PY " + * "z", "Z" , "pZ" , "Pz" , "PZ " + * + * @endcode + * + * @attention Named fields must be ordered <code>(x,y,z)</code> + * + * @param result (output) the parsed point + * @param input (input) the input string + * @return status code + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse(ROOT::Math::XYZPoint& result, const std::string& input); + + // ======================================================================== + /** parse 3D-vector + * + * Valid representations of 3D-vector: + * + * - a'la python tuple with 3 elements ("canonical") + * - a'la python list with 3 elements + * - tuple or list with named ordered fields + * + * @code + * + * " (1,2,3) " + * " [1,2,3] " + * " [ x : 1, 2, Z:3 ] " + * " [ pX : 1 , PY : 2, 3] " + * + * @endcode + * + * Valid keys for names fields: + * + * @code + * + * "x", "X" , "pX" , "Px" , "PX " + * "y", "Y" , "pY" , "Py" , "PY " + * "z", "Z" , "pZ" , "Pz" , "PZ " + * + * @endcode + * + * @attention Named fields must be ordered <code>(x,y,z)</code> + * + * @param result (output) the parsed vector + * @param input (input) the input string + * @return status code + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse( ROOT::Math::XYZVector& result, const std::string& input); + + // ======================================================================== + /** parse PxPyPzEVector + * + * Valid representations of Lorenzt vector + * + * - a'la python tuple with 4 elements ("canonical") + * - a'la python list with 4 elements + * - python/list with inner representation of 3D-point/vector + * - tuple or list with named ordered fields + * + * @code + * + * " (1,2,3,4) " + * " (1,2,3;4) " + * + * " [1,2,3,4] " + * " [1,2,3;4] " + * + * " [ x:1 ,2,3; e= 4] " + * " [ pX : 1 , PY : 2, 3 , T= 4] " + * + * " [ ( pX : 1 , PY : 2, 3 ) , 4] " + * " [ ( pX : 1 , PY : 2, 3 ) ; 4] " + * + * " [ 4 , ( pX : 1 , PY : 2, 3 ) ] " + * " [ 4 ; ( pX : 1 , PY : 2, 3 ) ] " + * + * " [ [ pX : 1 , PY : 2, 3 ] , 4] " + * " [ [ pX : 1 , PY : 2, 3 ] ; 4] " + * + * " [ 4 , [ pX : 1 , PY : 2, 3 ] ] " + * " [ 4 ; [ pX : 1 , PY : 2, 3 ] ] " + * + * " ( ( pX : 1 , PY : 2, 3 ) , 4 )" + * " ( ( pX : 1 , PY : 2, 3 ) ; 4 )" + * + * " ( 4 , ( pX : 1 , PY : 2, 3 ) )" + * " ( 4 ; ( pX : 1 , PY : 2, 3 ) )" + * + * " ( [ pX : 1 , PY : 2, 3 ] , 4 )" + * " ( [ pX : 1 , PY : 2, 3 ] ; 4 )" + * + * " ( 4 , [ pX : 1 , PY : 2, 3 ] )" + * " ( 4 ; [ pX : 1 , PY : 2, 3 ] )" + * + * + * @endcode + * + * Note that "eenrgy" element can be separated with semicolumn. + * + * Valid keys for names fields: + * + * @code + * + * "x", "X" , "pX" , "Px" , "PX " + * "y", "Y" , "pY" , "Py" , "PY " + * "z", "Z" , "pZ" , "Pz" , "PZ " + * "t", "T" , "e" , "E" + * + * @endcode + * + * @attention Named fields must be ordered <code>(x,y,z)</code> + * + * @param result (output) the parsed lorentz vector + * @param input (input) the input string + * @return status code + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse(ROOT::Math::PxPyPzEVector& result, const std::string& input); + // ======================================================================== + /** parse the vector of points + * @param resut (OUTPUT) the parser vector + * @param input (INPIUT) the string to be parsed + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse( std::vector<ROOT::Math::XYZPoint>& result, const std::string& input); + // ======================================================================== + /** parse the vector of vectors + * @param resut (OUTPUT) the parser vector + * @param input (INPIUT) the string to be parsed + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse( std::vector<ROOT::Math::XYZVector>& result, const std::string& input); + // ======================================================================== + /** parse the vector of vectors + * @param resut (OUTPUT) the parser vector + * @param input (INPIUT) the string to be parsed + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse(std::vector<ROOT::Math::PxPyPzEVector>& result, const std::string& input); + // ======================================================================== }// end of namespace Parsers // ========================================================================== diff --git a/DDG4/include/DDG4/ToStream.h b/DDG4/include/DDG4/ToStream.h index 2cbe3cef8a28a499af8c1e15dd77b7906536c690..4e59d803ed6fb3a483712a0537c9615557f1cbff 100644 --- a/DDG4/include/DDG4/ToStream.h +++ b/DDG4/include/DDG4/ToStream.h @@ -17,6 +17,9 @@ #include <list> #include <string> #include <sstream> +#include "Math/Point3D.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" // ============================================================================ /** @file DD4hepKernel/ToStream.h @@ -53,12 +56,12 @@ namespace DD4hep { * @date 2009-09-15 */ template <class ITERATOR> - inline std::ostream& toStream(ITERATOR first, // begin of the sequence - ITERATOR last, // end of the sequence - std::ostream& s, // the stream - const std::string& open, // opening - const std::string& close, // closing - const std::string& delim); // delimiter + inline std::ostream& toStream(ITERATOR first, // begin of the sequence + ITERATOR last, // end of the sequence + std::ostream& s, // the stream + const std::string& open, // opening + const std::string& close, // closing + const std::string& delim); // delimiter // ======================================================================== /** the printtout of the strings. * the string is printed a'la Python using the quotes @@ -244,12 +247,12 @@ namespace DD4hep { * @date 2009-09-15 */ template <class ITERATOR> - inline std::ostream& toStream(ITERATOR first, // begin of the sequence - ITERATOR last, // end of the sequence - std::ostream& s, // the stream - const std::string& open, // opening - const std::string& close, // closing - const std::string& delim) // delimiter + inline std::ostream& toStream(ITERATOR first, // begin of the sequence + ITERATOR last, // end of the sequence + std::ostream& s, // the stream + const std::string& open, // opening + const std::string& close, // closing + const std::string& delim) // delimiter { s << open; for (ITERATOR curr = first; curr != last; ++curr) { @@ -279,6 +282,14 @@ namespace DD4hep { s.flags(orig_flags); return s.str(); } + // ============================================================================ + /// print XYZ point + std::ostream& toStream(const ROOT::Math::XYZPoint& obj, std::ostream& s); + // print XYZ-vector + std::ostream& toStream(const ROOT::Math::XYZVector& obj, std::ostream& s); + /// print Lorentz vector + std::ostream& toStream(const ROOT::Math::PxPyPzEVector& obj, std::ostream& s); + // ======================================================================== }// end of namespace DD4hep::Utils // ========================================================================== diff --git a/DDG4/lcio/EventReader.cpp b/DDG4/lcio/EventReader.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bca719711682a5f4e34602648c83e791c6bd0ef1 --- /dev/null +++ b/DDG4/lcio/EventReader.cpp @@ -0,0 +1,303 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// @author P.Kostka (main author) +// @author M.Frank (code reshuffeling into new DDG4 scheme) +// +//==================================================================== + +// Framework include files +#include "EventReader.h" +#include "DD4hep/Plugins.h" +#include "G4Event.hh" +#include "G4PrimaryParticle.hh" +#include "G4ParticleTable.hh" +#include "IMPL/MCParticleImpl.h" +#include "IMPL/LCCollectionVec.h" +#include "Randomize.hh" + +using namespace std; + +/// Initializing constructor +DD4hep::lcio::EventReader::EventReader(const string& nam) : m_name(nam) { +} + +/// Default destructor +DD4hep::lcio::EventReader::~EventReader() { +} + +static void setCharge(IMPL::MCParticleImpl* p) { + // Initialise the particle charge looking at the actual + // Particle Table, or set it to -1000 to flag missing + // information, if PDG unknown in the actual physics list.. + G4ParticleTable* tab = G4ParticleTable::GetParticleTable(); + G4ParticleDefinition* def = tab->FindParticle(p->getPDG()); + double charge = def ? def->GetPDGCharge() : -1000; + p->setCharge(charge); +} + +/// Standard constructor +DD4hep::lcio::LcioGeneratorAction::LcioGeneratorAction(Simulation::Geant4Context* context, const string& name) + : Simulation::Geant4GeneratorAction(context,name) +{ + declareProperty("Input", m_input); + declareProperty("Sync", m_SYNCEVT=0); + declareProperty("zSpread", m_primaryVertexSpreadZ=0.0); + declareProperty("lorentzAngle",m_lorentzTransformationAngle=0.0); +} + +/// Default destructor +DD4hep::lcio::LcioGeneratorAction::~LcioGeneratorAction() { +} + +static string issue(int i) { + stringstream str; + str << "lcio::LcioGeneratorAction: Event " << i; + return str.str(); +} + +/// Read an event and return a LCCollectionVec of MCParticles. +DD4hep::lcio::LcioGeneratorAction::Particles *DD4hep::lcio::LcioGeneratorAction::readEvent(int which) { + int evid = which + m_SYNCEVT; + if ( 0 == m_reader ) { + if ( m_input.empty() ) { + throw runtime_error("LcioGeneratorAction: No inoput file declared!"); + } + string err; + Simulation::TypeName tn = Simulation::TypeName::split(m_input,"|"); + try { + m_reader = PluginService::Create<DD4hep::lcio::EventReader*>(tn.first,tn.second,int(0)); + if ( 0 == m_reader ) { + abortRun(issue(evid),"Error creating reader plugin.", + "Failed to create file reader of type %s. Cannot open dataset %s", + tn.first.c_str(),tn.second.c_str()); + return 0; + } + } + catch(const std::exception& e) { + err = e.what(); + } + if ( !err.empty() ) { + abortRun(issue(evid),err,"Error when reading file %s",m_input.c_str()); + return 0; + } + } + Particles* parts = m_reader->readEvent(); + if ( 0 == parts ) { + abortRun(issue(evid),"Error when reading file - may be end of file.", + "Error when reading file %s",m_input.c_str()); + } + return parts; +} + +/// Callback to generate primary particles +void DD4hep::lcio::LcioGeneratorAction::operator()(G4Event* evt) { + //********************************************************** + // Read in the event + //********************************************************** + m_g4ParticleMap.clear(); + m_convertedParticles.clear(); + Particles* primaries = readEvent(evt->GetEventID()); + if ( 0 == primaries ) return; + + vector<MCParticle*> col(primaries->size()) ; + int NHEP = col.size(); + + // parameters of the Lorentz transformation matrix + const double zspreadparameter = m_primaryVertexSpreadZ; + const double alpha = m_lorentzTransformationAngle; + const double gamma = sqrt(1 + sqr(tan(alpha))); + const double betagamma = tan(alpha); + + double zspread = ( zspreadparameter == 0.0 ? 0.0 : G4RandGauss::shoot(0,zspreadparameter/mm ) ); + //cout << " PrimaryVertexSpreadZ=" << zspreadparameter/mm <<"mm, zspread=" << zspread << "mm" << endl; + G4ThreeVector particle_position = G4ThreeVector(0,0,zspread); + + // build collection of MCParticles + for(int IHEP=0; IHEP<NHEP; IHEP++ ) { + Particle* mcp = dynamic_cast<Particle*> ( primaries->getElementAt( IHEP ) ); + setCharge(mcp); + //Boost to lab frame with crossing angle alpha + if (alpha != 0) { + + double c_light = 299.792; + const double* p = mcp->getMomentum(); + const double m = mcp->getMass(); + //double e = mcp->getEnergy(); + // after the transformation (boost in x-direction) + double pPrime[3]; + pPrime[0] = betagamma * sqrt(sqr(p[0])+sqr(p[1])+sqr(p[2])+sqr(m)) + gamma * p[0]; + pPrime[1] = p[1]; + pPrime[2] = p[2]; + + const double* v = mcp->getVertex(); + const double t = mcp->getTime(); + double vPrime[3]; + double tPrime; + tPrime = gamma * t + betagamma * v[0] / c_light; + vPrime[0] = gamma * v[0] + betagamma * c_light * t; + vPrime[1] = v[1]; + vPrime[2] = v[2]; + + // py and pz remain the same, E changes implicitly with px + mcp->setMomentum( pPrime ); + mcp->setVertex( vPrime ); + mcp->setTime( tPrime ); + //cout << IHEP << " alpha=" << alpha <<" gamma=" << gamma << " betagamma=" << betagamma << endl; + // cout << IHEP << " unboosted momentum: (" << e << "," << p[0] << "," << p[1] << "," << p[2] << ")" << endl; + // cout << IHEP << " boosted momentum: (" << mcp->getEnergy() << "," << betagamma * sqrt(sqr(p[0])+sqr(p[1])+sqr(p[2])+sqr(m)) + gamma * p[0] << "," << pPrime[1] << "," << pPrime[2] << ")" << endl; + // cout << IHEP << " (" << mcp->getEnergy()-e << "," << ( betagamma * sqrt(sqr(p[0])+sqr(p[1])+sqr(p[2])+sqr(m)) + gamma * p[0])-p[0] << "," << pPrime[1]-p[1] << "," << pPrime[2]-p[2] << ")" << endl; + //cout << IHEP << " unboosted vertex: (" << t << "," << v[0] << "," << v[1] << "," << v[2] << ")" << endl; + // cout << IHEP << " boosted vertex: (" << tPrime << "," << vPrime[0] << "," << vPrime[1] << "," << vPrime[2] << ")" << endl; + // cout << IHEP << " boosted vertex: (" << tPrime-t << "," << vPrime[0]-v[0] << "," << vPrime[1]-v[1] << "," << vPrime[2]-v[2] << ")" << endl; + } + // apply z spread + if (zspreadparameter != 0){ + const double* v = mcp->getVertex(); + double vPrime[3]; + vPrime[0] = v[0]; + vPrime[1] = v[1]; + vPrime[2] = v[2] + zspread; + mcp->setVertex(vPrime); + } + col[IHEP] = dynamic_cast<MCParticle*>( mcp ); + } + + // check if there is at least one particle + if( NHEP == 0 ) return; + + // create G4PrimaryVertex object + double particle_time = 0.0; + G4PrimaryVertex* vertex = new G4PrimaryVertex(particle_position,particle_time); + // put initial particles to the vertex + for( size_t i=0; i< col.size(); i++ ){ + MCParticle* mcp = col[i]; + if ( mcp->getParents().size()==0 ) { + G4PrimarySet g4set = getRelevantParticles(mcp); + for (G4PrimarySet::iterator setit=g4set.begin(); setit != g4set.end(); setit++ ){ + vertex->SetPrimary(*setit); + //cout << "G4PrimaryParticle ("<< (*setit)->GetPDGcode() << ") added to G4PrimaryVertex." << endl; + } + } + } + // Put the vertex to G4Event object + evt->AddPrimaryVertex(vertex); +} + +DD4hep::lcio::LcioGeneratorAction::G4PrimarySet +DD4hep::lcio::LcioGeneratorAction::getRelevantParticles(MCParticle* p) { + // log each particle which has been called, to avoid double counting and increase efficiency + m_convertedParticles.insert(p); + G4ParticleMap::iterator mcpIT; + G4PrimarySet relevantParticlesSet; //holds all relevant decay particles of p + + // logic starts here: + // CASE1: if p is a stable particle: search for it in G4ParticleMap + // if it is already there: get G4PrimaryParticle version of p from G4ParticleMap + // else: create G4PrimaryParticle version of p and write it to G4ParticleMap + // in the end: insert this single G4PrimaryParticle verion of p to the + // relevant particle list and return this "list". + if (p->getGeneratorStatus() == 1) { + G4PrimaryParticle* g4p; + mcpIT = m_g4ParticleMap.find( p ); + if( mcpIT != m_g4ParticleMap.end() ){ + g4p = mcpIT->second; + } + else { + int IDHEP = p->getPDG(); + double PHEP1 = p->getMomentum()[0]; + double PHEP2 = p->getMomentum()[1]; + double PHEP3 = p->getMomentum()[2]; + double PHEP5 = p->getMass(); + // create G4PrimaryParticle object + g4p = new G4PrimaryParticle(IDHEP, PHEP1*GeV, PHEP2*GeV, PHEP3*GeV); + g4p->SetMass(PHEP5*GeV); + m_g4ParticleMap[ p ] = g4p; + } + relevantParticlesSet.insert(g4p); + return relevantParticlesSet; + } + + // CASE2: if p is not stable: get first decay daughter and calculate the proper time of p + // if proper time is not zero: particle is "relevant", since it carries vertex information + // if p is already in G4ParticleMap: take it + // otherwise: create G4PrimaryParticle version of p and write it to G4ParticleMap + // now collect all relevant particles of all daughters and setup "relevant mother- + // daughter-relations" between relevant decay particles and G4PrimaryParticle version of p + // in the end: insert only the G4PrimaryParticle version of p to the + // relevant particle list and return this "list". The added particle has now the correct pre-assigned + // decay products and (hopefully) the right lifetime. + else if( p->getDaughters().size() > 0 ) { //fg: require that there is at least one daughter - otherwise forget the particle + cout << "case 2" << endl; + // calculate proper time + MCParticle* dp = p->getDaughters()[0]; + double proper_time = fabs((dp->getTime()-p->getTime()) * p->getMass()) / p->getEnergy(); + // fix by S.Morozov for real != 0 + double aPrecision = dp->getTime() * p->getMass() / p->getEnergy(); + double bPrecision = p->getTime() * p->getMass() / p->getEnergy(); + double ten = 10; + double proper_time_Precision = pow(ten,-DBL_DIG)*fmax(fabs(aPrecision),fabs(bPrecision)); + + bool isProperTimeZero = ( proper_time <= proper_time_Precision ); + + // -- remove original --- if (proper_time != 0) { + if ( isProperTimeZero == false ) { + G4PrimaryParticle* g4p; + mcpIT = m_g4ParticleMap.find( p ); + if( mcpIT != m_g4ParticleMap.end() ){ + g4p = mcpIT->second; + } + else { + int IDHEP = p->getPDG(); + double PHEP1 = p->getMomentum()[0]; + double PHEP2 = p->getMomentum()[1]; + double PHEP3 = p->getMomentum()[2]; + double PHEP5 = p->getMass(); + // create G4PrimaryParticle object + g4p = new G4PrimaryParticle( IDHEP, PHEP1*GeV, PHEP2*GeV, PHEP3*GeV ); + g4p->SetMass( PHEP5*GeV ); + g4p->SetProperTime( proper_time*ns ); + m_g4ParticleMap[ p ] = g4p; + G4PrimarySet vec3; + for (size_t i=0; i<p->getDaughters().size(); i++){ + if (m_convertedParticles.count(p->getDaughters()[i])==0) { + G4PrimarySet vec2 = getRelevantParticles(p->getDaughters()[i]); + G4PrimarySet::iterator setit; + for (setit=vec2.begin(); setit != vec2.end(); setit++ ){ + vec3.insert(*setit); + } + } + } + G4PrimarySet::iterator setit; + for (setit=vec3.begin(); setit != vec3.end(); setit++ ){ + g4p->SetDaughter(*setit); + } + } + relevantParticlesSet.insert(g4p); + return relevantParticlesSet; + } + + // CASE3: if p is not stable AND has zero lifetime: forget about p and retrieve all relevant + // decay particles of all daughters of p. In this case this step of recursion is just there for + // collecting all relevant decay products of daughters (and return them). + else { + for (size_t i=0; i<p->getDaughters().size(); i++){ + if (m_convertedParticles.count(p->getDaughters()[i])==0){ + G4PrimarySet vec2 = getRelevantParticles(p->getDaughters()[i]); + G4PrimarySet::iterator setit; + for (setit=vec2.begin(); setit != vec2.end(); setit++ ){ + relevantParticlesSet.insert(*setit); + } + } + } + return relevantParticlesSet; + } + } + //fg: add a default return + return relevantParticlesSet; +} + +#include "DDG4/Factories.h" +DECLARE_GEANT4ACTION_NS(DD4hep::lcio,LcioGeneratorAction) diff --git a/DDG4/lcio/EventReader.h b/DDG4/lcio/EventReader.h new file mode 100644 index 0000000000000000000000000000000000000000..00219a9332cc0a13715a4e36577ca6ad2fa5801e --- /dev/null +++ b/DDG4/lcio/EventReader.h @@ -0,0 +1,125 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +//==================================================================== +#ifndef DD4HEP_DDG4_LCIO_EVENTREADER_H +#define DD4HEP_DDG4_LCIO_EVENTREADER_H + +// Framework include files +#include "DD4hep/Primitives.h" +#include "DDG4/Geant4GeneratorAction.h" + +// C/C++ include files +#include <set> +#include <map> + +// Forward declarations +class G4Event; +class G4PrimaryParticle; +namespace IO { class LCReader; } +namespace UTIL { class LCStdHepRdr; } +namespace EVENT { class MCParticle; } +namespace IMPL { class MCParticleImpl; } +namespace IMPL { class LCCollectionVec; } + + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + /* + * lcio namespace declaration + */ + namespace lcio { + + /** @class EventReader EventReader.h DDG4/EventReader.h + * + * Base class to read lcio files + * + * @author P.Kostka (main author) + * @author M.Frank (code reshuffeling into new DDG4 scheme) + * @version 1.0 + */ + class EventReader { + public: + typedef IMPL::LCCollectionVec Particles; + protected: + /// File name + std::string m_name; + public: + /// Initializing constructor + EventReader(const std::string& nam); + /// Default destructor + virtual ~EventReader(); + /// File name + const std::string& name() const { return m_name; } + /// Read an event and return a LCCollectionVec of MCParticles. + virtual Particles *readEvent() = 0; + }; + + /** @class LcioGeneratorAction Geant4GeneratorAction.h DDG4/Geant4GeneratorAction.h + * + * Concrete implementation of the Geant4 generator action base class + * populating Geant4 primaries from LCIO and HepStd files. + * + * @author P.Kostka (main author) + * @author M.Frank (code reshuffeling into new DDG4 scheme) + * @version 1.0 + */ + class LcioGeneratorAction : public Simulation::Geant4GeneratorAction { + public: + typedef EVENT::MCParticle MCParticle; + typedef IMPL::MCParticleImpl Particle; + typedef IMPL::LCCollectionVec Particles; + typedef std::set<MCParticle*> LCPrimarySet; + typedef std::set<G4PrimaryParticle*> G4PrimarySet; + typedef std::map<MCParticle*,G4PrimaryParticle*> G4ParticleMap; + protected: + /// Property: input file + std::string m_input; + /// Property: spread of primary vertex in Z + double m_primaryVertexSpreadZ; + /// Property: lorentz transformation angle + double m_lorentzTransformationAngle; + /// Property: SYNCEVT + int m_SYNCEVT; + /// map of partciles already import to Geant4 event + G4ParticleMap m_g4ParticleMap; + /// Set of LCIO particles already converted + LCPrimarySet m_convertedParticles; + /// Event reader object + EventReader* m_reader; + + /// Read an event and return a LCCollectionVec of MCParticles. + Particles* readEvent(int which); + public: + /// Standard constructor + LcioGeneratorAction(Simulation::Geant4Context* context, const std::string& name); + /// Default destructor + virtual ~LcioGeneratorAction(); + /// Callback to generate primary particles + virtual void operator()(G4Event*); + G4PrimarySet getRelevantParticles(MCParticle* p); + }; + } /* End namespace lcio */ +} /* End namespace DD4hep */ + +#include "DD4hep/Plugins.h" +namespace { + /// Factory to create Geant4 physics constructions + template <typename P> class Factory<P, DD4hep::lcio::EventReader*(std::string, int)> { public: + static void Func(void *ret, void*, const std::vector<void*>& a, void*) + { *(DD4hep::lcio::EventReader**)ret = (DD4hep::lcio::EventReader*)new P(*(std::string*)a[0],*(int*)a[1]);} + }; +} +/// Plugin defintion to create event reader objects +#define DECLARE_LCIO_EVENT_READER(name) \ + PLUGINSVC_FACTORY_WITH_ID(name,std::string(#name),DD4hep::lcio::EventReader*(std::string,int)) + +/// Plugin defintion to create event reader objects +#define DECLARE_LCIO_EVENT_READER_NS(ns,name) typedef ns::name __##name##__; \ + PLUGINSVC_FACTORY_WITH_ID(__##name##__,std::string(#name),DD4hep::lcio::EventReader*(std::string,int)) + +#endif /* DD4HEP_DDG4_LCIO_EVENTREADER_H */ diff --git a/DDG4/lcio/HepEventReader.cpp b/DDG4/lcio/HepEventReader.cpp new file mode 100644 index 0000000000000000000000000000000000000000..3548c10dfe55cb75eb1aab0011962926d743e340 --- /dev/null +++ b/DDG4/lcio/HepEventReader.cpp @@ -0,0 +1,233 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +//==================================================================== + +// Framework include files +#include "EventReader.h" + +// C/C++ include files +#include <fstream> + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + /* + * lcio namespace declaration + */ + namespace lcio { + + /** @class HepEventReader HepEventReader.h DDG4/HepEventReader.h + * + * Class to populate Geant4 primary particles and vertices from a + * file in StdHep format (ASCII) + * + * @author P.Kostka (main author) + * @author M.Frank (code reshuffeling into new DDG4 scheme) + * @version 1.0 + */ + struct HepEventReader : public EventReader { + 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 Particles *readEvent(); + }; + } /* 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; +#define HEPEvt 1 + +// Factory entry +typedef DD4hep::lcio::HepEventReader HepEventReader; +DECLARE_LCIO_EVENT_READER(HepEventReader) + +/// Initializing constructor +HepEventReader::HepEventReader(const std::string& nam, int fmt) +: EventReader(nam), m_format(fmt) +{ +} + +/// 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. +HepEventReader::Particles *HepEventReader::readEvent() { + 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/LcioEventReader.cpp b/DDG4/lcio/LcioEventReader.cpp new file mode 100644 index 0000000000000000000000000000000000000000..25e4bbf890a9f76424c0aad90349afa030146dd4 --- /dev/null +++ b/DDG4/lcio/LcioEventReader.cpp @@ -0,0 +1,77 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +//==================================================================== + +// Framework include files +#include "EventReader.h" + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + /* + * lcio namespace declaration + */ + namespace lcio { + + /** @class LcioEventReader LcioEventReader.h DDG4/LcioEventReader.h + * + * Base class to read lcio event files + * + * @author P.Kostka (main author) + * @author M.Frank (code reshuffeling into new DDG4 scheme) + * @version 1.0 + */ + struct LcioEventReader : public EventReader { + protected: + /// Reference to reader object + IO::LCReader* m_reader; + public: + /// Initializing constructor + LcioEventReader(const std::string& nam, int); + /// Default destructor + virtual ~LcioEventReader(); + /// Read an event and return a LCCollectionVec of MCParticles. + virtual Particles *readEvent(); + }; + } +} + +#include "DD4hep/Printout.h" +#include "lcio.h" +#include "EVENT/LCIO.h" +#include "IO/LCReader.h" + +using namespace EVENT; + +// Factory entry +typedef DD4hep::lcio::LcioEventReader LcioEventReader; +DECLARE_LCIO_EVENT_READER(LcioEventReader) + +/// Initializing constructor +LcioEventReader::LcioEventReader(const std::string& nam, int /* arg */) +: EventReader(nam) +{ + m_reader = ::lcio::LCFactory::getInstance()->createLCReader(); + printout(INFO,"LcioEventReader","Created file reader. Try to open input %s",nam.c_str()); + m_reader->open(nam); +} + +/// Default destructor +LcioEventReader::~LcioEventReader() { + deletePtr(m_reader); +} + +/// Read an event and return a LCCollectionVec of MCParticles. +DD4hep::lcio::EventReader::Particles *LcioEventReader::readEvent() { + ::lcio::LCEvent* evt = m_reader->readNextEvent(); + if ( evt ) { + Particles* mcVec = (Particles*)evt->getCollection(LCIO::MCPARTICLE); + return mcVec; + } + return 0; +} + diff --git a/DDG4/lcio/LcioStdHepReader.cpp b/DDG4/lcio/LcioStdHepReader.cpp new file mode 100644 index 0000000000000000000000000000000000000000..0bf0b6b3124b2534e234c74d46ef5363fa7b6b94 --- /dev/null +++ b/DDG4/lcio/LcioStdHepReader.cpp @@ -0,0 +1,68 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +//==================================================================== + +// Framework include files +#include "EventReader.h" + +/* + * DD4hep namespace declaration + */ +namespace DD4hep { + /* + * lcio namespace declaration + */ + namespace lcio { + + /** @class StdHepReader StdHepReader.h DDG4/StdHepReader.h + * + * Base class to read StdHep files with lcio + * + * @author P.Kostka (main author) + * @author M.Frank (code reshuffeling into new DDG4 scheme) + * @version 1.0 + */ + struct StdHepReader : public EventReader { + protected: + /// Reference to reader object + UTIL::LCStdHepRdr* m_reader; + public: + /// Initializing constructor + StdHepReader(const std::string& nam, int); + /// Default destructor + virtual ~StdHepReader(); + /// Read an event and return a LCCollectionVec of MCParticles. + virtual Particles *readEvent(); + }; + } /* End namespace lcio */ +} /* End namespace DD4hep */ + + +#include "lcio.h" +#include "EVENT/LCIO.h" +#include "UTIL/LCStdHepRdr.h" + +// Factory entry +typedef DD4hep::lcio::StdHepReader LcioStdHepReader; +DECLARE_LCIO_EVENT_READER(LcioStdHepReader) + +/// Initializing constructor +DD4hep::lcio::StdHepReader::StdHepReader(const std::string& nam, int) + : EventReader(nam) +{ + m_reader = new ::lcio::LCStdHepRdr(m_name.c_str()); +} + +/// Default destructor +DD4hep::lcio::StdHepReader::~StdHepReader() { + deletePtr(m_reader); +} + +/// Read an event and return a LCCollectionVec of MCParticles. +DD4hep::lcio::EventReader::Particles *DD4hep::lcio::StdHepReader::readEvent() { + return m_reader->readEvent(); +} + diff --git a/DDG4/parsers/Grammars.h b/DDG4/parsers/Grammars.h index 81140ff5464758857924d5017b0d0a68ade5f607..37513082bd99b064d956551a603c01800441436b 100644 --- a/DDG4/parsers/Grammars.h +++ b/DDG4/parsers/Grammars.h @@ -3,7 +3,7 @@ #define DD4HEPKERNEL_GRAMMARS_H 1 #ifdef __GNUC__ #warning \ - The headers GaudiKernel/Grammars.h and GaudiKernel/Parsers.icpp are deprecated \ + The headers Grammars.h and Parsers.icpp are deprecated \ and will be removed from the next release of Gaudi. You should migrate your \ code the new pasers based on Boost.Spirit 2. #endif diff --git a/DDG4/parsers/GrammarsV2.h b/DDG4/parsers/GrammarsV2.h index 419de5d6529461d98187844c1f69e28560658db2..d9b3f3c48fac3b47797f32f3b18a5ad75afe9142 100644 --- a/DDG4/parsers/GrammarsV2.h +++ b/DDG4/parsers/GrammarsV2.h @@ -29,6 +29,10 @@ #include <boost/spirit/repository/include/qi_confix.hpp> +#include "Math/Point3D.h" +#include "Math/Vector3D.h" +#include "Math/Vector4D.h" + //============================================================================== namespace DD4hep { namespace Parsers { //============================================================================== @@ -39,14 +43,16 @@ namespace DD4hep { namespace Parsers { namespace qi = sp::qi; namespace enc = sp::ascii; namespace rep = sp::repository; + + template <typename T> T evaluate_string(const std::string& value); + //============================================================================== // Grammars //============================================================================== typedef std::string::const_iterator DefaultIterator; typedef enc::space_type DefaultSkipper; //============================================================================== - template <typename Iterator, typename T, typename Skipper, - class Enable = void> + template <typename Iterator, typename T, typename Skipper, class Enable=void> struct Grammar_ { /* READ THIS IF YOUR COMPILE BREAKS ON THE FOLLOWING LINE * @@ -64,8 +70,7 @@ namespace DD4hep { namespace Parsers { typedef GrammarName<Iterator, Skipper> Grammar; \ } //============================================================================== - template< typename Iterator> - struct SkipperGrammar : qi::grammar<Iterator> + template< typename Iterator> struct SkipperGrammar : qi::grammar<Iterator> { SkipperGrammar() : SkipperGrammar::base_type(comments) { comments = enc::space | rep::confix("/*", "*/")[*(qi::char_ - "*/")] @@ -76,41 +81,39 @@ namespace DD4hep { namespace Parsers { }; //============================================================================== template< typename Iterator, typename Skipper> - struct StringGrammar : qi::grammar<Iterator, std::string(), qi::locals<char>, - Skipper> - { - //------------------------------------------------------------------------------ - typedef std::string ResultT; - //------------------------------------------------------------------------------ - StringGrammar() : StringGrammar::base_type( str ) { - begin_quote = enc::char_("\"'"); - quote = enc::char_(qi::_r1); + struct StringGrammar : qi::grammar<Iterator, std::string(), qi::locals<char>,Skipper> + { + //------------------------------------------------------------------------------ + typedef std::string ResultT; + //------------------------------------------------------------------------------ + StringGrammar() : StringGrammar::base_type( str ) { + begin_quote = enc::char_("\"'"); + quote = enc::char_(qi::_r1); - str = qi::lexeme[begin_quote[qi::_a = qi::_1] - > *( (enc::char_('\\') >> quote(qi::_a))[qi::_val += qi::_a] - | (enc::char_[qi::_val += qi::_1] - quote(qi::_a))) > - quote(qi::_a)] - ; - } - //------------------------------------------------------------------------------ - qi::rule<Iterator, std::string(), qi::locals<char>, Skipper> str; - qi::rule<Iterator, char()> begin_quote; - qi::rule<Iterator, void(char)> quote; - //------------------------------------------------------------------------------ - }; + str = qi::lexeme[begin_quote[qi::_a = qi::_1] + > *( (enc::char_('\\') >> quote(qi::_a))[qi::_val += qi::_a] + | (enc::char_[qi::_val += qi::_1] - quote(qi::_a))) > + quote(qi::_a)] + ; + } + //------------------------------------------------------------------------------ + qi::rule<Iterator, std::string(), qi::locals<char>, Skipper> str; + qi::rule<Iterator, char()> begin_quote; + qi::rule<Iterator, void(char)> quote; + //------------------------------------------------------------------------------ + }; REGISTER_GRAMMAR(std::string, StringGrammar); //============================================================================== template< typename Iterator, typename Skipper> - struct CharGrammar : qi::grammar<Iterator, char(), Skipper> - { - typedef char ResultT; - CharGrammar() : CharGrammar::base_type( ch ) { - ch = qi::int_parser<char>() - | - '\'' >> (qi::char_-'\'') >> '\''; - } - qi::rule<Iterator, char(), Skipper> ch; - }; + struct CharGrammar : qi::grammar<Iterator, char(), Skipper> { + typedef char ResultT; + CharGrammar() : CharGrammar::base_type( ch ) { + ch = qi::int_parser<char>() + | + '\'' >> (qi::char_-'\'') >> '\''; + } + qi::rule<Iterator, char(), Skipper> ch; + }; REGISTER_GRAMMAR(char, CharGrammar); //============================================================================== template< typename Iterator, typename Skipper> @@ -128,8 +131,7 @@ namespace DD4hep { namespace Parsers { REGISTER_GRAMMAR(bool, BoolGrammar); //============================================================================== template< typename Iterator, typename RT , typename Skipper> - struct IntGrammar : qi::grammar<Iterator, RT(), Skipper> - { + struct IntGrammar : qi::grammar<Iterator, RT(), Skipper> { typedef RT ResultT; IntGrammar() : IntGrammar::base_type( integer ) { integer = qi::int_parser<RT>()[qi::_val = qi::_1] @@ -148,8 +150,7 @@ namespace DD4hep { namespace Parsers { }; //============================================================================== template< typename Iterator, typename RT, typename Skipper> - struct RealGrammar : qi::grammar<Iterator, RT(), Skipper> - { + struct RealGrammar : qi::grammar<Iterator, RT(), Skipper> { typedef RT ResultT; RealGrammar() : RealGrammar::base_type(real) { real = qi::real_parser<RT>(); @@ -161,14 +162,12 @@ namespace DD4hep { namespace Parsers { // ---------------------------------------------------------------------------- template <typename Iterator, typename T, typename Skipper > struct Grammar_<Iterator, T, Skipper, - typename boost::enable_if<boost::is_floating_point<T> >::type > - { + typename boost::enable_if<boost::is_floating_point<T> >::type > { typedef RealGrammar<Iterator, T, Skipper> Grammar; }; //============================================================================== template< typename Iterator, typename VectorT, typename Skipper> - struct VectorGrammar : qi::grammar<Iterator, - VectorT(), qi::locals<char>,Skipper> + struct VectorGrammar : qi::grammar<Iterator, VectorT(), qi::locals<char>,Skipper> { //------------------------------------------------------------------------------ typedef VectorT ResultT; @@ -194,36 +193,30 @@ namespace DD4hep { namespace Parsers { // ---------------------------------------------------------------------------- // Register VectorGrammar for std::vector: // ---------------------------------------------------------------------------- - template <typename Iterator, typename InnerT, typename AllocatorT, - typename Skipper> - struct Grammar_<Iterator, std::vector<InnerT, AllocatorT>, Skipper > - { - typedef - VectorGrammar<Iterator, std::vector<InnerT, AllocatorT>,Skipper> - Grammar; - }; + template <typename Iterator,typename InnerT,typename AllocatorT,typename Skipper> + struct Grammar_<Iterator, std::vector<InnerT, AllocatorT>, Skipper > { + typedef + VectorGrammar<Iterator, std::vector<InnerT, AllocatorT>,Skipper> + Grammar; + }; // ---------------------------------------------------------------------------- // Register VectorGrammar for std::list: // ---------------------------------------------------------------------------- - template <typename Iterator, typename InnerT, typename AllocatorT, - typename Skipper> - struct Grammar_<Iterator, std::list<InnerT, AllocatorT>, Skipper > - { - typedef - VectorGrammar<Iterator, std::list<InnerT, AllocatorT>,Skipper> - Grammar; - }; + template <typename Iterator, typename InnerT, typename AllocatorT,typename Skipper> + struct Grammar_<Iterator, std::list<InnerT, AllocatorT>, Skipper > { + typedef + VectorGrammar<Iterator, std::list<InnerT, AllocatorT>,Skipper> + Grammar; + }; // ---------------------------------------------------------------------------- // Register VectorGrammar for std::set: // ---------------------------------------------------------------------------- - template <typename Iterator, typename InnerT, typename CompareT, - typename AllocatorT, typename Skipper> - struct Grammar_<Iterator, std::set<InnerT, CompareT, AllocatorT>, Skipper > - { - typedef - VectorGrammar<Iterator, std::set<InnerT, CompareT, AllocatorT>,Skipper> - Grammar; - }; + template <typename Iterator, typename InnerT, typename CompareT,typename AllocatorT, typename Skipper> + struct Grammar_<Iterator, std::set<InnerT, CompareT, AllocatorT>, Skipper > { + typedef + VectorGrammar<Iterator, std::set<InnerT, CompareT, AllocatorT>,Skipper> + Grammar; + }; //============================================================================== template< typename Iterator, typename PairT, typename Skipper> @@ -276,58 +269,56 @@ namespace DD4hep { namespace Parsers { }; // ============================================================================ template< typename Iterator, typename MapT, typename Skipper> - struct MapGrammar : qi::grammar<Iterator,MapT(), Skipper> - { - //------------------------------------------------------------------------------ - typedef MapT ResultT; - typedef typename MapT::key_type KeyT; - typedef typename MapT::mapped_type MappedT; - typedef std::pair<KeyT, MappedT> PairT; + struct MapGrammar : qi::grammar<Iterator,MapT(), Skipper> { + //------------------------------------------------------------------------------ + typedef MapT ResultT; + typedef typename MapT::key_type KeyT; + typedef typename MapT::mapped_type MappedT; + typedef std::pair<KeyT, MappedT> PairT; - typedef std::vector<PairT> VectorPairT; - //------------------------------------------------------------------------------ - struct tag_key{}; - struct tag_mapped{}; - struct Operations - { - template <typename A, typename B = boost::fusion::unused_type, - typename C = boost::fusion::unused_type, - typename D = boost::fusion::unused_type> - struct result { typedef void type; }; - //---------------------------------------------------------------------- - void operator()(ResultT& res, const VectorPairT& vec) const{ - for(typename VectorPairT::const_iterator cur = vec.begin(); - cur != vec.end(); cur++){ - res.insert(*cur); - } - } - void operator()(PairT& res, const KeyT& key, tag_key) const{ - res.first = key; - } - void operator()(PairT& res, const MappedT& value, tag_mapped) const{ - res.second = value; + typedef std::vector<PairT> VectorPairT; + //------------------------------------------------------------------------------ + struct tag_key{}; + struct tag_mapped{}; + struct Operations { + template <typename A, typename B = boost::fusion::unused_type, + typename C = boost::fusion::unused_type, + typename D = boost::fusion::unused_type> + struct result { typedef void type; }; + //---------------------------------------------------------------------- + void operator()(ResultT& res, const VectorPairT& vec) const{ + for(typename VectorPairT::const_iterator cur = vec.begin(); + cur != vec.end(); cur++){ + res.insert(*cur); } - //---------------------------------------------------------------------- - }; - //------------------------------------------------------------------------------ - MapGrammar() : MapGrammar::base_type(map) { - pair = key[op(qi::_val,qi::_1, tag_key())] > (qi::lit(':') | '=') > - value[op(qi::_val,qi::_1, tag_mapped())]; - list = -(pair % enc::char_(',')); - map = (('[' >> list >> ']') - | ('{' >> list >> '}'))[op(qi::_val,qi::_1)]; } - // ---------------------------------------------------------------------------- - typename - Grammar_<Iterator, typename MapT::key_type, Skipper>::Grammar key; - typename - Grammar_<Iterator, typename MapT::mapped_type, Skipper>::Grammar value; - qi::rule<Iterator, PairT(), Skipper> pair; - qi::rule<Iterator, VectorPairT(), Skipper> list; - qi::rule<Iterator, ResultT(), Skipper> map; - ph::function<Operations> op; - // ---------------------------------------------------------------------------- + void operator()(PairT& res, const KeyT& key, tag_key) const{ + res.first = key; + } + void operator()(PairT& res, const MappedT& value, tag_mapped) const{ + res.second = value; + } + //---------------------------------------------------------------------- }; + //------------------------------------------------------------------------------ + MapGrammar() : MapGrammar::base_type(map) { + pair = key[op(qi::_val,qi::_1, tag_key())] > (qi::lit(':') | '=') > + value[op(qi::_val,qi::_1, tag_mapped())]; + list = -(pair % enc::char_(',')); + map = (('[' >> list >> ']') + | ('{' >> list >> '}'))[op(qi::_val,qi::_1)]; + } + // ---------------------------------------------------------------------------- + typename + Grammar_<Iterator, typename MapT::key_type, Skipper>::Grammar key; + typename + Grammar_<Iterator, typename MapT::mapped_type, Skipper>::Grammar value; + qi::rule<Iterator, PairT(), Skipper> pair; + qi::rule<Iterator, VectorPairT(), Skipper> list; + qi::rule<Iterator, ResultT(), Skipper> map; + ph::function<Operations> op; + // ---------------------------------------------------------------------------- + }; // ---------------------------------------------------------------------------- // Register MapGrammar for std::map: // ---------------------------------------------------------------------------- @@ -359,6 +350,121 @@ namespace DD4hep { namespace Parsers { // ---------------------------------------------------------------------------- }; // END KeyValueGrammar // We don't register KeyalueGrammar because it's a special parser + + + // ============================================================================ + template< typename Iterator, typename PointT, typename Skipper> + struct Pnt3DGrammar : qi::grammar<Iterator, PointT(), Skipper> { + typedef PointT ResultT; + typedef std::string Scalar; + // ---------------------------------------------------------------------------- + struct Operations { + template <typename A, typename B = boost::fusion::unused_type, + typename C = boost::fusion::unused_type, + typename D = boost::fusion::unused_type> + struct result { typedef void type; }; + void operator()(ResultT& res, const Scalar& value,const char xyz) const{ + typename PointT::Scalar val = evaluate_string<typename PointT::Scalar>(value); + switch(xyz) { + case 'x': res.SetX(val); break; + case 'y': res.SetY(val); break; + case 'z': res.SetZ(val); break; + default: break; + } + } + }; // Operations + // ---------------------------------------------------------------------------- + Pnt3DGrammar() : Pnt3DGrammar::base_type(point) { + point = list | ('(' >> list >> ')') | ('[' >> list >> ']'); + list = -(enc::no_case[qi::lit("x") | qi::lit("px")] >> ':') + >> scalar[op(qi::_val,qi::_1,'x')] >> + ',' >> -(enc::no_case[qi::lit("y") | qi::lit("py")] >> ':') + >> scalar[op(qi::_val,qi::_1,'y')] >> + ',' >> -(enc::no_case[qi::lit("z") | qi::lit("pz")] >> ':') + >> scalar[op(qi::_val,qi::_1,'z')]; + } + // ---------------------------------------------------------------------------- + qi::rule<Iterator, ResultT(), Skipper> point, list; + typename Grammar_<Iterator, Scalar, Skipper>::Grammar scalar; + ph::function<Operations> op; + // ---------------------------------------------------------------------------- + }; // Pnt3DGrammar + // ---------------------------------------------------------------------------- + // Register Pnt3DGrammar for ROOT::Math::PositionVector3D: + // ---------------------------------------------------------------------------- + template <typename Iterator, typename T1, typename T2, typename Skipper> + struct Grammar_<Iterator, ROOT::Math::PositionVector3D<T1,T2>, Skipper>{ + typedef Pnt3DGrammar<Iterator, ROOT::Math::PositionVector3D<T1,T2>, Skipper> Grammar; + }; + // ---------------------------------------------------------------------------- + // Register Pnt3DGrammar for ROOT::Math::DisplacementVector3D: + // ---------------------------------------------------------------------------- + template <typename Iterator, typename T1, typename T2, typename Skipper> + struct Grammar_<Iterator, ROOT::Math::DisplacementVector3D<T1,T2>, Skipper>{ + typedef Pnt3DGrammar<Iterator,ROOT::Math::DisplacementVector3D<T1,T2>, Skipper> Grammar; + }; + // ============================================================================ + template< typename Iterator, typename PointT, typename Skipper> + struct Pnt4DGrammar : qi::grammar<Iterator, PointT(), Skipper> { + typedef PointT ResultT; + typedef std::string ScalarT; + //----------------------------------------------------------------------------- + struct Operations { + template <typename A, typename B = boost::fusion::unused_type, + typename C = boost::fusion::unused_type, + typename D = boost::fusion::unused_type> + struct result { typedef void type; }; + + void operator()(ResultT& res, const ScalarT& value,const char xyz) const{ + typename PointT::Scalar val = evaluate_string<typename PointT::Scalar>(value); + switch(xyz){ + case 'x': res.SetPx(val); break; + case 'y': res.SetPy(val); break; + case 'z': res.SetPz(val); break; + case 'e': res.SetE(val); break; + default: break; + } + } + void operator()(ResultT& res, const ResultT& xyz) const{ + res.SetPx(xyz.Px()); + res.SetPy(xyz.Py()); + res.SetPz(xyz.Pz()); + } + }; // Operations + // ---------------------------------------------------------------------------- + Pnt4DGrammar() : Pnt4DGrammar::base_type(point4d) { + point4d = list4d | ('(' >> list4d >> ')') | ('[' >> list4d >> ']'); + list4d = (point3d[op(qi::_val,qi::_1)] >> enc::char_(";,") + >> e[op(qi::_val, qi::_1, 'e')]) + | + (e[op(qi::_val,qi::_1, 'e')] >> enc::char_(";,") + >> point3d[op(qi::_val, qi::_1)]); + e = -(enc::no_case[enc::char_("te")] >> ':') + >> scalar[qi::_val = qi::_1]; + + point3d = list3d | ('(' >> list3d >> ')') | ('[' >> list3d >> ']'); + list3d = -(enc::no_case[qi::lit("x") | qi::lit("px")] >> ':') + >> scalar[op(qi::_val, qi::_1,'x')] >> + ',' >> -(enc::no_case[qi::lit("y") | qi::lit("py")] >> ':') + >> scalar[op(qi::_val, qi::_1,'y')] >> + ',' >> -(enc::no_case[qi::lit("z") | qi::lit("pz")] >> ':') + >> scalar[op(qi::_val, qi::_1,'z')]; + } + // ---------------------------------------------------------------------------- + qi::rule<Iterator, ResultT(), Skipper> point3d, point4d, list3d, + list4d; + qi::rule<Iterator, ScalarT(), Skipper> e; + typename Grammar_<Iterator, ScalarT, Skipper>::Grammar scalar; + ph::function<Operations> op; + // ---------------------------------------------------------------------------- + }; // Pnt4DGrammar + // ---------------------------------------------------------------------------- + // Register Pnt4DGrammar for ROOT::Math::LorentzVector: + // ---------------------------------------------------------------------------- + template <typename Iterator, typename T1, typename Skipper> + struct Grammar_<Iterator, ROOT::Math::LorentzVector<T1>, Skipper > { + typedef Pnt4DGrammar<Iterator, ROOT::Math::LorentzVector<T1>, Skipper> Grammar; + }; // ============================================================================ }} // DD4hep::Parsers //============================================================================ diff --git a/DDG4/parsers/Parsers.icpp b/DDG4/parsers/Parsers.icpp index a4e1421e60e5741f503864ae27dd6d90f85ad76b..ec4ab1ef090d2092ad899877e0058ccd931be694 100644 --- a/DDG4/parsers/Parsers.icpp +++ b/DDG4/parsers/Parsers.icpp @@ -1,7 +1,7 @@ // $Id: Parsers.icpp,v 1.5 2008/10/28 14:02:18 marcocle Exp $ // ============================================================================ -#ifndef GAUDIKERNEL_PARSERS_ICPP -#define GAUDIKERNEL_PARSERS_ICPP 1 +#ifndef DD4HEP_PARSERS_ICPP +#define DD4HEP_PARSERS_ICPP 1 // ============================================================================ // Include files // ============================================================================ diff --git a/DDG4/parsers/ParsersObjects.cpp b/DDG4/parsers/ParsersObjects.cpp new file mode 100644 index 0000000000000000000000000000000000000000..dc84c595e86c674aa99c20a6d5339fbeef76977f --- /dev/null +++ b/DDG4/parsers/ParsersObjects.cpp @@ -0,0 +1,84 @@ +// ============================================================================ +// Include files +// ============================================================================ +#include "ParsersFactory.h" +#include "DDG4/ToStream.h" + +using namespace std; + +// ============================================================================ +namespace DD4hep { + namespace Parsers { + // ========================================================================== + template<typename T1, typename T2> inline int + parse_(ROOT::Math::PositionVector3D<T1,T2>& result, const string& input){ + Skipper skipper; + typename Grammar_<IteratorT,ROOT::Math::PositionVector3D<T1,T2>,Skipper>::Grammar g; + IteratorT iter = input.begin(), end = input.end(); + if (qi::phrase_parse( iter, end, g, skipper, result)){ + return 1; + } + return 0; + } + // ========================================================================== + int parse(ROOT::Math::XYZPoint& result,const string& input) { + return parse_(result, input); + } + // ========================================================================== + /* parse 3D-vector + * @param result (output) the parsed vector + * @param input (input) the input string + * @return status code + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + int parse(ROOT::Math::XYZVector& result,const string& input) { + ROOT::Math::XYZPoint point; + int sc = parse(point,input); + if ( 0 == sc ){ return sc; } // RETURN + result = point; + return 1; + } + int parse(ROOT::Math::PxPyPzEVector& result, const string& input) { + return parse_(result, input); + } + + // ========================================================================== + /* parse the vector of points + * @param resut (OUTPUT) the parser vector + * @param input (INPIUT) the string to be parsed + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + // ========================================================================== + int parse(vector<ROOT::Math::XYZPoint>& result, const string& input) { + result.clear(); + return parse_(result, input); + } + // ========================================================================== + /* parse the vector of vectors + * @param resut (OUTPUT) the parser vector + * @param input (INPIUT) the string to be parsed + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + // ========================================================================== + int parse(vector<ROOT::Math::XYZVector>& result, const string& input) { + result.clear(); + return parse_(result, input); + } + // ========================================================================== + /* parse the vector of vectors + * @param resut (OUTPUT) the parser vector + * @param input (INPIUT) the string to be parsed + * @author Vanya BELYAEV Ivan.Belyaev@nikhef.nl + * @date 2009-09-05 + */ + // ========================================================================== + int parse(vector<ROOT::Math::PxPyPzEVector>& result,const string& input) { + result.clear(); + return parse_(result, input); + } + + } +} diff --git a/DDG4/plugins/Geant4Factories.cpp b/DDG4/plugins/Geant4Factories.cpp index 46fa29f499ced2b50dac297968ff7381a6217c8d..624eb7acb7bde2ce0d6adfc98cb59d8018c61363 100644 --- a/DDG4/plugins/Geant4Factories.cpp +++ b/DDG4/plugins/Geant4Factories.cpp @@ -34,6 +34,14 @@ DECLARE_GEANT4ACTION(Geant4PhysicsListActionSequence) #include "DDG4/Geant4SensDetAction.h" DECLARE_GEANT4ACTION(Geant4SensDetActionSequence) +#include "DDG4/Geant4UIManager.h" +DECLARE_GEANT4ACTION(Geant4UIManager) + +#include "DDG4/Geant4TrackPersistency.h" +DECLARE_GEANT4ACTION(Geant4TrackPersistency) +#include "DDG4/Geant4MonteCarloRecordManager.h" +DECLARE_GEANT4ACTION(Geant4MonteCarloRecordManager) + //============================= #include "DDG4/Geant4TrackingPreAction.h" DECLARE_GEANT4ACTION(Geant4TrackingPreAction) diff --git a/DDG4/plugins/Geant4SDActions.cpp b/DDG4/plugins/Geant4SDActions.cpp index a7372d17a92b4c1562966b7e468b1e4937bb99f8..d529b5d8eb162b661e157b5d6c615288ad58b5eb 100644 --- a/DDG4/plugins/Geant4SDActions.cpp +++ b/DDG4/plugins/Geant4SDActions.cpp @@ -8,6 +8,7 @@ //==================================================================== // Framework include files #include "DDG4/Geant4SensDetAction.h" +#include "DDG4/Geant4MonteCarloTruth.h" #include "DDG4/Geant4Data.h" #include "DD4hep/Printout.h" @@ -145,6 +146,7 @@ namespace DD4hep { hit->momentum = direction; hit->length = hit_len; collection(m_collectionID)->add(hit); + mcTruthMgr().mark(h.track,true); return hit != 0; } typedef Geant4SensitiveAction<SimpleTracker> Geant4SimpleTrackerAction; @@ -178,6 +180,7 @@ namespace DD4hep { } hit->truth.push_back(contrib); hit->energyDeposit += contrib.deposit; + mcTruthMgr().mark(h.track,true); return true; } typedef Geant4SensitiveAction<SimpleCalorimeter> Geant4SimpleCalorimeterAction; @@ -215,6 +218,7 @@ namespace DD4hep { hit->energyDeposit += contrib.deposit; hit->truth.push_back(contrib); track->SetTrackStatus(fStopAndKill); // don't step photon any further + mcTruthMgr().mark(h.track,true); return true; } } @@ -267,6 +271,7 @@ namespace DD4hep { hit->length = path_len; clear(); c->insert(hit); + mcTruthMgr().mark(h.track,true); return hit; } }; diff --git a/DDG4/plugins/Geant4SensDet.cpp b/DDG4/plugins/Geant4SensDet.cpp index c73b0ee4bacc82ec69435a37a9eb0f9d94ca7651..66c29fb130e4bbb15357c347e2a62f04d592bd50 100644 --- a/DDG4/plugins/Geant4SensDet.cpp +++ b/DDG4/plugins/Geant4SensDet.cpp @@ -25,29 +25,27 @@ namespace DD4hep { */ namespace Simulation { - namespace { - template <typename T> struct _Seq { - typedef _Seq<T> Base; - T* m_sequence; - _Seq() : m_sequence(0) { } - _Seq(T* seq) { _aquire(seq); } - virtual ~_Seq() { _release(); } - void _aquire(T* s) { - InstanceCount::increment(this); - m_sequence = s; - m_sequence->addRef(); - } - void _release() { - releasePtr(m_sequence); - InstanceCount::decrement(this); - } - }; - } + template <typename T> struct RefCountedSequence { + typedef RefCountedSequence<T> Base; + T* m_sequence; + RefCountedSequence() : m_sequence(0) { } + RefCountedSequence(T* seq) { _aquire(seq); } + virtual ~RefCountedSequence() { _release(); } + void _aquire(T* s) { + InstanceCount::increment(this); + m_sequence = s; + m_sequence->addRef(); + } + void _release() { + releasePtr(m_sequence); + InstanceCount::decrement(this); + } + }; struct Geant4SensDet : virtual public G4VSensitiveDetector, virtual public G4VSDFilter, virtual public Geant4ActionSD, - virtual public _Seq<Geant4SensDetActionSequence> + virtual public RefCountedSequence<Geant4SensDetActionSequence> { /// Constructor. The detector element is identified by the name Geant4SensDet(const std::string& nam, Geometry::LCDD& lcdd) diff --git a/DDG4/plugins/Geant4SensDetFilters.cpp b/DDG4/plugins/Geant4SensDetFilters.cpp index d283da00e0f9dced88b9bfdc2487719e87456c67..2fa2a4d96c2034f020c5ff1142ad2fa8b53b21ab 100644 --- a/DDG4/plugins/Geant4SensDetFilters.cpp +++ b/DDG4/plugins/Geant4SensDetFilters.cpp @@ -140,7 +140,7 @@ DECLARE_GEANT4ACTION(EnergyDepositMinimumCut) /// Constructor. ParticleFilter::ParticleFilter(Geant4Context* context, const std::string& name) - : Geant4Filter(context,name) +: Geant4Filter(context,name), m_definition(0) { declareProperty("particle",m_particle); InstanceCount::increment(this); diff --git a/DDG4/plugins/Geant4XMLSetup.cpp b/DDG4/plugins/Geant4XMLSetup.cpp index 49dfd5322b3da8f2187190674d027e609719c709..4b612987b9a7d78df54475f51d091df699dea634 100644 --- a/DDG4/plugins/Geant4XMLSetup.cpp +++ b/DDG4/plugins/Geant4XMLSetup.cpp @@ -86,7 +86,7 @@ namespace DD4hep { Kernel& kernel = Kernel::access(lcdd); Action action((what==FILTER) ? (Geant4Action*)kernel.globalFilter(typ.second,false) : (what==ACTION) ? kernel.globalAction(typ.second,false) - : (what==FILTER) ? kernel.globalAction(typ.second,false) + /// : (what==FILTER) ? kernel.globalAction(typ.second,false) : 0); // Create the object using the factory method if ( !action ) { diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py index 012d5df2e8ef297bf0bb6fb2d5b11d2fbc8acb32..79e41c61db0ad6ad271d29f6c42c0590ca9a091b 100644 --- a/DDG4/python/DDG4.py +++ b/DDG4/python/DDG4.py @@ -6,6 +6,7 @@ import ROOT # We compile the DDG4 plugin on the fly if it does not exist using the AClick mechanism: def compileAClick(dictionary,g4=True): from ROOT import gInterpreter, gSystem + import os.path dd4hep = os.environ['DD4hep_DIR'] inc = ' -I'+os.environ['ROOTSYS']+'/include -I'+dd4hep+'/include ' lib = ' -L'+dd4hep+'/lib -lDD4hepCore -lDD4hepG4 -lDDSegmentation ' @@ -17,7 +18,10 @@ def compileAClick(dictionary,g4=True): gSystem.AddIncludePath(inc) gSystem.AddLinkedLibs(lib) #####print "Includes: ",gSystem.GetIncludePath(),"\n","Linked libs:",gSystem.GetLinkedLibs() - gInterpreter.ProcessLine('.L '+dictionary+'+') + package = imp.find_module('DDG4') + dic = os.path.dirname(package[1])+os.sep+dictionary + ###print dic + gInterpreter.ProcessLine('.L '+dic+'+') #####gInterpreter.Load('DDG4Dict_C.so') from ROOT import DD4hep as module return module @@ -29,7 +33,7 @@ def _import_class(ns,nam): setattr(current,nam,getattr(scope,nam)) #--------------------------------------------------------------------------- -DD4hep = compileAClick(dictionary='./DDG4Dict.C',g4=True) +DD4hep = compileAClick(dictionary='DDG4Dict.C',g4=True) Sim = DD4hep.Simulation Simulation = DD4hep.Simulation @@ -52,13 +56,39 @@ def _registerGlobalAction(self,action): self.get().registerGlobalAction(Interface.toAction(action)) def _registerGlobalFilter(self,filter): self.get().registerGlobalFilter(Interface.toAction(filter)) +#--------------------------------------------------------------------------- +def _getKernelProperty(self, name): + print '_getKernelProperty:',str(type(self)),name + ret = Interface.getPropertyKernel(self.get(),name) + if ret.status > 0: + return ret.data + elif hasattr(self.get(),name): + return getattr(self.get(),name) + elif hasattr(self,name): + return getattr(self,name) + msg = 'Geant4Kernel::GetProperty [Unhandled]: Cannot access Kernel.'+name + raise exceptions.KeyError(msg) + +#--------------------------------------------------------------------------- +def _setKernelProperty(self, name, value): + #print '_setKernelProperty:',name,value + if Interface.setPropertyKernel(self.get(),name,str(value)): + return + msg = 'Geant4Kernel::SetProperty [Unhandled]: Cannot set Kernel.'+name+' = '+str(value) + raise exceptions.KeyError(msg) Kernel.registerGlobalAction = _registerGlobalAction Kernel.registerGlobalFilter = _registerGlobalFilter +Kernel.__getattr__ = _getKernelProperty +Kernel.__setattr__ = _setKernelProperty + + ActionHandle = Sim.ActionHandle #--------------------------------------------------------------------------- def SensitiveAction(kernel,nam,det): return Interface.createSensitive(kernel,nam,det) #--------------------------------------------------------------------------- +def Action(kernel,nam): return Interface.createAction(kernel,nam) +#--------------------------------------------------------------------------- def Filter(kernel,nam): return Interface.createFilter(kernel,nam) #--------------------------------------------------------------------------- def RunAction(kernel,nam): return Interface.createRunAction(kernel,nam) diff --git a/DDG4/python/DDG4Dict.C b/DDG4/python/DDG4Dict.C index 9f712bb5eb233e5a53c9bf2cab30a14f07a3f1d7..e2c87bb45f1dd1d8c48342b846bb727f14cf5b98 100644 --- a/DDG4/python/DDG4Dict.C +++ b/DDG4/python/DDG4Dict.C @@ -48,7 +48,8 @@ namespace DD4hep { Geant4##x* action; \ x##Handle(Geant4##x* a) : action(a) { if ( action ) action->addRef();} \ x##Handle(const x##Handle& h) : action(h.action) { if ( action ) action->addRef();} \ - ~x##Handle() { if ( action) action->release(); } \ + ~x##Handle() { if ( action) action->release(); } \ + Geant4##x* release() { Geant4##x* tmp = action; action=0; return tmp; } \ operator Geant4##x* () const { return action; } \ Geant4##x* operator->() const { return action; } \ Geant4##x* get() const { return action; } \ @@ -162,6 +163,19 @@ namespace DD4hep { } return 0; } + static PropertyResult getPropertyKernel(Geant4Kernel* kernel, const string& name) { + if ( kernel->hasProperty(name) ) { + return PropertyResult(kernel->property(name).str(),1); + } + return PropertyResult("",0); + } + static int setPropertyKernel(Geant4Kernel* kernel, const string& name, const string& value) { + if ( kernel->hasProperty(name) ) { + kernel->property(name).str(value); + return 1; + } + return 0; + } }; } } diff --git a/DDG4/src/ComponentProperties.cpp b/DDG4/src/ComponentProperties.cpp index eae34d214c1bc276650530b1a5f6107967a60336..816fc35db4f55622e50f9f60fc71c6054ccd5de8 100644 --- a/DDG4/src/ComponentProperties.cpp +++ b/DDG4/src/ComponentProperties.cpp @@ -10,10 +10,12 @@ // Framework include files #include "DD4hep/Printout.h" #include "DDG4/ComponentProperties_inl.h" +#include "DDG4/ToStream.h" // C/C++ include files #include <stdexcept> #include <cstring> +#include <set> using namespace std; using namespace DD4hep; @@ -89,7 +91,7 @@ string Property::str() const { /// Conversion from string value Property& Property::str(const std::string& input) { - if (m_hdl && m_par ) { + if (m_hdl && m_par ) { m_hdl->fromString(m_par,input); return *this; } @@ -192,27 +194,168 @@ void PropertyManager::dump() const { } } -#define DD4HEP_INSTANTIATE_PROPERTY_TYPE1(x) DD4HEP_INSTANTIATE_PROPERTY_TYPE(x); \ - DD4HEP_INSTANTIATE_PROPERTY_TYPE(std::vector<x>); \ - DD4HEP_INSTANTIATE_PROPERTY_TYPE(std::list<x>); \ +#include "XML/Evaluator.h" +namespace DD4hep { XmlTools::Evaluator& g4Evaluator(); } +namespace { static XmlTools::Evaluator& s__eval(DD4hep::g4Evaluator()); } + +static string pre_parse_obj(const string& in) { + string res = ""; + res.reserve(1024); + for(const char* c = in.c_str(); *c; ++c) { + switch(*c) { + case '\'': + return "Bad object representation"; + case ',': + res += "','"; + break; + case '(': + case '[': + res += "['"; + break; + case ')': + case ']': + res += "']"; + break; + default: + res += *c; + break; + } + } + //cout << "Pre-parsed:" << res << endl; + return res; +} + +template <typename TYPE> static int fill(std::vector<TYPE>* p,const std::vector<string>& temp) { + const Grammar<TYPE>& g = Grammar<TYPE>::instance(); + TYPE val; + for(std::vector<string>::const_iterator i=temp.begin(); i != temp.end(); ++i) { + if ( !g.fromString(&val,*i) ) + return 0; + p->push_back(val); + } + return 1; +} +template <typename TYPE> static int fill(std::list<TYPE>* p,const std::vector<string>& temp) { + const Grammar<TYPE>& g = Grammar<TYPE>::instance(); + TYPE val; + for(std::vector<string>::const_iterator i=temp.begin(); i != temp.end(); ++i) { + if ( !g.fromString(&val,*i) ) + return 0; + p->push_back(val); + } + return 1; +} +template <typename TYPE> static int fill(std::set<TYPE>* p,const std::vector<string>& temp) { + const Grammar<TYPE>& g = Grammar<TYPE>::instance(); + TYPE val; + for(std::vector<string>::const_iterator i=temp.begin(); i != temp.end(); ++i) { + if ( !g.fromString(&val,*i) ) + return 0; + p->insert(val); + } + return 1; +} + +template <typename TYPE> static int eval(TYPE* p, const string& str) { +#ifdef DD4HEP_USE_BOOST + std::vector<string> buff; + int sc = Parsers::parse(buff,str); + if ( sc ) { + return fill(p,buff); + } + else { + TYPE temp; + std::string temp_str = pre_parse_obj(str); + sc = Parsers::parse(temp,temp_str); + if ( sc ) { + *p = temp; + return 1; + } + buff.clear(); + sc = Parsers::parse(buff,temp_str); + if ( sc ) { + return fill(p,buff); + } + } +#else + if ( p && str.empty() ) return 0; +#endif + return 0; +} + +template <typename T> static int _eval(T* p, string s) { + size_t idx = s.find("(int)"); + if (idx != string::npos) + s.erase(idx, 5); + while (s[0] == ' ') + s.erase(0, 1); + double result = s__eval.evaluate(s.c_str()); + if (s__eval.status() != XmlTools::Evaluator::OK) { + return 0; + } + *p = (T)result; + return 1; +} + +static int _eval(ROOT::Math::XYZPoint* p, const string& str) +{ return Grammar<ROOT::Math::XYZPoint>::instance().fromString(p,pre_parse_obj(str)); } +static int _eval(ROOT::Math::XYZVector* p, const string& str) +{ return Grammar<ROOT::Math::XYZVector>::instance().fromString(p,pre_parse_obj(str)); } +static int _eval(ROOT::Math::PxPyPzEVector* p, const string& str) +{ return Grammar<ROOT::Math::PxPyPzEVector>::instance().fromString(p,pre_parse_obj(str)); } + +#define DD4HEP_INSTANTIATE_PROPERTY_EVALUATE(x) \ + template<> int Grammar<x >::evaluate(void* p, const string& v) const { return _eval((x*)p,v); } + +#define DD4HEP_INSTANTIATE_PROPERTY_EVALUATE1(x) DD4HEP_INSTANTIATE_PROPERTY_EVALUATE(x)\ + template<> int Grammar<vector<x> >::evaluate(void* p,const string& v)const{return eval((std::vector<x>*)p,v);} \ + template<> int Grammar<list<x> >::evaluate(void* p,const string& v)const{return eval((std::list<x>*)p,v);} \ + template<> int Grammar<set<x> >::evaluate(void* p,const string& v)const{return eval((std::set<x>*)p,v);} + +#define DD4HEP_INSTANTIATE_PROPERTY_TYPE1(x) \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE(x); \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE(std::vector<x>); \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE(std::list<x>); \ DD4HEP_INSTANTIATE_PROPERTY_TYPE(std::set<x>) #define DD4HEP_INSTANTIATE_PROPERTY_TYPE2(x) \ DD4HEP_INSTANTIATE_PROPERTY_TYPE1(x); \ DD4HEP_INSTANTIATE_PROPERTY_TYPE1(unsigned x) +// Macros for evaluated properties: + +#define DD4HEP_INSTANTIATE_PROPERTY_TYPE0E(x) \ + DD4HEP_INSTANTIATE_PROPERTY_EVALUATE(x) \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE(x) + +#define DD4HEP_INSTANTIATE_PROPERTY_TYPE1E(x) \ + DD4HEP_INSTANTIATE_PROPERTY_EVALUATE1(x) \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE1(x) + +#define DD4HEP_INSTANTIATE_PROPERTY_TYPE2E(x) \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE1E(x); \ + DD4HEP_INSTANTIATE_PROPERTY_TYPE1E(unsigned x) + namespace DD4hep { - DD4HEP_INSTANTIATE_PROPERTY_TYPE2(char); - DD4HEP_INSTANTIATE_PROPERTY_TYPE2(short); - DD4HEP_INSTANTIATE_PROPERTY_TYPE2(int); - DD4HEP_INSTANTIATE_PROPERTY_TYPE2(long); - DD4HEP_INSTANTIATE_PROPERTY_TYPE2(long long); - - DD4HEP_INSTANTIATE_PROPERTY_TYPE1(bool); - DD4HEP_INSTANTIATE_PROPERTY_TYPE1(float); - DD4HEP_INSTANTIATE_PROPERTY_TYPE1(double); + + DD4HEP_INSTANTIATE_PROPERTY_TYPE2E(char); + DD4HEP_INSTANTIATE_PROPERTY_TYPE2E(short); + DD4HEP_INSTANTIATE_PROPERTY_TYPE2E(int); + DD4HEP_INSTANTIATE_PROPERTY_TYPE2E(long); + DD4HEP_INSTANTIATE_PROPERTY_TYPE2E(long long); + + DD4HEP_INSTANTIATE_PROPERTY_TYPE1E(bool); + DD4HEP_INSTANTIATE_PROPERTY_TYPE1E(float); + DD4HEP_INSTANTIATE_PROPERTY_TYPE1E(double); + + // STL objects DD4HEP_INSTANTIATE_PROPERTY_TYPE1(string); typedef map<string, int> map_string_int; - DD4HEP_INSTANTIATE_PROPERTY_TYPE (map_string_int); + DD4HEP_INSTANTIATE_PROPERTY_TYPE(map_string_int); + + // ROOT::Math Object instances + DD4HEP_INSTANTIATE_PROPERTY_TYPE0E(ROOT::Math::XYZPoint); + DD4HEP_INSTANTIATE_PROPERTY_TYPE0E(ROOT::Math::XYZVector); + DD4HEP_INSTANTIATE_PROPERTY_TYPE0E(ROOT::Math::PxPyPzEVector); } diff --git a/DDG4/src/Geant4Action.cpp b/DDG4/src/Geant4Action.cpp index f2c8617bb315a3e4844494dfa1407fc2ff7d4f27..ebcf23464ada782359f5c5e53be96c63fe656a14 100644 --- a/DDG4/src/Geant4Action.cpp +++ b/DDG4/src/Geant4Action.cpp @@ -21,8 +21,8 @@ using namespace std; using namespace DD4hep; using namespace DD4hep::Simulation; -TypeName TypeName::split(const string& type_name) { - size_t idx = type_name.find("/"); +TypeName TypeName::split(const string& type_name, const std::string& delim) { + size_t idx = type_name.find(delim); string typ = type_name, nam = type_name; if (idx != string::npos) { typ = type_name.substr(0, idx); @@ -31,6 +31,10 @@ TypeName TypeName::split(const string& type_name) { return TypeName(typ, nam); } +TypeName TypeName::split(const string& type_name) { + return split(type_name,"/"); +} + /// Default constructor Geant4TrackInformation::Geant4TrackInformation() : G4VUserTrackInformation(), m_flags(0) { @@ -197,32 +201,55 @@ void Geant4Action::except(const string& fmt, ...) const { throw runtime_error(err); } +/// Abort Geant4 Run by throwing a G4Exception with type RunMustBeAborted +void Geant4Action::abortRun(const string& exception, const string& fmt, ...) const { + string desc, typ = typeinfoName(typeid(*this)); + string issuer = name()+" ["+typ+"]"; + va_list args; + va_start(args, fmt); + desc = DD4hep::format("*** Geant4Action:", fmt, args); + va_end(args); + G4Exception(issuer.c_str(),exception.c_str(),RunMustBeAborted,desc.c_str()); + //throw runtime_error(issuer+"> "+desc); +} + /// Access to the main run action sequence from the kernel object Geant4RunActionSequence& Geant4Action::runAction() const { - return m_context->runAction(); + return m_context->kernel().runAction(); } /// Access to the main event action sequence from the kernel object Geant4EventActionSequence& Geant4Action::eventAction() const { - return m_context->eventAction(); + return m_context->kernel().eventAction(); } /// Access to the main stepping action sequence from the kernel object Geant4SteppingActionSequence& Geant4Action::steppingAction() const { - return m_context->steppingAction(); + return m_context->kernel().steppingAction(); } /// Access to the main tracking action sequence from the kernel object Geant4TrackingActionSequence& Geant4Action::trackingAction() const { - return m_context->trackingAction(); + return m_context->kernel().trackingAction(); } /// Access to the main stacking action sequence from the kernel object Geant4StackingActionSequence& Geant4Action::stackingAction() const { - return m_context->stackingAction(); + return m_context->kernel().stackingAction(); } /// Access to the main generator action sequence from the kernel object Geant4GeneratorActionSequence& Geant4Action::generatorAction() const { - return m_context->generatorAction(); + return m_context->kernel().generatorAction(); +} + +/// Access to the Track Manager from the kernel object +Geant4MonteCarloTruth& Geant4Action::mcTruthMgr() const { + return m_context->kernel().mcTruthMgr(); } + +/// Access to the MC record manager from the kernel object +Geant4MonteCarloRecordManager& Geant4Action::mcRecordMgr() const { + return m_context->kernel().mcRecordMgr(); +} + diff --git a/DDG4/src/Geant4Call.cpp b/DDG4/src/Geant4Call.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6b53c9a322b8516c46b7f1c50fdf0454572bf9c7 --- /dev/null +++ b/DDG4/src/Geant4Call.cpp @@ -0,0 +1,15 @@ +// $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/Geant4Call.h" + +/// Default destructor (keep here to avoid weak linkage to vtable) +DD4hep::Simulation::Geant4Call::~Geant4Call() { +} + diff --git a/DDG4/src/Geant4Context.cpp b/DDG4/src/Geant4Context.cpp index 876bfc8a00fb2e6c17f00a6fdc34c81ec7fb8e8c..3f6c668bd3da213dd5b29603c89f69aabb8badd6 100644 --- a/DDG4/src/Geant4Context.cpp +++ b/DDG4/src/Geant4Context.cpp @@ -18,7 +18,7 @@ using namespace DD4hep::Simulation; /// Default constructor Geant4Context::Geant4Context(Geant4Kernel* kernel) - : m_kernel(kernel), m_trackMgr(0) { + : m_kernel(kernel) { } /// Default destructor @@ -37,6 +37,11 @@ G4VTrajectory* Geant4Context::createTrajectory(const G4Track* /* track */) const throw runtime_error(err); } +/// Access the tracking manager +G4TrackingManager* Geant4Context::trackMgr() const { + return m_kernel->trackMgr(); +} + /// Access to the main run action sequence from the kernel object Geant4RunActionSequence& Geant4Context::runAction() const { return m_kernel->runAction(); @@ -71,3 +76,14 @@ Geant4GeneratorActionSequence& Geant4Context::generatorAction() const { Geant4SensDetSequences& Geant4Context::sensitiveActions() const { return m_kernel->sensitiveActions(); } + +/// Access to the Track Manager from the kernel object +Geant4MonteCarloTruth& Geant4Context::mcTruthMgr() const { + return m_kernel->mcTruthMgr(); +} + +/// Access to the MC record manager from the kernel object +Geant4MonteCarloRecordManager& Geant4Context::mcRecordMgr() const { + return m_kernel->mcRecordMgr(); +} + diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp index de6c5bbe2b1ac33cde977cef94a0618571ee28d8..d09f88eb0db555907294864dcd7042a319ceed66 100644 --- a/DDG4/src/Geant4Converter.cpp +++ b/DDG4/src/Geant4Converter.cpp @@ -481,7 +481,7 @@ void* Geant4Converter::handleVolume(const string& name, const TGeoVolume* volume throw runtime_error("G4Converter: No Geant4 material present for volume:" + n); } if (user_limits) { - printout(DEBUG, "++ Volume + Apply LIMITS settings:%-24s to volume %s.", lim.name(), _v.name()); + printout(DEBUG, "Geant4Converter", "++ Volume + Apply LIMITS settings:%-24s to volume %s.", lim.name(), _v.name()); } vol = new G4LogicalVolume(solid, medium, n, 0, sd, user_limits); if (region) { @@ -869,6 +869,8 @@ Geant4Converter& Geant4Converter::create(DetElement top) { mat = m_lcdd.material("Silicon"); handleMaterial(mat.name(), mat.ptr()); + //setPrintLevel(VERBOSE); + handle(this, geo.volumes, &Geant4Converter::collectVolume); handle(this, geo.solids, &Geant4Converter::handleSolid); printout(INFO, "Geant4Converter", "++ Handled %ld solids.", geo.solids.size()); diff --git a/DDG4/src/Geant4EventAction.cpp b/DDG4/src/Geant4EventAction.cpp index d5897a9fdd63eab562a9d10f202579b453a7a431..006b9655b2dab55af36a66955f36c4f63dc59400 100644 --- a/DDG4/src/Geant4EventAction.cpp +++ b/DDG4/src/Geant4EventAction.cpp @@ -47,6 +47,7 @@ Geant4EventActionSequence::~Geant4EventActionSequence() { m_actors.clear(); m_begin.clear(); m_end.clear(); + m_final.clear(); InstanceCount::decrement(this); } @@ -70,4 +71,5 @@ void Geant4EventActionSequence::begin(const G4Event* event) { void Geant4EventActionSequence::end(const G4Event* event) { m_end(event); m_actors(&Geant4EventAction::end, event); + m_final(event); } diff --git a/DDG4/src/Geant4Exec.cpp b/DDG4/src/Geant4Exec.cpp index d388b7b136a529691635cf7cf5754ec1b309a848..9e85ef88a2fe6778edc79e5ea0b284e392158ee9 100644 --- a/DDG4/src/Geant4Exec.cpp +++ b/DDG4/src/Geant4Exec.cpp @@ -29,8 +29,11 @@ #include "DDG4/Geant4GeneratorAction.h" #include "DDG4/Geant4PhysicsList.h" #include "DDG4/Geant4Kernel.h" +#include "DDG4/Geant4UIManager.h" #include <memory> +#include <stdexcept> + /* * DD4hep namespace declaration */ @@ -40,30 +43,29 @@ namespace DD4hep { * Simulation namespace declaration */ namespace Simulation { - namespace { - template <typename T> struct _Seq { - typedef _Seq<T> Base; - T* m_sequence; - _Seq() - : m_sequence(0) { - } - _Seq(T* seq) { - _aquire(seq); - } - virtual ~_Seq() { - _release(); - } - void _aquire(T* s) { - InstanceCount::increment(this); - m_sequence = s; - m_sequence->addRef(); - } - void _release() { - releasePtr(m_sequence); - InstanceCount::decrement(this); - } - }; - } + template <typename T> struct SequenceHdl { + typedef SequenceHdl<T> Base; + T* m_sequence; + SequenceHdl() + : m_sequence(0) { + } + SequenceHdl(T* seq) { + _aquire(seq); + } + virtual ~SequenceHdl() { + _release(); + } + void _aquire(T* s) { + InstanceCount::increment(this); + m_sequence = s; + m_sequence->addRef(); + } + void _release() { + releasePtr(m_sequence); + InstanceCount::decrement(this); + } + }; + /** @class Geant4UserRunAction * * Concrete implementation of the Geant4 run action @@ -71,7 +73,7 @@ namespace DD4hep { * @author M.Frank * @version 1.0 */ - struct Geant4UserRunAction : public G4UserRunAction, public _Seq<Geant4RunActionSequence> { + struct Geant4UserRunAction : public G4UserRunAction, public SequenceHdl<Geant4RunActionSequence> { /// Standard constructor Geant4UserRunAction(Geant4RunActionSequence* seq) : Base(seq) { @@ -96,7 +98,7 @@ namespace DD4hep { * @author M.Frank * @version 1.0 */ - struct Geant4UserEventAction : public G4UserEventAction, public _Seq<Geant4EventActionSequence> { + struct Geant4UserEventAction : public G4UserEventAction, public SequenceHdl<Geant4EventActionSequence> { /// Standard constructor Geant4UserEventAction(Geant4EventActionSequence* seq) : Base(seq) { @@ -121,7 +123,7 @@ namespace DD4hep { * @author M.Frank * @version 1.0 */ - struct Geant4UserTrackingAction : public G4UserTrackingAction, public _Seq<Geant4TrackingActionSequence> { + struct Geant4UserTrackingAction : public G4UserTrackingAction, public SequenceHdl<Geant4TrackingActionSequence> { /// Standard constructor Geant4UserTrackingAction(Geant4TrackingActionSequence* seq) : Base(seq) { @@ -131,11 +133,13 @@ namespace DD4hep { } /// Pre-track action callback virtual void PreUserTrackingAction(const G4Track* trk) { + m_sequence->context()->kernel().setTrackMgr(fpTrackingManager); m_sequence->begin(trk); } /// Post-track action callback virtual void PostUserTrackingAction(const G4Track* trk) { m_sequence->end(trk); + m_sequence->context()->kernel().setTrackMgr(0); } }; @@ -146,7 +150,7 @@ namespace DD4hep { * @author M.Frank * @version 1.0 */ - struct Geant4UserStackingAction : public G4UserStackingAction, public _Seq<Geant4StackingActionSequence> { + struct Geant4UserStackingAction : public G4UserStackingAction, public SequenceHdl<Geant4StackingActionSequence> { /// Standard constructor Geant4UserStackingAction(Geant4StackingActionSequence* seq) : Base(seq) { @@ -171,7 +175,7 @@ namespace DD4hep { * @author M.Frank * @version 1.0 */ - struct Geant4UserGeneratorAction : public G4VUserPrimaryGeneratorAction, public _Seq<Geant4GeneratorActionSequence> { + struct Geant4UserGeneratorAction : public G4VUserPrimaryGeneratorAction, public SequenceHdl<Geant4GeneratorActionSequence> { /// Standard constructor Geant4UserGeneratorAction(Geant4GeneratorActionSequence* seq) : G4VUserPrimaryGeneratorAction(), Base(seq) { @@ -192,7 +196,7 @@ namespace DD4hep { * @author M.Frank * @version 1.0 */ - struct Geant4UserSteppingAction : public G4UserSteppingAction, public _Seq<Geant4SteppingActionSequence> { + struct Geant4UserSteppingAction : public G4UserSteppingAction, public SequenceHdl<Geant4SteppingActionSequence> { /// Standard constructor Geant4UserSteppingAction(Geant4SteppingActionSequence* seq) : Base(seq) { @@ -222,18 +226,6 @@ using namespace DD4hep::Simulation; // Geant4 include files #include "G4RunManager.hh" -//-- -#define G4VIS_USE_OPENGL -#include "G4UImanager.hh" -#include "G4UIsession.hh" -#ifdef G4VIS_USE_OPENGLX -#include "G4OpenGLImmediateX.hh" -#include "G4OpenGLStoredX.hh" -#endif -#include "G4VisManager.hh" -#include "G4VisExecutive.hh" -#include "G4UIExecutive.hh" - /// Configure the simulation int Geant4Exec::configure(Geant4Kernel& kernel) { CLHEP::HepRandom::setTheEngine(new CLHEP::RanecuEngine); @@ -282,6 +274,11 @@ int Geant4Exec::configure(Geant4Kernel& kernel) { Geant4UserEventAction* action = new Geant4UserEventAction(kernel.eventAction(false)); runManager.SetUserAction(action); } + // Set the tracking action sequence + if (kernel.trackingAction(false)) { + Geant4UserTrackingAction* action = new Geant4UserTrackingAction(kernel.trackingAction(false)); + runManager.SetUserAction(action); + } // Set the stepping action sequence if (kernel.steppingAction(false)) { Geant4UserSteppingAction* action = new Geant4UserSteppingAction(kernel.steppingAction(false)); @@ -307,50 +304,23 @@ int Geant4Exec::initialize(Geant4Kernel& kernel) { } /// Run the simulation -int Geant4Exec::run(Geant4Kernel&) { - bool is_visual = false; - string gui_setup, vis_setup, gui_type = "tcsh"; - - // Initialize visualization - G4VisManager* visManager = 0; - if (is_visual) { - visManager = new G4VisExecutive(); -#ifdef G4VIS_USE_OPENGLX - //visManager->RegisterGraphicsSystem (new G4OpenGLImmediateX()); - //visManager->RegisterGraphicsSystem (new G4OpenGLStoredX()); -#endif - visManager->Initialize(); - } - - // Get the pointer to the User Interface manager - G4UImanager* UImanager = G4UImanager::GetUIpointer(); - G4String command = "/control/execute run.mac"; - cout << "++++++++++++++++++++++++++++ executing command:" << command << endl; - UImanager->ApplyCommand(command); - - G4UIExecutive* ui = 0; - if (!gui_type.empty()) { // interactive mode : define UI session - const char* args[] = { "cmd" }; - ui = new G4UIExecutive(1, (char**) args); //,gui_type); - if (is_visual && !vis_setup.empty()) { - UImanager->ApplyCommand("/control/execute vis.mac"); - cout << "++++++++++++++++++++++++++++ executed vis.mac" << endl; - } - if (ui->IsGUI()) { - if (!gui_setup.empty()) { - UImanager->ApplyCommand("/control/execute " + gui_setup); - cout << "++++++++++++++++++++++++++++ executed gui.mac" << endl; +int Geant4Exec::run(Geant4Kernel& kernel) { + Property& p = kernel.property("UI"); + string value = p.value<string>(); + if ( !value.empty() ) { + Geant4Action* ui = kernel.globalAction(value); + if ( ui ) { + Geant4Call* c = dynamic_cast<Geant4Call*>(ui); + if ( c ) { + (*c)(0); + return 1; } + ui->except("++ Geant4Exec: Failed to start UI interface."); } - ui->SessionStart(); - delete ui; + throw runtime_error(format("Geant4Exec","++ Failed to locate UI interface %s.",value.c_str())); } - // Job termination - // Free the store: user actions, physics_list and detector_description are - // owned and deleted by the run manager, so they should not - // be deleted in the main() program ! - if (visManager) - delete visManager; + long nevt = kernel.property("NumEvent").value<long>(); + kernel.runManager().BeamOn(nevt); return 1; } diff --git a/DDG4/src/Geant4Field.cpp b/DDG4/src/Geant4Field.cpp index 6e114cf786661dce5805112aaf83b797ba4a6b42..6a03a4a4eb5d5ec632126f50296f6f721c7f8254 100644 --- a/DDG4/src/Geant4Field.cpp +++ b/DDG4/src/Geant4Field.cpp @@ -8,6 +8,8 @@ //==================================================================== // Framework include files #include "DDG4/Geant4Field.h" +#include "DD4hep/TGeoUnits.h" +#include "CLHEP/Units/SystemOfUnits.h" using namespace DD4hep::Simulation; @@ -16,7 +18,13 @@ G4bool Geant4Field::DoesFieldChangeEnergy() const { } void Geant4Field::GetFieldValue(const double pos[4], double *field) const { - field[0] = field[1] = field[2] = 0.0; - return m_field.magneticField(pos, field); + static const double fac1 = tgeo::mm/CLHEP::mm; + static const double fac2 = CLHEP::tesla/tgeo::tesla; + double p[3] = {pos[0]*fac1, pos[1]*fac1, pos[2]*fac1}; // Convert from CLHEP units to tgeo units + field[0] = field[1] = field[2] = 0.0; // Reset field vector + m_field.magneticField(p, field); + field[0] *= fac2; // Convert from tgeo units to CLHEP units + field[1] *= fac2; + field[2] *= fac2; + //::printf("Pos: %7.4f %7.4f %7.4f --> %9g %9g %9g\n",p[0],p[1],p[2],field[0],field[1],field[2]); } - diff --git a/DDG4/src/Geant4Kernel.cpp b/DDG4/src/Geant4Kernel.cpp index 08016760d5631562e86245268dfc2567e77074bd..3cab6f7e74e8a9ec6e79a971e0553b613da8460a 100644 --- a/DDG4/src/Geant4Kernel.cpp +++ b/DDG4/src/Geant4Kernel.cpp @@ -25,6 +25,8 @@ #include "DDG4/Geant4GeneratorAction.h" #include "DDG4/Geant4SensDetAction.h" #include "DDG4/Geant4DetectorConstruction.h" +#include "DDG4/Geant4MonteCarloRecordManager.h" +#include "DDG4/Geant4TrackPersistency.h" // Geant4 include files #include "G4RunManager.hh" @@ -47,6 +49,14 @@ Geant4Kernel::PhaseSelector::PhaseSelector(const PhaseSelector& c) : m_kernel(c.m_kernel) { } +/// Assignment operator +Geant4Kernel::PhaseSelector& Geant4Kernel::PhaseSelector::operator=(const PhaseSelector& c) { + if ( this != &c ) { + m_kernel = c.m_kernel; + } + return *this; +} + /// Phase access to the map Geant4ActionPhase& Geant4Kernel::PhaseSelector::operator[](const std::string& name) const { Geant4ActionPhase* phase = m_kernel->getPhase(name); @@ -58,8 +68,10 @@ Geant4ActionPhase& Geant4Kernel::PhaseSelector::operator[](const std::string& na /// Standard constructor Geant4Kernel::Geant4Kernel(LCDD& lcdd) - : m_context(0), m_runManager(0), m_generatorAction(0), m_runAction(0), m_eventAction(0), m_trackingAction(0), m_steppingAction( - 0), m_stackingAction(0), m_sensDetActions(0), m_physicsList(0), m_lcdd(lcdd), phase(this) { + : m_context(0), m_runManager(0), m_generatorAction(0), m_runAction(0), m_eventAction(0), + m_trackingAction(0), m_steppingAction(0), m_stackingAction(0), m_sensDetActions(0), + m_physicsList(0), m_mcTruthMgr(0), m_mcRecordMgr(0), + m_lcdd(lcdd), phase(this) { #if 0 registerSequence(m_runAction, "RunAction"); registerSequence(m_eventAction, "EventAction"); @@ -72,6 +84,8 @@ Geant4Kernel::Geant4Kernel(LCDD& lcdd) m_context = new Geant4Context(this); m_lcdd.addExtension < Geant4Kernel > (this); + declareProperty("UI",m_uiName); + declareProperty("NumEvents",m_numEvent = 10); m_controlName = "/ddg4/"; m_control = new G4UIdirectory(m_controlName.c_str()); m_control->SetGuidance("Control for all named Geant4 actions"); @@ -83,7 +97,9 @@ Geant4Kernel::~Geant4Kernel() { destroyPhases(); for_each(m_globalFilters.begin(), m_globalFilters.end(), releaseObjects(m_globalFilters)); for_each(m_globalActions.begin(), m_globalActions.end(), releaseObjects(m_globalActions)); - deletePtr (m_runManager); + deletePtr (m_runManager); + deletePtr (m_mcTruthMgr); + releasePtr (m_mcRecordMgr); releasePtr (m_physicsList); releasePtr (m_stackingAction); releasePtr (m_steppingAction); @@ -91,8 +107,8 @@ Geant4Kernel::~Geant4Kernel() { releasePtr (m_eventAction); releasePtr (m_generatorAction); releasePtr (m_runAction); - deletePtr (m_sensDetActions); - deletePtr (m_context); + deletePtr (m_sensDetActions); + deletePtr (m_context); m_lcdd.destroyInstance(); InstanceCount::decrement(this); } @@ -103,6 +119,16 @@ Geant4Kernel& Geant4Kernel::instance(LCDD& lcdd) { return obj; } +/// Check property for existence +bool Geant4Kernel::hasProperty(const std::string& name) const { + return m_properties.exists(name); +} + +/// Access single property +DD4hep::Property& Geant4Kernel::property(const std::string& name) { + return properties()[name]; +} + /// Accessof the Geant4Kernel object from the LCDD reference extension (if present and registered) Geant4Kernel& Geant4Kernel::access(LCDD& lcdd) { Geant4Kernel* kernel = lcdd.extension<Geant4Kernel>(); @@ -154,7 +180,9 @@ void Geant4Kernel::terminate() { destroyPhases(); for_each(m_globalFilters.begin(), m_globalFilters.end(), releaseObjects(m_globalFilters)); for_each(m_globalActions.begin(), m_globalActions.end(), releaseObjects(m_globalActions)); - deletePtr (m_runManager); + deletePtr (m_runManager); + deletePtr (m_mcTruthMgr); + releasePtr (m_mcRecordMgr); releasePtr (m_physicsList); releasePtr (m_stackingAction); releasePtr (m_steppingAction); @@ -162,8 +190,8 @@ void Geant4Kernel::terminate() { releasePtr (m_eventAction); releasePtr (m_generatorAction); releasePtr (m_runAction); - deletePtr (m_sensDetActions); - deletePtr (m_context); + deletePtr (m_sensDetActions); + deletePtr (m_context); //return *this; } @@ -188,6 +216,8 @@ Geant4Kernel& Geant4Kernel::registerGlobalAction(Geant4Action* action) { if (i == m_globalActions.end()) { action->addRef(); m_globalActions[nam] = action; + printout(INFO,"Geant4Kernel","++ Registered global action %s of type %s", + nam.c_str(),typeinfoName(typeid(*action)).c_str()); return *this; } throw runtime_error(format("Geant4Kernel", "DDG4: The action '%s' is already globally " @@ -350,3 +380,31 @@ Geant4PhysicsListActionSequence* Geant4Kernel::physicsList(bool create) { return m_physicsList; } +/// Access to the Track Manager from the kernel object +Geant4MonteCarloTruth& Geant4Kernel::mcTruthMgr() { + if ( m_mcTruthMgr ) return *m_mcTruthMgr; + // If not present, check if the action is registered. + Geant4Action* a = globalAction("MonteCarloTruthHandler",false); + if ( 0 != a ) { + m_mcTruthMgr = dynamic_cast<Geant4MonteCarloTruth*>(a); + if ( m_mcTruthMgr ) return *m_mcTruthMgr; + } + // No action registered to handle monte carlo truth. This is fatal + throw runtime_error(format("Geant4Kernel", "DDG4: No Geant4MonteCarloTruth defined. " + "Geant4 monte carlo information cannot be saved!")); +} + +/// Access to the MC record manager from the kernel object +Geant4MonteCarloRecordManager& Geant4Kernel::mcRecordMgr() { + if ( m_mcRecordMgr ) return *m_mcRecordMgr; + // If not present, check if the action is registered. + Geant4Action* a = globalAction("MonteCarloRecordManager",false); + if ( 0 != a ) { + m_mcRecordMgr = dynamic_cast<Geant4MonteCarloRecordManager*>(a); + if ( m_mcRecordMgr ) return *m_mcRecordMgr; + } + // No action registered to save tracks. This is fatal + throw runtime_error(format("Geant4Kernel", "DDG4: No MonteCarloRecordManager defined. " + "Geant4 tracks cannot be saved!")); +} + diff --git a/DDG4/src/Geant4MonteCarloRecordManager.cpp b/DDG4/src/Geant4MonteCarloRecordManager.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a98a613e492425a15ca87c0f7bf03967db3a17f8 --- /dev/null +++ b/DDG4/src/Geant4MonteCarloRecordManager.cpp @@ -0,0 +1,54 @@ +// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +// Framework include files +#include "DD4hep/Printout.h" +#include "DDG4/Geant4TrackPersistency.h" +#include "DDG4/Geant4TrackingAction.h" +#include "DDG4/Geant4MonteCarloRecordManager.h" +#include "DDG4/Geant4TrackHandler.h" +#include "DDG4/Geant4StepHandler.h" + +#include "G4TrackingManager.hh" +#include "G4Track.hh" + +using namespace DD4hep::Simulation; + +/// Standard constructor +Geant4MonteCarloRecordManager::Geant4MonteCarloRecordManager(Geant4Context* context, const std::string& nam) + : Geant4Action(context,nam) +{ + declareProperty("Collect",m_collectInfo=true); +} + +/// Default destructor +Geant4MonteCarloRecordManager::~Geant4MonteCarloRecordManager() { +} + +/// Save G4Track data +void Geant4MonteCarloRecordManager::save(const Geant4TrackPersistency::TrackInfo& trk) { + if ( m_collectInfo ) { + Geant4TrackHandler track(trk.track); +#if 0 + // const G4Step* step = track.step(); + std::string proc = track.creatorName(); + float ene = trk.initialEnergy; + float loss = ene-track.energy(); + G4ThreeVector pos = track.position(); + G4ThreeVector org = track.vertex(); + + info("Track ID=%4d Parent:%4d E:%7.2f MeV Loss:%7.2f MeV " + "%-6s %-8s L:%7.2f Org:(%8.2f %8.2f %8.2f) Pos:(%8.2f %8.2f %8.2f)", + track.id(), track.parent(), ene, loss, + track.name().c_str(), ("["+proc+"]").c_str(), track.length(), + org.x(),org.y(),org.z(), + pos.x(),pos.y(),pos.z()); +#endif + } +} + diff --git a/DDG4/src/Geant4MonteCarloTruth.cpp b/DDG4/src/Geant4MonteCarloTruth.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b597e40c3f201bb92c5d50051b374b09fd0fc405 --- /dev/null +++ b/DDG4/src/Geant4MonteCarloTruth.cpp @@ -0,0 +1,24 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== + +// Framework include files +#include "DD4hep/InstanceCount.h" +#include "DDG4/Geant4MonteCarloTruth.h" + +using namespace DD4hep::Simulation; + +/// Standard action constructor +Geant4MonteCarloTruth::Geant4MonteCarloTruth() { + InstanceCount::increment(this); +} + +/// Default destructor +Geant4MonteCarloTruth::~Geant4MonteCarloTruth() { + InstanceCount::decrement(this); +} diff --git a/DDG4/src/Geant4ParticleGun.cpp b/DDG4/src/Geant4ParticleGun.cpp index 0dfd09317424fdd20aa64421f04eb83b7a3d91e2..da647f25c062da7d1b82e93701db8aa52fad78c9 100644 --- a/DDG4/src/Geant4ParticleGun.cpp +++ b/DDG4/src/Geant4ParticleGun.cpp @@ -17,23 +17,24 @@ // C/C++ include files #include <stdexcept> +#include <limits> +#include <cmath> +using namespace std; using namespace DD4hep::Simulation; /// Standard constructor -Geant4ParticleGun::Geant4ParticleGun(Geant4Context* context, const std::string& name) - : Geant4GeneratorAction(context, name), m_particle(0), m_gun(0) { +Geant4ParticleGun::Geant4ParticleGun(Geant4Context* context, const string& name) + : Geant4GeneratorAction(context, name), m_position(0,0,0), m_direction(1,1,0.3), + m_particle(0), m_gun(0) +{ InstanceCount::increment(this); m_needsControl = true; + declareProperty("particle", m_particleName = "e-"); declareProperty("energy", m_energy = 50 * MeV); declareProperty("multiplicity", m_multiplicity = 1); - declareProperty("pos_x", m_position.x = 0); - declareProperty("pos_y", m_position.y = 0); - declareProperty("pos_z", m_position.z = 0); - declareProperty("direction_x", m_direction.x = 1); - declareProperty("direction_y", m_direction.y = 1); - declareProperty("direction_z", m_direction.z = 0.3); - declareProperty("particle", m_particleName = "e-"); + declareProperty("position", m_position); + declareProperty("direction", m_direction); } /// Default destructor @@ -52,12 +53,16 @@ void Geant4ParticleGun::operator()(G4Event* event) { G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); m_particle = particleTable->FindParticle(m_particleName); if (0 == m_particle) { - throw std::runtime_error("Bad particle type!"); + throw runtime_error("Bad particle type!"); } } + double r = m_direction.R(), eps = numeric_limits<float>::epsilon(); + if ( r > eps && std::fabs(r-1.0) > eps ) { + m_direction.SetXYZ(m_direction.X()/r, m_direction.Y()/r, m_direction.Z()/r); + } m_gun->SetParticleDefinition(m_particle); m_gun->SetParticleEnergy(m_energy); - m_gun->SetParticleMomentumDirection(G4ThreeVector(m_direction.x, m_direction.y, m_direction.z)); - m_gun->SetParticlePosition(G4ThreeVector(m_position.x, m_position.y, m_position.z)); + m_gun->SetParticleMomentumDirection(G4ThreeVector(m_direction.X(), m_direction.Y(), m_direction.Z())); + m_gun->SetParticlePosition(G4ThreeVector(m_position.X(), m_position.Y(), m_position.Z())); m_gun->GeneratePrimaryVertex(event); } diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp index 0350b3e855b033f553f9a5be076babccdc1938a9..37adf5d6bacd32b1d0a211ebb08f918bb7fe9f98 100644 --- a/DDG4/src/Geant4SensDetAction.cpp +++ b/DDG4/src/Geant4SensDetAction.cpp @@ -16,7 +16,9 @@ #include "DDG4/Geant4SensDetAction.h" #include "DDG4/Geant4StepHandler.h" #include "DDG4/Geant4VolumeManager.h" +#include "DDG4/Geant4MonteCarloTruth.h" +#include <G4Step.hh> #include <G4SDManager.hh> #include <G4VSensitiveDetector.hh> @@ -168,6 +170,16 @@ bool Geant4Sensitive::process(G4Step* /* step */, G4TouchableHistory* /* history void Geant4Sensitive::clear(G4HCofThisEvent* /* HCE */) { } +/// Mark the track to be kept for MC truth propagation during hit processing +void Geant4Sensitive::mark(const G4Track* track) const { + mcTruthMgr().mark(track); +} + +/// Mark the track of this step to be kept for MC truth propagation during hit processing +void Geant4Sensitive::mark(const G4Step* step) const { + mcTruthMgr().mark(step); +} + /// Returns the volumeID of the sensitive volume corresponding to the step long long int Geant4Sensitive::volumeID(G4Step* s) { Geant4StepHandler step(s); diff --git a/DDG4/src/Geant4TestActions.cpp b/DDG4/src/Geant4TestActions.cpp index fb8c2b77b34d7d05d2f5810e0a5f5d25a3cbd00c..4db4b11c7ff043102a06032c697cd7552b464fe9 100644 --- a/DDG4/src/Geant4TestActions.cpp +++ b/DDG4/src/Geant4TestActions.cpp @@ -59,7 +59,8 @@ void Geant4TestRunAction::begin(const G4Run* run) { } /// End-of-run callback void Geant4TestRunAction::end(const G4Run* run) { - printout(INFO, name(), "%s> calling end(run_id=%d, num_event=%d)", m_type.c_str(), run->GetRunID(), run->GetNumberOfEvent()); + printout(INFO, name(), "%s> calling end(run_id=%d, num_event=%d)", + m_type.c_str(), run->GetRunID(), run->GetNumberOfEvent()); } /// begin-of-event callback void Geant4TestRunAction::beginEvent(const G4Event* evt) { @@ -89,14 +90,16 @@ void Geant4TestEventAction::end(const G4Event* evt) { /// begin-of-run callback void Geant4TestEventAction::beginRun(const G4Run* run) { - printout(INFO, name(), "%s> calling beginRun(run_id=%d,num_event=%d)", m_type.c_str(), run->GetRunID(), - run->GetNumberOfEventToBeProcessed()); + printout(INFO, name(), "%s> calling beginRun(run_id=%d,num_event=%d)", + m_type.c_str(), run->GetRunID(), + run->GetNumberOfEventToBeProcessed()); } /// End-of-run callback void Geant4TestEventAction::endRun(const G4Run* run) { - printout(INFO, name(), "%s> calling endRun(run_id=%d, num_event=%d)", m_type.c_str(), run->GetRunID(), - run->GetNumberOfEvent()); + printout(INFO, name(), "%s> calling endRun(run_id=%d, num_event=%d)", + m_type.c_str(), run->GetRunID(), + run->GetNumberOfEvent()); } Geant4TestTrackAction::Geant4TestTrackAction(Geant4Context* c, const std::string& n) @@ -108,14 +111,16 @@ Geant4TestTrackAction::~Geant4TestTrackAction() { } /// Begin-of-tracking callback void Geant4TestTrackAction::begin(const G4Track* trk) { - printout(INFO, name(), "%s> calling begin(track=%d, parent=%d, position=(%f,%f,%f))", m_type.c_str(), trk->GetTrackID(), - trk->GetParentID(), trk->GetPosition().x(), trk->GetPosition().y(), trk->GetPosition().z()); + printout(INFO, name(), "%s> calling begin(track=%d, parent=%d, position=(%f,%f,%f))", + m_type.c_str(), trk->GetTrackID(), + trk->GetParentID(), trk->GetPosition().x(), trk->GetPosition().y(), trk->GetPosition().z()); } /// End-of-tracking callback void Geant4TestTrackAction::end(const G4Track* trk) { - printout(INFO, name(), "%s> calling end(track=%d, parent=%d, position=(%f,%f,%f))", m_type.c_str(), trk->GetTrackID(), - trk->GetParentID(), trk->GetPosition().x(), trk->GetPosition().y(), trk->GetPosition().z()); + printout(INFO, name(), "%s> calling end(track=%d, parent=%d, position=(%f,%f,%f))", + m_type.c_str(), trk->GetTrackID(), + trk->GetParentID(), trk->GetPosition().x(), trk->GetPosition().y(), trk->GetPosition().z()); } Geant4TestStepAction::Geant4TestStepAction(Geant4Context* c, const std::string& n) @@ -143,23 +148,27 @@ Geant4TestSensitive::~Geant4TestSensitive() { /// Begin-of-tracking callback void Geant4TestSensitive::begin(G4HCofThisEvent* hce) { Geant4HitCollection* c = collectionByID(m_collectionID); - printout(INFO, name(), "%s> calling begin(num_coll=%d, coll=%s)", m_type.c_str(), hce->GetNumberOfCollections(), - c ? c->GetName().c_str() : "None"); + printout(INFO, name(), "%s> calling begin(num_coll=%d, coll=%s)", + m_type.c_str(), hce->GetNumberOfCollections(), + c ? c->GetName().c_str() : "None"); } /// End-of-tracking callback void Geant4TestSensitive::end(G4HCofThisEvent* hce) { Geant4HitCollection* c = collection(m_collectionID); - printout(INFO, name(), "%s> calling end(num_coll=%d, coll=%s)", m_type.c_str(), hce->GetNumberOfCollections(), - c ? c->GetName().c_str() : "None"); + printout(INFO, name(), "%s> calling end(num_coll=%d, coll=%s)", + m_type.c_str(), hce->GetNumberOfCollections(), + c ? c->GetName().c_str() : "None"); } /// Method for generating hit(s) using the information of G4Step object. bool Geant4TestSensitive::process(G4Step* step, G4TouchableHistory*) { Geant4HitCollection* c = collection(m_collectionID); printout(INFO, name(), "%s> calling process(track=%d, dE=%f, dT=%f len=%f, First,last in Vol=(%c,%c), coll=%s)", - m_type.c_str(), step->GetTrack()->GetTrackID(), step->GetTotalEnergyDeposit(), step->GetDeltaTime(), - step->GetStepLength(), step->IsFirstStepInVolume() ? 'Y' : 'N', step->IsLastStepInVolume() ? 'Y' : 'N', - c ? c->GetName().c_str() : "None"); + m_type.c_str(), step->GetTrack()->GetTrackID(), + step->GetTotalEnergyDeposit(), step->GetDeltaTime(), + step->GetStepLength(), step->IsFirstStepInVolume() ? 'Y' : 'N', + step->IsLastStepInVolume() ? 'Y' : 'N', + c ? c->GetName().c_str() : "None"); return true; } diff --git a/DDG4/src/Geant4TrackPersistency.cpp b/DDG4/src/Geant4TrackPersistency.cpp new file mode 100644 index 0000000000000000000000000000000000000000..91edb34f68780e5481ef876dd7fc2dfd53358536 --- /dev/null +++ b/DDG4/src/Geant4TrackPersistency.cpp @@ -0,0 +1,538 @@ +// $Id: Geant4Field.cpp 888 2013-11-14 15:54:56Z markus.frank@cern.ch $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +// Framework include files +#include "DD4hep/Printout.h" +#include "DD4hep/Primitives.h" +#include "DD4hep/InstanceCount.h" +#include "DDG4/Geant4TrackPersistency.h" +#include "DDG4/Geant4EventAction.h" +#include "DDG4/Geant4TrackingAction.h" +#include "DDG4/Geant4SteppingAction.h" +#include "DDG4/Geant4MonteCarloRecordManager.h" +#include "DDG4/Geant4TrackHandler.h" +#include "DDG4/Geant4StepHandler.h" + +#include "G4TrackingManager.hh" +#include "G4ParticleDefinition.hh" +#include "G4Track.hh" +#include "G4Step.hh" +#include "G4TrackStatus.hh" +#include <set> +#include <stdexcept> + +using namespace std; +using DD4hep::InstanceCount; +using namespace DD4hep; +using namespace DD4hep::Simulation; + + + +#include <map> +namespace DD4hep { + namespace Simulation { + + struct Track; + struct Vertex; + typedef std::vector<int> TrackIDs; + typedef std::vector<Track*> Tracks; + typedef std::vector<Vertex*> Vertices; + typedef Geant4TrackPersistency::FourMomentum FourMomentum; + + struct Vertex { + G4ThreeVector position; + Tracks outgoing; + Track* incoming; + int refCount; + Vertex() : incoming(0), refCount(0) { + InstanceCount::increment(this); + } + virtual ~Vertex() { + InstanceCount::decrement(this); + } + void addRef() { + ++refCount; + } + int release() { + int tmp = --refCount; + if ( tmp < 1 ) { + delete this; + } + return tmp; + } + }; + + typedef std::map<int,void*> TrackMap; + typedef std::map<const void*,Track*> TrackMap2; + + struct Track { + enum State { + KEEP=0x1, + CREATED_HIT=0x2, + VALID_G4TRACK=0x4, + TRACKED_G4TRACK=0x8 + }; + + TrackMap2 secondaries; + TrackIDs contained; + FourMomentum momentum; + Vertices daughters; + Vertex* vertex; + Track* parent; + int id; + int refCount; + int flag; + const G4VProcess* process; + const G4ParticleDefinition* def; + + /// Standard constructor + Track() : vertex(0), parent(0), id(0), refCount(0), flag(0), process(0), def(0) + { + InstanceCount::increment(this); + } + virtual ~Track() + { + InstanceCount::decrement(this); + } + bool toBeKept() const { return (flag&KEEP) == KEEP; } + bool createdHit() const { return (flag&CREATED_HIT) == CREATED_HIT; } + bool validG4Track() const { return (flag&VALID_G4TRACK) == VALID_G4TRACK; } + bool wasTracked() const { return (flag&TRACKED_G4TRACK) == TRACKED_G4TRACK; } + }; + } +} + + +void collectContained(const Track* trk, std::map<int,int>& equiv) { + if ( !trk->contained.empty() ) { + for(TrackIDs::const_iterator i=trk->contained.begin(); i!=trk->contained.end(); ++i) + equiv[*i] = trk->id; + } + for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) + collectContained((*j).second,equiv); +} + +void printTrack(int flg, const char* msg); +void printTree(const Track* trk,const char* msg, int cnt); +void printProperties(const Track* trk); +void printSecondaries(const Track* trk); +void printTrack(const Track* trk, int flg=0, const char* msg="Interesting"); +/// Release track object +void releaseTrack(Geant4TrackPersistency* p, Track* trk); +bool cleanTree(Geant4TrackPersistency* p, const string& proc, Track* trk); + +bool keepTrackTree(const string& proc, Track* trk) { + // We have to check if any of the daughters has the KEEP flag ON. + // If YES, this track may not be deleted, otherwise, we can release the track + string pnam = trk->process ? trk->process->GetProcessName().c_str() : "---"; + if ( pnam == proc ) trk->flag &= ~(Track::KEEP); + if ( (trk->flag&Track::KEEP) == Track::KEEP ) return true; + for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) + if ( keepTrackTree(proc,(*j).second) ) return true; + return false; +} + +bool cleanTree(Geant4TrackPersistency* p, const string& proc, Track* trk) { + bool remove = false; + for(TrackMap2::iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) + remove |= cleanTree(p,proc,(*j).second); + + if ( remove ) { + string pnam = trk->process ? trk->process->GetProcessName().c_str() : "---"; + if ( remove || pnam == proc ) { + bool keep = keepTrackTree(proc,trk); + if ( !keep ) { + releaseTrack(p,trk); + return true; + } + } + } + return false; +} + +void printTree(const Track* trk,const char* msg, int cnt) { + char text[64]; + ::sprintf(text," %-8d | ",cnt); + printf(" %-6d ",cnt); + ::sprintf(text,"| "); + std::string m(text); + m += msg; + printTrack(trk, 1, msg); + int c = 0; + for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) + printTree((*j).second,m.c_str(),++c); +} + +void printProperties(const Track* trk) { + const char* pnam = trk->process ? trk->process->GetProcessName().c_str() : "---"; + int pid = trk->parent ? trk->parent->id : 0; + int sec = int(trk->secondaries.size()); + int vtx = 0; + for(Vertices::const_iterator i=trk->daughters.begin(); i!=trk->daughters.end();++i) + vtx += (*i)->outgoing.size(); + if ( vtx != sec ) { + ::printf(" ++++ ERROR! ++++ "); + } + ::printf("Track %4d Parent:%4d %-8s %-6s Sec:%3d [%3d] Cont:%3d Keep:%s Tracked:%s Hit:%s", + trk->id, pid, pnam, trk->def->GetParticleName().c_str(), sec, vtx, + int(trk->contained.size()), + yes_no(trk->toBeKept()), yes_no(trk->wasTracked()), + yes_no(trk->createdHit())); +} + +void printContained(const Track* trk) { + const TrackIDs& ids = trk->contained; + if ( !ids.empty() ) { + ::printf(" Equiv(%d):",int(ids.size())); + for(TrackIDs::const_iterator i=ids.begin(); i!=ids.end(); ++i) ::printf(" %d",*i); + } +} + +void printSecondaries(const Track* trk) { + char text[256]; + int ndau = 1; + for(TrackMap2::const_iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j, ++ndau) { + ::sprintf(text," ---> Daughter [%3d]",ndau); + ::printTrack((*j).second,2,text); + } +} +void printParents(const Track* trk) { + char text[256]; + int npar = 1; + for(Track* p=trk->parent; p; p=p->parent, ++npar) { + ::sprintf(text," ---> Parent [%3d]",npar); + printTrack(p,4,text); + } +} + +void printTrack(const Track* trk, int flg, const char* msg) { + printf("%s> ",msg); + printProperties(trk); + if ( (flg&1)==1 ) printContained(trk); + printf("\n"); + if ( (flg&2)==2 ) printSecondaries(trk); + if ( (flg&4)==4 ) printParents(trk); +} + +void referenceTrack(Track* trk) { + if ( trk ) { + ++trk->refCount; + referenceTrack(trk->parent); + } +} + +/// Release track object +void releaseTrack(Geant4TrackPersistency* p, Track* trk) { + int cnt = --trk->refCount; + Track* par = trk->parent; + if ( cnt < 1 ) { + int id = trk->id; + Vertex* v = trk->vertex; + if ( v ) { + if ( par ) { + Tracks& trks = v->outgoing; + for(Tracks::iterator it=trks.begin(); it != trks.end(); ++it) { + if ( *it == trk ) { trks.erase(it); break; } + } + if ( trks.size() == 0 ) { + Vertices::iterator iv=find(par->daughters.begin(),par->daughters.end(),v); + if ( iv != par->daughters.end() ) { + par->daughters.erase(iv); + } + } + for(TrackMap2::iterator is=par->secondaries.begin(); is!=par->secondaries.end(); ++is) { + if ( (*is).second == trk ) { + par->secondaries.erase(is); + break; + } + } + par->contained.push_back(id); + } + v->release(); + } + std::for_each(trk->daughters.begin(),trk->daughters.end(),releasePtr<Vertex>); + delete trk; + + TrackMap::iterator m = p->m_trackMap.find(id); + if ( m != p->m_trackMap.end() ) p->m_trackMap.erase(m); + //printf(" ---> DT: REMOVE track %p [%s] %d size:%d\n", (void*)trk, pn, id, int(p->m_trackMap.size())); + } + if ( par ) { + releaseTrack(p, par); + } +} + +/// Initialize track data +Track* initTrack(const G4Track* track, Track* trk, Track* par) { + G4ThreeVector mom = track->GetMomentum(); + trk->momentum.SetPxPyPzE(mom.x(),mom.y(),mom.z(),track->GetTotalEnergy()); + trk->process = track->GetCreatorProcess(); + trk->def = track->GetParticleDefinition(); + trk->id = track->GetTrackID(); + trk->parent = par; + if ( par ) { + par->secondaries[track] = trk; + } + referenceTrack(trk); + return trk; +} + +void deleteTrackTree(Track* trk) { + Vertex* v = trk->vertex; + if ( v ) { + //if ( v->refCount == 2 ) --v->refCount; + v->release(); + } + for(TrackMap2::iterator j=trk->secondaries.begin(); j!=trk->secondaries.end(); ++j) + deleteTrackTree((*j).second); + delete trk; +} + +/// Standard constructor +Geant4TrackPersistency::TrackInfo::TrackInfo() : G4VUserTrackInformation() { + this->manager = 0; + this->track = 0; + this->store = false; + this->createdHit = false; + this->initialEnergy = 0; + InstanceCount::increment(this); +} + +/// Default destructor +Geant4TrackPersistency::TrackInfo::~TrackInfo() { + InstanceCount::decrement(this); +} + +/// Clear alkl stored information for this track +void Geant4TrackPersistency::TrackInfo::set(const G4Track* trk, void* inf) { + this->track = trk; + if ( trk ) { + G4ThreeVector mom = trk->GetMomentum(); + this->initialMomentum.SetPxPyPzE(mom.x(),mom.y(),mom.z(),trk->GetTotalEnergy()); + this->initialEnergy = trk->GetTotalEnergy(); + this->info = inf; + } + else { + this->initialMomentum.SetPxPyPzE(0e0,0e0,0e0,0e0); + this->initialEnergy = 0; + this->createdHit = false; + this->store = false; + this->info = 0; + } +} + +/// Standard constructor +Geant4TrackPersistency::Geant4TrackPersistency(Geant4Context* context, const std::string& nam) + : Geant4Action(context,nam), Geant4MonteCarloTruth(), m_current() +{ + m_current.manager = this; + eventAction().callAtBegin(this,&Geant4TrackPersistency::beginEvent); + eventAction().callAtFinal(this,&Geant4TrackPersistency::endEvent); + trackingAction().callAtFinal(this,&Geant4TrackPersistency::end,CallbackSequence::FRONT); + trackingAction().callUpFront(this,&Geant4TrackPersistency::begin,CallbackSequence::FRONT); + steppingAction().call(this,&Geant4TrackPersistency::step); + + declareProperty("PrintEndTracking",m_printEndTracking = false); + declareProperty("PrintStartTracking",m_printStartTracking = false); + + InstanceCount::increment(this); +} + +/// Default destructor +Geant4TrackPersistency::~Geant4TrackPersistency() { + InstanceCount::decrement(this); +} + +/// Mark a Geant4 track to be kept for later MC truth analysis +void Geant4TrackPersistency::mark(const G4Track* track, bool created_hit) { + mark(track); // Does all the checking... + if ( created_hit ) { + Track* trk = (Track*)m_current.info; + trk->flag |= Track::CREATED_HIT; + m_current.createdHit = created_hit; + } +} + +/// Store a track produced in a step to be kept for later MC truth analysis +void Geant4TrackPersistency::mark(const G4Step* step, bool created_hit) { + if ( step ) { + mark(step->GetTrack(),created_hit); + return; + } + except("Cannot mark the G4Track if the step-pointer is invalid!", c_name()); +} + +/// Mark a Geant4 track of the step to be kept for later MC truth analysis +void Geant4TrackPersistency::mark(const G4Step* step) { + if ( step ) { + mark(step->GetTrack()); + return; + } + except("Cannot mark the G4Track if the step-pointer is invalid!", c_name()); +} + +/// Mark a Geant4 track of the step to be kept for later MC truth analysis +void Geant4TrackPersistency::mark(const G4Track* track) { + if ( 0 == m_current.track ) { + except("Miconfigured program. The tracking preaction was never called for the Geant4TrackPersistency"); + } + else if ( track == m_current.track ) { + m_current.store = true; + Track* trk = (Track*)m_current.info; + trk->flag |= Track::KEEP; + referenceTrack(trk); + return; + } + except("Call to Geant4TrackPersistency with inconsistent G4Track pointer: %p <> %p", + track,m_current.track); +} + +/// Pre-event action callback +void Geant4TrackPersistency::beginEvent(const G4Event* ) { + m_trackMap.clear(); +} + +/// Pre-event action callback +void Geant4TrackPersistency::endEvent(const G4Event* ) { + int cnt = 1; + printf("---------------------------- End of Event: %d tracks kept ----------------------------\n", + int(m_trackMap.size())); + for (TrackMap::iterator m=m_trackMap.begin(); m != m_trackMap.end(); ++m, ++cnt) { + Track* trk = (Track*)(*m).second; + if ( trk->parent == 0 ) { + printTree(trk,"+--",1); + cleanTree(this,"eBrem",trk); + break; + } + } + + printf("---------------------------- End of Event: %d tracks kept ----------------------------\n", + int(m_trackMap.size())); + for(TrackMap::iterator m=m_trackMap.begin(); m != m_trackMap.end(); ++m, ++cnt) { + Track* trk = (Track*)(*m).second; + if ( trk->parent == 0 ) { + printTree(trk,"+--",1); + cleanTree(this,"eIoni",trk); + break; + } + } + + printf("---------------------------- End of Event: %d tracks kept ----------------------------\n", + int(m_trackMap.size())); + for(TrackMap::iterator m=m_trackMap.begin(); m != m_trackMap.end(); ++m, ++cnt) { + Track* trk = (Track*)(*m).second; + if ( trk->parent == 0 ) { + printTree(trk,"+--",1); + deleteTrackTree(trk); + break; + } + + //char text[64]; + //::sprintf(text,"Intersting [%3d]",cnt); + //trk->print(7,text); + } + m_trackMap.clear(); +} + +/// User stepping callback +void Geant4TrackPersistency::step(const G4Step* step, G4SteppingManager* /* mgr */) { + typedef std::vector<const G4Track*> _Sec; + const _Sec* sec=step->GetSecondaryInCurrentStep(); + if ( sec->size() > 0 ) { + Track* t = (Track*)m_current.info; + Vertex* v = new Vertex(); + v->incoming = t; + t->daughters.push_back(v); + + for(_Sec::const_iterator i=sec->begin(); i!=sec->end(); ++i) { + Track* secondary = initTrack(*i,new Track(),t); + secondary->vertex = v; + secondary->vertex->addRef(); + v->outgoing.push_back(secondary); + } + } +} + +/// Pre-track action callback +void Geant4TrackPersistency::begin(const G4Track* track) { + G4VUserTrackInformation* info = &m_current; + int pid = track->GetParentID(); + Track* trk = 0; + + if ( pid == 0 ) { + trk = (Track*)initTrack(track,new Track(),0); + trk->vertex = new Vertex(); + trk->vertex->incoming = 0; + trk->vertex->position = track->GetPosition(); + m_trackMap.insert(make_pair(trk->id,trk)); + } + else { + Track* par = (Track*)m_trackMap[pid]; + trk = par->secondaries[track]; + trk->id = track->GetTrackID(); + trk->vertex->incoming = par; + trk->vertex->position = track->GetPosition(); + m_trackMap.insert(make_pair(trk->id,trk)); + } + + if ( m_printStartTracking ) { + Geant4TrackHandler t(track); + std::string proc = t.creatorName(); + G4ThreeVector pos = t.position(); + + this->info("START Track %p ID=%4d Parent:%4d E:%7.2f MeV " + "%-6s %-8s L:%7.2f Pos:(%8.2f %8.2f %8.2f) Ref:%d",track, + t.id(), t.parent(), t.energy(), + t.name().c_str(), ("["+proc+"]").c_str(), t.length(), + pos.x(),pos.y(),pos.z(),trk->refCount); + } + + trk->flag |= Track::VALID_G4TRACK; + m_current.set(track,trk); + trackMgr()->SetUserTrackInformation(info); + + Geant4TrackHandler t(track); + std::string proc = t.creatorName(); + if ( proc == "conv" ) { + mark(track); + } + if ( proc == "eBrem" ) { + mark(track); + } + else if ( trk->id < 5 ) + mark(track); +} + +/// Post-track action callback +void Geant4TrackPersistency::end(const G4Track* track) { + Track* trk = (Track*)m_current.info; + + trk->flag |= Track::TRACKED_G4TRACK; + trk->flag &= ~(Track::VALID_G4TRACK); + if ( m_printEndTracking ) { + Geant4TrackHandler t(track); + std::string proc = t.creatorName(); + G4ThreeVector pos = t.position(); + + this->info("END Track %p ID=%4d Parent:%4d E:%7.2f MeV " + "%-6s %-8s L:%7.2f Pos:(%8.2f %8.2f %8.2f) Ref:%d Keep:%s Tracked:%s Hit:%s", + track,t.id(), t.parent(), t.energy(), + t.name().c_str(), ("["+proc+"]").c_str(), t.length(), + pos.x(),pos.y(),pos.z(),trk->refCount, + yes_no(trk->toBeKept()), yes_no(trk->wasTracked()), + yes_no(trk->createdHit())); + } + /// If required save Track record... + if ( m_current.store ) { + mcRecordMgr().save(m_current); + } + m_current.set(0,0); + trackMgr()->SetUserTrackInformation(0); + // All done. Check if the track can be released/deleted + releaseTrack(this,trk); +} diff --git a/DDG4/src/Geant4TrackingAction.cpp b/DDG4/src/Geant4TrackingAction.cpp index df74d8657bd7e5bb1396073b38b32d938e61117f..ec703b5a63d653bf40c552b14983fd7707a99053 100644 --- a/DDG4/src/Geant4TrackingAction.cpp +++ b/DDG4/src/Geant4TrackingAction.cpp @@ -10,6 +10,7 @@ // Framework include files #include "DD4hep/InstanceCount.h" #include "DDG4/Geant4TrackingAction.h" +#include "DDG4/Geant4MonteCarloTruth.h" // Geant4 include files #include "G4Track.hh" @@ -36,6 +37,8 @@ Geant4TrackingActionSequence::Geant4TrackingActionSequence(Geant4Context* contex Geant4TrackingActionSequence::~Geant4TrackingActionSequence() { m_actors(&Geant4TrackingAction::release); m_actors.clear(); + m_front.clear(); + m_final.clear(); m_begin.clear(); m_end.clear(); InstanceCount::decrement(this); @@ -53,6 +56,7 @@ void Geant4TrackingActionSequence::adopt(Geant4TrackingAction* action) { /// Pre-track action callback void Geant4TrackingActionSequence::begin(const G4Track* track) { + m_front(track); m_actors(&Geant4TrackingAction::begin, track); m_begin(track); } @@ -61,6 +65,7 @@ void Geant4TrackingActionSequence::begin(const G4Track* track) { void Geant4TrackingActionSequence::end(const G4Track* track) { m_end(track); m_actors(&Geant4TrackingAction::end, track); + m_final(track); } /// Standard constructor @@ -82,6 +87,11 @@ void Geant4TrackingAction::begin(const G4Track*) { void Geant4TrackingAction::end(const G4Track*) { } +/// Mark the track to be kept for MC truth propagation +void Geant4TrackingAction::mark(const G4Track* track) const { + mcTruthMgr().mark(track); +} + /// Get the valid Geant4 tarck information Geant4TrackInformation* Geant4TrackingAction::trackInfo(G4Track* track) const { if (track) { diff --git a/DDG4/src/Geant4UI.cpp b/DDG4/src/Geant4UI.cpp deleted file mode 100644 index 85c2ba1bfb2ee9dae96d931092d354864330c4b5..0000000000000000000000000000000000000000 --- a/DDG4/src/Geant4UI.cpp +++ /dev/null @@ -1,65 +0,0 @@ -// $Id: Geant4Field.cpp 875 2013-11-04 16:15:14Z markus.frank@cern.ch $ -//==================================================================== -// AIDA Detector description implementation for LCD -//-------------------------------------------------------------------- -// -// Author : M.Frank -// -//==================================================================== -// Framework include files -#include "DDG4/Geant4Field.h" - -using namespace DD4hep::Simulation; -using namespace DD4hep::Geometry; -using namespace DD4hep; -using namespace std; - -#if 0 -#ifdef G4VIS_USE_OPENGLX -#include "G4OpenGLImmediateX.hh" -#include "G4OpenGLStoredX.hh" -#endif - -void startX() { - // Initialize visualization - G4VisManager* visManager = 0; - if ( is_visual ) { - visManager = new G4VisExecutive(); -#ifdef G4VIS_USE_OPENGLX - //visManager->RegisterGraphicsSystem (new G4OpenGLImmediateX()); - //visManager->RegisterGraphicsSystem (new G4OpenGLStoredX()); -#endif - visManager->Initialize(); - } -} -#endif - -#if 0 -Geant4UIManager::Geant4UIManager() -{ - declareProperty(""); -} - -Geant4UIManager::operator() { - // Get the pointer to the User Interface manager - G4UImanager* mgr = G4UImanager::GetUIpointer(); - - if ( !gui_type.empty() ) { // interactive mode : define UI session - G4UIExecutive* ui = 0; - const char* args[] = {"cmd"}; - ui = new G4UIExecutive(1,(char**)args); - if ( is_visual && !vis_setup.empty() ) { - mgr->ApplyCommand("/control/execute vis.mac"); - cout << "++++++++++++++++++++++++++++ executed vis.mac" << endl; - } - - if (ui->IsGUI()) { - if ( !gui_setup.empty() ) { - mgr->ApplyCommand("/control/execute "+gui_setup); - cout << "++++++++++++++++++++++++++++ executed gui.mac" << endl; - } - } - ui->SessionStart(); - delete ui; - } -#endif diff --git a/DDG4/src/Geant4UIManager.cpp b/DDG4/src/Geant4UIManager.cpp new file mode 100644 index 0000000000000000000000000000000000000000..89be211a958d9116757fbb16507abf888e0c1d20 --- /dev/null +++ b/DDG4/src/Geant4UIManager.cpp @@ -0,0 +1,103 @@ +// $Id: Geant4Field.cpp 875 2013-11-04 16:15:14Z markus.frank@cern.ch $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== + +#include "DDG4/Geant4UIManager.h" +#include "DDG4/Geant4Kernel.h" +#include "DD4hep/Primitives.h" +#include "DD4hep/Printout.h" +#include "G4VisExecutive.hh" +#include "G4UImanager.hh" +#include "G4UIsession.hh" +#include "G4VisExecutive.hh" +#include "G4UIExecutive.hh" +#include "G4RunManager.hh" + + +using namespace DD4hep::Simulation; +using namespace std; + +/// Initializing constructor +Geant4UIManager::Geant4UIManager(Geant4Context* context, const std::string& name) + : Geant4Action(context,name), m_vis(0), m_ui(0) +{ + declareProperty("SetupUI", m_uiSetup=""); + declareProperty("SetupVIS", m_visSetup=""); + declareProperty("SessionType", m_sessionType="cmd"); + declareProperty("HaveVIS", m_haveVis=false); + declareProperty("HaveUI", m_haveUI=true); +} + +/// Default destructor +Geant4UIManager::~Geant4UIManager() { +} + +/// Start visualization +G4VisManager* Geant4UIManager::startVis() { + // Initialize visualization + G4VisManager* vis = new G4VisExecutive(); + vis->Initialize(); + return vis; +} + +/// Start UI +G4UIExecutive* Geant4UIManager::startUI() { + G4UIExecutive* ui = 0; + const char* args[] = {"DDG4","",""}; + ui = new G4UIExecutive(1,(char**)args,m_sessionType.c_str()); + return ui; +} + +/// Run UI +void Geant4UIManager::operator()(void* ) { + start(); + stop(); +} + +/// Start manager & session +void Geant4UIManager::start() { + // Get the pointer to the User Interface manager + G4UImanager* mgr = G4UImanager::GetUIpointer(); + + // Start visualization + if ( m_haveVis || !m_visSetup.empty() ) { + m_vis = startVis(); + m_haveUI = true; // No graphics without UI! + } + // Start UI instance + if ( m_haveUI || !m_uiSetup.empty() ) { + m_ui = startUI(); + } + // Configure visualization instance + if ( !m_visSetup.empty() ) { + string cmd = "/control/execute "+m_visSetup; + mgr->ApplyCommand(cmd.c_str()); + printout(INFO,"Geant4UIManager","++ Executed visualization setup:%s",m_visSetup.c_str()); + } + // Configure UI instance + if ( !m_uiSetup.empty() ) { + string cmd = "/control/execute "+m_uiSetup; + mgr->ApplyCommand(cmd.c_str()); + printout(INFO,"Geant4UIManager","++ Executed UI setup:%s",m_uiSetup.c_str()); + } + // Start UI session if present + if ( m_ui ) { + m_ui->SessionStart(); + return; + } + // No UI. Simply execute requested number of events + long numEvent = context()->kernel().property("NumEvents").value<long>(); + printout(INFO,"Geant4UIManager","++ Start run with %d events.",numEvent); + context()->kernel().runManager().BeamOn(numEvent); +} + +/// Stop and release resources +void Geant4UIManager::stop() { + deletePtr(m_vis); + deletePtr(m_ui); +} diff --git a/DDG4/src/Geant4UIMessenger.cpp b/DDG4/src/Geant4UIMessenger.cpp index 906ce22db7740a8d2dc153b2f1e043787200aca8..555f9d1be4c2bce46a062f22522d0b0a635d9208 100644 --- a/DDG4/src/Geant4UIMessenger.cpp +++ b/DDG4/src/Geant4UIMessenger.cpp @@ -65,7 +65,8 @@ void Geant4UIMessenger::addCall(const std::string& name, const std::string& desc void Geant4UIMessenger::exportProperties(PropertyManager& mgr) { InstallProperties installer(m_propertyCmd, m_path, this); m_properties = &mgr; - addCall("show", "Show all properties of Geant4 component:" + m_name, Callback(m_properties).make(&PropertyManager::dump)); + addCall("show", "Show all properties of Geant4 component:" + m_name, + Callback(m_properties).make(&PropertyManager::dump)); m_properties->for_each(installer); } @@ -76,7 +77,8 @@ G4String Geant4UIMessenger::GetCurrentValue(G4UIcommand * c) { const string& n = (*i).second; return (*m_properties)[n].str(); } - printout(INFO, "Geant4UIMessenger", "+++ %s> Failed to access property value.", m_name.c_str()); + printout(INFO, "Geant4UIMessenger", + "+++ %s> Failed to access property value.", m_name.c_str()); return ""; } @@ -86,13 +88,16 @@ void Geant4UIMessenger::SetNewValue(G4UIcommand *c, G4String v) { if (m_properties && i != m_propertyCmd.end()) { const string& n = (*i).second; if (!v.empty()) { - printout(INFO, "Geant4UIMessenger", "+++ %s> Setting new property value %s = %s.", m_name.c_str(), n.c_str(), v.c_str()); - (*m_properties)[n].str(v); + Property& p = (*m_properties)[n]; + p.str(v); + printout(INFO, "Geant4UIMessenger", + "+++ %s> Setting property value %s = %s native:%s.", + m_name.c_str(), n.c_str(), v.c_str(), p.str().c_str()); } else { string value = (*m_properties)[n].str(); - printout(INFO, "Geant4UIMessenger", "+++ %s> Unchanged property value %s = %s.", m_name.c_str(), n.c_str(), - value.c_str()); + printout(INFO, "Geant4UIMessenger", "+++ %s> Unchanged property value %s = %s.", + m_name.c_str(), n.c_str(), value.c_str()); } return; } diff --git a/DDG4/src/ToStream.cpp b/DDG4/src/ToStream.cpp new file mode 100644 index 0000000000000000000000000000000000000000..407c1870b59f87fa86c82dd1528f081278d09c13 --- /dev/null +++ b/DDG4/src/ToStream.cpp @@ -0,0 +1,89 @@ +// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== + +// Framework include files +#include "DDG4/ToStream.h" +#include "XML/Evaluator.h" + +// C/C++ include files +#include <stdexcept> + +namespace DD4hep { + XmlTools::Evaluator& g4Evaluator(); +} +namespace { + XmlTools::Evaluator& eval(DD4hep::g4Evaluator()); +} + +using namespace std; +//============================================================================== +namespace DD4hep { namespace Parsers { + template <typename T> T evaluate_string(const string& value) { + throw "Bad undefined call"; + } + + template <> double evaluate_string<double>(const string& value) { + double result = eval.evaluate(value.c_str()); + if (eval.status() != XmlTools::Evaluator::OK) { + cerr << value << ": "; + eval.print_error(); + throw runtime_error("DD4hep::Properties: Severe error during expression evaluation of " + value); + } + return result; + } + template <> float evaluate_string<float>(const string& value) { + double result = eval.evaluate(value.c_str()); + if (eval.status() != XmlTools::Evaluator::OK) { + cerr << value << ": "; + eval.print_error(); + throw runtime_error("DD4hep::Properties: Severe error during expression evaluation of " + value); + } + return (float) result; + } + } +} + +// ============================================================================ +// print XYZ-point +ostream& DD4hep::Utils::toStream(const ROOT::Math::XYZPoint& obj, ostream& s) { + s << "( "; + toStream(obj.X () , s ); + s << " , "; + toStream(obj.Y () , s ); + s << " , "; + toStream(obj.Z () , s ); + s << " )"; + return s; +} +// ============================================================================ +// print XYZ-vector +ostream& DD4hep::Utils::toStream(const ROOT::Math::XYZVector& obj, ostream& s) { + s << "( "; + toStream(obj.X () , s ); + s << " , "; + toStream(obj.Y () , s ); + s << " , "; + toStream(obj.Z () , s ); + s << " )"; + return s; +} +// ============================================================================ +// print LorentzVector +ostream& DD4hep::Utils::toStream(const ROOT::Math::PxPyPzEVector& obj, ostream& s){ + s << "( "; + toStream(obj.Px () , s , 12 ); + s << " , "; + toStream(obj.Py () , s , 12 ); + s << " , "; + toStream(obj.Pz () , s , 13 ); + s << " , "; + toStream(obj.E () , s , 14 ); + s << " )"; + return s; +}