diff options
| -rw-r--r-- | chroma/cuda/daq.cu | 117 | ||||
| -rw-r--r-- | chroma/cuda/pdf.cu | 125 | ||||
| -rw-r--r-- | chroma/gpu/pdf.py | 4 |
3 files changed, 127 insertions, 119 deletions
diff --git a/chroma/cuda/daq.cu b/chroma/cuda/daq.cu index 5a68846..f80fcaa 100644 --- a/chroma/cuda/daq.cu +++ b/chroma/cuda/daq.cu @@ -77,122 +77,5 @@ convert_sortable_int_to_float(int n, unsigned int *sortable_ints, float_output[id] = sortable_int_to_float(sortable_ints[id]); } -__global__ void -bin_hits(int nchannels, unsigned int *channel_q, float *channel_time, - unsigned int *hitcount, int tbins, float tmin, float tmax, int qbins, - float qmin, float qmax, unsigned int *pdf) -{ - int id = threadIdx.x + blockDim.x * blockIdx.x; - - if (id >= nchannels) - return; - - unsigned int q = channel_q[id]; - float t = channel_time[id]; - - if (t < 1e8 && t >= tmin && t < tmax && q >= qmin && q < qmax) { - hitcount[id] += 1; - - int tbin = (t - tmin) / (tmax - tmin) * tbins; - int qbin = (q - qmin) / (qmax - qmin) * qbins; - - // row major order (channel, t, q) - int bin = id * (tbins * qbins) + tbin * qbins + qbin; - pdf[bin] += 1; - } -} - -__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! - unsigned int *hitcount, unsigned int *bincount, - float min_twidth, float tmin, float tmax, - float min_qwidth, float qmin, float qmax, - float *nearest_mc, int min_bin_content) -{ - int id = threadIdx.x + blockDim.x * blockIdx.x; - - if (id >= nchannels) - return; - - // Was this channel hit in the Monte Carlo? - float channel_mc_time = mc_time[id]; - if (channel_mc_time >= 1e8) - return; // Nothing else to do - - // Is this channel inside the range of the PDF? - float distance; - int channel_bincount = bincount[id]; - if (time_only) { - if (channel_mc_time < tmin || channel_mc_time > tmax) - return; // Nothing else to do - - // This MC information is contained in the PDF - hitcount[id] += 1; - - // Was this channel hit in the event-of-interest? - int channel_event_hit = event_hit[id]; - if (!channel_event_hit) - return; // No need to update PDF value for unhit channel - - // Are we inside the minimum size bin? - float channel_event_time = event_time[id]; - distance = fabsf(channel_mc_time - channel_event_time); - if (distance < min_twidth/2.0) - bincount[id] = channel_bincount + 1; - - } - else { // time and charge PDF - float channel_mc_charge = mc_charge[id]; // int->float conversion because DAQ just returns an integer - - if (channel_mc_time < tmin || channel_mc_time > tmax || - channel_mc_charge < qmin || channel_mc_charge > qmax) - return; // Nothing else to do - - // This MC information is contained in the PDF - hitcount[id] += 1; - - // Was this channel hit in the event-of-interest? - int channel_event_hit = event_hit[id]; - if (!channel_event_hit) - return; // No need to update PDF value for unhit channel - - // Are we inside the minimum size bin? - float channel_event_time = event_time[id]; - float channel_event_charge = event_charge[id]; - float normalized_time_distance = fabsf(channel_event_time - channel_mc_time)/min_twidth/2.0; - float normalized_charge_distance = fabsf(channel_event_charge - channel_mc_charge)/min_qwidth/2.0; - distance = sqrt(normalized_time_distance*normalized_time_distance + normalized_charge_distance*normalized_charge_distance); - - if (distance < 1.0f) - bincount[id] = channel_bincount + 1; - } - - // Do we need to keep updating the nearest_mc list? - if (channel_bincount + 1 >= min_bin_content) - return; // No need to perform insertion into nearest_mc because bincount is a better estimate of the PDF value - - // insertion sort the distance into the array of nearest MC points - int offset = min_bin_content * id; - int entry = min_bin_content - 1; - - // If last entry less than new entry, nothing to update - if (distance > nearest_mc[offset + entry]) - return; - - // Find where to insert the new entry while sliding the rest - // to the right - entry--; - while (entry >= 0) { - if (nearest_mc[offset+entry] >= distance) - nearest_mc[offset+entry+1] = nearest_mc[offset+entry]; - else - break; - entry--; - } - - nearest_mc[offset+entry+1] = distance; -} } // extern "C" diff --git a/chroma/cuda/pdf.cu b/chroma/cuda/pdf.cu new file mode 100644 index 0000000..c9c2ea4 --- /dev/null +++ b/chroma/cuda/pdf.cu @@ -0,0 +1,125 @@ +// -*-c++-*- +#include <curand_kernel.h> + +extern "C" +{ + +__global__ void +bin_hits(int nchannels, unsigned int *channel_q, float *channel_time, + unsigned int *hitcount, int tbins, float tmin, float tmax, int qbins, + float qmin, float qmax, unsigned int *pdf) +{ + int id = threadIdx.x + blockDim.x * blockIdx.x; + + if (id >= nchannels) + return; + + unsigned int q = channel_q[id]; + float t = channel_time[id]; + + if (t < 1e8 && t >= tmin && t < tmax && q >= qmin && q < qmax) { + hitcount[id] += 1; + + int tbin = (t - tmin) / (tmax - tmin) * tbins; + int qbin = (q - qmin) / (qmax - qmin) * qbins; + + // row major order (channel, t, q) + int bin = id * (tbins * qbins) + tbin * qbins + qbin; + pdf[bin] += 1; + } +} + +__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! + unsigned int *hitcount, unsigned int *bincount, + float min_twidth, float tmin, float tmax, + float min_qwidth, float qmin, float qmax, + float *nearest_mc, int min_bin_content) +{ + int id = threadIdx.x + blockDim.x * blockIdx.x; + + if (id >= nchannels) + return; + + // Was this channel hit in the Monte Carlo? + float channel_mc_time = mc_time[id]; + if (channel_mc_time >= 1e8) + return; // Nothing else to do + + // Is this channel inside the range of the PDF? + float distance; + int channel_bincount = bincount[id]; + if (time_only) { + if (channel_mc_time < tmin || channel_mc_time > tmax) + return; // Nothing else to do + + // This MC information is contained in the PDF + hitcount[id] += 1; + + // Was this channel hit in the event-of-interest? + int channel_event_hit = event_hit[id]; + if (!channel_event_hit) + return; // No need to update PDF value for unhit channel + + // Are we inside the minimum size bin? + float channel_event_time = event_time[id]; + distance = fabsf(channel_mc_time - channel_event_time); + if (distance < min_twidth/2.0) + bincount[id] = channel_bincount + 1; + + } + else { // time and charge PDF + float channel_mc_charge = mc_charge[id]; // int->float conversion because DAQ just returns an integer + + if (channel_mc_time < tmin || channel_mc_time > tmax || + channel_mc_charge < qmin || channel_mc_charge > qmax) + return; // Nothing else to do + + // This MC information is contained in the PDF + hitcount[id] += 1; + + // Was this channel hit in the event-of-interest? + int channel_event_hit = event_hit[id]; + if (!channel_event_hit) + return; // No need to update PDF value for unhit channel + + // Are we inside the minimum size bin? + float channel_event_time = event_time[id]; + float channel_event_charge = event_charge[id]; + float normalized_time_distance = fabsf(channel_event_time - channel_mc_time)/min_twidth/2.0; + float normalized_charge_distance = fabsf(channel_event_charge - channel_mc_charge)/min_qwidth/2.0; + distance = sqrt(normalized_time_distance*normalized_time_distance + normalized_charge_distance*normalized_charge_distance); + + if (distance < 1.0f) + bincount[id] = channel_bincount + 1; + } + + // Do we need to keep updating the nearest_mc list? + if (channel_bincount + 1 >= min_bin_content) + return; // No need to perform insertion into nearest_mc because bincount is a better estimate of the PDF value + + // insertion sort the distance into the array of nearest MC points + int offset = min_bin_content * id; + int entry = min_bin_content - 1; + + // If last entry less than new entry, nothing to update + if (distance > nearest_mc[offset + entry]) + return; + + // Find where to insert the new entry while sliding the rest + // to the right + entry--; + while (entry >= 0) { + if (nearest_mc[offset+entry] >= distance) + nearest_mc[offset+entry+1] = nearest_mc[offset+entry]; + else + break; + entry--; + } + + nearest_mc[offset+entry+1] = distance; +} + +} // extern "C" diff --git a/chroma/gpu/pdf.py b/chroma/gpu/pdf.py index 053d5f3..0ad78d2 100644 --- a/chroma/gpu/pdf.py +++ b/chroma/gpu/pdf.py @@ -5,7 +5,7 @@ from chroma.gpu.tools import get_cu_module, cuda_options, GPUFuncs class GPUPDF(object): def __init__(self): - self.module = get_cu_module('daq.cu', options=cuda_options, + self.module = get_cu_module('pdf.cu', options=cuda_options, include_source_directory=False) self.gpu_funcs = GPUFuncs(self.module) @@ -142,7 +142,7 @@ class GPUPDF(object): pdf_frac_uncert = np.zeros_like(pdf_value) # PDF value for high stats bins - high_stats = bincount >= self.min_bin_content + high_stats = (bincount >= self.min_bin_content) if high_stats.any(): if self.time_only: pdf_value[high_stats] = bincount[high_stats].astype(float) / hitcount[high_stats] / self.min_twidth |