source: fact/tools/pyscripts/pyfact/calfactfits.h@ 17691

Last change on this file since 17691 was 17691, checked in by kraehenb, 10 years ago
The file extern_Mars_mcore/factfits.h was imported instead of factfits.h. Now corrected.
File size: 9.5 KB
Line 
1//********************************
2//
3// Class CalFactFits
4// Wrapper class for factfits.h (or fits.h, 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// Try /usr/lib/x86_64-linux-gnu/libz.so at ETH for the libz.so file.
20//
21// Note: the produced izstream_h.so must be in LD_LIBRARY_PATH, e.g.
22// export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/guest/kraehenb/software/PythonIncludes/
23//
24// Usage in Python:
25// See pyscripts/examples/CalFitsTest.py
26//
27//********************************
28
29//ToDo: shared library creation debuggen
30
31#ifndef CALFACTFITS_H
32#define CALFACTFITS_H
33
34#include <cstdio>
35#include <string>
36
37#ifndef __MARS__
38#include <vector>
39#include <iomanip>
40#include <iostream>
41#define gLog cerr
42#define ___err___ ""
43#define ___all___ ""
44#else
45#include "MLog.h"
46#include "MLogManip.h"
47#define ___err___ err
48#define ___all___ all
49#endif
50
51#ifdef __EXCEPTIONS
52#include <stdexcept>
53#endif
54
55#include "factfits.h"
56
57class CalFactFits
58{
59public:
60 //No standard constructor CalFactFits()!
61
62 //Direct handlers of the files
63 factfits datafile, calibfile; //Class name should be PyFits or better FactPyFits...
64
65 //Basic file parameters
66 UInt_t calib_nroi, calib_npix;
67 UInt_t calib_blm_size, calib_gm_size, calib_tom_size;
68 UInt_t data_nroi, data_npix, data_ndata;
69
70 //Common variables
71 UInt_t nroi, npix, nevents;
72
73 //Calibration variables
74 float* calib_baselinemean;
75 float* calib_gainmean;
76 float* calib_triggeroffsetmean;
77 // The same variables as above, but this time
78 // they are doubled in memory, so the modulo operation is not necessary.
79 // idea by TPK again
80 // since we like to work with doubles for some reason, I don't know yet
81 // I cast the float to double already in this stage.
82
83 double* baseline;
84 double* gain;
85 double* trigger_offset;
86 //Using <vector> instead of arrays makes no visible difference
87 //ToDo: use arrays of size 1440x1024 (x2 for wrap-arounds) and read all variables into those
88
89 //Event variables
90 UInt_t event_id;
91 UShort_t event_triggertype;
92 int16_t* event_data;
93 int16_t* event_offset;
94 int32_t* event_boardtimes;
95 double* npcaldata;
96
97 CalFactFits(const string &datafilename, const string &calibfilename) //Constructor with two filenames
98 : datafile(datafilename),
99 calibfile(calibfilename),
100 npcaldata(NULL)
101 {
102 //cout << "Constructor called" << endl;
103 //Read basic parameters of the two files
104// std::cout << "...Reading basic file parameters..." << std::endl;
105 try {
106 calib_nroi = calibfile.GetUInt("NROI");
107 }
108 catch(runtime_error &err) {
109 std::cout << err.what() << std::endl;
110 std::cout << "Standard value NROI=1024 assumed." << std::endl;
111 calib_nroi = 1024;
112 }
113 try {
114 calib_npix = calibfile.GetUInt("NPIX");
115 }
116 catch(runtime_error &err) {
117 std::cout << err.what() << std::endl;
118 std::cout << "Standard value NPIX=1440 assumed." << std::endl;
119 calib_npix = 1440;
120 }
121
122 data_nroi = datafile.GetUInt("NROI");
123 data_npix = datafile.GetUInt("NPIX");
124 data_ndata = datafile.GetN("Data");
125
126 calib_blm_size = calibfile.GetN("BaselineMean")/calib_npix;
127 calib_gm_size = calibfile.GetN("GainMean")/calib_npix;
128 calib_tom_size = calibfile.GetN("TriggerOffsetMean")/calib_npix;
129
130// std::cout << "Column sizes: " << calib_blm_size << " " << calib_gm_size << " " << calib_tom_size << std::endl;
131
132 //Define the common variables
133 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)) {
134 nroi = data_nroi;
135 npix = data_npix;
136 }
137 else {
138 ostringstream str;
139 str << "Data/calib file error: NROI mismatch, NPIX mismatch, data column size wrong or calib columns mismatch.";
140#ifdef __EXCEPTIONS
141 throw runtime_error(str.str());
142#else
143 gLog << ___err___ << "ERROR - " << str.str() << endl;
144 return;
145#endif
146 }
147 nevents = datafile.GetNumRows();
148
149 //Read the calibration data
150// std::cout << "...Reading calibration data..." << std::endl;
151 calib_baselinemean = new float[calibfile.GetN("BaselineMean")];
152 calibfile.SetPtrAddress("BaselineMean", calib_baselinemean, calibfile.GetN("BaselineMean"));
153 baseline = new double[calibfile.GetN("BaselineMean")*2];
154
155 calib_gainmean = new float[calibfile.GetN("GainMean")];
156 calibfile.SetPtrAddress("GainMean", calib_gainmean, calibfile.GetN("GainMean"));
157 gain = new double[calibfile.GetN("GainMean")*2];
158
159 calib_triggeroffsetmean = new float[calibfile.GetN("TriggerOffsetMean")];
160 calibfile.SetPtrAddress("TriggerOffsetMean", calib_triggeroffsetmean, calibfile.GetN("TriggerOffsetMean"));
161 trigger_offset = new double[calibfile.GetN("TriggerOffsetMean")*2];
162
163 calibfile.GetRow(0);
164
165int orig_index, new_index1, new_index2;
166 for (int pix=0 ; pix < (int)calib_npix; ++pix)
167 {
168 for (int sl = 0; sl < (int)calib_blm_size; ++sl)
169 {
170 orig_index = pix * (int)calib_blm_size + sl;
171 new_index1 = 2*pix * (int)calib_blm_size + sl;
172 new_index2 = (2*pix+1) * (int)calib_blm_size + sl;
173
174 baseline[new_index2] = baseline[new_index1] = double(calib_baselinemean[orig_index]) *4096./2000.;
175 }
176 }
177 for (int pix=0 ; pix < (int)calib_npix; ++pix)
178 {
179 for (int sl = 0; sl < (int)calib_gm_size; ++sl)
180 {
181 orig_index = pix * (int)calib_gm_size + sl;
182 new_index1 = 2*pix * (int)calib_gm_size + sl;
183 new_index2 = (2*pix+1) * (int)calib_gm_size + sl;
184
185 gain[new_index2] = gain[new_index1] = double(calib_gainmean[orig_index]) /( 1907.35 /4096. *2000.);
186 }
187 }
188
189 for (int pix=0 ; pix < (int)calib_npix; ++pix)
190 {
191 for (int sl = 0; sl < (int)calib_tom_size; ++sl)
192 {
193 orig_index = pix * (int)calib_tom_size + sl;
194 new_index1 = 2*pix * (int)calib_tom_size + sl;
195 new_index2 = (2*pix+1) * (int)calib_tom_size + sl;
196
197 trigger_offset[new_index2] = trigger_offset[new_index1] = double(calib_triggeroffsetmean[orig_index]) *4096. /2000.;
198 }
199 }
200
201
202
203 //Set the event pointers
204// std::cout << "...Setting event pointers..." << std::endl;
205 datafile.SetRefAddress("EventNum", event_id);
206 datafile.SetRefAddress("TriggerType", event_triggertype);
207
208 event_data = new int16_t[data_ndata];
209 datafile.SetPtrAddress("Data", event_data, data_ndata);
210
211 event_offset = new int16_t[datafile.GetN("StartCellData")];
212 datafile.SetPtrAddress("StartCellData", event_offset, datafile.GetN("StartCellData"));
213
214 event_boardtimes = new int32_t[datafile.GetN("BoardTime")];
215 datafile.SetPtrAddress("BoardTime", event_boardtimes, datafile.GetN("BoardTime"));
216 }
217
218 ~CalFactFits() //Standard destructor
219 {
220 delete[] calib_baselinemean;
221 delete[] calib_gainmean;
222 delete[] calib_triggeroffsetmean;
223 delete[] event_data;
224 delete[] event_offset;
225 delete[] event_boardtimes;
226 delete[] baseline;
227 delete[] gain;
228 delete[] trigger_offset;
229 }
230
231 bool GetCalEvent() //Read calibrated event into the event variables
232 {
233 if(!npcaldata) {
234 ostringstream str;
235 str << "Pointer to the calibrated data not initialized!";
236#ifdef __EXCEPTIONS
237 throw runtime_error(str.str());
238#else
239 gLog << ___err___ << "ERROR - " << str.str() << endl;
240 return false;
241#endif
242 }
243 if(datafile.GetNextRow() == false) {
244// std::cout << "Last event reached..." << std::endl;
245 return false;
246 }
247 else
248 {
249 UInt_t drs_calib_offset;
250 double raw, raw_bsl, shifted;
251 for(UInt_t pixel = 0; pixel < data_npix; pixel++)
252 {
253 for(UInt_t slice = 0; slice < data_nroi; slice++)
254 {
255 drs_calib_offset = slice + event_offset[pixel];
256 raw = (double)event_data[pixel*data_nroi+slice];
257 raw_bsl = raw - baseline[pixel*calib_blm_size*2+drs_calib_offset];
258 shifted = raw_bsl - trigger_offset[pixel*data_nroi*2+slice];
259 npcaldata[pixel*data_nroi+slice] = shifted / gain[pixel*calib_blm_size*2+drs_calib_offset];
260 //Note: data_nroi=calib_nroi, calib_blm_size=calib_gm_size
261 }
262
263 // stolen from TBs MARS ... DrsCalib.h or something like that
264 double *p = npcaldata+(pixel*data_nroi);
265
266 // second loop for despiking
267 for (size_t i=1; i<data_nroi-2; i++)
268 {
269 if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
270 {
271 //std::cout << "single cand:\t" << pixel << "\t" << i << std::endl;
272 p[i] = (p[i-1]+p[i+1])/2;
273 }
274 if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
275 {
276 //std::cout << "double cand:\t" << pixel << "\t" << i << std::endl;
277 p[i] = (p[i-1]+p[i+2])/2;
278 p[i+1] = p[i];
279 }
280 }
281
282 }
283 }
284 return true;
285 }
286
287 void SetNpcaldataPtr(double *numpyptr) //Set the pointer for the calibrated data to the numpy array
288 {
289 npcaldata = numpyptr;
290 return;
291 }
292};
293#endif /* CALFACTFITS_H */
Note: See TracBrowser for help on using the repository browser.