From 70374fc61d1479f142433fe630533c6c40ca211d Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Mon, 26 Jan 2015 17:46:12 +0000 Subject: [PATCH] Add component to allow the setup of the magnetic field tracking in Geant4 from python. See DDG4/examples/CLICSidSimu.py for the usage. --- DDG4/examples/CLICSidSimu.py | 18 +- DDG4/examples/CLICSidSimuMarkus.py | 96 ++++++--- DDG4/include/DDG4/Geant4ActionPhase.h | 38 ++-- DDG4/include/DDG4/Geant4Config.h | 2 + DDG4/include/DDG4/Geant4Kernel.h | 7 + DDG4/plugins/Geant4Factories.cpp | 3 + DDG4/plugins/Geant4FieldSetup.cpp | 98 ---------- DDG4/plugins/Geant4FieldTrackingSetup.cpp | 226 ++++++++++++++++++++++ DDG4/python/DDG4.py | 95 +++++++-- DDG4/python/DDG4Dict.C | 10 + DDG4/src/Geant4ActionPhase.cpp | 37 +++- DDG4/src/Geant4Exec.cpp | 21 +- DDG4/src/Geant4Handle.cpp | 1 + DDG4/src/Geant4Kernel.cpp | 15 ++ 14 files changed, 502 insertions(+), 165 deletions(-) delete mode 100644 DDG4/plugins/Geant4FieldSetup.cpp create mode 100644 DDG4/plugins/Geant4FieldTrackingSetup.cpp diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py index f05e9845d..793a1d10f 100644 --- a/DDG4/examples/CLICSidSimu.py +++ b/DDG4/examples/CLICSidSimu.py @@ -19,7 +19,6 @@ def run(): install_dir = os.environ['DD4hepINSTALL'] example_dir = install_dir+'/examples/DDG4/examples'; kernel.loadGeometry("file:"+install_dir+"/DDDetectors/compact/SiD.xml") - kernel.loadXML("file:"+example_dir+"/DDG4_field.xml") DDG4.importConstants(lcdd) simple = DDG4.Simple(kernel,tracker='Geant4TrackerCombineAction') @@ -27,6 +26,23 @@ def run(): # Configure UI simple.setupCshUI() + # Configure G4 magnetic field tracking + field = geant4.addConfig('Geant4FieldTrackingSetupAction/MagFieldTrackingSetup') + field.stepper = "HelixSimpleRunge" + field.equation = "Mag_UsualEqRhs" + field.eps_min = 5e-05 * mm + field.eps_max = 0.001 * mm + field.min_chord_step = 0.01 * mm + field.delta_chord = 0.25 * mm + field.delta_intersection = 1e-05 * mm + field.delta_one_step = 0.001 * mm + print '+++++> ',field.name,'-> stepper = ',field.stepper + print '+++++> ',field.name,'-> equation = ',field.equation + print '+++++> ',field.name,'-> eps_min = ',field.eps_min + print '+++++> ',field.name,'-> eps_max = ',field.eps_max + print '+++++> ',field.name,'-> delta_one_step = ',field.delta_one_step + + # Configure Run actions run1 = DDG4.RunAction(kernel,'Geant4TestRunAction/RunInit') run1.Property_int = 12345 diff --git a/DDG4/examples/CLICSidSimuMarkus.py b/DDG4/examples/CLICSidSimuMarkus.py index 15c822c49..05ce087d7 100644 --- a/DDG4/examples/CLICSidSimuMarkus.py +++ b/DDG4/examples/CLICSidSimuMarkus.py @@ -12,6 +12,28 @@ from SystemOfUnits import * @author M.Frank @version 1.0 +import os, time, DDG4 +from DDG4 import OutputLevel as Output +from SystemOfUnits import * + +kernel = DDG4.Kernel() +lcdd = kernel.lcdd() +install_dir = os.environ['DD4hepINSTALL'] +example_dir = install_dir+'/examples/DDG4/examples'; +kernel.loadGeometry("file:"+install_dir+"/DDDetectors/compact/SiD_Markus.xml") +kernel.loadXML("file:"+example_dir+"/DDG4_field.xml") +DDG4.importConstants(lcdd,debug=False) +geant4 = DDG4.Geant4(kernel,tracker='Geant4TrackerCombineAction') +geant4.printDetectors() +# Configure UI +geant4.setupCshUI() + +a = DDG4.PhaseAction(kernel,'Geant4FieldPhaseAction/Geant4FieldPhaseAction_1') +kernel.phase('configure').add(a) + +kernel.configure() + + """ def run(): kernel = DDG4.Kernel() @@ -19,12 +41,29 @@ def run(): install_dir = os.environ['DD4hepINSTALL'] example_dir = install_dir+'/examples/DDG4/examples'; kernel.loadGeometry("file:"+install_dir+"/DDDetectors/compact/SiD_Markus.xml") - kernel.loadXML("file:"+example_dir+"/DDG4_field.xml") + ##kernel.loadXML("file:"+example_dir+"/DDG4_field.xml") DDG4.importConstants(lcdd,debug=False) - simple = DDG4.Simple(kernel,tracker='Geant4TrackerCombineAction') - simple.printDetectors() + geant4 = DDG4.Geant4(kernel,tracker='Geant4TrackerCombineAction') + geant4.printDetectors() # Configure UI - simple.setupCshUI() + geant4.setupCshUI() + + field = geant4.addConfig('Geant4FieldTrackingSetupAction/MagFieldTrackingSetup') + field.stepper = "HelixSimpleRunge" + field.equation = "Mag_UsualEqRhs" + field.eps_min = 5e-05*mm + field.eps_max = 0.001*mm + field.min_chord_step = 0.01*mm + field.delta_chord = 0.25*mm + field.delta_intersection = 1e-05*mm + field.delta_one_step = 0.001*mm + print '+++++> ',field.name,'-> stepper = ',field.stepper + print '+++++> ',field.name,'-> equation = ',field.equation + print '+++++> ',field.name,'-> eps_min = ',field.eps_min + print '+++++> ',field.name,'-> eps_max = ',field.eps_max + print '+++++> ',field.name,'-> delta_one_step = ',field.delta_one_step + + # Configure Run actions run1 = DDG4.RunAction(kernel,'Geant4TestRunAction/RunInit') @@ -41,15 +80,15 @@ def run(): generator_output_level = Output.WARNING # Configure I/O - evt_lcio = simple.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) + evt_lcio = geant4.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) ##evt_lcio.OutputLevel = generator_output_level - evt_root = simple.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) + evt_root = geant4.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV """ Generation of primary particles from LCIO input files """ - + """ # First particle file reader gen = DDG4.GeneratorAction(kernel,"LCIOInputAction/LCIO1"); #gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/e2e2nn_gen_1343_1.stdhep" @@ -66,13 +105,14 @@ def run(): #gen.Input = "LCIOFileReader|/afs/cern.ch/user/n/nikiforo/public/Markus/geantinos.slcio" gen.MomentumScale = 1.0 gen.Mask = 1 - simple.buildInputStage([gen],output_level=generator_output_level) - #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + geant4.buildInputStage([gen],output_level=generator_output_level) """ - gen = simple.setupGun("Gun",particle='geantino',energy=20*GeV,position=(0*mm,0*mm,0*cm),multiplicity=3) + #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + gen = geant4.setupGun("Gun",particle='geantino',energy=20*GeV,position=(0*mm,0*mm,0*cm),multiplicity=3) gen.isotrop = False gen.direction = (1,0,0) gen.OutputLevel = generator_output_level + """ # And handle the simulation particles. part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler") @@ -88,7 +128,7 @@ def run(): user.enableUI() part.adopt(user) """ - """ + """ rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader") rdr.zSpread = 0.0 @@ -99,24 +139,26 @@ def run(): kernel.generatorAction().adopt(rdr) """ - #seq,act = simple.setupTracker('SiTrackerBarrel') + """ + #seq,act = geant4.setupTracker('SiTrackerBarrel') # First the tracking detectors - seq,act = simple.setupTracker('SiVertexBarrel') - seq,act = simple.setupTracker('SiVertexEndcap') - seq,act = simple.setupTracker('SiTrackerBarrel') - seq,act = simple.setupTracker('SiTrackerEndcap') - seq,act = simple.setupTracker('SiTrackerForward') + seq,act = geant4.setupTracker('SiVertexBarrel') + seq,act = geant4.setupTracker('SiVertexEndcap') + seq,act = geant4.setupTracker('SiTrackerBarrel') + seq,act = geant4.setupTracker('SiTrackerEndcap') + seq,act = geant4.setupTracker('SiTrackerForward') # Now the calorimeters - seq,act = simple.setupCalorimeter('EcalBarrel') - seq,act = simple.setupCalorimeter('EcalEndcap') - seq,act = simple.setupCalorimeter('HcalBarrel') - seq,act = simple.setupCalorimeter('HcalEndcap') - seq,act = simple.setupCalorimeter('HcalPlug') - seq,act = simple.setupCalorimeter('MuonBarrel') - seq,act = simple.setupCalorimeter('MuonEndcap') - seq,act = simple.setupCalorimeter('LumiCal') - seq,act = simple.setupCalorimeter('BeamCal') + seq,act = geant4.setupCalorimeter('EcalBarrel') + seq,act = geant4.setupCalorimeter('EcalEndcap') + seq,act = geant4.setupCalorimeter('HcalBarrel') + seq,act = geant4.setupCalorimeter('HcalEndcap') + seq,act = geant4.setupCalorimeter('HcalPlug') + seq,act = geant4.setupCalorimeter('MuonBarrel') + seq,act = geant4.setupCalorimeter('MuonEndcap') + seq,act = geant4.setupCalorimeter('LumiCal') + seq,act = geant4.setupCalorimeter('BeamCal') + """ """ scan = DDG4.SteppingAction(kernel,'Geant4MaterialScanner/MaterialScan') kernel.steppingAction().adopt(scan) @@ -124,7 +166,7 @@ def run(): # Now build the physics list: - phys = simple.setupPhysics('QGSP_BERT') + phys = geant4.setupPhysics('QGSP_BERT') ph = DDG4.PhysicsList(kernel,'Geant4PhysicsList/Myphysics') ph.addParticleConstructor('G4Geantino') ph.addParticleConstructor('G4BosonConstructor') diff --git a/DDG4/include/DDG4/Geant4ActionPhase.h b/DDG4/include/DDG4/Geant4ActionPhase.h index 959b47b9f..3ba94568e 100644 --- a/DDG4/include/DDG4/Geant4ActionPhase.h +++ b/DDG4/include/DDG4/Geant4ActionPhase.h @@ -21,8 +21,26 @@ namespace DD4hep { /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit namespace Simulation { - /* + /// Generic action for Geant4 phases + /** + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4PhaseAction : public Geant4Action { + public: + /// Standard constructor + Geant4PhaseAction(Geant4Context* context, const std::string& name); + /// Default destructor + virtual ~Geant4PhaseAction(); + /// Callback to generate primary particles + virtual void operator()(); + /// Create bound callback to operator()() + virtual Callback callback(); + }; + /* Geant4Phase, G4EventGenerator --> G4VUserPrimaryGeneratorAction::GeneratePrimaries Geant4Begin, G4Run --> G4UserRunAction::BeginOfRunAction Geant4End, G4Run --> G4UserRunAction::EndOfRunAction @@ -44,7 +62,7 @@ namespace DD4hep { */ class Geant4ActionPhase : public Geant4Action { public: - typedef std::vector<Callback> Members; + typedef std::vector<std::pair<Geant4Action*, Callback> > Members; protected: /// Phase members (actions) being called for a particular phase Members m_members; @@ -68,16 +86,15 @@ namespace DD4hep { /// Execute all members in the phase context void execute(void* argument); /// Add a new member to the phase - virtual bool add(Callback callback); + virtual bool add(Geant4Action* action, Callback callback); /// Remove an existing member from the phase. If not existing returns false - virtual bool remove(Callback callback); + virtual bool remove(Geant4Action* action, Callback callback); /// Add a new member to the phase template <typename TYPE, typename IF_TYPE, typename A0, typename R> bool add(TYPE* member, R (IF_TYPE::*callback)(A0 arg)) { typeinfoCheck(typeid(A0), *m_argTypes[0], "Invalid ARG0 type. Failed to add phase callback."); if (dynamic_cast<IF_TYPE*>(member)) { - //member->addRef(); - return add(Callback(member).make(callback)); + return add(member,Callback(member).make(callback)); } throw unrelated_type_error(typeid(TYPE), typeid(IF_TYPE), "Failed to add phase callback."); } @@ -87,8 +104,7 @@ namespace DD4hep { typeinfoCheck(typeid(A0), *m_argTypes[0], "Invalid ARG0 type. Failed to add phase callback."); typeinfoCheck(typeid(A1), *m_argTypes[1], "Invalid ARG1 type. Failed to add phase callback."); if (dynamic_cast<IF_TYPE*>(member)) { - //member->addRef(); - return add(Callback(member).make(callback)); + return add(member,Callback(member).make(callback)); } throw unrelated_type_error(typeid(TYPE), typeid(IF_TYPE), "Failed to add phase callback."); } @@ -100,18 +116,18 @@ namespace DD4hep { typeinfoCheck(typeid(A2), *m_argTypes[2], "Invalid ARG2 type. Failed to add phase callback."); if (dynamic_cast<IF_TYPE*>(member)) { //member->addRef(); - return add(Callback(member).make(callback)); + return add(member,Callback(member).make(callback)); } throw unrelated_type_error(typeid(TYPE), typeid(IF_TYPE), "Failed to add phase callback."); } /// Remove all member callbacks from the phase. If not existing returns false template <typename TYPE, typename PMF> bool remove(TYPE* member) { - return remove(Callback(member)); + return remove(member,Callback(member)); } /// Remove an existing member callback from the phase. If not existing returns false template <typename TYPE, typename PMF> bool remove(TYPE* member, PMF callback) { Callback cb(member); - return remove(cb.make(callback)); + return remove(member,cb.make(callback)); } /// Create action to execute phase members void call(); diff --git a/DDG4/include/DDG4/Geant4Config.h b/DDG4/include/DDG4/Geant4Config.h index d91288dd7..2690d9ec8 100644 --- a/DDG4/include/DDG4/Geant4Config.h +++ b/DDG4/include/DDG4/Geant4Config.h @@ -23,6 +23,7 @@ namespace DD4hep { class Geant4Kernel; class Geant4Action; class Geant4Filter; + class Geant4PhaseAction; class Geant4RunAction; class Geant4EventAction; class Geant4TrackingAction; @@ -50,6 +51,7 @@ namespace DD4hep { // Actions typedef Geant4Handle<Geant4Action> Action; typedef Geant4Handle<Geant4Filter> Filter; + typedef Geant4Handle<Geant4PhaseAction> PhaseAction; typedef Geant4Handle<Geant4GeneratorAction> GenAction; typedef Geant4Handle<Geant4RunAction> RunAction; typedef Geant4Handle<Geant4EventAction> EventAction; diff --git a/DDG4/include/DDG4/Geant4Kernel.h b/DDG4/include/DDG4/Geant4Kernel.h index 4cc10552a..74f36eaf9 100644 --- a/DDG4/include/DDG4/Geant4Kernel.h +++ b/DDG4/include/DDG4/Geant4Kernel.h @@ -212,6 +212,10 @@ namespace DD4hep { /// Access phase by name Geant4ActionPhase* getPhase(const std::string& name); + + /// Add a new phase to the phase + virtual Geant4ActionPhase* addSimplePhase(const std::string& name, bool throw_on_exist); + /// Add a new phase to the phase virtual Geant4ActionPhase* addPhase(const std::string& name, const std::type_info& arg1, const std::type_info& arg2, const std::type_info& arg3, bool throw_on_exist); @@ -234,6 +238,9 @@ namespace DD4hep { /// Destroy all phases. To be called only at shutdown. virtual void destroyPhases(); + /// Execute phase action if it exists + virtual bool executePhase(const std::string& name, const void** args) const; + /// Access generator action sequence Geant4GeneratorActionSequence* generatorAction(bool create); /// Access generator action sequence diff --git a/DDG4/plugins/Geant4Factories.cpp b/DDG4/plugins/Geant4Factories.cpp index fa91f307e..1aaee887c 100644 --- a/DDG4/plugins/Geant4Factories.cpp +++ b/DDG4/plugins/Geant4Factories.cpp @@ -9,6 +9,9 @@ #include "DDG4/Factories.h" using namespace DD4hep::Simulation; +#include "DDG4/Geant4ActionPhase.h" +DECLARE_GEANT4ACTION(Geant4PhaseAction) + #include "DDG4/Geant4RunAction.h" DECLARE_GEANT4ACTION(Geant4RunActionSequence) diff --git a/DDG4/plugins/Geant4FieldSetup.cpp b/DDG4/plugins/Geant4FieldSetup.cpp deleted file mode 100644 index fd275aacc..000000000 --- a/DDG4/plugins/Geant4FieldSetup.cpp +++ /dev/null @@ -1,98 +0,0 @@ -// $Id: Geant4Setup.cpp 578 2013-05-17 22:33:09Z markus.frank $ -//==================================================================== -// AIDA Detector description implementation for LCD -//-------------------------------------------------------------------- -// -// Author : M.Frank -// -//==================================================================== -// Framework include files -#include "DD4hep/Handle.h" -#include "DD4hep/Fields.h" -#include "DDG4/Factories.h" -#include "DDG4/Geant4Field.h" -#include "DDG4/Geant4Converter.h" - -#include "G4TransportationManager.hh" -#include "G4MagIntegratorStepper.hh" -#include "G4Mag_EqRhs.hh" -#include "G4ChordFinder.hh" -#include "G4PropagatorInField.hh" - -using namespace std; -using namespace DD4hep; -using namespace DD4hep::Simulation; -typedef DD4hep::Geometry::LCDD lcdd_t; - -/// Local declaration in anonymous namespace -namespace { - - struct Geant4FieldSetup; - struct Geant4SetupPropertyMap { - const map<string,string>& vals; - Geant4SetupPropertyMap(const map<string,string>& v) : vals(v) {} - string value(const string& key) const; - double toDouble(const string& key) const; - bool operator[](const string& key) const { return vals.find(key) != vals.end(); } - }; - - string Geant4SetupPropertyMap::value(const string& key) const { - lcdd_t::PropertyValues::const_iterator iV = vals.find(key); - return iV == vals.end() ? "" : (*iV).second; - } - - double Geant4SetupPropertyMap::toDouble(const string& key) const { - return Geometry::_toDouble(this->value(key)); - } -} - -static long setup_fields(lcdd_t& lcdd, const DD4hep::Simulation::Geant4Converter& /* cnv */, const map<string,string>& vals) { - Geant4SetupPropertyMap pm(vals); - DD4hep::Geometry::OverlayedField fld = lcdd.field(); - G4MagIntegratorStepper* stepper = 0; - G4MagneticField* field = 0; - G4FieldManager* fieldMgr = 0; - G4TransportationManager* tr = 0; - G4PropagatorInField* propagator = 0; - G4ChordFinder* chordFinder = 0; - - lcdd_t::PropertyValues::const_iterator iV; - - string eq_typ = pm.value("equation"); - string stepper_typ = pm.value("stepper"); - double value; - - field = new Simulation::Geant4Field(fld); - tr = G4TransportationManager::GetTransportationManager(); - fieldMgr = tr->GetFieldManager(); - - fieldMgr->SetFieldChangesEnergy(fld.changesEnergy()); - fieldMgr->SetDetectorField(field); - - G4Mag_EqRhs* equation = ROOT::Reflex::PluginService::Create<G4Mag_EqRhs*>(eq_typ,field); - stepper = ROOT::Reflex::PluginService::Create<G4MagIntegratorStepper*>(stepper_typ,equation); - - iV = vals.find("min_chord_step"); - value = Geometry::_toDouble((iV == vals.end()) ? string("1.0e-2 * mm") : (*iV).second); - chordFinder = new G4ChordFinder(field,value,stepper); - propagator = tr->GetPropagatorInField(); - fieldMgr->SetChordFinder(chordFinder); - - if ( pm["delta_chord"] ) { - chordFinder->SetDeltaChord(pm.toDouble("delta_chord")); - } - if ( pm["delta_one_step"] ) { - fieldMgr->SetAccuraciesWithDeltaOneStep(pm.toDouble("delta_one_step")); - } - if ( pm["delta_intersection"] ) { - fieldMgr->SetDeltaIntersection(pm.toDouble("delta_intersection")); - } - if ( pm["eps_min"] ) { - propagator->SetMinimumEpsilonStep(pm.toDouble("eps_min")); - } - if ( pm["eps_max"] ) { - propagator->SetMaximumEpsilonStep(pm.toDouble("eps_max")); - } - return 1; -} -DECLARE_GEANT4_SETUP(Geant4FieldSetup,setup_fields) diff --git a/DDG4/plugins/Geant4FieldTrackingSetup.cpp b/DDG4/plugins/Geant4FieldTrackingSetup.cpp new file mode 100644 index 000000000..a373565cb --- /dev/null +++ b/DDG4/plugins/Geant4FieldTrackingSetup.cpp @@ -0,0 +1,226 @@ +// $Id: Geant4Setup.cpp 578 2013-05-17 22:33:09Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +#ifndef DD4HEP_DDG4_GEANT4FIELDTRACKINGSETUP_H +#define DD4HEP_DDG4_GEANT4FIELDTRACKINGSETUP_H 1 + +// Framework include files +#include "DD4hep/LCDD.h" +#include "DDG4/Geant4ActionPhase.h" + + +/// Namespace for the AIDA detector description toolkit +namespace DD4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace Simulation { + + /// Generic Setup component to perform the magnetic field tracking in Geant4 + /** Geant4FieldTrackingSetup. + * + * This base class is use jointly between the XML setup and the + * phase action used by the python setup. + * + * Note: + * Negative parameters are not passed to Geant4 objects, but ignored -- if possible. + * + * @author M.Frank + * @version 1.0 + */ + struct Geant4FieldTrackingSetup { + protected: + /** Variables to be filled before calling execute */ + /// Name of the G4Mag_EqRhs class + std::string eq_typ; + /// Name of the G4MagIntegratorStepper class + std::string stepper_typ; + /// G4ChordFinder parameter: min_chord_step + double min_chord_step; + /// G4ChordFinder parameter: delta + double delta_chord; + /// G4FieldManager parameter: delta_one_step + double delta_one_step; + /// G4FieldManager parameter: delta_intersection + double delta_intersection; + /// G4PropagatorInField parameter: eps_min + double eps_min; + /// G4PropagatorInField parameter: eps_min + double eps_max; + + public: + /// Default constructor + Geant4FieldTrackingSetup(); + /// Default destructor + virtual ~Geant4FieldTrackingSetup(); + /// Perform the setup of the magnetic field tracking in Geant4 + virtual int execute(Geometry::LCDD& lcdd); + }; + + /// Phase action to perform the setup of the Geant4 tracking in magnetic fields + /** Geant4FieldTrackingSetupAction. + * + * The phase action configures the Geant4FieldTrackingSetup base class using properties + * and then configures the Geant4 tracking in magnetic fields. + * + * @author M.Frank + * @version 1.0 + */ + class Geant4FieldTrackingSetupAction : public Geant4PhaseAction, public Geant4FieldTrackingSetup { + protected: + public: + /// Standard constructor + Geant4FieldTrackingSetupAction(Geant4Context* context, const std::string& nam); + + /// Default destructor + virtual ~Geant4FieldTrackingSetupAction() {} + + /// Phase action callback + void operator()(); + + }; + } // End namespace Simulation +} // End namespace DD4hep +#endif // DD4HEP_DDG4_GEANT4FIELDTRACKINGSETUP_H + + +// $Id: Geant4Setup.cpp 578 2013-05-17 22:33:09Z markus.frank $ +//==================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------- +// +// Author : M.Frank +// +//==================================================================== +// Framework include files +#include "DD4hep/Handle.h" +#include "DD4hep/Fields.h" +#include "DDG4/Factories.h" +#include "DDG4/Geant4Field.h" +#include "DDG4/Geant4Converter.h" + +#include "G4TransportationManager.hh" +#include "G4MagIntegratorStepper.hh" +#include "G4Mag_EqRhs.hh" +#include "G4ChordFinder.hh" +#include "G4PropagatorInField.hh" +#include <limits> + +using namespace std; +using namespace DD4hep; +using namespace DD4hep::Simulation; +typedef DD4hep::Geometry::LCDD lcdd_t; + +/// Local declaration in anonymous namespace +namespace { + + struct Geant4SetupPropertyMap { + const map<string,string>& vals; + Geant4SetupPropertyMap(const map<string,string>& v) : vals(v) {} + string value(const string& key) const; + double toDouble(const string& key) const; + bool operator[](const string& key) const { return vals.find(key) != vals.end(); } + }; + + string Geant4SetupPropertyMap::value(const string& key) const { + lcdd_t::PropertyValues::const_iterator iV = vals.find(key); + return iV == vals.end() ? "" : (*iV).second; + } + + double Geant4SetupPropertyMap::toDouble(const string& key) const { + return Geometry::_toDouble(this->value(key)); + } + +} + +/// Default constructor +Geant4FieldTrackingSetup::Geant4FieldTrackingSetup() : eq_typ(), stepper_typ() { + eps_min = -1.0; + eps_max = -1.0; + min_chord_step = 1.0e-2 *mm; + delta_chord = -1.0; + delta_one_step = -1.0; + delta_intersection = -1.0; +} + +/// Default destructor +Geant4FieldTrackingSetup::~Geant4FieldTrackingSetup() { +} + +/// Perform the setup of the magnetic field tracking in Geant4 +int Geant4FieldTrackingSetup::execute(Geometry::LCDD& lcdd) { + Geometry::OverlayedField fld = lcdd.field(); + G4MagneticField* mag_field = new Simulation::Geant4Field(fld); + G4Mag_EqRhs* mag_equation = ROOT::Reflex::PluginService::Create<G4Mag_EqRhs*>(eq_typ,mag_field); + G4MagIntegratorStepper* fld_stepper = ROOT::Reflex::PluginService::Create<G4MagIntegratorStepper*>(stepper_typ,mag_equation); + G4ChordFinder* chordFinder = new G4ChordFinder(mag_field,min_chord_step,fld_stepper); + G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager(); + G4PropagatorInField* propagator = transportMgr->GetPropagatorInField(); + G4FieldManager* fieldManager = transportMgr->GetFieldManager(); + + fieldManager->SetFieldChangesEnergy(fld.changesEnergy()); + fieldManager->SetDetectorField(mag_field); + fieldManager->SetChordFinder(chordFinder); + + if ( delta_chord >= 0e0 ) + chordFinder->SetDeltaChord(delta_chord); + if ( delta_one_step >= 0e0 ) + fieldManager->SetAccuraciesWithDeltaOneStep(delta_one_step); + if ( delta_intersection >= 0e0 ) + fieldManager->SetDeltaIntersection(delta_intersection); + if ( eps_min >= 0e0 ) + propagator->SetMinimumEpsilonStep(eps_min); + if ( eps_max >= 0e0 ) + propagator->SetMaximumEpsilonStep(eps_max); + return 1; +} + +static long setup_fields(lcdd_t& lcdd, const DD4hep::Simulation::Geant4Converter& /* cnv */, const map<string,string>& vals) { + struct XMLFieldTrackingSetup : public Geant4FieldTrackingSetup { + XMLFieldTrackingSetup(const map<string,string>& vals) : Geant4FieldTrackingSetup() { + Geant4SetupPropertyMap pm(vals); + lcdd_t::PropertyValues::const_iterator iV = vals.find("min_chord_step"); + eq_typ = pm.value("equation"); + stepper_typ = pm.value("stepper"); + min_chord_step = Geometry::_toDouble((iV==vals.end()) ? string("1e-2 * mm") : (*iV).second); + if ( pm["eps_min"] ) eps_min = pm.toDouble("eps_min"); + if ( pm["eps_max"] ) eps_max = pm.toDouble("eps_max"); + if ( pm["delta_chord"] ) delta_chord = pm.toDouble("delta_chord"); + if ( pm["delta_one_step"] ) delta_one_step = pm.toDouble("delta_one_step"); + if ( pm["delta_intersection"] ) delta_intersection = pm.toDouble("delta_intersection"); + } + virtual ~XMLFieldTrackingSetup() {} + } setup(vals); + return setup.execute(lcdd); +} + +/// Standard constructor +Geant4FieldTrackingSetupAction::Geant4FieldTrackingSetupAction(Geant4Context* context, const std::string& nam) +: Geant4PhaseAction(context,nam), Geant4FieldTrackingSetup() +{ + declareProperty("equation", eq_typ); + declareProperty("stepper", stepper_typ); + declareProperty("min_chord_step", min_chord_step = 1.0e-2); + declareProperty("delta_chord", delta_chord = -1.0); + declareProperty("delta_one_step", delta_one_step = -1.0); + declareProperty("delta_intersection", delta_intersection = -1.0); + declareProperty("eps_min", eps_min = -1.0); + declareProperty("eps_max", eps_max = -1.0); +} + +/// Post-track action callback +void Geant4FieldTrackingSetupAction::operator()() { + execute(context()->lcdd()); + print("Geant4 magnetic field tracking configured. G4MagIntegratorStepper:%s G4Mag_EqRhs:%s " + "Epsilon:[min:%f mm max:%f mm] Delta:[chord:%f 1-step:%f intersect:%f]", + stepper_typ.c_str(),eq_typ.c_str(),eps_min, eps_max, + delta_chord,delta_one_step,delta_intersection); +} + +DECLARE_GEANT4_SETUP(Geant4FieldSetup,setup_fields) +DECLARE_GEANT4ACTION(Geant4FieldTrackingSetupAction) + diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py index c2233f3ee..dbc9890ab 100644 --- a/DDG4/python/DDG4.py +++ b/DDG4/python/DDG4.py @@ -115,11 +115,16 @@ def _setKernelProperty(self, name, value): msg = 'Geant4Kernel::SetProperty [Unhandled]: Cannot set Kernel.'+name+' = '+str(value) raise exceptions.KeyError(msg) +#--------------------------------------------------------------------------- +def _kernel_phase(self,name): return self.addSimplePhase(name,False) +#--------------------------------------------------------------------------- +Kernel.phase = _kernel_phase 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) @@ -128,6 +133,8 @@ def Action(kernel,nam): return Interface.createAction(kernel,nam) #--------------------------------------------------------------------------- def Filter(kernel,nam): return Interface.createFilter(kernel,nam) #--------------------------------------------------------------------------- +def PhaseAction(kernel,nam): return Interface.createPhaseAction(kernel,nam) +#--------------------------------------------------------------------------- def RunAction(kernel,nam): return Interface.createRunAction(kernel,nam) #--------------------------------------------------------------------------- def EventAction(kernel,nam): return Interface.createEventAction(kernel,nam) @@ -152,7 +159,15 @@ def _setup(obj): setattr(o,'adopt',_adopt) setattr(o,'add',_adopt) +def _setup_callback(obj): + def _adopt(self,action): self.__adopt(action.get(),action.callback()) + _import_class('Sim',obj) + o = getattr(current,obj) + setattr(o,'__adopt',getattr(o,'add')) + setattr(o,'add',_adopt) + #--------------------------------------------------------------------------- +_setup_callback('Geant4ActionPhase') _setup('Geant4RunActionSequence') _setup('Geant4EventActionSequence') _setup('Geant4GeneratorActionSequence') @@ -165,6 +180,7 @@ _setup('Geant4Sensitive') _setup('Geant4ParticleHandler') _import_class('Sim','Geant4Filter') _import_class('Sim','Geant4RunAction') +_import_class('Sim','Geant4PhaseAction') _import_class('Sim','Geant4UserParticleHandler') #--------------------------------------------------------------------------- @@ -202,6 +218,7 @@ def _props(obj): _props('FilterHandle') _props('ActionHandle') +_props('PhaseActionHandle') _props('RunActionHandle') _props('EventActionHandle') _props('GeneratorActionHandle') @@ -224,29 +241,12 @@ _props('SensDetActionSequenceHandle') _props('Geant4PhysicsListActionSequence') - - -#class Iter(): -# def Iterator(self): -# ''' Fix for map iteration on macos ''' -# n = self.m.size() -# it = self.m.begin() -# for i in range(0,n): -# yield it -# it.__preinc__() -# def __init__(self,m): -# self.m = m -# def __iter__(self): -# return self.Iterator() - - - """ Helper object to perform stuff, which occurs very often. I am sick of typing the same over and over again. -@author M.Frank -@version 1.0 +\author M.Frank +\version 1.0 """ class Simple: @@ -261,6 +261,61 @@ class Simple: self.sensitive_types['tracker'] = tracker self.sensitive_types['calorimeter'] = calo + def addPhaseAction(self,phase_name,factory_specification,ui=True): + action = PhaseAction(self.kernel,factory_specification) + self.kernel.phase('configure').add(action) + if ui: action.enableUI() + return action + + """ + Add a new phase action to the 'configure' step. + Called at the beginning of Geant4Exec::configure. + The factory specification is the typical string "<factory_name>/<instance name>". + If no instance name is specified it defaults to the factory name. + + \author M.Frank + """ + def addConfig(self, factory_specification): + return self.addPhaseAction('configure',factory_specification) + + """ + Add a new phase action to the 'initialize' step. + Called at the beginning of Geant4Exec::initialize. + The factory specification is the typical string "<factory_name>/<instance name>". + If no instance name is specified it defaults to the factory name. + + \author M.Frank + """ + def addInit(self, factory_specification): + return self.addPhaseAction('initialize',factory_specification) + + """ + Add a new phase action to the 'start' step. + Called at the beginning of Geant4Exec::run. + The factory specification is the typical string "<factory_name>/<instance name>". + If no instance name is specified it defaults to the factory name. + + \author M.Frank + """ + def addStart(self, factory_specification): + return self.addPhaseAction('start',factory_specification) + + """ + Add a new phase action to the 'stop' step. + Called at the end of Geant4Exec::run. + The factory specification is the typical string "<factory_name>/<instance name>". + If no instance name is specified it defaults to the factory name. + + \author M.Frank + """ + def addStop(self, factory_specification): + return self.addPhaseAction('stop',factory_specification) + + """ + Execute the geant 4 program with all steps. + + \author M.Frank + """ def execute(self): self.kernel.configure() self.kernel.initialize() @@ -269,7 +324,7 @@ class Simple: return self def printDetectors(self): - print '+++ List of detectors:' + print '+++ List of sensitive detectors:' for i in self.lcdd.detectors(): o = DetElement(i.second) sd = self.lcdd.sensitiveDetector(o.name()) diff --git a/DDG4/python/DDG4Dict.C b/DDG4/python/DDG4Dict.C index 22c665d51..baba4aed4 100644 --- a/DDG4/python/DDG4Dict.C +++ b/DDG4/python/DDG4Dict.C @@ -59,6 +59,7 @@ namespace DD4hep { ACTIONHANDLE(Filter); ACTIONHANDLE(Action); + ACTIONHANDLE(PhaseAction); ACTIONHANDLE(RunAction); ACTIONHANDLE(EventAction); ACTIONHANDLE(GeneratorAction); @@ -97,6 +98,8 @@ namespace DD4hep { { return cr<ActionHandle,Setup::Action>(kernel,name_type); } static FilterHandle createFilter(KernelHandle& kernel, const string& name_type) { return cr<FilterHandle,Setup::Filter>(kernel,name_type); } + static PhaseActionHandle createPhaseAction(KernelHandle& kernel, const string& name_type) + { return cr<PhaseActionHandle,Setup::PhaseAction>(kernel,name_type); } static PhysicsListHandle createPhysicsList(KernelHandle& kernel, const string& name_type) { return cr<PhysicsListHandle,Setup::PhysicsList>(kernel,name_type); } static RunActionHandle createRunAction(KernelHandle& kernel, const string& name_type) @@ -118,6 +121,7 @@ namespace DD4hep { static Geant4Action* toAction(Geant4Filter* f) { return f; } static Geant4Action* toAction(Geant4Action* f) { return f; } + static Geant4Action* toAction(Geant4PhaseAction* f) { return f; } static Geant4Action* toAction(Geant4Sensitive* f) { return f; } static Geant4Action* toAction(Geant4PhysicsList* f) { return f; } static Geant4Action* toAction(Geant4RunAction* f) { return f; } @@ -137,6 +141,7 @@ namespace DD4hep { static Geant4Action* toAction(FilterHandle f) { return f.action; } static Geant4Action* toAction(ActionHandle f) { return f.action; } + static Geant4Action* toAction(PhaseActionHandle f) { return f.action; } static Geant4Action* toAction(SensitiveHandle f) { return f.action; } static Geant4Action* toAction(PhysicsListHandle f) { return f.action; } static Geant4Action* toAction(RunActionHandle f) { return f.action; } @@ -195,6 +200,7 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation; #pragma link C++ class ActionHandle; #pragma link C++ class FilterHandle; #pragma link C++ class RunActionHandle; +#pragma link C++ class PhaseActionHandle; #pragma link C++ class GeneratorActionHandle; #pragma link C++ class EventActionHandle; #pragma link C++ class PhysicsListHandle; @@ -249,7 +255,11 @@ typedef DD4hep::Simulation::Geant4ActionCreation Geant4ActionCreation; #pragma link C++ class Geant4ActionSD; #pragma link C++ class Geant4Sensitive; #pragma link C++ class Geant4SensDetActionSequence; + #pragma link C++ class Geant4ActionPhase; +#pragma link C++ class Geant4PhaseAction; +#pragma link C++ class Callback; +#pragma link C++ class Callback::mfunc_t; // Work around CINT bug: // somehow the symbol Geometry moved into global namespace. Redeclare it here diff --git a/DDG4/src/Geant4ActionPhase.cpp b/DDG4/src/Geant4ActionPhase.cpp index 9efcf941c..418f6322a 100644 --- a/DDG4/src/Geant4ActionPhase.cpp +++ b/DDG4/src/Geant4ActionPhase.cpp @@ -14,6 +14,24 @@ using namespace std; using namespace DD4hep::Simulation; +/// Standard constructor +Geant4PhaseAction::Geant4PhaseAction(Geant4Context* context, const std::string& name) +: Geant4Action(context,name) +{ +} + +/// Default destructor +Geant4PhaseAction::~Geant4PhaseAction() { +} + +/// Callback to generate primary particles +void Geant4PhaseAction::operator()() { +} + +DD4hep::Callback Geant4PhaseAction::callback() { + return Callback(this).make(&Geant4PhaseAction::operator()); +} + /// Standard constructor Geant4ActionPhase::Geant4ActionPhase(Geant4Context* context, const string& nam, const type_info& arg_type0, const type_info& arg_type1, const type_info& arg_type2) @@ -26,21 +44,25 @@ Geant4ActionPhase::Geant4ActionPhase(Geant4Context* context, const string& nam, /// Default destructor Geant4ActionPhase::~Geant4ActionPhase() { + for (Members::iterator i = m_members.begin(); i != m_members.end(); ++i) + (*i).first->release(); m_members.clear(); InstanceCount::decrement(this); } /// Add a new member to the phase -bool Geant4ActionPhase::add(Callback callback) { - m_members.push_back(callback); +bool Geant4ActionPhase::add(Geant4Action* action, Callback callback) { + action->addRef(); + m_members.push_back(make_pair(action,callback)); return true; } /// Remove an existing member from the phase. If not existing returns false -bool Geant4ActionPhase::remove(Callback callback) { - if (callback.func.first) { - Members::iterator i = find(m_members.begin(), m_members.end(), callback); +bool Geant4ActionPhase::remove(Geant4Action* action, Callback callback) { + if (action && callback.func.first) { + Members::iterator i = find(m_members.begin(), m_members.end(), make_pair(action,callback)); if (i != m_members.end()) { + (*i).first->release(); m_members.erase(i); return true; } @@ -48,7 +70,8 @@ bool Geant4ActionPhase::remove(Callback callback) { } size_t len = m_members.size(); for (Members::iterator i = m_members.begin(); i != m_members.end(); ++i) { - if ((*i).par == callback.par) { + if ( (*i).first == action && (*i).second.par == callback.par) { + (*i).first->release(); m_members.erase(i); i = m_members.begin(); } @@ -59,7 +82,7 @@ bool Geant4ActionPhase::remove(Callback callback) { /// Execute all members in the phase context void Geant4ActionPhase::execute(void* argument) { for (Members::iterator i = m_members.begin(); i != m_members.end(); ++i) { - (*i).execute((const void**) &argument); + (*i).second.execute((const void**) &argument); } } diff --git a/DDG4/src/Geant4Exec.cpp b/DDG4/src/Geant4Exec.cpp index 040c8386f..217bd3f18 100644 --- a/DDG4/src/Geant4Exec.cpp +++ b/DDG4/src/Geant4Exec.cpp @@ -75,6 +75,13 @@ namespace DD4hep { releasePtr(m_sequence); InstanceCount::decrement(this); } + Geant4Context* context() const { + return m_activeContext; + } + Geant4Kernel& kernel() const { + return context()->kernel(); + } + void setContextToClients() { (Geant4Action::ContextUpdate(m_activeContext))(m_sequence); } @@ -263,6 +270,7 @@ namespace DD4hep { /// Begin-of-run callback void Geant4UserRunAction::BeginOfRunAction(const G4Run* run) { createClientContext(run); + kernel().executePhase("begin-run",(const void**)&run); eventAction->setContextToClients(); m_sequence->begin(run); } @@ -270,6 +278,7 @@ namespace DD4hep { /// End-of-run callback void Geant4UserRunAction::EndOfRunAction(const G4Run* run) { m_sequence->end(run); + kernel().executePhase("end-run",(const void**)&run); eventAction->releaseContextFromClients(); destroyClientContext(run); } @@ -278,6 +287,7 @@ namespace DD4hep { void Geant4UserEventAction::BeginOfEventAction(const G4Event* evt) { runAction->setContextToClients(); setContextToClients(); + kernel().executePhase("begin-event",(const void**)&evt); m_sequence->begin(evt); } @@ -285,6 +295,7 @@ namespace DD4hep { void Geant4UserEventAction::EndOfEventAction(const G4Event* evt) { m_sequence->end(evt); runAction->releaseContextFromClients(); + kernel().executePhase("end-event",(const void**)&evt); destroyClientContext(evt); } } @@ -311,6 +322,8 @@ int Geant4Exec::configure(Geant4Kernel& kernel) { // For now do this: /* Geant4Random* rnd = */ s_globalRandom = new Geant4Random(); + kernel.executePhase("configure",0); + // Construct the default run manager G4RunManager& runManager = kernel.runManager(); @@ -379,6 +392,7 @@ int Geant4Exec::initialize(Geant4Kernel& kernel) { // // Initialize G4 engine // + kernel.executePhase("initialize",0); runManager.Initialize(); return 1; } @@ -387,12 +401,15 @@ int Geant4Exec::initialize(Geant4Kernel& kernel) { int Geant4Exec::run(Geant4Kernel& kernel) { Property& p = kernel.property("UI"); string value = p.value<string>(); + + kernel.executePhase("start",0); if ( !value.empty() ) { Geant4Action* ui = kernel.globalAction(value); if ( ui ) { Geant4Call* c = dynamic_cast<Geant4Call*>(ui); if ( c ) { (*c)(0); + kernel.executePhase("stop",0); return 1; } ui->except("++ Geant4Exec: Failed to start UI interface."); @@ -401,11 +418,13 @@ int Geant4Exec::run(Geant4Kernel& kernel) { } long nevt = kernel.property("NumEvents").value<long>(); kernel.runManager().BeamOn(nevt); + kernel.executePhase("stop",0); return 1; } /// Run the simulation -int Geant4Exec::terminate(Geant4Kernel&) { +int Geant4Exec::terminate(Geant4Kernel& kernel) { + kernel.executePhase("terminate",0); deletePtr(s_globalContext); return 1; } diff --git a/DDG4/src/Geant4Handle.cpp b/DDG4/src/Geant4Handle.cpp index f69bcc684..64bd95f8a 100644 --- a/DDG4/src/Geant4Handle.cpp +++ b/DDG4/src/Geant4Handle.cpp @@ -223,6 +223,7 @@ template class Geant4Handle<Geant4Action> ; template class Geant4Handle<Geant4Filter> ; template class Geant4Handle<Geant4Sensitive> ; template class Geant4Handle<Geant4ActionPhase> ; +template class Geant4Handle<Geant4PhaseAction> ; template class Geant4Handle<Geant4GeneratorAction> ; template class Geant4Handle<Geant4RunAction> ; template class Geant4Handle<Geant4EventAction> ; diff --git a/DDG4/src/Geant4Kernel.cpp b/DDG4/src/Geant4Kernel.cpp index 4d652ca34..08369aacb 100644 --- a/DDG4/src/Geant4Kernel.cpp +++ b/DDG4/src/Geant4Kernel.cpp @@ -303,6 +303,16 @@ Geant4Action* Geant4Kernel::globalFilter(const std::string& filter_name, bool th return (*i).second; } +/// Execute phase action if it exists +bool Geant4Kernel::executePhase(const std::string& nam, const void** arguments) const { + Phases::const_iterator i = m_phases.find(nam); + if (i != m_phases.end()) { + (*i).second->execute(arguments); + return true; + } + return false; +} + /// Access phase by name Geant4ActionPhase* Geant4Kernel::getPhase(const std::string& nam) { Phases::const_iterator i = m_phases.find(nam); @@ -313,6 +323,11 @@ Geant4ActionPhase* Geant4Kernel::getPhase(const std::string& nam) { "does not exist. [No-Entry]", nam.c_str())); } +/// Add a new phase to the phase +Geant4ActionPhase* Geant4Kernel::addSimplePhase(const std::string& name, bool throw_on_exist) { + return addPhase(name,typeid(void),typeid(void),typeid(void),throw_on_exist); +} + /// Add a new phase Geant4ActionPhase* Geant4Kernel::addPhase(const std::string& nam, const type_info& arg0, const type_info& arg1, const type_info& arg2, bool throw_on_exist) { -- GitLab