1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
// -*-c++-*-
#include <curand_kernel.h>
__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"
|
