diff options
Diffstat (limited to 'src/daq.cu')
| -rw-r--r-- | src/daq.cu | 298 |
1 files changed, 144 insertions, 154 deletions
@@ -1,208 +1,198 @@ // -*-c++-*- #include <curand_kernel.h> -__device__ unsigned int float_to_sortable_int(float f) +__device__ unsigned int +float_to_sortable_int(float f) { - return __float_as_int(f); - //int i = __float_as_int(f); - //unsigned int mask = -(int)(i >> 31) | 0x80000000; - //return i ^ mask; + return __float_as_int(f); + //int i = __float_as_int(f); + //unsigned int mask = -(int)(i >> 31) | 0x80000000; + //return i ^ mask; } -__device__ float sortable_int_to_float(unsigned int i) +__device__ float +sortable_int_to_float(unsigned int i) { - return __int_as_float(i); - //unsigned int mask = ((i >> 31) - 1) | 0x80000000; - //return __int_as_float(i ^ mask); + return __int_as_float(i); + //unsigned int mask = ((i >> 31) - 1) | 0x80000000; + //return __int_as_float(i ^ mask); } extern "C" { -__global__ void reset_earliest_time_int(float maxtime, int ntime_ints, unsigned int *time_ints) +__global__ void +reset_earliest_time_int(float maxtime, int ntime_ints, unsigned int *time_ints) { - int id = threadIdx.x + blockDim.x * blockIdx.x; - if (id < ntime_ints) { - unsigned int maxtime_int = float_to_sortable_int(maxtime); - time_ints[id] = maxtime_int; - } + int id = threadIdx.x + blockDim.x * blockIdx.x; + if (id < ntime_ints) { + unsigned int maxtime_int = float_to_sortable_int(maxtime); + time_ints[id] = maxtime_int; + } } -__global__ void run_daq(curandState *s, unsigned int detection_state, - float time_rms, - 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) +__global__ void +run_daq(curandState *s, unsigned int detection_state, float time_rms, + 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 id = threadIdx.x + blockDim.x * blockIdx.x; + int id = threadIdx.x + blockDim.x * blockIdx.x; - if (id < nphotons) - { - curandState rng = s[id]; - int photon_id = id + first_photon; - int triangle_id = last_hit_triangles[photon_id]; + if (id < nphotons) { + curandState rng = s[id]; + int photon_id = id + first_photon; + int triangle_id = last_hit_triangles[photon_id]; - if (triangle_id > -1) - { - int solid_id = solid_map[triangle_id]; - unsigned int history = photon_histories[photon_id]; - - if (solid_id < nsolids && (history & detection_state)) - { - float time = photon_times[photon_id] + curand_normal(&rng) * time_rms; - 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); - } - - } + if (triangle_id > -1) { + int solid_id = solid_map[triangle_id]; + unsigned int history = photon_histories[photon_id]; - s[id] = rng; + if (solid_id < nsolids && (history & detection_state)) { + float time = photon_times[photon_id] + curand_normal(&rng) * time_rms; + 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); + } } + s[id] = rng; + + } + } -__global__ void convert_sortable_int_to_float(int n, - unsigned int *sortable_ints, - float *float_output) +__global__ void +convert_sortable_int_to_float(int n, unsigned int *sortable_ints, + float *float_output) { - int id = threadIdx.x + blockDim.x * blockIdx.x; + int id = threadIdx.x + blockDim.x * blockIdx.x; - if (id < n) - float_output[id] = sortable_int_to_float(sortable_ints[id]); + if (id < n) + 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) +__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; + int id = threadIdx.x + blockDim.x * blockIdx.x; - if (id >= nchannels) - return; + if (id >= nchannels) + return; - unsigned int q = channel_q[id]; - float t = channel_time[id]; + 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; + 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; + 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; - } + // 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) +__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; + int id = threadIdx.x + blockDim.x * blockIdx.x; - if (id >= nchannels) - return; + 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 + // 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 + // 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; + // 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 + // 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 + // 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 + 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; + // 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 + // 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; + // 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 + // 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--; - 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; + nearest_mc[offset+entry+1] = distance; } - } // extern "C" |