Set the GEANT4 and CUDA RNG seeds using current time and process ID if not set on command line.

2 files changed, 25 insertions, 8 deletions

diff --git a/g4gen.py b/g4gen.py
index 409570f..718df6e 100644
--- a/g4gen.py
+++ b/g4gen.py
@@ -18,7 +18,7 @@ except:
     import G4chroma
 
 class G4Generator(object):
-    def __init__(self, material):
+    def __init__(self, material, seed=None):
         '''Create generator to produce photons inside the specified material.
 
            material: chroma.geometry.Material object with density, 
@@ -26,7 +26,12 @@ class G4Generator(object):
 
                      composition dictionary should be 
                         { element_symbol : fraction_by_weight, ... }.
+           seed: Random number generator seed for HepRandom. If None,
+                 generator is not seeded.
         '''
+        if seed is not None:
+            HepRandom.setTheSeed(seed)
+
         g4py.NISTmaterials.Construct()
         g4py.ezgeom.Construct()
         self.physics_list = G4chroma.ChromaPhysicsList()
diff --git a/sim.py b/sim.py
index 9f89cd1..26a6b68 100755
--- a/sim.py
+++ b/sim.py
@@ -2,6 +2,7 @@
 import sys
 import optparse
 import time
+import os
 import multiprocessing
 
 import detectors
@@ -13,6 +14,11 @@ import numpy as np
 import math
 import ROOT
 
+def pick_seed():
+    '''Returns a seed for a random number generator selected using
+    a mixture of the current time and the current process ID.'''
+    return int(time.time()) ^ (os.getpid() << 16)
+
 def info(type, value, tb):
     if hasattr(sys, 'ps1') or not sys.stderr.isatty():
        # we are in interactive mode or we don't have a tty-like
@@ -26,10 +32,8 @@ def info(type, value, tb):
        # ...then start the debugger in post-mortem mode.
        pdb.pm()
 
-sys.excepthook = info
-
 class GeneratorProcess(multiprocessing.Process):
-    def __init__(self, particle, energy, position, direction, nevents, material):
+    def __init__(self, particle, energy, position, direction, nevents, material, seed=None):
         multiprocessing.Process.__init__(self)
 
         self.particle = particle
@@ -38,11 +42,12 @@ class GeneratorProcess(multiprocessing.Process):
         self.direction = direction
         self.nevents = nevents
         self.material = material
+        self.seed = seed
         self.queue = multiprocessing.Queue()
 
     def run(self):
         print >>sys.stderr, 'Starting generator thread...'
-        generator = g4gen.G4Generator(self.material)
+        generator = g4gen.G4Generator(self.material, seed=self.seed)
 
         for i in xrange(self.nevents):
             photons = generator.generate_photons(particle_name=self.particle, 
@@ -87,7 +92,8 @@ def main():
     parser.add_option('-b', type='int', dest='nbits', default=10)
     parser.add_option('-j', type='int', dest='device', default=None)
     parser.add_option('-n', type='int', dest='nblocks', default=64)
-
+    parser.add_option('-s', type='int', dest='seed', default=None,
+                      help='Set random number generator seed')
     parser.add_option('--detector', type='string', dest='detector', default='microlbne')
     parser.add_option('--nevents', type='int', dest='nevents', default=100)    
     parser.add_option('--particle', type='string', dest='particle', default='e-')
@@ -115,6 +121,10 @@ def main():
     position = np.array(eval(options.pos), dtype=float)
     direction = np.array(eval(options.dir), dtype=float)
     detector = detectors.find(options.detector)
+    if options.seed is None:
+        options.seed = pick_seed()
+
+    print >>sys.stderr, 'RNG seed:', options.seed
 
     print >>sys.stderr, 'Creating BVH for detector "%s" with %d bits...' % (options.detector, options.nbits)
     detector.build(bits=options.nbits)
@@ -126,7 +136,8 @@ def main():
                                         position=position,
                                         direction=direction,
                                         nevents=options.nevents,
-                                        material=detector_material)
+                                        material=detector_material,
+                                        seed=options.seed)
     print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WCSIM WATER!!'
 
     # Do this now so we can get ahead of the photon propagation
@@ -140,7 +151,7 @@ def main():
     gpu_worker.load_geometry(detector)
     
     print >>sys.stderr, 'Initializing random numbers generators...'
-    gpu_worker.setup_propagate()
+    gpu_worker.setup_propagate(seed=options.seed)
     gpu_worker.setup_daq(max(detector.pmtids))
 
     # Create output file
@@ -191,4 +202,5 @@ def main():
     print >>sys.stderr, 'Done. %1.1f events/sec, %1.0f photons/sec.' % (options.nevents/(end_sim - start_sim), nphotons/(end_sim - start_sim))
 
 if __name__ == '__main__':
+    sys.excepthook = info
     main()