Update FiberTubeCalorimeter example

DDG4/src/Geant4Output2ROOT.cpp

@@ -29,14 +29,16 @@
 #include <TFile.h>
 #include <TTree.h>
 #include <TBranch.h>
+#include <TSystem.h>
+
 
 using namespace dd4hep::sim;
 using namespace dd4hep;
 using namespace std;
 
 /// Standard constructor
 Geant4Output2ROOT::Geant4Output2ROOT(Geant4Context* ctxt, const string& nam)
-  : Geant4OutputAction(ctxt, nam), m_file(0), m_tree(0) {
+  : Geant4OutputAction(ctxt, nam), m_file(nullptr), m_tree(nullptr) {
   declareProperty("Section",              m_section = "EVENT");
   declareProperty("HandleMCTruth",        m_handleMCTruth = true);
   declareProperty("DisabledCollections",  m_disabledCollections);
@@ -94,11 +96,21 @@ void Geant4Output2ROOT::beginRun(const G4Run* run) {
   }
   if ( !m_file && !fname.empty() ) {
     TDirectory::TContext ctxt(TDirectory::CurrentDirectory());
-    m_file = TFile::Open(fname.c_str(), "RECREATE", "dd4hep Simulation data");
-    if (m_file->IsZombie()) {
+    if ( !gSystem->AccessPathName(fname.c_str()) )  {
+      gSystem->Unlink(fname.c_str());
+    }
+    std::unique_ptr<TFile> file(TFile::Open(fname.c_str(), "RECREATE", "dd4hep Simulation data"));
+    if ( !file )  {
+      file.reset(TFile::Open((fname+".1").c_str(), "RECREATE", "dd4hep Simulation data"));
+    }
+    if ( !file )  {
+      except("Failed to create ROOT output file:'%s'", fname.c_str());
+    }
+    if (file->IsZombie()) {
       detail::deletePtr (m_file);
       except("Failed to open ROOT output file:'%s'", fname.c_str());
     }
+    m_file = file.release();
     m_tree = section(m_section);
   }
   Geant4OutputAction::beginRun(run);

examples/ClientTests/compact/FiberTubeCalorimeter.xml

@@ -70,11 +70,12 @@
   </define>
 
   <display>
-    <vis name="AbsVis" alpha="1.0" r="0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
+    <vis name="AbsVis"   alpha="0.5" r="0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
     <vis name="CerenVis" alpha="0.5" r="0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
     <vis name="ScintVis" alpha="0.5" r="1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/>
     <vis name="phdetVis" alpha="1.0" r="1.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
     <vis name="towerVis" alpha="0.1" r="0.3" g="0.3" b="0.3" showDaughters="true" visible="true"/>
+    <vis name="layerVis" alpha="0.5" r="0.1" g="0.1" b="0.9" showDaughters="true" visible="true"/>
   </display>
 
   <readouts>
@@ -92,18 +93,18 @@
          zmin is where the dtector goes
     -->
     <dimensions numsides="DRFiberNsizeSCE"
-                thickness="DRFiberAbswidthSCE"
+                thickness="DRFiberAbswidthSCE+holeoverSCE"
                 z_length="DRFiberlengthSCE"
                 gap="gapSCE"
                 zmin="-world_side/2.+2*killthick+edgeoffset+DRcrystallength+EcalHcalgap"
                 z1="killthick"/>
       <structure>
         <core1  name="scintillator" rmax="DRFiberFibwidthSCE" rmin="0.0" material="Polystyrene" vis="ScintVis" sensitive="yes"/> 
-        <core2  name="quartz" rmax="DRFiberFibwidthSCE" rmin="0.0" material="Quartz" vis="CerenVis" sensitive="yes"/> 
-        <hole   name="hole" rmax="(DRFiberFibwidthSCE+holeoverSCE)" rmin="0.0" material="Air" vis="holeVis" sensitive="yes"/> 
-        <absorb material="Brass" vis="AbsVis" sensitive="yes"/>
-        <phdet1 name="phdet1" rmax="DRFiberFibwidthSCE" rmin="0.0" material="killMedia2" vis="phdetVis" sensitive="yes"/> 
-        <phdet2 name="phdet2" rmax="DRFiberFibwidthSCE" rmin="0.0" material="killMedia3" vis="phdetVis" sensitive="yes"/> 
+        <core2  name="quartz"       rmax="DRFiberFibwidthSCE" rmin="0.0" material="Quartz" vis="CerenVis" sensitive="yes"/> 
+        <hole   name="hole"         rmax="DRFiberFibwidthSCE+holeoverSCE" rmin="DRFiberFibwidthSCE" material="Air" vis="holeVis" sensitive="yes"/> 
+        <absorb material="Brass"    vis="AbsVis" sensitive="yes"/>
+        <phdet1 name="phdet1"       rmax="DRFiberFibwidthSCE" rmin="0.0" material="killMedia2" vis="phdetVis" sensitive="yes"/> 
+        <phdet2 name="phdet2"       rmax="DRFiberFibwidthSCE" rmin="0.0" material="killMedia3" vis="phdetVis" sensitive="yes"/> 
       </structure>
     </detector>
   </detectors>

examples/ClientTests/scripts/FiberTubeCalorimeter.py

@@ -0,0 +1,95 @@
+# ==========================================================================
+#  AIDA Detector description implementation
+# --------------------------------------------------------------------------
+# 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.
+#
+# ==========================================================================
+#
+#
+from __future__ import absolute_import, unicode_literals
+import os
+import time
+import DDG4
+from DDG4 import OutputLevel as Output
+from g4units import GeV, MeV, m
+#
+#
+"""
+
+   dd4hep simulation example setup using the python configuration
+
+   @author  M.Frank
+   @version 1.0
+
+"""
+
+
+def run():
+  args = DDG4.CommandLine()
+  kernel = DDG4.Kernel()
+  install_dir = os.environ['DD4hepExamplesINSTALL']
+  kernel.loadGeometry(str("file:" + install_dir + "/examples/ClientTests/compact/FiberTubeCalorimeter.xml"))
+
+  DDG4.importConstants(kernel.detectorDescription(), debug=False)
+  geant4 = DDG4.Geant4(kernel, tracker='Geant4TrackerCombineAction')
+  geant4.printDetectors()
+  # Configure UI
+  if args.macro:
+    ui = geant4.setupCshUI(macro=args.macro)
+  else:
+    ui = geant4.setupCshUI()
+  if args.batch:
+    ui.Commands = ['/run/beamOn ' + str(args.events), '/ddg4/UI/terminate']
+
+  # 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', 'FiberTubeCalorimeter_' + time.strftime('%Y-%m-%d_%H-%M'))
+
+  # Setup particle gun
+  gun = geant4.setupGun("Gun", particle='e+', energy=20 * GeV, multiplicity=1, position=(0.0, 0.0, -369.0))
+  gun.OutputLevel = generator_output_level
+
+  # And handle the simulation particles.
+  part = DDG4.GeneratorAction(kernel, "Geant4ParticleHandler/ParticleHandler")
+  kernel.generatorAction().adopt(part)
+  part.SaveProcesses = ['Decay']
+  part.MinimalKineticEnergy = 100 * MeV
+  part.OutputLevel = Output.INFO  # generator_output_level
+  part.enableUI()
+  user = DDG4.Action(kernel, "Geant4TCUserParticleHandler/UserParticleHandler")
+  user.TrackingVolume_Zmax = 3.0 * m
+  user.TrackingVolume_Rmax = 3.0 * m
+  user.enableUI()
+  part.adopt(user)
+
+  geant4.setupCalorimeter('FiberTubeCalorimeter')
+
+  # Now build the physics list:
+  phys = geant4.setupPhysics('QGSP_BERT')
+  phys.dump()
+  geant4.execute()
+
+
+if __name__ == "__main__":
+  run()

examples/ClientTests/src/FiberTubeCalorimeter_geo.cpp

@@ -68,27 +68,17 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
   // these refer to different fields in the xml file for this detector
   xml_comp_t fX_struct( x_det.child( _Unicode(structure) ) );
   xml_comp_t fX_absorb( fX_struct.child( _Unicode(absorb) ) );
-  xml_comp_t fX_core1( fX_struct.child( _Unicode(core1) ) );
-  xml_comp_t fX_core2( fX_struct.child( _Unicode(core2) ) );
-  xml_comp_t fX_hole( fX_struct.child( _Unicode(hole) ) );
+  xml_comp_t fX_core1(  fX_struct.child( _Unicode(core1) ) );
+  xml_comp_t fX_core2(  fX_struct.child( _Unicode(core2) ) );
+  xml_comp_t fX_hole(   fX_struct.child( _Unicode(hole) ) );
   xml_comp_t fX_phdet1( fX_struct.child( _Unicode(phdet1) ) );
   xml_comp_t fX_phdet2( fX_struct.child( _Unicode(phdet2) ) );
-
-  // detector element for entire detector.  
-  DetElement    sdet      (det_name, det_id);
-  Volume        motherVol = description.pickMotherVolume(sdet);
-  Box           env_box   ((2*Ncount+1)*(hthick+agap+tol),(2*Ncount+1)*(hthick+agap+tol), (hzlength+hzph+tol));
-  Volume        envelopeVol  (det_name, env_box, air);
-  envelopeVol.setAttributes(description,x_det.regionStr(),x_det.limitsStr(),x_det.visStr());
 
   Material mat;
   Transform3D trafo;
   PlacedVolume pv;
   Solid sol;
-
-  pv = motherVol.placeVolume(envelopeVol, Position(0.,0.,azmin+hzlength+hzph+tol));
-  pv.addPhysVolID("system", det_id);
-  sdet.setPlacement(pv);  // associate the placed volume to the detector element
+
   sens.setType("calorimeter");
 
   // scint fiber
@@ -182,24 +172,40 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
   }
 
   // setup the volumes with the shapes and properties in one horixontal layer
-  Volume tube_row_vol("layer", Box(hthick,hthick,hzlength+hzph), air);
+  double dx = 2*(Ncount + Ncount+1)/2e0 * (hthick+agap) + tol;
+  double dy = hthick + tol;
+  double dz = hzlength+hzph + tol;
+  Box    tube_row_box(dx, dy, dz);
+  Volume tube_row_vol("layer", tube_row_box, air);
   tube_row_vol.setVisAttributes(description, x_det.visStr());
   tube_row_vol.setSensitiveDetector(sens);
-
+  cout << tube_row_vol.name()
+       << " dx: " << tube_row_box.x()
+       << " dy: " << tube_row_box.y()
+       << " dz: " << tube_row_box.z() << endl;
+  tube_row_vol.setVisAttributes(description, "layerVis");
+
   for (int ijk=-Ncount; ijk<Ncount+1; ijk++) {
-    double mod_x_off = (ijk)*2*(hthick+agap);
-    Transform3D tr(RotationZYX(0.,0.,0.), Position(mod_x_off,0.,0.));
+    double mod_x_off = (ijk) * 2 * (hthick + agap);
     int towernum = Ncount + ijk + 1;
-    pv = tube_row_vol.placeVolume((towernum%2 == 0) ? quartz_abs_vol : scint_abs_vol, tr);
+    pv = tube_row_vol.placeVolume((towernum%2 == 0) ? quartz_abs_vol : scint_abs_vol, Position(mod_x_off,0.,0.));
     pv.addPhysVolID("tube", towernum);
-    cout << "Placing row  " << setw(5) << right << ijk
-         << " x-offset: "   << setw(7) << right<< mod_x_off
+    //Box bounding_box = pv.volume().solid().GetBoundingBox();
+    cout << "Placing row  "    << setw(5) << right << ijk
+         << " x-offset: "      << setw(7) << right << mod_x_off
          << " volume of type " << pv.volume().name()
          << endl;
   }
+
+  dy = 2*(Ncount + Ncount+1)/2e0 * (hthick+agap) + tol;
+  DetElement    sdet      (det_name, det_id);
+  Box           env_box   (dx+tol, dy+tol, dz+tol);
+  Volume        envelopeVol  (det_name, env_box, air);
+  envelopeVol.setAttributes(description, x_det.regionStr(), x_det.limitsStr(), x_det.visStr());
+
   // Now stack multiple horizontal layers to form the final box
   for (int ijk=-Ncount; ijk<Ncount+1; ijk++) {
-    double mod_y_off = (ijk)*2*(hthick+agap);
+    double mod_y_off = (ijk) * 2 * (tube_row_box.y() + agap);
     Transform3D tr(RotationZYX(0.,0.,0.), Position(0.,mod_y_off,0.));
     pv = envelopeVol.placeVolume(tube_row_vol, tr);
     pv.addPhysVolID("layer", Ncount+ijk+1);
@@ -212,6 +218,12 @@ static Ref_t create_detector(Detector& description, xml_h e, SensitiveDetector s
          << " volume of type " << pv.volume().name()
          << endl;
   }
+
+  // detector element for entire detector.  
+  Volume        motherVol = description.pickMotherVolume(sdet);
+  pv = motherVol.placeVolume(envelopeVol, Position(0.,0.,azmin+hzlength+hzph+tol));
+  pv.addPhysVolID("system", det_id);
+  sdet.setPlacement(pv);  // associate the placed volume to the detector element
   return sdet;
 }