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, 79 insertions, 5 deletions

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)