Allow for smartless configuration of volumes in Geant4

DDG4/examples/SiDSim.py

@@ -79,8 +79,9 @@ def run():
   if args.batch:
     cmds.append('/ddg4/UI/terminate')
 
-  ui.Commands = cmds
-
+  if len(cmds) > 0:
+    ui.Commands = cmds
+
   logger.info("#  Configure G4 magnetic field tracking")
   geant4.setupTrackingField()
 

DDG4/src/Geant4Converter.cpp

@@ -785,7 +785,7 @@ void* Geant4Converter::handleVolume(const std::string& name, const TGeoVolume* v
     }
 
     G4LogicalVolume* g4vol = nullptr;
-    if ( _v.hasProperties() && !_v.getProperty(GEANT4_TAG_PLUGIN,"").empty() )   {
+    if( _v.hasProperties() && !_v.getProperty(GEANT4_TAG_PLUGIN,"").empty() )   {
       Detector* det = const_cast<Detector*>(&m_detDesc); 
       std::string plugin = _v.getProperty(GEANT4_TAG_PLUGIN,"");
       g4vol = PluginService::Create<G4LogicalVolume*>(plugin, det, _v, g4solid, g4medium);
@@ -796,17 +796,19 @@ void* Geant4Converter::handleVolume(const std::string& name, const TGeoVolume* v
     else  {
       g4vol = new G4LogicalVolume(g4solid, g4medium, vnam, nullptr, nullptr, nullptr);
     }
+    PrintLevel plevel = (debugVolumes||debugRegions||debugLimits) ? ALWAYS : outputLevel;
     /// Set smartless optimization
     unsigned char smart_less_value = _v.smartlessValue();
-    if ( smart_less_value != Volume::NO_SMARTLESS_OPTIMIZATION )  {
-      g4vol->SetSmartless(smart_less_value);
+    if( smart_less_value != Volume::NO_SMARTLESS_OPTIMIZATION )  {
+      printout(ALWAYS, "Geant4Converter", "++ Volume %s Set Smartless value to %d",
+               vnam, int(smart_less_value));
+      g4vol->SetSmartless( smart_less_value );
     }
     /// Assign limits if necessary
-    if ( g4limits )   {
+    if( g4limits )   {
       g4vol->SetUserLimits(g4limits);
     }
-    if ( g4region )   {
-      PrintLevel plevel = (debugVolumes||debugRegions||debugLimits) ? ALWAYS : outputLevel;
+    if( g4region )   {
       printout(plevel, "Geant4Converter", "++ Volume     + Apply REGION settings: %-24s to volume %s.",
                reg.name(), vnam);
       // Handle the region settings for the world volume seperately.

examples/ClientTests/CMakeLists.txt

@@ -641,7 +641,8 @@ if (DD4HEP_USE_GEANT4)
       REGEX_PASS "ResourcesAfterConstruction       ConstructSD:   VmRSS"
       REGEX_FAIL "Error;ERROR; Exception"
   )
-
+  #
+  #
   dd4hep_add_test_reg(ClientTests_ddsim_setup_BoxOfStraws
       COMMAND    "${CMAKE_INSTALL_PREFIX}/bin/run_test_ClientTests.sh"
       EXEC_ARGS  ddsim
@@ -650,9 +651,18 @@ if (DD4HEP_USE_GEANT4)
                  --enableGun
                  --numberOfEvents 1
                  --outputFile regex.slcio
-
       REGEX_PASS "BoxOfStrawsDet Handled 1 nodes with 1 sensitive volume type"
       REGEX_FAIL "Error;ERROR; Exception"
   )
-
+  #
+  #  Test Changing Geant4 voxelization
+  dd4hep_add_test_reg( ClientTests_geant4_change_voxelization
+    COMMAND    "${CMAKE_INSTALL_PREFIX}/bin/run_test_ClientTests.sh"
+    EXEC_ARGS  ${Python_EXECUTABLE} ${ClientTestsEx_INSTALL}/scripts/DriftChamber.py
+    	       -verbose 2 -events 1
+    REGEX_PASS "100.00         0...         0...            0k DriftChamber_vol"
+    REGEX_FAIL "Exception;EXCEPTION;ERROR;Error;FATAL"
+  )
+  #
+  #
 endif(DD4HEP_USE_GEANT4)

examples/ClientTests/compact/DriftChamber.xml

@@ -22,8 +22,8 @@
     <constant name="world_x" value="world_size"/>
     <constant name="world_y" value="world_size"/>
     <constant name="world_z" value="world_size"/>
-    <constant name="DriftChamber_rmin" value="30*cm"/>
-    <constant name="DriftChamber_rmax" value="100*cm"/>
+    <constant name="DriftChamber_rmin" value="20*cm"/>
+    <constant name="DriftChamber_rmax" value="800*cm"/>
     <constant name="DriftChamber_zmax" value="100*cm"/>
   </define>
 
@@ -36,9 +36,11 @@
 
   <detectors>
     <detector id="1" name="DriftChamber" type="DD4hep_DriftChamber" readout="DriftChamberHits" vis="VisibleBlue">
-      <comment>A barrel hadronic calorimeter inspired on the ATLAS Tile hadronic calorimeter</comment>
-      <dimensions rmin="DriftChamber_rmin" rmax="DriftChamber_rmax" dz="DriftChamber_zmax" layers="20" />
-      <wire       count="60" thickness="3*mm" material="Cu" vis="VisibleRed"/>
+<!--  Geant4 default:
+      <dimensions rmin="DriftChamber_rmin" rmax="DriftChamber_rmax" dz="DriftChamber_zmax" layers="200" option="2"/>
+-->
+      <dimensions rmin="DriftChamber_rmin" rmax="DriftChamber_rmax" dz="DriftChamber_zmax" layers="200" option="0"/>
+      <wire       count="360" thickness="0.03*mm" material="Cu" vis="VisibleRed"/>
     </detector>
   </detectors>
 

examples/ClientTests/scripts/DriftChamber.py

@@ -14,8 +14,11 @@
 import os
 import time
 import DDG4
+import logging
 from DDG4 import OutputLevel as Output
 from g4units import GeV, MeV, m
+logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.INFO)
+logger = logging.getLogger(__name__)
 #
 #
 """
@@ -47,62 +50,91 @@ def run():
   cmds = []
   if args.verbose:
     cmds.append('/run/verbose ' + str(args.verbose))
+
   if args.events:
     cmds.append('/run/beamOn ' + str(args.events))
     cmds.append('/ddg4/UI/terminate')
 
   if len(cmds) > 0:
     ui.Commands = cmds
 
-  # Configure field
-  geant4.setupTrackingField(prt=True)
-  # 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)
-
-  generator_output_level = Output.INFO
-
-  # Configure G4 geometry setup
-  seq, act = geant4.addDetectorConstruction("Geant4DetectorGeometryConstruction/ConstructGeo")
-  act.DebugMaterials = True
-  act.DebugElements = False
-  act.DebugVolumes = True
-  act.DebugShapes = True
-  seq, act = geant4.addDetectorConstruction("Geant4DetectorSensitivesConstruction/ConstructSD")
-
-  # Configure I/O
-  geant4.setupROOTOutput('RootOutput', 'DriftChamber_' + time.strftime('%Y-%m-%d_%H-%M'))
 
-  # Setup particle gun
-  gun = geant4.setupGun("Gun", particle='e+', energy=20 * GeV, multiplicity=1)
-  gun.OutputLevel = generator_output_level
+  logger.info("#  Configure G4 magnetic field tracking")
+  geant4.setupTrackingField()
 
-  # And handle the simulation particles.
+  logger.info("#  Setup random generator")
+  rndm = DDG4.Action(kernel, 'Geant4Random/Random')
+  rndm.Seed = 987654321
+  if args.seed_time:
+    rndm.Seed = int(time.time())
+  rndm.initialize()
+
+  logger.info("""
+  Configure I/O
+  """)
+  geant4.setupROOTOutput('RootOutput', 'DriftChamber_' + time.strftime('%Y-%m-%d_%H-%M'))
+
+  gen = DDG4.GeneratorAction(kernel, "Geant4GeneratorActionInit/GenerationInit")
+  kernel.generatorAction().adopt(gen)
+
+  logger.info("""
+  Generation of isotrope tracks of a given multiplicity with overlay:
+  """)
+  logger.info("#  First particle generator: pi+")
+  gen = DDG4.GeneratorAction(kernel, "Geant4IsotropeGenerator/IsotropPi+")
+  gen.Mask = 1
+  gen.Particle = 'pi+'
+  gen.Energy = 100 * GeV
+  gen.Multiplicity = 2
+  gen.Distribution = 'cos(theta)'
+  kernel.generatorAction().adopt(gen)
+  logger.info("#  Install vertex smearing for this interaction")
+
+  logger.info("#  Merge all existing interaction records")
+  gen = DDG4.GeneratorAction(kernel, "Geant4InteractionMerger/InteractionMerger")
+  gen.OutputLevel = 4  # generator_output_level
+  gen.enableUI()
+  kernel.generatorAction().adopt(gen)
+  #
+  logger.info("#  Finally generate Geant4 primaries")
+  gen = DDG4.GeneratorAction(kernel, "Geant4PrimaryHandler/PrimaryHandler")
+  gen.OutputLevel = 4  # generator_output_level
+  gen.enableUI()
+  kernel.generatorAction().adopt(gen)
+  #
+  logger.info("#  ....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.INFO  # generator_output_level
+  part.OutputLevel = 5  # generator_output_level
   part.enableUI()
   user = DDG4.Action(kernel, "Geant4TCUserParticleHandler/UserParticleHandler")
-  user.TrackingVolume_Zmax = 3.0 * m
-  user.TrackingVolume_Rmax = 3.0 * m
+  user.TrackingVolume_Zmax = 2*m
+  user.TrackingVolume_Rmax = 2*m
   user.enableUI()
   part.adopt(user)
-
-  geant4.setupTracker('DriftChamber')
-
-  # Now build the physics list:
+  #
+  seq, act = geant4.setupTracker('DriftChamber')
+  #
+  logger.info("#  Now build the physics list:")
   phys = geant4.setupPhysics('QGSP_BERT')
-  ph = DDG4.PhysicsList(kernel, str('Geant4PhysicsList/Myphysics'))
-  ph.addParticleConstructor(str('G4BosonConstructor'))
-  ph.enableUI()
-  phys.adopt(ph)
+  ph = geant4.addPhysics(str('Geant4PhysicsList/Myphysics'))
+  ph.addPhysicsConstructor(str('G4StepLimiterPhysics'))
+  #
+  # Add special particle types from specialized physics constructor
+  part = geant4.addPhysics(str('Geant4ExtraParticles/ExtraParticles'))
+  part.pdgfile = os.path.join(install_dir, 'examples/DDG4/examples/particle.tbl')
+  #
   phys.dump()
+  #
+  kernel.configure()
+  kernel.initialize()
 
-  geant4.execute()
+  # DDG4.setPrintLevel(Output.DEBUG)
+  kernel.run()
+  kernel.terminate()
 
 
 if __name__ == "__main__":

examples/ClientTests/src/DriftChamber_geo.cpp

@@ -54,12 +54,17 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
   Volume        wire_vol("Wire_vol", wire_cyl, wire_mat);
   PlacedVolume  pv;
 
-  sdet_vol.setSmartlessValue(2);
+  /// The Geant4 voxelization change must be applied to the parent volume
+  if ( x_dim.hasAttr(_U(option)) )  {
+    int value = x_dim.attr<int>(_U(option));
+    printout(ALWAYS, "DriftChamber", "+++ Setting smartlessValue to %d for %s",
+             value, sdet_vol.name());
+    sdet_vol.setSmartlessValue(value);
+  }
+  /// Set volume attributes
   sdet_vol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
-
-  wire_vol.setSmartlessValue(2);
   wire_vol.setVisAttributes(description.visAttributes(x_wire.visStr()));
-
+  /// Place the wires in layers around the origin
   for( std::size_t l=0; l<layer_cnt; ++l )  {
     double radius = cyl_rmin + (0.5+double(l)) * layer_thickness;
     double phi_s  = phi_start + delta_phi * 0.5 * ((l%2) - 1.0);