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.

4 files changed, 67 insertions, 14 deletions

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