11 files changed, 75 insertions, 9 deletions

diff --git a/color/__init__.py b/color/__init__.py
new file mode 100644
index 0000000..c5fda6a
--- /dev/null
+++ b/color/__init__.py
@@ -0,0 +1 @@
+from chromaticity import map_wavelength
diff --git a/color/chromaticity.py b/color/chromaticity.py
index 6229b51..fbd812a 100644
--- a/color/chromaticity.py
+++ b/color/chromaticity.py
@@ -1,6 +1,9 @@
 import numpy as np
+import os
 
-f = open('sbrgb10w.csv')
+dir = os.path.split(os.path.realpath(__file__))[0]
+
+f = open(dir + '/sbrgb10w.csv')
 
 color_map = []
 for line in f:
diff --git a/models/hamamatsu_12inch_inner.stl b/models/hamamatsu_12inch_inner.stl
index a324f2d..fc805ad 100644
--- a/models/hamamatsu_12inch_inner.stl
+++ b/models/hamamatsu_12inch_inner.stl
diff --git a/models/hamamatsu_12inch_outer.stl b/models/hamamatsu_12inch_outer.stl
index 77dfc14..23b2c55 100644
--- a/models/hamamatsu_12inch_outer.stl
+++ b/models/hamamatsu_12inch_outer.stl
diff --git a/models/snopmt_inner.stl b/models/snopmt_inner.stl
index 35d6ab8..24929ca 100644
--- a/models/snopmt_inner.stl
+++ b/models/snopmt_inner.stl
diff --git a/models/snopmt_outer.stl b/models/snopmt_outer.stl
index de8a4d8..d498559 100644
--- a/models/snopmt_outer.stl
+++ b/models/snopmt_outer.stl
diff --git a/src/intersect.h b/src/intersect.h
index bfde782..0b713c8 100644
--- a/src/intersect.h
+++ b/src/intersect.h
@@ -60,7 +60,7 @@ __device__ bool intersect_triangle(const float3 &origin, const float3 &direction
    `direction` must be normalized. */
 __device__ int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2, const int base_color=0xFFFFFFFF)
 {
-	float scale = dot(normalize(cross(v1-v0,v2-v0)),-direction);
+	float scale = dot(normalize(cross(v1-v0,v2-v1)),-direction);
 
 	unsigned int r = 0xFF & (base_color >> 16);
 	unsigned int g = 0xFF & (base_color >> 8);
diff --git a/src/kernel.cu b/src/kernel.cu
index b7dadfa..6c3ef1b 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -267,7 +267,7 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f
 		int material_out_index = tex1Dfetch(material2_lookup, last_hit_triangle);
 		int surface_index = tex1Dfetch(surface_lookup, last_hit_triangle);
 
-		float3 surface_normal = cross(v1-v0,v2-v0);
+		float3 surface_normal = cross(v1-v0,v2-v1);
 		surface_normal /= norm(surface_normal);
 		
 		Material material1, material2;
diff --git a/stl.py b/stl.py
index 9bd3a38..f88c3fa 100644
--- a/stl.py
+++ b/stl.py
@@ -60,7 +60,7 @@ def mesh_from_binary_stl(filename):
     ntriangles = struct.unpack('<I', f.read(4))[0]
 
     for i in range(ntriangles):
-        f.read(12)
+        normal = tuple(struct.unpack('<fff', f.read(12)))
 
         triangle = [None]*3
         for j in range(3):
diff --git a/threadtest.py b/threadtest.py
index 81b6d62..91a1857 100644
--- a/threadtest.py
+++ b/threadtest.py
@@ -11,7 +11,7 @@ from histogram import Histogram
 jobs = Queue()
 output = Queue()
 
-def create_job(position=(0,0,0), nphotons=10000, max_steps=10):
+def create_job(position=(0,0,0), nphotons=1000, max_steps=10):
     positions = np.tile(position, nphotons).reshape((nphotons, 3))
     directions = uniform_sphere(nphotons)
     polarizations = uniform_sphere(nphotons)
@@ -23,7 +23,7 @@ def create_job(position=(0,0,0), nphotons=10000, max_steps=10):
     return Job(positions, directions, wavelengths, polarizations, times,
                states, last_hit_triangles, max_steps)
 
-def generate_event(position=(0,0,0), nphotons=10000):
+def generate_event(position=(0,0,0), nphotons=1000):
     jobs.put(create_job(position, nphotons))
     jobs.join()
 
@@ -34,6 +34,9 @@ def generate_event(position=(0,0,0), nphotons=10000):
     solids[last_hit_triangles == -1] = -1
     surface_absorbed = job.states == 2
 
+    for i in np.unique(job.states):
+        print 'state %2i %i' % (i, len(job.states[job.states == i]))
+
     event_times = []
     for i in lbne.pmtids:
         photons = np.logical_and(solids == i, surface_absorbed)
@@ -43,12 +46,11 @@ def generate_event(position=(0,0,0), nphotons=10000):
         if len(hit_times) > 0:
             event_times.append((i, np.min(hit_times)))
 
-    print len(event_times)
-    print event_times
+    print '%i hit pmts' % len(event_times)
 
     return event_times
 
-def likelihood(event_times, position=(0,0,0), nphotons=10000, neval=100):
+def likelihood(event_times, position=(0,0,0), nphotons=1000, neval=100):
     for i in range(neval):
         jobs.put(create_job(position, nphotons))
     jobs.join()
diff --git a/track.py b/track.py
new file mode 100644
index 0000000..f73f2ed
--- /dev/null
+++ b/track.py
@@ -0,0 +1,60 @@
+import numpy as np
+import pycuda.driver as cuda
+from gputhread import GPUThread
+from detectors import LBNE
+from Queue import Queue
+from threadtest import create_job
+import matplotlib.pyplot as plt
+from itertoolset import roundrobin
+from color import map_wavelength
+from solids import r7081
+from geometry import Geometry
+
+nphotons = 100000
+
+jobs = Queue()
+output = Queue()
+
+#geometry = LBNE()
+geometry = Geometry()
+geometry.add_solid(r7081, displacement=(0,-1,0))
+geometry.build(bits=8)
+
+cuda.init()
+
+try:
+    gputhread = GPUThread(0, geometry, jobs, output, 64)
+    gputhread.start()
+
+    job = create_job((0,0,0), nphotons)
+
+    x = np.empty((nphotons, 10, 3))
+    for i in range(10):
+        print '%i' % i
+        x[:,i,0] = job.positions['x']
+        x[:,i,1] = job.positions['y']
+        x[:,i,2] = job.positions['z']
+
+        jobs.put(job)
+        jobs.join()
+
+        job = output.get()
+
+        for j in np.unique(job.states):
+            print 'state %2i, %i' % (j, len(job.states[job.states == j]))
+finally:
+    gputhread.stop()
+    gputhread.join()
+
+mask = job.states != 0
+
+rgb = (map_wavelength(job.wavelengths[mask])*255).astype(np.uint32)
+
+def format_hex_string(s):
+    return '#' + s.rstrip('L')[2:].zfill(6)
+
+colors = map(format_hex_string, map(hex, rgb[:,0] << 16 | rgb[:,1] << 8 | rgb[:,2]))
+
+plt.figure()
+plt.plot(*roundrobin(x[mask,:,0], x[mask,:,1], colors))
+plt.show()