Create a Detector class to hold information about the PMTs in a

geometry, like the mapping from solid IDs to channels, and the time
and charge distributions.

Detector is a subclass of Geometry, so that a Detector can be used
wherever a Geometry is used.  Only code (like the DAQ stuff) that
needs to know how PMT solids map to channels should look for a
Detector object.

There is a corresponding GPUDetector class as well, with its own
device side struct to hold PMT channel information.

The GPU code now can sample an arbitrary time and charge PDF, but
on the host side, the only interface exposed right now creates
a Gaussian distribution.

3 files changed, 61 insertions, 13 deletions

diff --git a/chroma/cuda/daq.cu b/chroma/cuda/daq.cu
index f80fcaa..b995bfb 100644
--- a/chroma/cuda/daq.cu
+++ b/chroma/cuda/daq.cu
@@ -1,5 +1,6 @@
 // -*-c++-*-
-#include <curand_kernel.h>
+#include "detector.h"
+#include "random.h"
 
 __device__ unsigned int
 float_to_sortable_int(float f)
@@ -32,11 +33,13 @@ reset_earliest_time_int(float maxtime, int ntime_ints, unsigned int *time_ints)
 }
 
 __global__ void
-run_daq(curandState *s, unsigned int detection_state, float time_rms,
+run_daq(curandState *s, unsigned int detection_state,
 	int first_photon, int nphotons, float *photon_times,
 	unsigned int *photon_histories, int *last_hit_triangles,
-	int *solid_map, int nsolids, unsigned int *earliest_time_int,
-	unsigned int *channel_q, unsigned int *channel_histories)
+	int *solid_map,
+	Detector *detector,
+	unsigned int *earliest_time_int,
+	unsigned int *channel_q_int, unsigned int *channel_histories)
 {
 
     int id = threadIdx.x + blockDim.x * blockIdx.x;
@@ -49,14 +52,22 @@ run_daq(curandState *s, unsigned int detection_state, float time_rms,
 	if (triangle_id > -1) {
 	    int solid_id = solid_map[triangle_id];
 	    unsigned int history = photon_histories[photon_id];
+	    int channel_index = detector->solid_id_to_channel_index[solid_id];
 
-	    if (solid_id < nsolids && (history & detection_state)) {
-		float time = photon_times[photon_id] + curand_normal(&rng) * time_rms;
+	    if (channel_index >= 0 && (history & detection_state)) {
+		float time = photon_times[photon_id] + 
+		    sample_cdf(&rng, detector->time_cdf_len,
+			       detector->time_cdf_x, detector->time_cdf_y);
 		unsigned int time_int = float_to_sortable_int(time);
 			  
-		atomicMin(earliest_time_int + solid_id, time_int);
-		atomicAdd(channel_q + solid_id, 1);
-		atomicOr(channel_histories + solid_id, history);
+		float charge = sample_cdf(&rng, detector->charge_cdf_len,
+					  detector->charge_cdf_x,
+					  detector->charge_cdf_y);
+		unsigned int charge_int = roundf(charge / detector->charge_unit);
+
+		atomicMin(earliest_time_int + channel_index, time_int);
+		atomicAdd(channel_q_int + channel_index, charge_int);
+		atomicOr(channel_histories + channel_index, history);
 	    }
 
 	}
@@ -78,4 +89,16 @@ convert_sortable_int_to_float(int n, unsigned int *sortable_ints,
 }
 
 
+__global__ void
+convert_charge_int_to_float(Detector *detector, 
+			    unsigned int *charge_int,
+			    float *charge_float)
+{
+    int id = threadIdx.x + blockDim.x * blockIdx.x;
+	
+    if (id < detector->nchannels)
+	charge_float[id] = charge_int[id] * detector->charge_unit;
+}
+
+
 } // extern "C"
diff --git a/chroma/cuda/detector.h b/chroma/cuda/detector.h
new file mode 100644
index 0000000..16bd164
--- /dev/null
+++ b/chroma/cuda/detector.h
@@ -0,0 +1,25 @@
+#ifndef __DETECTOR_H__
+#define __DETECTOR_H__
+
+struct Detector
+{
+    // Order in decreasing size to avoid alignment problems
+    int *solid_id_to_channel_index;
+
+    float *time_cdf_x;
+    float *time_cdf_y;
+
+    float *charge_cdf_x;
+    float *charge_cdf_y;
+
+    int nchannels;
+    int time_cdf_len;
+    int charge_cdf_len;
+    float charge_unit; 
+    // Convert charges to/from quantized integers with
+    // q_int = (int) roundf(q / charge_unit )
+    // q = q_int * charge_unit
+};
+
+
+#endif // __DETECTOR_H__
diff --git a/chroma/cuda/pdf.cu b/chroma/cuda/pdf.cu
index 9f547d0..0d82e3a 100644
--- a/chroma/cuda/pdf.cu
+++ b/chroma/cuda/pdf.cu
@@ -5,7 +5,7 @@ extern "C"
 {
 
 __global__ void
-bin_hits(int nchannels, unsigned int *channel_q, float *channel_time,
+bin_hits(int nchannels, float *channel_q, float *channel_time,
 	 unsigned int *hitcount, int tbins, float tmin, float tmax, int qbins,
 	 float qmin, float qmax, unsigned int *pdf)
 {
@@ -32,7 +32,7 @@ bin_hits(int nchannels, unsigned int *channel_q, float *channel_time,
 __global__ void
 accumulate_pdf_eval(int time_only, int nchannels, unsigned int *event_hit,
 		    float *event_time, float *event_charge, float *mc_time,
-		    unsigned int *mc_charge, // quirk of DAQ!
+		    float *mc_charge, // quirk of DAQ!
 		    unsigned int *hitcount, unsigned int *bincount,
 		    float min_twidth, float tmin, float tmax,
 		    float min_qwidth, float qmin, float qmax,
@@ -126,7 +126,7 @@ accumulate_pdf_eval(int time_only, int nchannels, unsigned int *event_hit,
 __global__ void
 accumulate_moments(int time_only, int nchannels,
 		   float *mc_time,
-		   unsigned int *mc_charge, // quirk of DAQ!
+		   float *mc_charge,
 		   float tmin, float tmax,
 		   float qmin, float qmax,
 		   unsigned int *mom0,
@@ -174,7 +174,7 @@ static const float rootPiBy2 = 1.2533141373155001f; // sqrt(M_PI/2)
 __global__ void
 accumulate_kernel_eval(int time_only, int nchannels, unsigned int *event_hit,
 		       float *event_time, float *event_charge, float *mc_time,
-		       unsigned int *mc_charge, // quirk of DAQ!
+		       float *mc_charge,
 		       float tmin, float tmax,
 		       float qmin, float qmax,
 		       float *inv_time_bandwidths,