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
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
|
// -*-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;
}
__global__ void
accumulate_moments(int time_only, int nchannels,
float *mc_time,
unsigned int *mc_charge, // quirk of DAQ!
float tmin, float tmax,
float qmin, float qmax,
unsigned int *mom0,
float *t_mom1,
float *t_mom2,
float *q_mom1,
float *q_mom2)
{
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];
// Is this channel inside the range of the PDF?
if (time_only) {
if (channel_mc_time < tmin || channel_mc_time > tmax)
return; // Nothing else to do
mom0[id] += 1;
t_mom1[id] += channel_mc_time;
t_mom2[id] += channel_mc_time*channel_mc_time;
}
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
mom0[id] += 1;
t_mom1[id] += channel_mc_time;
t_mom2[id] += channel_mc_time*channel_mc_time;
q_mom1[id] += channel_mc_charge;
q_mom2[id] += channel_mc_charge*channel_mc_charge;
}
}
static const float invroot2 = 0.70710678118654746f; // 1/sqrt(2)
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 tmin, float tmax,
float qmin, float qmax,
float *inv_time_bandwidths,
float *inv_charge_bandwidths,
unsigned int *hitcount,
float *pdf_values)
{
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];
// Is this channel inside the range of the PDF?
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
// Kernel argument
float channel_event_time = event_time[id];
float inv_bandwidth = inv_time_bandwidths[id];
float arg = (channel_mc_time - channel_event_time) * inv_bandwidth;
// evaluate 1D Gaussian normalized within time window
float term = expf(-0.5f * arg * arg) * inv_bandwidth;
float norm = tmax - tmin;
if (inv_bandwidth > 0.0f) {
float loarg = (tmin - channel_mc_time)*inv_bandwidth*invroot2;
float hiarg = (tmax - channel_mc_time)*inv_bandwidth*invroot2;
norm = (erff(hiarg) - erff(loarg)) * rootPiBy2;
}
pdf_values[id] += term / norm;
}
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
// Kernel argument: time dim
float channel_event_obs = event_time[id];
float inv_bandwidth = inv_time_bandwidths[id];
float arg = (channel_mc_time - channel_event_obs) * inv_bandwidth;
float loarg = (tmin - channel_mc_time)*inv_bandwidth*invroot2;
float hiarg = (tmax - channel_mc_time)*inv_bandwidth*invroot2;
float norm = (erff(hiarg) - erff(loarg)) * rootPiBy2;
// Kernel argument: charge dim
channel_event_obs = event_time[id];
inv_bandwidth = inv_time_bandwidths[id];
arg *= (channel_mc_charge - channel_event_obs) * inv_bandwidth;
loarg = (tmin - channel_mc_charge)*inv_bandwidth*invroot2;
hiarg = (tmax - channel_mc_charge)*inv_bandwidth*invroot2;
norm *= (erff(hiarg) - erff(loarg)) * rootPiBy2;
// evaluate 2D Gaussian normalized within time window
float term = expf(-0.5f * arg * arg);
pdf_values[id] += term / norm;
}
}
} // extern "C"
|
