diff --git a/Simulation/G4Extensions/RHadrons/CMakeLists.txt b/Simulation/G4Extensions/RHadrons/CMakeLists.txt
index a60a2bb20bdd0a6d21bab873813489f4906229a3..f47e493bc2279a939575b7e917271721d32f0a5c 100644
--- a/Simulation/G4Extensions/RHadrons/CMakeLists.txt
+++ b/Simulation/G4Extensions/RHadrons/CMakeLists.txt
@@ -20,7 +20,7 @@ if( NOT GENERATIONBASE )
                        OBJECT
                        NO_PUBLIC_HEADERS
                        PRIVATE_INCLUDE_DIRS ${GEANT4_INCLUDE_DIRS} ${XERCESC_INCLUDE_DIRS} ${CLHEP_INCLUDE_DIRS} ${PYTHIA8_INCLUDE_DIRS}
-                       PRIVATE_LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} GaudiKernel AthenaBaseComps G4AtlasInterfaces G4AtlasToolsLib G4ExternalDecay SimHelpers Pythia8_iLib )
+                       PRIVATE_LINK_LIBRARIES ${GEANT4_LIBRARIES} ${XERCESC_LIBRARIES} ${CLHEP_LIBRARIES} ${PYTHIA8_LIBRARIES} GaudiKernel AthenaBaseComps CxxUtils G4AtlasInterfaces G4AtlasToolsLib G4ExternalDecay SimHelpers Pythia8_iLib )
 endif()
 
 # Install files from the package:
diff --git a/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.cxx b/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.cxx
index a8d21adaf4620493121b098545709fe9ada8eacb..b68ab4a39de53b3d7d0329a7417942215a53f2fb 100644
--- a/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.cxx
+++ b/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 #include <fstream>
@@ -17,18 +17,22 @@
 #include "G4DecayTable.hh"
 #include "G4PhaseSpaceDecayChannel.hh"
 
-bool CustomParticleFactory::loaded = false;
-std::set<G4ParticleDefinition *> CustomParticleFactory::m_particles;
-
 bool CustomParticleFactory::isCustomParticle(G4ParticleDefinition *particle)
 {
-  return (m_particles.find(particle)!=m_particles.end());
+  static const std::set<G4ParticleDefinition *> particles = load();
+  return (particle!=nullptr && particles.find(particle)!=particles.end());
 }
 
 void CustomParticleFactory::loadCustomParticles()
 {
-  if(loaded) return;
-  loaded = true;
+  // tickle the loading of the particles if it wasn't done yet
+  isCustomParticle(nullptr);
+}
+
+std::set<G4ParticleDefinition *> CustomParticleFactory::load()
+{
+  std::set<G4ParticleDefinition *> particles;
+
   std::ifstream configFile("particles.txt");
   G4String pType="custom";
   G4String pSubType="";
@@ -117,7 +121,7 @@ void CustomParticleFactory::loadCustomParticles()
             <<G4endl;
     }
 
-    m_particles.insert(particle);
+    particles.insert(particle);
   }
   configFile.close();
 
@@ -132,7 +136,7 @@ void CustomParticleFactory::loadCustomParticles()
   decayFile.close();
 
   // Looping over custom particles to add decays
-  for (std::set<G4ParticleDefinition *>::iterator part=m_particles.begin();part!=m_particles.end();part++) {
+  for (std::set<G4ParticleDefinition *>::iterator part=particles.begin();part!=particles.end();part++) {
     name=(*part)->GetParticleName();
     std::vector<std::vector<std::string> > mydecays;
     for (unsigned int i = 0; i!= decays.size(); i++){
@@ -163,5 +167,5 @@ void CustomParticleFactory::loadCustomParticles()
       (*part)->SetDecayTable(table);
     }
   }
-  return;
+  return particles;
 }
diff --git a/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.h b/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.h
index 9d429fd9590ac998bb04b2d620ea1152cad3c17a..dcb0fec10a8b1d1c2484293e6e8a0ed07a057faf 100644
--- a/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.h
+++ b/Simulation/G4Extensions/RHadrons/src/CustomParticleFactory.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef CustomParticleFactory_h
@@ -15,15 +15,13 @@
 
 class CustomParticleFactory
 {
-private:
-  static bool loaded;
-  static std::set<G4ParticleDefinition *> m_particles;
 
 public:
   static void loadCustomParticles();
   static bool isCustomParticle(G4ParticleDefinition *particle);
 
-
+private:
+  static std::set<G4ParticleDefinition *> load();
 };
 
 
diff --git a/Simulation/G4Extensions/RHadrons/src/FullModelHadronicProcess.cxx b/Simulation/G4Extensions/RHadrons/src/FullModelHadronicProcess.cxx
index 6a23da42ebc86ef8db7c0c008d9348e4714a5caf..9c2ae20a812c0dfd3faccc699e0398bba2b54823 100644
--- a/Simulation/G4Extensions/RHadrons/src/FullModelHadronicProcess.cxx
+++ b/Simulation/G4Extensions/RHadrons/src/FullModelHadronicProcess.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "FullModelHadronicProcess.hh"
@@ -82,7 +82,7 @@ G4VParticleChange* FullModelHadronicProcess::PostStepDoIt(const G4Track& aTrack,
   // A little setting up
   aParticleChange.Initialize(aTrack);
   //  G4DynamicParticle* OrgPart = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
-  G4DynamicParticle* IncidentRhadron = const_cast<G4DynamicParticle*>(aTrack.GetDynamicParticle());
+  const G4DynamicParticle* IncidentRhadron = aTrack.GetDynamicParticle();
   CustomParticle* CustomIncident = static_cast<CustomParticle*>(IncidentRhadron->GetDefinition());
   const G4ThreeVector& aPosition = aTrack.GetPosition();
   const G4int theIncidentPDG = IncidentRhadron->GetDefinition()->GetPDGEncoding();
diff --git a/Simulation/G4Extensions/RHadrons/src/FullModelReactionDynamics.cxx b/Simulation/G4Extensions/RHadrons/src/FullModelReactionDynamics.cxx
index 67b8e06d4d5f20b5d4abb4b9c0d93f5742967a9c..6a40c30f1ca3a339d8885d9c43ec2334eef740b9 100644
--- a/Simulation/G4Extensions/RHadrons/src/FullModelReactionDynamics.cxx
+++ b/Simulation/G4Extensions/RHadrons/src/FullModelReactionDynamics.cxx
@@ -1021,7 +1021,8 @@ G4bool FullModelReactionDynamics::GenerateXandPt(
           G4cerr << "currentParticle side = " << currentParticle.GetSide() << G4endl;
           for( i=0; i<vecLen; ++i )
             G4cerr << "particle #" << i << " side = " << vec[i]->GetSide() << G4endl;
-          exit( EXIT_FAILURE );
+          throw std::runtime_error("FullModelReactionDynamics::GenerateXandPt: "
+                                   "tempLen is not the same as backwardNucleonCount");
         }
       constantCrossSection = true;
       // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
diff --git a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx
index 4b3e490764a8cdf8be348a0a36bcdee881bc12f2..6aa4cd3d6f8387cd03fae83120ce3170e7baa7e0 100644
--- a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx
+++ b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "G4ProcessHelper.hh"
@@ -13,6 +13,8 @@
 #include <fstream>
 #include <stdexcept>
 
+#include "CxxUtils/checker_macros.h"
+
 G4ProcessHelper::G4ProcessHelper()
   : theTarget(0)
   , theRmesoncloud(0)
@@ -160,25 +162,23 @@ G4ProcessHelper::G4ProcessHelper()
   return;
 }
 
-G4ProcessHelper* G4ProcessHelper::pinstance = 0;
-
 G4ProcessHelper* G4ProcessHelper::Instance()
 {
-  if (pinstance == 0)
-    {
-      pinstance = new G4ProcessHelper();
-    }
-  return pinstance;
+  static G4ProcessHelper instance;
+  return &instance;
 }
 
-G4bool G4ProcessHelper::ApplicabilityTester(const G4ParticleDefinition& aPart){
-  const G4ParticleDefinition* aP = &aPart;
-  if (known_particles[aP]) return true;
-  return false;
+G4bool G4ProcessHelper::ApplicabilityTester(const G4ParticleDefinition& aPart) const {
+  try {
+    return known_particles.at(&aPart);
+  }
+  catch (const std::out_of_range& e) {
+    return false;
+  }
 }
 
 G4double G4ProcessHelper::GetInclusiveCrossSection(const G4DynamicParticle *aParticle,
-                                                   const G4Element *anElement){
+                                                   const G4Element *anElement) const {
   //We really do need a dedicated class to handle the cross sections. They might not always be constant
 
   //Disassemble the PDG-code
@@ -484,7 +484,8 @@ G4double G4ProcessHelper::ReactionProductMass(const ReactionProduct& aReaction,c
   G4double M_after = 0;
   for (ReactionProduct::const_iterator r_it = aReaction.begin(); r_it !=aReaction.end(); r_it++){
     //G4cout<<"Mass contrib: "<<(particleTable->FindParticle(*r_it)->GetPDGMass())/CLHEP::MeV<<" MeV"<<G4endl;
-    M_after += particleTable->FindParticle(*r_it)->GetPDGMass();
+    auto table ATLAS_THREAD_SAFE = particleTable;  // safe because table has been loaded by now
+    M_after += table->FindParticle(*r_it)->GetPDGMass();
   }
   //G4cout<<"Intending to return this ReactionProductMass: " << sqrts << " - " <<  M_after << " MeV"<<G4endl;
   return sqrts - M_after;
diff --git a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh
index ee30bd1aff719be2cdd556b5b935ab28b461651c..29c3b69c20c75241b86cd935958fccccc32b2818 100644
--- a/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh
+++ b/Simulation/G4Extensions/RHadrons/src/G4ProcessHelper.hh
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef RHADRONS_G4PROCESSHELPER_HH
@@ -25,10 +25,10 @@ class G4ProcessHelper {
 public:
   static G4ProcessHelper* Instance();
 
-  G4bool ApplicabilityTester(const G4ParticleDefinition& aPart);
+  G4bool ApplicabilityTester(const G4ParticleDefinition& aPart) const;
 
   G4double GetInclusiveCrossSection(const G4DynamicParticle *aParticle,
-                                    const G4Element *anElement);
+                                    const G4Element *anElement) const;
 
   //Make sure the element is known (for n/p-decision)
   ReactionProduct GetFinalState(const G4Track& aTrack,G4ParticleDefinition*& aTarget);
@@ -52,8 +52,6 @@ private:
 
   ReactionMap* theReactionMap;
 
-  static G4ProcessHelper* pinstance;
-
   // Version where we know if we baryonize already
   ReactionProduct GetFinalStateInternal(const G4Track& aTrack,G4ParticleDefinition*& aTarget, const bool baryonize_failed);