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

Last change on this file since 19832 was 17730, checked in by dneise, 12 years ago
Removed the dependency of SlowData class from the special Python Accessors in fits.h So now we can use the C++ header files from Mars/mcore. makelibs.C is not needed anymore. In order to make the so-files just call pyfact.py For testing one can do: python -i pyfact.py <filename> After a lot of compiler warnings one should get: len(test_m1) 0 len(test_m2) 0 <ROOT.pair<string,fits::Entry> object at 0x...> <ROOT.pair<string,fits::Table::Column> object at 0x...> Well ... for now.
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 "extern_Mars_mcore/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.