A GPU-side "DAQ" implementation that identifies the first photon on each channel, which is presumed to trigger that channel. Major speed up in conversion of detection times to time PDFs.

2 files changed, 81 insertions, 0 deletions

diff --git a/src/__init__.py b/src/__init__.py
index d2958f1..865d612 100644
--- a/src/__init__.py
+++ b/src/__init__.py
@@ -3,3 +3,4 @@ import os
 dir = os.path.split(os.path.realpath(__file__))[0]
 
 kernel = open(dir + '/kernel.cu').read()
+daq = open(dir + '/daq.cu').read()
diff --git a/src/daq.cu b/src/daq.cu
new file mode 100644
index 0000000..c79401c
--- /dev/null
+++ b/src/daq.cu
@@ -0,0 +1,80 @@
+// -*-c++-*-
+#include <curand_kernel.h>
+
+__device__ unsigned int float_to_sortable_int(float f)
+{
+  return __float_as_int(f);
+  //int i = __float_as_int(f);
+  //unsigned int mask = -(int)(i >> 31) | 0x80000000;
+  //return i ^ mask;
+}
+
+__device__ float sortable_int_to_float(unsigned int i)
+{
+  return __int_as_float(i);
+  //unsigned int mask = ((i >> 31) - 1) | 0x80000000;
+  //return __int_as_float(i ^ mask);
+}
+
+
+__device__ curandState daq_rng_states[100000];
+
+extern "C" {
+
+  __global__ void init_daq_rng(int nthreads, 
+			  unsigned long long seed, unsigned long long offset)
+  {
+    int id = blockIdx.x*blockDim.x + threadIdx.x;
+    
+    if (id >= nthreads)
+      return;
+    
+    curand_init(seed, id, offset, daq_rng_states+id);
+  }
+
+  __global__ void reset_earliest_time_int(float maxtime, 
+				     int ntime_ints, unsigned int *time_ints)
+  {
+    int id = threadIdx.x + blockDim.x * blockIdx.x;
+    if (id < ntime_ints) {
+      unsigned int maxtime_int = float_to_sortable_int(maxtime);
+      time_ints[id] = maxtime_int;
+    }
+  }
+
+  __global__ void run_daq(int detection_state, float time_rms,
+		     int nphotons, float *photon_times, int *photon_states,
+		     int *last_hit_triangles, int *solid_map,
+		     int nsolids, unsigned int *earliest_time_int)
+  {
+    int id = threadIdx.x + blockDim.x * blockIdx.x;
+    
+    curandState_t rng = daq_rng_states[id];
+
+    if (id < nphotons) {
+      int triangle_id = last_hit_triangles[id];
+      
+      if (triangle_id > -1) {
+	int solid_id = solid_map[triangle_id];
+	int state = photon_states[id];
+	float time = photon_times[id];// + curand_normal(&rng) * time_rms;
+	unsigned int time_int = float_to_sortable_int(time);
+	if (solid_id < nsolids && state == detection_state)
+	  atomicMin(earliest_time_int + solid_id, time_int);
+      }
+    }
+
+    daq_rng_states[id] = rng;
+  }
+
+  __global__ void convert_sortable_int_to_float(int n, 
+						unsigned int *sortable_ints, 
+						float *float_output)
+  {
+    int id = threadIdx.x + blockDim.x * blockIdx.x;
+    
+    if (id < n)
+      float_output[id] = sortable_int_to_float(sortable_ints[id]);
+  }
+
+} // extern "C"