diff --git a/include/SimCore/GammaPhysics.h b/include/SimCore/GammaPhysics.h index c2dcd65..0b112af 100644 --- a/include/SimCore/GammaPhysics.h +++ b/include/SimCore/GammaPhysics.h @@ -3,10 +3,12 @@ * @brief Class used to enhanced the gamma physics list. * @author Jeremy McCormick, SLAC National Accelerator Laboratory * @author Omar Moreno, SLAC National Accelerator Laboratory + * @author Einar Elen, Lund University + * @author Tom Eichlersmith, University of Minnesota */ #ifndef SIMCORE_GAMMAPHYSICS_H_ -#define SIMCORE_GAMMAPHYSICS_H 1_ +#define SIMCORE_GAMMAPHYSICS_H_ //------------// // Geant4 // @@ -16,6 +18,7 @@ #include "G4VPhysicsConstructor.hh" #include "G4VProcess.hh" #include "SimCore/PhotoNuclearModel.h" + namespace simcore { /** @@ -38,11 +41,7 @@ class GammaPhysics : public G4VPhysicsConstructor { * @param parameters The python configuration */ GammaPhysics(const G4String& name, - const framework::config::Parameters& parameters) - : G4VPhysicsConstructor(name), - modelParameters{parameters.getParameter( - "photonuclear_model")} {} - GammaPhysics(const G4String& name = "GammaPhysics"); + const framework::config::Parameters& parameters); /** * Class destructor. @@ -50,34 +49,31 @@ class GammaPhysics : public G4VPhysicsConstructor { virtual ~GammaPhysics() = default; /** - * Construct particles (no-op). - */ - void ConstructParticle() {} - - /** - * Construct the process (gamma to muons). + * Construct particles + * + * We don't do anything here since we are just attaching/updating + * the photon physics. */ - void ConstructProcess(); + void ConstructParticle() final; /** - * Search the list for the process "photoNuclear". When found, - * change the calling order so photonNuclear is called before - * EM processes. The biasing operator needs the photonNuclear - * process to be called first because the cross-section is - * needed to bias down other process. + * We do two things for this call back during initialization. + * + * 1. We add the muon-conversion process to the photon's process + * table, enabling it to be one of the processes that can happen. + * 2. We move the photonNuclear process to be the first process in + * the photon's ordering so that the bias operator has a chance + * to lower its cross-section before any EM process is called + * (if need be). */ - void SetPhotonNuclearAsFirstProcess() const; + void ConstructProcess() final; private: - /* - * Returns the process manager object for the G4Gamma class from the list of - * Geant4 particles. - */ - G4ProcessManager* GetGammaProcessManager() const; /** * The gamma to muons process. */ G4GammaConversionToMuons gammaConvProcess; + /** * Parameters from the configuration to pass along to the photonuclear model. */ diff --git a/src/SimCore/GammaPhysics.cxx b/src/SimCore/GammaPhysics.cxx index 0f26d16..ed65361 100644 --- a/src/SimCore/GammaPhysics.cxx +++ b/src/SimCore/GammaPhysics.cxx @@ -9,60 +9,54 @@ namespace simcore { -// needed for GEANT4 10.3.0 and later -// -// TODO: I (Einar) don't think this really is necessary if we change -// GetGammaProcessManager below -#ifndef aParticleIterator -#define aParticleIterator \ - ((subInstanceManager.offset[g4vpcInstanceID])._aParticleIterator) -#endif +GammaPhysics::GammaPhysics(const G4String& name, + const framework::config::Parameters& parameters) + : G4VPhysicsConstructor(name), + modelParameters{parameters.getParameter( + "photonuclear_model")} {} -G4ProcessManager* GammaPhysics::GetGammaProcessManager() const { - /** - * Not entirely clear that there isn't a simpler way to do this, but to keep - * this function extraction to a pure refactoring, I'm leaving it here for - * now. /Einar - */ - aParticleIterator->reset(); - - // Loop through all of the particles and find the "gamma". - while ((*aParticleIterator)()) { - auto particle{aParticleIterator->value()}; - auto particleName{particle->GetParticleName()}; - if (particleName == "gamma") { - return particle->GetProcessManager(); - } - } - EXCEPTION_RAISE("GammaPhysics", - "Was unable to access the process manager for photons, " - "something is very wrong!"); -} +void GammaPhysics::ConstructParticle() {} -void GammaPhysics::SetPhotonNuclearAsFirstProcess() const { - auto processManager{GetGammaProcessManager()}; - // Get the process list associated with the gamma. - const auto processes{processManager->GetProcessList()}; - - for (int i{0}; i < processes->size(); ++i) { - const auto process{(*processes)[i]}; - const auto processName{process->GetProcessName()}; - if (processName == "photonNuclear") { - processManager->SetProcessOrderingToFirst( - process, G4ProcessVectorDoItIndex::idxAll); - } - } -} void GammaPhysics::ConstructProcess() { - G4ProcessManager* processManager = GetGammaProcessManager(); - + G4ProcessManager* processManager = G4Gamma::Gamma()->GetProcessManager(); + if (processManager == nullptr) { + EXCEPTION_RAISE("GammaPhysics", + "Was unable to access the process manager for photons, " + "something is very wrong!"); + } + // configure our PN model based on runtime parameters auto pn = PhotoNuclearModel::Factory::get().make( modelParameters.getParameter("class_name"), modelParameters.getParameter("instance_name"), modelParameters); pn->removeExistingModel(processManager); pn->ConstructGammaProcess(processManager); - SetPhotonNuclearAsFirstProcess(); + /** + * Put the PN process first in the ordering in case PN biasing is happening. + * + * Process ordering is a complicated concept and unfortunately + * its affect on biasing is poorly documented. What has been said is + * that some processes _need_ to be called _last_[1]. In addition, + * the practical experience of working on defining a custom + * process for + * [G4DarkBreM](https://github.com/LDMX-Software/G4DarkBreM) + * showed that sometimes Geant4 does + * not get through the full list of processes but it always starts + * at the beginning. For these reasons, we put the PN process first + * in the ordering so that we can insure it is always check by Geant4 + * before continuing. + * + * [1] + * https://indico.cern.ch/event/58317/contributions/2047449/attachments/992300/1411062/PartDef-ProcMan-AS.pdf + */ + const auto processes{processManager->GetProcessList()}; + for (int i{0}; i < processes->size(); i++) { + const auto process{(*processes)[i]}; + if (process->GetProcessName() == "photonNuclear") { + processManager->SetProcessOrderingToFirst( + process, G4ProcessVectorDoItIndex::idxAll); + } + } // Add the gamma -> mumu to the physics list. processManager->AddDiscreteProcess(&gammaConvProcess); }