get process manager from particle singleton

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,46 +41,39 @@ class GammaPhysics : public G4VPhysicsConstructor {
    * @param parameters The python configuration
    */
   GammaPhysics(const G4String& name,
-               const framework::config::Parameters& parameters)
-      : G4VPhysicsConstructor(name),
-        modelParameters{parameters.getParameter<framework::config::Parameters>(
-            "photonuclear_model")} {}
-  GammaPhysics(const G4String& name = "GammaPhysics");
+               const framework::config::Parameters& parameters);
 
   /**
    * Class destructor.
    */
   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.
    */

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<framework::config::Parameters>(
+          "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<std::string>("class_name"),
       modelParameters.getParameter<std::string>("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);
 }