From d8f3a3720a5f04156f71da445fc969e6ca76a943 Mon Sep 17 00:00:00 2001
From: Jan Kieseler <jkiesele@cern.ch>
Date: Thu, 22 Aug 2024 15:50:24 +0200
Subject: [PATCH] proper seedings

---
 bind/G4Calo.py      | 50 ++++++++++++++++++++-------------------------
 bind/example.py     |  2 +-
 include/G4System.hh |  2 +-
 src/G4System.cc     |  6 ++++--
 4 files changed, 28 insertions(+), 32 deletions(-)

diff --git a/bind/G4Calo.py b/bind/G4Calo.py
index a57e0b0..d8f5d4e 100644
--- a/bind/G4Calo.py
+++ b/bind/G4Calo.py
@@ -1,5 +1,7 @@
 import multiprocessing
 import warnings
+import contextlib
+import sys
 
 def set_start_method():
     try:
@@ -312,43 +314,22 @@ def _run_mini_batch(
         particleSpec: str,
         minEnergy_GeV: float,
         maxEnergy_GeV: float = -1.0,
+        batch_seed : int = 0,
         counter : int = 0):
 
     
-    print(f"Running mini batch {counter} with {nEvents} events")
-    import time
-    time.sleep(counter/20)
-    print(f'done sleeping {counter}')
-
     #this is now encapsuled
     from  G4Calo import __G4System
     G4System = __G4System()
-    G4System.init(cw)
+    seed = int(batch_seed + counter)
+    G4System.init(cw, seed) #counter gives random seed offset
 
     df = G4System.run_batch(nEvents, particleSpec, minEnergy_GeV, maxEnergy_GeV,"")
     return df
 
 def _run_mini_batch_silent(*args, **kwargs):
-    import sys
-    # Redirect stdout and stderr to silence the output
-    devnull = open(os.devnull, 'w')
-    original_stdout = sys.stdout
-    original_stderr = sys.stderr
-    sys.stdout = devnull
-    sys.stderr = devnull
-
-    try:
-        return _run_mini_batch(*args, **kwargs)  # Call your actual function
-    except Exception as e:
-        # Restore original stdout and stderr if an exception occurs
-        sys.stdout = original_stdout
-        sys.stderr = original_stderr
-        raise e
-    finally:
-        # Restore original stdout and stderr
-        sys.stdout = original_stdout
-        sys.stderr = original_stderr
-        devnull.close()
+    #I tried
+    return _run_mini_batch(*args, **kwargs)
 
 
 def run_batch(
@@ -375,11 +356,24 @@ def run_batch(
     
     nevents = [nEventsPerCore if i < nCores - 1 else nEventsLastCore for i in range(nCores)]
 
+    batch_seed = int(time.time()*1000)
+    #truncate to 32 bit, lower 32 bit are used in G4
+    batch_seed = batch_seed & 0xFFFFFFFF
+
+    print(f"Batch seed: {batch_seed}")
     #use a multiprocessing pool to run the mini batches in parallel
     with multiprocessing.Pool(nCores) as pool:
-        dfs = pool.starmap(_run_mini_batch_silent, [(gd, nevents[i], particleSpec, minEnergy_GeV, maxEnergy_GeV, i) for i in range(nCores)])
+        dfs = pool.starmap(_run_mini_batch_silent, [(gd, nevents[i], particleSpec, minEnergy_GeV, maxEnergy_GeV, batch_seed, i) for i in range(nCores)])
 
-    return pd.concat(dfs)
+    alldf = pd.concat(dfs)
+    #convert alldf['total_dep_energy'] to numpy array and check if all non-zero entries are unique
+    non_zero = alldf['total_dep_energy'][alldf['total_dep_energy'] != 0.]
+    n_unique = np.unique(non_zero.to_numpy()).shape[0]
+    n_ex_non_zero = non_zero.shape[0]
+
+    assert (n_ex_non_zero==n_unique), f"only {n_unique} unique non-zero entries in total_dep_energy out of {n_ex_non_zero}"
+
+    return alldf
 
 def _fill_event(gd : GeometryDescriptor, 
                  particleSpec: str,
diff --git a/bind/example.py b/bind/example.py
index 766d551..236f827 100644
--- a/bind/example.py
+++ b/bind/example.py
@@ -18,4 +18,4 @@ if __name__ == '__main__':
         gd.addLayer(1.,"G4_POLYSTYRENE",True,1)
     
     df = run_batch(gd,10, 'gamma', 1)
-    display_event(gd,"gamma", 2)
\ No newline at end of file
+    #display_event(gd,"gamma", 2)
\ No newline at end of file
diff --git a/include/G4System.hh b/include/G4System.hh
index 1d332d5..2b305fe 100644
--- a/include/G4System.hh
+++ b/include/G4System.hh
@@ -34,7 +34,7 @@ G4System(bool Gui):gui(Gui){};
 };
 
 
-void init(GeometryDescriptor &cw);
+void init(GeometryDescriptor &cw,  int seed=0);
 //will use dawn for visualization, also wrap more in python
 void run_visualize(const std::string& partSpecies, double minEnergy_GeV, double maxEnergy_GeV);
 //runs the whole gui if available
diff --git a/src/G4System.cc b/src/G4System.cc
index 9137f8f..72ac38c 100644
--- a/src/G4System.cc
+++ b/src/G4System.cc
@@ -9,7 +9,7 @@
 
 #include <ctime>
 
-void G4System::init(GeometryDescriptor &CW){
+void G4System::init(GeometryDescriptor &CW, int seed){
 
   if(CW.isAssigned()){
     throw std::runtime_error("GeometryDescriptor already assigned");
@@ -28,7 +28,9 @@ void G4System::init(GeometryDescriptor &CW){
 
   bool first_init=true;
 
-  int seed=std::time(NULL) ;
+  if(seed < 0){
+    seed = std::time(NULL);
+  }
   CLHEP::HepRandom::setTheSeed(seed); G4Random::setTheSeed(seed);