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
|
// -*-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;
}
}
} // extern "C"
|
