multi-threading working

common/framework/include/OMSim.hh

@@ -19,7 +19,7 @@
 #include "OMSimSteppingAction.hh"
 #include "OMSimUIinterface.hh"
 
-#include <G4RunManager.hh>
+#include <G4MTRunManager.hh>
 #include <G4VisExecutive.hh>
 #include <G4UIExecutive.hh>
 
@@ -53,28 +53,24 @@ public:
     void startVisualisationIfRequested();
     // OMSimDetectorConstruction* getDetectorConstruction();
 
-    G4Navigator *getNavigator() { return mNavigator; };
+    G4Navigator *getNavigator() { return mNavigator.get(); };
     void extendOptions(po::options_description pNewOptions);
     po::options_description mGeneralOptions;
 
 private:
     void initialLoggerConfiguration();
+    int determineNumberOfThreads();
     po::variables_map parseArguments(int pArgumentCount, char *pArgumentVector[]);
     void setUserArgumentsToArgTable(po::variables_map pVariablesMap);
     void setGeneralOptions();
 
-    G4RunManager *mRunManager = nullptr;
-    G4VisExecutive *mVisManager = nullptr;
-    G4Navigator *mNavigator = nullptr;
-    G4VUserPhysicsList *mPhysics = nullptr;
-    G4VUserPrimaryGeneratorAction *mGenAction = nullptr;
-    G4UserRunAction *mRunAction = nullptr;
-    G4UserEventAction *mEventAction = nullptr;
-    G4UserTrackingAction *mTracking = nullptr;
-    G4UserSteppingAction *mStepping = nullptr;
-    G4TouchableHistory *mHistory = nullptr;
+    std::unique_ptr<G4MTRunManager> mRunManager;
+    std::unique_ptr<G4VisExecutive> mVisManager;
+    std::unique_ptr<G4VUserPhysicsList> mPhysics;
+    std::unique_ptr<G4TouchableHistory> mHistory;
+    std::unique_ptr<G4Navigator> mNavigator;
 
-    G4double mStartingTime = 0;
+    std::chrono::high_resolution_clock::time_point mStartingTime;
 };
 
 #endif // OMSIM_H

common/framework/include/OMSimActionInitialization.hh

@@ -0,0 +1,19 @@
+#ifndef OMSimActionInitialization_h
+#define OMSimActionInitialization_h 1
+
+#include "G4VUserActionInitialization.hh"
+
+class OMSimActionInitialization : public G4VUserActionInitialization
+{
+public:
+    OMSimActionInitialization(long pSeed);
+    virtual ~OMSimActionInitialization();
+
+    virtual void BuildForMaster() const;
+    virtual void Build() const;
+
+private:
+    long mMasterSeed;
+};
+
+#endif

common/framework/include/OMSimDataFileTypes.hh

@@ -41,6 +41,7 @@ class abcDataFile
 {
 public:
     abcDataFile(G4String pFileName);
+    ~abcDataFile();
     G4String mFileName;
     G4String mObjectName;
 
@@ -66,7 +67,6 @@ public:
     abcMaterialData(G4String pFileName) : abcDataFile(pFileName){};
     G4Material *mMaterial;
     G4MaterialPropertiesTable *mMPT;
-    G4NistManager *mMatDatBase;
 
     void createMaterial();
     void extractAbsorptionLength();

common/framework/include/OMSimHitManager.hh

@@ -100,7 +100,6 @@ private:
         std::map<G4int, HitStats> moduleHits;
     };
     G4ThreadLocal static ThreadLocalData* mThreadData;
-    G4bool mDataWasMerged = false;
 };
 
 #endif

common/framework/include/OMSimRunAction.hh

@@ -2,7 +2,8 @@
 #define OMSimRunAction_h 1
 
 #include <G4UserRunAction.hh>
-
+#include "OMSimHitManager.hh"
+#include "OMSimLogger.hh"
 class G4Run;
 
 class OMSimRunAction : public G4UserRunAction
@@ -12,8 +13,13 @@ public:
   ~OMSimRunAction(){};
 
 public:
-	void BeginOfRunAction(const G4Run*){};
-	void EndOfRunAction(const G4Run*){};
+  void BeginOfRunAction(const G4Run *) {};
+  void EndOfRunAction(const G4Run *run) override
+  {
+    log_debug("EndOfRunAction called");
+    // Merge thread data here
+    OMSimHitManager::getInstance().mergeThreadData();
+  }
 };
 
 #endif

common/framework/include/OMSimSensitiveDetector.hh

@@ -51,7 +51,7 @@ class OMSimSensitiveDetector : public G4VSensitiveDetector
 {
 public:
     OMSimSensitiveDetector(G4String pName, DetectorType pDetectorType);
-    virtual ~OMSimSensitiveDetector() {};
+    ~OMSimSensitiveDetector() {};
 
     G4bool ProcessHits(G4Step *pStep, G4TouchableHistory *pTouchableHistory) override;
     void setPMTResponse(OMSimPMTResponse *pResponse);

common/framework/include/OMSimUIinterface.hh

@@ -57,7 +57,7 @@ public:
      */
     void runBeamOn(G4int pNumberOfEvents = -1)
     {
-        log_debug("Running beamOn command");
+        log_trace("Running beamOn command");
         G4int lNumEvents = pNumberOfEvents >= 0 ? pNumberOfEvents : OMSimCommandArgsTable::getInstance().get<G4int>("numevents");
         applyCommand("/run/beamOn ", lNumEvents);
     }

common/framework/src/OMSim.cc

@@ -10,12 +10,18 @@
 
 #include "OMSim.hh"
 #include "OMSimLogger.hh"
+#include "OMSimActionInitialization.hh"
 #include <boost/program_options.hpp>
 
 namespace po = boost::program_options;
 extern std::shared_ptr<spdlog::logger> globalLogger;
 
-OMSim::OMSim() : mStartingTime(clock() / CLOCKS_PER_SEC), mGeneralOptions("General options"), mRunManager(new G4RunManager()), mVisManager(new G4VisExecutive()), mNavigator(new G4Navigator())
+OMSim::OMSim() : 
+mStartingTime(std::chrono::high_resolution_clock::now()), 
+mGeneralOptions("General options"), 
+mRunManager(std::make_unique<G4MTRunManager>()),
+mVisManager(std::make_unique<G4VisExecutive>()),
+mNavigator(std::make_unique<G4Navigator>())
 {
     setGeneralOptions();
     initialLoggerConfiguration();
@@ -41,7 +47,8 @@ void OMSim::setGeneralOptions()
     ("glass", po::value<G4int>()->default_value(0), "DEPRECATED. Index to select glass type [VITROVEX = 0, Chiba = 1, Kopp = 2, myVitroVex = 3, myChiba = 4, WOMQuartz = 5, fusedSilica = 6]")
     ("gel", po::value<G4int>()->default_value(1), "DEPRECATED. Index to select gel type [Wacker = 0, Chiba = 1, IceCube = 2, Wacker_company = 3]")
     ("reflective_surface", po::value<G4int>()->default_value(0), "DEPRECATED. Index to select reflective surface type [Refl_V95Gel = 0, Refl_V98Gel = 1, Refl_Aluminium = 2, Refl_Total98 = 3]")
-    ("pmt_model", po::value<G4int>()->default_value(0), "DEPRECATED. R15458 (mDOM) = 0,  R7081 (DOM) = 1, 4inch (LOM) = 2, R5912_20_100 (D-Egg)= 3");
+    ("pmt_model", po::value<G4int>()->default_value(0), "DEPRECATED. R15458 (mDOM) = 0,  R7081 (DOM) = 1, 4inch (LOM) = 2, R5912_20_100 (D-Egg)= 3")
+    ("threads", po::value<int>()->default_value(1), "number of threads to use.");
 }
 
 void OMSim::initialLoggerConfiguration()
@@ -117,6 +124,21 @@ OMSimDetectorConstruction* OMSim::getDetectorConstruction()
 }
 */
 
+int OMSim::determineNumberOfThreads()
+{
+    int requestedThreads = OMSimCommandArgsTable::getInstance().get<int>("threads");
+    int availableThreads = G4Threading::G4GetNumberOfCores();
+
+    if (requestedThreads <= 0) {
+        log_info("Auto-detected {} available cores", availableThreads);
+        return availableThreads;
+    } else {
+        int threadsToUse = std::min(requestedThreads, availableThreads);
+        log_info("Using {} out of {} available cores", threadsToUse, availableThreads);
+        return threadsToUse;
+    }
+}
+
 /**
  * @brief Initialize the simulation constructing all Geant instances.
  */
@@ -131,41 +153,35 @@ void OMSim::initialiseSimulation(OMSimDetectorConstruction* pDetectorConstructio
     if (lArgs.get<bool>("save_args"))
         lArgs.writeToJson(lFileName);
 
-    CLHEP::HepRandom::setTheEngine(new CLHEP::RanluxEngine(lArgs.get<long>("seed"), 3));
+    //CLHEP::HepRandom::setTheEngine(new CLHEP::RanluxEngine(lArgs.get<long>("seed"), 3));
+    //CLHEP::HepRandom::setTheEngine(new CLHEP::MixMaxRng(lArgs.get<long>("seed")));
+    //G4Random::setTheEngine(new CLHEP::RanluxEngine(lArgs.get<long>("seed"), 3));
+
 
     mRunManager->SetUserInitialization(pDetectorConstruction);
 
-    mPhysics = new OMSimPhysicsList;
-    mRunManager->SetUserInitialization(mPhysics);
+    mPhysics = std::make_unique<OMSimPhysicsList>();
+    mRunManager->SetUserInitialization(mPhysics.get());
+    mPhysics.release();
 
     mVisManager->Initialize();
 
-    mGenAction = new OMSimPrimaryGeneratorAction();
-    mRunManager->SetUserAction(mGenAction);
-
-    mRunAction = new OMSimRunAction();
-    mRunManager->SetUserAction(mRunAction);
-
-    mEventAction = new OMSimEventAction();
-    mRunManager->SetUserAction(mEventAction);
-
-    mTracking = new OMSimTrackingAction();
-    mRunManager->SetUserAction(mTracking);
+    OMSimActionInitialization* actionInitialization = new OMSimActionInitialization(lArgs.get<long>("seed"));
+    mRunManager->SetUserInitialization(actionInitialization);
 
-    mStepping = new OMSimSteppingAction();
-    mRunManager->SetUserAction(mStepping);
+    // Set number of threads
+    int nThreads = determineNumberOfThreads();
+    mRunManager->SetNumberOfThreads(nThreads);
 
     mRunManager->Initialize();
 
     OMSimUIinterface &lUIinterface = OMSimUIinterface::getInstance();
     lUIinterface.setUI(G4UImanager::GetUIpointer());
 
-    mNavigator->SetWorldVolume(pDetectorConstruction->mWorldPhysical);
-    mNavigator->LocateGlobalPointAndSetup(G4ThreeVector(0., 0., 0.));
+    mNavigator.get()->SetWorldVolume(pDetectorConstruction->mWorldPhysical);
+    mNavigator.get()->LocateGlobalPointAndSetup(G4ThreeVector(0., 0., 0.));
 
-    mHistory = mNavigator->CreateTouchableHistory();
-
-    lUIinterface.applyCommand("/control/execute ", lArgs.get<bool>("visual"));
+    mHistory = std::unique_ptr<G4TouchableHistory>(mNavigator->CreateTouchableHistory());
 }
 
 /**
@@ -238,23 +254,45 @@ bool OMSim::handleArguments(int pArgumentCount, char *pArgumentVector[])
 OMSim::~OMSim()
 {
     log_trace("OMSim destructor");
-    if (mRunManager) {
-        log_trace("Deleting RunManager");
-        delete mRunManager;
-        mRunManager = nullptr;
-    }
-
-    if (mVisManager) {
-        log_trace("Deleting VisManager");
-        delete mVisManager;
-        mVisManager = nullptr;
-    }
-
-    if (mNavigator) {
-        log_trace("Deleting Navigator");
-        delete mNavigator;
-        mNavigator = nullptr;
-    }
-    double lFinishtime = clock() / CLOCKS_PER_SEC;
-    log_info("Computation time: {} {}", lFinishtime - mStartingTime, " seconds.");
+    // if (mRunManager) {
+    //     log_trace("Deleting RunManager");
+    //     delete mRunManager;
+    //     mRunManager = nullptr;
+    // }
+
+    // if (mVisManager) {
+    //     log_trace("Deleting VisManager");
+    //     delete mVisManager;
+    //     mVisManager = nullptr;
+    // }
+
+    // if (mNavigator) {
+    //     log_trace("Deleting Navigator");
+    //     delete mNavigator;
+    //     mNavigator = nullptr;
+    // }
+
+    log_trace("OMSim destructor started");
+
+    log_trace("Resetting mHistory");
+    mHistory.reset();
+
+    log_trace("Resetting mNavigator");
+    mNavigator.reset();
+
+    log_trace("Resetting mPhysics");
+    mPhysics.reset();
+
+    log_trace("Resetting mVisManager");
+    mVisManager.reset();
+
+    log_trace("Resetting mRunManager");
+    mRunManager.reset();
+
+    log_trace("OMSim destructor finished");
+
+
+    std::chrono::high_resolution_clock::time_point lFinishtime = std::chrono::high_resolution_clock::now();
+    const std::chrono::duration<double> lDiff = lFinishtime - mStartingTime;
+    log_info("Computation time: {} {}", lDiff.count(), " seconds.");
 }

common/framework/src/OMSimActionInitialization.cc

@@ -0,0 +1,36 @@
+#include "OMSimActionInitialization.hh"
+#include "OMSimPrimaryGeneratorAction.hh"
+#include "OMSimRunAction.hh"
+#include "OMSimEventAction.hh"
+#include "OMSimTrackingAction.hh"
+#include "OMSimSteppingAction.hh"
+#include <G4RandomTools.hh>
+
+OMSimActionInitialization::OMSimActionInitialization(long pSeed)
+    : G4VUserActionInitialization(), mMasterSeed(pSeed)
+{
+}
+
+OMSimActionInitialization::~OMSimActionInitialization()
+{
+}
+
+void OMSimActionInitialization::BuildForMaster() const
+{
+    SetUserAction(new OMSimRunAction);
+}
+
+void OMSimActionInitialization::Build() const
+{
+    SetUserAction(new OMSimPrimaryGeneratorAction);
+    SetUserAction(new OMSimRunAction);
+    SetUserAction(new OMSimEventAction);
+    SetUserAction(new OMSimTrackingAction);
+    SetUserAction(new OMSimSteppingAction);
+    const long lPrime = 2147483647; 
+    long lSeed = (mMasterSeed + (G4Threading::G4GetThreadId()+1) * lPrime);
+    lSeed = lSeed % std::numeric_limits<long>::max();
+    log_debug("Random engine of thread {} was assigned seed {}", G4Threading::G4GetThreadId(), lSeed);
+    G4Random::setTheSeed(lSeed);
+    G4Random::setTheEngine(new CLHEP::RanluxEngine);
+}