source: branches/Corsika7500Compatibility/fact/processing/threshold.C@ 20115

Last change on this file since 20115 was 18309, checked in by Daniela Dorner, 9 years ago
macros to read the thresholds from the auxiliary files
File size: 1.6 KB
Line 
1#include <algorithm>
2
3void threshold(const char *fname, double beg=0, double end=100000)
4{
5 fits file(fname);
6
7 //file.PrintColumns();
8 //file.PrintKeys();
9
10 Double_t time;
11 UShort_t th[160];
12 file.SetPtrAddress("Time", &time);
13 file.SetPtrAddress("PatchThresh", th);
14
15 UInt_t offset = file.GetUInt("MJDREF");
16 if (beg < 30000)
17 beg+=offset;
18 if (end < 30000)
19 end+=offset;
20
21 double average = 0;
22 double med_avg = 0;
23 double med_rms = 0;
24 UShort_t max_tot = 0;
25 int cnt = 0;
26
27 double avg_last = -1;
28 double med_last = -1;
29 UShort_t max_last = -1;
30 double diff = -1;
31
32 while (file.GetNextRow())
33 {
34 time += offset;
35
36 if (time>end)
37 break;
38
39 sort(th, th+160);
40
41 Float_t med = (th[80] + th[81]) / 2;
42 UShort_t max = th[159];
43
44 Double_t avg = 0;
45 for (int i=0; i<160; i++)
46 avg += th[i];
47 avg /= 160;
48
49 if (time<beg)
50 {
51 avg_last = avg;
52 med_last = med;
53 max_last = max;
54 diff = beg-time;
55 continue;
56 }
57
58 med_avg += med;
59 med_rms += med*med;
60 average += avg;
61
62 if (max>max_tot)
63 max_tot = max;
64
65 cnt ++;
66 }
67
68 if (cnt==0)
69 {
70 if (diff<5./24/3600)
71 return;
72
73 cout << "result " << med_last << " 0 " << max_last << " " << avg_last << endl;
74 return;
75 }
76
77 average /= cnt;
78 med_avg /= cnt;
79 med_rms /= cnt;
80 med_rms = sqrt(med_rms-med_avg*med_avg);
81
82 cout << "result " << med_avg << " " << med_rms << " " << max_tot << " " << average << endl;
83}
Note: See TracBrowser for help on using the repository browser.