From 06ece999c3866f2d19acfd0f23b2f62d02b50577 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Tue, 9 Aug 2011 15:43:36 -0400 Subject: Store a photon history for each hit channel. If multiple photons hit the same channel, their history bits are OR'ed together. --- gpu.py | 10 ++++++++-- io/root.C | 13 ++++++++----- sim.py | 7 ++++--- src/daq.cu | 6 ++++-- 4 files changed, 24 insertions(+), 12 deletions(-) diff --git a/gpu.py b/gpu.py index a21bc38..7b08473 100644 --- a/gpu.py +++ b/gpu.py @@ -278,6 +278,7 @@ class GPU(object): self.earliest_time_gpu = gpuarray.GPUArray(shape=(max_pmt_id+1,), dtype=np.float32) self.earliest_time_int_gpu = gpuarray.GPUArray(shape=self.earliest_time_gpu.shape, dtype=np.uint32) + self.channel_history_gpu = gpuarray.zeros_like(self.earliest_time_int_gpu) self.daq_pmt_rms = pmt_rms def run_daq(self): @@ -286,6 +287,8 @@ class GPU(object): self.earliest_time_int_gpu, block=(self.nthread_per_block,1,1), grid=(len(self.earliest_time_int_gpu)//self.nthread_per_block+1,1)) + self.channel_history_gpu.fill(0) + #self.context.synchronize() for first_photon, photons_this_round, blocks in chunk_iterator(self.nphotons, self.nthread_per_block, self.max_blocks): self.daq_funcs.run_daq(self.rng_states_gpu, np.uint32(0x1 << 2), @@ -296,7 +299,8 @@ class GPU(object): self.last_hit_triangles_gpu, self.solid_id_map_gpu, np.int32(len(self.earliest_time_int_gpu)), - self.earliest_time_int_gpu, + self.earliest_time_int_gpu, + self.channel_history_gpu, block=(self.nthread_per_block,1,1), grid=(blocks,1)) #self.context.synchronize() self.daq_funcs.convert_sortable_int_to_float(np.int32(len(self.earliest_time_int_gpu)), @@ -308,7 +312,9 @@ class GPU(object): def get_hits(self): - return self.earliest_time_gpu.get() + return { 't': self.earliest_time_gpu.get(), + 'q': np.ones(shape=len(self.earliest_time_int_gpu), dtype=np.float32), + 'history': self.channel_history_gpu.get()} def __del__(self): self.context.pop() diff --git a/io/root.C b/io/root.C index 43f4e73..096811a 100644 --- a/io/root.C +++ b/io/root.C @@ -9,7 +9,7 @@ struct Photon { TVector3 dir; TVector3 pol; double wavelength; // nm - int history; + unsigned int history; int last_hit_triangle; }; @@ -29,6 +29,7 @@ struct Channel { int channel_id; double t; double q; + unsigned int mc_history; }; struct Event { @@ -67,17 +68,19 @@ void fill_photons(Event *ev, bool start, } -void fill_hits(Event *ev, unsigned int nchannels, float *channel_times) +void fill_hits(Event *ev, unsigned int nchannels, float *t, + float *q, unsigned int *history) { ev->channel.resize(0); ev->nhit = 0; Channel ch; for (unsigned int i=0; i < nchannels; i++) { - if (channel_times[i] < 1e8) { + if (t[i] < 1e8) { ev->nhit++; ch.channel_id = i; - ch.t = channel_times[i]; - ch.q = 1.0; + ch.t = t[i]; + ch.q = q[i]; + ch.mc_history = history[i]; ev->channel.push_back(ch); } } diff --git a/sim.py b/sim.py index 87993ab..c8f9fdb 100755 --- a/sim.py +++ b/sim.py @@ -72,9 +72,10 @@ def write_event(T, ev, event_id, hits, photon_start=None, photon_stop=None): photons['wavelength'], photons['t0'], photons['histories'], photons['last_hit_triangles']) - root.fill_hits(ev, len(hits), hits) + root.fill_hits(ev, len(hits), hits['t'], hits['q'], hits['history']) T.Fill() +#@profile def main(): parser = optparse.OptionParser('%prog') parser.add_option('-b', type='int', dest='nbits', default=10) @@ -159,7 +160,7 @@ def main(): gpu_worker.load_photons(**photons) gpu_worker.propagate() gpu_worker.run_daq() - hit_times = gpu_worker.get_hits() + hits = gpu_worker.get_hits() if options.save_photon_start: photon_start = photons @@ -170,7 +171,7 @@ def main(): else: photon_stop = None - write_event(T, ev, i, hit_times, + write_event(T, ev, i, hits, photon_start=photon_start, photon_stop=photon_stop) if i % 10 == 0: print >>sys.stderr, "\rEvent:", i, diff --git a/src/daq.cu b/src/daq.cu index 7c5e6a5..9824d7d 100644 --- a/src/daq.cu +++ b/src/daq.cu @@ -34,7 +34,8 @@ __global__ void run_daq(curandState *s, unsigned int detection_state, int nphotons, float *photon_times, unsigned int *photon_histories, int *last_hit_triangles, int *solid_map, - int nsolids, unsigned int *earliest_time_int) + int nsolids, unsigned int *earliest_time_int, + unsigned int *channel_histories) { int id = threadIdx.x + blockDim.x * blockIdx.x; @@ -48,7 +49,7 @@ __global__ void run_daq(curandState *s, unsigned int detection_state, if (triangle_id > -1) { int solid_id = solid_map[triangle_id]; - int history = photon_histories[photon_id]; + unsigned int history = photon_histories[photon_id]; if (solid_id < nsolids && (history & detection_state)) { @@ -56,6 +57,7 @@ __global__ void run_daq(curandState *s, unsigned int detection_state, unsigned int time_int = float_to_sortable_int(time); atomicMin(earliest_time_int + solid_id, time_int); + atomicOr(channel_histories + solid_id, history); } } -- cgit From 5b478fe72e600e06cd7b2e8a05a600f30c44d5c0 Mon Sep 17 00:00:00 2001 From: Stan Seibert Date: Tue, 9 Aug 2011 15:53:35 -0400 Subject: Put number of detected photons into charge value for channel. --- gpu.py | 9 ++++++--- src/daq.cu | 2 ++ 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/gpu.py b/gpu.py index 7b08473..e4eecbc 100644 --- a/gpu.py +++ b/gpu.py @@ -72,13 +72,13 @@ class GPU(object): device = cuda.Device(device_id) self.context = device.make_context() print 'device %s' % self.context.get_device().name() - self.module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True, cache_dir=False) + self.module = SourceModule(src.kernel, options=['-I' + src.dir], no_extern_c=True) self.geo_funcs = CUDAFuncs(self.module, ['set_wavelength_range', 'set_material', 'set_surface', 'set_global_mesh_variables', 'color_solids']) self.prop_funcs = CUDAFuncs(self.module, ['init_rng', 'propagate']) self.nthread_per_block = 64 self.max_blocks = 1024 - self.daq_module = SourceModule(src.daq, options=['-I' + src.dir], no_extern_c=True, cache_dir=False) + self.daq_module = SourceModule(src.daq, options=['-I' + src.dir], no_extern_c=True) self.daq_funcs = CUDAFuncs(self.daq_module, ['reset_earliest_time_int', 'run_daq', 'convert_sortable_int_to_float']) @@ -279,6 +279,7 @@ class GPU(object): self.earliest_time_int_gpu = gpuarray.GPUArray(shape=self.earliest_time_gpu.shape, dtype=np.uint32) self.channel_history_gpu = gpuarray.zeros_like(self.earliest_time_int_gpu) + self.channel_q_gpu = gpuarray.zeros_like(self.earliest_time_int_gpu) self.daq_pmt_rms = pmt_rms def run_daq(self): @@ -287,6 +288,7 @@ class GPU(object): self.earliest_time_int_gpu, block=(self.nthread_per_block,1,1), grid=(len(self.earliest_time_int_gpu)//self.nthread_per_block+1,1)) + self.channel_q_gpu.fill(0) self.channel_history_gpu.fill(0) #self.context.synchronize() @@ -300,6 +302,7 @@ class GPU(object): 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=(self.nthread_per_block,1,1), grid=(blocks,1)) #self.context.synchronize() @@ -313,7 +316,7 @@ class GPU(object): def get_hits(self): return { 't': self.earliest_time_gpu.get(), - 'q': np.ones(shape=len(self.earliest_time_int_gpu), dtype=np.float32), + 'q': self.channel_q_gpu.get().astype(np.float32), 'history': self.channel_history_gpu.get()} def __del__(self): diff --git a/src/daq.cu b/src/daq.cu index 9824d7d..2b95560 100644 --- a/src/daq.cu +++ b/src/daq.cu @@ -35,6 +35,7 @@ __global__ void run_daq(curandState *s, unsigned int detection_state, unsigned int *photon_histories, int *last_hit_triangles, int *solid_map, int nsolids, unsigned int *earliest_time_int, + unsigned int *channel_q, unsigned int *channel_histories) { @@ -57,6 +58,7 @@ __global__ void run_daq(curandState *s, unsigned int detection_state, unsigned int time_int = float_to_sortable_int(time); atomicMin(earliest_time_int + solid_id, time_int); + atomicAdd(channel_q + solid_id, 1); atomicOr(channel_histories + solid_id, history); } -- cgit