diff options
Diffstat (limited to 'threadtest.py')
| -rw-r--r-- | threadtest.py | 97 |
1 files changed, 68 insertions, 29 deletions
diff --git a/threadtest.py b/threadtest.py index 4d41580..08a1b42 100644 --- a/threadtest.py +++ b/threadtest.py @@ -1,7 +1,7 @@ from gputhread import * from Queue import Queue from detectors import LBNE -from photon import uniform_sphere +from sample import uniform_sphere import numpy as np from pycuda import gpuarray import pycuda.driver as cuda @@ -11,47 +11,86 @@ from histogram import Histogram jobs = Queue() output = Queue() -def generate_event(pos=(0,0,0), nphotons=1000): - origins_float3 = np.tile(gpuarray.vec.make_float3(*pos), (nphotons,1)) +def create_job(position=(0,0,0), nphotons=10000, max_steps=10): + positions = np.tile(gpuarray.vec.make_float3(*position), (nphotons,1)) directions = uniform_sphere(nphotons) - directions_float3 = np.empty(nphotons, dtype=gpuarray.vec.float3) - directions_float3['x'] = directions[:,0] - directions_float3['y'] = directions[:,1] - directions_float3['z'] = directions[:,2] + polarizations = uniform_sphere(nphotons) + wavelengths = np.random.uniform(200, 800, size=nphotons) + times = np.zeros(nphotons) + states = -np.ones(nphotons) + last_hit_triangles = -np.ones(nphotons) + + return Job(positions, directions, wavelengths, polarizations, times, + states, last_hit_triangles, max_steps) + +#@profile +def generate_event(position=(0,0,0), nphotons=10000): + jobs.put(create_job(position, nphotons)) + jobs.join() - jobs.put(Job(origins_float3, directions_float3)) + job = output.get() - jobs.join() + last_hit_triangles = job.last_hit_triangles + solids = lbne.solid_id[last_hit_triangles] + solids[last_hit_triangles == -1] = -1 + surface_absorbed = job.states == 2 - return output.get() + event_times = [] + for i in lbne.pmtids: + photons = np.logical_and(solids == i, surface_absorbed) -def likelihood(event_bincount, pos=(0,0,0), nphotons=1000, neval=1000): - for i in range(neval): - origins_float3 = np.tile(gpuarray.vec.make_float3(*pos), (nphotons,1)) - directions = uniform_sphere(nphotons) - directions_float3 = np.empty(nphotons, dtype=gpuarray.vec.float3) - directions_float3['x'] = directions[:,0] - directions_float3['y'] = directions[:,1] - directions_float3['z'] = directions[:,2] + hit_times = job.times[photons] - jobs.put(Job(origins_float3, directions_float3)) + if len(hit_times) > 0: + event_times.append((i, np.min(hit_times))) + print event_times + + return np.array(event_times) + +#@profile +def likelihood(event_times, position=(0,0,0), nphotons=10000, neval=100): + for i in range(neval): + jobs.put(create_job(position, nphotons)) jobs.join() - bincount = np.zeros((neval, len(lbne.solids))) + t = {} for i in range(neval): - bincount[i] = output.get() + job = output.get() + + last_hit_triangles = job.last_hit_triangles + solids = lbne.solid_id[job.last_hit_triangles] + solids[last_hit_triangles == -1] = -1 + surface_absorbed = job.states == 2 + + for j in lbne.pmtids: + pmt_photons = solids == j + photons = np.logical_and(pmt_photons, surface_absorbed) + + if j not in t: + t[j] = [] + + hit_times = job.times[photons] + + if len(hit_times) > 0: + t[j].append(np.min(hit_times)) log_likelihood = ufloat((0,0)) - for i in range(len(lbne.solids)): - h = Histogram(100, (-0.5, 99.5)) - h.fill(bincount[:,i]) - h.normalize() + for i, time in event_times: + h = Histogram(100, (-0.5e-9, 99.5e-9)) + h.fill(t[i]) + + if h.nentries > 0: + h.normalize() - probability = h.ueval(event_bincount[i]) + probability = h.ueval(time) if probability.nominal_value == 0.0: - probability = ufloat((0.5/h.nentries, 0.5/h.nentries)) + if h.nentries > 0: + probability = ufloat((0.5/h.nentries, 0.5/h.nentries)) + else: + probability = \ + ufloat((1.0/len(h.bincenters), 1.0/len(h.bincenters))) log_likelihood += umath.log(probability) @@ -79,11 +118,11 @@ if __name__ == '__main__': gputhreads[-1].start() try: - event_bincount = generate_event() + event_times = generate_event() for z in np.linspace(-1.0, 1.0, 100): t0 = time.time() - log_likelihood = likelihood(event_bincount, (z,0,0)) + log_likelihood = likelihood(event_times, (z,0,0)) elapsed = time.time() - t0 print 'z = %5.2f, likelihood = %s, elapsed %.2f sec' % (z, log_likelihood, elapsed) finally: |