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
|
// -*-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"
|
