merge heads

10 files changed, 481 insertions, 14 deletions

diff --git a/.hgignore b/.hgignore
index ea5e656..a2dec2a 100644
--- a/.hgignore
+++ b/.hgignore
@@ -1,3 +1,5 @@
 syntax:glob
 *.pyc
-*~
\ No newline at end of file
+*~
+*.so
+*.d
diff --git a/.rootlogon.C b/.rootlogon.C
new file mode 100644
index 0000000..f1b4cf6
--- /dev/null
+++ b/.rootlogon.C
@@ -0,0 +1,3 @@
+{
+  gROOT->ProcessLine(".L io/root.C+g");
+}
diff --git a/detectors/__init__.py b/detectors/__init__.py
index d7583a6..043bdbd 100644
--- a/detectors/__init__.py
+++ b/detectors/__init__.py
@@ -3,6 +3,7 @@ from sno import build_sno as build_sno_detector
 
 import os
 import sys
+import inspect
 
 dir = os.path.split(os.path.realpath(__file__))[0]
 sys.path.append(dir + '/..')
@@ -35,3 +36,15 @@ def build_sno():
 @buildable('real_sno')
 def build_real_sno():
     return build_sno_detector(real_av=True)
+
+
+def find(detector_name):
+    members = globals()
+    buildable_lookup = {}
+    for member in members.values():
+        if inspect.isfunction(member) and \
+                hasattr(member, 'buildable') and member.buildable == True:
+            buildable_lookup[member.identifier] = member
+
+    if detector_name in buildable_lookup:
+        return buildable_lookup[detector_name]()
diff --git a/gpu.py b/gpu.py
index e5fd5f1..11a83ee 100644
--- a/gpu.py
+++ b/gpu.py
@@ -2,6 +2,7 @@ import numpy as np
 import numpy.ma as ma
 from pycuda.tools import make_default_context
 from pycuda.compiler import SourceModule
+from pycuda.characterize import sizeof
 import pycuda.driver as cuda
 from pycuda import gpuarray
 from copy import copy
@@ -27,11 +28,39 @@ def format_array(name, array):
     return '%-15s %6s %6s' % \
         (name, format_size(len(array)), format_size(array.nbytes))
 
+class CUDAFuncs(object):
+    '''simple container class for GPU functions as attributes'''
+    def __init__(self, cuda_module, func_names):
+        '''Extract __global__ functions listed in func_names from
+        the PyCUDA module object.  The functions are assigned
+        to attributes of this object with the same name.'''
+        for name in func_names:
+            setattr(self, name, cuda_module.get_function(name))
+
 class GPU(object):
-    def __init__(self):
-        self.context = make_default_context()
+    def __init__(self, device_id=None):
+        '''Initialize a GPU context on the specified device.  If device_id==None,
+        the default device is used.'''
+
+        if device_id == None:
+            self.context = make_default_context()
+        else:
+            device = pycuda.driver.Device(device_id)
+            self.context = device.make_context()
         print 'device %s' % self.context.get_device().name()
         self.module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
+        self.geo_funcs = CUDAFuncs(self.module, ['set_wavelength_range', 'set_material',
+                                                 'set_surface', 
+                                                 'set_global_mesh_variables'])
+                                                 
+                                                 
+        self.prop_funcs = CUDAFuncs(self.module, ['init_rng', 'propagate'])
+        self.nblocks = 64
+        self.max_nthreads = 200000
+        self.daq_module = SourceModule(src.daq, options=['-I' + src.dir], no_extern_c=True, cache_dir=False)
+        self.daq_funcs = CUDAFuncs(self.daq_module, 
+                              ['reset_earliest_time_int', 'run_daq',
+                               'convert_sortable_int_to_float'])
 
     def print_device_usage(self):
         print 'device usage:'
@@ -47,10 +76,8 @@ class GPU(object):
         if not hasattr(geometry, 'mesh'):
             geometry.build(bits=8)
 
-        set_wavelength_range = self.module.get_function('set_wavelength_range')
-        set_wavelength_range(np.float32(standard_wavelengths[0]), np.float32(standard_wavelengths[-1]), np.float32(standard_wavelengths[1]-standard_wavelengths[0]), np.uint32(standard_wavelengths.size), block=(1,1,1), grid=(1,1))
+        self.geo_funcs.set_wavelength_range(np.float32(standard_wavelengths[0]), np.float32(standard_wavelengths[-1]), np.float32(standard_wavelengths[1]-standard_wavelengths[0]), np.uint32(standard_wavelengths.size), block=(1,1,1), grid=(1,1))
 
-        set_material = self.module.get_function('set_material')
         self.materials = []
         for i in range(len(geometry.unique_materials)):
             material = copy(geometry.unique_materials[i])
@@ -62,11 +89,10 @@ class GPU(object):
             material.absorption_length_gpu = gpuarray.to_gpu(np.interp(standard_wavelengths, material.absorption_length[:,0], material.absorption_length[:,1]).astype(np.float32))
             material.scattering_length_gpu = gpuarray.to_gpu(np.interp(standard_wavelengths, material.scattering_length[:,0], material.scattering_length[:,1]).astype(np.float32))
 
-            set_material(np.int32(i), material.refractive_index_gpu, material.absorption_length_gpu, material.scattering_length_gpu, block=(1,1,1), grid=(1,1))
+            self.geo_funcs.set_material(np.int32(i), material.refractive_index_gpu, material.absorption_length_gpu, material.scattering_length_gpu, block=(1,1,1), grid=(1,1))
 
             self.materials.append(material)
 
-        set_surface = self.module.get_function('set_surface')
         self.surfaces = []
         for i in range(len(geometry.unique_surfaces)):
             surface = copy(geometry.unique_surfaces[i])
@@ -79,7 +105,7 @@ class GPU(object):
             surface.reflect_diffuse_gpu = gpuarray.to_gpu(np.interp(standard_wavelengths, surface.reflect_diffuse[:,0], surface.reflect_diffuse[:,1]).astype(np.float32))
             surface.reflect_specular_gpu = gpuarray.to_gpu(np.interp(standard_wavelengths, surface.reflect_specular[:,0], surface.reflect_specular[:,1]).astype(np.float32))
 
-            set_surface(np.int32(i), surface.detect_gpu, surface.absorb_gpu, surface.reflect_diffuse_gpu, surface.reflect_specular_gpu, block=(1,1,1), grid=(1,1))
+            self.geo_funcs.set_surface(np.int32(i), surface.detect_gpu, surface.absorb_gpu, surface.reflect_diffuse_gpu, surface.reflect_specular_gpu, block=(1,1,1), grid=(1,1))
 
             self.surfaces.append(surface)
 
@@ -110,8 +136,7 @@ class GPU(object):
         self.node_length_gpu = gpuarray.to_gpu(geometry.node_length.astype(np.uint32))
         self.solid_id_map_gpu = gpuarray.to_gpu(geometry.solid_id.astype(np.uint32))
 
-        set_global_mesh_variables = self.module.get_function('set_global_mesh_variables')
-        set_global_mesh_variables(self.triangles_gpu, self.vertices_gpu, self.colors_gpu, np.uint32(geometry.node_map.size-1), np.uint32(geometry.first_node), block=(1,1,1), grid=(1,1))
+        self.geo_funcs.set_global_mesh_variables(self.triangles_gpu, self.vertices_gpu, self.colors_gpu, np.uint32(geometry.node_map.size-1), np.uint32(geometry.first_node), block=(1,1,1), grid=(1,1))
 
         self.lower_bounds_tex = self.module.get_texref('lower_bounds')
         self.upper_bounds_tex = self.module.get_texref('upper_bounds')
@@ -138,3 +163,118 @@ class GPU(object):
             triangle_colors[solid_id_map == i] = color
 
         self.colors_gpu.set(triangle_colors)
+
+    def setup_propagate(self, seed=1):
+        self.rng_states_gpu = cuda.mem_alloc(self.max_nthreads*sizeof('curandStateXORWOW', '#include <curand_kernel.h>'))
+        self.prop_funcs.init_rng(np.int32(self.max_nthreads), self.rng_states_gpu, 
+                                 np.uint64(seed), np.uint64(0), 
+                                 block=(self.nblocks,1,1), 
+                                 grid=(self.max_nthreads//self.nblocks+1,1))
+        #self.context.synchronize()
+
+    def load_photons(self, pos, dir, pol, wavelength, t0, 
+                     histories=None, last_hit_triangles=None):
+        '''Load N photons onto the GPU.
+
+           pos: numpy.array(shape=(N, 3)) of photon starting positions (meters)
+           dir: numpy.array(shape=(N, 3)) of photon starting directions (unit vectors)
+           pol: numpy.array(shape=(N, 3)) of photon polarization directions (unit vectors)
+           wavelength: numpy.array(shape=N) of photon wavelengths (nm)
+           t0: numpy.array(shape=N) of photon start times (s)
+           
+           Optional args will be loaded with defaults on GPU if not set:
+           histories: Bitmask of interactions that have occurred over history of photon
+           last_hit_triangles: The triangle ID number that the photon last interacted with,
+                               if any.  -1 if no triangle was hit in the last step
+        '''
+        self.nphotons = len(pos)
+        assert self.nphotons < self.max_nthreads
+        assert len(dir) == self.nphotons
+        assert len(pol) == self.nphotons
+        assert len(wavelength) == self.nphotons
+        assert len(t0) == self.nphotons
+
+        self.positions_gpu = gpuarray.to_gpu(pos.astype(np.float32).view(gpuarray.vec.float3))
+        self.directions_gpu = gpuarray.to_gpu(dir.astype(np.float32).view(gpuarray.vec.float3))
+        self.polarizations_gpu = gpuarray.to_gpu(pol.astype(np.float32).view(gpuarray.vec.float3))
+        self.wavelengths_gpu = gpuarray.to_gpu(wavelength.astype(np.float32))
+        self.times_gpu = gpuarray.to_gpu(t0.astype(np.float32))
+
+        if histories:
+            self.histories_gpu = gpuarray.to_gpu(histories.astype(np.uint32))
+        else:
+            self.histories_gpu = gpuarray.zeros(self.nphotons, dtype=np.uint32)
+
+        if last_hit_triangles:
+            self.last_hit_triangles_gpu = gpuarray.to_gpu(last_hit_triangles.astype(np.int32))
+        else:
+            self.last_hit_triangles_gpu = gpuarray.empty(self.nphotons, dtype=np.int32)
+            self.last_hit_triangles_gpu.fill(-1)
+
+        self.per_photon_kernel_config = {'block': (self.nblocks,1,1), 
+                                         'grid': (self.nphotons//self.nblocks+1,1)}
+
+    def propagate(self, max_steps=10):
+        '''Propagate photons on GPU to termination or max_steps, whichever comes first.
+
+        May be called repeatedly without reloading photon information if single-stepping
+        through photon history.
+        '''
+        self.prop_funcs.propagate(np.int32(self.nphotons), 
+                                 self.rng_states_gpu, self.positions_gpu, 
+                                 self.directions_gpu, 
+                                 self.wavelengths_gpu, self.polarizations_gpu, 
+                                 self.times_gpu, self.histories_gpu, 
+                                 self.last_hit_triangles_gpu, 
+                                 np.int32(max_steps), **self.per_photon_kernel_config)
+        #self.context.synchronize()
+
+    def get_photons(self):
+        '''Returns a dictionary of current photon state information.
+
+        Contents of dictionary have the same names as the parameters to load_photons().
+        '''
+        return { 'pos' : self.positions_gpu.get().view(np.float32).reshape(self.positions_gpu.size, 3),
+                 'dir' : self.directions_gpu.get().view(np.float32).reshape(self.directions_gpu.size, 3),
+                 'pol' : self.polarizations_gpu.get().view(np.float32).reshape(self.polarization_gpu.size, 3),
+                 'wavelength' : self.wavelengths_gpu.get(),
+                 't0' : self.times_gpu.get(),
+                 'histories' : self.histories_gpu.get(),
+                 'last_hit_triangles' : self.last_hit_triangles_gpu.get()}
+    
+    def setup_daq(self, max_pmt_id, pmt_rms=1.2e-9):
+        self.earliest_time_gpu = gpuarray.GPUArray(shape=(max_pmt_id+1,), dtype=np.float32)
+        self.earliest_time_int_gpu = gpuarray.GPUArray(shape=self.earliest_time_gpu.shape, 
+                                                       dtype=np.uint32)
+        self.daq_pmt_rms = pmt_rms
+
+    def run_daq(self):
+        self.daq_funcs.reset_earliest_time_int(np.float32(1e9), 
+                                               np.int32(len(self.earliest_time_int_gpu)), 
+                                               self.earliest_time_int_gpu, 
+                                               block=(self.nblocks,1,1), 
+                                               grid=(len(self.earliest_time_int_gpu)//self.nblocks+1,1))
+        #self.context.synchronize()
+        self.daq_funcs.run_daq(self.rng_states_gpu, np.uint32(0x1 << 2), 
+                               np.float32(self.daq_pmt_rms), 
+                               np.int32(len(self.positions_gpu)), 
+                               self.times_gpu, self.histories_gpu,
+                               self.last_hit_triangles_gpu, 
+                               self.solid_id_map_gpu, 
+                               np.int32(len(self.earliest_time_int_gpu)), 
+                               self.earliest_time_int_gpu, 
+                               block=(self.nblocks,1,1), grid=(self.nphotons//self.nblocks+1,1))
+        #self.context.synchronize()
+        self.daq_funcs.convert_sortable_int_to_float(np.int32(len(self.earliest_time_int_gpu)), 
+                                                     self.earliest_time_int_gpu, 
+                                                     self.earliest_time_gpu, 
+                                                     block=(self.nblocks,1,1), 
+                                                     grid=(len(self.earliest_time_int_gpu)//self.nblocks+1,1))
+        #self.context.synchronize()
+        
+
+    def get_hits(self):
+        return self.earliest_time_gpu.get()
+
+    def shutdown(self):
+        self.context.pop()
diff --git a/io/__init__.py b/io/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/io/__init__.py
diff --git a/io/root.C b/io/root.C
new file mode 100644
index 0000000..dfcb3c0
--- /dev/null
+++ b/io/root.C
@@ -0,0 +1,78 @@
+#include <TVector3.h>
+#include <vector>
+#include <TTree.h>
+#include <string>
+
+struct Photon {
+  double t;
+  TVector3 pos;
+  TVector3 dir;
+  TVector3 pol;
+  double wavelength; // nm
+  int state;
+  int last_triangle;
+  int last_solid;
+};
+
+struct MC {
+  std::string particle;
+  TVector3 gen_pos;
+  TVector3 gen_dir;
+  double gen_ke;
+
+  std::vector<Photon> photon;
+
+  void clear() {
+    particle = "none";
+    photon.clear();
+  }
+};
+
+struct Channel {
+  Channel() : channel_id(-1), t(-9999.0), q(-9999.0) { };
+  int channel_id;
+  double t;
+  double q;
+};
+  
+struct Event {
+  int event_id;
+  MC mc;
+  int nhit;
+  std::vector<Channel> channel;
+
+  void clear() {
+    mc.clear();
+    channel.clear();
+  }
+};
+
+void fill_event(TTree *T,
+		Event *ev, int event_id, double *gen_pos, int nchannels, 
+		float *channel_times)
+{
+  ev->clear();
+  ev->event_id = event_id;
+  MC &mc = ev->mc;
+  mc.gen_pos.SetXYZ(gen_pos[0], gen_pos[1], gen_pos[2]);
+  ev->nhit = 0;
+  Channel ch;
+  for (int i=0; i < nchannels; i++) {
+    if (channel_times[i] < 1e8) {
+      ev->nhit++;
+      ch.channel_id = i;
+      ch.t = channel_times[i];
+      ch.q = 1.0;
+      ev->channel.push_back(ch);
+    }
+  }
+  T->Fill();
+}
+
+
+#ifdef __MAKECINT__
+#pragma link C++ class vector<Photon>;
+#pragma link C++ class vector<Channel>;
+#endif
+
+
diff --git a/io/root.py b/io/root.py
new file mode 100644
index 0000000..d970d1d
--- /dev/null
+++ b/io/root.py
@@ -0,0 +1,18 @@
+import ROOT
+import os.path
+
+ROOT.gROOT.ProcessLine('.L '+os.path.join(os.path.dirname(__file__), 'root.C+g'))
+
+from ROOT import Event
+
+fill_event = ROOT.fill_event
+
+def make_tree(name, desc=''):
+    '''Create a ROOT tree for holding event information.
+
+    Returns tuple of Event object for filling and TTree.
+    '''
+    tree = ROOT.TTree(name, desc)
+    ev = ROOT.Event()
+    tree.Branch('ev', ev)
+    return ev, tree
diff --git a/sim.py b/sim.py
new file mode 100755
index 0000000..de219c8
--- /dev/null
+++ b/sim.py
@@ -0,0 +1,132 @@
+#!/usr/bin/env python
+import sys
+import optparse
+import time
+
+import detectors
+import optics
+import gpu
+import g4gen
+from io import root
+import numpy as np
+import math
+import ROOT
+
+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
+       # device, so we call the default hook
+       sys.__excepthook__(type, value, tb)
+    else:
+       import traceback, pdb
+       # we are NOT in interactive mode, print the exception...
+       traceback.print_exception(type, value, tb)
+       print
+       # ...then start the debugger in post-mortem mode.
+       pdb.pm()
+
+sys.excepthook = info
+
+if __name__ == '__main__':
+
+    parser = optparse.OptionParser('%prog')
+    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('--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-')
+    parser.add_option('--energy', type='float', dest='energy', default=100.0)
+    parser.add_option('--pos', type='string', dest='pos', default='(0,0,0)')
+    parser.add_option('--dir', type='string', dest='pos', default='(1,0,0)')
+
+    options, args = parser.parse_args()
+    if len(args) != 1:
+        print 'Must specify output filename!'
+        sys.exit(1)
+    else:
+        output_filename = args[0]
+
+    if options.nevents <= 0:
+        print '--nevents must be greater than 0!'
+        sys.exit(1)
+
+    position = np.array(eval(options.pos), dtype=float)
+    direction = np.array(eval(options.pos), dtype=float)
+    detector = detectors.find(options.detector)
+
+    print >>sys.stderr, 'Creating detector...'
+    detector.build(bits=options.nbits)
+
+    print >>sys.stderr, 'Initializing generator...'
+    detector_material = optics.water
+    generator = g4gen.G4Generator(detector_material)
+    print >>sys.stderr, 'WARNING: ASSUMING DETECTOR IS WATER!!'
+
+    print >>sys.stderr, 'Initializing GPU...'
+    gpu_worker = gpu.GPU(options.device)
+    
+    print >>sys.stderr, 'Loading detector onto GPU...'
+    gpu_worker.load_geometry(detector)
+    
+    print >>sys.stderr, 'Initializing random numbers generators...'
+    gpu_worker.setup_propagate()
+    gpu_worker.setup_daq(max(detector.pmtids))
+
+    # Create output file
+    f = ROOT.TFile(output_filename, 'RECREATE')
+    ev, T = root.make_tree('T')
+
+    print >>sys.stderr, 'Starting simulation...'
+    start_sim = time.time()
+    #### Do the generation and writing in this offset order to ensure
+    #### that propagation on the GPU overlaps with CPU work to create 
+    #### and save photons. 
+    #### WARNING: THIS MAKES THE CODE LOOK A LITTLE CRAZY.  I AM SORRY
+    
+    # Do first event
+    photons = generator.generate_photons(particle_name=options.particle, 
+                                         total_energy=options.energy,
+                                         position=position,
+                                         direction=direction)
+    nphotons = len(photons['pos'])
+    gpu_worker.load_photons(**photons)
+    gpu_worker.propagate()
+    gpu_worker.run_daq()
+
+    for i in xrange(1, options.nevents-1):
+        # photons for next event while previous event propagates on GPU
+        photons = generator.generate_photons(particle_name=options.particle, 
+                                             total_energy=options.energy,
+                                             position=position,
+                                             direction=direction)
+        nphotons += len(photons['pos'])
+
+        # this will stop and wait for event to finish
+        hit_times = gpu_worker.get_hits()
+
+        # turn around next event
+        gpu_worker.load_photons(**photons)
+        gpu_worker.propagate()
+        gpu_worker.run_daq()
+
+        # write out this event
+        root.fill_event(T, ev, i-1, position, len(hit_times), hit_times)
+        
+        if i % 10 == 0:
+            print >>sys.stderr, "\rEvent:", i,
+
+    # Get results for last event and write it out
+    hit_times = gpu_worker.get_hits()
+    root.fill_event(T, ev, options.nevents - 1, position, len(hit_times), hit_times)
+
+    end_sim = time.time()
+    print >>sys.stderr, "\rEvent:", options.nevents - 1
+
+    T.Write()
+    f.Close()
+
+    print >>sys.stderr, 'Done. %1.1f events/sec, %1.0f photons/sec.' % (options.nevents/(end_sim - start_sim), nphotons/(end_sim - start_sim))
+
+    gpu_worker.shutdown()
diff --git a/spnav.py b/spnav.py
new file mode 100644
index 0000000..9079114
--- /dev/null
+++ b/spnav.py
@@ -0,0 +1,80 @@
+from ctypes import cdll, c_int, c_uint, c_void_p, byref, Structure, Union
+
+libspnav = cdll.LoadLibrary('libspnav.so')
+
+SPNAV_EVENT_ANY = 0
+SPNAV_EVENT_MOTION = 1
+SPNAV_EVENT_BUTTON = 2
+
+class spnav_event_motion(Structure):
+    _fields_ = [('type', c_int),
+                ('x', c_int),
+                ('y', c_int),
+                ('z', c_int),
+                ('rx', c_int),
+                ('ry', c_int),
+                ('rz', c_int),
+                ('period', c_uint),
+                ('data', c_void_p)]
+
+class spnav_event_button(Structure):
+    _fields_ = [('type', c_int),
+                ('press', c_int),
+                ('bnum', c_int)]
+
+class spnav_event(Union):
+    _fields_ = [('type', c_int),
+                ('motion', spnav_event_motion),
+                ('button', spnav_event_button) ]
+
+    def __str__(self):
+        if self.type == SPNAV_EVENT_ANY:
+            return 'SPNAV_EVENT_ANY'
+        elif self.type == SPNAV_EVENT_MOTION:
+            m = self.motion
+            return 'SPNAV_EVENT_MOTION t(%d,%d,%d) r(%d,%d,%d)' % (m.x, m.y, m.z, m.rx, m.ry, m.rz)
+        elif self.type == SPNAV_EVENT_BUTTON:
+            if self.button.press:
+                state = 'down'
+            else:
+                state = 'up'
+            return 'SPNAV_EVENT_BUTTON %d %s' % (self.button.bnum, state)
+    
+def spnav_open():
+    if libspnav.spnav_open() == -1:
+        raise Exception('failed to connect to the space navigator daemon')
+
+def spnav_close():
+    libspnav.spnav_close()
+
+def spnav_fd():
+    return libspnav.spnav_fd()
+
+def spnav_wait_event():
+    ev = spnav_event()
+    ret = libspnav.spnav_wait_event(byref(ev))
+    if ret:
+        return ev
+    else:
+        raise Exception('non-zero return code from spnav_wait_event()')
+
+def spnav_poll_event():
+    ev = spnav_event()
+    ret = libspnav.spnav_poll_event(byref(ev))
+    if ret == 0:
+        return None
+    else:
+        return ev
+
+def spnav_remove_events(event_type):
+    return libspnav.spnav_remove_events(event_type)
+
+if __name__ == '__main__':
+    spnav_open()
+    try:
+        while True:
+            ev = spnav_wait_event()
+            print ev
+    finally:
+        spnav_close()
+
diff --git a/src/daq.cu b/src/daq.cu
index 07e2d5e..2b5f9b4 100644
--- a/src/daq.cu
+++ b/src/daq.cu
@@ -48,12 +48,13 @@ __global__ void run_daq(curandState *s, unsigned int detection_state,
 		{
 			int solid_id = solid_map[triangle_id];
 			int history = photon_histories[id];
-			float time = photon_times[id] + curand_normal(&rng) * time_rms;
-			unsigned int time_int = float_to_sortable_int(time);
 
 			if (solid_id < nsolids && (history & detection_state)) 
 			{
-			    atomicMin(earliest_time_int + solid_id, time_int);
+			  float time = photon_times[id] + curand_normal(&rng) * time_rms;
+			  unsigned int time_int = float_to_sortable_int(time);
+			  
+			  atomicMin(earliest_time_int + solid_id, time_int);
 			}
 
 		}