added early stop option for decays chains

CMakeLists.txt

@@ -38,8 +38,8 @@ include(${ROOT_USE_FILE})
 
 # Include the subdirectories
 add_subdirectory(common)
-add_subdirectory(effective_area)
-#add_subdirectory(radioactive_decays)
+#add_subdirectory(effective_area)
+add_subdirectory(radioactive_decays)
 #add_subdirectory(supernova)
 
 # Copy auxiliary files from source directory to binary directory

radioactive_decays/OMSim_radioactive_decays.cc

@@ -21,14 +21,14 @@ void decaySimulation(OMSimRadDecaysDetector *pDetector)
 
 	lAnalysisManager.setOutputFileName(lArgs.get<std::string>("output_file"));
 
-	IsotopeDecays *lDecays = new IsotopeDecays(280);
+	OMSimDecaysGPS &lDecays = OMSimDecaysGPS::getInstance();
 
-	lDecays->setOpticalModule(pDetector->mOpticalModule);
+	lDecays.setOpticalModule(pDetector->mOpticalModule);
 
 	for (int i = 0; i < (int)lArgs.get<G4int>("numevents"); i++)
 	{
 
-		lDecays->simulateDecaysInOpticalModule(lArgs.get<G4double>("time_window"));
+		lDecays.simulateDecaysInOpticalModule(lArgs.get<G4double>("time_window"));
 
 		if (lArgs.get<bool>("multiplicity_study"))
 		{

radioactive_decays/include/OMSimDecaysGPS.hh

@@ -1,6 +1,6 @@
 /**
  * @file
- * @brief Defines the IsotopeDecays class for the radioactive decays simulation.
+ * @brief Defines the OMSimDecaysGPS class for the radioactive decays simulation.
  * @ingroup radioactive
  */
 
@@ -9,40 +9,46 @@
 #include "OMSimOpticalModule.hh"
 #include <globals.hh>
 /**
- * @class IsotopeDecays
+ * @class OMSimDecaysGPS
  * @brief A class for simulating isotope decays inside the pressure vessel and PMT glass.
  * @ingroup radioactive
  */
-class IsotopeDecays
+class OMSimDecaysGPS
 {
 public:
-  IsotopeDecays(G4double pProductionRadius);
-  ~IsotopeDecays();
+
+  static OMSimDecaysGPS &getInstance()
+  {
+    static OMSimDecaysGPS instance;
+    return instance;
+  }
   void simulateDecaysInOpticalModule(G4double pTimeWindow);
 
   /**
    * @brief Set the optical module to be used.
    * @param pOpticalModule Pointer to the optical module.
    */
-  void setOpticalModule(OMSimOpticalModule* pOpticalModule){mOM = pOpticalModule;};
+  void setOpticalModule(OMSimOpticalModule *pOpticalModule) { mOM = pOpticalModule; };
+  void setProductionRadius(G4double pProductionRadius);
+  G4String getDecayTerminationNuclide();
 
 private:
-
   // Define a map that maps isotopes to their commands
   std::map<G4String, G4String> mIsotopeCommands = {
       {"U238", "/gps/ion 92 238 0"},
       {"U235", "/gps/ion 92 235 0"},
       {"Th232", "/gps/ion 90 232 0"},
-      {"K40", "/gps/ion 19 40 0"}
-  };
-  void generalGPS(); 
+      {"K40", "/gps/ion 19 40 0"}};
+  void generalGPS();
   void configureIsotopeGPS(G4String Isotope, G4String location, G4int optParam = -999);
-  std::map<G4String, G4int> calculateNumberOfDecays(G4MaterialPropertiesTable* pMPT, G4double pTimeWindow, G4double pMass);
-  OMSimOpticalModule* mOM;
+  std::map<G4String, G4int> calculateNumberOfDecays(G4MaterialPropertiesTable *pMPT, G4double pTimeWindow, G4double pMass);
+  OMSimOpticalModule *mOM;
   G4double mProductionRadius;
+  G4String mNuclideStopName;
+
+  OMSimDecaysGPS() = default;
+  ~OMSimDecaysGPS() = default;
+  OMSimDecaysGPS(const OMSimDecaysGPS &) = delete;
+  OMSimDecaysGPS &operator=(const OMSimDecaysGPS &) = delete;
 };
 #endif
-
-
-
-

radioactive_decays/include/OMSimG4RadioactiveDecay.hh

@@ -178,6 +178,8 @@ class G4RadioactiveDecay : public G4VRestDiscreteProcess
 
     void DecayAnalog(const G4Track& theTrack);
 
+    G4VParticleChange* TerminateDecay();
+
     G4DecayProducts* DoDecay(const G4ParticleDefinition& theParticleDef);
 
     // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)

radioactive_decays/src/OMSimDecaysGPS.cc

@@ -3,15 +3,19 @@
 #include <G4SystemOfUnits.hh>
 #include <G4Poisson.hh>
 
-IsotopeDecays::IsotopeDecays(G4double pProductionRadius)
+void OMSimDecaysGPS::setProductionRadius(G4double pProductionRadius)
 {
     mProductionRadius = pProductionRadius;
 }
 
+G4String OMSimDecaysGPS::getDecayTerminationNuclide()
+{
+    return mNuclideStopName;
+}
 /**
  * @brief Configures common GPS commands for the radioactive decays.
  */
-void IsotopeDecays::generalGPS()
+void OMSimDecaysGPS::generalGPS()
 {
     OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance();
     OMSimCommandArgsTable &lArgs = OMSimCommandArgsTable::getInstance();
@@ -45,7 +49,7 @@ void IsotopeDecays::generalGPS()
  * @param pVolumeName The volume name where the isotope decays.
  * @param optParam An optional parameter related to the location. For "PMT", use e.g. the PMT number.
  */
-void IsotopeDecays::configureIsotopeGPS(G4String Isotope, G4String pVolumeName, G4int optParam)
+void OMSimDecaysGPS::configureIsotopeGPS(G4String Isotope, G4String pVolumeName, G4int optParam)
 {
     OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance();
     generalGPS();
@@ -73,7 +77,7 @@ void IsotopeDecays::configureIsotopeGPS(G4String Isotope, G4String pVolumeName,
  * @param pMass Mass of the volume in which decays occur.
  * @return Map of isotopes and their respective number of decays.
  */
-std::map<G4String, G4int> IsotopeDecays::calculateNumberOfDecays(G4MaterialPropertiesTable *pMPT, G4double pTimeWindow, G4double pMass)
+std::map<G4String, G4int> OMSimDecaysGPS::calculateNumberOfDecays(G4MaterialPropertiesTable *pMPT, G4double pTimeWindow, G4double pMass)
 {
     std::map<G4String, G4int> mNumberDecays;
     for (auto &pair : mIsotopeCommands)
@@ -89,7 +93,7 @@ std::map<G4String, G4int> IsotopeDecays::calculateNumberOfDecays(G4MaterialPrope
  * @brief Simulates the decays in the optical module.
  * @param pTimeWindow The livetime that should be simulated.
  */
-void IsotopeDecays::simulateDecaysInOpticalModule(G4double pTimeWindow)
+void OMSimDecaysGPS::simulateDecaysInOpticalModule(G4double pTimeWindow)
 {
     OMSimCommandArgsTable &lArgs = OMSimCommandArgsTable::getInstance();
     OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance();

radioactive_decays/src/OMSimG4RadioactiveDecay.cc

@@ -38,6 +38,7 @@
 
 #include "OMSimDecaysAnalysis.hh"
 #include "OMSimCommandArgsTable.hh"
+#include "OMSimDecaysGPS.hh"
 
 #include <OMSimG4RadioactiveDecay.hh>
 #include <G4RadioactiveDecayMessenger.hh>
@@ -1017,6 +1018,19 @@ void G4RadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A, G4String filenam
 //                                                                            //
 ////////////////////////////////////////////////////////////////////////////////
 
+// Method to terminate the decay process without producing secondaries
+G4VParticleChange *G4RadioactiveDecay::TerminateDecay()
+{
+  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
+
+  // Kill the parent particle.
+  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
+  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
+  ClearNumberOfInteractionLengthLeft();
+
+  return &fParticleChangeForRadDecay;
+}
+
 G4VParticleChange *
 G4RadioactiveDecay::DecayIt(const G4Track &theTrack, const G4Step &)
 {
@@ -1044,13 +1058,8 @@ G4RadioactiveDecay::DecayIt(const G4Track &theTrack, const G4Step &)
           G4cout << ValidVolumes[i] << G4endl;
       }
 #endif
-      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
 
-      // Kill the parent particle.
-      fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
-      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
-      ClearNumberOfInteractionLengthLeft();
-      return &fParticleChangeForRadDecay;
+      return TerminateDecay();
     }
   }
 
@@ -1068,13 +1077,13 @@ G4RadioactiveDecay::DecayIt(const G4Track &theTrack, const G4Step &)
              << G4endl;
     }
 #endif
-    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
+    return TerminateDecay();
+  }
 
-    // Kill the parent particle
-    fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
-    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
-    ClearNumberOfInteractionLengthLeft();
-    return &fParticleChangeForRadDecay;
+  // Check if the current nucleus is the target stop nucleus
+  if (theParticleDef->GetParticleName() == OMSimDecaysGPS::getInstance().getDecayTerminationNuclide())
+  {
+    return TerminateDecay();
   }
 
   G4DecayTable *theDecayTable = GetDecayTable(theParticleDef);
@@ -1093,13 +1102,7 @@ G4RadioactiveDecay::DecayIt(const G4Track &theTrack, const G4Step &)
              << G4endl;
     }
 #endif
-    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
-
-    // Kill the parent particle.
-    fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
-    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
-    ClearNumberOfInteractionLengthLeft();
-    return &fParticleChangeForRadDecay;
+    return TerminateDecay();
   }
   else
   {
@@ -1239,7 +1242,7 @@ void G4RadioactiveDecay::DecayAnalog(const G4Track &theTrack)
     //////////////////////////////////////////////////////////////////////////////////////////////////
     if (temptime > 10e-1 * s)
     {
-      G4double lSimulatedTime = OMSimCommandArgsTable::getInstance().get<G4double>("time_window")*s;
+      G4double lSimulatedTime = OMSimCommandArgsTable::getInstance().get<G4double>("time_window") * s;
       temptime = G4UniformRand() * lSimulatedTime;
       finalGlobalTime = temptime;
       finalLocalTime = temptime;
@@ -1266,7 +1269,6 @@ void G4RadioactiveDecay::DecayAnalog(const G4Track &theTrack)
     // Note that the cut is not on the average, mean lifetime, but on the actual
     // sampled global decay time.
 
-
     if (finalGlobalTime > fThresholdForVeryLongDecayTime)
     {
       fParticleChangeForRadDecay.SetNumberOfSecondaries(0);