From 8e1257f5bb9c2682370ae26b9350094e1ef23c45 Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Thu, 29 Oct 2015 21:43:38 +0000 Subject: [PATCH] DDG4 improvements. See doc/release.notes --- DDDetectors/compact/SiD_Markus.xml | 123 ++---- DDDetectors/src/CaloFaceBarrel_surfaces.cpp | 46 +-- DDDetectors/src/CaloFaceEndcap_surfaces.cpp | 44 +-- DDG4/CMakeLists.txt | 5 +- DDG4/examples/CLICSidSimuMarkus.py | 24 +- DDG4/examples/SiD_Markus.py | 185 +++++++++ DDG4/include/DDG4/Geant4Action.h | 15 +- DDG4/include/DDG4/Geant4Data.h | 15 + DDG4/include/DDG4/Geant4IsotropeGenerator.h | 9 + DDG4/include/DDG4/Geant4PhysicsList.h | 3 + DDG4/include/DDG4/Geant4SensDetAction.h | 4 +- DDG4/include/DDG4/Geant4StepHandler.h | 27 ++ DDG4/plugins/Geant4PhysicsLists.cpp | 169 --------- DDG4/{src => plugins}/Geant4ROOTDump.cpp | 2 +- DDG4/plugins/Geant4TrackerWeightedSD.cpp | 393 ++++++++++++++++++++ DDG4/reco/Geant4SurfaceTest.cpp | 197 ++++++++++ DDG4/src/Geant4Action.cpp | 2 +- DDG4/src/Geant4Data.cpp | 2 +- DDG4/src/Geant4Exec.cpp | 15 +- DDG4/src/Geant4GeneratorActionInit.cpp | 2 +- DDG4/src/Geant4IsotropeGenerator.cpp | 8 +- DDG4/src/Geant4PhysicsList.cpp | 20 +- DDG4/src/Geant4SensDetAction.cpp | 4 +- DDG4/src/Geant4StepHandler.cpp | 32 +- DDG4/src/Geant4TrackingAction.cpp | 3 +- DDRec/include/DDRec/Surface.h | 12 +- cmake/DD4hepBuild.cmake | 21 +- doc/release.notes | 12 + 28 files changed, 1024 insertions(+), 370 deletions(-) create mode 100644 DDG4/examples/SiD_Markus.py delete mode 100644 DDG4/plugins/Geant4PhysicsLists.cpp rename DDG4/{src => plugins}/Geant4ROOTDump.cpp (99%) create mode 100644 DDG4/plugins/Geant4TrackerWeightedSD.cpp create mode 100644 DDG4/reco/Geant4SurfaceTest.cpp diff --git a/DDDetectors/compact/SiD_Markus.xml b/DDDetectors/compact/SiD_Markus.xml index 5aabb185d..b7185809d 100644 --- a/DDDetectors/compact/SiD_Markus.xml +++ b/DDDetectors/compact/SiD_Markus.xml @@ -12,13 +12,11 @@ </info> <includes> - <gdmlFile ref="./elements.xml"/> - <gdmlFile ref="./materials.xml"/> + <gdmlFile ref="elements.xml"/> + <gdmlFile ref="materials.xml"/> </includes> <define> - <constant name="materials_dir" value="." type="string"/>; - <constant name="SiD_dir" value="./SiD" type="string"/>; <constant name="world_side" value="30000*mm"/> <constant name="world_x" value="world_side"/> <constant name="world_y" value="world_side"/> @@ -210,101 +208,11 @@ </display> <comment>Additional design specific material definitions</comment> - <include ref="./SiD/SiD_Materials.xml"/> -<!-- - <display> - <vis name="SiTrackerLayerVis" alpha="1.0" r="1.0" g="1.0" b="0.6" showDaughters="true" visible="true"/> - <vis name="SiTrackerModuleVis" alpha="0.1" r="0.0" g="1.0" b="0.6" showDaughters="false" visible="true"/> - <vis name="SiTrackerForwardVis" alpha="1.0" r="0.8" g="0.1" b="0.1" showDaughters="false" visible="true"/> - </display> + <include ref="SiD/SiD_Materials.xml"/> - <readouts> - <readout name="SiTrackerBarrelHits"> - <id>system:8,barrel:3,layer:4,module:14,sensor:2,side:32:-2,strip:20</id> - </readout> - </readouts> ---> - <!-- Includes for sensitives and support --> - <!-- include ref="SiD/SiD_TrackerBarrel.xml" --> -<!-- ---> <comment>Tracking detectors</comment> - <include ref="./SiD/SiD_Tracker.xml"/> - <include ref="./SiD/SiD_Ecal.xml"/> - - <plugins> - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"> - <argument value="SiTrackerBarrel"/> - <argument value="dimension=1"/> - </plugin> - <plugin name="DD4hep_SurfaceExamplePlugin"> - <argument value="SiTrackerBarrel"/> - <argument value="aa=1*mm"/> - <argument value="bb=2*cm"/> - <argument value="cc=3*m"/> - </plugin> - </plugins> -<!-- - - <plugins> - <plugin name="TestSurfaces"><argument value="EcalEndcap"/></plugin> - <plugin name="DD4hep_PolyhedraEndcapCalorimeterSurfacePlugin"><argument value="EcalEndcap"/></plugin> - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"><argument value="EcalBarrel"/></plugin> - </plugins> - - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"> - <argument value="SiTrackerBarrel"/> - </plugin> - <plugin name="TestSurfaces"><argument value="SiTrackerBarrel"/></plugin> - - <plugin name="TestSurfaces"> <argument value="EcalEndcap"/></plugin> - <plugin name="DD4hep_SiTrackerEndcapSurfacePlugin"><argument value="EcalEndcap"/></plugin> - - <plugin name="SubdetectorExtensionPlugin"> <argument value="EcalEndcap"/></plugin> - <plugin name="LayeringExtensionPlugin"> <argument value="EcalEndcap"/></plugin> - - <include ref="SiD/SiD_Ecal.xml"/> - <plugin name="TestSurfaces"> <argument value="EcalBarrel"/></plugin> - <plugin name="SubdetectorExtensionPlugin"> <argument value="EcalBarrel"/></plugin> - <plugin name="LayeringExtensionPlugin"> <argument value="EcalBarrel"/></plugin> - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"><argument value="EcalBarrel"/></plugin> - - <include ref="SiD/SiD_HcalBarrel.xml"/> - <plugin name="TestSurfaces"> <argument value="HcalBarrel"/></plugin> - <plugin name="SubdetectorExtensionPlugin"> <argument value="HcalBarrel"/></plugin> - <plugin name="LayeringExtensionPlugin"> <argument value="HcalBarrel"/></plugin> - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"><argument value="HcalBarrel"/></plugin> - - <include ref="SiD/SiD_Muon.xml"/> - <plugins> - <plugin name="SubdetectorExtensionPlugin"> <argument value="MuonBarrel"/></plugin> - <plugin name="LayeringExtensionPlugin"> <argument value="MuonBarrel"/></plugin> - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"><argument value="MuonBarrel"/></plugin> - </plugins> ---> -<!-- - <plugins> - <plugin name="TestSurfaces"> - <argument value="SiTrackerEndcap"/> - </plugin> - </plugins> - - <include ref="SiD/SiD_Tracker.xml"/> - <plugins> - <plugin name="DD4hep_SiTrackerEndcapSurfacePlugin"> - <argument value="SiTrackerEndcap"/> - </plugin> - </plugins> - <plugins> - <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"> - <argument value="SiTrackerBarrel"/> - </plugin> - </plugins> - <include ref="SiD/SiD_TrackerEndcap.xml"/> - <include ref="SiD/SiD_TrackerBarrel.xml"/> ---> -<!-- <include ref="SiD/SiD_Vertex.xml"/> +<!-- <include ref="SiD/SiD_Tracker.xml"/> <comment>Calorimeters</comment> @@ -315,13 +223,14 @@ <include ref="SiD/SiD_Lumical.xml"/> <include ref="SiD/SiD_Beamcal.xml"/> - <comment>Dead material and supports</comment> - <include ref="SiD/SiD_Solenoid.xml"/> - <include ref="SiD/SiD_Shields.xml"/> - <comment>Beampipe</comment> <include ref="SiD/SiD_Beampipe.xml"/> + <comment>Dead material, supports and magnet</comment> + <include ref="SiD/SiD_Shields.xml"/> + <include ref="SiD/SiD_Solenoid.xml"/> +--> + <fields> <field name="GlobalSolenoid" type="solenoid" inner_field="5.0*tesla" @@ -330,6 +239,18 @@ outer_radius="SolenoidalFieldRadius"> </field> </fields> ---> + <plugins> + <plugin name="DD4hep_SiTrackerEndcapSurfacePlugin"> + <argument value="SiVertexEndcap"/> + <argument value="dimension=1"/> + </plugin> + </plugins> + <plugins> + <plugin name="DD4hep_SiTrackerBarrelSurfacePlugin"> + <argument value="SiVertexBarrel"/> + <argument value="dimension=1"/> + </plugin> + <plugin name="InstallSurfaceManager"/> + </plugins> </lccdd> diff --git a/DDDetectors/src/CaloFaceBarrel_surfaces.cpp b/DDDetectors/src/CaloFaceBarrel_surfaces.cpp index 2e9b7c2ce..252bd34c5 100644 --- a/DDDetectors/src/CaloFaceBarrel_surfaces.cpp +++ b/DDDetectors/src/CaloFaceBarrel_surfaces.cpp @@ -37,16 +37,17 @@ namespace { namespace{ /// helper class for a planar surface placed into a polyhedral calorimeter barrel - class CaloBarrelPlaneImpl : public DD4hep::DDRec::VolPlaneImpl{ + class CaloBarrelPlaneImpl : public DD4hep::DDRec::VolPlaneImpl { double _length, _width ; + public: /// standard c'tor with all necessary arguments - origin is (0,0,0) if not given. CaloBarrelPlaneImpl( DDSurfaces::SurfaceType typ, - double thickness_inner ,double thickness_outer, - DDSurfaces::Vector3D u_val ,DDSurfaces::Vector3D v_val , - DDSurfaces::Vector3D n_val , DDSurfaces::Vector3D o_val, - DD4hep::Geometry::Volume vol, int id ) : - DD4hep::DDRec::VolPlaneImpl( typ, thickness_inner,thickness_outer, u_val, v_val, n_val, o_val, vol, id), + double thickness_inner ,double thickness_outer, + DDSurfaces::Vector3D u_val ,DDSurfaces::Vector3D v_val , + DDSurfaces::Vector3D n_val , DDSurfaces::Vector3D o_val, + DD4hep::Geometry::Volume vol, int id_val ) : + DD4hep::DDRec::VolPlaneImpl( typ, thickness_inner,thickness_outer, u_val, v_val, n_val, o_val, vol, id_val), _length(0),_width(0) {} void setData( double length, double width){ @@ -54,15 +55,14 @@ namespace{ _width = width ; } - void setID( DD4hep::long64 id ) { _id = id ; } + void setID( DD4hep::long64 id_val ) { _id = id_val ; } // overwrite to include points inside the inner radius of the barrel bool insideBounds(const DDSurfaces::Vector3D& point, double epsilon) const { - DDSurfaces::Vector2D uvVec = globalToLocal( point ) ; return ( std::abs ( distance( point ) ) < epsilon ) && - std::abs( uvVec[0] ) < _width/2. && std::abs( uvVec[1] ) < _length/2. ; + std::abs( uvVec[0] ) < _width/2. && std::abs( uvVec[1] ) < _length/2. ; } /// create outer bounding lines for the given symmetry of the polyhedron @@ -93,15 +93,15 @@ namespace{ std::cout << "DD4hep_CaloFaceBarrelSurfacePlugin: argument[" << i << "] = " << name << " = " << value << std::endl; - if( name=="length" ) data.length = value ; - else if( name=="radius" ) data.radius = value ; - else if( name=="phi0" ) data.phi0 = value ; - else if( name=="symmetry") data.symmetry = value ; - else if( name=="systemID") data.systemID = value ; - else if( name=="encoding") data.encoding = value ; - else { - std::cout << "DD4hep_CaloFaceBarrelSurfacePlugin: WARNING unknown parameter: " << name << std::endl ; - } + if( name=="length" ) data.length = value ; + else if( name=="radius" ) data.radius = value ; + else if( name=="phi0" ) data.phi0 = value ; + else if( name=="symmetry") data.symmetry = value ; + else if( name=="systemID") data.systemID = value ; + else if( name=="encoding") data.encoding = value ; + else { + std::cout << "DD4hep_CaloFaceBarrelSurfacePlugin: WARNING unknown parameter: " << name << std::endl ; + } } } } @@ -126,7 +126,7 @@ namespace{ double outer_thickness = 1e-6 ; std::cout << "DD4hep_CaloFaceBarrelSurfacePlugin: install tracking surfaces for : " - << component.name() << std::endl ; + << component.name() << std::endl ; DD4hep::DDSegmentation::BitField64 bf( "system:5,side:-2,layer:9,module:8,sensor:8" ) ; @@ -144,10 +144,10 @@ namespace{ double gam = phi0 + alpha/2. + i*alpha; Vector3D - u( cos(gam+M_PI/2.), sin(gam+M_PI/2.), 0. ), - v( 0. , 0. , 1. ), - n( cos(gam) , sin(gam) , 0. ), - o( radius*cos(gam) , radius*sin(gam) , 0. ); + u( cos(gam+M_PI/2.), sin(gam+M_PI/2.), 0. ), + v( 0. , 0. , 1. ), + n( cos(gam) , sin(gam) , 0. ), + o( radius*cos(gam) , radius*sin(gam) , 0. ); CaloBarrelPlane surf(comp_vol,Type(Type::Helper,Type::Sensitive), inner_thickness, outer_thickness, u, v, n, o); diff --git a/DDDetectors/src/CaloFaceEndcap_surfaces.cpp b/DDDetectors/src/CaloFaceEndcap_surfaces.cpp index d4d40dda0..554e555e8 100644 --- a/DDDetectors/src/CaloFaceEndcap_surfaces.cpp +++ b/DDDetectors/src/CaloFaceEndcap_surfaces.cpp @@ -43,11 +43,11 @@ namespace{ public: /// standard c'tor with all necessary arguments - origin is (0,0,0) if not given. CaloEndcapPlaneImpl( DDSurfaces::SurfaceType typ, - double thickness_inner ,double thickness_outer, - DDSurfaces::Vector3D u_val ,DDSurfaces::Vector3D v_val , - DDSurfaces::Vector3D n_val , DDSurfaces::Vector3D o_val, - DD4hep::Geometry::Volume vol, int id ) : - DD4hep::DDRec::VolPlaneImpl( typ, thickness_inner,thickness_outer, u_val, v_val, n_val, o_val, vol, id), + double thickness_inner ,double thickness_outer, + DDSurfaces::Vector3D u_val ,DDSurfaces::Vector3D v_val , + DDSurfaces::Vector3D n_val , DDSurfaces::Vector3D o_val, + DD4hep::Geometry::Volume vol, int id_val ) : + DD4hep::DDRec::VolPlaneImpl( typ, thickness_inner,thickness_outer, u_val, v_val, n_val, o_val, vol, id_val), _r(0),_phi0(0),_sym(0) {} void setData( double radius, double phi0, unsigned symmetry){ @@ -55,7 +55,7 @@ namespace{ _phi0 = phi0 ; _sym = symmetry ; } - void setID( DD4hep::long64 id ) { _id = id ; } + void setID( DD4hep::long64 id_val ) { _id = id_val; } // overwrite to include points inside the inner radius of the endcap bool insideBounds(const DDSurfaces::Vector3D& point, double epsilon) const { @@ -63,8 +63,8 @@ namespace{ double ri = _r * cos( 2.* M_PI / _sym ) ; return ( std::abs ( distance( point ) ) < epsilon ) && - ( (point.rho() < ri ) || - ( volume()->GetShape()->Contains( const_cast<double*> (point.const_array() ) ) ) ) ; + ( (point.rho() < ri ) || + ( volume()->GetShape()->Contains( const_cast<double*> (point.const_array() ) ) ) ) ; } /// create outer bounding lines for the given symmetry of the polyhedron @@ -75,10 +75,10 @@ namespace{ double alpha = ( _sym ? 2.* M_PI / _sym : 0. ) ; for(unsigned i=0 ; i < _sym ; ++i){ - double gam0 = i * alpha + _phi0 ; - double gam1 = (i+1) * alpha + _phi0 ; - lines.push_back( std::make_pair( DDSurfaces::Vector3D( _r*cos(gam0), _r*sin(gam0), origin().z() ), - DDSurfaces::Vector3D( _r*cos(gam1), _r*sin(gam1), origin().z() ) ) ) ; + double gam0 = i * alpha + _phi0 ; + double gam1 = (i+1) * alpha + _phi0 ; + lines.push_back( std::make_pair( DDSurfaces::Vector3D( _r*cos(gam0), _r*sin(gam0), origin().z() ), + DDSurfaces::Vector3D( _r*cos(gam1), _r*sin(gam1), origin().z() ) ) ) ; } return lines; } @@ -98,15 +98,15 @@ namespace{ std::cout << "DD4hep_CaloFaceEndcapSurfacePlugin: argument[" << i << "] = " << name << " = " << value << std::endl; - if( name=="zpos" ) data.zpos = value ; - else if( name=="radius" ) data.radius = value ; - else if( name=="phi0" ) data.phi0 = value ; - else if( name=="symmetry") data.symmetry = value ; - else if( name=="systemID") data.systemID = value ; - else if( name=="encoding") data.encoding = value ; - else { - std::cout << "DD4hep_CaloFaceEndcapSurfacePlugin: WARNING unknown parameter: " << name << std::endl ; - } + if( name=="zpos" ) data.zpos = value ; + else if( name=="radius" ) data.radius = value ; + else if( name=="phi0" ) data.phi0 = value ; + else if( name=="symmetry") data.symmetry = value ; + else if( name=="systemID") data.systemID = value ; + else if( name=="encoding") data.encoding = value ; + else { + std::cout << "DD4hep_CaloFaceEndcapSurfacePlugin: WARNING unknown parameter: " << name << std::endl ; + } } } } @@ -130,7 +130,7 @@ namespace{ double outer_thickness = 1e-6 ; std::cout << "DD4hep_CaloFaceEndcapSurfacePlugin: install tracking surfaces for : " - << component.name() << std::endl ; + << component.name() << std::endl ; DD4hep::DDSegmentation::BitField64 bf( "system:5,side:-2,layer:9,module:8,sensor:8" ) ; diff --git a/DDG4/CMakeLists.txt b/DDG4/CMakeLists.txt index 6a1e2f92f..0383c15b3 100644 --- a/DDG4/CMakeLists.txt +++ b/DDG4/CMakeLists.txt @@ -28,12 +28,15 @@ dd4hep_add_plugin(DDG4Plugins #--------------------------- LCIO Plugins for new simulation framework ----------- dd4hep_add_plugin(DDG4LCIO OPTIONAL [LCIO REQUIRED SOURCES lcio/*.cpp] ) +#--------------------------- DDRec dependent Plugins ----------------------------- +dd4hep_add_plugin(DDG4Reco + OPTIONAL [DDRec REQUIRED SOURCES reco/*.cpp] ) #----------------------------------------------------------------------------------- dd4hep_add_executable(g4gdmlDisplay SOURCES g4gdmlDisplay.cpp) #----------------------------------------------------------------------------------- dd4hep_add_executable(g4FromXML SOURCES g4FromXML.cpp) #----------------------------------------------------------------------------------- dd4hep_add_executable(dd_sim SOURCES ddsim.cpp) -#---Package installation procedure(s) ------------------------------------- +#---Package installation procedure(s) ---------------------------------------------- dd4hep_install_dir(examples DESTINATION examples/DDG4) dd4hep_install_files(FILES python/*.py python/*.C DESTINATION python) diff --git a/DDG4/examples/CLICSidSimuMarkus.py b/DDG4/examples/CLICSidSimuMarkus.py index e948d9500..98f720047 100644 --- a/DDG4/examples/CLICSidSimuMarkus.py +++ b/DDG4/examples/CLICSidSimuMarkus.py @@ -46,7 +46,8 @@ def run(): geant4 = DDG4.Geant4(kernel,tracker='Geant4TrackerCombineAction') geant4.printDetectors() # Configure UI - geant4.setupCshUI(macro='run.mac',ui=None) + #geant4.setupCshUI(macro='run.mac',ui=None) + geant4.setupCshUI() field = geant4.addConfig('Geant4FieldTrackingSetupAction/MagFieldTrackingSetup') field.stepper = "HelixSimpleRunge" @@ -82,8 +83,9 @@ def run(): # Configure I/O evt_lcio = geant4.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) ##evt_lcio.OutputLevel = generator_output_level - evt_root = geant4.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) + #evt_root = geant4.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) + prim = DDG4.GeneratorAction(kernel,"Geant4GeneratorActionInit/GenerationInit") #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV """ Generation of primary particles from LCIO input files @@ -108,11 +110,11 @@ def run(): geant4.buildInputStage([gen],output_level=generator_output_level) """ #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - gen = geant4.setupGun("Gun",particle='geantino',energy=20*GeV,position=(0*mm,0*mm,0*cm),multiplicity=3) - gen.isotrop = False + gen = geant4.setupGun("Gun",particle='mu+',energy=20*GeV,position=(0*mm,0*mm,0*cm),multiplicity=3) + gen.isotrop = True gen.direction = (1,0,0) gen.OutputLevel = generator_output_level - + gen.Standalone = False """ # And handle the simulation particles. part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler") @@ -128,6 +130,10 @@ def run(): user.enableUI() part.adopt(user) """ + geant4.buildInputStage([prim,gen],Output.ERROR) + + """ + """ """ rdr = DDG4.GeneratorAction(kernel,"LcioGeneratorAction/Reader") @@ -139,16 +145,16 @@ def run(): kernel.generatorAction().adopt(rdr) """ - """ - #seq,act = geant4.setupTracker('SiTrackerBarrel') # First the tracking detectors - 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 = geant4.setupTracker('SiVertexBarrel') + seq,act = geant4.setupTracker('SiVertexEndcap') + seq,act = geant4.setupCalorimeter('EcalBarrel') seq,act = geant4.setupCalorimeter('EcalEndcap') seq,act = geant4.setupCalorimeter('HcalBarrel') diff --git a/DDG4/examples/SiD_Markus.py b/DDG4/examples/SiD_Markus.py new file mode 100644 index 000000000..a96b8f4b1 --- /dev/null +++ b/DDG4/examples/SiD_Markus.py @@ -0,0 +1,185 @@ +# +# +import os, time, DDG4 +from DDG4 import OutputLevel as Output +from SystemOfUnits import * +# +# +""" + + DD4hep simulation example setup using the python configuration + + @author M.Frank + @version 1.0 + +""" +def run(): + kernel = DDG4.Kernel() + lcdd = kernel.lcdd() + install_dir = os.environ['DD4hepINSTALL'] + kernel.loadGeometry("file:"+install_dir+"/DDDetectors/compact/SiD_Markus.xml") + DDG4.importConstants(lcdd) + DDG4.Core.setPrintLevel(Output.WARNING) + + geant4 = DDG4.Geant4(kernel,tracker='Geant4TrackerWeightedAction') + geant4.printDetectors() + # Configure UI + geant4.setupCshUI() + + # Configure G4 magnetic field tracking + field = geant4.addConfig('Geant4FieldTrackingSetupAction/MagFieldTrackingSetup') + field.stepper = "HelixGeant4Runge" + 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 + + # Setup random generator + rndm = DDG4.Action(kernel,'Geant4Random/Random') + rndm.Seed = 987654321 + rndm.initialize() + + # Configure Run actions + run1 = DDG4.RunAction(kernel,'Geant4TestRunAction/RunInit') + run1.Property_int = 12345 + run1.Property_double = -5e15*keV + run1.Property_string = 'Startrun: Hello_2' + print run1.Property_string, run1.Property_double, run1.Property_int + run1.enableUI() + kernel.registerGlobalAction(run1) + kernel.runAction().adopt(run1) + + # Configure Event actions + prt = DDG4.EventAction(kernel,'Geant4ParticlePrint/ParticlePrint') + prt.OutputLevel = Output.DEBUG + prt.OutputType = 3 # Print both: table and tree + kernel.eventAction().adopt(prt) + + # Configure Event actions + prt = DDG4.EventAction(kernel,'Geant4SurfaceTest/SurfaceTest') + prt.OutputLevel = Output.INFO + kernel.eventAction().adopt(prt) + + # Configure I/O + evt_lcio = geant4.setupLCIOOutput('LcioOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) + evt_lcio.OutputLevel = Output.DEBUG + + #evt_root = geant4.setupROOTOutput('RootOutput','CLICSiD_'+time.strftime('%Y-%m-%d_%H-%M')) + #generator_output_level = Output.INFO + + gen = DDG4.GeneratorAction(kernel,"Geant4GeneratorActionInit/GenerationInit") + kernel.generatorAction().adopt(gen) + + #VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV + """ + Generation of isotrope tracks of a given multiplicity with overlay: + """ + # First particle generator: pi+ + gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropPi+"); + gen.Particle = 'pi+' + gen.Energy = 100 * GeV + gen.Multiplicity = 2 + gen.Mask = 1 + gen.OutputLevel = Output.DEBUG + gen.PhiMin = 0 + gen.PhiMax = 0 + gen.ThetaMin = 1.61 + gen.ThetaMax = 1.61 + + kernel.generatorAction().adopt(gen) + # Install vertex smearing for this interaction + gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearPi+"); + gen.Mask = 1 + gen.Offset = (20*mm, 10*mm, 10*mm, 0*ns) + gen.Sigma = (4*mm, 1*mm, 1*mm, 0*ns) + kernel.generatorAction().adopt(gen) + """ + # Second particle generator: e- + gen = DDG4.GeneratorAction(kernel,"Geant4IsotropeGenerator/IsotropE-"); + gen.Particle = 'e-' + gen.Energy = 25 * GeV + gen.Multiplicity = 3 + gen.Mask = 2 + gen.OutputLevel = Output.DEBUG + kernel.generatorAction().adopt(gen) + # Install vertex smearing for this interaction + gen = DDG4.GeneratorAction(kernel,"Geant4InteractionVertexSmear/SmearE-"); + gen.Mask = 2 + gen.Offset = (-20*mm, -10*mm, -10*mm, 0*ns) + gen.Sigma = (12*mm, 8*mm, 8*mm, 0*ns) + kernel.generatorAction().adopt(gen) + #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + """ + # Merge all existing interaction records + gen = DDG4.GeneratorAction(kernel,"Geant4InteractionMerger/InteractionMerger") + gen.OutputLevel = 4 #generator_output_level + gen.enableUI() + kernel.generatorAction().adopt(gen) + + # Finally generate Geant4 primaries + gen = DDG4.GeneratorAction(kernel,"Geant4PrimaryHandler/PrimaryHandler") + gen.OutputLevel = Output.DEBUG #generator_output_level + gen.enableUI() + kernel.generatorAction().adopt(gen) + + # And handle the simulation particles. + part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler") + kernel.generatorAction().adopt(part) + #part.SaveProcesses = ['conv','Decay'] + part.SaveProcesses = ['Decay'] + part.MinimalKineticEnergy = 100*MeV + part.OutputLevel = Output.DEBUG # generator_output_level + part.enableUI() + user = DDG4.Action(kernel,"Geant4TCUserParticleHandler/UserParticleHandler") + user.TrackingVolume_Zmax = DDG4.EcalEndcap_zmin + user.TrackingVolume_Rmax = DDG4.EcalBarrel_rmin + user.enableUI() + part.adopt(user) + + # Setup global filters fur use in sensntive detectors + f1 = DDG4.Filter(kernel,'GeantinoRejectFilter/GeantinoRejector') + kernel.registerGlobalFilter(f1) + + # First the tracking detectors + seq,act = geant4.setupTracker('SiVertexBarrel') + act.OutputLevel = Output.ERROR + act.CollectSingleDeposits = False + seq,act = geant4.setupTracker('SiVertexEndcap') + act.OutputLevel = Output.ERROR + act.CollectSingleDeposits = False + + #seq,act = geant4.setupTracker('SiTrackerBarrel') + #seq,act = geant4.setupTracker('SiTrackerEndcap') + #seq,act = geant4.setupTracker('SiTrackerForward') + # Now the calorimeters + #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') + + # Now build the physics list: + phys = geant4.setupPhysics('QGSP_BERT') + phys.dump() + + kernel.configure() + kernel.initialize() + + #DDG4.setPrintLevel(Output.DEBUG) + kernel.run() + kernel.terminate() + +if __name__ == "__main__": + run() diff --git a/DDG4/include/DDG4/Geant4Action.h b/DDG4/include/DDG4/Geant4Action.h index ed5c20ac7..bc71dcacf 100644 --- a/DDG4/include/DDG4/Geant4Action.h +++ b/DDG4/include/DDG4/Geant4Action.h @@ -1,4 +1,3 @@ - #ifndef DD4HEP_DDG4_GEANT4ACTION_H #define DD4HEP_DDG4_GEANT4ACTION_H @@ -63,14 +62,14 @@ namespace DD4hep { class TypeName : public std::pair<std::string, std::string> { public: /// Default constructor - TypeName() - : std::pair<std::string, std::string>() { + TypeName() + : std::pair<std::string, std::string>() { } - TypeName(const std::pair<std::string, std::string>& c) - : std::pair<std::string, std::string>(c) { + TypeName(const std::pair<std::string, std::string>& c) + : std::pair<std::string, std::string>(c) { } - TypeName(const std::string& typ, const std::string& nam) - : std::pair<std::string, std::string>(typ, nam) { + 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); @@ -313,7 +312,7 @@ namespace DD4hep { /// Support of error messages. void error(const char* fmt, ...) const; /// Action to support error messages. - bool error(bool return_value, const char* fmt, ...) const; + bool return_error(bool return_value, const char* fmt, ...) const; /// Support of fatal messages. Throws exception void fatal(const char* fmt, ...) const; /// Support of exceptions: Print fatal message and throw runtime_error. diff --git a/DDG4/include/DDG4/Geant4Data.h b/DDG4/include/DDG4/Geant4Data.h index 4ef2994a0..8162cb9f2 100644 --- a/DDG4/include/DDG4/Geant4Data.h +++ b/DDG4/include/DDG4/Geant4Data.h @@ -107,8 +107,23 @@ namespace DD4hep { */ class Geant4HitData { public: + enum HitFlags { + HIT_KILLED_TRACK = 1<<0, + HIT_SECONDARY_TRACK = 1<<1, + HIT_NEW_TRACK = 1<<2, + HIT_STARTED_INSIDE = 1<<10, + HIT_STARTED_SURFACE = 1<<11, + HIT_STARTED_OUTSIDE = 1<<12, + HIT_ENDED_INSIDE = 1<<13, + HIT_ENDED_SURFACE = 1<<14, + HIT_ENDED_OUTSIDE = 1<<15 + }; /// cellID long long int cellID; + /// User flag to classify hits + long flag; + /// Original Geant 4 track identifier of the creating track (debugging) + long g4ID; /// User data extension if required dd4hep_ptr<DataExtension> extension; diff --git a/DDG4/include/DDG4/Geant4IsotropeGenerator.h b/DDG4/include/DDG4/Geant4IsotropeGenerator.h index e6a67e679..413f35003 100644 --- a/DDG4/include/DDG4/Geant4IsotropeGenerator.h +++ b/DDG4/include/DDG4/Geant4IsotropeGenerator.h @@ -31,6 +31,15 @@ namespace DD4hep { */ class Geant4IsotropeGenerator: public Geant4ParticleGenerator { protected: + /// Property: Minimal phi angular value + double m_phiMin; + /// Property: Maximal phi angular value + double m_phiMax; + /// Property: Minimal theta angular value + double m_thetaMin; + /// Property: Maximal theta angular value + double m_thetaMax; + /// Particle modification. Caller presets defaults to: ( direction = m_direction, momentum = m_energy) /** Use this function to implement isotrop guns, multiple guns etc. User must return a UNIT vector, which gets scaled with momentum. diff --git a/DDG4/include/DDG4/Geant4PhysicsList.h b/DDG4/include/DDG4/Geant4PhysicsList.h index 7d5adb070..34dccdcd5 100644 --- a/DDG4/include/DDG4/Geant4PhysicsList.h +++ b/DDG4/include/DDG4/Geant4PhysicsList.h @@ -24,6 +24,7 @@ // Forward declarations class G4VPhysicsConstructor; +class G4VUserPhysicsList; /// Namespace for the AIDA detector description toolkit namespace DD4hep { @@ -246,6 +247,8 @@ namespace DD4hep { virtual void constructProcess(Geant4UserPhysics* physics); /// begin-of-event callback virtual void constructParticles(Geant4UserPhysics* physics); + /// Extend physics list from factory: + G4VUserPhysicsList* extensionList() const; }; } // End namespace Simulation diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h index b607c336d..91c0faff2 100644 --- a/DDG4/include/DDG4/Geant4SensDetAction.h +++ b/DDG4/include/DDG4/Geant4SensDetAction.h @@ -251,14 +251,14 @@ namespace DD4hep { /** Combining the VolIDS of the complete geometry path (Geant4TouchableHistory) * from the current sensitive volume to the world volume */ - long long int volumeID(G4Step* step); + long long int volumeID(const G4Step* step); /// Returns the cellID of the sensitive volume corresponding to the step /** The CellID is the VolumeID + the local coordinates of the sensitive area. * Calculated by combining the VolIDS of the complete geometry path (Geant4TouchableHistory) * from the current sensitive volume to the world volume */ - long long int cellID(G4Step* step); + long long int cellID(const G4Step* step); }; diff --git a/DDG4/include/DDG4/Geant4StepHandler.h b/DDG4/include/DDG4/Geant4StepHandler.h index 40b6e080b..be86c7f1e 100644 --- a/DDG4/include/DDG4/Geant4StepHandler.h +++ b/DDG4/include/DDG4/Geant4StepHandler.h @@ -132,6 +132,9 @@ namespace DD4hep { int trkID() const { return track->GetTrackID(); } + int parentID() const { + return track->GetParentID(); + } double trkTime() const { return track->GetGlobalTime(); } @@ -141,6 +144,15 @@ namespace DD4hep { double trkKineEnergy() const { return track->GetKineticEnergy(); } + int trkStatus() const { + return track->GetTrackStatus(); + } + bool trkAlive() const { + return track->GetTrackStatus() == fAlive; + } + const G4VProcess* trkProcess() const { + return track->GetCreatorProcess(); + } const G4VTouchable* preTouchable() const { return pre->GetTouchable(); } @@ -154,6 +166,9 @@ namespace DD4hep { G4VPhysicalVolume* volume(const G4StepPoint* p) const { return p->GetTouchableHandle()->GetVolume(); } + G4VSolid* solid(const G4StepPoint* p) const { + return p->GetTouchableHandle()->GetSolid(); + } G4VPhysicalVolume* physvol(const G4StepPoint* p) const { return p->GetPhysicalVolume(); } @@ -206,6 +221,18 @@ namespace DD4hep { Position localToGlobal(const G4ThreeVector& local) const; /// Coordinate transformation to global coordinates in MM Position localToGlobal(double x, double y, double z) const; + + /// Coordinate transformation to local coordinates + Position globalToLocal(double x, double y, double z) const; + /// Coordinate transformation to local coordinates + Position globalToLocal(const Position& global) const; + /// Coordinate transformation to local coordinates + Position globalToLocal(const G4ThreeVector& global) const; + /// Coordinate transformation to local coordinates + G4ThreeVector globalToLocalG4(double x, double y, double z) const; + /// Coordinate transformation to local coordinates with G4 objects + G4ThreeVector globalToLocalG4(const G4ThreeVector& loc) const; + /// Apply BirksLaw double BirkAttenuation(const G4Step* aStep) const { diff --git a/DDG4/plugins/Geant4PhysicsLists.cpp b/DDG4/plugins/Geant4PhysicsLists.cpp deleted file mode 100644 index 5ee3af210..000000000 --- a/DDG4/plugins/Geant4PhysicsLists.cpp +++ /dev/null @@ -1,169 +0,0 @@ -// $Id: $ -//========================================================================== -// AIDA Detector description implementation for LCD -//-------------------------------------------------------------------------- -// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) -// All rights reserved. -// -// For the licensing terms see $DD4hepINSTALL/LICENSE. -// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. -// -// Author : M.Frank -// -//========================================================================== - -// Framework include files -#include "DDG4/Factories.h" - -// Geant4 physics lists -#include "G4Version.hh" -#include "G4DecayPhysics.hh" - -#if G4VERSION_NUMBER>=960 -#define GEANT4_9_6 -#endif - -// ====================================================================== -// Predefined physics lists -// -// Usage: -// -// <physicslist name="Geant4PhysicsList/MyPhysics.0"> -// <list name="TQGSP_FTFP_BERT_95"/> -// </physicslist> -// -// -// ====================================================================== -#include "DDG4/Geant4UserPhysicsList.h" -namespace { - struct EmptyPhysics : public G4VModularPhysicsList { - EmptyPhysics(int) {} - virtual ~EmptyPhysics() {} - virtual void ConstructProcess() {} - virtual void ConstructParticle() {} - }; -} -DECLARE_GEANT4_PHYSICS_LIST(EmptyPhysics) - -// Physics constructors from source/physics_lists - -#if G4VERSION_NUMBER>=960 and G4VERSION_NUMBER<1000 -#include "CHIPS.hh" -DECLARE_GEANT4_PHYSICS_LIST(CHIPS) -#endif - -#ifdef GEANT4_9_6 - -#if G4VERSION_NUMBER < 1000 -#include "CHIPS_HP.hh" -DECLARE_GEANT4_PHYSICS_LIST(CHIPS_HP) -#include "QGSP_BERT_95.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_95) -#include "QGSP_BERT_95XS.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_95XS) -#endif - -#include "FTFP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT) -#include "FTFP_BERT_HP.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT_HP) - -#if G4VERSION_NUMBER < 1000 - -// extra EM physics lists only via physicsListFactory -#include "FTFP_BERT_EMV.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT_EMV) -#include "FTFP_BERT_EMX.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT_EMX) -#include "LHEP_EMV.hh" -DECLARE_GEANT4_PHYSICS_LIST(LHEP_EMV) -#include "QGSP_BERT_EMV.hh" - DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_EMV) -#include "QGSP_BERT_EMX.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_EMX) -#include "QGSP_BIC_EMY.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BIC_EMY) -#include "LHEP.hh" -DECLARE_GEANT4_PHYSICS_LIST(LHEP) -#include "QGSC_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSC_BERT) -#include "QGSC_CHIPS.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSC_CHIPS) -#include "QGSP_BERT_CHIPS.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_CHIPS) -#include "QGSP_BERT_EMV.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_EMV) -#include "QGSP_BERT_EMX.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_EMX) - -#include "QGSP_BERT_NOLEP.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_NOLEP) -#include "QGSP_BERT_TRV.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_TRV) -#include "QGSP.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP) -#include "QGSP_QEL.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_QEL) -#include "HadronPhysicsCHIPS.hh" - -#endif - -#include "FTFP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT) -#include "FTFP_BERT_TRV.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT_TRV) -#include "LBE.hh" -//DECLARE_GEANT4_PHYSICS_LIST(LBE) takes no verbosity arg! -#include "QBBC.hh" -DECLARE_GEANT4_PHYSICS_LIST(QBBC) -#include "QGS_BIC.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGS_BIC) -#include "QGSP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT) -#endif - -#if G4VERSION_NUMBER>=1000 -#include "NuBeam.hh" -DECLARE_GEANT4_PHYSICS_LIST(NuBeam) -#include "FTFP_INCLXX.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_INCLXX) -#include "QGSP_BERT_HP.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT_HP) -#include "QGSP_BIC_AllHP.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BIC_AllHP) -#include "QGSP_INCLXX.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_INCLXX) -#endif - -#include "FTFP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT) -#include "FTFP_BERT_HP.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT_HP) -#include "QGSP_INCLXX.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_INCLXX) - -#include "FTF_BIC.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTF_BIC) -#include "FTFP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT) -#include "FTFP_BERT_TRV.hh" -DECLARE_GEANT4_PHYSICS_LIST(FTFP_BERT_TRV) -#include "LBE.hh" -//DECLARE_GEANT4_PHYSICS_LIST(LBE) takes no verbosity arg! -#include "QBBC.hh" -DECLARE_GEANT4_PHYSICS_LIST(QBBC) -#include "QGS_BIC.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGS_BIC) -#include "QGSP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BERT) -#include "QGSP_BIC_HP.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BIC_HP) -#include "QGSP_BIC.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_BIC) -#include "QGSP_FTFP_BERT.hh" -DECLARE_GEANT4_PHYSICS_LIST(QGSP_FTFP_BERT) - -#if 0 -#include ".hh" -DECLARE_GEANT4_PHYSICS_LIST() -#endif diff --git a/DDG4/src/Geant4ROOTDump.cpp b/DDG4/plugins/Geant4ROOTDump.cpp similarity index 99% rename from DDG4/src/Geant4ROOTDump.cpp rename to DDG4/plugins/Geant4ROOTDump.cpp index 0fefa91ba..e997b2dea 100644 --- a/DDG4/src/Geant4ROOTDump.cpp +++ b/DDG4/plugins/Geant4ROOTDump.cpp @@ -1,4 +1,4 @@ -// $Id: $ +// $Id$ //========================================================================== // AIDA Detector description implementation for LCD //-------------------------------------------------------------------------- diff --git a/DDG4/plugins/Geant4TrackerWeightedSD.cpp b/DDG4/plugins/Geant4TrackerWeightedSD.cpp new file mode 100644 index 000000000..b67f340b8 --- /dev/null +++ b/DDG4/plugins/Geant4TrackerWeightedSD.cpp @@ -0,0 +1,393 @@ +// $Id$ +//========================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include "DD4hep/LCDD.h" +#include "DD4hep/DD4hepUnits.h" +#include "DDG4/Geant4SensDetAction.inl" +#include "DDG4/Geant4SteppingAction.h" +#include "DDG4/Geant4TrackingAction.h" +#include "DDG4/Geant4EventAction.h" +#include "G4OpticalPhoton.hh" +#include "G4VProcess.hh" +#include "G4Event.hh" +#include "G4VSolid.hh" + +#include "DDRec/Surface.h" +#include "DDRec/DetectorSurfaces.h" +#include "DDRec/SurfaceManager.h" +#include "DDRec/SurfaceHelper.h" + +#include <map> +#include <limits> +#include <sstream> + +using namespace std; + +/// Namespace for the AIDA detector description toolkit +namespace DD4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace Simulation { + + using namespace DDRec; + using namespace Geometry; + + /// Geant4 sensitive detector combining all deposits of one G4Track within one sensitive element. + /** + * Geant4SensitiveAction<TrackerWeighted> + * + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + struct TrackerWeighted { + typedef Geant4HitCollection HitCollection; + + /// Enumeration to define the calculation of the position of the energy deposit + enum { + POSITION_WEIGHTED = 1, // Use energy weights to define the position of the energy deposit + POSITION_MIDDLE = 2, // Set the hit position between the step pre and post point + POSITION_PREPOINT = 3, // Set the hit position to the position of the step pre point + POSITION_POSTPOINT = 4 // Set the hit position to the position of the step post point + }; + + Geant4Tracker::Hit pre, post; + Position mean_pos; + Geant4Sensitive* sensitive; + G4VSensitiveDetector* thisSD; + double distance_to_inside; + double distance_to_outside; + double mean_time; + double e_cut; + int current, parent; + int combined; + int hit_position_type; + int hit_flag; + int g4ID; + EInside last_inside; + long long int cell; + bool single_deposit_mode; + TrackerWeighted() : pre(), post(), sensitive(0), thisSD(0), + distance_to_inside(0.0), distance_to_outside(0.0), mean_time(0.0), + e_cut(0.0), current(-1), parent(0), combined(0), + hit_position_type(POSITION_MIDDLE), hit_flag(0), g4ID(0), cell(0), + single_deposit_mode(false) + { + } + + /// Clear collected information and restart for new hit + TrackerWeighted& clear() { + mean_pos.SetXYZ(0,0,0); + distance_to_inside = 0; + distance_to_outside = 0; + mean_time = 0; + post.clear(); + pre.clear(); + current = -1; + parent = -1; + combined = 0; + cell = 0; + hit_flag = 0; + g4ID = 0; + last_inside = kOutside; + return *this; + } + + /// Start a new hit + TrackerWeighted& start(const G4Step* step, const G4StepPoint* point) { + clear(); + pre.storePoint(step,point); + pre.truth.deposit = 0.0; + post.truth.deposit = 0.0; + current = pre.truth.trackID; + sensitive->mark(step->GetTrack()); + post = pre; + parent = step->GetTrack()->GetParentID(); + g4ID = step->GetTrack()->GetTrackID(); + return *this; + } + + /// Update energy and track information during hit info accumulation + TrackerWeighted& update(const G4Step* step) { + post.storePoint(step,step->GetPostStepPoint()); + pre.truth.deposit += post.truth.deposit; + mean_pos.SetX(mean_pos.x()+post.position.x()*post.truth.deposit); + mean_pos.SetY(mean_pos.y()+post.position.y()*post.truth.deposit); + mean_pos.SetZ(mean_pos.z()+post.position.z()*post.truth.deposit); + mean_time += post.truth.time*post.truth.deposit; + if ( 0 == cell ) { + cell = sensitive->cellID(step); + if ( 0 == cell ) { + cell = sensitive->volumeID(step) ; + sensitive->except("+++ Invalid CELL ID for hit!"); + } + } + ++combined; + return *this; + } + + /// Helper function to decide if the hit has to be extracted and saved in the collection + bool mustSaveTrack(const G4Track* tr) const { + return current > 0 && current != tr->GetTrackID(); + } + + /// Extract hit information and add the created hit to the collection + void extractHit(EInside ended) { + Geant4HitCollection* collection = sensitive->collection(0); + extractHit(collection, ended); + } + + TrackerWeighted& calc_dist_out(const G4VSolid* solid) { + Position v(pre.momentum.unit()), &p=post.position; + double dist = solid->DistanceToOut(G4ThreeVector(p.X(),p.Y(),p.Z()), + G4ThreeVector(v.X(),v.Y(),v.Z())); + distance_to_outside = dist; + return *this; + } + + TrackerWeighted& calc_dist_in(const G4VSolid* solid) { + Position v(pre.momentum.unit()), &p=pre.position; + double dist = solid->DistanceToOut(G4ThreeVector(p.X(),p.Y(),p.Z()), + G4ThreeVector(v.X(),v.Y(),v.Z())); + distance_to_inside = dist; + return *this; + } + + void extractHit(Geant4HitCollection* collection, EInside ended) { + double deposit = pre.truth.deposit; + if ( current != -1 && deposit/CLHEP::keV > 0 ) { + Position pos; + Momentum mom = 0.5 * (pre.momentum + post.momentum); + double time = mean_time / deposit; + double path = (post.position - pre.position).R(); + char dist_in[64], dist_out[64]; + + switch(hit_position_type) { + case POSITION_WEIGHTED: + pos = mean_pos / deposit; + break; + case POSITION_PREPOINT: + pos = pre.position; + break; + case POSITION_POSTPOINT: + pos = post.position; + break; + case POSITION_MIDDLE: + default: + pos = (post.position + pre.position) / 2.0; + break; + } + + if ( ended == kSurface || distance_to_outside < numeric_limits<float>::epsilon() ) + hit_flag |= Geant4Tracker::Hit::HIT_ENDED_SURFACE; + else if ( ended == kInside ) + hit_flag |= Geant4Tracker::Hit::HIT_ENDED_INSIDE; + else if ( ended == kOutside ) + hit_flag |= Geant4Tracker::Hit::HIT_ENDED_OUTSIDE; + + Geant4Tracker::Hit* hit = new Geant4Tracker::Hit(pre.truth.trackID, + pre.truth.pdgID, + deposit,time); + hit->flag = hit_flag; + hit->position = pos; + hit->momentum = mom; + hit->length = path; + hit->cellID = cell; + hit->g4ID = g4ID; + + dist_in[0] = dist_out[0] = 0; + if ( !(hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE) ) + ::snprintf(dist_in,sizeof(dist_in)," [%.2e um]",distance_to_inside/CLHEP::um); + if ( !(hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE) ) + ::snprintf(dist_out,sizeof(dist_out)," [%.2e um]",distance_to_outside/CLHEP::um); + sensitive->print("+++ G4Track:%5d CREATE hit[%03d]:%3d deps E:" + " %.2e keV Pos:%7.2f %7.2f %7.2f [mm] Start:%s%s%s%s End:%s%s%s%s", + pre.truth.trackID,int(collection->GetSize()), + combined,pre.truth.deposit/CLHEP::keV, + pos.X()/CLHEP::mm,pos.Y()/CLHEP::mm,pos.Z()/CLHEP::mm, + hit_flag&Geant4Tracker::Hit::HIT_STARTED_SURFACE ? "SURFACE" : "", + hit_flag&Geant4Tracker::Hit::HIT_STARTED_OUTSIDE ? "OUTSIDE" : "", + hit_flag&Geant4Tracker::Hit::HIT_STARTED_INSIDE ? "INSIDE " : "", + dist_in, + hit_flag&Geant4Tracker::Hit::HIT_ENDED_SURFACE ? "SURFACE" : "", + hit_flag&Geant4Tracker::Hit::HIT_ENDED_OUTSIDE ? "OUTSIDE" : "", + hit_flag&Geant4Tracker::Hit::HIT_ENDED_INSIDE ? "INSIDE " : "", + dist_out); + collection->add(hit); + } + clear(); + } + + /// Method for generating hit(s) using the information of G4Step object. + G4bool process(const G4Step* step, G4TouchableHistory* ) { + Geant4StepHandler h(step); + G4VSolid* preSolid = h.solid(h.pre); + G4VSolid* postSolid = h.solid(h.post); + G4ThreeVector local_pre = h.globalToLocalG4(h.prePosG4()); + G4ThreeVector local_post = h.globalToLocalG4(h.postPosG4()); + EInside pre_inside = preSolid->Inside(local_pre); + EInside post_inside = postSolid->Inside(local_post); + + const void* postPV = h.postVolume(); + const void* prePV = h.preVolume(); + const void* postSD = h.postSD(); + const void* preSD = h.preSD(); + G4VSolid* solid = (preSD == thisSD) ? preSolid : postSolid; + // 1) Track killed inside SD: trace incomplete. This deposition must be added as well. + if ( current == h.trkID() && !h.trkAlive() ) { + hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK; + update(step).calc_dist_out(solid).extractHit(post_inside); + return true; + } + // 2) Track leaving SD volume. Sensitive detector changed. Store hit. + else if ( current == h.trkID() && postSD != thisSD ) { + update(step).calc_dist_out(solid).extractHit(kOutside); + return true; + } + // 3) Track leaving SD volume. Store hit. + else if ( current == h.trkID() && postSD == thisSD && post_inside == kSurface ) { + update(step).calc_dist_out(solid).extractHit(kSurface); + return true; + } + // 4) Track leaving SD volume. Store hit. + else if ( current == h.trkID() && postSD == thisSD && post_inside == kOutside ) { + update(step).calc_dist_out(solid).extractHit(post_inside); + return true; + } + // 5) Normal update: either intermediate deposition or track is going to be killed. + else if ( current == h.trkID() && postSD == thisSD && post_inside == kInside ) { + last_inside = post_inside; + update(step).calc_dist_out(solid); + return true; + } + // 6) Track from secondary created in SD volume. Store hit from previous. + // --> New hit will be created, to whom also this deposition belongs + else if ( current != h.trkID() && current >= 0 ) { + extractHit(last_inside); + } + + // If the hit got extracted, a new one must be set up + if ( current < 0 ) { + EInside inside = pre_inside; + // Track entering SD volume + if ( preSD != thisSD ) { + start(step, h.post); + inside = post_inside; + sensitive->print("++++++++++ Using POST step volume to start hit -- dubious ?"); + } + else { + start(step, h.pre); + } + calc_dist_in(solid); + if ( inside == kSurface ) + hit_flag |= Geant4Tracker::Hit::HIT_STARTED_SURFACE; + else if ( inside == kInside ) + hit_flag |= Geant4Tracker::Hit::HIT_STARTED_INSIDE; + else if ( inside == kOutside ) + hit_flag |= Geant4Tracker::Hit::HIT_STARTED_OUTSIDE; + + // New (secondary) track created by some process starting inside the volume + if ( inside == kInside ) { + hit_flag |= Geant4Tracker::Hit::HIT_SECONDARY_TRACK; + } + } + + // Update Data + last_inside = post_inside; + update(step); + calc_dist_out(solid); + + // Track killed inside SD: trace incomplete. This deposition must be added as well. + if ( !h.trkAlive() ) { + hit_flag |= Geant4Tracker::Hit::HIT_KILLED_TRACK; + extractHit(post_inside); + } + // Avoid danglich hits if the track leaves the sensitive volume + else if ( post_inside == kSurface ) { + extractHit(post_inside); + } + // Avoid danglich hits if the track leaves the sensitive volume + else if ( thisSD == preSD && (preSD != postSD || prePV != postPV) ) { + extractHit(post_inside); + } + // This should simply not happen! + else if ( thisSD == postSD && (preSD != postSD || prePV != postPV) ) { + sensitive->error("+++++ WRONG!!! Extract. How did we get here?"); + extractHit(post_inside); + } + // In single hit mode we MUST write the hit after update + else if ( single_deposit_mode ) { + extractHit(post_inside); + } + + return true; + } + + /// Post-event action callback + void endEvent(const G4Event*) { + // We need to add the possibly last added hit to the collection here. + // otherwise the last hit would be assigned to the next event and the + // MC truth would be screwed. + // + if ( current > 0 ) { + Geant4HitCollection* coll = sensitive->collection(0); + sensitive->print("++++++++++ Found dangling hit: Is the hit extraction logic correct?"); + extractHit(coll,last_inside); + } + } + /// Pre event action callback + void startEvent(const G4Event* event) { + thisSD = dynamic_cast<G4VSensitiveDetector*>(&sensitive->detector()); + if ( event ) { + sensitive->print("++++++++++++++++++++++++++ START EVENT %d",event->GetEventID()); + } + } + }; + + /// Initialization overload for specialization + template <> void Geant4SensitiveAction<TrackerWeighted>::initialize() { + declareProperty("HitPostionCombination", m_userData.hit_position_type); + declareProperty("CollectSingleDeposits", m_userData.single_deposit_mode); + eventAction().callAtBegin(&m_userData,&TrackerWeighted::startEvent); + eventAction().callAtEnd(&m_userData,&TrackerWeighted::endEvent); + m_userData.e_cut = m_sensitive.energyCutoff(); + m_userData.sensitive = this; + } + + /// Define collections created by this sensitivie action object + template <> void Geant4SensitiveAction<TrackerWeighted>::defineCollections() { + m_collectionID = defineCollection<Geant4Tracker::Hit>(m_sensitive.readout().name()); + } + + /// Method for generating hit(s) using the information of G4Step object. + template <> void Geant4SensitiveAction<TrackerWeighted>::clear(G4HCofThisEvent*) { + m_userData.clear(); + } + + /// Method for generating hit(s) using the information of G4Step object. + template <> G4bool + Geant4SensitiveAction<TrackerWeighted>::process(G4Step* step, G4TouchableHistory* history) { + return m_userData.process(step, history); + } + + typedef Geant4SensitiveAction<TrackerWeighted> Geant4TrackerWeightedAction; + } +} + +using namespace DD4hep::Simulation; + +#include "DDG4/Factories.h" +DECLARE_GEANT4SENSITIVE(Geant4TrackerWeightedAction) diff --git a/DDG4/reco/Geant4SurfaceTest.cpp b/DDG4/reco/Geant4SurfaceTest.cpp new file mode 100644 index 000000000..eab2b3b1f --- /dev/null +++ b/DDG4/reco/Geant4SurfaceTest.cpp @@ -0,0 +1,197 @@ +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// \author Markus Frank +// \date 2015-10-14 +// +//========================================================================== +// $Id$ +#ifndef DD4HEP_DDG4_GEANT4SURFACETEST_H +#define DD4HEP_DDG4_GEANT4SURFACETEST_H + +// Framework include files +//#include "Geant4SurfaceTest.h" +#include "DDG4/Geant4EventAction.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 { + + /** @class Geant4SurfaceTest Geant4SurfaceTest.h reco/Geant4SurfaceTest.h + */ + class Geant4SurfaceTest : public Geant4EventAction { + public: + /// Standard constructor + Geant4SurfaceTest(Geant4Context* context, const std::string& nam); + /// Destructor + virtual ~Geant4SurfaceTest( ); + + /// Begin-of-event callback + virtual void begin(const G4Event* event); + /// End-of-event callback + virtual void end(const G4Event* event); + + protected: + + private: + + }; + } // End namespace Simulation +} // End namespace DD4hep + +#endif // DD4HEP_DDG4_GEANT4SURFACETEST_H + +// Framework include files +#include "DD4hep/LCDD.h" +#include "DD4hep/DD4hepUnits.h" +#include "DD4hep/InstanceCount.h" + +#include "DDG4/Factories.h" +#include "DDG4/Geant4Data.h" +#include "DDG4/Geant4Context.h" +#include "DDG4/Geant4SensDetAction.h" +#include "DDG4/Geant4HitCollection.h" + +#include "DDRec/Surface.h" +#include "DDRec/DetectorSurfaces.h" +#include "DDRec/SurfaceManager.h" +#include "DDRec/SurfaceHelper.h" + +// Geant 4 includes +#include "G4VHitsCollection.hh" +#include "G4HCofThisEvent.hh" +#include "G4Event.hh" + +// C/C++ include files +#include <stdexcept> +#include <map> +#include <sstream> + +using namespace std; +using namespace DD4hep::DDRec; +using namespace DD4hep::Geometry; +using namespace DD4hep::Simulation; + +DECLARE_GEANT4ACTION(Geant4SurfaceTest) + +/// Standard constructor, initializes variables +Geant4SurfaceTest::Geant4SurfaceTest(Geant4Context* ctxt, const std::string& nam) +: Geant4EventAction(ctxt, nam) +{ + InstanceCount::increment(this); +} + +/// Default Destructor +Geant4SurfaceTest::~Geant4SurfaceTest() { + InstanceCount::decrement(this); +} + +/// Begin-of-event callback +void Geant4SurfaceTest::begin(const G4Event* event) { + if ( event ) { + } +} + +/// End-of-event callback +void Geant4SurfaceTest::end(const G4Event* evt) { + stringstream sst; + LCDD& lcdd = context()->lcdd(); + SurfaceManager& surfMan = *lcdd.extension< SurfaceManager >() ; + const SurfaceMap& surfMap = *surfMan.map( "world" ) ; + G4HCofThisEvent* hce = evt->GetHCofThisEvent(); + + if ( !hce ) { + error("Event:%d: No Geant4 event record (HCE) found.",evt->GetEventID()); + return; + } + + for(int j=0, nCol=hce->GetNumberOfCollections(); j<nCol; ++j) { + G4VHitsCollection* hc = hce->GetHC(j); + Geant4HitCollection* coll = dynamic_cast<Geant4HitCollection*>(hc); + + if ( !(coll && coll->type().type == typeid(Geant4Tracker::Hit)) ) { + continue; + } + + SensitiveDetector sd = coll->sensitive()->sensitiveDetector(); + if ( !sd.isValid() ) { + error("Event:%d: Collection:%s Invalid sensitive detector.", + evt->GetEventID(),hc->GetName().c_str()); + continue; + } + IDDescriptor idDesc = sd.readout().idSpec(); + const string idStr = idDesc.fieldDescription(); + BitField64 idDecoder(idStr); + info("Event: %d [%s: %d hits] Desc:%s", + evt->GetEventID(),hc->GetName().c_str(),coll->GetSize(),idStr.c_str()); + for(size_t i=0, nHit=coll->GetSize(); i<nHit; ++i) { + Geant4Tracker::Hit* h = coll->hit(i); + const Position& pos = h->position; + int trackID = h->truth.trackID; + idDecoder.setValue( h->cellID ) ; + SurfaceMap::const_iterator si = surfMap.find(h->cellID); + ISurface* surf = (si != surfMap.end() ? si->second : 0); + + // ===== test that we have a surface with the correct ID for every hit ====================== + if ( !surf ) { + error("FAILED: Hit:%ld Track:%d No surface found cell id: %016llX",i,trackID,h->cellID); + continue; + } + info("PASSED: Hit:%ld Track:%d Surface found for cell: %016llX",i,trackID,h->cellID); + + Vector3D hit_point(pos.x()*dd4hep::mm,pos.y()*dd4hep::mm,pos.z()*dd4hep::mm); + bool isInside = surf->insideBounds(hit_point); + string flag = ""; + if ( h->flag & Geant4Tracker::Hit::HIT_SECONDARY_TRACK ) flag += " TRK:SECONDARY"; + if ( h->flag & Geant4Tracker::Hit::HIT_KILLED_TRACK ) flag += " TRK:KILLED"; + if ( h->flag & Geant4Tracker::Hit::HIT_NEW_TRACK ) flag += " NEW_TRACK"; + if ( h->flag & Geant4Tracker::Hit::HIT_STARTED_OUTSIDE ) flag += " START:OUTSIDE"; + if ( h->flag & Geant4Tracker::Hit::HIT_STARTED_INSIDE ) flag += " START:INSIDE"; + if ( h->flag & Geant4Tracker::Hit::HIT_STARTED_SURFACE ) flag += " START:SURFACE"; + if ( h->flag & Geant4Tracker::Hit::HIT_ENDED_OUTSIDE ) flag += " END:OUTSIDE"; + if ( h->flag & Geant4Tracker::Hit::HIT_ENDED_INSIDE ) flag += " END:INSIDE"; + if ( h->flag & Geant4Tracker::Hit::HIT_ENDED_SURFACE ) flag += " END:SURFACE"; + + sst.str(""); + if ( isInside ) { + sst << hc->GetName() << " " << "PASSED: Track:" << trackID << " Hit:" << i + //<< " Point " << hit_point + << " is ON SURFACE. " + << "(FLAG:" << h->flag << flag << ")."; + print(sst.str().c_str()); + } + else { + double dist = surf->distance(hit_point)/dd4hep::mm; + sst.str(""); + sst << hc->GetName() << " FAILED: Track:" << trackID << " Hit:" << i + << " Point: " << hit_point; + error(sst.str().c_str()); + //sst.str(""); + //sst << -> Surface:" << *surf; + //error(sst.str().c_str()); + sst.str(""); + sst << " -> Hit is NOT ON SURFACE (Flag:" << h->flag << flag << ")" + << " Distance to surface:" << dist << " mm."; + error(sst.str().c_str()); + // No need for further tests, if this one failed.... + continue; + } + // ====== Test to show that moved hit points are outside their surface + hit_point = hit_point + 1e-3 * surf->normal() ; + sst.str("") ; + sst << "Track:" << trackID << " Hit:" << i << " Moved:" << hit_point; + if ( (isInside=surf->insideBounds(hit_point)) ) { + error("FAILED: %s is ON SURFACE (SHOULD NOT BE)",sst.str().c_str()); + continue; + } + info("PASSED: %s is NOT ON SURFACE (SHOULD NOT BE)",sst.str().c_str()); + } + } +} diff --git a/DDG4/src/Geant4Action.cpp b/DDG4/src/Geant4Action.cpp index 9e72d1cf7..b6bab5291 100644 --- a/DDG4/src/Geant4Action.cpp +++ b/DDG4/src/Geant4Action.cpp @@ -238,7 +238,7 @@ void Geant4Action::error(const char* fmt, ...) const { } /// Action to support error messages. -bool Geant4Action::error(bool return_value, const char* fmt, ...) const { +bool Geant4Action::return_error(bool return_value, const char* fmt, ...) const { va_list args; va_start(args, fmt); DD4hep::printout(DD4hep::ERROR, m_name, fmt, args); diff --git a/DDG4/src/Geant4Data.cpp b/DDG4/src/Geant4Data.cpp index 7e6d877a1..7e6229f88 100644 --- a/DDG4/src/Geant4Data.cpp +++ b/DDG4/src/Geant4Data.cpp @@ -55,7 +55,7 @@ DataExtension::~DataExtension() { /// Default constructor Geant4HitData::Geant4HitData() -: cellID(0), extension() { + : cellID(0), flag(0), g4ID(0), extension() { InstanceCount::increment(this); } diff --git a/DDG4/src/Geant4Exec.cpp b/DDG4/src/Geant4Exec.cpp index 4a9c59c9d..5c94ec606 100644 --- a/DDG4/src/Geant4Exec.cpp +++ b/DDG4/src/Geant4Exec.cpp @@ -320,6 +320,7 @@ using namespace DD4hep::Simulation; // Geant4 include files #include "G4RunManager.hh" +#include "G4PhysListFactory.hh" /// Configure the simulation int Geant4Exec::configure(Geant4Kernel& kernel) { @@ -348,18 +349,12 @@ int Geant4Exec::configure(Geant4Kernel& kernel) { Geant4DetectorConstruction* detector = Geant4DetectorConstruction::instance(kernel); runManager.SetUserInitialization(detector); - G4VUserPhysicsList* physics = 0; Geant4PhysicsListActionSequence* seq = kernel.physicsList(false); - if (seq) { - string nam = kernel.physicsList().properties()["extends"].value<string>(); - if (nam.empty()) - nam = "EmptyPhysics"; - physics = PluginService::Create<G4VUserPhysicsList*>(nam, seq, int(1)); - } - else { - printout(INFO, "Geant4Exec", "+++ Using Geant4 physics constructor QGSP_BERT"); - physics = PluginService::Create<G4VUserPhysicsList*>(string("QGSP_BERT"), seq, int(1)); + if ( 0 == seq ) { + seq = kernel.physicsList(true); + seq->property("Extends").set<string>("QGSP_BERT"); } + G4VUserPhysicsList* physics = seq->extensionList(); if (0 == physics) { throw runtime_error("Panic! No valid user physics list present!"); } diff --git a/DDG4/src/Geant4GeneratorActionInit.cpp b/DDG4/src/Geant4GeneratorActionInit.cpp index 05ccf9c68..107ccd949 100644 --- a/DDG4/src/Geant4GeneratorActionInit.cpp +++ b/DDG4/src/Geant4GeneratorActionInit.cpp @@ -46,7 +46,7 @@ void Geant4GeneratorActionInit::begin(const G4Run* run) { /// End-run action callback void Geant4GeneratorActionInit::end(const G4Run* /* run */) { - printP1("+++ Funished run %d after %d events (%d events in total)",m_run,m_evtRun,m_evtTotal); + printP1("+++ Finished run %d after %d events (%d events in total)",m_run,m_evtRun,m_evtTotal); m_evtRun = 0; m_run = 0; } diff --git a/DDG4/src/Geant4IsotropeGenerator.cpp b/DDG4/src/Geant4IsotropeGenerator.cpp index b98a3df6d..c49712656 100644 --- a/DDG4/src/Geant4IsotropeGenerator.cpp +++ b/DDG4/src/Geant4IsotropeGenerator.cpp @@ -26,6 +26,10 @@ Geant4IsotropeGenerator::Geant4IsotropeGenerator(Geant4Context* ctxt, const stri : Geant4ParticleGenerator(ctxt, nam) { InstanceCount::increment(this); + declareProperty("PhiMin", m_phiMin = 0.0); + declareProperty("PhiMax", m_phiMax = 2.0*M_PI); + declareProperty("ThetaMin", m_thetaMin = 0.0); + declareProperty("ThetaMax", m_thetaMax = M_PI); } /// Default destructor @@ -37,8 +41,8 @@ Geant4IsotropeGenerator::~Geant4IsotropeGenerator() { void Geant4IsotropeGenerator::getParticleDirection(int, ROOT::Math::XYZVector& direction, double& momentum) const { Geant4Event& evt = context()->event(); Geant4Random& rnd = evt.random(); - double phi = 2*M_PI*rnd.rndm(); - double theta = M_PI*rnd.rndm(); + double phi = m_phiMin+(m_phiMax-m_phiMin)*rnd.rndm(); + double theta = m_thetaMin+(m_thetaMax-m_thetaMin)*rnd.rndm(); double x1 = std::sin(theta)*std::cos(phi); double x2 = std::sin(theta)*std::sin(phi); double x3 = std::cos(theta); diff --git a/DDG4/src/Geant4PhysicsList.cpp b/DDG4/src/Geant4PhysicsList.cpp index daef90924..7504b04b4 100644 --- a/DDG4/src/Geant4PhysicsList.cpp +++ b/DDG4/src/Geant4PhysicsList.cpp @@ -15,14 +15,16 @@ // Framework include files #include "DDG4/Geant4PhysicsList.h" #include "DDG4/Geant4UIMessenger.h" +#include "DDG4/Geant4UserPhysicsList.h" #include "DD4hep/InstanceCount.h" #include "DD4hep/Printout.h" #include "DD4hep/Plugins.h" // Geant4 include files #include "G4VPhysicsConstructor.hh" -#include "G4ParticleTable.hh" +#include "G4PhysListFactory.hh" #include "G4ProcessManager.hh" +#include "G4ParticleTable.hh" #include "G4VProcess.hh" #include "G4Decay.hh" @@ -35,6 +37,14 @@ using namespace DD4hep; using namespace DD4hep::Simulation; namespace { + + struct EmptyPhysics : public G4VModularPhysicsList { + EmptyPhysics() {} + virtual ~EmptyPhysics() {} + virtual void ConstructProcess() {} + virtual void ConstructParticle() {} + }; + void _findDef(const string& expression, vector<G4ParticleDefinition*>& results) { string exp = expression; //'^'+expression+"$"; G4ParticleTable* pt = G4ParticleTable::GetParticleTable(); @@ -261,6 +271,14 @@ Geant4PhysicsListActionSequence::~Geant4PhysicsListActionSequence() { InstanceCount::decrement(this); } +/// Extend physics list from factory: +G4VUserPhysicsList* Geant4PhysicsListActionSequence::extensionList() const { + G4VModularPhysicsList* physics = ( m_extends.empty() ) + ? new EmptyPhysics() + : G4PhysListFactory().GetReferencePhysList(m_extends); + return physics; +} + /// Install command control messenger if wanted void Geant4PhysicsListActionSequence::installCommandMessenger() { control()->addCall("dump", "Dump content of " + name(), Callback(this).make(&Geant4PhysicsListActionSequence::dump)); diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp index b75d9b7da..04a9c098b 100644 --- a/DDG4/src/Geant4SensDetAction.cpp +++ b/DDG4/src/Geant4SensDetAction.cpp @@ -193,7 +193,7 @@ void Geant4Sensitive::mark(const G4Step* step) const { } /// Returns the volumeID of the sensitive volume corresponding to the step -long long int Geant4Sensitive::volumeID(G4Step* s) { +long long int Geant4Sensitive::volumeID(const G4Step* s) { Geant4StepHandler step(s); Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager(); VolumeID id = volMgr.volumeID(step.preTouchable()); @@ -201,7 +201,7 @@ long long int Geant4Sensitive::volumeID(G4Step* s) { } /// Returns the cellID(volumeID+local coordinate encoding) of the sensitive volume corresponding to the step -long long int Geant4Sensitive::cellID(G4Step* s) { +long long int Geant4Sensitive::cellID(const G4Step* s) { StepHandler h(s); Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager(); VolumeID volID = volMgr.volumeID(h.preTouchable()); diff --git a/DDG4/src/Geant4StepHandler.cpp b/DDG4/src/Geant4StepHandler.cpp index c5886dcbf..2410a0635 100644 --- a/DDG4/src/Geant4StepHandler.cpp +++ b/DDG4/src/Geant4StepHandler.cpp @@ -77,8 +77,36 @@ Position Geant4StepHandler::localToGlobal(double x, double y, double z) const /// Coordinate transformation to global coordinates Position Geant4StepHandler::localToGlobal(const G4ThreeVector& loc) const { - G4StepPoint* sp = step->GetPreStepPoint(); - G4TouchableHandle t = sp->GetTouchableHandle(); + G4TouchableHandle t = step->GetPreStepPoint()->GetTouchableHandle(); G4ThreeVector p = t->GetHistory()->GetTopTransform().Inverse().TransformPoint(loc); return Position(p.x(),p.y(),p.z()); } + +/// Coordinate transformation to local coordinates +Position Geant4StepHandler::globalToLocal(double x, double y, double z) const { + G4ThreeVector p = globalToLocalG4(G4ThreeVector(x,y,z)); + return Position(p.x(),p.y(),p.z()); +} + +/// Coordinate transformation to local coordinates +Position Geant4StepHandler::globalToLocal(const Position& global) const { + G4ThreeVector p = globalToLocalG4(G4ThreeVector(global.X(),global.Y(),global.Z())); + return Position(p.x(),p.y(),p.z()); +} + +/// Coordinate transformation to local coordinates +Position Geant4StepHandler::globalToLocal(const G4ThreeVector& global) const { + G4ThreeVector p = globalToLocalG4(global); + return Position(p.x(),p.y(),p.z()); +} + +/// Coordinate transformation to local coordinates +G4ThreeVector Geant4StepHandler::globalToLocalG4(double x, double y, double z) const { + return globalToLocalG4(G4ThreeVector(x,y,z)); +} + +/// Coordinate transformation to local coordinates +G4ThreeVector Geant4StepHandler::globalToLocalG4(const G4ThreeVector& global) const { + G4TouchableHandle t = step->GetPreStepPoint()->GetTouchableHandle(); + return t->GetHistory()->GetTopTransform().TransformPoint(global); +} diff --git a/DDG4/src/Geant4TrackingAction.cpp b/DDG4/src/Geant4TrackingAction.cpp index 1ea3c6d89..72767b028 100644 --- a/DDG4/src/Geant4TrackingAction.cpp +++ b/DDG4/src/Geant4TrackingAction.cpp @@ -144,5 +144,6 @@ bool Geant4TrackingAction::storeChild(Geant4TrackInformation* track_info) const } return true; } - return error(false, "storeChild: Geant4TrackInformation points to NULL!"); + error("storeChild: Geant4TrackInformation points to NULL!"); + return false; } diff --git a/DDRec/include/DDRec/Surface.h b/DDRec/include/DDRec/Surface.h index 2a3481c97..e273168bc 100644 --- a/DDRec/include/DDRec/Surface.h +++ b/DDRec/include/DDRec/Surface.h @@ -80,7 +80,7 @@ namespace DD4hep { VolSurfaceBase( SurfaceType typ, double thickness_inner ,double thickness_outer, Vector3D u_val ,Vector3D v_val , - Vector3D n ,Vector3D o, Geometry::Volume vol,int id ) : + Vector3D n ,Vector3D o, Geometry::Volume vol,int identifier ) : _type(typ ) , _u( u_val ) , _v( v_val ) , @@ -91,7 +91,7 @@ namespace DD4hep { _innerMat( MaterialData() ), _outerMat( MaterialData() ), _vol(vol) , - _id( id ), _refCount(0) { + _id( identifier ), _refCount(0) { } @@ -350,9 +350,9 @@ namespace DD4hep { /// standard c'tor with all necessary arguments - origin is (0,0,0) if not given. VolPlaneImpl( SurfaceType typ, double thickness_inner ,double thickness_outer, - Vector3D u_val ,Vector3D v_val ,Vector3D n_val , Vector3D o_val, Geometry::Volume vol, int id ) : + Vector3D u_val ,Vector3D v_val ,Vector3D n_val , Vector3D o_val, Geometry::Volume vol, int id_val ) : - VolSurfaceBase( typ, thickness_inner, thickness_outer, u_val,v_val, n_val, o_val, vol,id ) { + VolSurfaceBase( typ, thickness_inner, thickness_outer, u_val,v_val, n_val, o_val, vol, id_val ) { _type.setProperty( SurfaceType::Plane , true ) ; _type.setProperty( SurfaceType::Cylinder , false ) ; @@ -483,8 +483,8 @@ namespace DD4hep { class VolCylinder : public VolSurface{ public: - VolCylinder( Geometry::Volume vol, SurfaceType type, double thickness_inner ,double thickness_outer, Vector3D origin ) : - VolSurface( new VolCylinderImpl( vol, type, thickness_inner , thickness_outer, origin ) ) {} + VolCylinder( Geometry::Volume vol, SurfaceType type_val, double thickness_inner ,double thickness_outer, Vector3D origin_pos ) : + VolSurface( new VolCylinderImpl( vol, type_val, thickness_inner , thickness_outer, origin_pos ) ) {} } ; class VolCone : public VolSurface{ diff --git a/cmake/DD4hepBuild.cmake b/cmake/DD4hepBuild.cmake index 398b56991..ccdb8a6a9 100644 --- a/cmake/DD4hepBuild.cmake +++ b/cmake/DD4hepBuild.cmake @@ -1172,16 +1172,23 @@ function( dd4hep_add_dictionary dictionary ) dd4hep_debug ( "${tag} Unparsed:'${ARG_UNPARSED_ARGUMENTS}'" ) dd4hep_debug ( "${tag} Sources: '${CMAKE_CURRENT_SOURCE_DIR}'" ) # - add_custom_command(OUTPUT ${dictionary}.cxx ${dictionary}.h - COMMAND ${ROOTCINT_EXECUTABLE} -cint -f ${dictionary}.cxx - -s ${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary} -c -p ${ARG_OPTIONS} ${comp_defs} ${inc_dirs} ${headers} ${linkdefs} - DEPENDS ${headers} ${linkdefs} ) - # Install the binary to the destination directory if ( ${ROOT_VERSION_MAJOR} GREATER 5 ) - set_source_files_properties(${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary}_rdict.pcm PROPERTIES GENERATED TRUE ) + ## ${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary}_rdict.pcm + add_custom_command(OUTPUT ${dictionary}.cxx + COMMAND ${ROOTCINT_EXECUTABLE} -cint -f ${dictionary}.cxx + -s ${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary} -c -p ${ARG_OPTIONS} ${comp_defs} ${inc_dirs} ${headers} ${linkdefs} + DEPENDS ${headers} ${linkdefs} ) + # Install the binary to the destination directory + #set_source_files_properties(${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary}_rdict.pcm PROPERTIES GENERATED TRUE ) install(FILES ${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary}_rdict.pcm DESTINATION lib) + #set_source_files_properties( ${dictionary}.h ${dictionary}.cxx PROPERTIES GENERATED TRUE ) + else() + add_custom_command(OUTPUT ${dictionary}.h ${dictionary}.cxx + COMMAND ${ROOTCINT_EXECUTABLE} -cint -f ${dictionary}.cxx + -s ${CMAKE_CURRENT_BINARY_DIR}/../lib/${dictionary} -c -p ${ARG_OPTIONS} ${comp_defs} ${inc_dirs} ${headers} ${linkdefs} + DEPENDS ${headers} ${linkdefs} ) + #set_source_files_properties( ${dictionary}.h ${dictionary}.cxx PROPERTIES GENERATED TRUE ) endif() - set_source_files_properties( ${dictionary}.h ${dictionary}.cxx PROPERTIES GENERATED TRUE ) endif() endfunction() diff --git a/doc/release.notes b/doc/release.notes index 6b0e2321a..ea450104f 100644 --- a/doc/release.notes +++ b/doc/release.notes @@ -3,6 +3,18 @@ DD4hep ---- Release Notes ================================= +2015-10-13 M.Frank + DDG4 + - Remove explicit constructors for modular physics lists. + Use the native Geant4 provided G4PhysListFactory instead. + The physics list is instantiated as before with the "Extends" property. + - Implement angular ranges in the Geant4IsotropeGenerator for phi [0,2pi] and theta[0,pi] + - New sensitive detector Geant4TrackerWeightedAction + Attempt to properly process combined deposits in tracking detectors. + - Debug component Geant4SurfaceTest similar to the standalone program + test_surfaces, but to be appended as a DDG4 event action for event by event + tests of hits. + 2015-10-09 M.Frank DDG4 - Extend the functionality of the DDG4 plugins -- GitLab