diff --git a/SimG4Components/CMakeLists.txt b/SimG4Components/CMakeLists.txt index 72f0cae..80443fd 100644 --- a/SimG4Components/CMakeLists.txt +++ b/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) diff --git a/SimG4Components/src/SimG4OpticalPhysicsList.cpp b/SimG4Components/src/SimG4OpticalPhysicsList.cpp new file mode 100644 index 0000000..015e862 --- /dev/null +++ b/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(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; +} + diff --git a/SimG4Components/src/SimG4OpticalPhysicsList.h b/SimG4Components/src/SimG4OpticalPhysicsList.h new file mode 100644 index 0000000..24f1f23 --- /dev/null +++ b/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 SetCerenkov{this, "SetCerenkov", true, "Bool variable that enables Cerenkov process. Default true."}; + Gaudi::Property SetScintillation{this, "SetScintillation", true, "Bool variable that enables Scintillation process. Default true."}; + Gaudi::Property SetTransitionRadiation{this, "SetTransitionRadiation", false, "Bool variable that enables transition_radiation process. Default false."}; +private: + /// Handle for the full physics list tool + ToolHandle m_physicsListTool{"SimG4FtfpBert", this, true}; +}; + +#endif diff --git a/SimG4Components/tests/options/optical_physics_test.py b/SimG4Components/tests/options/optical_physics_test.py new file mode 100644 index 0000000..91b944b --- /dev/null +++ b/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, + )