summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rwxr-xr-xchroma/benchmark.py38
-rw-r--r--chroma/cuda/daq.cu41
-rw-r--r--chroma/cuda/detector.h25
-rw-r--r--chroma/cuda/pdf.cu8
-rw-r--r--chroma/demo/__init__.py24
-rw-r--r--chroma/detector.py129
-rw-r--r--chroma/gpu/__init__.py1
-rw-r--r--chroma/gpu/daq.py30
-rw-r--r--chroma/gpu/detector.py40
-rw-r--r--chroma/gpu/pdf.py10
-rw-r--r--chroma/sim.py16
-rw-r--r--test/test_detector.py76
-rw-r--r--test/test_pdf.py6
-rw-r--r--test/test_propagation.py1
-rw-r--r--test/test_rayleigh.py1
15 files changed, 379 insertions, 67 deletions
diff --git a/chroma/benchmark.py b/chroma/benchmark.py
index af422dd..44cf44b 100755
--- a/chroma/benchmark.py
+++ b/chroma/benchmark.py
@@ -64,7 +64,7 @@ def load_photons(number=100, nphotons=500000):
return nphotons/ufloat((np.mean(run_times),np.std(run_times)))
-def propagate(gpu_geometry, number=10, nphotons=500000, nthreads_per_block=64,
+def propagate(gpu_detector, number=10, nphotons=500000, nthreads_per_block=64,
max_blocks=1024):
"Returns the average number of photons propagated on the GPU per second."
rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks)
@@ -79,7 +79,7 @@ def propagate(gpu_geometry, number=10, nphotons=500000, nthreads_per_block=64,
gpu_photons = gpu.GPUPhotons(photons)
t0 = time.time()
- gpu_photons.propagate(gpu_geometry, rng_states, nthreads_per_block,
+ gpu_photons.propagate(gpu_detector, rng_states, nthreads_per_block,
max_blocks)
cuda.Context.get_current().synchronize()
elapsed = time.time() - t0
@@ -91,7 +91,7 @@ def propagate(gpu_geometry, number=10, nphotons=500000, nthreads_per_block=64,
return nphotons/ufloat((np.mean(run_times),np.std(run_times)))
@tools.profile_if_possible
-def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
+def pdf(gpu_detector, npdfs=10, nevents=100, nreps=16, ndaq=1,
nthreads_per_block=64, max_blocks=1024):
"""
Returns the average number of 100 MeV events per second that can be
@@ -99,9 +99,7 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
Args:
- gpu_instance, chroma.gpu.GPU
- The GPU instance passed to the GPUGeometry constructor.
- - max_pmt_id, int
- The channel number of the highest PMT
+ The GPU instance passed to the GPUDetector constructor.
- npdfs, int
The number of pdf generations to average.
- nevents, int
@@ -114,9 +112,9 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
"""
rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks)
- gpu_daq = gpu.GPUDaq(gpu_geometry, max_pmt_id)
+ gpu_daq = gpu.GPUDaq(gpu_detector)
gpu_pdf = gpu.GPUPDF()
- gpu_pdf.setup_pdf(max_pmt_id, 100, (-0.5, 999.5), 10, (-0.5, 9.5))
+ gpu_pdf.setup_pdf(gpu_detector.nchannels, 100, (-0.5, 999.5), 10, (-0.5, 9.5))
run_times = []
for i in tools.progress(range(npdfs)):
@@ -130,7 +128,7 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
for ev in g4generator.generate_events(vertex_iter):
gpu_photons = gpu.GPUPhotons(ev.photons_beg, ncopies=nreps)
- gpu_photons.propagate(gpu_geometry, rng_states,
+ gpu_photons.propagate(gpu_detector, rng_states,
nthreads_per_block, max_blocks)
for gpu_photon_slice in gpu_photons.iterate_copies():
for idaq in xrange(ndaq):
@@ -148,7 +146,7 @@ def pdf(gpu_geometry, max_pmt_id, npdfs=10, nevents=100, nreps=16, ndaq=1,
return nevents*nreps*ndaq/ufloat((np.mean(run_times),np.std(run_times)))
-def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
+def pdf_eval(gpu_detector, npdfs=10, nevents=25, nreps=16, ndaq=128,
nthreads_per_block=64, max_blocks=1024):
"""
Returns the average number of 100 MeV events that can be
@@ -156,9 +154,7 @@ def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
Args:
- gpu_instance, chroma.gpu.GPU
- The GPU instance passed to the GPUGeometry constructor.
- - max_pmt_id, int
- The channel number of the highest PMT
+ The GPU instance passed to the GPUDetector constructor.
- npdfs, int
The number of pdf generations to average.
- nevents, int
@@ -177,9 +173,9 @@ def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
1)).next()
gpu_photons = gpu.GPUPhotons(data_ev.photons_beg)
- gpu_photons.propagate(gpu_geometry, rng_states,
+ gpu_photons.propagate(gpu_detector, rng_states,
nthreads_per_block, max_blocks)
- gpu_daq = gpu.GPUDaq(gpu_geometry, max_pmt_id)
+ gpu_daq = gpu.GPUDaq(gpu_detector)
data_ev_channels = gpu_daq.acquire(gpu_photons, rng_states, nthreads_per_block, max_blocks).get()
# Setup PDF evaluation
@@ -206,7 +202,7 @@ def pdf_eval(gpu_geometry, max_pmt_id, npdfs=10, nevents=25, nreps=16, ndaq=128,
for ev in g4generator.generate_events(vertex_iter):
gpu_photons = gpu.GPUPhotons(ev.photons_beg, ncopies=nreps)
- gpu_photons.propagate(gpu_geometry, rng_states,
+ gpu_photons.propagate(gpu_detector, rng_states,
nthreads_per_block, max_blocks)
for gpu_photon_slice in gpu_photons.iterate_copies():
for idaq in xrange(ndaq):
@@ -237,11 +233,11 @@ if __name__ == '__main__':
detector.build(bits=11)
context = gpu.create_cuda_context()
- gpu_geometry = gpu.GPUGeometry(detector)
+ gpu_detector = gpu.GPUDetector(detector)
if 'ray' in tests:
print '%s ray intersections/sec.' % \
- tools.ufloat_to_str(intersect(gpu_geometry))
+ tools.ufloat_to_str(intersect(gpu_detector))
# run garbage collection since there is a reference loop
# in the GPUArray class.
gc.collect()
@@ -252,16 +248,16 @@ if __name__ == '__main__':
if 'propagate' in tests:
print '%s photons propagated/sec.' % \
- tools.ufloat_to_str(propagate(gpu_geometry))
+ tools.ufloat_to_str(propagate(gpu_detector))
gc.collect()
if 'pdf' in tests:
print '%s 100 MeV events histogrammed/s' % \
- tools.ufloat_to_str(pdf(gpu_geometry, max(detector.pmtids)))
+ tools.ufloat_to_str(pdf(gpu_detector))
gc.collect()
if 'pdf_eval' in tests:
print '%s 100 MeV events/s accumulated in PDF evaluation data structure (100 GEANT4 x 16 Chroma x 128 DAQ)' % \
- tools.ufloat_to_str(pdf_eval(gpu_geometry, max(detector.pmtids)))
+ tools.ufloat_to_str(pdf_eval(gpu_detector))
context.pop()
diff --git a/chroma/cuda/daq.cu b/chroma/cuda/daq.cu
index f80fcaa..b995bfb 100644
--- a/chroma/cuda/daq.cu
+++ b/chroma/cuda/daq.cu
@@ -1,5 +1,6 @@
// -*-c++-*-
-#include <curand_kernel.h>
+#include "detector.h"
+#include "random.h"
__device__ unsigned int
float_to_sortable_int(float f)
@@ -32,11 +33,13 @@ reset_earliest_time_int(float maxtime, int ntime_ints, unsigned int *time_ints)
}
__global__ void
-run_daq(curandState *s, unsigned int detection_state, float time_rms,
+run_daq(curandState *s, unsigned int detection_state,
int first_photon, int nphotons, float *photon_times,
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)
+ int *solid_map,
+ Detector *detector,
+ unsigned int *earliest_time_int,
+ unsigned int *channel_q_int, unsigned int *channel_histories)
{
int id = threadIdx.x + blockDim.x * blockIdx.x;
@@ -49,14 +52,22 @@ run_daq(curandState *s, unsigned int detection_state, float time_rms,
if (triangle_id > -1) {
int solid_id = solid_map[triangle_id];
unsigned int history = photon_histories[photon_id];
+ int channel_index = detector->solid_id_to_channel_index[solid_id];
- if (solid_id < nsolids && (history & detection_state)) {
- float time = photon_times[photon_id] + curand_normal(&rng) * time_rms;
+ if (channel_index >= 0 && (history & detection_state)) {
+ float time = photon_times[photon_id] +
+ sample_cdf(&rng, detector->time_cdf_len,
+ detector->time_cdf_x, detector->time_cdf_y);
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);
+ float charge = sample_cdf(&rng, detector->charge_cdf_len,
+ detector->charge_cdf_x,
+ detector->charge_cdf_y);
+ unsigned int charge_int = roundf(charge / detector->charge_unit);
+
+ atomicMin(earliest_time_int + channel_index, time_int);
+ atomicAdd(channel_q_int + channel_index, charge_int);
+ atomicOr(channel_histories + channel_index, history);
}
}
@@ -78,4 +89,16 @@ convert_sortable_int_to_float(int n, unsigned int *sortable_ints,
}
+__global__ void
+convert_charge_int_to_float(Detector *detector,
+ unsigned int *charge_int,
+ float *charge_float)
+{
+ int id = threadIdx.x + blockDim.x * blockIdx.x;
+
+ if (id < detector->nchannels)
+ charge_float[id] = charge_int[id] * detector->charge_unit;
+}
+
+
} // extern "C"
diff --git a/chroma/cuda/detector.h b/chroma/cuda/detector.h
new file mode 100644
index 0000000..16bd164
--- /dev/null
+++ b/chroma/cuda/detector.h
@@ -0,0 +1,25 @@
+#ifndef __DETECTOR_H__
+#define __DETECTOR_H__
+
+struct Detector
+{
+ // Order in decreasing size to avoid alignment problems
+ int *solid_id_to_channel_index;
+
+ float *time_cdf_x;
+ float *time_cdf_y;
+
+ float *charge_cdf_x;
+ float *charge_cdf_y;
+
+ int nchannels;
+ int time_cdf_len;
+ int charge_cdf_len;
+ float charge_unit;
+ // Convert charges to/from quantized integers with
+ // q_int = (int) roundf(q / charge_unit )
+ // q = q_int * charge_unit
+};
+
+
+#endif // __DETECTOR_H__
diff --git a/chroma/cuda/pdf.cu b/chroma/cuda/pdf.cu
index 9f547d0..0d82e3a 100644
--- a/chroma/cuda/pdf.cu
+++ b/chroma/cuda/pdf.cu
@@ -5,7 +5,7 @@ extern "C"
{
__global__ void
-bin_hits(int nchannels, unsigned int *channel_q, float *channel_time,
+bin_hits(int nchannels, float *channel_q, float *channel_time,
unsigned int *hitcount, int tbins, float tmin, float tmax, int qbins,
float qmin, float qmax, unsigned int *pdf)
{
@@ -32,7 +32,7 @@ bin_hits(int nchannels, unsigned int *channel_q, float *channel_time,
__global__ void
accumulate_pdf_eval(int time_only, int nchannels, unsigned int *event_hit,
float *event_time, float *event_charge, float *mc_time,
- unsigned int *mc_charge, // quirk of DAQ!
+ float *mc_charge, // quirk of DAQ!
unsigned int *hitcount, unsigned int *bincount,
float min_twidth, float tmin, float tmax,
float min_qwidth, float qmin, float qmax,
@@ -126,7 +126,7 @@ accumulate_pdf_eval(int time_only, int nchannels, unsigned int *event_hit,
__global__ void
accumulate_moments(int time_only, int nchannels,
float *mc_time,
- unsigned int *mc_charge, // quirk of DAQ!
+ float *mc_charge,
float tmin, float tmax,
float qmin, float qmax,
unsigned int *mom0,
@@ -174,7 +174,7 @@ static const float rootPiBy2 = 1.2533141373155001f; // sqrt(M_PI/2)
__global__ void
accumulate_kernel_eval(int time_only, int nchannels, unsigned int *event_hit,
float *event_time, float *event_charge, float *mc_time,
- unsigned int *mc_charge, // quirk of DAQ!
+ float *mc_charge,
float tmin, float tmax,
float qmin, float qmax,
float *inv_time_bandwidths,
diff --git a/chroma/demo/__init__.py b/chroma/demo/__init__.py
index 7df7df7..46f4529 100644
--- a/chroma/demo/__init__.py
+++ b/chroma/demo/__init__.py
@@ -6,13 +6,15 @@ import numpy.linalg
from chroma.make import sphere
from chroma.stl import mesh_from_stl
-from chroma.geometry import *
+from chroma.geometry import Solid
+from chroma.detector import Detector
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, black_surface
from chroma.demo.checkerboard import build_checkerboard_scene as checkerboard_scene
+from chroma.log import logger
def spherical_spiral(radius, spacing):
'''Returns iterator generating points on a spiral wrapping the
@@ -29,16 +31,13 @@ def spherical_spiral(radius, spacing):
def detector(pmt_radius=14000.0, sphere_radius=14500.0, spiral_step=350.0):
pmt = build_8inch_pmt_with_lc()
- geo = Geometry(water)
+ geo = Detector(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 = -normalize(position)
@@ -49,10 +48,19 @@ def detector(pmt_radius=14000.0, sphere_radius=14500.0, spiral_step=350.0):
rotation = make_rotation_matrix(angle, axis)
# Place PMT (note that position is front face of PMT)
- chroma_id = geo.add_solid(pmt, rotation, position)
- geo.pmtids.append(chroma_id)
+ geo.add_pmt(pmt, rotation, position)
- print 'Demo detector: %d PMTs' % len(geo.pmtids)
+
+ time_rms = 1.5 # ns
+ charge_mean = 1.0
+ charge_rms = 0.1 # Don't I wish!
+
+ geo.set_time_dist_gaussian(time_rms, -5 * time_rms, 5*time_rms)
+ geo.set_charge_dist_gaussian(charge_mean, charge_rms, 0.0, charge_mean + 5*charge_rms)
+
+ logger.info('Demo detector: %d PMTs' % geo.num_channels())
+ logger.info(' %1.1f ns time RMS' % time_rms)
+ logger.info(' %1.1f%% charge RMS' % (100.0*charge_rms/charge_mean))
return geo
def tiny():
diff --git a/chroma/detector.py b/chroma/detector.py
new file mode 100644
index 0000000..f3f2d3b
--- /dev/null
+++ b/chroma/detector.py
@@ -0,0 +1,129 @@
+import numpy as np
+
+from chroma.geometry import Geometry
+
+class Detector(Geometry):
+ '''A Detector is a subclass of Geometry that allows some Solids
+ to be marked as photon detectors, which we will suggestively call
+ "PMTs." Each detector is imagined to be connected to an electronics
+ channel that records a hit time and charge.
+
+ Each PMT has two integers identifying it: a channel index and a
+ channel ID. When all of the channels in the detector are stored
+ in a Numpy array, they will be stored in index order. Channel
+ indices star from zero and have no gaps. Channel ID numbers are
+ arbitrary integers that identify a PMT uniquely in a stable way.
+ They are written out to disk when using the Chroma ROOT format,
+ and are used when reading events back in to map channels back
+ into the correct array index.
+
+ For now, all the PMTs share a single set of time and charge
+ distributions. In the future, this will be generalized to
+ allow per-channel distributions.
+ '''
+
+ def __init__(self, detector_material=None):
+ Geometry.__init__(self, detector_material=detector_material)
+
+ # Using numpy arrays here to allow for fancy indexing
+ self.solid_id_to_channel_index = np.zeros(0, dtype=np.int32)
+ self.channel_index_to_solid_id = np.zeros(0, dtype=np.int32)
+
+ self.channel_index_to_channel_id = np.zeros(0, dtype=np.int32)
+
+ # If the ID numbers are arbitrary, we can't treat them
+ # as array indices, so have to use a dictionary
+ self.channel_id_to_channel_index = {}
+
+ # zero time and unit charge distributions
+ self.time_cdf = (np.array([-0.01, 0.01]), np.array([0.0, 1.0]))
+ self.charge_cdf = (np.array([0.99, 1.00]), np.array([0.0, 1.0]))
+
+ def add_solid(self, solid, rotation=None, displacement=None):
+ """
+ Add the solid `solid` to the geometry. When building the final triangle
+ mesh, `solid` will be placed by rotating it with the rotation matrix
+ `rotation` and displacing it by the vector `displacement`.
+ """
+ solid_id = Geometry.add_solid(self, solid=solid, rotation=rotation,
+ displacement=displacement)
+ self.solid_id_to_channel_index.resize(solid_id+1)
+ self.solid_id_to_channel_index[solid_id] = -1 # solid maps to no channel
+ return solid_id
+
+ def add_pmt(self, pmt, rotation=None, displacement=None, channel_id=None):
+ """Add the PMT `pmt` to the geometry. When building the final triangle
+ mesh, `solid` will be placed by rotating it with the rotation matrix
+ `rotation` and displacing it by the vector `displacement`, just like
+ add_solid().
+
+ `pmt``: instance of chroma.Solid
+ Solid representing a PMT.
+ `rotation`: numpy.matrix (3x3)
+ Rotation to apply to PMT mesh before displacement. Defaults to
+ identity rotation.
+ `displacement`: numpy.ndarray (shape=3)
+ 3-vector displacement to apply to PMT mesh after rotation.
+ Defaults to zero vector.
+ `channel_id`: int
+ Integer identifier for this PMT. May be any integer, with no
+ requirement for consective numbering. Defaults to None,
+ where the ID number will be set to the generated channel index.
+ The channel_id must be representable as a 32-bit integer.
+
+ Returns: dictionary { 'solid_id' : solid_id,
+ 'channel_index' : channel_index,
+ 'channel_id' : channel_id }
+ """
+
+ solid_id = self.add_solid(solid=pmt, rotation=rotation,
+ displacement=displacement)
+
+ channel_index = len(self.channel_index_to_solid_id)
+ if channel_id is None:
+ channel_id = channel_index
+
+ # add_solid resized this array already
+ self.solid_id_to_channel_index[solid_id] = channel_index
+
+ # resize channel_index arrays before filling
+ self.channel_index_to_solid_id.resize(channel_index+1)
+ self.channel_index_to_solid_id[channel_index] = solid_id
+ self.channel_index_to_channel_id.resize(channel_index+1)
+ self.channel_index_to_channel_id[channel_index] = channel_id
+
+ # dictionary does not need resizing
+ self.channel_id_to_channel_index[channel_id] = channel_index
+
+ return { 'solid_id' : solid_id,
+ 'channel_index' : channel_index,
+ 'channel_id' : channel_id }
+
+ def _pdf_to_cdf(self, bin_edges, bin_contents):
+ '''Returns tuple of arrays `(cdf_x, cdf_y)` corresponding to
+ the cumulative distribution function for a binned PDF with the
+ given bin_edges and bin_contents.
+
+ `bin_edges`: array
+ bin edges of PDF. length is len(bin_contents) + 1
+ `bin_contents`: array
+ The contents of each bin. The array should NOT be
+ normalized for bin width if variable binning is present.
+ '''
+ cdf_x = np.copy(bin_edges)
+ cdf_y = np.array([0.0]+bin_contents.cumsum())
+ cdf_y /= cdf_y[-1] # normalize
+ return (cdf_x, cdf_y)
+
+ def set_time_dist_gaussian(self, rms, lo, hi, nsamples=50):
+ pdf_x = np.linspace(lo, hi, nsamples+1, endpoint=True)
+ pdf_y = np.exp(-0.5 * (pdf_x[1:]/rms)**2)
+ self.time_cdf = self._pdf_to_cdf(pdf_x, pdf_y)
+
+ def set_charge_dist_gaussian(self, mean, rms, lo, hi, nsamples=50):
+ pdf_x = np.linspace(lo, hi, nsamples+1, endpoint=True)
+ pdf_y = np.exp(-0.5 * ((pdf_x[1:] - mean)/rms)**2)
+ self.charge_cdf = self._pdf_to_cdf(pdf_x, pdf_y)
+
+ def num_channels(self):
+ return len(self.channel_index_to_channel_id)
diff --git a/chroma/gpu/__init__.py b/chroma/gpu/__init__.py
index b5a1ff3..74dc9f2 100644
--- a/chroma/gpu/__init__.py
+++ b/chroma/gpu/__init__.py
@@ -4,3 +4,4 @@ from chroma.gpu.render import *
from chroma.gpu.photon import *
from chroma.gpu.daq import *
from chroma.gpu.pdf import *
+from chroma.gpu.detector import *
diff --git a/chroma/gpu/daq.py b/chroma/gpu/daq.py
index 6ec1a79..571feba 100644
--- a/chroma/gpu/daq.py
+++ b/chroma/gpu/daq.py
@@ -13,27 +13,30 @@ class GPUChannels(object):
def get(self):
t = self.t.get()
- q = self.q.get().astype(np.float32)
+ q = self.q.get()
# For now, assume all channels with small
# enough hit time were hit.
return event.Channels(t<1e8, t, q, self.flags.get())
class GPUDaq(object):
- 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)
+ def __init__(self, gpu_detector):
+ self.earliest_time_gpu = ga.empty(gpu_detector.nchannels, dtype=np.float32)
+ self.earliest_time_int_gpu = ga.empty(gpu_detector.nchannels, dtype=np.uint32)
self.channel_history_gpu = ga.zeros_like(self.earliest_time_int_gpu)
- self.channel_q_gpu = ga.zeros_like(self.earliest_time_int_gpu)
- self.daq_pmt_rms = pmt_rms
- self.solid_id_map_gpu = gpu_geometry.solid_id_map
+ self.channel_q_int_gpu = ga.zeros_like(self.earliest_time_int_gpu)
+ self.channel_q_gpu = ga.zeros(len(self.earliest_time_int_gpu), dtype=np.float32)
+ self.detector_gpu = gpu_detector.detector_gpu
+ self.solid_id_map_gpu = gpu_detector.solid_id_map
+ self.solid_id_to_channel_index_gpu = gpu_detector.solid_id_to_channel_index_gpu
self.module = get_cu_module('daq.cu', options=cuda_options,
- include_source_directory=False)
+ include_source_directory=True)
self.gpu_funcs = GPUFuncs(self.module)
def acquire(self, gpuphotons, rng_states, nthreads_per_block=64, max_blocks=1024):
self.gpu_funcs.reset_earliest_time_int(np.float32(1e9), np.int32(len(self.earliest_time_int_gpu)), self.earliest_time_int_gpu, block=(nthreads_per_block,1,1), grid=(len(self.earliest_time_int_gpu)//nthreads_per_block+1,1))
+ self.channel_q_int_gpu.fill(0)
self.channel_q_gpu.fill(0)
self.channel_history_gpu.fill(0)
@@ -41,9 +44,18 @@ class GPUDaq(object):
for first_photon, photons_this_round, blocks in \
chunk_iterator(n, nthreads_per_block, max_blocks):
- self.gpu_funcs.run_daq(rng_states, np.uint32(0x1 << 2), np.float32(self.daq_pmt_rms), np.int32(first_photon), np.int32(photons_this_round), gpuphotons.t, gpuphotons.flags, gpuphotons.last_hit_triangles, 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=(nthreads_per_block,1,1), grid=(blocks,1))
+ self.gpu_funcs.run_daq(rng_states, np.uint32(0x1 << 2),
+ np.int32(first_photon), np.int32(photons_this_round), gpuphotons.t,
+ gpuphotons.flags, gpuphotons.last_hit_triangles,
+ self.solid_id_map_gpu,
+ self.detector_gpu,
+ self.earliest_time_int_gpu,
+ self.channel_q_int_gpu, self.channel_history_gpu,
+ block=(nthreads_per_block,1,1), grid=(blocks,1))
self.gpu_funcs.convert_sortable_int_to_float(np.int32(len(self.earliest_time_int_gpu)), self.earliest_time_int_gpu, self.earliest_time_gpu, block=(nthreads_per_block,1,1), grid=(len(self.earliest_time_int_gpu)//nthreads_per_block+1,1))
+ self.gpu_funcs.convert_charge_int_to_float(self.detector_gpu, self.channel_q_int_gpu, self.channel_q_gpu, block=(nthreads_per_block,1,1), grid=(len(self.channel_q_int_gpu)//nthreads_per_block+1,1))
+
return GPUChannels(self.earliest_time_gpu, self.channel_q_gpu, self.channel_history_gpu)
diff --git a/chroma/gpu/detector.py b/chroma/gpu/detector.py
new file mode 100644
index 0000000..ac0092e
--- /dev/null
+++ b/chroma/gpu/detector.py
@@ -0,0 +1,40 @@
+import numpy as np
+import pycuda.driver as cuda
+from pycuda import gpuarray as ga
+from pycuda import characterize
+
+from chroma.geometry import standard_wavelengths
+from chroma.gpu.tools import get_cu_module, get_cu_source, cuda_options, \
+ chunk_iterator, format_array, format_size, to_uint3, to_float3, \
+ make_gpu_struct
+from chroma.log import logger
+
+from chroma.gpu.geometry import GPUGeometry
+
+class GPUDetector(GPUGeometry):
+ def __init__(self, detector, wavelengths=None, print_usage=False):
+ GPUGeometry.__init__(self, detector, wavelengths=wavelengths, print_usage=False)
+ self.solid_id_to_channel_index_gpu = \
+ ga.to_gpu(detector.solid_id_to_channel_index.astype(np.int32))
+ self.nchannels = detector.num_channels()
+
+
+ self.time_cdf_x_gpu = ga.to_gpu(detector.time_cdf[0].astype(np.float32))
+ self.time_cdf_y_gpu = ga.to_gpu(detector.time_cdf[1].astype(np.float32))
+
+ self.charge_cdf_x_gpu = ga.to_gpu(detector.charge_cdf[0].astype(np.float32))
+ self.charge_cdf_y_gpu = ga.to_gpu(detector.charge_cdf[1].astype(np.float32))
+
+ detector_source = get_cu_source('detector.h')
+ detector_struct_size = characterize.sizeof('Detector', detector_source)
+ self.detector_gpu = make_gpu_struct(detector_struct_size,
+ [self.solid_id_to_channel_index_gpu,
+ self.time_cdf_x_gpu,
+ self.time_cdf_y_gpu,
+ self.charge_cdf_x_gpu,
+ self.charge_cdf_y_gpu,
+ np.int32(self.nchannels),
+ np.int32(len(detector.time_cdf[0])),
+ np.int32(len(detector.charge_cdf[0])),
+ np.float32(detector.charge_cdf[0][-1] / 2**16)])
+
diff --git a/chroma/gpu/pdf.py b/chroma/gpu/pdf.py
index ba0d2af..5fb4b1a 100644
--- a/chroma/gpu/pdf.py
+++ b/chroma/gpu/pdf.py
@@ -76,7 +76,7 @@ class GPUKernelPDF(object):
d = 2
dimensionality_factor = ((4.0/(d+2)) / (mom0/scale_factor))**(-1.0/(d+4))
- time_bandwidths = dimensionality_factor * trms
+ time_bandwidths = dimensionality_factor * trms
inv_time_bandwidths = np.zeros_like(time_bandwidths)
inv_time_bandwidths[time_bandwidths > 0] = time_bandwidths[time_bandwidths > 0] ** -1
#inv_time_bandwidths /= 4.0
@@ -164,18 +164,18 @@ class GPUPDF(object):
include_source_directory=False)
self.gpu_funcs = GPUFuncs(self.module)
- def setup_pdf(self, max_pmt_id, tbins, trange, qbins, qrange):
+ def setup_pdf(self, nchannels, tbins, trange, qbins, qrange):
"""Setup GPU arrays to hold PDF information.
- max_pmt_id: int, largest PMT id #
+ nchannels: int, number of channels
tbins: number of time bins
trange: tuple of (min, max) time in PDF
qbins: number of charge bins
qrange: tuple of (min, max) charge in PDF
"""
self.events_in_histogram = 0
- self.hitcount_gpu = ga.zeros(max_pmt_id+1, dtype=np.uint32)
- self.pdf_gpu = ga.zeros(shape=(max_pmt_id+1, tbins, qbins),
+ self.hitcount_gpu = ga.zeros(nchannels, dtype=np.uint32)
+ self.pdf_gpu = ga.zeros(shape=(nchannels, tbins, qbins),
dtype=np.uint32)
self.tbins = tbins
self.trange = trange
diff --git a/chroma/sim.py b/chroma/sim.py
index d15c02d..f00f536 100644
--- a/chroma/sim.py
+++ b/chroma/sim.py
@@ -41,10 +41,13 @@ class Simulation(object):
# need to build geometry
detector.build()
- self.gpu_geometry = gpu.GPUGeometry(detector)
- self.gpu_daq = gpu.GPUDaq(self.gpu_geometry, max(self.detector.pmtids))
- self.gpu_pdf = gpu.GPUPDF()
- self.gpu_pdf_kernel = gpu.GPUKernelPDF()
+ if hasattr(detector, 'num_channels'):
+ self.gpu_geometry = gpu.GPUDetector(detector)
+ self.gpu_daq = gpu.GPUDaq(self.gpu_geometry)
+ self.gpu_pdf = gpu.GPUPDF()
+ self.gpu_pdf_kernel = gpu.GPUKernelPDF()
+ else:
+ self.gpu_geometry = gpu.GPUGeometry(detector)
self.rng_states = gpu.get_rng_states(self.nthreads_per_block*self.max_blocks, seed=self.seed)
@@ -81,7 +84,8 @@ class Simulation(object):
if keep_photons_end:
ev.photons_end = gpu_photons.get()
- if run_daq:
+ # Skip running DAQ if we don't have one
+ if hasattr(self, 'gpu_daq') and run_daq:
gpu_channels = self.gpu_daq.acquire(gpu_photons, self.rng_states, nthreads_per_block=self.nthreads_per_block, max_blocks=self.max_blocks)
ev.channels = gpu_channels.get()
@@ -98,7 +102,7 @@ class Simulation(object):
pdf_config = (tbins, trange, qbins, qrange)
if pdf_config != self.pdf_config:
self.pdf_config = pdf_config
- self.gpu_pdf.setup_pdf(max(self.detector.pmtids), tbins, trange,
+ self.gpu_pdf.setup_pdf(self.detector.num_channels(), tbins, trange,
qbins, qrange)
else:
self.gpu_pdf.clear_pdf()
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)
generated by cgit (git 2.34.1) at 2025-05-07 22:22:02 +0000