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.

15 files changed, 379 insertions, 67 deletions

diff --git a/chroma/benchmark.py b/chroma/benchmark.py
index af422dd..44cf44b 100755
--- a/chroma/benchmark.py
+++ b/chroma/benchmark.py
@@ -64,7 +64,7 @@ def load_photons(number=100, nphotons=500000):
 
     return nphotons/ufloat((np.mean(run_times),np.std(run_times)))
 
-def propagate(gpu_geometry, number=10, nphotons=500000, nthreads_per_block=64,
+def propagate(gpu_detector, number=10, nphotons=500000, nthreads_per_block=64,
               max_blocks=1024):
     "Returns the average number of photons propagated on the GPU per second."
     rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks)
@@ -79,7 +79,7 @@ def propagate(gpu_geometry, number=10, nphotons=500000, nthreads_per_block=64,
         gpu_photons = gpu.GPUPhotons(photons)
 
         t0 = time.time()
-        gpu_photons.propagate(gpu_geometry, rng_states, nthreads_per_block,
+        gpu_photons.propagate(gpu_detector, rng_states, nthreads_per_block,
                               max_blocks)
         cuda.Context.get_current().synchronize()
         elapsed = time.time() - t0
@@ -91,7 +91,7 @@ def propagate(gpu_geometry, number=10, nphotons=500000, nthreads_per_block=64,
     return nphotons/ufloat((np.mean(run_times),np.std(run_times)))
 
 @tools.profile_if_possible
-def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
+def pdf(gpu_detector, npdfs=10, nevents=100, nreps=16, ndaq=1,
         nthreads_per_block=64, max_blocks=1024):
     """
     Returns the average number of 100 MeV events per second that can be
@@ -99,9 +99,7 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
 
     Args:
         - gpu_instance, chroma.gpu.GPU
-            The GPU instance passed to the GPUGeometry constructor.
-        - max_pmt_id, int
-            The channel number of the highest PMT
+            The GPU instance passed to the GPUDetector constructor.
         - npdfs, int
             The number of pdf generations to average.
         - nevents, int
@@ -114,9 +112,9 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
     """
     rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks)
 
-    gpu_daq = gpu.GPUDaq(gpu_geometry, max_pmt_id)
+    gpu_daq = gpu.GPUDaq(gpu_detector)
     gpu_pdf = gpu.GPUPDF()
-    gpu_pdf.setup_pdf(max_pmt_id, 100, (-0.5, 999.5), 10, (-0.5, 9.5))
+    gpu_pdf.setup_pdf(gpu_detector.nchannels, 100, (-0.5, 999.5), 10, (-0.5, 9.5))
 
     run_times = []
     for i in tools.progress(range(npdfs)):
@@ -130,7 +128,7 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
         for ev in g4generator.generate_events(vertex_iter):
             gpu_photons = gpu.GPUPhotons(ev.photons_beg, ncopies=nreps)
 
-            gpu_photons.propagate(gpu_geometry, rng_states,
+            gpu_photons.propagate(gpu_detector, rng_states,
                                   nthreads_per_block, max_blocks)
             for gpu_photon_slice in gpu_photons.iterate_copies():
                 for idaq in xrange(ndaq):
@@ -148,7 +146,7 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
 
     return nevents*nreps*ndaq/ufloat((np.mean(run_times),np.std(run_times)))
 
-def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
+def pdf_eval(gpu_detector, npdfs=10, nevents=25, nreps=16, ndaq=128,
              nthreads_per_block=64, max_blocks=1024):
     """
     Returns the average number of 100 MeV events that can be
@@ -156,9 +154,7 @@ def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
 
     Args:
         - gpu_instance, chroma.gpu.GPU
-            The GPU instance passed to the GPUGeometry constructor.
-        - max_pmt_id, int
-            The channel number of the highest PMT
+            The GPU instance passed to the GPUDetector constructor.
         - npdfs, int
             The number of pdf generations to average.
         - nevents, int
@@ -177,9 +173,9 @@ def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
                                                            1)).next()
     gpu_photons = gpu.GPUPhotons(data_ev.photons_beg)
 
-    gpu_photons.propagate(gpu_geometry, rng_states,
+    gpu_photons.propagate(gpu_detector, rng_states,
                           nthreads_per_block, max_blocks)
-    gpu_daq = gpu.GPUDaq(gpu_geometry, max_pmt_id)
+    gpu_daq = gpu.GPUDaq(gpu_detector)
     data_ev_channels = gpu_daq.acquire(gpu_photons, rng_states, nthreads_per_block, max_blocks).get()
     
     # Setup PDF evaluation
@@ -206,7 +202,7 @@ def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
         for ev in g4generator.generate_events(vertex_iter):
             gpu_photons = gpu.GPUPhotons(ev.photons_beg, ncopies=nreps)
 
-            gpu_photons.propagate(gpu_geometry, rng_states,
+            gpu_photons.propagate(gpu_detector, rng_states,
                                   nthreads_per_block, max_blocks)
             for gpu_photon_slice in gpu_photons.iterate_copies():
                 for idaq in xrange(ndaq):
@@ -237,11 +233,11 @@ if __name__ == '__main__':
     detector.build(bits=11)
 
     context = gpu.create_cuda_context()
-    gpu_geometry = gpu.GPUGeometry(detector)
+    gpu_detector = gpu.GPUDetector(detector)
     
     if 'ray' in tests:
         print '%s ray intersections/sec.' % \
-            tools.ufloat_to_str(intersect(gpu_geometry))
+            tools.ufloat_to_str(intersect(gpu_detector))
         # run garbage collection since there is a reference loop
         # in the GPUArray class.
         gc.collect()
@@ -252,16 +248,16 @@ if __name__ == '__main__':
 
     if 'propagate' in tests:
         print '%s photons propagated/sec.' % \
-            tools.ufloat_to_str(propagate(gpu_geometry))
+            tools.ufloat_to_str(propagate(gpu_detector))
         gc.collect()
 
     if 'pdf' in tests:
         print '%s 100 MeV events histogrammed/s' % \
-            tools.ufloat_to_str(pdf(gpu_geometry, max(detector.pmtids)))
+            tools.ufloat_to_str(pdf(gpu_detector))
         gc.collect()
 
     if 'pdf_eval' in tests:
         print '%s 100 MeV events/s accumulated in PDF evaluation data structure (100 GEANT4 x 16 Chroma x 128 DAQ)' % \
-            tools.ufloat_to_str(pdf_eval(gpu_geometry, max(detector.pmtids)))
+            tools.ufloat_to_str(pdf_eval(gpu_detector))
 
     context.pop()
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,
diff --git a/chroma/demo/__init__.py b/chroma/demo/__init__.py
index 7df7df7..46f4529 100644
--- a/chroma/demo/__init__.py
+++ b/chroma/demo/__init__.py
@@ -6,13 +6,15 @@ import numpy.linalg
 
 from chroma.make import sphere
 from chroma.stl import mesh_from_stl
-from chroma.geometry import *
+from chroma.geometry import Solid
+from chroma.detector import Detector
 from chroma.transform import rotate, make_rotation_matrix, normalize
 
 from chroma.demo.pmt import build_8inch_pmt_with_lc
 from chroma.demo.optics import water, black_surface
 
 from chroma.demo.checkerboard import build_checkerboard_scene as checkerboard_scene
+from chroma.log import logger
 
 def spherical_spiral(radius, spacing):
     '''Returns iterator generating points on a spiral wrapping the
@@ -29,16 +31,13 @@ def spherical_spiral(radius, spacing):
 
 def detector(pmt_radius=14000.0, sphere_radius=14500.0, spiral_step=350.0):
     pmt = build_8inch_pmt_with_lc()
-    geo = Geometry(water)
+    geo = Detector(water)
 
     geo.add_solid(Solid(sphere(sphere_radius,nsteps=200), 
                         water, water, 
                         surface=black_surface,
                         color=0xBBFFFFFF))
 
-    geo.pmtids = []
-
-
     for position in spherical_spiral(pmt_radius, spiral_step):
         direction = -normalize(position)
 
@@ -49,10 +48,19 @@ def detector(pmt_radius=14000.0, sphere_radius=14500.0, spiral_step=350.0):
         rotation = make_rotation_matrix(angle, axis)
 
         # Place PMT (note that position is front face of PMT)
-        chroma_id = geo.add_solid(pmt, rotation, position)
-        geo.pmtids.append(chroma_id)
+        geo.add_pmt(pmt, rotation, position)
         
-    print 'Demo detector: %d PMTs' % len(geo.pmtids)
+    
+    time_rms = 1.5 # ns
+    charge_mean = 1.0
+    charge_rms = 0.1 # Don't I wish!
+    
+    geo.set_time_dist_gaussian(time_rms, -5 * time_rms, 5*time_rms)
+    geo.set_charge_dist_gaussian(charge_mean, charge_rms, 0.0, charge_mean + 5*charge_rms)
+
+    logger.info('Demo detector: %d PMTs' % geo.num_channels())
+    logger.info('               %1.1f ns time RMS' % time_rms)
+    logger.info('               %1.1f%% charge RMS' % (100.0*charge_rms/charge_mean))
     return geo
 
 def tiny():
diff --git a/chroma/detector.py b/chroma/detector.py
new file mode 100644
index 0000000..f3f2d3b
--- /dev/null
+++ b/chroma/detector.py
@@ -0,0 +1,129 @@
+import numpy as np
+
+from chroma.geometry import Geometry
+
+class Detector(Geometry):
+    '''A Detector is a subclass of Geometry that allows some Solids
+    to be marked as photon detectors, which we will suggestively call
+    "PMTs."  Each detector is imagined to be connected to an electronics
+    channel that records a hit time and charge.
+
+    Each PMT has two integers identifying it: a channel index and a
+    channel ID.  When all of the channels in the detector are stored
+    in a Numpy array, they will be stored in index order.  Channel
+    indices star from zero and have no gaps.  Channel ID numbers are
+    arbitrary integers that identify a PMT uniquely in a stable way.
+    They are written out to disk when using the Chroma ROOT format,
+    and are used when reading events back in to map channels back
+    into the correct array index.
+
+    For now, all the PMTs share a single set of time and charge
+    distributions.  In the future, this will be generalized to
+    allow per-channel distributions.
+    '''
+
+    def __init__(self, detector_material=None):
+        Geometry.__init__(self, detector_material=detector_material)
+
+        # Using numpy arrays here to allow for fancy indexing
+        self.solid_id_to_channel_index = np.zeros(0, dtype=np.int32)
+        self.channel_index_to_solid_id = np.zeros(0, dtype=np.int32)
+
+        self.channel_index_to_channel_id = np.zeros(0, dtype=np.int32)
+
+        # If the ID numbers are arbitrary, we can't treat them
+        # as array indices, so have to use a dictionary
+        self.channel_id_to_channel_index = {}
+
+        # zero time and unit charge distributions
+        self.time_cdf = (np.array([-0.01, 0.01]), np.array([0.0, 1.0]))
+        self.charge_cdf = (np.array([0.99, 1.00]), np.array([0.0, 1.0]))
+
+    def add_solid(self, solid, rotation=None, displacement=None):
+        """
+        Add the solid `solid` to the geometry. When building the final triangle
+        mesh, `solid` will be placed by rotating it with the rotation matrix
+        `rotation` and displacing it by the vector `displacement`.
+        """
+        solid_id = Geometry.add_solid(self, solid=solid, rotation=rotation, 
+                                      displacement=displacement)
+        self.solid_id_to_channel_index.resize(solid_id+1)
+        self.solid_id_to_channel_index[solid_id] = -1 # solid maps to no channel
+        return solid_id
+
+    def add_pmt(self, pmt, rotation=None, displacement=None, channel_id=None):
+        """Add the PMT `pmt` to the geometry. When building the final triangle
+        mesh, `solid` will be placed by rotating it with the rotation matrix
+        `rotation` and displacing it by the vector `displacement`, just like
+        add_solid().
+
+            `pmt``: instance of chroma.Solid
+                Solid representing a PMT.
+            `rotation`: numpy.matrix (3x3)
+                Rotation to apply to PMT mesh before displacement.  Defaults to
+                identity rotation.
+            `displacement`: numpy.ndarray (shape=3)
+                3-vector displacement to apply to PMT mesh after rotation.
+                Defaults to zero vector.
+            `channel_id`: int
+                Integer identifier for this PMT.  May be any integer, with no
+                requirement for consective numbering.  Defaults to None,
+                where the ID number will be set to the generated channel index.
+                The channel_id must be representable as a 32-bit integer.
+        
+            Returns: dictionary { 'solid_id' : solid_id, 
+                                  'channel_index' : channel_index,
+                                  'channel_id' : channel_id }
+        """
+
+        solid_id = self.add_solid(solid=pmt, rotation=rotation, 
+                                  displacement=displacement)
+
+        channel_index = len(self.channel_index_to_solid_id)
+        if channel_id is None:
+            channel_id = channel_index
+
+        # add_solid resized this array already
+        self.solid_id_to_channel_index[solid_id] = channel_index
+
+        # resize channel_index arrays before filling
+        self.channel_index_to_solid_id.resize(channel_index+1)
+        self.channel_index_to_solid_id[channel_index] = solid_id
+        self.channel_index_to_channel_id.resize(channel_index+1)
+        self.channel_index_to_channel_id[channel_index] = channel_id
+
+        # dictionary does not need resizing
+        self.channel_id_to_channel_index[channel_id] = channel_index
+        
+        return { 'solid_id' : solid_id, 
+                 'channel_index' : channel_index,
+                 'channel_id' : channel_id }
+                           
+    def _pdf_to_cdf(self, bin_edges, bin_contents):
+        '''Returns tuple of arrays `(cdf_x, cdf_y)` corresponding to
+        the cumulative distribution function for a binned PDF with the
+        given bin_edges and bin_contents.
+
+           `bin_edges`: array
+               bin edges of PDF.  length is len(bin_contents) + 1
+           `bin_contents`: array
+               The contents of each bin.  The array should NOT be
+               normalized for bin width if variable binning is present.
+        '''
+        cdf_x = np.copy(bin_edges)
+        cdf_y = np.array([0.0]+bin_contents.cumsum())
+        cdf_y /= cdf_y[-1] # normalize
+        return (cdf_x, cdf_y)
+
+    def set_time_dist_gaussian(self, rms, lo, hi, nsamples=50):
+        pdf_x = np.linspace(lo, hi, nsamples+1, endpoint=True)
+        pdf_y = np.exp(-0.5 * (pdf_x[1:]/rms)**2)
+        self.time_cdf = self._pdf_to_cdf(pdf_x, pdf_y)
+
+    def set_charge_dist_gaussian(self, mean, rms, lo, hi, nsamples=50):
+        pdf_x = np.linspace(lo, hi, nsamples+1, endpoint=True)
+        pdf_y = np.exp(-0.5 * ((pdf_x[1:] - mean)/rms)**2)
+        self.charge_cdf = self._pdf_to_cdf(pdf_x, pdf_y)
+
+    def num_channels(self):
+        return len(self.channel_index_to_channel_id)
diff --git a/chroma/gpu/__init__.py b/chroma/gpu/__init__.py
index b5a1ff3..74dc9f2 100644
--- a/chroma/gpu/__init__.py
+++ b/chroma/gpu/__init__.py
@@ -4,3 +4,4 @@ from chroma.gpu.render import *
 from chroma.gpu.photon import *
 from chroma.gpu.daq import *
 from chroma.gpu.pdf import *
+from chroma.gpu.detector import *
diff --git a/chroma/gpu/daq.py b/chroma/gpu/daq.py
index 6ec1a79..571feba 100644
--- a/chroma/gpu/daq.py
+++ b/chroma/gpu/daq.py
@@ -13,27 +13,30 @@ class GPUChannels(object):
 
     def get(self):
         t = self.t.get()
-        q = self.q.get().astype(np.float32)
+        q = self.q.get()
 
         # For now, assume all channels with small
         # enough hit time were hit.
         return event.Channels(t<1e8, t, q, self.flags.get())
 
 class GPUDaq(object):
-    def __init__(self, gpu_geometry, max_pmt_id, pmt_rms=1.2):
-        self.earliest_time_gpu = ga.empty(max_pmt_id+1, dtype=np.float32)
-        self.earliest_time_int_gpu = ga.empty(max_pmt_id+1, dtype=np.uint32)
+    def __init__(self, gpu_detector):
+        self.earliest_time_gpu = ga.empty(gpu_detector.nchannels, dtype=np.float32)
+        self.earliest_time_int_gpu = ga.empty(gpu_detector.nchannels, dtype=np.uint32)
         self.channel_history_gpu = ga.zeros_like(self.earliest_time_int_gpu)
-        self.channel_q_gpu = ga.zeros_like(self.earliest_time_int_gpu)
-        self.daq_pmt_rms = pmt_rms
-        self.solid_id_map_gpu = gpu_geometry.solid_id_map
+        self.channel_q_int_gpu = ga.zeros_like(self.earliest_time_int_gpu)
+        self.channel_q_gpu = ga.zeros(len(self.earliest_time_int_gpu), dtype=np.float32)
+        self.detector_gpu = gpu_detector.detector_gpu
+        self.solid_id_map_gpu = gpu_detector.solid_id_map
+        self.solid_id_to_channel_index_gpu = gpu_detector.solid_id_to_channel_index_gpu
 
         self.module = get_cu_module('daq.cu', options=cuda_options, 
-                                    include_source_directory=False)
+                                    include_source_directory=True)
         self.gpu_funcs = GPUFuncs(self.module)
 
     def acquire(self, gpuphotons, rng_states, nthreads_per_block=64, max_blocks=1024):
         self.gpu_funcs.reset_earliest_time_int(np.float32(1e9), np.int32(len(self.earliest_time_int_gpu)), self.earliest_time_int_gpu, block=(nthreads_per_block,1,1), grid=(len(self.earliest_time_int_gpu)//nthreads_per_block+1,1))
+        self.channel_q_int_gpu.fill(0)
         self.channel_q_gpu.fill(0)
         self.channel_history_gpu.fill(0)
 
@@ -41,9 +44,18 @@ class GPUDaq(object):
 
         for first_photon, photons_this_round, blocks in \
                 chunk_iterator(n, nthreads_per_block, max_blocks):
-            self.gpu_funcs.run_daq(rng_states, np.uint32(0x1 << 2), np.float32(self.daq_pmt_rms), np.int32(first_photon), np.int32(photons_this_round), gpuphotons.t, gpuphotons.flags, gpuphotons.last_hit_triangles, self.solid_id_map_gpu, np.int32(len(self.earliest_time_int_gpu)), self.earliest_time_int_gpu, self.channel_q_gpu, self.channel_history_gpu, block=(nthreads_per_block,1,1), grid=(blocks,1))
+            self.gpu_funcs.run_daq(rng_states, np.uint32(0x1 << 2), 
+                                   np.int32(first_photon), np.int32(photons_this_round), gpuphotons.t, 
+                                   gpuphotons.flags, gpuphotons.last_hit_triangles, 
+                                   self.solid_id_map_gpu,
+                                   self.detector_gpu,
+                                   self.earliest_time_int_gpu, 
+                                   self.channel_q_int_gpu, self.channel_history_gpu, 
+                                   block=(nthreads_per_block,1,1), grid=(blocks,1))
 
         self.gpu_funcs.convert_sortable_int_to_float(np.int32(len(self.earliest_time_int_gpu)), self.earliest_time_int_gpu, self.earliest_time_gpu, block=(nthreads_per_block,1,1), grid=(len(self.earliest_time_int_gpu)//nthreads_per_block+1,1))
 
+        self.gpu_funcs.convert_charge_int_to_float(self.detector_gpu, self.channel_q_int_gpu, self.channel_q_gpu, block=(nthreads_per_block,1,1), grid=(len(self.channel_q_int_gpu)//nthreads_per_block+1,1))
+
         return GPUChannels(self.earliest_time_gpu, self.channel_q_gpu, self.channel_history_gpu)
 
diff --git a/chroma/gpu/detector.py b/chroma/gpu/detector.py
new file mode 100644
index 0000000..ac0092e
--- /dev/null
+++ b/chroma/gpu/detector.py
@@ -0,0 +1,40 @@
+import numpy as np
+import pycuda.driver as cuda
+from pycuda import gpuarray as ga
+from pycuda import characterize
+
+from chroma.geometry import standard_wavelengths
+from chroma.gpu.tools import get_cu_module, get_cu_source, cuda_options, \
+    chunk_iterator, format_array, format_size, to_uint3, to_float3, \
+    make_gpu_struct
+from chroma.log import logger
+
+from chroma.gpu.geometry import GPUGeometry
+
+class GPUDetector(GPUGeometry):
+    def __init__(self, detector, wavelengths=None, print_usage=False):
+        GPUGeometry.__init__(self, detector, wavelengths=wavelengths, print_usage=False)
+        self.solid_id_to_channel_index_gpu = \
+            ga.to_gpu(detector.solid_id_to_channel_index.astype(np.int32))
+        self.nchannels = detector.num_channels()
+
+
+        self.time_cdf_x_gpu = ga.to_gpu(detector.time_cdf[0].astype(np.float32))
+        self.time_cdf_y_gpu = ga.to_gpu(detector.time_cdf[1].astype(np.float32))
+
+        self.charge_cdf_x_gpu = ga.to_gpu(detector.charge_cdf[0].astype(np.float32))
+        self.charge_cdf_y_gpu = ga.to_gpu(detector.charge_cdf[1].astype(np.float32))
+
+        detector_source = get_cu_source('detector.h')
+        detector_struct_size = characterize.sizeof('Detector', detector_source)
+        self.detector_gpu = make_gpu_struct(detector_struct_size,
+                                            [self.solid_id_to_channel_index_gpu,
+                                             self.time_cdf_x_gpu,
+                                             self.time_cdf_y_gpu,
+                                             self.charge_cdf_x_gpu,
+                                             self.charge_cdf_y_gpu,
+                                             np.int32(self.nchannels),
+                                             np.int32(len(detector.time_cdf[0])),
+                                             np.int32(len(detector.charge_cdf[0])),
+                                             np.float32(detector.charge_cdf[0][-1] / 2**16)])
+                                             
diff --git a/chroma/gpu/pdf.py b/chroma/gpu/pdf.py
index ba0d2af..5fb4b1a 100644
--- a/chroma/gpu/pdf.py
+++ b/chroma/gpu/pdf.py
@@ -76,7 +76,7 @@ class GPUKernelPDF(object):
             d = 2
         dimensionality_factor = ((4.0/(d+2)) / (mom0/scale_factor))**(-1.0/(d+4))
                                  
-        time_bandwidths = dimensionality_factor * trms 
+        time_bandwidths = dimensionality_factor * trms
         inv_time_bandwidths = np.zeros_like(time_bandwidths)
         inv_time_bandwidths[time_bandwidths  > 0] = time_bandwidths[time_bandwidths > 0] ** -1
         #inv_time_bandwidths /= 4.0
@@ -164,18 +164,18 @@ class GPUPDF(object):
                                     include_source_directory=False)
         self.gpu_funcs = GPUFuncs(self.module)
 
-    def setup_pdf(self, max_pmt_id, tbins, trange, qbins, qrange):
+    def setup_pdf(self, nchannels, tbins, trange, qbins, qrange):
         """Setup GPU arrays to hold PDF information.
 
-           max_pmt_id: int, largest PMT id #
+           nchannels: int, number of channels
            tbins: number of time bins
            trange: tuple of (min, max) time in PDF
            qbins: number of charge bins
            qrange: tuple of (min, max) charge in PDF
         """
         self.events_in_histogram = 0
-        self.hitcount_gpu = ga.zeros(max_pmt_id+1, dtype=np.uint32)
-        self.pdf_gpu = ga.zeros(shape=(max_pmt_id+1, tbins, qbins), 
+        self.hitcount_gpu = ga.zeros(nchannels, dtype=np.uint32)
+        self.pdf_gpu = ga.zeros(shape=(nchannels, tbins, qbins), 
                                       dtype=np.uint32)
         self.tbins = tbins
         self.trange = trange
diff --git a/chroma/sim.py b/chroma/sim.py
index d15c02d..f00f536 100644
--- a/chroma/sim.py
+++ b/chroma/sim.py
@@ -41,10 +41,13 @@ class Simulation(object):
             # need to build geometry
             detector.build()
 
-        self.gpu_geometry = gpu.GPUGeometry(detector)
-        self.gpu_daq = gpu.GPUDaq(self.gpu_geometry, max(self.detector.pmtids))
-        self.gpu_pdf = gpu.GPUPDF()
-        self.gpu_pdf_kernel = gpu.GPUKernelPDF()
+        if hasattr(detector, 'num_channels'):
+            self.gpu_geometry = gpu.GPUDetector(detector)
+            self.gpu_daq = gpu.GPUDaq(self.gpu_geometry)
+            self.gpu_pdf = gpu.GPUPDF()
+            self.gpu_pdf_kernel = gpu.GPUKernelPDF()
+        else:
+            self.gpu_geometry = gpu.GPUGeometry(detector)
 
         self.rng_states = gpu.get_rng_states(self.nthreads_per_block*self.max_blocks, seed=self.seed)
 
@@ -81,7 +84,8 @@ class Simulation(object):
             if keep_photons_end:
                 ev.photons_end = gpu_photons.get()
 
-            if run_daq:
+            # Skip running DAQ if we don't have one
+            if hasattr(self, 'gpu_daq') and run_daq:
                 gpu_channels = self.gpu_daq.acquire(gpu_photons, self.rng_states, nthreads_per_block=self.nthreads_per_block, max_blocks=self.max_blocks)
                 ev.channels = gpu_channels.get()
 
@@ -98,7 +102,7 @@ class Simulation(object):
         pdf_config = (tbins, trange, qbins, qrange)
         if pdf_config != self.pdf_config:
             self.pdf_config = pdf_config
-            self.gpu_pdf.setup_pdf(max(self.detector.pmtids), tbins, trange,
+            self.gpu_pdf.setup_pdf(self.detector.num_channels(), tbins, trange,
                                    qbins, qrange)
         else:
             self.gpu_pdf.clear_pdf()
diff --git a/test/test_detector.py b/test/test_detector.py
new file mode 100644
index 0000000..9660e59
--- /dev/null
+++ b/test/test_detector.py
@@ -0,0 +1,76 @@
+import unittest
+import numpy as np
+
+from chroma.geometry import Solid, Geometry, vacuum
+from chroma.detector import Detector
+from chroma.make import box
+from chroma.sim import Simulation
+from chroma.event import Photons
+
+from chroma.demo.optics import r7081hqe_photocathode
+
+class TestDetector(unittest.TestCase):
+    def setUp(self):
+        # Setup geometry
+        cube = Detector(vacuum)
+        cube.add_pmt(Solid(box(10.0,10,10), vacuum, vacuum, surface=r7081hqe_photocathode))
+        cube.set_time_dist_gaussian(1.2, -6.0, 6.0)
+        cube.set_charge_dist_gaussian(1.0, 0.1, 0.5, 1.5)
+
+        cube.build(use_cache=False)
+
+        self.cube = cube
+        self.sim = Simulation(cube, geant4_processes=0)
+        
+    def testTime(self):
+        '''Test PMT time distribution'''
+
+        # Run only one photon at a time
+        nphotons = 1
+        pos = np.tile([0,0,0], (nphotons,1)).astype(np.float32)
+        dir = np.tile([0,0,1], (nphotons,1)).astype(np.float32)
+        pol = np.zeros_like(pos)
+        phi = np.random.uniform(0, 2*np.pi, nphotons).astype(np.float32)
+        pol[:,0] = np.cos(phi)
+        pol[:,1] = np.sin(phi)
+        t = np.zeros(nphotons, dtype=np.float32) + 100.0 # Avoid negative photon times
+        wavelengths = np.empty(nphotons, np.float32)
+        wavelengths.fill(400.0)
+
+        photons = Photons(pos=pos, dir=dir, pol=pol, t=t,
+                          wavelengths=wavelengths)
+
+        hit_times = []
+        for ev in self.sim.simulate(photons for i in xrange(10000)):
+            if ev.channels.hit[0]:
+                hit_times.append(ev.channels.t[0])
+        hit_times = np.array(hit_times)
+
+        self.assertAlmostEqual(hit_times.std(),  1.2, delta=1e-1)
+
+
+    def testCharge(self):
+        '''Test PMT charge distribution'''
+
+        # Run only one photon at a time
+        nphotons = 1
+        pos = np.tile([0,0,0], (nphotons,1)).astype(np.float32)
+        dir = np.tile([0,0,1], (nphotons,1)).astype(np.float32)
+        pol = np.zeros_like(pos)
+        phi = np.random.uniform(0, 2*np.pi, nphotons).astype(np.float32)
+        pol[:,0] = np.cos(phi)
+        pol[:,1] = np.sin(phi)
+        t = np.zeros(nphotons, dtype=np.float32)
+        wavelengths = np.empty(nphotons, np.float32)
+        wavelengths.fill(400.0)
+
+        photons = Photons(pos=pos, dir=dir, pol=pol, t=t,
+                          wavelengths=wavelengths)
+
+        hit_charges = []
+        for ev in self.sim.simulate(photons for i in xrange(1000)):
+            if ev.channels.hit[0]:
+                hit_charges.append(ev.channels.q[0])
+        hit_charges = np.array(hit_charges)
+        self.assertAlmostEqual(hit_charges.mean(),  1.0, delta=1e-1)
+        self.assertAlmostEqual(hit_charges.std(), 0.1, delta=1e-1)
diff --git a/test/test_pdf.py b/test/test_pdf.py
index 2eafd67..df13a2a 100644
--- a/test/test_pdf.py
+++ b/test/test_pdf.py
@@ -22,15 +22,15 @@ class TestPDF(unittest.TestCase):
 
         context = gpu.create_cuda_context()
 
-        gpu_geometry = gpu.GPUGeometry(self.detector)
+        gpu_geometry = gpu.GPUDetector(self.detector)
 
         nthreads_per_block, max_blocks = 64, 1024
 
         rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks)
 
-        gpu_daq = gpu.GPUDaq(gpu_geometry, max(self.detector.pmtids))
+        gpu_daq = gpu.GPUDaq(gpu_geometry)
         gpu_pdf = gpu.GPUPDF()
-        gpu_pdf.setup_pdf(max(self.detector.pmtids), 100, (-0.5, 999.5), 10, (-0.5, 9.5))
+        gpu_pdf.setup_pdf(self.detector.num_channels(), 100, (-0.5, 999.5), 10, (-0.5, 9.5))
         
         gpu_pdf.clear_pdf()
 
diff --git a/test/test_propagation.py b/test/test_propagation.py
index d3dc2ed..34e64f5 100644
--- a/test/test_propagation.py
+++ b/test/test_propagation.py
@@ -18,7 +18,6 @@ class TestPropagation(unittest.TestCase):
         # Setup geometry
         cube = Geometry(vacuum)
         cube.add_solid(Solid(box(100,100,100), vacuum, vacuum))
-        cube.pmtids = [0]
         cube.build(use_cache=False)
         sim = Simulation(cube, geant4_processes=0)
 
diff --git a/test/test_rayleigh.py b/test/test_rayleigh.py
index 487184f..7657fc1 100644
--- a/test/test_rayleigh.py
+++ b/test/test_rayleigh.py
@@ -15,7 +15,6 @@ class TestRayleigh(unittest.TestCase):
     def setUp(self):
         self.cube = Geometry(water)
         self.cube.add_solid(Solid(box(100,100,100), water, water))
-        self.cube.pmtids = [0]
         self.cube.build(use_cache=False)
         self.sim = Simulation(self.cube, geant4_processes=0)