Store a photon history for each hit channel. If multiple photons hit the

same channel, their history bits are OR'ed together.

4 files changed, 24 insertions, 12 deletions

diff --git a/gpu.py b/gpu.py
index a21bc38..7b08473 100644
--- a/gpu.py
+++ b/gpu.py
@@ -278,6 +278,7 @@ class GPU(object):
         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.channel_history_gpu = gpuarray.zeros_like(self.earliest_time_int_gpu)
         self.daq_pmt_rms = pmt_rms
 
     def run_daq(self):
@@ -286,6 +287,8 @@ class GPU(object):
                                                self.earliest_time_int_gpu, 
                                                block=(self.nthread_per_block,1,1), 
                                                grid=(len(self.earliest_time_int_gpu)//self.nthread_per_block+1,1))
+        self.channel_history_gpu.fill(0)
+
         #self.context.synchronize()
         for first_photon, photons_this_round, blocks in chunk_iterator(self.nphotons, self.nthread_per_block, self.max_blocks):
             self.daq_funcs.run_daq(self.rng_states_gpu, np.uint32(0x1 << 2), 
@@ -296,7 +299,8 @@ class GPU(object):
                                    self.last_hit_triangles_gpu, 
                                    self.solid_id_map_gpu, 
                                    np.int32(len(self.earliest_time_int_gpu)), 
-                                   self.earliest_time_int_gpu, 
+                                   self.earliest_time_int_gpu,
+                                   self.channel_history_gpu,
                                    block=(self.nthread_per_block,1,1), grid=(blocks,1))
         #self.context.synchronize()
         self.daq_funcs.convert_sortable_int_to_float(np.int32(len(self.earliest_time_int_gpu)), 
@@ -308,7 +312,9 @@ class GPU(object):
         
 
     def get_hits(self):
-        return self.earliest_time_gpu.get()
+        return { 't': self.earliest_time_gpu.get(), 
+                 'q': np.ones(shape=len(self.earliest_time_int_gpu), dtype=np.float32),
+                 'history': self.channel_history_gpu.get()}
 
     def __del__(self):
         self.context.pop()
diff --git a/io/root.C b/io/root.C
index 43f4e73..096811a 100644
--- a/io/root.C
+++ b/io/root.C
@@ -9,7 +9,7 @@ struct Photon {
   TVector3 dir;
   TVector3 pol;
   double wavelength; // nm
-  int history;
+  unsigned int history;
   int last_hit_triangle;
 };
 
@@ -29,6 +29,7 @@ struct Channel {
   int channel_id;
   double t;
   double q;
+  unsigned int mc_history;
 };
   
 struct Event {
@@ -67,17 +68,19 @@ void fill_photons(Event *ev, bool start,
 }
 		  
 
-void fill_hits(Event *ev, unsigned int nchannels, float *channel_times)
+void fill_hits(Event *ev, unsigned int nchannels, float *t,
+	       float *q, unsigned int *history)
 {
   ev->channel.resize(0);
   ev->nhit = 0;
   Channel ch;
   for (unsigned int i=0; i < nchannels; i++) {
-    if (channel_times[i] < 1e8) {
+    if (t[i] < 1e8) {
       ev->nhit++;
       ch.channel_id = i;
-      ch.t = channel_times[i];
-      ch.q = 1.0;
+      ch.t = t[i];
+      ch.q = q[i];
+      ch.mc_history = history[i];
       ev->channel.push_back(ch);
     }
   }
diff --git a/sim.py b/sim.py
index 87993ab..c8f9fdb 100755
--- a/sim.py
+++ b/sim.py
@@ -72,9 +72,10 @@ def write_event(T, ev, event_id, hits, photon_start=None, photon_stop=None):
                           photons['wavelength'], photons['t0'],
                           photons['histories'], photons['last_hit_triangles'])
 
-    root.fill_hits(ev, len(hits), hits)
+    root.fill_hits(ev, len(hits), hits['t'], hits['q'], hits['history'])
     T.Fill()
 
+#@profile
 def main():
     parser = optparse.OptionParser('%prog')
     parser.add_option('-b', type='int', dest='nbits', default=10)
@@ -159,7 +160,7 @@ def main():
         gpu_worker.load_photons(**photons)
         gpu_worker.propagate()
         gpu_worker.run_daq()
-        hit_times = gpu_worker.get_hits()
+        hits = gpu_worker.get_hits()
 
         if options.save_photon_start:
             photon_start = photons
@@ -170,7 +171,7 @@ def main():
         else:
             photon_stop = None
             
-        write_event(T, ev, i, hit_times, 
+        write_event(T, ev, i, hits,
                     photon_start=photon_start, photon_stop=photon_stop)
         if i % 10 == 0:
             print >>sys.stderr, "\rEvent:", i,
diff --git a/src/daq.cu b/src/daq.cu
index 7c5e6a5..9824d7d 100644
--- a/src/daq.cu
+++ b/src/daq.cu
@@ -34,7 +34,8 @@ __global__ void run_daq(curandState *s, unsigned int detection_state,
 			int nphotons, float *photon_times,
 			unsigned int *photon_histories,
 			int *last_hit_triangles, int *solid_map,
-			int nsolids, unsigned int *earliest_time_int)
+			int nsolids, unsigned int *earliest_time_int,
+			unsigned int *channel_histories)
 {
 
 	int id = threadIdx.x + blockDim.x * blockIdx.x;
@@ -48,7 +49,7 @@ __global__ void run_daq(curandState *s, unsigned int detection_state,
 		if (triangle_id > -1)
 		{
 			int solid_id = solid_map[triangle_id];
-			int history = photon_histories[photon_id];
+			unsigned int history = photon_histories[photon_id];
 
 			if (solid_id < nsolids && (history & detection_state)) 
 			{
@@ -56,6 +57,7 @@ __global__ void run_daq(curandState *s, unsigned int detection_state,
 			  unsigned int time_int = float_to_sortable_int(time);
 			  
 			  atomicMin(earliest_time_int + solid_id, time_int);
+			  atomicOr(channel_histories + solid_id, history);
 			}
 
 		}