// -*-c++-*- #include __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; } __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); } extern "C" { __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; } } __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; 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); } } s[id] = rng; } } __global__ void convert_sortable_int_to_float(int n, unsigned int *sortable_ints, float *float_output) { int id = threadIdx.x + blockDim.x * blockIdx.x; 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) { 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"