Merge pull request #44 from atolosadelgado/optical_physicslist_tool

Optical physicslist tool

SimG4Components/CMakeLists.txt

@@ -31,6 +31,12 @@ add_test(NAME MagFieldFromDD4hep
          WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
          COMMAND bash -c "source k4simgeant4env.sh;  k4run ${CMAKE_CURRENT_LIST_DIR}/tests/options/magFieldTool.py"
 )
+add_test(NAME OpticalPhysicsTest
+         WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
+         COMMAND bash -c "source k4simgeant4env.sh;  k4run ${CMAKE_CURRENT_LIST_DIR}/tests/options/optical_physics_test.py"
+)
+SET_TESTS_PROPERTIES( OpticalPhysicsTest PROPERTIES FAIL_REGULAR_EXPRESSION " Exception;EXCEPTION;ERROR;Error" )
+SET_TESTS_PROPERTIES( OpticalPhysicsTest PROPERTIES PASS_REGULAR_EXPRESSION " Cerenkov process active:               1;  Scintillation process active:          1;  Rayleigh process active:               1;  Absorption process active:             1")
 #gaudi_add_test(GeantFullSimGdml
 #               WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
 #               FRAMEWORK tests/options/geant_fullsim_gdml.py)

SimG4Components/src/SimG4OpticalPhysicsList.cpp

@@ -0,0 +1,47 @@
+#include "SimG4OpticalPhysicsList.h"
+
+// Geant4
+#include "G4VModularPhysicsList.hh"
+#include "G4OpticalPhysics.hh"
+#include "G4OpticalParameters.hh"
+
+DECLARE_COMPONENT(SimG4OpticalPhysicsList)
+
+SimG4OpticalPhysicsList::SimG4OpticalPhysicsList(const std::string& aType, const std::string& aName, const IInterface* aParent)
+: AlgTool(aType, aName, aParent) {
+  declareInterface<ISimG4PhysicsList>(this);
+  declareProperty("fullphysics", m_physicsListTool, "Handle for the full physics list tool");
+}
+
+SimG4OpticalPhysicsList::~SimG4OpticalPhysicsList() {}
+
+StatusCode SimG4OpticalPhysicsList::initialize() {
+  if (AlgTool::initialize().isFailure())
+    return StatusCode::FAILURE;
+
+  m_physicsListTool.retrieve().ignore();
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode SimG4OpticalPhysicsList::finalize() { return AlgTool::finalize(); }
+
+G4VModularPhysicsList* SimG4OpticalPhysicsList::physicsList() {
+  // ownership passed to SimG4Svc which will register it in G4RunManager. To be deleted in ~G4RunManager()
+  G4VModularPhysicsList* physicsList = m_physicsListTool->physicsList();
+
+  G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics(); // deleted in ~G4VModularPhysicsList()
+  opticalPhysics->SetVerboseLevel (2);
+  physicsList->RegisterPhysics(opticalPhysics);
+
+  auto* opticalParams = G4OpticalParameters::Instance();
+  opticalParams->SetBoundaryInvokeSD(true);
+  opticalParams->SetProcessActivation("Cerenkov",SetCerenkov);
+  opticalParams->SetProcessActivation("Scintillation",SetScintillation);
+  opticalParams->SetProcessActivation("TransitionRadiation",SetTransitionRadiation);
+  opticalParams->SetCerenkovTrackSecondariesFirst(true);
+  opticalParams->SetScintTrackSecondariesFirst(true);
+
+  return physicsList;
+}
+

SimG4Components/src/SimG4OpticalPhysicsList.h

@@ -0,0 +1,39 @@
+#ifndef SimG4OpticalPhysicsList_h
+#define SimG4OpticalPhysicsList_h 1
+
+// Gaudi
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+// FCCSW
+#include "k4Interface/ISimG4PhysicsList.h"
+
+/** @class SimG4OpticalPhysicsList SimG4Components/src/SimG4OpticalPhysicsList.h SimG4OpticalPhysicsList.h
+ *
+ *  FTFP_BERT physics list + Optical photons physics lists tool.
+ *  
+ *  When instantiating the tool three booleans can be passed as arguments:
+ *  --SetCerenkov, default true, to enable Cerenkov process
+ *  --SetScintillation, default true, to enable Scintillation process
+ *  --SetTransitionRadiation, default true, to enable Transition Radiation process
+ *
+ *  @author Alvaro Tolosa-Delgado
+ */
+
+class SimG4OpticalPhysicsList : public AlgTool, virtual public ISimG4PhysicsList {
+public:
+  explicit SimG4OpticalPhysicsList(const std::string& aType, const std::string& aName, const IInterface* aParent);
+  virtual ~SimG4OpticalPhysicsList();
+
+  virtual StatusCode initialize();
+  virtual StatusCode finalize();
+  virtual G4VModularPhysicsList* physicsList();
+  Gaudi::Property<bool> SetCerenkov{this, "SetCerenkov", true, "Bool variable that enables Cerenkov process. Default true."};
+  Gaudi::Property<bool> SetScintillation{this, "SetScintillation", true, "Bool variable that enables Scintillation process. Default true."};
+  Gaudi::Property<bool> SetTransitionRadiation{this, "SetTransitionRadiation", false, "Bool variable that enables transition_radiation process. Default false."};
+private:
+  /// Handle for the full physics list tool
+  ToolHandle<ISimG4PhysicsList> m_physicsListTool{"SimG4FtfpBert", this, true};
+};
+
+#endif

SimG4Components/tests/options/optical_physics_test.py

@@ -0,0 +1,109 @@
+# This script launch a simulation of 1 event and a given detector
+# The physics invoked is the FTFpBert baseline + Optical physics
+
+import os
+import copy
+
+from GaudiKernel.SystemOfUnits import MeV, GeV, tesla
+
+# Input for simulations (momentum is expected in GeV!)
+# Parameters for the particle gun simulations, dummy if use_pythia is set to True
+# theta from 80 to 100 degrees corresponds to -0.17 < eta < 0.17 
+momentum = 1 # in GeV
+#thetaMin = 90.25 # degrees
+#thetaMax = 90.25 # degrees
+thetaMin = 50 # degrees
+thetaMax = 130 # degrees
+pdgCode = 11 # 11 electron, 13 muon, 22 photon, 111 pi0, 211 pi+
+
+from Gaudi.Configuration import *
+
+from Configurables import FCCDataSvc
+podioevent  = FCCDataSvc("EventDataSvc")
+
+################## Particle gun setup
+_pi = 3.14159
+
+from Configurables import GenAlg
+genAlg = GenAlg()
+
+from Configurables import  MomentumRangeParticleGun
+pgun = MomentumRangeParticleGun("ParticleGun_Electron")
+pgun.PdgCodes = [pdgCode]
+pgun.MomentumMin = momentum * GeV
+pgun.MomentumMax = momentum * GeV
+pgun.PhiMin = 0
+#pgun.PhiMax = 0
+pgun.PhiMax = 2 * _pi
+pgun.ThetaMin = thetaMin * _pi / 180.
+pgun.ThetaMax = thetaMax * _pi / 180.
+genAlg.SignalProvider = pgun
+genAlg.hepmc.Path = "hepmc"
+
+from Configurables import HepMCToEDMConverter
+hepmc_converter = HepMCToEDMConverter()
+hepmc_converter.hepmc.Path="hepmc"
+genParticlesOutputName = "genParticles"
+hepmc_converter.GenParticles.Path = genParticlesOutputName
+hepmc_converter.hepmcStatusList = []
+
+################## Simulation setup
+# Detector geometry
+from Configurables import GeoSvc
+geoservice = GeoSvc("GeoSvc")
+# if FCC_DETECTORS is empty, this should use relative path to working directory
+path_to_detector = os.environ.get("FCCDETECTORS", "")
+print(path_to_detector)
+detectors_to_use=[
+                    'Detector/DetFCCeeIDEA-LAr/compact/FCCee_DectMaster.xml',
+                  ]
+# prefix all xmls with path_to_detector
+geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use]
+geoservice.OutputLevel = INFO
+
+
+# Geant4 service
+# Configures the Geant simulation: geometry, physics list and user actions
+
+from Configurables import SimG4FullSimActions, SimG4Alg, SimG4PrimariesFromEdmTool, SimG4SaveParticleHistory
+actions = SimG4FullSimActions()
+# Uncomment if history from Geant4 decays is needed (e.g. to get the photons from pi0) and set actions=actions in SimG4Svc + uncomment saveHistTool in SimG4Alg
+#actions.enableHistory=True
+#actions.energyCut = 0.2 * GeV
+#saveHistTool = SimG4SaveParticleHistory("saveHistory")
+
+from Configurables import SimG4OpticalPhysicsList
+physicslisttool = SimG4OpticalPhysicsList("Physics", fullphysics="SimG4FtfpBert")
+
+from Configurables import SimG4Svc
+geantservice = SimG4Svc("SimG4Svc", detector='SimG4DD4hepDetector', physicslist=physicslisttool, actions=actions)
+
+# Fixed seed to have reproducible results, change it for each job if you split one production into several jobs
+# Mind that if you leave Gaudi handle random seed and some job start within the same second (very likely) you will have duplicates
+geantservice.randomNumbersFromGaudi = False
+geantservice.seedValue = 4242
+
+# Range cut
+geantservice.g4PreInitCommands += ["/tracking/verbose 0"]
+# next, create the G4 algorithm, giving the list of names of tools ("XX/YY")
+from Configurables import SimG4PrimariesFromEdmTool
+particle_converter = SimG4PrimariesFromEdmTool("EdmConverter")
+particle_converter.GenParticles.Path = genParticlesOutputName
+
+from Configurables import SimG4Alg
+geantsim = SimG4Alg("SimG4Alg",
+                       eventProvider=particle_converter,
+                       OutputLevel=INFO)
+
+from Configurables import ApplicationMgr
+ApplicationMgr(
+    TopAlg = [
+              genAlg,
+              hepmc_converter,
+              geantsim,
+              ],
+    EvtSel = 'NONE',
+    EvtMax   = 1,
+    ExtSvc = [geoservice, geantservice],
+    StopOnSignal = True,
+ )