ToolDAQFramework
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
library_daq.h
Go to the documentation of this file.
1 
2 #ifndef LIBRARY_DAQ_H
3 #define LIBRARY_DAQ_H
4 
5 #pragma once
6 
7 #include "build.h"
8 
9 #include <unistd.h>
10 #include <stdlib.h>
11 #include <stdio.h>
12 #include <string.h>
13 #include <thrust/extrema.h>
14 #include <limits>
15 #include <limits.h>
16 #include <thrust/sort.h>
17 #include <thrust/device_ptr.h>
18 
19 typedef unsigned short offset_t;
20 
22 // define variable types
24 #if defined __HISTOGRAM_UCHAR__
25 typedef unsigned char histogram_t;
26 #elif defined __HISTOGRAM_USHORT__
27 typedef unsigned short histogram_t;
28 #elif defined __HISTOGRAM_UINT__
29 typedef unsigned int histogram_t;
30 #endif
31 
32 #if defined __TIME_OF_FLIGHT_UCHAR__
33 typedef unsigned char time_of_flight_t;
34 #elif defined __TIME_OF_FLIGHT_USHORT__
35 typedef unsigned short time_of_flight_t;
36 #elif defined __TIME_OF_FLIGHT_UINT__
37 typedef unsigned int time_of_flight_t;
38 #endif
39 
41 // define global variables //
44 double distance_between_vertices; // linear distance between test vertices
45 double wall_like_distance; // distance from wall (in units of distance_between_vertices) to define wall-like events
46 unsigned int time_step_size; // time binning for the trigger
47 __constant__ unsigned int constant_time_step_size;
48 unsigned int water_like_threshold_number_of_pmts; // number of pmts above which a trigger is possible for water-like events
49 unsigned int wall_like_threshold_number_of_pmts; // number of pmts above which a trigger is possible for wall-like events
50 unsigned int nhits_threshold_min, nhits_threshold_max;
51 double coalesce_time; // time such that if two triggers are closer than this they are coalesced into a single trigger
52 double trigger_gate_up; // duration to be saved after the trigger time
53 double trigger_gate_down; // duration to be saved before the trigger time
54 unsigned int max_n_hits_per_job; // max n of pmt hits per job
55 double dark_rate;
56 float costheta_cone_cut; // max distance between measured costheta and cerenkov costheta
57 __constant__ float constant_costheta_cone_cut;
58 bool select_based_on_cone; // for mode == 10, set to 0 to select based on vertices, to 1 to select based on cone
59 __constant__ bool constant_select_based_on_cone;
60 bool return_vertex; // return vertex position
61 bool return_direction; // return cone direction
62 __constant__ bool constant_return_direction;
64 double detector_height; // detector height
65 double detector_radius; // detector radius
67 unsigned int n_PMTs; // number of pmts in the detector
68 __constant__ unsigned int constant_n_PMTs;
69 double * PMT_x, *PMT_y, *PMT_z; // coordinates of the pmts in the detector
71 unsigned int n_test_vertices; // number of test vertices
72 unsigned int n_water_like_test_vertices; // number of test vertices
73 __constant__ unsigned int constant_n_test_vertices;
74 __constant__ unsigned int constant_n_water_like_test_vertices;
75 double * vertex_x, * vertex_y, * vertex_z; // coordinates of test vertices
77 unsigned int number_of_kernel_blocks; // number of cores to be used
78 dim3 number_of_kernel_blocks_3d;
79 unsigned int number_of_threads_per_block; // number of threads per core to be used
80 dim3 number_of_threads_per_block_3d;
81 unsigned int grid_size; // grid = (n cores) X (n threads / core)
83 offset_t time_offset; // ns, offset to make times positive
84 __constant__ offset_t constant_time_offset;
85 unsigned int n_time_bins; // number of time bins
86 __constant__ unsigned int constant_n_time_bins;
87 unsigned int n_direction_bins_theta; // number of direction bins
88 __constant__ unsigned int constant_n_direction_bins_theta;
89 unsigned int n_direction_bins_phi; // number of direction bins
90 __constant__ unsigned int constant_n_direction_bins_phi;
91 unsigned int n_direction_bins; // number of direction bins
92 __constant__ unsigned int constant_n_direction_bins;
93 unsigned int n_hits; // number of input hits from the detector
94 __constant__ unsigned int constant_n_hits;
95 unsigned int * host_ids; // pmt id of a hit
96 unsigned int *device_ids;
97 texture<unsigned int, 1, cudaReadModeElementType> tex_ids;
98 unsigned int * host_times; // time of a hit
99 unsigned int *device_times;
100 texture<unsigned int, 1, cudaReadModeElementType> tex_times;
101 // corrected tim bin of each hit (for each vertex)
102 unsigned int * host_time_bin_of_hit;
103 unsigned int * device_time_bin_of_hit;
104 // npmts per time bin
105 histogram_t * device_n_pmts_per_time_bin; // number of active pmts in a time bin
106 unsigned int * host_n_pmts_per_time_bin;
107 unsigned int * device_n_pmts_nhits; // number of active pmts
108 unsigned int * host_n_pmts_nhits;
109 unsigned int * device_n_pmts_per_time_bin_and_direction_bin; // number of active pmts in a time bin and direction bin
110 float * device_dx_per_time_bin; // dx from vertex to hit in a time bin
111 float * device_dy_per_time_bin; // dy from vertex to hit in a time bin
112 float * device_dz_per_time_bin; // dz from vertex to hit in a time bin
113 float * device_meanx_per_time_bin; // sum of x positions of chosen PMTs
114 float * device_meany_per_time_bin; // sum of y positions of chosen PMTs
115 float * device_meanz_per_time_bin; // sum of z positions of chosen PMTs
116 //unsigned int * device_time_nhits; // trigger time
117 //unsigned int * host_time_nhits;
118 // tof
119 double speed_light_water;
120 double cerenkov_angle_water;
121 float cerenkov_costheta;
122 __constant__ float constant_cerenkov_costheta;
123 double twopi;
124 bool cylindrical_grid;
125 time_of_flight_t *device_times_of_flight; // time of flight between a vertex and a pmt
126 time_of_flight_t *host_times_of_flight;
127 float *device_light_dx; // x distance between a vertex and a pmt
128 float *host_light_dx;
129 float *device_light_dy; // y distance between a vertex and a pmt
130 float *host_light_dy;
131 float *device_light_dz; // z distance between a vertex and a pmt
132 float *host_light_dz;
133 float *device_light_dr; // distance between a vertex and a pmt
134 float *host_light_dr;
135 float *device_PMTx; // PMT x
136 float *host_PMTx;
137 float *device_PMTy; // PMT y
138 float *host_PMTy;
139 float *device_PMTz; // PMT z
140 float *host_PMTz;
141 bool *device_directions_for_vertex_and_pmt; // test directions for vertex and pmt
142 bool *host_directions_for_vertex_and_pmt;
143 texture<time_of_flight_t, 1, cudaReadModeElementType> tex_times_of_flight;
144 texture<float, 1, cudaReadModeElementType> tex_light_dx;
145 texture<float, 1, cudaReadModeElementType> tex_light_dy;
146 texture<float, 1, cudaReadModeElementType> tex_light_dz;
147 texture<float, 1, cudaReadModeElementType> tex_light_dr;
148 texture<float, 1, cudaReadModeElementType> tex_PMTx;
149 texture<float, 1, cudaReadModeElementType> tex_PMTy;
150 texture<float, 1, cudaReadModeElementType> tex_PMTz;
151 //texture<bool, 1, cudaReadModeElementType> tex_directions_for_vertex_and_pmt;
152 // triggers
153 std::vector<std::pair<unsigned int,unsigned int> > candidate_trigger_pair_vertex_time; // pair = (v, t) = (a vertex, a time at the end of the 2nd of two coalesced bins)
154 std::vector<unsigned int> candidate_trigger_npmts_in_time_bin; // npmts in time bin
155 std::vector<unsigned int> candidate_trigger_npmts_in_cone_in_time_bin;
156 std::vector<std::pair<unsigned int,unsigned int> > trigger_pair_vertex_time;
157 std::vector<unsigned int> trigger_npmts_in_time_bin;
158 std::vector<unsigned int> trigger_npmts_in_cone_in_time_bin;
159 std::vector<std::pair<unsigned int,unsigned int> > final_trigger_pair_vertex_time;
160 std::vector<double> output_trigger_information;
161 std::vector<float> candidate_trigger_meanx_per_time_bin; // sum of pmts x in time bin
162 std::vector<float> candidate_trigger_meany_per_time_bin; // sum of pmts y in time bin
163 std::vector<float> candidate_trigger_meanz_per_time_bin; // sum of pmts z in time bin
164 std::vector<float> trigger_meanx_per_time_bin;
165 std::vector<float> trigger_meany_per_time_bin;
166 std::vector<float> trigger_meanz_per_time_bin;
167 // C timing
168 struct timeval t0;
169 struct timeval t1;
170 // CUDA timing
171 cudaEvent_t start, stop, total_start, total_stop;
172 // make output txt file for plotting?
173 bool output_txt;
174 unsigned int correct_mode;
175 unsigned int write_output_mode;
176 // find candidates
177 histogram_t * host_max_number_of_pmts_in_time_bin;
178 histogram_t * device_max_number_of_pmts_in_time_bin;
179 unsigned int * host_vertex_with_max_n_pmts;
180 unsigned int * device_vertex_with_max_n_pmts;
181 unsigned int * device_number_of_pmts_in_cone_in_time_bin;
182 unsigned int * host_max_number_of_pmts_in_cone_in_time_bin;
183 unsigned int * device_max_number_of_pmts_in_cone_in_time_bin;
184 float * device_meanx_for_vertex_with_max_n_pmts_per_time_bin;
185 float * device_meany_for_vertex_with_max_n_pmts_per_time_bin;
186 float * device_meanz_for_vertex_with_max_n_pmts_per_time_bin;
187 float * host_meanx_for_vertex_with_max_n_pmts_per_time_bin;
188 float * host_meany_for_vertex_with_max_n_pmts_per_time_bin;
189 float * host_meanz_for_vertex_with_max_n_pmts_per_time_bin;
190 // gpu properties
191 int max_n_threads_per_block;
192 int max_n_blocks;
193 // verbosity level
194 bool use_verbose;
195 bool use_timing;
196 // files
197 std::string event_file;
198 std::string detector_file;
199 std::string pmts_file;
200 std::string output_file;
201 std::string event_file_base;
202 std::string event_file_suffix;
203 std::string output_file_base;
204 float elapsed_parameters, elapsed_pmts, elapsed_detector, elapsed_vertices,
205  elapsed_threads, elapsed_tof, elapsed_directions, elapsed_memory_tofs_dev, elapsed_memory_directions_dev, elapsed_memory_candidates_host, elapsed_tofs_copy_dev, elapsed_directions_copy_dev,
206  elapsed_input, elapsed_memory_dev, elapsed_copy_dev, elapsed_kernel_correct_times_and_get_n_pmts_per_time_bin,
207  elapsed_threads_candidates, elapsed_candidates_memory_dev, elapsed_candidates_kernel,
208  elapsed_candidates_copy_host, choose_candidates, elapsed_coalesce, elapsed_gates, elapsed_free, elapsed_total,
209  elapsed_tofs_free, elapsed_reset, elapsed_write_output;
210 unsigned int greatest_divisor;
211 unsigned int the_max_time;
212 unsigned int nhits_window;
213 int n_events;
214 
215 
216 __global__ void kernel_find_vertex_with_max_npmts_in_timebin(histogram_t * np, histogram_t * mnp, unsigned int * vmnp);
217 __global__ void kernel_find_vertex_with_max_npmts_in_timebin(histogram_t * np, histogram_t * mnp, unsigned int * vmnp, float * meanxin, float * meanyin, float * meanzin, float * meanxout, float * meanyout, float * meanzout);
218 __global__ void kernel_find_vertex_with_max_npmts_and_center_of_mass_in_timebin(histogram_t * np, histogram_t * mnp, unsigned int * vmnp, unsigned int *nc, unsigned int *mnc);
219 __global__ void kernel_find_vertex_with_max_npmts_in_timebin_and_directionbin(unsigned int * np, histogram_t * mnp, unsigned int * vmnp);
220 
221 unsigned int read_number_of_input_hits();
222 bool read_input();
223 void print_parameters();
224 void print_parameters_2d();
225 void print_input();
226 void print_times_of_flight();
227 void print_directions();
228 void print_pmts();
229 bool read_detector();
230 bool read_the_pmts();
231 bool read_the_detector();
232 void make_test_vertices();
233 bool setup_threads_for_tof();
234 bool setup_threads_for_tof_biparallel();
235 bool setup_threads_for_tof_2d(unsigned int A, unsigned int B);
236 bool setup_threads_to_find_candidates();
237 bool setup_threads_nhits();
238 bool read_the_input();
239 bool read_the_input_ToolDAQ(std::vector<int> PMTid, std::vector<int> time);
240 void allocate_tofs_memory_on_device();
241 void allocate_directions_memory_on_device();
242 void allocate_correct_memory_on_device();
243 void allocate_correct_memory_on_device_nhits();
244 void allocate_candidates_memory_on_host();
245 void allocate_candidates_memory_on_device();
246 void make_table_of_tofs();
247 void make_table_of_directions();
248 void fill_correct_memory_on_device();
249 void fill_tofs_memory_on_device();
250 void fill_directions_memory_on_device();
251 void fill_tofs_memory_on_device_nhits();
252 void coalesce_triggers();
253 void separate_triggers_into_gates();
254 void separate_triggers_into_gates(std::vector<int> * trigger_ns, std::vector<int> * trigger_ts, std::vector<double> * trigger_vtx_xs, std::vector<double> * trigger_vtx_ys, std::vector<double> * trigger_vtx_zs, std::vector<double> * trigger_dir_xs, std::vector<double> * trigger_dir_ys, std::vector<double> * trigger_dir_zs);
255 float timedifference_msec(struct timeval t0, struct timeval t1);
256 void start_c_clock();
257 double stop_c_clock();
258 void start_cuda_clock();
259 double stop_cuda_clock();
260 void start_total_cuda_clock();
261 double stop_total_cuda_clock();
262 unsigned int get_distance_index(unsigned int pmt_id, unsigned int vertex_block);
263 unsigned int get_time_index(unsigned int hit_index, unsigned int vertex_block);
264 unsigned int get_direction_index_at_pmt(unsigned int pmt_id, unsigned int vertex_index, unsigned int direction_index);
265 unsigned int get_direction_index_at_angles(unsigned int iphi, unsigned int itheta);
266 unsigned int get_direction_index_at_time(unsigned int time_bin, unsigned int vertex_index, unsigned int direction_index);
267 __device__ unsigned int device_get_distance_index(unsigned int pmt_id, unsigned int vertex_block);
268 __device__ unsigned int device_get_time_index(unsigned int hit_index, unsigned int vertex_block);
269 __device__ unsigned int device_get_direction_index_at_pmt(unsigned int pmt_id, unsigned int vertex_index, unsigned int direction_index);
270 __device__ unsigned int device_get_direction_index_at_angles(unsigned int iphi, unsigned int itheta);
271 __device__ unsigned int device_get_direction_index_at_time(unsigned int time_bin, unsigned int vertex_index, unsigned int direction_index);
272 void print_gpu_properties();
273 unsigned int read_number_of_pmts();
274 bool read_pmts();
275 void free_event_memories();
276 void free_event_memories_nhits();
277 void free_global_memories();
278 void copy_candidates_from_device_to_host();
279 void choose_candidates_above_threshold();
280 bool set_input_file_for_event(int n);
281 void set_output_file();
282 void set_output_file_nhits(unsigned int threshold);
283 void write_output();
284 void write_output_nhits(unsigned int n);
285 void initialize_output();
286 void initialize_output_nhits();
287 float read_value_from_file(std::string paramname, std::string filename);
288 void read_user_parameters();
289 void read_user_parameters_nhits();
290 void check_cudamalloc_float(unsigned int size);
291 void check_cudamalloc_int(unsigned int size);
292 void check_cudamalloc_unsigned_short(unsigned int size);
293 void check_cudamalloc_unsigned_char(unsigned int size);
294 void check_cudamalloc_unsigned_int(unsigned int size);
295 void check_cudamalloc_bool(unsigned int size);
296 void setup_threads_for_histo(unsigned int n);
297 unsigned int find_greatest_divisor(unsigned int n, unsigned int max);
298 void setup_threads_for_histo();
299 void setup_threads_for_histo_iterated(bool last);
300 void setup_threads_for_histo_per(unsigned int n);
301 
302 
303 unsigned int read_number_of_input_hits(){
304 
305  FILE *f=fopen(event_file.c_str(), "r");
306  if (f == NULL){
307  printf(" [2] cannot read input file \n");
308  fclose(f);
309  return 0;
310  }
311 
312  unsigned int n_hits = 0;
313 
314  for (char c = getc(f); c != EOF; c = getc(f))
315  if (c == '\n')
316  n_hits ++;
317 
318  fclose(f);
319  return n_hits;
320 
321 }
322 
323 bool read_input(){
324 
325  FILE *f=fopen(event_file.c_str(), "r");
326 
327  int time;
328  double timef;
329  unsigned int id;
330  int min = INT_MAX;
331  int max = INT_MIN;
332  for( unsigned int i=0; i<n_hits; i++){
333  if( fscanf(f, "%d %lf", &id, &timef) != 2 ){
334  printf(" [2] problem scanning hit %d from input \n", i);
335  fclose(f);
336  return false;
337  }
338  time = int(floor(timef));
339  host_times[i] = time;
340  host_ids[i] = id;
341  if( time > max ) max = time;
342  if( time < min ) min = time;
343  }
344 
345  if( min < 0 ){
346  for(int i=0; i<n_hits; i++){
347  host_times[i] -= min;
348  }
349  max -= min;
350  }
351 
352 
353  the_max_time = max;
354 
355  fclose(f);
356  return true;
357 
358 }
359 
360 
361 bool read_detector(){
362 
363  FILE *f=fopen(detector_file.c_str(), "r");
364  double pmt_radius;
365  if( fscanf(f, "%lf %lf %lf", &detector_height, &detector_radius, &pmt_radius) != 3 ){
366  printf(" [2] problem scanning detector \n");
367  fclose(f);
368  return false;
369  }
370 
371  fclose(f);
372  return true;
373 
374 }
375 
376 
377 
378 void print_parameters(){
379 
380  printf(" [2] n_test_vertices = %d \n", n_test_vertices);
381  printf(" [2] n_water_like_test_vertices = %d \n", n_water_like_test_vertices);
382  printf(" [2] n_PMTs = %d \n", n_PMTs);
383  printf(" [2] number_of_kernel_blocks = %d \n", number_of_kernel_blocks);
384  printf(" [2] number_of_threads_per_block = %d \n", number_of_threads_per_block);
385  printf(" [2] grid size = %d -> %d \n", number_of_kernel_blocks*number_of_threads_per_block, grid_size);
386 
387 }
388 
389 void print_parameters_2d(){
390 
391  printf(" [2] n_test_vertices = %d \n", n_test_vertices);
392  printf(" [2] n_water_like_test_vertices = %d \n", n_water_like_test_vertices);
393  printf(" [2] n_PMTs = %d \n", n_PMTs);
394  printf(" [2] number_of_kernel_blocks = (%d, %d) = %d \n", number_of_kernel_blocks_3d.x, number_of_kernel_blocks_3d.y, number_of_kernel_blocks_3d.x * number_of_kernel_blocks_3d.y);
395  printf(" [2] number_of_threads_per_block = (%d, %d) = %d \n", number_of_threads_per_block_3d.x, number_of_threads_per_block_3d.y, number_of_threads_per_block_3d.x * number_of_threads_per_block_3d.y);
396  printf(" [2] grid size = %d -> %d \n", number_of_kernel_blocks_3d.x*number_of_kernel_blocks_3d.y*number_of_threads_per_block_3d.x*number_of_threads_per_block_3d.y, grid_size);
397 
398 }
399 
400 void print_input(){
401 
402  for(unsigned int i=0; i<n_hits; i++)
403  printf(" [2] hit %d time %d id %d \n", i, host_times[i], host_ids[i]);
404 
405 }
406 
407 void print_pmts(){
408 
409  for(unsigned int i=0; i<n_PMTs; i++)
410  printf(" [2] pmt %d x %f y %f z %f \n", i, PMT_x[i], PMT_y[i], PMT_z[i]);
411 
412 }
413 
414 void print_times_of_flight(){
415 
416  printf(" [2] times_of_flight: (vertex, PMT) \n");
417  unsigned int distance_index;
418  for(unsigned int iv=0; iv<n_test_vertices; iv++){
419  printf(" ( ");
420  for(unsigned int ip=0; ip<n_PMTs; ip++){
421  distance_index = get_distance_index(ip + 1, n_PMTs*iv);
422  printf(" %d ", host_times_of_flight[distance_index]);
423  }
424  printf(" ) \n");
425  }
426 }
427 
428 
429 void print_directions(){
430 
431  printf(" [2] directions: (vertex, PMT) \n");
432  for(unsigned int iv=0; iv<n_test_vertices; iv++){
433  printf(" [ ");
434  for(unsigned int ip=0; ip<n_PMTs; ip++){
435  printf(" ( ");
436  for(unsigned int id=0; id<n_direction_bins; id++){
437  printf("%d ", host_directions_for_vertex_and_pmt[get_direction_index_at_pmt(ip, iv, id)]);
438  }
439  printf(" ) ");
440  }
441  printf(" ] \n");
442  }
443 }
444 
445 
446 bool read_the_pmts(){
447 
448  printf(" [2] --- read pmts \n");
449  n_PMTs = read_number_of_pmts();
450  if( !n_PMTs ) return false;
451  printf(" [2] detector contains %d PMTs \n", n_PMTs);
452  PMT_x = (double *)malloc(n_PMTs*sizeof(double));
453  PMT_y = (double *)malloc(n_PMTs*sizeof(double));
454  PMT_z = (double *)malloc(n_PMTs*sizeof(double));
455  if( !read_pmts() ) return false;
456  //print_pmts();
457  return true;
458 
459 }
460 
461 bool read_the_detector(){
462 
463  printf(" [2] --- read detector \n");
464  if( !read_detector() ) return false;
465  printf(" [2] detector height %f cm, radius %f cm \n", detector_height, detector_radius);
466  return true;
467 
468 }
469 
470 void make_test_vertices(){
471 
472  printf(" [2] --- make test vertices \n");
473  float semiheight = detector_height/2.;
474  n_test_vertices = 0;
475 
476 
477  if( !cylindrical_grid ){
478 
479  // 1: count number of test vertices
480  for(int i=-1*semiheight; i <= semiheight; i+=distance_between_vertices) {
481  for(int j=-1*detector_radius; j<=detector_radius; j+=distance_between_vertices) {
482  for(int k=-1*detector_radius; k<=detector_radius; k+=distance_between_vertices) {
483  if(pow(j,2)+pow(k,2) > pow(detector_radius,2))
484  continue;
485  n_test_vertices++;
486  }
487  }
488  }
489  vertex_x = (double *)malloc(n_test_vertices*sizeof(double));
490  vertex_y = (double *)malloc(n_test_vertices*sizeof(double));
491  vertex_z = (double *)malloc(n_test_vertices*sizeof(double));
492 
493  // 2: assign coordinates to test vertices
494  // water-like events
495  n_test_vertices = 0;
496  for(int i=-1*semiheight; i <= semiheight; i+=distance_between_vertices) {
497  for(int j=-1*detector_radius; j<=detector_radius; j+=distance_between_vertices) {
498  for(int k=-1*detector_radius; k<=detector_radius; k+=distance_between_vertices) {
499 
500 
501  if(
502  // skip endcap region
503  abs(i) > semiheight - wall_like_distance*distance_between_vertices ||
504  // skip sidewall region
505  pow(j,2)+pow(k,2) > pow(detector_radius - wall_like_distance*distance_between_vertices,2)
506  ) continue;
507 
508  vertex_x[n_test_vertices] = j*1.;
509  vertex_y[n_test_vertices] = k*1.;
510  vertex_z[n_test_vertices] = i*1.;
511  n_test_vertices++;
512  }
513  }
514  }
515  n_water_like_test_vertices = n_test_vertices;
516 
517  // wall-like events
518  for(int i=-1*semiheight; i <= semiheight; i+=distance_between_vertices) {
519  for(int j=-1*detector_radius; j<=detector_radius; j+=distance_between_vertices) {
520  for(int k=-1*detector_radius; k<=detector_radius; k+=distance_between_vertices) {
521 
522  if(
523  abs(i) > semiheight - wall_like_distance*distance_between_vertices ||
524  pow(j,2)+pow(k,2) > pow(detector_radius - wall_like_distance*distance_between_vertices,2)
525  ){
526 
527  if(pow(j,2)+pow(k,2) > pow(detector_radius,2)) continue;
528 
529  vertex_x[n_test_vertices] = j*1.;
530  vertex_y[n_test_vertices] = k*1.;
531  vertex_z[n_test_vertices] = i*1.;
532  n_test_vertices++;
533  }
534  }
535  }
536  }
537 
538  }else{ // cylindrical grid
539 
540  int n_vertical = detector_height/distance_between_vertices;
541  double distance_vertical = detector_height/n_vertical;
542  int n_radial = 2.*detector_radius/distance_between_vertices;
543  double distance_radial = 2.*detector_radius/n_radial;
544  int n_angular;
545  double distance_angular;
546 
547  printf(" [2] distance_between_vertices %f, distance_vertical %f, distance_radial %f \n",
548  distance_between_vertices, distance_vertical, distance_radial);
549 
550  double the_r, the_z, the_phi;
551  bool first = false; // true: add extra layer near wall
552  // false: regular spacing
553 
554  bool add_extra_layer = first;
555 
556  // 1: count number of test vertices
557  the_r = detector_radius;
558  while( the_r >= 0. ){
559  n_angular = twopi*the_r / distance_between_vertices;
560  distance_angular = twopi/n_angular;
561 
562  the_z = -semiheight;
563 
564  while( the_z <= semiheight){
565 
566  the_phi = 0.;
567  while( the_phi < twopi - distance_angular/2. ){
568 
569  n_test_vertices ++;
570 
571  if( the_r == 0. ) break;
572  the_phi += distance_angular;
573  }
574 
575 
576  if( add_extra_layer ){
577  if( the_z + semiheight < 0.3*distance_vertical ) // only true at bottom endcap
578  the_z += distance_vertical/2.;
579  else if( semiheight - the_z < 0.7*distance_vertical ) // only true near top endcap
580  the_z += distance_vertical/2.;
581  else
582  the_z += distance_vertical;
583  }else{
584  the_z += distance_vertical;
585  }
586  }
587  if( first ){
588  the_r -= distance_radial/2.;
589  first = false;
590  }
591  else{
592  the_r -= distance_radial;
593  }
594  }
595 
596  vertex_x = (double *)malloc(n_test_vertices*sizeof(double));
597  vertex_y = (double *)malloc(n_test_vertices*sizeof(double));
598  vertex_z = (double *)malloc(n_test_vertices*sizeof(double));
599 
600  first = add_extra_layer;
601  // 2: assign coordinates to test vertices
602  // water-like events
603  n_test_vertices = 0;
604 
605  the_r = detector_radius;
606  while( the_r >= 0. ){
607 
608  // skip sidewall region
609  if(the_r <= detector_radius - wall_like_distance*distance_between_vertices ){
610 
611  n_angular = twopi*the_r / distance_between_vertices;
612  distance_angular = twopi/n_angular;
613 
614  the_z = -semiheight;
615 
616  while( the_z <= semiheight){
617 
618  // skip endcap region
619  if( fabs(the_z) <= semiheight - wall_like_distance*distance_between_vertices ){
620 
621  the_phi = 0.;
622  while( the_phi < twopi - distance_angular/2. ){
623 
624  vertex_x[n_test_vertices] = the_r*cos(the_phi);
625  vertex_y[n_test_vertices] = the_r*sin(the_phi);
626  vertex_z[n_test_vertices] = the_z;
627  n_test_vertices ++;
628 
629  if( the_r == 0. ) break;
630  the_phi += distance_angular;
631  }
632  }
633 
634  if( add_extra_layer ){
635  if( the_z + semiheight < 0.3*distance_vertical ) // only true at bottom endcap
636  the_z += distance_vertical/2.;
637  else if( semiheight - the_z < 0.7*distance_vertical ) // only true near top endcap
638  the_z += distance_vertical/2.;
639  else
640  the_z += distance_vertical;
641  }else{
642  the_z += distance_vertical;
643  }
644  }
645 
646  }
647  if( first ){
648  the_r -= distance_radial/2.;
649  first = false;
650  }
651  else{
652  the_r -= distance_radial;
653  }
654  }
655 
656 
657  n_water_like_test_vertices = n_test_vertices;
658 
659  first = add_extra_layer;
660  // wall-like events
661  the_r = detector_radius;
662  while( the_r >= 0. ){
663 
664  n_angular = twopi*the_r / distance_between_vertices;
665  distance_angular = twopi/n_angular;
666 
667  the_z = -semiheight;
668 
669  while( the_z <= semiheight){
670 
671  if( fabs(the_z) > semiheight - wall_like_distance*distance_between_vertices ||
672  the_r > detector_radius - wall_like_distance*distance_between_vertices ){
673 
674  the_phi = 0.;
675  while( the_phi < twopi - distance_angular/2. ){
676 
677  vertex_x[n_test_vertices] = the_r*cos(the_phi);
678  vertex_y[n_test_vertices] = the_r*sin(the_phi);
679  vertex_z[n_test_vertices] = the_z;
680  n_test_vertices ++;
681 
682  if( the_r == 0. ) break;
683  the_phi += distance_angular;
684  }
685  }
686  if( add_extra_layer ){
687  if( the_z + semiheight < 0.3*distance_vertical ) // only true at bottom endcap
688  the_z += distance_vertical/2.;
689  else if( semiheight - the_z < 0.7*distance_vertical ) // only true near top endcap
690  the_z += distance_vertical/2.;
691  else
692  the_z += distance_vertical;
693  }else{
694  the_z += distance_vertical;
695  }
696  }
697  if( first ){
698  the_r -= distance_radial/2.;
699  first = false;
700  }
701  else{
702  the_r -= distance_radial;
703  }
704  }
705 
706 
707  }
708 
709  printf(" [2] made %d test vertices \n", n_test_vertices);
710 
711  return;
712 
713 }
714 
715 bool setup_threads_for_tof(){
716 
717  grid_size = n_test_vertices;
718 
719  number_of_kernel_blocks = grid_size / max_n_threads_per_block + 1;
720  number_of_threads_per_block = ( number_of_kernel_blocks > 1 ? max_n_threads_per_block : grid_size);
721 
722  print_parameters();
723 
724  if( number_of_threads_per_block > max_n_threads_per_block ){
725  printf(" [2] warning: number_of_threads_per_block = %d cannot exceed max value %d \n", number_of_threads_per_block, max_n_threads_per_block );
726  return false;
727  }
728 
729  if( number_of_kernel_blocks > max_n_blocks ){
730  printf(" [2] warning: number_of_kernel_blocks = %d cannot exceed max value %d \n", number_of_kernel_blocks, max_n_blocks );
731  return false;
732  }
733 
734  return true;
735 }
736 
737 
738 bool setup_threads_for_tof_biparallel(){
739 
740  grid_size = n_test_vertices * n_hits;
741 
742  number_of_kernel_blocks = grid_size / max_n_threads_per_block + 1;
743  number_of_threads_per_block = ( number_of_kernel_blocks > 1 ? max_n_threads_per_block : grid_size);
744 
745  print_parameters();
746 
747  if( number_of_threads_per_block > max_n_threads_per_block ){
748  printf(" [2] --------------------- warning: number_of_threads_per_block = %d cannot exceed max value %d \n", number_of_threads_per_block, max_n_threads_per_block );
749  return false;
750  }
751 
752  if( number_of_kernel_blocks > max_n_blocks ){
753  printf(" [2] warning: number_of_kernel_blocks = %d cannot exceed max value %d \n", number_of_kernel_blocks, max_n_blocks );
754  return false;
755  }
756 
757  return true;
758 
759 }
760 
761 bool setup_threads_for_tof_2d(unsigned int A, unsigned int B){
762 
763  if( std::numeric_limits<unsigned int>::max() / B < A ){
764  printf(" [2] --------------------- warning: B = %d times A = %d cannot exceed max value %u \n", B, A, std::numeric_limits<unsigned int>::max() );
765  return false;
766  }
767 
768  grid_size = A * B;
769  unsigned int max_n_threads_per_block_2d = sqrt(max_n_threads_per_block);
770 
771  number_of_kernel_blocks_3d.x = A / max_n_threads_per_block_2d + 1;
772  number_of_kernel_blocks_3d.y = B / max_n_threads_per_block_2d + 1;
773 
774  number_of_threads_per_block_3d.x = ( number_of_kernel_blocks_3d.x > 1 ? max_n_threads_per_block_2d : A);
775  number_of_threads_per_block_3d.y = ( number_of_kernel_blocks_3d.y > 1 ? max_n_threads_per_block_2d : B);
776 
777  print_parameters_2d();
778 
779  if( number_of_threads_per_block_3d.x > max_n_threads_per_block_2d ){
780  printf(" [2] --------------------- warning: number_of_threads_per_block x = %d cannot exceed max value %d \n", number_of_threads_per_block_3d.x, max_n_threads_per_block_2d );
781  return false;
782  }
783 
784  if( number_of_threads_per_block_3d.y > max_n_threads_per_block_2d ){
785  printf(" [2] --------------------- warning: number_of_threads_per_block y = %d cannot exceed max value %d \n", number_of_threads_per_block_3d.y, max_n_threads_per_block_2d );
786  return false;
787  }
788 
789  if( number_of_kernel_blocks_3d.x > max_n_blocks ){
790  printf(" [2] warning: number_of_kernel_blocks x = %d cannot exceed max value %d \n", number_of_kernel_blocks_3d.x, max_n_blocks );
791  return false;
792  }
793 
794  if( number_of_kernel_blocks_3d.y > max_n_blocks ){
795  printf(" [2] warning: number_of_kernel_blocks y = %d cannot exceed max value %d \n", number_of_kernel_blocks_3d.y, max_n_blocks );
796  return false;
797  }
798 
799  if( std::numeric_limits<int>::max() / (number_of_kernel_blocks_3d.x*number_of_kernel_blocks_3d.y) < number_of_threads_per_block_3d.x*number_of_threads_per_block_3d.y ){
800  printf(" [2] --------------------- warning: grid size cannot exceed max value %u \n", std::numeric_limits<int>::max() );
801  return false;
802  }
803 
804 
805  return true;
806 
807 }
808 
809 bool setup_threads_to_find_candidates(){
810 
811  number_of_kernel_blocks = n_time_bins / max_n_threads_per_block + 1;
812  number_of_threads_per_block = ( number_of_kernel_blocks > 1 ? max_n_threads_per_block : n_time_bins);
813 
814  if( number_of_threads_per_block > max_n_threads_per_block ){
815  printf(" [2] warning: number_of_threads_per_block = %d cannot exceed max value %d \n", number_of_threads_per_block, max_n_threads_per_block );
816  return false;
817  }
818 
819  return true;
820 }
821 
822 bool setup_threads_nhits(){
823 
824  number_of_kernel_blocks_3d.x = 100;
825  number_of_kernel_blocks_3d.y = 1;
826 
827  number_of_threads_per_block_3d.x = 1024;
828  number_of_threads_per_block_3d.y = 1;
829 
830  print_parameters_2d();
831 
832  return true;
833 
834 }
835 
836 
837 
838 bool read_the_input(){
839 
840  printf(" [2] --- read input \n");
841  n_hits = read_number_of_input_hits();
842  if( !n_hits ) return false;
843  printf(" [2] input contains %d hits \n", n_hits);
844  host_ids = (unsigned int *)malloc(n_hits*sizeof(unsigned int));
845  host_times = (unsigned int *)malloc(n_hits*sizeof(unsigned int));
846  if( !read_input() ) return false;
847  //time_offset = 600.; // set to constant to match trevor running
848  n_time_bins = int(floor((the_max_time + time_offset)/time_step_size))+1; // floor returns the integer below
849  printf(" [2] input max_time %d, n_time_bins %d \n", the_max_time, n_time_bins);
850  printf(" [2] time_offset = %d ns \n", time_offset);
851  //print_input();
852 
853  checkCudaErrors( cudaMemcpyToSymbol(constant_n_time_bins, &n_time_bins, sizeof(n_time_bins)) );
854  checkCudaErrors( cudaMemcpyToSymbol(constant_n_hits, &n_hits, sizeof(n_hits)) );
855 
856  return true;
857 }
858 
859 bool read_the_input_ToolDAQ(std::vector<int> PMTids, std::vector<int> times, int * earliest_time){
860 
861  printf(" [2] --- read input \n");
862  n_hits = PMTids.size();
863  if( !n_hits ) return false;
864  if( n_hits != times.size() ){
865  printf(" [2] n PMT ids %d but n times %d \n", n_hits, times.size());
866  return false;
867  }
868  printf(" [2] input contains %d hits \n", n_hits);
869  host_ids = (unsigned int *)malloc(n_hits*sizeof(unsigned int));
870  host_times = (unsigned int *)malloc(n_hits*sizeof(unsigned int));
871 
872  // if( !read_input() ) return false;
873  // read_input()
874  {
875  int min = INT_MAX;
876  int max = INT_MIN;
877  int time;
878  for(int i=0; i<PMTids.size(); i++){
879  time = int(floor(times[i]));
880  host_times[i] = time;
881  host_ids[i] = PMTids[i];
882  // printf(" [2] input %d PMT %d time %d \n", i, host_ids[i], host_times[i]);
883  if( time > max ) max = time;
884  if( time < min ) min = time;
885  }
886  if( min < 0 ){
887  for(int i=0; i<PMTids.size(); i++){
888  host_times[i] -= min;
889  }
890  max -= min;
891  min -= min;
892  }
893  the_max_time = max;
894  *earliest_time = min - min % time_step_size;
895  }
896 
897 
898  //time_offset = 600.; // set to constant to match trevor running
899  n_time_bins = int(floor((the_max_time + time_offset)/time_step_size))+1; // floor returns the integer below
900  printf(" [2] input max_time %d, n_time_bins %d \n", the_max_time, n_time_bins);
901  printf(" [2] time_offset = %d ns \n", time_offset);
902  //print_input();
903 
904  checkCudaErrors( cudaMemcpyToSymbol(constant_n_time_bins, &n_time_bins, sizeof(n_time_bins)) );
905  checkCudaErrors( cudaMemcpyToSymbol(constant_n_hits, &n_hits, sizeof(n_hits)) );
906 
907  return true;
908 }
909 
910 void allocate_tofs_memory_on_device(){
911 
912  printf(" [2] --- allocate memory tofs \n");
913 #if defined __TIME_OF_FLIGHT_UCHAR__
914  check_cudamalloc_unsigned_char(n_test_vertices*n_PMTs);
915 #elif defined __TIME_OF_FLIGHT_USHORT__
916  check_cudamalloc_unsigned_short(n_test_vertices*n_PMTs);
917 #elif defined __TIME_OF_FLIGHT_UINT__
918  check_cudamalloc_unsigned_int(n_test_vertices*n_PMTs);
919 #endif
920  checkCudaErrors(cudaMalloc((void **)&device_times_of_flight, n_test_vertices*n_PMTs*sizeof(time_of_flight_t)));
921 
922  if( correct_mode == 8 && return_direction ){
923  check_cudamalloc_float(3*n_PMTs);
924  checkCudaErrors(cudaMalloc((void **)&device_PMTx, n_PMTs*sizeof(float)));
925  checkCudaErrors(cudaMalloc((void **)&device_PMTy, n_PMTs*sizeof(float)));
926  checkCudaErrors(cudaMalloc((void **)&device_PMTz, n_PMTs*sizeof(float)));
927  }
928  if( correct_mode == 10 ){
929 
930  check_cudamalloc_float(3*n_test_vertices*n_PMTs);
931  checkCudaErrors(cudaMalloc((void **)&device_light_dx, n_test_vertices*n_PMTs*sizeof(float)));
932  checkCudaErrors(cudaMalloc((void **)&device_light_dy, n_test_vertices*n_PMTs*sizeof(float)));
933  checkCudaErrors(cudaMalloc((void **)&device_light_dz, n_test_vertices*n_PMTs*sizeof(float)));
934  checkCudaErrors(cudaMalloc((void **)&device_light_dr, n_test_vertices*n_PMTs*sizeof(float)));
935 
936  }
937  /*
938  if( n_hits*n_test_vertices > available_memory ){
939  printf(" [2] cannot allocate vector of %d, available_memory %d \n", n_hits*n_test_vertices, available_memory);
940  return 0;
941  }
942  */
943 
944 
945  if( correct_mode == 1 ){
946 
947  unsigned int max = max_n_threads_per_block;
948  greatest_divisor = find_greatest_divisor ( n_test_vertices , max);
949  printf(" [2] greatest divisor of %d below %d is %d \n", n_test_vertices, max, greatest_divisor);
950 
951  }
952 
953  return;
954 
955 }
956 
957 void allocate_directions_memory_on_device(){
958 
959  printf(" [2] --- allocate memory directions \n");
960  check_cudamalloc_bool( n_test_vertices*n_direction_bins*n_PMTs);
961  checkCudaErrors(cudaMalloc((void **)&device_directions_for_vertex_and_pmt, n_test_vertices*n_direction_bins*n_PMTs*sizeof(bool)));
962  /*
963  if( n_hits*n_test_vertices > available_memory ){
964  printf(" [2] cannot allocate vector of %d, available_memory %d \n", n_hits*n_test_vertices, available_memory);
965  return 0;
966  }
967  */
968 
969  return;
970 
971 }
972 
973 void allocate_correct_memory_on_device(){
974 
975  printf(" [2] --- allocate memory \n");
976  /*
977  if( n_hits > available_memory ){
978  printf(" [2] cannot allocate vector of %d, available_memory %d \n", n_hits, available_memory);
979  return 0;
980  }
981  */
982  check_cudamalloc_unsigned_int(n_hits);
983  checkCudaErrors(cudaMalloc((void **)&device_ids, n_hits*sizeof(unsigned int)));
984 
985  check_cudamalloc_unsigned_int(n_hits);
986  checkCudaErrors(cudaMalloc((void **)&device_times, n_hits*sizeof(unsigned int)));
987  /*
988  if( n_test_vertices*n_PMTs > available_memory ){
989  printf(" [2] cannot allocate vector of %d, available_memory %d \n", n_test_vertices*n_PMTs, available_memory);
990  return 0;
991  }
992  */
993 
994  if( correct_mode != 9 ){
995  check_cudamalloc_unsigned_int(n_time_bins*n_test_vertices);
996  checkCudaErrors(cudaMalloc((void **)&device_n_pmts_per_time_bin, n_time_bins*n_test_vertices*sizeof(unsigned int)));
997  if( correct_mode == 8 && return_direction ){
998  checkCudaErrors(cudaMalloc((void **)&device_meanx_per_time_bin, n_time_bins*n_test_vertices*sizeof(float)));
999  checkCudaErrors(cudaMalloc((void **)&device_meany_per_time_bin, n_time_bins*n_test_vertices*sizeof(float)));
1000  checkCudaErrors(cudaMalloc((void **)&device_meanz_per_time_bin, n_time_bins*n_test_vertices*sizeof(float)));
1001  }
1002  if( correct_mode == 10 ){
1003  check_cudamalloc_float(3*n_time_bins*n_test_vertices);
1004  checkCudaErrors(cudaMalloc((void **)&device_dx_per_time_bin, n_time_bins*n_test_vertices*sizeof(float)));
1005  checkCudaErrors(cudaMalloc((void **)&device_dy_per_time_bin, n_time_bins*n_test_vertices*sizeof(float)));
1006  checkCudaErrors(cudaMalloc((void **)&device_dz_per_time_bin, n_time_bins*n_test_vertices*sizeof(float)));
1007  check_cudamalloc_unsigned_int(n_time_bins*n_test_vertices);
1008  checkCudaErrors(cudaMalloc((void **)&device_number_of_pmts_in_cone_in_time_bin, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1009  }
1010  }
1011 
1012  if( correct_mode == 0 ){
1013  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1014  } else if( correct_mode == 1 ){
1015  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1016 
1017  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1018  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1019  // checkCudaErrors(cudaMemset(device_time_bin_of_hit, 0, n_hits*n_test_vertices*sizeof(unsigned int)));
1020  }
1021  else if( correct_mode == 2 ){
1022  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1023  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1024  //checkCudaErrors(cudaMemset(device_time_bin_of_hit, 0, n_hits*n_test_vertices*sizeof(unsigned int)));
1025 
1026  host_time_bin_of_hit = (unsigned int *)malloc(n_hits*n_test_vertices*sizeof(unsigned int));
1027 
1028  host_n_pmts_per_time_bin = (unsigned int *)malloc(n_time_bins*n_test_vertices*sizeof(unsigned int));
1029  memset(host_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int));
1030  } else if( correct_mode == 3 ){
1031  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1032 
1033  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1034  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1035  //checkCudaErrors(cudaMemset(device_time_bin_of_hit, 0, n_hits*n_test_vertices*sizeof(unsigned int)));
1036  } else if( correct_mode == 4 ){
1037  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1038 
1039  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1040  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1041  } else if( correct_mode == 5 ){
1042  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1043 
1044  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1045  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1046  } else if( correct_mode == 6 ){
1047  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1048 
1049  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1050  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1051  } else if( correct_mode == 7 ){
1052  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1053 
1054  check_cudamalloc_unsigned_int(n_hits*n_test_vertices);
1055  checkCudaErrors(cudaMalloc((void **)&device_time_bin_of_hit, n_hits*n_test_vertices*sizeof(unsigned int)));
1056  } else if( correct_mode == 8 ){
1057  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1058  if( return_direction ){
1059  checkCudaErrors(cudaMemset(device_meanx_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(float)));
1060  checkCudaErrors(cudaMemset(device_meany_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(float)));
1061  checkCudaErrors(cudaMemset(device_meanz_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(float)));
1062  }
1063  } else if( correct_mode == 9 ){
1064 
1065  check_cudamalloc_unsigned_int(n_time_bins*n_direction_bins*n_test_vertices);
1066  checkCudaErrors(cudaMalloc((void **)&device_n_pmts_per_time_bin_and_direction_bin, n_time_bins*n_direction_bins*n_test_vertices*sizeof(unsigned int)));
1067  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin_and_direction_bin, 0, n_time_bins*n_direction_bins*n_test_vertices*sizeof(unsigned int)));
1068  } else if( correct_mode == 10 ){
1069  checkCudaErrors(cudaMemset(device_n_pmts_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(unsigned int)));
1070  checkCudaErrors(cudaMemset(device_dx_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(float)));
1071  checkCudaErrors(cudaMemset(device_dy_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(float)));
1072  checkCudaErrors(cudaMemset(device_dz_per_time_bin, 0, n_time_bins*n_test_vertices*sizeof(float)));
1073  }
1074 
1075 
1076  return;
1077 
1078 }
1079 
1080 
1081 void allocate_correct_memory_on_device_nhits(){
1082 
1083  printf(" [2] --- allocate memory \n");
1084  /*
1085  if( n_hits > available_memory ){
1086  printf(" [2] cannot allocate vector of %d, available_memory %d \n", n_hits, available_memory);
1087  return 0;
1088  }
1089  */
1090  check_cudamalloc_unsigned_int(n_hits);
1091  checkCudaErrors(cudaMalloc((void **)&device_ids, n_hits*sizeof(unsigned int)));
1092 
1093  check_cudamalloc_unsigned_int(n_hits);
1094  checkCudaErrors(cudaMalloc((void **)&device_times, n_hits*sizeof(unsigned int)));
1095  /*
1096  if( n_test_vertices*n_PMTs > available_memory ){
1097  printf(" [2] cannot allocate vector of %d, available_memory %d \n", n_test_vertices*n_PMTs, available_memory);
1098  return 0;
1099  }
1100  */
1101 
1102  check_cudamalloc_unsigned_int(1);
1103  checkCudaErrors(cudaMalloc((void **)&device_n_pmts_nhits, 1*sizeof(unsigned int)));
1104  // checkCudaErrors(cudaMalloc((void **)&device_time_nhits, (nhits_window/time_step_size + 1)*sizeof(unsigned int)));
1105 
1106  host_n_pmts_nhits = (unsigned int *)malloc(1*sizeof(unsigned int));
1107  // host_time_nhits = (unsigned int *)malloc((nhits_window/time_step_size + 1)*sizeof(unsigned int));
1108 
1109  return;
1110 
1111 }
1112 
1113 void allocate_candidates_memory_on_host(){
1114 
1115  printf(" [2] --- allocate candidates memory on host \n");
1116 
1117  host_max_number_of_pmts_in_time_bin = (histogram_t *)malloc(n_time_bins*sizeof(histogram_t));
1118  host_vertex_with_max_n_pmts = (unsigned int *)malloc(n_time_bins*sizeof(unsigned int));
1119 
1120  if( correct_mode == 8 && return_direction ){
1121  host_meanx_for_vertex_with_max_n_pmts_per_time_bin = (float *)malloc(n_time_bins*sizeof(float));
1122  host_meany_for_vertex_with_max_n_pmts_per_time_bin = (float *)malloc(n_time_bins*sizeof(float));
1123  host_meanz_for_vertex_with_max_n_pmts_per_time_bin = (float *)malloc(n_time_bins*sizeof(float));
1124  }
1125 
1126  if( correct_mode == 10 ){
1127  host_max_number_of_pmts_in_cone_in_time_bin = (unsigned int *)malloc(n_time_bins*sizeof(unsigned int));
1128  }
1129 
1130  return;
1131 
1132 }
1133 
1134 void allocate_candidates_memory_on_device(){
1135 
1136  printf(" [2] --- allocate candidates memory on device \n");
1137 
1138 #if defined __HISTOGRAM_UCHAR__
1139  check_cudamalloc_unsigned_char(n_time_bins);
1140 #elif defined __HISTOGRAM_USHORT__
1141  check_cudamalloc_unsigned_short(n_time_bins);
1142 #elif defined __HISTOGRAM_UINT__
1143  check_cudamalloc_unsigned_int(n_time_bins);
1144 #endif
1145  checkCudaErrors(cudaMalloc((void **)&device_max_number_of_pmts_in_time_bin, n_time_bins*sizeof(histogram_t)));
1146 
1147  if( correct_mode == 8 && return_direction ){
1148  checkCudaErrors(cudaMalloc((void **)&device_meanx_for_vertex_with_max_n_pmts_per_time_bin, n_time_bins*sizeof(float)));
1149  checkCudaErrors(cudaMalloc((void **)&device_meany_for_vertex_with_max_n_pmts_per_time_bin, n_time_bins*sizeof(float)));
1150  checkCudaErrors(cudaMalloc((void **)&device_meanz_for_vertex_with_max_n_pmts_per_time_bin, n_time_bins*sizeof(float)));
1151  }
1152 
1153  check_cudamalloc_unsigned_int(n_time_bins);
1154  checkCudaErrors(cudaMalloc((void **)&device_vertex_with_max_n_pmts, n_time_bins*sizeof(unsigned int)));
1155 
1156  if( correct_mode == 10 ){
1157  check_cudamalloc_unsigned_int(n_time_bins);
1158  checkCudaErrors(cudaMalloc((void **)&device_max_number_of_pmts_in_cone_in_time_bin, n_time_bins*sizeof(unsigned int)));
1159  }
1160 
1161  return;
1162 
1163 }
1164 
1165 void make_table_of_tofs(){
1166 
1167  printf(" [2] --- fill times_of_flight \n");
1168  host_times_of_flight = (time_of_flight_t*)malloc(n_test_vertices*n_PMTs * sizeof(time_of_flight_t));
1169  printf(" [2] speed_light_water %f \n", speed_light_water);
1170  if( correct_mode == 8 && return_direction ){
1171  host_PMTx = (float*)malloc(n_PMTs * sizeof(double));
1172  host_PMTy = (float*)malloc(n_PMTs * sizeof(double));
1173  host_PMTz = (float*)malloc(n_PMTs * sizeof(double));
1174  }
1175  if( correct_mode == 10 ){
1176  host_light_dx = (float*)malloc(n_test_vertices*n_PMTs * sizeof(double));
1177  host_light_dy = (float*)malloc(n_test_vertices*n_PMTs * sizeof(double));
1178  host_light_dz = (float*)malloc(n_test_vertices*n_PMTs * sizeof(double));
1179  host_light_dr = (float*)malloc(n_test_vertices*n_PMTs * sizeof(double));
1180  }
1181  unsigned int distance_index;
1182  time_offset = 0.;
1183  for(unsigned int ip=0; ip<n_PMTs; ip++){
1184  if( correct_mode == 8 && return_direction ){
1185  host_PMTx[ip] = PMT_x[ip];
1186  host_PMTy[ip] = PMT_y[ip];
1187  host_PMTz[ip] = PMT_z[ip];
1188  }
1189  for(unsigned int iv=0; iv<n_test_vertices; iv++){
1190  distance_index = get_distance_index(ip + 1, n_PMTs*iv);
1191  host_times_of_flight[distance_index] = sqrt(pow(vertex_x[iv] - PMT_x[ip],2) + pow(vertex_y[iv] - PMT_y[ip],2) + pow(vertex_z[iv] - PMT_z[ip],2))/speed_light_water;
1192  if( correct_mode == 10 ){
1193  host_light_dx[distance_index] = PMT_x[ip] - vertex_x[iv];
1194  host_light_dy[distance_index] = PMT_y[ip] - vertex_y[iv];
1195  host_light_dz[distance_index] = PMT_z[ip] - vertex_z[iv];
1196  host_light_dr[distance_index] = sqrt(pow(host_light_dx[distance_index],2) + pow(host_light_dy[distance_index],2) + pow(host_light_dz[distance_index],2));
1197  }
1198  if( host_times_of_flight[distance_index] > time_offset )
1199  time_offset = host_times_of_flight[distance_index];
1200 
1201  }
1202  }
1203  //print_times_of_flight();
1204 
1205  return;
1206 }
1207 
1208 
1209 void make_table_of_directions(){
1210 
1211  printf(" [2] --- fill directions \n");
1212  printf(" [2] cerenkov_angle_water %f \n", cerenkov_angle_water);
1213  host_directions_for_vertex_and_pmt = (bool*)malloc(n_test_vertices*n_PMTs*n_direction_bins * sizeof(bool));
1214  float dx, dy, dz, dr, phi, cos_theta, sin_theta;
1215  float phi2, cos_theta2, angle;
1216  unsigned int dir_index_at_angles;
1217  unsigned int dir_index_at_pmt;
1218  for(unsigned int ip=0; ip<n_PMTs; ip++){
1219  for(unsigned int iv=0; iv<n_test_vertices; iv++){
1220  dx = PMT_x[ip] - vertex_x[iv];
1221  dy = PMT_y[ip] - vertex_y[iv];
1222  dz = PMT_z[ip] - vertex_z[iv];
1223  dr = sqrt(pow(dx,2) + pow(dy,2) + pow(dz,2));
1224  phi = atan2(dy,dx);
1225  // light direction
1226  cos_theta = dz/dr;
1227  sin_theta = sqrt(1. - pow(cos_theta,2));
1228  // particle direction
1229  for(unsigned int itheta = 0; itheta < n_direction_bins_theta; itheta++){
1230  cos_theta2 = -1. + 2.*itheta/(n_direction_bins_theta - 1);
1231  for(unsigned int iphi = 0; iphi < n_direction_bins_phi; iphi++){
1232  phi2 = 0. + twopi*iphi/n_direction_bins_phi;
1233 
1234  if( (itheta == 0 || itheta + 1 == n_direction_bins_theta ) && iphi != 0 ) break;
1235 
1236  // angle between light direction and particle direction
1237  angle = acos( sin_theta*sqrt(1 - pow(cos_theta2,2)) * cos(phi - phi2) + cos_theta*cos_theta2 );
1238 
1239  dir_index_at_angles = get_direction_index_at_angles(iphi, itheta);
1240  dir_index_at_pmt = get_direction_index_at_pmt(ip, iv, dir_index_at_angles);
1241 
1242  //printf(" [2] phi %f ctheta %f phi' %f ctheta' %f angle %f dir_index_at_angles %d dir_index_at_pmt %d \n",
1243  // phi, cos_theta, phi2, cos_theta2, angle, dir_index_at_angles, dir_index_at_pmt);
1244 
1245  host_directions_for_vertex_and_pmt[dir_index_at_pmt]
1246  = (bool)(fabs(angle - cerenkov_angle_water) < twopi/(2.*n_direction_bins_phi));
1247  }
1248  }
1249  }
1250  }
1251  //print_directions();
1252 
1253  return;
1254 }
1255 
1256 
1257 void fill_tofs_memory_on_device(){
1258 
1259  printf(" [2] --- copy tofs from host to device \n");
1260  checkCudaErrors(cudaMemcpy(device_times_of_flight,
1261  host_times_of_flight,
1262  n_test_vertices*n_PMTs*sizeof(time_of_flight_t),
1263  cudaMemcpyHostToDevice));
1264  if( correct_mode == 8 && return_direction ){
1265  checkCudaErrors(cudaMemcpy(device_PMTx,host_PMTx,n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1266  checkCudaErrors(cudaMemcpy(device_PMTy,host_PMTy,n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1267  checkCudaErrors(cudaMemcpy(device_PMTz,host_PMTz,n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1268  }
1269  if( correct_mode == 10 ){
1270  checkCudaErrors(cudaMemcpy(device_light_dx,host_light_dx,n_test_vertices*n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1271  checkCudaErrors(cudaMemcpy(device_light_dy,host_light_dy,n_test_vertices*n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1272  checkCudaErrors(cudaMemcpy(device_light_dz,host_light_dz,n_test_vertices*n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1273  checkCudaErrors(cudaMemcpy(device_light_dr,host_light_dr,n_test_vertices*n_PMTs*sizeof(float),cudaMemcpyHostToDevice));
1274  }
1275  checkCudaErrors( cudaMemcpyToSymbol(constant_time_step_size, &time_step_size, sizeof(time_step_size)) );
1276  checkCudaErrors( cudaMemcpyToSymbol(constant_n_test_vertices, &n_test_vertices, sizeof(n_test_vertices)) );
1277  checkCudaErrors( cudaMemcpyToSymbol(constant_n_water_like_test_vertices, &n_water_like_test_vertices, sizeof(n_water_like_test_vertices)) );
1278  checkCudaErrors( cudaMemcpyToSymbol(constant_n_PMTs, &n_PMTs, sizeof(n_PMTs)) );
1279 
1280  // Bind the array to the texture
1281  checkCudaErrors(cudaBindTexture(0,tex_times_of_flight, device_times_of_flight, n_test_vertices*n_PMTs*sizeof(time_of_flight_t)));
1282  if( correct_mode == 8 && return_direction ){
1283  checkCudaErrors(cudaBindTexture(0,tex_PMTx, device_PMTx, n_PMTs*sizeof(float)));
1284  checkCudaErrors(cudaBindTexture(0,tex_PMTy, device_PMTy, n_PMTs*sizeof(float)));
1285  checkCudaErrors(cudaBindTexture(0,tex_PMTz, device_PMTz, n_PMTs*sizeof(float)));
1286  }
1287  if( correct_mode == 10 ){
1288  checkCudaErrors(cudaBindTexture(0,tex_light_dx, device_light_dx, n_test_vertices*n_PMTs*sizeof(float)));
1289  checkCudaErrors(cudaBindTexture(0,tex_light_dy, device_light_dy, n_test_vertices*n_PMTs*sizeof(float)));
1290  checkCudaErrors(cudaBindTexture(0,tex_light_dz, device_light_dz, n_test_vertices*n_PMTs*sizeof(float)));
1291  checkCudaErrors(cudaBindTexture(0,tex_light_dr, device_light_dr, n_test_vertices*n_PMTs*sizeof(float)));
1292  checkCudaErrors( cudaMemcpyToSymbol(constant_cerenkov_costheta, &cerenkov_costheta, sizeof(float)) );
1293  checkCudaErrors( cudaMemcpyToSymbol(constant_costheta_cone_cut, &costheta_cone_cut, sizeof(float)) );
1294  checkCudaErrors( cudaMemcpyToSymbol(constant_select_based_on_cone, &select_based_on_cone, sizeof(bool)) );
1295  checkCudaErrors( cudaMemcpyToSymbol(constant_return_direction, &return_direction, sizeof(bool)) );
1296  }
1297 
1298 
1299  return;
1300 }
1301 
1302 
1303 void fill_directions_memory_on_device(){
1304 
1305  printf(" [2] --- copy directions from host to device \n");
1306  checkCudaErrors(cudaMemcpy(device_directions_for_vertex_and_pmt,
1307  host_directions_for_vertex_and_pmt,
1308  n_test_vertices*n_PMTs*n_direction_bins*sizeof(bool),
1309  cudaMemcpyHostToDevice));
1310  checkCudaErrors( cudaMemcpyToSymbol(constant_n_direction_bins_theta, &n_direction_bins_theta, sizeof(n_direction_bins_theta)) );
1311  checkCudaErrors( cudaMemcpyToSymbol(constant_n_direction_bins_phi, &n_direction_bins_phi, sizeof(n_direction_bins_phi)) );
1312  checkCudaErrors( cudaMemcpyToSymbol(constant_n_direction_bins, &n_direction_bins, sizeof(n_direction_bins)) );
1313 
1314  // Bind the array to the texture
1315  // checkCudaErrors(cudaBindTexture(0,tex_directions_for_vertex_and_pmt, device_directions_for_vertex_and_pmt, n_test_vertices*n_PMTs*n_direction_bins_theta*n_direction_bins_theta*sizeof(bool)));
1316 
1317 
1318 
1319  return;
1320 }
1321 
1322 
1323 void fill_tofs_memory_on_device_nhits(){
1324 
1325  printf(" [2] --- copy tofs from host to device \n");
1326  checkCudaErrors( cudaMemcpyToSymbol(constant_time_step_size, &time_step_size, sizeof(time_step_size)) );
1327  checkCudaErrors( cudaMemcpyToSymbol(constant_n_PMTs, &n_PMTs, sizeof(n_PMTs)) );
1328 
1329 
1330  return;
1331 }
1332 
1333 
1334 void fill_correct_memory_on_device(){
1335 
1336  printf(" [2] --- copy from host to device \n");
1337  checkCudaErrors(cudaMemcpy(device_ids,
1338  host_ids,
1339  n_hits*sizeof(unsigned int),
1340  cudaMemcpyHostToDevice));
1341  checkCudaErrors(cudaMemcpy(device_times,
1342  host_times,
1343  n_hits*sizeof(unsigned int),
1344  cudaMemcpyHostToDevice));
1345  checkCudaErrors( cudaMemcpyToSymbol(constant_time_offset, &time_offset, sizeof(time_offset)) );
1346 
1347  checkCudaErrors(cudaBindTexture(0,tex_ids, device_ids, n_hits*sizeof(unsigned int)));
1348  checkCudaErrors(cudaBindTexture(0,tex_times, device_times, n_hits*sizeof(unsigned int)));
1349 
1350 
1351  return;
1352 }
1353 
1354 
1355 
1356 
1357 
1358 unsigned int read_number_of_pmts(){
1359 
1360  FILE *f=fopen(pmts_file.c_str(), "r");
1361  if (f == NULL){
1362  printf(" [2] cannot read pmts file %s \n", pmts_file.c_str());
1363  fclose(f);
1364  return 0;
1365  }
1366 
1367  unsigned int n_pmts = 0;
1368 
1369  for (char c = getc(f); c != EOF; c = getc(f))
1370  if (c == '\n')
1371  n_pmts ++;
1372 
1373  fclose(f);
1374  return n_pmts;
1375 
1376 }
1377 
1378 bool read_pmts(){
1379 
1380  FILE *f=fopen(pmts_file.c_str(), "r");
1381 
1382  double x, y, z;
1383  unsigned int id;
1384  for( unsigned int i=0; i<n_PMTs; i++){
1385  if( fscanf(f, "%d %lf %lf %lf", &id, &x, &y, &z) != 4 ){
1386  printf(" [2] problem scanning pmt %d \n", i);
1387  fclose(f);
1388  return false;
1389  }
1390  PMT_x[id-1] = x;
1391  PMT_y[id-1] = y;
1392  PMT_z[id-1] = z;
1393  }
1394 
1395  fclose(f);
1396 
1397 
1398  return true;
1399 
1400 }
1401 
1402 
1403 void coalesce_triggers(){
1404 
1405  trigger_pair_vertex_time.clear();
1406  trigger_npmts_in_time_bin.clear();
1407  if( correct_mode == 8 && return_direction ){
1408  trigger_meanx_per_time_bin.clear();
1409  trigger_meany_per_time_bin.clear();
1410  trigger_meanz_per_time_bin.clear();
1411  }
1412  if( correct_mode == 10 ){
1413  trigger_npmts_in_cone_in_time_bin.clear();
1414  }
1415 
1416  unsigned int vertex_index, time_upper, number_of_pmts_in_time_bin, number_of_pmts_in_cone_in_time_bin;
1417  unsigned int max_pmt=0,max_vertex_index=0,max_time=0,max_pmt_in_cone=0;
1418  bool first_trigger, last_trigger, coalesce_triggers;
1419  unsigned int trigger_index;
1420  float local_trigger_meanx_per_time_bin, local_trigger_meany_per_time_bin, local_trigger_meanz_per_time_bin;
1421  float max_trigger_meanx_per_time_bin = 0, max_trigger_meany_per_time_bin = 0, max_trigger_meanz_per_time_bin = 0;
1422  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=candidate_trigger_pair_vertex_time.begin(); itrigger != candidate_trigger_pair_vertex_time.end(); ++itrigger){
1423 
1424  vertex_index = itrigger->first;
1425  time_upper = itrigger->second;
1426  trigger_index = itrigger - candidate_trigger_pair_vertex_time.begin();
1427  number_of_pmts_in_time_bin = candidate_trigger_npmts_in_time_bin.at(trigger_index);
1428  if( correct_mode == 8 && return_direction ){
1429  local_trigger_meanx_per_time_bin = candidate_trigger_meanx_per_time_bin.at(trigger_index);
1430  local_trigger_meany_per_time_bin = candidate_trigger_meany_per_time_bin.at(trigger_index);
1431  local_trigger_meanz_per_time_bin = candidate_trigger_meanz_per_time_bin.at(trigger_index);
1432  }
1433  if( correct_mode == 10 ){
1434  number_of_pmts_in_cone_in_time_bin = candidate_trigger_npmts_in_cone_in_time_bin.at(trigger_index);
1435  }
1436 
1437  first_trigger = (trigger_index == 0);
1438  last_trigger = (trigger_index == candidate_trigger_pair_vertex_time.size()-1);
1439 
1440  if( first_trigger ){
1441  max_pmt = number_of_pmts_in_time_bin;
1442  max_vertex_index = vertex_index;
1443  max_time = time_upper;
1444  if( correct_mode == 8 && return_direction ){
1445  max_trigger_meanx_per_time_bin = local_trigger_meanx_per_time_bin;
1446  max_trigger_meany_per_time_bin = local_trigger_meany_per_time_bin;
1447  max_trigger_meanz_per_time_bin = local_trigger_meanz_per_time_bin;
1448  }
1449  if( correct_mode == 10 ){
1450  max_pmt_in_cone = number_of_pmts_in_cone_in_time_bin;
1451  }
1452  }
1453  else{
1454  coalesce_triggers = (std::abs((int)(max_time - time_upper)) < coalesce_time/time_step_size);
1455 
1456  if( coalesce_triggers ){
1457  if( number_of_pmts_in_time_bin >= max_pmt) {
1458  max_pmt = number_of_pmts_in_time_bin;
1459  max_vertex_index = vertex_index;
1460  max_time = time_upper;
1461  if( correct_mode == 8 && return_direction ){
1462  max_trigger_meanx_per_time_bin = local_trigger_meanx_per_time_bin;
1463  max_trigger_meany_per_time_bin = local_trigger_meany_per_time_bin;
1464  max_trigger_meanz_per_time_bin = local_trigger_meanz_per_time_bin;
1465  }
1466  if( correct_mode == 10 ){
1467  max_pmt_in_cone = number_of_pmts_in_cone_in_time_bin;
1468  }
1469  }
1470  }else{
1471  trigger_pair_vertex_time.push_back(std::make_pair(max_vertex_index,max_time));
1472  trigger_npmts_in_time_bin.push_back(max_pmt);
1473  max_pmt = number_of_pmts_in_time_bin;
1474  max_vertex_index = vertex_index;
1475  max_time = time_upper;
1476  if( correct_mode == 8 && return_direction ){
1477  trigger_meanx_per_time_bin.push_back(max_trigger_meanx_per_time_bin);
1478  trigger_meany_per_time_bin.push_back(max_trigger_meany_per_time_bin);
1479  trigger_meanz_per_time_bin.push_back(max_trigger_meanz_per_time_bin);
1480  max_trigger_meanx_per_time_bin = local_trigger_meanx_per_time_bin;
1481  max_trigger_meany_per_time_bin = local_trigger_meany_per_time_bin;
1482  max_trigger_meanz_per_time_bin = local_trigger_meanz_per_time_bin;
1483  }
1484  if( correct_mode == 10 ){
1485  trigger_npmts_in_cone_in_time_bin.push_back(max_pmt_in_cone);
1486  max_pmt_in_cone = number_of_pmts_in_cone_in_time_bin;
1487  }
1488  }
1489  }
1490 
1491  if(last_trigger){
1492  trigger_pair_vertex_time.push_back(std::make_pair(max_vertex_index,max_time));
1493  trigger_npmts_in_time_bin.push_back(max_pmt);
1494  if( correct_mode == 8 && return_direction ){
1495  trigger_meanx_per_time_bin.push_back(max_trigger_meanx_per_time_bin);
1496  trigger_meany_per_time_bin.push_back(max_trigger_meany_per_time_bin);
1497  trigger_meanz_per_time_bin.push_back(max_trigger_meanz_per_time_bin);
1498  }
1499  if( correct_mode == 10 ){
1500  trigger_npmts_in_cone_in_time_bin.push_back(max_pmt_in_cone);
1501  }
1502  }
1503 
1504  }
1505 
1506  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=trigger_pair_vertex_time.begin(); itrigger != trigger_pair_vertex_time.end(); ++itrigger)
1507  printf(" [2] coalesced vertex index %d trigger timebin %d \n", itrigger->first, itrigger->second);
1508 
1509  return;
1510 
1511 }
1512 
1513 
1514 void separate_triggers_into_gates(){
1515 
1516  final_trigger_pair_vertex_time.clear();
1517  unsigned int trigger_index;
1518 
1519  unsigned int time_start=0;
1520  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=trigger_pair_vertex_time.begin(); itrigger != trigger_pair_vertex_time.end(); ++itrigger){
1521  //once a trigger is found, we must jump in the future before searching for the next
1522  if(itrigger->second > time_start) {
1523  unsigned int triggertime = itrigger->second*time_step_size - time_offset;
1524  final_trigger_pair_vertex_time.push_back(std::make_pair(itrigger->first,triggertime));
1525  time_start = itrigger->second + trigger_gate_up/time_step_size;
1526  trigger_index = itrigger - trigger_pair_vertex_time.begin();
1527  output_trigger_information.clear();
1528  output_trigger_information.push_back(vertex_x[itrigger->first]);
1529  output_trigger_information.push_back(vertex_y[itrigger->first]);
1530  output_trigger_information.push_back(vertex_z[itrigger->first]);
1531  output_trigger_information.push_back(trigger_npmts_in_time_bin.at(trigger_index));
1532  output_trigger_information.push_back(triggertime);
1533 
1534  printf(" [2] triggertime: %d, vertex index: %d, npmts: %d, x: %f, y: %f, z: %f \n", triggertime, itrigger->first, trigger_npmts_in_time_bin.at(trigger_index), vertex_x[itrigger->first], vertex_y[itrigger->first], vertex_z[itrigger->first]);
1535 
1536  /* if( output_txt ){ */
1537  /* FILE *of=fopen(output_file.c_str(), "w"); */
1538 
1539  /* unsigned int distance_index; */
1540  /* double tof; */
1541  /* double corrected_time; */
1542 
1543  /* for(unsigned int i=0; i<n_hits; i++){ */
1544 
1545  /* distance_index = get_distance_index(host_ids[i], n_PMTs*(itrigger->first)); */
1546  /* tof = host_times_of_flight[distance_index]; */
1547 
1548  /* corrected_time = host_times[i]-tof; */
1549 
1550  /* //fprintf(of, " %d %d %f \n", host_ids[i], host_times[i], corrected_time); */
1551  /* fprintf(of, " %d %f \n", host_ids[i], corrected_time); */
1552  /* } */
1553 
1554  /* fclose(of); */
1555  /* } */
1556 
1557  }
1558  }
1559 
1560 
1561  return;
1562 }
1563 
1564 void separate_triggers_into_gates(std::vector<int> * trigger_ns, std::vector<int> * trigger_ts, std::vector<double> * trigger_vtx_xs, std::vector<double> * trigger_vtx_ys, std::vector<double> * trigger_vtx_zs, std::vector<double> * trigger_dir_xs, std::vector<double> * trigger_dir_ys, std::vector<double> * trigger_dir_zs){
1565 
1566  final_trigger_pair_vertex_time.clear();
1567  unsigned int trigger_index;
1568 
1569  unsigned int time_start=0;
1570  float sum_x, mean_x, sum_y, mean_y, sum_z, mean_z;
1571  float dir_x, dir_y, dir_z;
1572  float dir;
1573  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=trigger_pair_vertex_time.begin(); itrigger != trigger_pair_vertex_time.end(); ++itrigger){
1574  //once a trigger is found, we must jump in the future before searching for the next
1575  if(itrigger->second > time_start) {
1576  unsigned int triggertime = itrigger->second*time_step_size - time_offset;
1577  final_trigger_pair_vertex_time.push_back(std::make_pair(itrigger->first,triggertime));
1578  time_start = itrigger->second + trigger_gate_up/time_step_size;
1579  trigger_index = itrigger - trigger_pair_vertex_time.begin();
1580  output_trigger_information.clear();
1581  output_trigger_information.push_back(vertex_x[itrigger->first]);
1582  output_trigger_information.push_back(vertex_y[itrigger->first]);
1583  output_trigger_information.push_back(vertex_z[itrigger->first]);
1584  output_trigger_information.push_back(trigger_npmts_in_time_bin.at(trigger_index));
1585  output_trigger_information.push_back(triggertime);
1586 
1587  trigger_ns->push_back(trigger_npmts_in_time_bin.at(trigger_index));
1588  trigger_ts->push_back(triggertime);
1589  if( return_vertex ){
1590  trigger_vtx_xs->push_back(vertex_x[itrigger->first]);
1591  trigger_vtx_ys->push_back(vertex_y[itrigger->first]);
1592  trigger_vtx_zs->push_back(vertex_z[itrigger->first]);
1593  }
1594  if( correct_mode == 8 && return_direction ){
1595  sum_x = trigger_meanx_per_time_bin.at(trigger_index);
1596  mean_x = sum_x/trigger_npmts_in_time_bin.at(trigger_index);
1597  dir_x = mean_x - vertex_x[itrigger->first];
1598  sum_y = trigger_meany_per_time_bin.at(trigger_index);
1599  mean_y = sum_y/trigger_npmts_in_time_bin.at(trigger_index);
1600  dir_y = mean_y - vertex_y[itrigger->first];
1601  sum_z = trigger_meanz_per_time_bin.at(trigger_index);
1602  mean_z = sum_z/trigger_npmts_in_time_bin.at(trigger_index);
1603  dir_z = mean_z - vertex_z[itrigger->first];
1604  dir = sqrt(pow(dir_x,2) + pow(dir_y,2) + pow(dir_z,2));
1605  trigger_dir_xs->push_back(dir_x/dir);
1606  trigger_dir_ys->push_back(dir_y/dir);
1607  trigger_dir_zs->push_back(dir_z/dir);
1608  }
1609 
1610  printf(" [2] triggertime: %d, vertex index: %d, npmts: %d, x: %f, y: %f, z: %f \n", triggertime, itrigger->first, trigger_npmts_in_time_bin.at(trigger_index), vertex_x[itrigger->first], vertex_y[itrigger->first], vertex_z[itrigger->first]);
1611 
1612  /* if( output_txt ){ */
1613  /* FILE *of=fopen(output_file.c_str(), "w"); */
1614 
1615  /* unsigned int distance_index; */
1616  /* double tof; */
1617  /* double corrected_time; */
1618 
1619  /* for(unsigned int i=0; i<n_hits; i++){ */
1620 
1621  /* distance_index = get_distance_index(host_ids[i], n_PMTs*(itrigger->first)); */
1622  /* tof = host_times_of_flight[distance_index]; */
1623 
1624  /* corrected_time = host_times[i]-tof; */
1625 
1626  /* //fprintf(of, " %d %d %f \n", host_ids[i], host_times[i], corrected_time); */
1627  /* fprintf(of, " %d %f \n", host_ids[i], corrected_time); */
1628  /* } */
1629 
1630  /* fclose(of); */
1631  /* } */
1632 
1633  }
1634  }
1635 
1636 
1637  return;
1638 }
1639 
1640 
1641 float timedifference_msec(struct timeval t0, struct timeval t1){
1642  return (t1.tv_sec - t0.tv_sec) * 1000.0f + (t1.tv_usec - t0.tv_usec) / 1000.0f;
1643 }
1644 
1645 
1646 
1647 void start_c_clock(){
1648  gettimeofday(&t0,0);
1649 
1650 }
1651 double stop_c_clock(){
1652  gettimeofday(&t1,0);
1653  return timedifference_msec(t0, t1);
1654 }
1655 void start_cuda_clock(){
1656  cudaEventRecord(start);
1657 
1658 }
1659 double stop_cuda_clock(){
1660  float milli;
1661  cudaEventRecord(stop);
1662  cudaEventSynchronize(stop);
1663  cudaEventElapsedTime(&milli, start, stop);
1664  return milli;
1665 }
1666 void start_total_cuda_clock(){
1667  cudaEventRecord(total_start);
1668 
1669 }
1670 double stop_total_cuda_clock(){
1671  float milli;
1672  cudaEventRecord(total_stop);
1673  cudaEventSynchronize(total_stop);
1674  cudaEventElapsedTime(&milli, total_start, total_stop);
1675  return milli;
1676 }
1677 
1678 unsigned int get_distance_index(unsigned int pmt_id, unsigned int vertex_block){
1679  // block = (npmts) * (vertex index)
1680 
1681  return pmt_id - 1 + vertex_block;
1682 
1683 }
1684 
1685 unsigned int get_time_index(unsigned int hit_index, unsigned int vertex_block){
1686  // block = (n time bins) * (vertex index)
1687 
1688  return hit_index + vertex_block;
1689 
1690 }
1691 
1692  unsigned int get_direction_index_at_angles(unsigned int iphi, unsigned int itheta){
1693 
1694  if( itheta == 0 ) return 0;
1695  if( itheta + 1 == n_direction_bins_theta ) return n_direction_bins - 1;
1696 
1697  return 1 + (itheta - 1) * n_direction_bins_phi + iphi;
1698 
1699 }
1700 
1701 unsigned int get_direction_index_at_pmt(unsigned int pmt_id, unsigned int vertex_index, unsigned int direction_index){
1702 
1703  // pmt id 1 ... pmt id p
1704  // [ vertex 1 vertex 2 ... vertex m] ... [vertex 1 ... vertex m]
1705  // [(dir 1 ... dir n) (dir 1 ... dir n) ... (dir 1 ... dir n)] ...
1706 
1707  return n_direction_bins * (pmt_id * n_test_vertices + vertex_index) + direction_index ;
1708 
1709 }
1710 
1711 unsigned int get_direction_index_at_time(unsigned int time_bin, unsigned int vertex_index, unsigned int direction_index){
1712 
1713  // time 1 ... time p
1714  // [ vertex 1 vertex 2 ... vertex m] ... [vertex 1 ... vertex m]
1715  // [(dir 1 ... dir n) (dir 1 ... dir n) ... (dir 1 ... dir n)] ...
1716 
1717  return n_direction_bins* (time_bin * n_test_vertices + vertex_index ) + direction_index ;
1718 
1719 }
1720 
1721 
1722 __device__ unsigned int device_get_distance_index(unsigned int pmt_id, unsigned int vertex_block){
1723  // block = (npmts) * (vertex index)
1724 
1725  return pmt_id - 1 + vertex_block;
1726 
1727 }
1728 
1729 __device__ unsigned int device_get_time_index(unsigned int hit_index, unsigned int vertex_block){
1730  // block = (n time bins) * (vertex index)
1731 
1732  return hit_index + vertex_block;
1733 
1734 }
1735 
1736 __device__ unsigned int device_get_direction_index_at_pmt(unsigned int pmt_id, unsigned int vertex_index, unsigned int direction_index){
1737 
1738  // pmt id 1 ... pmt id p
1739  // [ vertex 1 vertex 2 ... vertex m] ... [vertex 1 ... vertex m]
1740  // [(dir 1 ... dir n) (dir 1 ... dir n) ... (dir 1 ... dir n)] ...
1741 
1742  return constant_n_direction_bins * (pmt_id * constant_n_test_vertices + vertex_index) + direction_index ;
1743 
1744 }
1745 
1746 __device__ unsigned int device_get_direction_index_at_angles(unsigned int iphi, unsigned int itheta){
1747 
1748  if( itheta == 0 ) return 0;
1749  if( itheta + 1 == constant_n_direction_bins_theta ) return constant_n_direction_bins - 1;
1750 
1751  return 1 + (itheta - 1) * constant_n_direction_bins_phi + iphi;
1752 
1753 }
1754 
1755 __device__ unsigned int device_get_direction_index_at_time(unsigned int time_bin, unsigned int vertex_index, unsigned int direction_index){
1756 
1757  // time 1 ... time p
1758  // [ vertex 1 vertex 2 ... vertex m] ... [vertex 1 ... vertex m]
1759  // [(dir 1 ... dir n) (dir 1 ... dir n) ... (dir 1 ... dir n)] ...
1760 
1761  return constant_n_direction_bins* (time_bin * constant_n_test_vertices + vertex_index ) + direction_index ;
1762 
1763 }
1764 
1765 // Print device properties
1766 void print_gpu_properties(){
1767 
1768  int devCount;
1769  cudaGetDeviceCount(&devCount);
1770  printf(" [2] CUDA Device Query...\n");
1771  printf(" [2] There are %d CUDA devices.\n", devCount);
1772  cudaDeviceProp devProp;
1773  for (int i = 0; i < devCount; ++i){
1774  // Get device properties
1775  printf(" [2] CUDA Device #%d\n", i);
1776  cudaGetDeviceProperties(&devProp, i);
1777  printf("Major revision number: %d\n", devProp.major);
1778  printf("Minor revision number: %d\n", devProp.minor);
1779  printf("Name: %s\n", devProp.name);
1780  printf("Total global memory: %lu\n", devProp.totalGlobalMem);
1781  printf("Total shared memory per block: %lu\n", devProp.sharedMemPerBlock);
1782  printf("Total registers per block: %d\n", devProp.regsPerBlock);
1783  printf("Warp size: %d\n", devProp.warpSize);
1784  printf("Maximum memory pitch: %lu\n", devProp.memPitch);
1785  max_n_threads_per_block = devProp.maxThreadsPerBlock;
1786  printf("Maximum threads per block: %d\n", max_n_threads_per_block);
1787  for (int i = 0; i < 3; ++i)
1788  printf("Maximum dimension %d of block: %d\n", i, devProp.maxThreadsDim[i]);
1789  max_n_blocks = devProp.maxGridSize[0];
1790  for (int i = 0; i < 3; ++i)
1791  printf("Maximum dimension %d of grid: %d\n", i, devProp.maxGridSize[i]);
1792  printf("Clock rate: %d\n", devProp.clockRate);
1793  printf("Total constant memory: %lu\n", devProp.totalConstMem);
1794  printf("Texture alignment: %lu\n", devProp.textureAlignment);
1795  printf("Concurrent copy and execution: %s\n", (devProp.deviceOverlap ? "Yes" : "No"));
1796  printf("Number of multiprocessors: %d\n", devProp.multiProcessorCount);
1797  printf("Kernel execution timeout: %s\n", (devProp.kernelExecTimeoutEnabled ?"Yes" : "No"));
1798  printf("Memory Clock Rate (KHz): %d\n", devProp.memoryClockRate);
1799  printf("Memory Bus Width (bits): %d\n", devProp.memoryBusWidth);
1800  printf("Peak Memory Bandwidth (GB/s): %f\n", 2.0*devProp.memoryClockRate*(devProp.memoryBusWidth/8)/1.0e6);
1801  printf("Concurrent kernels: %d\n", devProp.concurrentKernels);
1802  }
1803  size_t available_memory, total_memory;
1804  cudaMemGetInfo(&available_memory, &total_memory);
1805  size_t stack_memory;
1806  cudaDeviceGetLimit(&stack_memory, cudaLimitStackSize);
1807  size_t fifo_memory;
1808  cudaDeviceGetLimit(&fifo_memory, cudaLimitPrintfFifoSize);
1809  size_t heap_memory;
1810  cudaDeviceGetLimit(&heap_memory, cudaLimitMallocHeapSize);
1811  printf(" [2] memgetinfo: available_memory %f MB, total_memory %f MB, stack_memory %f MB, fifo_memory %f MB, heap_memory %f MB \n", (double)available_memory/1.e6, (double)total_memory/1.e6, (double)stack_memory/1.e6, (double)fifo_memory/1.e6, (double)heap_memory/1.e6);
1812 
1813 
1814  return;
1815 }
1816 
1817 
1818 __global__ void kernel_find_vertex_with_max_npmts_in_timebin(histogram_t * np, histogram_t * mnp, unsigned int * vmnp){
1819 
1820 
1821  // get unique id for each thread in each block == time bin
1822  unsigned int time_bin_index = threadIdx.x + blockDim.x*blockIdx.x;
1823 
1824  // skip if thread is assigned to nonexistent time bin
1825  if( time_bin_index >= constant_n_time_bins ) return;
1826 
1827 
1828  unsigned int number_of_pmts_in_time_bin = 0;
1829  unsigned int time_index;
1830  histogram_t max_number_of_pmts_in_time_bin=0;
1831  unsigned int vertex_with_max_n_pmts = 0;
1832 
1833  for(unsigned int iv=0;iv<constant_n_test_vertices;iv++) { // loop over test vertices
1834  // sum the number of hit PMTs in this time window and the next
1835 
1836  time_index = time_bin_index + constant_n_time_bins*iv;
1837  if( time_index >= constant_n_time_bins*constant_n_test_vertices - 1 ) continue;
1838  number_of_pmts_in_time_bin = np[time_index] + np[time_index+1];
1839  if( number_of_pmts_in_time_bin >= max_number_of_pmts_in_time_bin ){
1840  max_number_of_pmts_in_time_bin = number_of_pmts_in_time_bin;
1841  vertex_with_max_n_pmts = iv;
1842  }
1843  }
1844 
1845  mnp[time_bin_index] = max_number_of_pmts_in_time_bin;
1846  vmnp[time_bin_index] = vertex_with_max_n_pmts;
1847 
1848  return;
1849 
1850 }
1851 
1852 
1853 __global__ void kernel_find_vertex_with_max_npmts_in_timebin(histogram_t * np, histogram_t * mnp, unsigned int * vmnp, float * meanxin, float * meanyin, float * meanzin, float * meanxout, float * meanyout, float * meanzout){
1854 
1855 
1856  // get unique id for each thread in each block == time bin
1857  unsigned int time_bin_index = threadIdx.x + blockDim.x*blockIdx.x;
1858 
1859  // skip if thread is assigned to nonexistent time bin
1860  if( time_bin_index >= constant_n_time_bins ) return;
1861 
1862 
1863  unsigned int number_of_pmts_in_time_bin = 0;
1864  unsigned int time_index;
1865  histogram_t max_number_of_pmts_in_time_bin=0;
1866  unsigned int vertex_with_max_n_pmts = 0;
1867  float meanx_for_vertex_with_max_n_pmts = 0., meany_for_vertex_with_max_n_pmts = 0., meanz_for_vertex_with_max_n_pmts = 0.;
1868 
1869  for(unsigned int iv=0;iv<constant_n_test_vertices;iv++) { // loop over test vertices
1870  // sum the number of hit PMTs in this time window and the next
1871 
1872  time_index = time_bin_index + constant_n_time_bins*iv;
1873  if( time_index >= constant_n_time_bins*constant_n_test_vertices - 1 ) continue;
1874  number_of_pmts_in_time_bin = np[time_index] + np[time_index+1];
1875  if( number_of_pmts_in_time_bin >= max_number_of_pmts_in_time_bin ){
1876  max_number_of_pmts_in_time_bin = number_of_pmts_in_time_bin;
1877  meanx_for_vertex_with_max_n_pmts = meanxin[time_index];
1878  meany_for_vertex_with_max_n_pmts = meanyin[time_index];
1879  meanz_for_vertex_with_max_n_pmts = meanzin[time_index];
1880  vertex_with_max_n_pmts = iv;
1881  }
1882  }
1883 
1884  mnp[time_bin_index] = max_number_of_pmts_in_time_bin;
1885  vmnp[time_bin_index] = vertex_with_max_n_pmts;
1886  meanxout[time_bin_index] = meanx_for_vertex_with_max_n_pmts;
1887  meanyout[time_bin_index] = meany_for_vertex_with_max_n_pmts;
1888  meanzout[time_bin_index] = meanz_for_vertex_with_max_n_pmts;
1889 
1890  return;
1891 
1892 }
1893 
1894 
1895 __global__ void kernel_find_vertex_with_max_npmts_and_center_of_mass_in_timebin(histogram_t * np, histogram_t * mnp, unsigned int * vmnp, unsigned int *nc, unsigned int *mnc){
1896 
1897 
1898  // get unique id for each thread in each block == time bin
1899  unsigned int time_bin_index = threadIdx.x + blockDim.x*blockIdx.x;
1900 
1901  // skip if thread is assigned to nonexistent time bin
1902  if( time_bin_index >= constant_n_time_bins ) return;
1903 
1904 
1905  unsigned int number_of_pmts_in_time_bin = 0;
1906  unsigned int number_of_pmts_in_cone_in_time_bin = 0;
1907  unsigned int time_index;
1908  unsigned int max_number_of_pmts_in_time_bin=0;
1909  unsigned int max_number_of_pmts_in_cone_in_time_bin=0;
1910  unsigned int vertex_with_max_n_pmts = 0;
1911 
1912  for(unsigned int iv=0;iv<constant_n_test_vertices;iv++) { // loop over test vertices
1913  // sum the number of hit PMTs in this time window and the next
1914 
1915  time_index = time_bin_index + constant_n_time_bins*iv;
1916  if( time_index >= constant_n_time_bins*constant_n_test_vertices - 1 ) continue;
1917  number_of_pmts_in_time_bin = np[time_index] + np[time_index+1];
1918  number_of_pmts_in_cone_in_time_bin = nc[time_index] + nc[time_index+1];
1919  if( !constant_select_based_on_cone ){
1920  // use this part to select events based on test vertices
1921  if( number_of_pmts_in_time_bin >= max_number_of_pmts_in_time_bin ){
1922  max_number_of_pmts_in_time_bin = number_of_pmts_in_time_bin;
1923  vertex_with_max_n_pmts = iv;
1924  max_number_of_pmts_in_cone_in_time_bin = number_of_pmts_in_cone_in_time_bin;
1925  }
1926  }else{
1927  // use this part to select events based on test vertices and cone
1928  if( number_of_pmts_in_cone_in_time_bin >= max_number_of_pmts_in_cone_in_time_bin ){
1929  max_number_of_pmts_in_cone_in_time_bin = number_of_pmts_in_cone_in_time_bin;
1930  vertex_with_max_n_pmts = iv;
1931  max_number_of_pmts_in_time_bin = number_of_pmts_in_time_bin;
1932  }
1933  }
1934  }
1935 
1936  mnp[time_bin_index] = max_number_of_pmts_in_time_bin;
1937  vmnp[time_bin_index] = vertex_with_max_n_pmts;
1938  mnc[time_bin_index] = max_number_of_pmts_in_cone_in_time_bin;
1939 
1940  return;
1941 
1942 }
1943 
1944 __global__ void kernel_find_vertex_with_max_npmts_in_timebin_and_directionbin(unsigned int * np, histogram_t * mnp, unsigned int * vmnp){
1945 
1946 
1947  // get unique id for each thread in each block == time bin
1948  unsigned int time_bin_index = threadIdx.x + blockDim.x*blockIdx.x;
1949 
1950  // skip if thread is assigned to nonexistent time bin
1951  if( time_bin_index >= constant_n_time_bins - 1 ) return;
1952 
1953 
1954  unsigned int number_of_pmts_in_time_bin = 0;
1955  unsigned int max_number_of_pmts_in_time_bin=0;
1956  unsigned int vertex_with_max_n_pmts = 0;
1957  unsigned int dir_index_1, dir_index_2;
1958 
1959  for(unsigned int iv=0;iv<constant_n_test_vertices;iv++) { // loop over test vertices
1960  // sum the number of hit PMTs in this time window
1961 
1962  for(unsigned int idir = 0; idir < constant_n_direction_bins; idir++){
1963 
1964  dir_index_1 = device_get_direction_index_at_time(time_bin_index, iv, idir);
1965  dir_index_2 = device_get_direction_index_at_time(time_bin_index + 1, iv, idir);
1966 
1967  number_of_pmts_in_time_bin = np[dir_index_1]
1968  + np[dir_index_2];
1969  if( number_of_pmts_in_time_bin > max_number_of_pmts_in_time_bin ){
1970  max_number_of_pmts_in_time_bin = number_of_pmts_in_time_bin;
1971  vertex_with_max_n_pmts = iv;
1972  }
1973  }
1974  }
1975 
1976  mnp[time_bin_index] = max_number_of_pmts_in_time_bin;
1977  vmnp[time_bin_index] = vertex_with_max_n_pmts;
1978 
1979  return;
1980 
1981 }
1982 
1983 void free_event_memories(){
1984 
1985  checkCudaErrors(cudaUnbindTexture(tex_ids));
1986  checkCudaErrors(cudaUnbindTexture(tex_times));
1987  free(host_ids);
1988  free(host_times);
1989  checkCudaErrors(cudaFree(device_ids));
1990  checkCudaErrors(cudaFree(device_times));
1991  if( correct_mode == 1 ){
1992  checkCudaErrors(cudaFree(device_time_bin_of_hit));
1993  } else if( correct_mode == 2 ){
1994  checkCudaErrors(cudaFree(device_time_bin_of_hit));
1995  free(host_time_bin_of_hit);
1996  free(host_n_pmts_per_time_bin);
1997  } else if( correct_mode == 3 ){
1998  checkCudaErrors(cudaFree(device_time_bin_of_hit));
1999  } else if( correct_mode == 4 ){
2000  checkCudaErrors(cudaFree(device_time_bin_of_hit));
2001  } else if( correct_mode == 5 ){
2002  checkCudaErrors(cudaFree(device_time_bin_of_hit));
2003  } else if( correct_mode == 6 ){
2004  checkCudaErrors(cudaFree(device_time_bin_of_hit));
2005  } else if( correct_mode == 7 ){
2006  checkCudaErrors(cudaFree(device_time_bin_of_hit));
2007  }
2008  if( correct_mode != 9 ){
2009  checkCudaErrors(cudaFree(device_n_pmts_per_time_bin));
2010  if( correct_mode == 8 && return_direction ){
2011  checkCudaErrors(cudaFree(device_meanx_per_time_bin));
2012  checkCudaErrors(cudaFree(device_meany_per_time_bin));
2013  checkCudaErrors(cudaFree(device_meanz_per_time_bin));
2014  }
2015  if( correct_mode == 10 ){
2016  checkCudaErrors(cudaFree(device_dx_per_time_bin));
2017  checkCudaErrors(cudaFree(device_dy_per_time_bin));
2018  checkCudaErrors(cudaFree(device_dz_per_time_bin));
2019  }
2020  }else{
2021  checkCudaErrors(cudaFree(device_n_pmts_per_time_bin_and_direction_bin));
2022  }
2023  free(host_max_number_of_pmts_in_time_bin);
2024  free(host_vertex_with_max_n_pmts);
2025  checkCudaErrors(cudaFree(device_max_number_of_pmts_in_time_bin));
2026  checkCudaErrors(cudaFree(device_vertex_with_max_n_pmts));
2027 
2028  if( correct_mode == 8 && return_direction ){
2029  free(host_meanx_for_vertex_with_max_n_pmts_per_time_bin);
2030  free(host_meany_for_vertex_with_max_n_pmts_per_time_bin);
2031  free(host_meanz_for_vertex_with_max_n_pmts_per_time_bin);
2032  checkCudaErrors(cudaFree(device_meanx_for_vertex_with_max_n_pmts_per_time_bin));
2033  checkCudaErrors(cudaFree(device_meany_for_vertex_with_max_n_pmts_per_time_bin));
2034  checkCudaErrors(cudaFree(device_meanz_for_vertex_with_max_n_pmts_per_time_bin));
2035  }
2036  if( correct_mode == 10 ){
2037  free(host_max_number_of_pmts_in_cone_in_time_bin);
2038  checkCudaErrors(cudaFree(device_max_number_of_pmts_in_cone_in_time_bin));
2039  checkCudaErrors(cudaFree(device_number_of_pmts_in_cone_in_time_bin));
2040  }
2041 
2042  return;
2043 }
2044 
2045 
2046 void free_event_memories_nhits(){
2047 
2048  checkCudaErrors(cudaUnbindTexture(tex_ids));
2049  checkCudaErrors(cudaUnbindTexture(tex_times));
2050  free(host_ids);
2051  free(host_times);
2052  checkCudaErrors(cudaFree(device_ids));
2053  checkCudaErrors(cudaFree(device_times));
2054  checkCudaErrors(cudaFree(device_n_pmts_nhits));
2055  // checkCudaErrors(cudaFree(device_time_nhits));
2056  free(host_n_pmts_nhits);
2057  // free(host_time_nhits);
2058 
2059  return;
2060 }
2061 
2062 void free_global_memories(){
2063 
2064  //unbind texture reference to free resource
2065  checkCudaErrors(cudaUnbindTexture(tex_times_of_flight));
2066  if( correct_mode == 8 && return_direction ){
2067  checkCudaErrors(cudaUnbindTexture(tex_PMTx));
2068  checkCudaErrors(cudaUnbindTexture(tex_PMTy));
2069  checkCudaErrors(cudaUnbindTexture(tex_PMTz));
2070  }
2071  if( correct_mode == 10 ){
2072  checkCudaErrors(cudaUnbindTexture(tex_light_dx));
2073  checkCudaErrors(cudaUnbindTexture(tex_light_dy));
2074  checkCudaErrors(cudaUnbindTexture(tex_light_dz));
2075  checkCudaErrors(cudaUnbindTexture(tex_light_dr));
2076  }
2077 
2078  if( correct_mode == 9 ){
2079  // checkCudaErrors(cudaUnbindTexture(tex_directions_for_vertex_and_pmt));
2080  checkCudaErrors(cudaFree(device_directions_for_vertex_and_pmt));
2081  free(host_directions_for_vertex_and_pmt);
2082  }
2083 
2084  free(PMT_x);
2085  free(PMT_y);
2086  free(PMT_z);
2087  free(vertex_x);
2088  free(vertex_y);
2089  free(vertex_z);
2090  checkCudaErrors(cudaFree(device_times_of_flight));
2091  free(host_times_of_flight);
2092  if( correct_mode == 8 && return_direction ){
2093  checkCudaErrors(cudaFree(device_PMTx));
2094  free(host_PMTx);
2095  checkCudaErrors(cudaFree(device_PMTy));
2096  free(host_PMTy);
2097  checkCudaErrors(cudaFree(device_PMTz));
2098  free(host_PMTz);
2099  }
2100  if( correct_mode == 10 ){
2101  checkCudaErrors(cudaFree(device_light_dx));
2102  free(host_light_dx);
2103  checkCudaErrors(cudaFree(device_light_dy));
2104  free(host_light_dy);
2105  checkCudaErrors(cudaFree(device_light_dz));
2106  free(host_light_dz);
2107  checkCudaErrors(cudaFree(device_light_dr));
2108  free(host_light_dr);
2109  }
2110 
2111  return;
2112 }
2113 
2114 void copy_candidates_from_device_to_host(){
2115 
2116  checkCudaErrors(cudaMemcpy(host_max_number_of_pmts_in_time_bin,
2117  device_max_number_of_pmts_in_time_bin,
2118  n_time_bins*sizeof(histogram_t),
2119  cudaMemcpyDeviceToHost));
2120  checkCudaErrors(cudaMemcpy(host_vertex_with_max_n_pmts,
2121  device_vertex_with_max_n_pmts,
2122  n_time_bins*sizeof(unsigned int),
2123  cudaMemcpyDeviceToHost));
2124  if( correct_mode == 8 && return_direction ){
2125  checkCudaErrors(cudaMemcpy(host_meanx_for_vertex_with_max_n_pmts_per_time_bin,
2126  device_meanx_for_vertex_with_max_n_pmts_per_time_bin,
2127  n_time_bins*sizeof(float),
2128  cudaMemcpyDeviceToHost));
2129  checkCudaErrors(cudaMemcpy(host_meany_for_vertex_with_max_n_pmts_per_time_bin,
2130  device_meany_for_vertex_with_max_n_pmts_per_time_bin,
2131  n_time_bins*sizeof(float),
2132  cudaMemcpyDeviceToHost));
2133  checkCudaErrors(cudaMemcpy(host_meanz_for_vertex_with_max_n_pmts_per_time_bin,
2134  device_meanz_for_vertex_with_max_n_pmts_per_time_bin,
2135  n_time_bins*sizeof(float),
2136  cudaMemcpyDeviceToHost));
2137  }
2138  if( correct_mode == 10 ){
2139  checkCudaErrors(cudaMemcpy(host_max_number_of_pmts_in_cone_in_time_bin,
2140  device_max_number_of_pmts_in_cone_in_time_bin,
2141  n_time_bins*sizeof(unsigned int),
2142  cudaMemcpyDeviceToHost));
2143 
2144  }
2145 
2146 }
2147 
2148 
2149 void choose_candidates_above_threshold(){
2150 
2151  candidate_trigger_pair_vertex_time.clear();
2152  candidate_trigger_npmts_in_time_bin.clear();
2153  if( correct_mode == 8 && return_direction ){
2154  candidate_trigger_meanx_per_time_bin.clear();
2155  candidate_trigger_meany_per_time_bin.clear();
2156  candidate_trigger_meanz_per_time_bin.clear();
2157  }
2158  if( correct_mode == 10 ){
2159  candidate_trigger_npmts_in_cone_in_time_bin.clear();
2160  }
2161 
2162  unsigned int the_threshold;
2163  unsigned int number_of_pmts_to_cut_on;
2164 
2165  for(unsigned int time_bin = 0; time_bin<n_time_bins - 1; time_bin++){ // loop over time bins
2166  // n_time_bins - 1 as we are checking the i and i+1 at the same time
2167 
2168  if(host_vertex_with_max_n_pmts[time_bin] < n_water_like_test_vertices )
2169  the_threshold = water_like_threshold_number_of_pmts;
2170  else
2171  the_threshold = wall_like_threshold_number_of_pmts;
2172 
2173  number_of_pmts_to_cut_on = host_max_number_of_pmts_in_time_bin[time_bin];
2174  if( correct_mode == 10 ){
2175  if( select_based_on_cone ){
2176  number_of_pmts_to_cut_on = host_max_number_of_pmts_in_cone_in_time_bin[time_bin];
2177  }
2178  }
2179 
2180  if(number_of_pmts_to_cut_on > the_threshold) {
2181 
2182  if( use_verbose ){
2183  printf(" [2] time %f vertex (%f, %f, %f) npmts %d \n", (time_bin + 2)*time_step_size - time_offset, vertex_x[host_vertex_with_max_n_pmts[time_bin]], vertex_y[host_vertex_with_max_n_pmts[time_bin]], vertex_z[host_vertex_with_max_n_pmts[time_bin]], number_of_pmts_to_cut_on);
2184  }
2185 
2186  candidate_trigger_pair_vertex_time.push_back(std::make_pair(host_vertex_with_max_n_pmts[time_bin],time_bin+1));
2187  candidate_trigger_npmts_in_time_bin.push_back(host_max_number_of_pmts_in_time_bin[time_bin]);
2188  if( correct_mode == 8 && return_direction ){
2189  candidate_trigger_meanx_per_time_bin.push_back(host_meanx_for_vertex_with_max_n_pmts_per_time_bin[time_bin]);
2190  candidate_trigger_meany_per_time_bin.push_back(host_meany_for_vertex_with_max_n_pmts_per_time_bin[time_bin]);
2191  candidate_trigger_meanz_per_time_bin.push_back(host_meanz_for_vertex_with_max_n_pmts_per_time_bin[time_bin]);
2192  }
2193  if( correct_mode == 10 ){
2194  candidate_trigger_npmts_in_cone_in_time_bin.push_back(host_max_number_of_pmts_in_cone_in_time_bin[time_bin]);
2195  }
2196  }
2197 
2198  }
2199 
2200  if( use_verbose )
2201  printf(" [2] n candidates: %d \n", candidate_trigger_pair_vertex_time.size());
2202 }
2203 
2204 bool set_input_file_for_event(int n){
2205 
2206  int nchar = (ceil(log10(n+1))+1);
2207  char * num = (char*)malloc(sizeof(char)*nchar);
2208  sprintf(num, "%d", n+1);
2209  event_file = event_file_base + num + event_file_suffix;
2210 
2211  bool file_exists = (access( event_file.c_str(), F_OK ) != -1);
2212 
2213  free(num);
2214 
2215  return file_exists;
2216 
2217 }
2218 
2219 void set_output_file(){
2220 
2221  int nchar = (ceil(log10(water_like_threshold_number_of_pmts))+1);
2222  char * num = (char*)malloc(sizeof(char)*nchar);
2223  sprintf(num, "%d", water_like_threshold_number_of_pmts);
2224  output_file = output_file_base + num + event_file_suffix;
2225 
2226  free(num);
2227 
2228  return ;
2229 
2230 }
2231 
2232 void set_output_file_nhits(unsigned int threshold){
2233 
2234  int nchar = (ceil(log10(threshold))+1);
2235  char * num = (char*)malloc(sizeof(char)*nchar);
2236  sprintf(num, "%d", threshold);
2237  output_file = output_file_base + num + event_file_suffix;
2238 
2239  free(num);
2240 
2241  return ;
2242 
2243 }
2244 
2245 void write_output_nhits(unsigned int * ntriggers){
2246 
2247  if( output_txt ){
2248 
2249  for(unsigned int u=nhits_threshold_min; u<=nhits_threshold_max; u++){
2250  set_output_file_nhits(u);
2251  FILE *of=fopen(output_file.c_str(), "a");
2252 
2253  // int trigger = (ntriggers[u - nhits_threshold_min] > 0 ? 1 : 0);
2254  int trigger = ntriggers[u - nhits_threshold_min];
2255  fprintf(of, " %d \n", trigger);
2256 
2257  fclose(of);
2258  }
2259 
2260  }
2261 
2262 
2263 }
2264 
2265 
2266 void write_output(){
2267 
2268  if( output_txt ){
2269  FILE *of=fopen(output_file.c_str(), "a");
2270 
2271  int trigger;
2272  if( write_output_mode == 0 ){
2273  // output 1 if there is a trigger, 0 otherwise
2274  trigger = (trigger_pair_vertex_time.size() > 0 ? 1 : 0);
2275  }
2276 
2277  if( write_output_mode == 1 ){
2278  // output the n of triggers
2279  trigger = trigger_pair_vertex_time.size();
2280  }
2281 
2282  if( write_output_mode == 2 ){
2283  // output the n of water-like triggers
2284  int trigger = 0;
2285  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=trigger_pair_vertex_time.begin(); itrigger != trigger_pair_vertex_time.end(); ++itrigger){
2286  if( itrigger->first < n_water_like_test_vertices )
2287  trigger ++;
2288  }
2289  }
2290 
2291  if( write_output_mode == 0 || write_output_mode == 1 || write_output_mode == 2 ){
2292  fprintf(of, " %d \n", trigger);
2293  }
2294 
2295  if( write_output_mode == 3 ){
2296  unsigned int triggertime, trigger_index;
2297  // output reconstructed vertices
2298  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=trigger_pair_vertex_time.begin(); itrigger != trigger_pair_vertex_time.end(); ++itrigger){
2299  triggertime = itrigger->second*time_step_size - time_offset;
2300  if( correct_mode == 10 ){
2301  trigger_index = itrigger - trigger_pair_vertex_time.begin();
2302  fprintf(of, " %d %f %f %f %d %d %d \n", n_events, vertex_x[itrigger->first], vertex_y[itrigger->first], vertex_z[itrigger->first], triggertime, trigger_npmts_in_time_bin.at(trigger_index), trigger_npmts_in_cone_in_time_bin.at(trigger_index));
2303  }else{
2304  fprintf(of, " %d %f %f %f %d \n", n_events, vertex_x[itrigger->first], vertex_y[itrigger->first], vertex_z[itrigger->first], triggertime);
2305  }
2306  }
2307  }
2308 
2309  if( write_output_mode == 4 ){
2310  // output non-corrected and corrected times for best vertex
2311  int max_n_pmts = 0;
2312  unsigned int best_vertex;
2313  for(std::vector<std::pair<unsigned int,unsigned int> >::const_iterator itrigger=trigger_pair_vertex_time.begin(); itrigger != trigger_pair_vertex_time.end(); ++itrigger){
2314  unsigned int vertex_index = itrigger->first;
2315  unsigned int time_index = itrigger->second;
2316  unsigned int local_n_pmts = host_max_number_of_pmts_in_time_bin[itrigger->second];
2317  if( local_n_pmts > max_n_pmts ){
2318  max_n_pmts = local_n_pmts;
2319  best_vertex = vertex_index;
2320  }
2321  }
2322  unsigned int distance_index;
2323  double tof;
2324  double corrected_time;
2325 
2326  for(unsigned int i=0; i<n_hits; i++){
2327 
2328  distance_index = get_distance_index(host_ids[i], n_PMTs*best_vertex);
2329  tof = host_times_of_flight[distance_index];
2330  corrected_time = host_times[i]-tof;
2331 
2332  fprintf(of, " %d %d %f \n", host_ids[i], host_times[i], corrected_time);
2333  //fprintf(of, " %d %f \n", host_ids[i], corrected_time);
2334  }
2335  }
2336 
2337  fclose(of);
2338  }
2339 
2340 
2341 }
2342 
2343 
2344 void initialize_output(){
2345 
2346  if( output_txt )
2347  remove( output_file.c_str() );
2348 
2349 }
2350 
2351 void initialize_output_nhits(){
2352 
2353  if( output_txt )
2354  for(unsigned int u=nhits_threshold_min; u<=nhits_threshold_max; u++){
2355  set_output_file_nhits(u);
2356  remove( output_file.c_str() );
2357  }
2358 }
2359 
2360 float read_value_from_file(std::string paramname, std::string filename){
2361 
2362  FILE * pFile = fopen (filename.c_str(),"r");
2363  if(pFile == NULL) {
2364  printf("Error: file %s could not be opened\n", filename.c_str());
2365  }
2366 
2367  char name[256];
2368  float value;
2369 
2370  while( EOF != fscanf(pFile, "%s %e", name, &value) ){
2371  if( paramname.compare(name) == 0 ){
2372  fclose(pFile);
2373  return value;
2374  }
2375  }
2376 
2377  printf(" [2] warning: could not find parameter %s in file %s \n", paramname.c_str(), filename.c_str());
2378 
2379  fclose(pFile);
2380  exit(0);
2381 
2382  return 0.;
2383 
2384 }
2385 
2386 void read_user_parameters(){
2387 
2388  std::string parameter_file = "parameters.txt";
2389 
2390  twopi = 2.*acos(-1.);
2391  speed_light_water = 29.9792/1.3330; // speed of light in water, cm/ns
2392  //double speed_light_water = 22.490023;
2393 
2394  cerenkov_costheta =1./1.3330;
2395  cerenkov_angle_water = acos(cerenkov_costheta);
2396  costheta_cone_cut = read_value_from_file("costheta_cone_cut", parameter_file);
2397  select_based_on_cone = (bool)read_value_from_file("select_based_on_cone", parameter_file);
2398  return_vertex = (bool)read_value_from_file("return_vertex", parameter_file);
2399  return_direction = (bool)read_value_from_file("return_direction", parameter_file);
2400 
2401  dark_rate = read_value_from_file("dark_rate", parameter_file); // Hz
2402  cylindrical_grid = (bool)read_value_from_file("cylindrical_grid", parameter_file);
2403  distance_between_vertices = read_value_from_file("distance_between_vertices", parameter_file); // cm
2404  wall_like_distance = read_value_from_file("wall_like_distance", parameter_file); // units of distance between vertices
2405  time_step_size = (unsigned int)(sqrt(3.)*distance_between_vertices/(4.*speed_light_water)); // ns
2406  int extra_threshold = (int)(dark_rate*n_PMTs*2.*time_step_size*1.e-9); // to account for dark current occupancy
2407  water_like_threshold_number_of_pmts = read_value_from_file("water_like_threshold_number_of_pmts", parameter_file) + extra_threshold;
2408  wall_like_threshold_number_of_pmts = read_value_from_file("wall_like_threshold_number_of_pmts", parameter_file) + extra_threshold;
2409  coalesce_time = read_value_from_file("coalesce_time", parameter_file); // ns
2410  trigger_gate_up = read_value_from_file("trigger_gate_up", parameter_file); // ns
2411  trigger_gate_down = read_value_from_file("trigger_gate_down", parameter_file); // ns
2412  max_n_hits_per_job = read_value_from_file("max_n_hits_per_job", parameter_file);
2413  output_txt = (bool)read_value_from_file("output_txt", parameter_file);
2414  correct_mode = read_value_from_file("correct_mode", parameter_file);
2415  write_output_mode = read_value_from_file("write_output_mode", parameter_file);
2416  number_of_kernel_blocks_3d.y = read_value_from_file("num_blocks_y", parameter_file);
2417  number_of_threads_per_block_3d.y = read_value_from_file("num_threads_per_block_y", parameter_file);
2418  number_of_threads_per_block_3d.x = read_value_from_file("num_threads_per_block_x", parameter_file);
2419 
2420  n_direction_bins_theta = read_value_from_file("n_direction_bins_theta", parameter_file);
2421  n_direction_bins_phi = 2*(n_direction_bins_theta - 1);
2422  n_direction_bins = n_direction_bins_phi*n_direction_bins_theta - 2*(n_direction_bins_phi - 1);
2423 
2424 
2425 }
2426 
2427 
2428 void read_user_parameters_nhits(){
2429 
2430  std::string parameter_file = "parameters.txt";
2431 
2432  twopi = 2.*acos(-1.);
2433  speed_light_water = 29.9792/1.3330; // speed of light in water, cm/ns
2434  //double speed_light_water = 22.490023;
2435 
2436  double dark_rate = read_value_from_file("dark_rate", parameter_file); // Hz
2437  distance_between_vertices = read_value_from_file("distance_between_vertices", parameter_file); // cm
2438  wall_like_distance = read_value_from_file("wall_like_distance", parameter_file); // units of distance between vertices
2439  time_step_size = read_value_from_file("nhits_step_size", parameter_file); // ns
2440  nhits_window = read_value_from_file("nhits_window", parameter_file); // ns
2441  int extra_threshold = (int)(dark_rate*n_PMTs*nhits_window*1.e-9); // to account for dark current occupancy
2442  extra_threshold = 0;
2443  water_like_threshold_number_of_pmts = read_value_from_file("water_like_threshold_number_of_pmts", parameter_file) + extra_threshold;
2444  wall_like_threshold_number_of_pmts = read_value_from_file("wall_like_threshold_number_of_pmts", parameter_file) + extra_threshold;
2445  coalesce_time = read_value_from_file("coalesce_time", parameter_file); // ns
2446  trigger_gate_up = read_value_from_file("trigger_gate_up", parameter_file); // ns
2447  trigger_gate_down = read_value_from_file("trigger_gate_down", parameter_file); // ns
2448  max_n_hits_per_job = read_value_from_file("max_n_hits_per_job", parameter_file);
2449  output_txt = (bool)read_value_from_file("output_txt", parameter_file);
2450  correct_mode = read_value_from_file("correct_mode", parameter_file);
2451  number_of_kernel_blocks_3d.y = read_value_from_file("num_blocks_y", parameter_file);
2452  number_of_threads_per_block_3d.y = read_value_from_file("num_threads_per_block_y", parameter_file);
2453  number_of_threads_per_block_3d.x = read_value_from_file("num_threads_per_block_x", parameter_file);
2454  nhits_threshold_min = read_value_from_file("nhits_threshold_min", parameter_file);
2455  nhits_threshold_max = read_value_from_file("nhits_threshold_max", parameter_file);
2456 
2457 }
2458 
2459 
2460 void check_cudamalloc_float(unsigned int size){
2461 
2462  unsigned int bytes_per_float = 4;
2463  size_t available_memory, total_memory;
2464  cudaMemGetInfo(&available_memory, &total_memory);
2465  if( size*bytes_per_float > available_memory*1000/1024 ){
2466  printf(" [2] cannot allocate %d floats, or %d B, available %d B \n",
2467  size, size*bytes_per_float, available_memory*1000/1024);
2468  }
2469 
2470 }
2471 
2472 void check_cudamalloc_int(unsigned int size){
2473 
2474  unsigned int bytes_per_int = 4;
2475  size_t available_memory, total_memory;
2476  cudaMemGetInfo(&available_memory, &total_memory);
2477  if( size*bytes_per_int > available_memory*1000/1024 ){
2478  printf(" [2] cannot allocate %d ints, or %d B, available %d B \n",
2479  size, size*bytes_per_int, available_memory*1000/1024);
2480  }
2481 
2482 }
2483 
2484 void check_cudamalloc_unsigned_int(unsigned int size){
2485 
2486  unsigned int bytes_per_unsigned_int = 4;
2487  size_t available_memory, total_memory;
2488  cudaMemGetInfo(&available_memory, &total_memory);
2489  if( size*bytes_per_unsigned_int > available_memory*1000/1024 ){
2490  printf(" [2] cannot allocate %d unsigned_ints, or %d B, available %d B \n",
2491  size, size*bytes_per_unsigned_int, available_memory*1000/1024);
2492  }
2493 
2494 }
2495 
2496 
2497 void check_cudamalloc_unsigned_short(unsigned int size){
2498 
2499  unsigned int bytes_per_unsigned_short = 2;
2500  size_t available_memory, total_memory;
2501  cudaMemGetInfo(&available_memory, &total_memory);
2502  if( size*bytes_per_unsigned_short > available_memory*1000/1024 ){
2503  printf(" cannot allocate %d unsigned_shorts, or %d B, available %d B \n",
2504  size, size*bytes_per_unsigned_short, available_memory*1000/1024);
2505  }
2506 
2507 }
2508 
2509 void check_cudamalloc_unsigned_char(unsigned int size){
2510 
2511  unsigned int bytes_per_unsigned_char = 1;
2512  size_t available_memory, total_memory;
2513  cudaMemGetInfo(&available_memory, &total_memory);
2514  if( size*bytes_per_unsigned_char > available_memory*1000/1024 ){
2515  printf(" cannot allocate %d unsigned_chars, or %d B, available %d B \n",
2516  size, size*bytes_per_unsigned_char, available_memory*1000/1024);
2517  }
2518 
2519 }
2520 
2521 void check_cudamalloc_bool(unsigned int size){
2522 
2523  unsigned int bytes_per_bool = 1;
2524  size_t available_memory, total_memory;
2525  cudaMemGetInfo(&available_memory, &total_memory);
2526  if( size*bytes_per_bool > available_memory*1000/1024 ){
2527  printf(" [2] cannot allocate %d unsigned_ints, or %d B, available %d B \n",
2528  size, size*bytes_per_bool, available_memory*1000/1024);
2529  }
2530 
2531 }
2532 
2533 unsigned int find_greatest_divisor(unsigned int n, unsigned int max){
2534 
2535  if( n == 1 ){
2536  return 1;
2537  }
2538 
2539  if (n % 2 == 0){
2540  if( n <= 2*max ){
2541  return n / 2;
2542  }
2543  else{
2544  float sqrtN = sqrt(n); // square root of n in float precision.
2545  unsigned int start = ceil(std::max((double)2,(double)n/(double)max));
2546  for(unsigned int d = start; d <= n; d += 1)
2547  if (n % d == 0)
2548  return n/d;
2549  return 1;
2550  }
2551  }
2552 
2553  // Now, the least prime divisor of n is odd.
2554  // So, we increment the counter by 2 in the loop, by starting in 3.
2555 
2556  float sqrtN = sqrt(n); // square root of n in float precision.
2557  unsigned int start = ceil(std::max((double)3,(double)n/(double)max));
2558  for(unsigned int d = start; d <= n; d += 2)
2559  if (n % d == 0)
2560  return n/d;
2561 
2562  // If the loop has reached its end normally,
2563  // it means that N is prime.
2564 
2565  return 1;
2566 
2567 }
2568 
2569 
2570 void setup_threads_for_histo(unsigned int n){
2571 
2572  number_of_kernel_blocks_3d.x = n/greatest_divisor;
2573  number_of_kernel_blocks_3d.y = 1;
2574 
2575  number_of_threads_per_block_3d.x = greatest_divisor;
2576  number_of_threads_per_block_3d.y = 1;
2577 
2578 }
2579 
2580 void setup_threads_for_histo(){
2581 
2582  number_of_kernel_blocks_3d.x = 1000;
2583  number_of_kernel_blocks_3d.y = 1;
2584 
2585  number_of_threads_per_block_3d.x = max_n_threads_per_block;
2586  number_of_threads_per_block_3d.y = 1;
2587 
2588 }
2589 
2590 void setup_threads_for_histo_iterated(bool last){
2591 
2592  number_of_kernel_blocks_3d.x = 1000;
2593  number_of_kernel_blocks_3d.y = 1;
2594 
2595  unsigned int size = n_time_bins*n_test_vertices;
2596  number_of_threads_per_block_3d.x = (last ? size - size/max_n_threads_per_block*max_n_threads_per_block : max_n_threads_per_block);
2597  number_of_threads_per_block_3d.y = 1;
2598 
2599 }
2600 
2601 void setup_threads_for_histo_per(unsigned int n){
2602 
2603  number_of_kernel_blocks_3d.x = n;
2604 
2605  print_parameters_2d();
2606 
2607  if( number_of_threads_per_block_3d.x * number_of_threads_per_block_3d.y > max_n_threads_per_block ){
2608  printf(" [2] --------------------- warning: number_of_threads_per_block (x*y) = %d cannot exceed max value %d \n", number_of_threads_per_block_3d.x * number_of_threads_per_block_3d.y, max_n_threads_per_block );
2609  }
2610 
2611  if( number_of_kernel_blocks_3d.x > max_n_blocks ){
2612  printf(" [2] warning: number_of_kernel_blocks x = %d cannot exceed max value %d \n", number_of_kernel_blocks_3d.x, max_n_blocks );
2613  }
2614 
2615  if( number_of_kernel_blocks_3d.y > max_n_blocks ){
2616  printf(" [2] warning: number_of_kernel_blocks y = %d cannot exceed max value %d \n", number_of_kernel_blocks_3d.y, max_n_blocks );
2617  }
2618 
2619  if( std::numeric_limits<int>::max() / (number_of_kernel_blocks_3d.x*number_of_kernel_blocks_3d.y) < number_of_threads_per_block_3d.x*number_of_threads_per_block_3d.y ){
2620  printf(" [2] --------------------- warning: grid size cannot exceed max value %u \n", std::numeric_limits<int>::max() );
2621  }
2622 
2623 }
2624 
2625 
2626 #endif
2627 
candidate_trigger_npmts_in_time_bin
std::vector< unsigned int > candidate_trigger_npmts_in_time_bin
Definition: library_daq.h:154
constant_n_hits
__constant__ unsigned int constant_n_hits
Definition: library_daq.h:94
set_output_file_nhits
void set_output_file_nhits(unsigned int threshold)
Definition: library_daq.h:2232
device_PMTz
float * device_PMTz
Definition: library_daq.h:139
device_meanx_for_vertex_with_max_n_pmts_per_time_bin
float * device_meanx_for_vertex_with_max_n_pmts_per_time_bin
Definition: library_daq.h:184
allocate_candidates_memory_on_device
void allocate_candidates_memory_on_device()
Definition: library_daq.h:1134
free_global_memories
void free_global_memories()
Definition: library_daq.h:2062
device_directions_for_vertex_and_pmt
bool * device_directions_for_vertex_and_pmt
Definition: library_daq.h:141
the_max_time
unsigned int the_max_time
Definition: library_daq.h:211
setup_threads_for_histo
void setup_threads_for_histo(unsigned int n)
Definition: library_daq.h:2570
set_output_file
void set_output_file()
Definition: library_daq.h:2219
tex_PMTx
texture< float, 1, cudaReadModeElementType > tex_PMTx
Definition: library_daq.h:148
read_user_parameters_nhits
void read_user_parameters_nhits()
Definition: library_daq.h:2428
device_get_direction_index_at_time
__device__ unsigned int device_get_direction_index_at_time(unsigned int time_bin, unsigned int vertex_index, unsigned int direction_index)
Definition: library_daq.h:1755
allocate_directions_memory_on_device
void allocate_directions_memory_on_device()
Definition: library_daq.h:957
histogram_t
unsigned short histogram_t
Definition: library_daq.h:27
elapsed_candidates_copy_host
float elapsed_candidates_copy_host
Definition: library_daq.h:204
allocate_candidates_memory_on_host
void allocate_candidates_memory_on_host()
Definition: library_daq.h:1113
nhits_threshold_max
unsigned int nhits_threshold_max
Definition: library_daq.h:50
stop_total_cuda_clock
double stop_total_cuda_clock()
Definition: library_daq.h:1670
n_test_vertices
unsigned int n_test_vertices
vertices
Definition: library_daq.h:71
distance_between_vertices
double distance_between_vertices
parameters
Definition: library_daq.h:44
nhits_window
unsigned int nhits_window
Definition: library_daq.h:212
device_dy_per_time_bin
float * device_dy_per_time_bin
Definition: library_daq.h:111
PMT_x
double * PMT_x
Definition: library_daq.h:69
wall_like_threshold_number_of_pmts
unsigned int wall_like_threshold_number_of_pmts
Definition: library_daq.h:49
device_times_of_flight
time_of_flight_t * device_times_of_flight
Definition: library_daq.h:125
setup_threads_for_tof_biparallel
bool setup_threads_for_tof_biparallel()
Definition: library_daq.h:738
tex_times_of_flight
texture< time_of_flight_t, 1, cudaReadModeElementType > tex_times_of_flight
Definition: library_daq.h:143
constant_n_direction_bins_phi
__constant__ unsigned int constant_n_direction_bins_phi
Definition: library_daq.h:90
device_dz_per_time_bin
float * device_dz_per_time_bin
Definition: library_daq.h:112
setup_threads_nhits
bool setup_threads_nhits()
Definition: library_daq.h:822
setup_threads_for_tof
bool setup_threads_for_tof()
Definition: library_daq.h:715
kernel_find_vertex_with_max_npmts_in_timebin_and_directionbin
__global__ void kernel_find_vertex_with_max_npmts_in_timebin_and_directionbin(unsigned int *np, histogram_t *mnp, unsigned int *vmnp)
Definition: library_daq.h:1944
tex_PMTz
texture< float, 1, cudaReadModeElementType > tex_PMTz
Definition: library_daq.h:150
vertex_z
double * vertex_z
Definition: library_daq.h:75
constant_n_direction_bins_theta
__constant__ unsigned int constant_n_direction_bins_theta
Definition: library_daq.h:88
host_max_number_of_pmts_in_cone_in_time_bin
unsigned int * host_max_number_of_pmts_in_cone_in_time_bin
Definition: library_daq.h:182
coalesce_triggers
void coalesce_triggers()
Definition: library_daq.h:1403
device_number_of_pmts_in_cone_in_time_bin
unsigned int * device_number_of_pmts_in_cone_in_time_bin
Definition: library_daq.h:181
dark_rate
double dark_rate
Definition: library_daq.h:55
setup_threads_for_histo_iterated
void setup_threads_for_histo_iterated(bool last)
Definition: library_daq.h:2590
elapsed_tofs_copy_dev
float elapsed_tofs_copy_dev
Definition: library_daq.h:204
check_cudamalloc_unsigned_short
void check_cudamalloc_unsigned_short(unsigned int size)
Definition: library_daq.h:2497
elapsed_detector
float elapsed_detector
Definition: library_daq.h:204
constant_time_step_size
__constant__ unsigned int constant_time_step_size
Definition: library_daq.h:47
total_start
cudaEvent_t total_start
Definition: library_daq.h:171
constant_n_PMTs
__constant__ unsigned int constant_n_PMTs
Definition: library_daq.h:68
device_light_dx
float * device_light_dx
Definition: library_daq.h:127
candidate_trigger_meany_per_time_bin
std::vector< float > candidate_trigger_meany_per_time_bin
Definition: library_daq.h:162
fill_tofs_memory_on_device_nhits
void fill_tofs_memory_on_device_nhits()
Definition: library_daq.h:1323
tex_light_dz
texture< float, 1, cudaReadModeElementType > tex_light_dz
Definition: library_daq.h:146
set_input_file_for_event
bool set_input_file_for_event(int n)
Definition: library_daq.h:2204
write_output_nhits
void write_output_nhits(unsigned int n)
read_the_detector
bool read_the_detector()
Definition: library_daq.h:461
elapsed_parameters
float elapsed_parameters
Definition: library_daq.h:204
free_event_memories_nhits
void free_event_memories_nhits()
Definition: library_daq.h:2046
stop_cuda_clock
double stop_cuda_clock()
Definition: library_daq.h:1659
setup_threads_for_histo_per
void setup_threads_for_histo_per(unsigned int n)
Definition: library_daq.h:2601
number_of_kernel_blocks
unsigned int number_of_kernel_blocks
threads
Definition: library_daq.h:77
n_PMTs
unsigned int n_PMTs
pmts
Definition: library_daq.h:67
cerenkov_angle_water
double cerenkov_angle_water
Definition: library_daq.h:120
trigger_meanz_per_time_bin
std::vector< float > trigger_meanz_per_time_bin
Definition: library_daq.h:166
choose_candidates_above_threshold
void choose_candidates_above_threshold()
Definition: library_daq.h:2149
device_PMTy
float * device_PMTy
Definition: library_daq.h:137
elapsed_memory_tofs_dev
float elapsed_memory_tofs_dev
Definition: library_daq.h:204
read_value_from_file
float read_value_from_file(std::string paramname, std::string filename)
Definition: library_daq.h:2360
device_meany_for_vertex_with_max_n_pmts_per_time_bin
float * device_meany_for_vertex_with_max_n_pmts_per_time_bin
Definition: library_daq.h:185
print_parameters_2d
void print_parameters_2d()
Definition: library_daq.h:389
elapsed_copy_dev
float elapsed_copy_dev
Definition: library_daq.h:204
device_meanx_per_time_bin
float * device_meanx_per_time_bin
Definition: library_daq.h:113
elapsed_memory_candidates_host
float elapsed_memory_candidates_host
Definition: library_daq.h:204
correct_mode
unsigned int correct_mode
Definition: library_daq.h:174
print_times_of_flight
void print_times_of_flight()
Definition: library_daq.h:414
host_light_dx
float * host_light_dx
Definition: library_daq.h:128
free_event_memories
void free_event_memories()
Definition: library_daq.h:1983
elapsed_reset
float elapsed_reset
Definition: library_daq.h:204
device_max_number_of_pmts_in_cone_in_time_bin
unsigned int * device_max_number_of_pmts_in_cone_in_time_bin
Definition: library_daq.h:183
get_distance_index
unsigned int get_distance_index(unsigned int pmt_id, unsigned int vertex_block)
Definition: library_daq.h:1678
device_light_dr
float * device_light_dr
Definition: library_daq.h:133
vertex_y
double * vertex_y
Definition: library_daq.h:75
speed_light_water
double speed_light_water
Definition: library_daq.h:119
t0
struct timeval t0
Definition: library_daq.h:168
number_of_threads_per_block
unsigned int number_of_threads_per_block
Definition: library_daq.h:79
elapsed_directions
float elapsed_directions
Definition: library_daq.h:204
output_file
std::string output_file
Definition: library_daq.h:200
n_time_bins
unsigned int n_time_bins
Definition: library_daq.h:85
output_txt
bool output_txt
Definition: library_daq.h:173
stop_c_clock
double stop_c_clock()
Definition: library_daq.h:1651
elapsed_candidates_memory_dev
float elapsed_candidates_memory_dev
Definition: library_daq.h:204
t1
struct timeval t1
Definition: library_daq.h:169
setup_threads_to_find_candidates
bool setup_threads_to_find_candidates()
Definition: library_daq.h:809
make_table_of_directions
void make_table_of_directions()
Definition: library_daq.h:1209
host_n_pmts_nhits
unsigned int * host_n_pmts_nhits
Definition: library_daq.h:108
max_n_threads_per_block
int max_n_threads_per_block
Definition: library_daq.h:191
tex_times
texture< unsigned int, 1, cudaReadModeElementType > tex_times
Definition: library_daq.h:100
select_based_on_cone
bool select_based_on_cone
Definition: library_daq.h:58
twopi
double twopi
Definition: library_daq.h:123
timedifference_msec
float timedifference_msec(struct timeval t0, struct timeval t1)
Definition: library_daq.h:1641
read_number_of_input_hits
unsigned int read_number_of_input_hits()
Definition: library_daq.h:303
device_max_number_of_pmts_in_time_bin
histogram_t * device_max_number_of_pmts_in_time_bin
Definition: library_daq.h:178
constant_n_water_like_test_vertices
__constant__ unsigned int constant_n_water_like_test_vertices
Definition: library_daq.h:74
device_light_dz
float * device_light_dz
Definition: library_daq.h:131
time_step_size
unsigned int time_step_size
Definition: library_daq.h:46
trigger_meany_per_time_bin
std::vector< float > trigger_meany_per_time_bin
Definition: library_daq.h:165
allocate_correct_memory_on_device
void allocate_correct_memory_on_device()
Definition: library_daq.h:973
tex_light_dy
texture< float, 1, cudaReadModeElementType > tex_light_dy
Definition: library_daq.h:145
host_n_pmts_per_time_bin
unsigned int * host_n_pmts_per_time_bin
Definition: library_daq.h:106
elapsed_memory_directions_dev
float elapsed_memory_directions_dev
Definition: library_daq.h:204
max_n_hits_per_job
unsigned int max_n_hits_per_job
Definition: library_daq.h:54
setup_threads_for_tof_2d
bool setup_threads_for_tof_2d(unsigned int A, unsigned int B)
Definition: library_daq.h:761
read_input
bool read_input()
Definition: library_daq.h:323
read_the_pmts
bool read_the_pmts()
Definition: library_daq.h:446
host_vertex_with_max_n_pmts
unsigned int * host_vertex_with_max_n_pmts
Definition: library_daq.h:179
trigger_pair_vertex_time
std::vector< std::pair< unsigned int, unsigned int > > trigger_pair_vertex_time
Definition: library_daq.h:156
allocate_correct_memory_on_device_nhits
void allocate_correct_memory_on_device_nhits()
Definition: library_daq.h:1081
separate_triggers_into_gates
void separate_triggers_into_gates()
Definition: library_daq.h:1514
candidate_trigger_meanz_per_time_bin
std::vector< float > candidate_trigger_meanz_per_time_bin
Definition: library_daq.h:163
device_light_dy
float * device_light_dy
Definition: library_daq.h:129
elapsed_gates
float elapsed_gates
Definition: library_daq.h:204
host_light_dz
float * host_light_dz
Definition: library_daq.h:132
device_n_pmts_per_time_bin_and_direction_bin
unsigned int * device_n_pmts_per_time_bin_and_direction_bin
Definition: library_daq.h:109
read_the_input_ToolDAQ
bool read_the_input_ToolDAQ(std::vector< int > PMTid, std::vector< int > time)
trigger_npmts_in_time_bin
std::vector< unsigned int > trigger_npmts_in_time_bin
Definition: library_daq.h:157
candidate_trigger_meanx_per_time_bin
std::vector< float > candidate_trigger_meanx_per_time_bin
Definition: library_daq.h:161
elapsed_pmts
float elapsed_pmts
Definition: library_daq.h:204
host_meany_for_vertex_with_max_n_pmts_per_time_bin
float * host_meany_for_vertex_with_max_n_pmts_per_time_bin
Definition: library_daq.h:188
print_input
void print_input()
Definition: library_daq.h:400
cylindrical_grid
bool cylindrical_grid
Definition: library_daq.h:124
initialize_output_nhits
void initialize_output_nhits()
Definition: library_daq.h:2351
get_direction_index_at_time
unsigned int get_direction_index_at_time(unsigned int time_bin, unsigned int vertex_index, unsigned int direction_index)
Definition: library_daq.h:1711
costheta_cone_cut
float costheta_cone_cut
Definition: library_daq.h:56
constant_cerenkov_costheta
__constant__ float constant_cerenkov_costheta
Definition: library_daq.h:122
number_of_threads_per_block_3d
dim3 number_of_threads_per_block_3d
Definition: library_daq.h:80
device_n_pmts_per_time_bin
histogram_t * device_n_pmts_per_time_bin
Definition: library_daq.h:105
host_meanx_for_vertex_with_max_n_pmts_per_time_bin
float * host_meanx_for_vertex_with_max_n_pmts_per_time_bin
Definition: library_daq.h:187
offset_t
unsigned short offset_t
Definition: library_daq.h:19
choose_candidates
float choose_candidates
Definition: library_daq.h:204
detector_radius
double detector_radius
Definition: library_daq.h:65
water_like_threshold_number_of_pmts
unsigned int water_like_threshold_number_of_pmts
Definition: library_daq.h:48
candidate_trigger_pair_vertex_time
std::vector< std::pair< unsigned int, unsigned int > > candidate_trigger_pair_vertex_time
Definition: library_daq.h:153
nhits_threshold_min
unsigned int nhits_threshold_min
Definition: library_daq.h:50
read_the_input
bool read_the_input()
Definition: library_daq.h:838
constant_return_direction
__constant__ bool constant_return_direction
Definition: library_daq.h:62
tex_ids
texture< unsigned int, 1, cudaReadModeElementType > tex_ids
Definition: library_daq.h:97
event_file_base
std::string event_file_base
Definition: library_daq.h:201
fill_tofs_memory_on_device
void fill_tofs_memory_on_device()
Definition: library_daq.h:1257
return_vertex
bool return_vertex
Definition: library_daq.h:60
trigger_gate_down
double trigger_gate_down
Definition: library_daq.h:53
get_time_index
unsigned int get_time_index(unsigned int hit_index, unsigned int vertex_block)
Definition: library_daq.h:1685
n_direction_bins_theta
unsigned int n_direction_bins_theta
Definition: library_daq.h:87
constant_n_time_bins
__constant__ unsigned int constant_n_time_bins
Definition: library_daq.h:86
start_c_clock
void start_c_clock()
Definition: library_daq.h:1647
device_ids
unsigned int * device_ids
Definition: library_daq.h:96
greatest_divisor
unsigned int greatest_divisor
Definition: library_daq.h:210
trigger_gate_up
double trigger_gate_up
Definition: library_daq.h:52
device_time_bin_of_hit
unsigned int * device_time_bin_of_hit
Definition: library_daq.h:103
event_file
std::string event_file
Definition: library_daq.h:197
device_get_direction_index_at_angles
__device__ unsigned int device_get_direction_index_at_angles(unsigned int iphi, unsigned int itheta)
Definition: library_daq.h:1746
check_cudamalloc_float
void check_cudamalloc_float(unsigned int size)
Definition: library_daq.h:2460
detector_height
double detector_height
detector
Definition: library_daq.h:64
n_hits
unsigned int n_hits
Definition: library_daq.h:93
copy_candidates_from_device_to_host
void copy_candidates_from_device_to_host()
Definition: library_daq.h:2114
return_direction
bool return_direction
Definition: library_daq.h:61
number_of_kernel_blocks_3d
dim3 number_of_kernel_blocks_3d
Definition: library_daq.h:78
host_times_of_flight
time_of_flight_t * host_times_of_flight
Definition: library_daq.h:126
output_trigger_information
std::vector< double > output_trigger_information
Definition: library_daq.h:160
elapsed_tofs_free
float elapsed_tofs_free
Definition: library_daq.h:204
stop
cudaEvent_t stop
Definition: library_daq.h:171
host_ids
unsigned int * host_ids
Definition: library_daq.h:95
host_max_number_of_pmts_in_time_bin
histogram_t * host_max_number_of_pmts_in_time_bin
Definition: library_daq.h:177
device_times
unsigned int * device_times
Definition: library_daq.h:99
start_total_cuda_clock
void start_total_cuda_clock()
Definition: library_daq.h:1666
host_light_dy
float * host_light_dy
Definition: library_daq.h:130
constant_n_direction_bins
__constant__ unsigned int constant_n_direction_bins
Definition: library_daq.h:92
event_file_suffix
std::string event_file_suffix
Definition: library_daq.h:202
fill_correct_memory_on_device
void fill_correct_memory_on_device()
Definition: library_daq.h:1334
use_timing
bool use_timing
Definition: library_daq.h:195
get_direction_index_at_pmt
unsigned int get_direction_index_at_pmt(unsigned int pmt_id, unsigned int vertex_index, unsigned int direction_index)
Definition: library_daq.h:1701
output_file_base
std::string output_file_base
Definition: library_daq.h:203
start
cudaEvent_t start
Definition: library_daq.h:171
check_cudamalloc_unsigned_char
void check_cudamalloc_unsigned_char(unsigned int size)
Definition: library_daq.h:2509
host_time_bin_of_hit
unsigned int * host_time_bin_of_hit
Definition: library_daq.h:102
time_offset
offset_t time_offset
hits
Definition: library_daq.h:83
print_pmts
void print_pmts()
Definition: library_daq.h:407
elapsed_input
float elapsed_input
Definition: library_daq.h:204
elapsed_vertices
float elapsed_vertices
Definition: library_daq.h:204
allocate_tofs_memory_on_device
void allocate_tofs_memory_on_device()
Definition: library_daq.h:910
final_trigger_pair_vertex_time
std::vector< std::pair< unsigned int, unsigned int > > final_trigger_pair_vertex_time
Definition: library_daq.h:159
start_cuda_clock
void start_cuda_clock()
Definition: library_daq.h:1655
cerenkov_costheta
float cerenkov_costheta
Definition: library_daq.h:121
read_user_parameters
void read_user_parameters()
Definition: library_daq.h:2386
PMT_z
double * PMT_z
Definition: library_daq.h:69
read_detector
bool read_detector()
Definition: library_daq.h:361
device_n_pmts_nhits
unsigned int * device_n_pmts_nhits
Definition: library_daq.h:107
host_PMTy
float * host_PMTy
Definition: library_daq.h:138
device_get_time_index
__device__ unsigned int device_get_time_index(unsigned int hit_index, unsigned int vertex_block)
Definition: library_daq.h:1729
elapsed_coalesce
float elapsed_coalesce
Definition: library_daq.h:204
device_vertex_with_max_n_pmts
unsigned int * device_vertex_with_max_n_pmts
Definition: library_daq.h:180
host_times
unsigned int * host_times
Definition: library_daq.h:98
host_PMTz
float * host_PMTz
Definition: library_daq.h:140
kernel_find_vertex_with_max_npmts_and_center_of_mass_in_timebin
__global__ void kernel_find_vertex_with_max_npmts_and_center_of_mass_in_timebin(histogram_t *np, histogram_t *mnp, unsigned int *vmnp, unsigned int *nc, unsigned int *mnc)
Definition: library_daq.h:1895
n_direction_bins
unsigned int n_direction_bins
Definition: library_daq.h:91
constant_costheta_cone_cut
__constant__ float constant_costheta_cone_cut
Definition: library_daq.h:57
check_cudamalloc_bool
void check_cudamalloc_bool(unsigned int size)
Definition: library_daq.h:2521
elapsed_threads
float elapsed_threads
Definition: library_daq.h:204
check_cudamalloc_unsigned_int
void check_cudamalloc_unsigned_int(unsigned int size)
Definition: library_daq.h:2484
make_test_vertices
void make_test_vertices()
Definition: library_daq.h:470
make_table_of_tofs
void make_table_of_tofs()
Definition: library_daq.h:1165
tex_light_dx
texture< float, 1, cudaReadModeElementType > tex_light_dx
Definition: library_daq.h:144
find_greatest_divisor
unsigned int find_greatest_divisor(unsigned int n, unsigned int max)
Definition: library_daq.h:2533
host_meanz_for_vertex_with_max_n_pmts_per_time_bin
float * host_meanz_for_vertex_with_max_n_pmts_per_time_bin
Definition: library_daq.h:189
PMT_y
double * PMT_y
Definition: library_daq.h:69
device_meanz_for_vertex_with_max_n_pmts_per_time_bin
float * device_meanz_for_vertex_with_max_n_pmts_per_time_bin
Definition: library_daq.h:186
read_pmts
bool read_pmts()
Definition: library_daq.h:1378
n_direction_bins_phi
unsigned int n_direction_bins_phi
Definition: library_daq.h:89
trigger_meanx_per_time_bin
std::vector< float > trigger_meanx_per_time_bin
Definition: library_daq.h:164
write_output_mode
unsigned int write_output_mode
Definition: library_daq.h:175
host_PMTx
float * host_PMTx
Definition: library_daq.h:136
initialize_output
void initialize_output()
Definition: library_daq.h:2344
trigger_npmts_in_cone_in_time_bin
std::vector< unsigned int > trigger_npmts_in_cone_in_time_bin
Definition: library_daq.h:158
elapsed_memory_dev
float elapsed_memory_dev
Definition: library_daq.h:204
vertex_x
double * vertex_x
Definition: library_daq.h:75
elapsed_kernel_correct_times_and_get_n_pmts_per_time_bin
float elapsed_kernel_correct_times_and_get_n_pmts_per_time_bin
Definition: library_daq.h:204
device_PMTx
float * device_PMTx
Definition: library_daq.h:135
print_directions
void print_directions()
Definition: library_daq.h:429
device_meanz_per_time_bin
float * device_meanz_per_time_bin
Definition: library_daq.h:115
grid_size
unsigned int grid_size
Definition: library_daq.h:81
fill_directions_memory_on_device
void fill_directions_memory_on_device()
Definition: library_daq.h:1303
device_meany_per_time_bin
float * device_meany_per_time_bin
Definition: library_daq.h:114
use_verbose
bool use_verbose
Definition: library_daq.h:194
host_light_dr
float * host_light_dr
Definition: library_daq.h:134
get_direction_index_at_angles
unsigned int get_direction_index_at_angles(unsigned int iphi, unsigned int itheta)
Definition: library_daq.h:1692
wall_like_distance
double wall_like_distance
Definition: library_daq.h:45
pmts_file
std::string pmts_file
Definition: library_daq.h:199
detector_file
std::string detector_file
Definition: library_daq.h:198
time_of_flight_t
unsigned short time_of_flight_t
Definition: library_daq.h:35
constant_select_based_on_cone
__constant__ bool constant_select_based_on_cone
Definition: library_daq.h:59
n_water_like_test_vertices
unsigned int n_water_like_test_vertices
Definition: library_daq.h:72
device_get_direction_index_at_pmt
__device__ unsigned int device_get_direction_index_at_pmt(unsigned int pmt_id, unsigned int vertex_index, unsigned int direction_index)
Definition: library_daq.h:1736
host_directions_for_vertex_and_pmt
bool * host_directions_for_vertex_and_pmt
Definition: library_daq.h:142
device_get_distance_index
__device__ unsigned int device_get_distance_index(unsigned int pmt_id, unsigned int vertex_block)
Definition: library_daq.h:1722
coalesce_time
double coalesce_time
Definition: library_daq.h:51
total_stop
cudaEvent_t total_stop
Definition: library_daq.h:171
max_n_blocks
int max_n_blocks
Definition: library_daq.h:192
device_dx_per_time_bin
float * device_dx_per_time_bin
Definition: library_daq.h:110
elapsed_directions_copy_dev
float elapsed_directions_copy_dev
Definition: library_daq.h:204
tex_PMTy
texture< float, 1, cudaReadModeElementType > tex_PMTy
Definition: library_daq.h:149
n_events
int n_events
Definition: library_daq.h:213
elapsed_tof
float elapsed_tof
Definition: library_daq.h:204
read_number_of_pmts
unsigned int read_number_of_pmts()
Definition: library_daq.h:1358
constant_n_test_vertices
__constant__ unsigned int constant_n_test_vertices
Definition: library_daq.h:73
kernel_find_vertex_with_max_npmts_in_timebin
__global__ void kernel_find_vertex_with_max_npmts_in_timebin(histogram_t *np, histogram_t *mnp, unsigned int *vmnp)
Definition: library_daq.h:1818
elapsed_total
float elapsed_total
Definition: library_daq.h:204
constant_time_offset
__constant__ offset_t constant_time_offset
Definition: library_daq.h:84
tex_light_dr
texture< float, 1, cudaReadModeElementType > tex_light_dr
Definition: library_daq.h:147
print_gpu_properties
void print_gpu_properties()
Definition: library_daq.h:1766
elapsed_write_output
float elapsed_write_output
Definition: library_daq.h:204
print_parameters
void print_parameters()
Definition: library_daq.h:378
elapsed_free
float elapsed_free
Definition: library_daq.h:204
write_output
void write_output()
Definition: library_daq.h:2266
check_cudamalloc_int
void check_cudamalloc_int(unsigned int size)
Definition: library_daq.h:2472
candidate_trigger_npmts_in_cone_in_time_bin
std::vector< unsigned int > candidate_trigger_npmts_in_cone_in_time_bin
Definition: library_daq.h:155
elapsed_candidates_kernel
float elapsed_candidates_kernel
Definition: library_daq.h:204
elapsed_threads_candidates
float elapsed_threads_candidates
Definition: library_daq.h:204
Generated by   doxygen 1.8.5