1 | //********************************
|
---|
2 | //
|
---|
3 | // Class CalFactFits
|
---|
4 | // Wrapper class for fits.h or pyfits.h
|
---|
5 | // Provides fast access to calibrated events of FACT raw files
|
---|
6 | //
|
---|
7 | // written by Thomas Kraehenbuehl, ETH Zurich
|
---|
8 | // tpk@phys.ethz.ch, +41 44 633 3973
|
---|
9 | // April 2012
|
---|
10 | //
|
---|
11 | //********************************
|
---|
12 | //
|
---|
13 | // Compilation (root in pyfact directory)
|
---|
14 | // root [0] gSystem->Load("/usr/lib64/libz.so");
|
---|
15 | // root [1] .L izstream.h++
|
---|
16 | // root [2] .L factfits.h++
|
---|
17 | // root [3] .L calfactfits.h++
|
---|
18 | //
|
---|
19 | // Usage in Python:
|
---|
20 | // See pyscripts/examples/CalFitsTest.py
|
---|
21 | //
|
---|
22 | //********************************
|
---|
23 |
|
---|
24 | //ToDo: shared library creation debuggen
|
---|
25 |
|
---|
26 | #ifndef CALFACTFITS_H
|
---|
27 | #define CALFACTFITS_H
|
---|
28 |
|
---|
29 | #include <cstdio>
|
---|
30 | #include <string>
|
---|
31 |
|
---|
32 | #ifndef __MARS__
|
---|
33 | #include <vector>
|
---|
34 | #include <iomanip>
|
---|
35 | #include <iostream>
|
---|
36 | #define gLog cerr
|
---|
37 | #define ___err___ ""
|
---|
38 | #define ___all___ ""
|
---|
39 | #else
|
---|
40 | #include "MLog.h"
|
---|
41 | #include "MLogManip.h"
|
---|
42 | #define ___err___ err
|
---|
43 | #define ___all___ all
|
---|
44 | #endif
|
---|
45 |
|
---|
46 | #ifdef __EXCEPTIONS
|
---|
47 | #include <stdexcept>
|
---|
48 | #endif
|
---|
49 |
|
---|
50 | #include "factfits.h"
|
---|
51 |
|
---|
52 | class CalFactFits
|
---|
53 | {
|
---|
54 | public:
|
---|
55 | //No standard constructor CalFactFits()!
|
---|
56 |
|
---|
57 | //Direct handlers of the files
|
---|
58 | FactFits datafile, calibfile; //Class name should be PyFits or better FactPyFits...
|
---|
59 |
|
---|
60 | //Basic file parameters
|
---|
61 | UInt_t calib_nroi, calib_npix;
|
---|
62 | UInt_t calib_blm_size, calib_gm_size, calib_tom_size;
|
---|
63 | UInt_t data_nroi, data_npix, data_ndata;
|
---|
64 |
|
---|
65 | //Common variables
|
---|
66 | UInt_t nroi, npix, nevents;
|
---|
67 |
|
---|
68 | //Calibration variables
|
---|
69 | float* calib_baselinemean;
|
---|
70 | float* calib_gainmean;
|
---|
71 | float* calib_triggeroffsetmean;
|
---|
72 | // The same variables as above, but this time
|
---|
73 | // they are doubled in memory, so the modulo operation is not necessary.
|
---|
74 | // idea by TPK again
|
---|
75 | // since we like to work with doubles for some reason, I don't know yet
|
---|
76 | // I cast the float to double already in this stage.
|
---|
77 |
|
---|
78 | double* baseline;
|
---|
79 | double* gain;
|
---|
80 | double* trigger_offset;
|
---|
81 | //Using <vector> instead of arrays makes no visible difference
|
---|
82 | //ToDo: use arrays of size 1440x1024 (x2 for wrap-arounds) and read all variables into those
|
---|
83 |
|
---|
84 | //Event variables
|
---|
85 | UInt_t event_id;
|
---|
86 | UShort_t event_triggertype;
|
---|
87 | int16_t* event_data;
|
---|
88 | int16_t* event_offset;
|
---|
89 | int32_t* event_boardtimes;
|
---|
90 | double* npcaldata;
|
---|
91 |
|
---|
92 | CalFactFits(const string &datafilename, const string &calibfilename) //Constructor with two filenames
|
---|
93 | : datafile(datafilename),
|
---|
94 | calibfile(calibfilename),
|
---|
95 | npcaldata(NULL)
|
---|
96 | {
|
---|
97 | //cout << "Constructor called" << endl;
|
---|
98 | //Read basic parameters of the two files
|
---|
99 | // std::cout << "...Reading basic file parameters..." << std::endl;
|
---|
100 | calib_nroi = calibfile.GetUInt("NROI");
|
---|
101 | calib_npix = calibfile.GetUInt("NPIX");
|
---|
102 | data_nroi = datafile.GetUInt("NROI");
|
---|
103 | data_npix = datafile.GetUInt("NPIX");
|
---|
104 | data_ndata = datafile.GetN("Data");
|
---|
105 |
|
---|
106 | calib_blm_size = calibfile.GetN("BaselineMean")/calib_npix;
|
---|
107 | calib_gm_size = calibfile.GetN("GainMean")/calib_npix;
|
---|
108 | calib_tom_size = calibfile.GetN("TriggerOffsetMean")/calib_npix;
|
---|
109 |
|
---|
110 | // std::cout << "Column sizes: " << calib_blm_size << " " << calib_gm_size << " " << calib_tom_size << std::endl;
|
---|
111 |
|
---|
112 | //Define the common variables
|
---|
113 | if((calib_nroi==data_nroi)&&(calib_npix==data_npix)&&(data_nroi*data_npix==data_ndata)&&(calib_blm_size==calib_gm_size)&&(calib_tom_size==calib_nroi)) {
|
---|
114 | nroi = data_nroi;
|
---|
115 | npix = data_npix;
|
---|
116 | }
|
---|
117 | else {
|
---|
118 | ostringstream str;
|
---|
119 | str << "Data/calib file error: NROI mismatch, NPIX mismatch, data column size wrong or calib columns mismatch.";
|
---|
120 | #ifdef __EXCEPTIONS
|
---|
121 | throw runtime_error(str.str());
|
---|
122 | #else
|
---|
123 | gLog << ___err___ << "ERROR - " << str.str() << endl;
|
---|
124 | return;
|
---|
125 | #endif
|
---|
126 | }
|
---|
127 | nevents = datafile.GetNumRows();
|
---|
128 |
|
---|
129 | //Read the calibration data
|
---|
130 | // std::cout << "...Reading calibration data..." << std::endl;
|
---|
131 | calib_baselinemean = new float[calibfile.GetN("BaselineMean")];
|
---|
132 | calibfile.SetPtrAddress("BaselineMean", calib_baselinemean, calibfile.GetN("BaselineMean"));
|
---|
133 | baseline = new double[calibfile.GetN("BaselineMean")*2];
|
---|
134 |
|
---|
135 | calib_gainmean = new float[calibfile.GetN("GainMean")];
|
---|
136 | calibfile.SetPtrAddress("GainMean", calib_gainmean, calibfile.GetN("GainMean"));
|
---|
137 | gain = new double[calibfile.GetN("GainMean")*2];
|
---|
138 |
|
---|
139 | calib_triggeroffsetmean = new float[calibfile.GetN("TriggerOffsetMean")];
|
---|
140 | calibfile.SetPtrAddress("TriggerOffsetMean", calib_triggeroffsetmean, calibfile.GetN("TriggerOffsetMean"));
|
---|
141 | trigger_offset = new double[calibfile.GetN("TriggerOffsetMean")];
|
---|
142 |
|
---|
143 | calibfile.GetRow(0);
|
---|
144 |
|
---|
145 | for (int i = 0; i < (int)calibfile.GetN("BaselineMean"); ++i)
|
---|
146 | {
|
---|
147 | baseline[i] = double(calib_baselinemean[i]) *4096./2000.;
|
---|
148 | baseline[i+(int)calibfile.GetN("BaselineMean")] = double(calib_baselinemean[i]) *4096./2000.;
|
---|
149 | }
|
---|
150 |
|
---|
151 | for (int i = 0; i < (int)calibfile.GetN("GainMean"); ++i)
|
---|
152 | {
|
---|
153 | gain[i] = double(calib_gainmean[i]) * 1907.35;
|
---|
154 | gain[i+(int)calibfile.GetN("GainMean")] = double(calib_gainmean[i]) * 1907.35 *4096. /2000.;
|
---|
155 | }
|
---|
156 |
|
---|
157 | for (int i = 0; i < (int)calibfile.GetN("TriggerOffsetMean"); ++i)
|
---|
158 | {
|
---|
159 | trigger_offset[i] = double(calib_triggeroffsetmean[i])*4096./2000.;
|
---|
160 | trigger_offset[i+(int)calibfile.GetN("TriggerOffsetMean")] = double(calib_triggeroffsetmean[i])*4096./2000.;
|
---|
161 | }
|
---|
162 |
|
---|
163 |
|
---|
164 |
|
---|
165 | //Set the event pointers
|
---|
166 | // std::cout << "...Setting event pointers..." << std::endl;
|
---|
167 | datafile.SetRefAddress("EventNum", event_id);
|
---|
168 | datafile.SetRefAddress("TriggerType", event_triggertype);
|
---|
169 |
|
---|
170 | event_data = new int16_t[data_ndata];
|
---|
171 | datafile.SetPtrAddress("Data", event_data, data_ndata);
|
---|
172 |
|
---|
173 | event_offset = new int16_t[datafile.GetN("StartCellData")];
|
---|
174 | datafile.SetPtrAddress("StartCellData", event_offset, datafile.GetN("StartCellData"));
|
---|
175 |
|
---|
176 | event_boardtimes = new int32_t[datafile.GetN("BoardTime")];
|
---|
177 | datafile.SetPtrAddress("BoardTime", event_boardtimes, datafile.GetN("BoardTime"));
|
---|
178 | }
|
---|
179 |
|
---|
180 | ~CalFactFits() //Standard destructor
|
---|
181 | {
|
---|
182 | //cout << "Destructor called " << endl;
|
---|
183 | //cout << "calib_baselinemean " << calib_baselinemean << endl;
|
---|
184 | //cout << "calib_gainmean " << calib_gainmean << endl;
|
---|
185 | //cout << "calib_triggeroffsetmean " << calib_triggeroffsetmean << endl;
|
---|
186 | //cout << "event_data " << event_data << endl;
|
---|
187 | //cout << "event_offset " << event_offset << endl;
|
---|
188 | //cout << "event_boardtimes " << event_boardtimes << endl;
|
---|
189 | //cout << "baseline " << baseline << endl;
|
---|
190 | //cout << "gain " << gain << endl;
|
---|
191 | //cout << "trigger_offset " << trigger_offset << endl;
|
---|
192 |
|
---|
193 | delete[] calib_baselinemean;
|
---|
194 | delete[] calib_gainmean;
|
---|
195 | //delete[] calib_triggeroffsetmean;
|
---|
196 | delete[] event_data;
|
---|
197 | delete[] event_offset;
|
---|
198 | delete[] event_boardtimes;
|
---|
199 | delete[] baseline;
|
---|
200 | //delete[] gain;
|
---|
201 | delete[] trigger_offset;
|
---|
202 | }
|
---|
203 |
|
---|
204 | bool GetCalEvent() //Read calibrated event into the event variables
|
---|
205 | {
|
---|
206 | if(!npcaldata) {
|
---|
207 | ostringstream str;
|
---|
208 | str << "Pointer to the calibrated data not initialized!";
|
---|
209 | #ifdef __EXCEPTIONS
|
---|
210 | throw runtime_error(str.str());
|
---|
211 | #else
|
---|
212 | gLog << ___err___ << "ERROR - " << str.str() << endl;
|
---|
213 | return false;
|
---|
214 | #endif
|
---|
215 | }
|
---|
216 | if(datafile.GetNextRow() == false) {
|
---|
217 | // std::cout << "Last event reached..." << std::endl;
|
---|
218 | return false;
|
---|
219 | }
|
---|
220 | else {
|
---|
221 | UInt_t drs_calib_offset;
|
---|
222 | double raw, raw_bsl, shifted;
|
---|
223 | for(UInt_t pixel = 0; pixel < data_npix; pixel++)
|
---|
224 | {
|
---|
225 | for(UInt_t slice = 0; slice < data_nroi; slice++)
|
---|
226 | {
|
---|
227 | drs_calib_offset = slice + event_offset[pixel];
|
---|
228 | raw = (double)event_data[pixel*data_nroi+slice];
|
---|
229 | raw_bsl = raw - baseline[pixel*calib_blm_size+drs_calib_offset];
|
---|
230 | shifted = raw_bsl - trigger_offset[pixel*data_nroi+slice];
|
---|
231 | npcaldata[pixel*data_nroi+slice] = shifted / gain[pixel*calib_blm_size+drs_calib_offset];
|
---|
232 | //npcaldata[pixel*data_nroi+slice] = double((event_data[pixel*data_nroi+slice]*2000./4096.-calib_baselinemean[pixel*calib_blm_size+drs_calib_offset]-calib_triggeroffsetmean[pixel*data_nroi+slice])/calib_gainmean[pixel*calib_blm_size+drs_calib_offset]*1907.35);
|
---|
233 | //Note: data_nroi=calib_nroi, calib_blm_size=calib_gm_size
|
---|
234 | }
|
---|
235 | }
|
---|
236 | }
|
---|
237 | return true;
|
---|
238 | }
|
---|
239 |
|
---|
240 | void SetNpcaldataPtr(double *numpyptr) //Set the pointer for the calibrated data to the numpy array
|
---|
241 | {
|
---|
242 | npcaldata = numpyptr;
|
---|
243 | return;
|
---|
244 | }
|
---|
245 | };
|
---|
246 | #endif /* CALFACTFITS_H */
|
---|