Epic port of Chroma from units of meters/seconds/MeV to

millimeters/nanoseconds/MeV in order to match GEANT4, and also avoid
huge discrepancies in magnitude caused by values like 10e-9 sec.

Along the way, cleaned up a few things:
 * Switch the PI and SPEED_OF_LIGHT constants from double to single
   precision.  This avoid some unnecessary double precision calculations
   in the GPU code.

 * Fixed a silly problem in the definition of the spherical spiral.  Now
   the demo detector looks totally awesome.  Also wrapped it in a
   black surface.  Demo detector now has 10055 PMTs.

 * Updated the test_ray_intersection data file to reflect the new units.

 * Fix a missing import in chroma.gpu.tools

16 files changed, 39 insertions, 39 deletions

diff --git a/chroma/camera.py b/chroma/camera.py
index 32c7e88..628307a 100644
--- a/chroma/camera.py
+++ b/chroma/camera.py
@@ -143,7 +143,7 @@ class Camera(multiprocessing.Process):
         self.axis1 = np.array([0,0,1], float)
         self.axis2 = np.array([1,0,0], float)
 
-        self.film_width = 0.035
+        self.film_width = 35.0 # mm
 
         pos, dir = from_film(self.point, axis1=self.axis1, axis2=self.axis2,
                              size=self.size, width=self.film_width)
diff --git a/chroma/cuda/physical_constants.h b/chroma/cuda/physical_constants.h
index 2ff87cd..0082503 100644
--- a/chroma/cuda/physical_constants.h
+++ b/chroma/cuda/physical_constants.h
@@ -1,7 +1,9 @@
 #ifndef __PHYSICAL_CONSTANTS_H__
 #define __PHYSICAL_CONSTANTS_H__
 
-#define SPEED_OF_LIGHT 2.99792458e8
-#define PI 3.141592653589793
+// mm/ns
+#define SPEED_OF_LIGHT 299.792458f
+
+#define PI 3.141592653589793f
 
 #endif
diff --git a/chroma/demo/__init__.py b/chroma/demo/__init__.py
index edd344a..7df7df7 100644
--- a/chroma/demo/__init__.py
+++ b/chroma/demo/__init__.py
@@ -7,10 +7,10 @@ import numpy.linalg
 from chroma.make import sphere
 from chroma.stl import mesh_from_stl
 from chroma.geometry import *
-from chroma.transform import rotate, make_rotation_matrix
+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
+from chroma.demo.optics import water, black_surface
 
 from chroma.demo.checkerboard import build_checkerboard_scene as checkerboard_scene
 
@@ -24,22 +24,23 @@ def spherical_spiral(radius, spacing):
 
     while t < np.pi:
         yield np.array([sin(t) * sin(a*t), sin(t) * cos(a*t), cos(t)])*radius
-        dt = spacing / sqrt(1 + a**2 * sin(t) ** 2)
+        dt = dl / sqrt(1 + a**2 * sin(t) ** 2)
         t += dt
 
-def detector(pmt_radius=30.0, sphere_radius=30.5, spiral_step=0.14):
+def detector(pmt_radius=14000.0, sphere_radius=14500.0, spiral_step=350.0):
     pmt = build_8inch_pmt_with_lc()
     geo = Geometry(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 = -position/numpy.linalg.norm(position)
+        direction = -normalize(position)
 
         # Orient PMT that starts facing Y axis
         y_axis = np.array((0.0,1.0,0.0))
@@ -55,5 +56,5 @@ def detector(pmt_radius=30.0, sphere_radius=30.5, spiral_step=0.14):
     return geo
 
 def tiny():
-    return detector(1.0, 1.5, 0.4)
+    return detector(2000.0, 2500.0, 700.0)
 
diff --git a/chroma/demo/checkerboard.py b/chroma/demo/checkerboard.py
index dd2e3e1..801f9da 100644
--- a/chroma/demo/checkerboard.py
+++ b/chroma/demo/checkerboard.py
@@ -6,8 +6,8 @@ from chroma.make import sphere
 from chroma.demo.optics import *
 
 def build_checkerboard_scene(checkers_per_side=10, squares_per_checker=50):
-    x = np.linspace(-5.0, 5.0, checkers_per_side*squares_per_checker+1)
-    y = np.linspace(-5.0, 5.0, checkers_per_side*squares_per_checker+1)
+    x = np.linspace(-5000.0, 5000.0, checkers_per_side*squares_per_checker+1)
+    y = np.linspace(-5000.0, 5000.0, checkers_per_side*squares_per_checker+1)
 
     vertices = np.array(tuple(product(x,y,[0])))
 
@@ -29,14 +29,15 @@ def build_checkerboard_scene(checkers_per_side=10, squares_per_checker=50):
 
     checkerboard = Solid(checkerboard_mesh, vacuum, vacuum, surface=checkerboard_surface, color=checkerboard_color)
 
-    sphere1 = Solid(sphere(nsteps=512), water, vacuum)
-    sphere2 = Solid(sphere(nsteps=512), vacuum, vacuum, surface=shiny_surface)
-    sphere3 = Solid(sphere(nsteps=512), vacuum, vacuum, surface=lambertian_surface)
+    sphere1 = Solid(sphere(1000.0, nsteps=512), water, vacuum)
+    sphere2 = Solid(sphere(1000.0, nsteps=512), vacuum, vacuum, 
+                    surface=shiny_surface)
+    sphere3 = Solid(sphere(1000.0, nsteps=512), vacuum, vacuum, surface=lambertian_surface)
 
     checkerboard_scene = Geometry()
-    checkerboard_scene.add_solid(checkerboard, displacement=(0,0,-1.5))
-    checkerboard_scene.add_solid(sphere1, displacement=(2.0,-2.0,0))
-    checkerboard_scene.add_solid(sphere2, displacement=(-2.0,-2.0,0))
-    checkerboard_scene.add_solid(sphere3, displacement=(0.0,2.0,0))
+    checkerboard_scene.add_solid(checkerboard, displacement=(0,0,-1500.0))
+    checkerboard_scene.add_solid(sphere1, displacement=(2000.0,-2000.0,0))
+    checkerboard_scene.add_solid(sphere2, displacement=(-2000.0,-2000.0,0))
+    checkerboard_scene.add_solid(sphere3, displacement=(0.0,2000.0,0))
 
     return checkerboard_scene
diff --git a/chroma/demo/optics.py b/chroma/demo/optics.py
index 2fd34d0..db14032 100644
--- a/chroma/demo/optics.py
+++ b/chroma/demo/optics.py
@@ -51,7 +51,7 @@ r7081hqe_photocathode.set('reflect_diffuse', 1.0 - r7081hqe_photocathode.detect[
 glass = Material('glass')
 glass.set('refractive_index', 1.49)
 glass.absorption_length = \
-    np.array([(200, 0.1e-6), (300, 0.1e-6), (330, 1.0), (500, 2.0), (600, 1.0), (770, 0.5), (800, 0.1e-6)])
+    np.array([(200, 0.1e-6), (300, 0.1e-6), (330, 1000.0), (500, 2000.0), (600, 1000.0), (770, 500.0), (800, 0.1e-6)])
 glass.set('scattering_length', 1e6)
 
 # From WCSim
@@ -99,7 +99,7 @@ water.set('absorption_length',
                  11751,   8795.34, 8741.23, 7102.37, 6060.68,
                  4498.56, 3039.56, 2232.2,  1938,    1811.58,
                  1610.32, 1338.7,  1095.3,  977.525, 965.258,
-                 1082.86, 876.434, 633.723, 389.87,  142.011])[::-1] / 100.0 # reversed, cm->m
+                 1082.86, 876.434, 633.723, 389.87,  142.011])[::-1] * 10.0 # reversed, cm->mm
                 )
       
 water.set('scattering_length',
@@ -123,5 +123,5 @@ water.set('scattering_length',
                           3960.210, 3473.413, 3032.937,
                           2635.746, 2278.907, 1959.588,
                           1675.064, 1422.710, 1200.004,
-                          1004.528, 833.9666, 686.1063])[::-1] / 100.0 * 0.625 # reversed, cm -> m, * magic tuning constant
+                          1004.528, 833.9666, 686.1063])[::-1] * 10.0 * 0.625 # reversed, cm -> mm, * magic tuning constant
           )
diff --git a/chroma/demo/pmt.py b/chroma/demo/pmt.py
index 17af87a..9b0953f 100644
--- a/chroma/demo/pmt.py
+++ b/chroma/demo/pmt.py
@@ -5,7 +5,7 @@ from chroma.pmt import build_pmt, build_light_collector_from_file
 from chroma.demo.optics import water, glass, vacuum, shiny_surface, r7081hqe_photocathode
 
 def build_8inch_pmt(outer_material=water, nsteps=24):
-    return build_pmt(dirname(__file__) + '/sno_pmt.txt', 0.003,
+    return build_pmt(dirname(__file__) + '/sno_pmt.txt', 3.0, # 3 mm
                      outer_material=outer_material,
                      glass=glass, vacuum=vacuum,
                      photocathode_surface=r7081hqe_photocathode,
diff --git a/chroma/generator/g4gen.py b/chroma/generator/g4gen.py
index 0ed5433..0132381 100644
--- a/chroma/generator/g4gen.py
+++ b/chroma/generator/g4gen.py
@@ -116,9 +116,9 @@ class G4Generator(object):
             mass = G4ParticleTable.GetParticleTable().FindParticle(vertex.particle_name).GetPDGMass()
             total_energy = vertex.ke*MeV + mass
             self.particle_gun.SetParticleEnergy(total_energy)
-            self.particle_gun.SetParticlePosition(G4ThreeVector(*vertex.pos)*m)
+            self.particle_gun.SetParticlePosition(G4ThreeVector(*vertex.pos)*mm)
             self.particle_gun.SetParticleMomentumDirection(G4ThreeVector(*vertex.dir).unit())
-            self.particle_gun.SetParticleTime(vertex.t0*s)
+            self.particle_gun.SetParticleTime(vertex.t0*ns)
             
             if vertex.pol is not None:
                 self.particle_gun.SetParticlePolarization(G4ThreeVector(*vertex.pol).unit())
diff --git a/chroma/gpu/daq.py b/chroma/gpu/daq.py
index 82958dc..6ec1a79 100644
--- a/chroma/gpu/daq.py
+++ b/chroma/gpu/daq.py
@@ -20,7 +20,7 @@ class GPUChannels(object):
         return event.Channels(t<1e8, t, q, self.flags.get())
 
 class GPUDaq(object):
-    def __init__(self, gpu_geometry, max_pmt_id, pmt_rms=1.2e-9):
+    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)
         self.channel_history_gpu = ga.zeros_like(self.earliest_time_int_gpu)
diff --git a/chroma/gpu/tools.py b/chroma/gpu/tools.py
index a11a561..38d67a1 100644
--- a/chroma/gpu/tools.py
+++ b/chroma/gpu/tools.py
@@ -2,6 +2,7 @@ import numpy as np
 import pycuda.tools
 from pycuda import characterize
 import pycuda.driver as cuda
+import pycuda.compiler
 from pycuda import gpuarray as ga
 
 from chroma.cuda import srcdir
diff --git a/chroma/likelihood.py b/chroma/likelihood.py
index c69dc72..93a95c3 100644
--- a/chroma/likelihood.py
+++ b/chroma/likelihood.py
@@ -6,7 +6,7 @@ from chroma.tools import profile_if_possible
 
 class Likelihood(object):
     "Class to evaluate likelihoods for detector events."
-    def __init__(self, sim, event=None, tbins=100, trange=(-0.5e-9, 999.5e-9), 
+    def __init__(self, sim, event=None, tbins=100, trange=(-0.5, 999.5), 
                  qbins=10, qrange=(-0.5, 9.5), time_only=True):
         """
         Args:
@@ -56,7 +56,7 @@ class Likelihood(object):
         hitcount, pdf_prob, pdf_prob_uncert = \
             self.sim.eval_pdf(self.event.channels,
                               vertex_generator,
-                              0.1e-9, self.trange, 
+                              0.1, self.trange, 
                               1, self.qrange,
                               nreps=nreps,
                               ndaq=ndaq,
diff --git a/chroma/pmt.py b/chroma/pmt.py
index 6ebf120..6f9e0df 100644
--- a/chroma/pmt.py
+++ b/chroma/pmt.py
@@ -34,8 +34,6 @@ def build_pmt_shell(filename, outer_material, nsteps=16):
     # so that the mesh is closed
     profile[0,0] = 0.0
     profile[-1,0] = 0.0
-    # convert mm -> m
-    profile /= 1000.0
 
     return Solid(rotate_extrude(profile[:,0], profile[:,1], nsteps), glass, outer_material, color=0xeeffffff)
 
@@ -52,8 +50,6 @@ def build_pmt(filename, glass_thickness, outer_material, glass,
     # so that the mesh is closed
     profile[0,0] = 0.0
     profile[-1,0] = 0.0
-    # convert mm -> m
-    profile /= 1000.0
 
     offset_profile = offset(profile, -glass_thickness)
 
@@ -80,9 +76,6 @@ def build_light_collector_from_file(filename, outer_material,
                                     surface, nsteps=48):
     profile = read_csv(filename)
     
-    # Convert mm to m
-    profile /= 1000.0
-
     mesh = rotate_extrude(profile[:,0], profile[:,1], nsteps)
     solid = Solid(mesh, outer_material, outer_material, surface=surface)
     return solid
diff --git a/chroma/project.py b/chroma/project.py
index 51fba1e..05e4d2e 100644
--- a/chroma/project.py
+++ b/chroma/project.py
@@ -2,7 +2,7 @@ import numpy as np
 from chroma.transform import normalize
 
 def from_film(position=(0,0,0), axis1=(0,0,1), axis2=(1,0,0), size=(800,600),
-              width=0.035, focal_length=0.018):
+              width=35.0, focal_length=18.0):
     """Project rays from a piece of film located at whose focal point is
     located at `position`. `axis1` and `axis2` specify the vectors pointing
     along the length and height of the film respectively.
diff --git a/src/G4chroma.cc b/src/G4chroma.cc
index 340cc4c..c0817a1 100644
--- a/src/G4chroma.cc
+++ b/src/G4chroma.cc
@@ -107,11 +107,11 @@ void PhotonTrackingAction::PreUserTrackingAction(const G4Track *track)
 {
   G4ParticleDefinition *particle = track->GetDefinition();
   if (particle->GetParticleName() == "opticalphoton") {
-    pos.push_back(track->GetPosition()/m);
+    pos.push_back(track->GetPosition()/mm);
     dir.push_back(track->GetMomentumDirection());
     pol.push_back(track->GetPolarization());
     wavelength.push_back( (h_Planck * c_light / track->GetKineticEnergy()) / nanometer );
-    t0.push_back(track->GetGlobalTime() / s);
+    t0.push_back(track->GetGlobalTime() / ns);
     const_cast<G4Track *>(track)->SetTrackStatus(fStopAndKill);
   }
 }
diff --git a/test/data/ray_intersection.npy b/test/data/ray_intersection.npy
new file mode 100644
index 0000000..43f51c2
--- /dev/null
+++ b/test/data/ray_intersection.npy
diff --git a/test/data/ray_intersection.npz b/test/data/ray_intersection.npz
deleted file mode 100644
index fba0e40..0000000
--- a/test/data/ray_intersection.npz
+++ /dev/null
diff --git a/test/test_ray_intersection.py b/test/test_ray_intersection.py
index ce364e0..8a33ac3 100644
--- a/test/test_ray_intersection.py
+++ b/test/test_ray_intersection.py
@@ -9,7 +9,7 @@ class TestRayIntersection(unittest.TestCase):
         self.context = chroma.gpu.create_cuda_context()
         self.module = chroma.gpu.get_cu_module('mesh.h')
         self.gpu_funcs = chroma.gpu.GPUFuncs(self.module)
-        self.box = chroma.gpu.GPUGeometry(chroma.build(chroma.make.cube()))
+        self.box = chroma.gpu.GPUGeometry(chroma.build(chroma.make.cube(size=1000.0)))
 
         pos, dir = chroma.project.from_film()
         self.pos_gpu = ga.to_gpu(chroma.gpu.to_float3(pos))
@@ -17,10 +17,12 @@ class TestRayIntersection(unittest.TestCase):
 
         testdir = os.path.dirname(os.path.abspath(__file__))
         self.dx_standard = np.load(os.path.join(testdir,
-                                                'data/ray_intersection.npz'))
+                                                'data/ray_intersection.npy'))
+
     def test_intersection_distance(self):
         dx = ga.zeros(self.pos_gpu.size, dtype=np.float32)
         self.gpu_funcs.distance_to_mesh(np.int32(self.pos_gpu.size), self.pos_gpu, self.dir_gpu, self.box.gpudata, dx, block=(64,1,1), grid=(self.pos_gpu.size//64+1,1))
+
         self.assertTrue((dx.get() == self.dx_standard).all())
 
     def tearDown(self):