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

Last change on this file since 16858 was 14487, checked in by neise, 12 years ago
debugged pointer f**k in implementation of despiker in calfactfits
File size: 8.9 KB
Line 
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
52class CalFactFits
53{
54public:
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")*2];
142
143 calibfile.GetRow(0);
144
145int orig_index, new_index1, new_index2;
146 for (int pix=0 ; pix < (int)calib_npix; ++pix)
147 {
148 for (int sl = 0; sl < (int)calib_blm_size; ++sl)
149 {
150 orig_index = pix * (int)calib_blm_size + sl;
151 new_index1 = 2*pix * (int)calib_blm_size + sl;
152 new_index2 = (2*pix+1) * (int)calib_blm_size + sl;
153
154 baseline[new_index2] = baseline[new_index1] = double(calib_baselinemean[orig_index]) *4096./2000.;
155 }
156 }
157 for (int pix=0 ; pix < (int)calib_npix; ++pix)
158 {
159 for (int sl = 0; sl < (int)calib_gm_size; ++sl)
160 {
161 orig_index = pix * (int)calib_gm_size + sl;
162 new_index1 = 2*pix * (int)calib_gm_size + sl;
163 new_index2 = (2*pix+1) * (int)calib_gm_size + sl;
164
165 gain[new_index2] = gain[new_index1] = double(calib_gainmean[orig_index]) /( 1907.35 /4096. *2000.);
166 }
167 }
168
169 for (int pix=0 ; pix < (int)calib_npix; ++pix)
170 {
171 for (int sl = 0; sl < (int)calib_tom_size; ++sl)
172 {
173 orig_index = pix * (int)calib_tom_size + sl;
174 new_index1 = 2*pix * (int)calib_tom_size + sl;
175 new_index2 = (2*pix+1) * (int)calib_tom_size + sl;
176
177 trigger_offset[new_index2] = trigger_offset[new_index1] = double(calib_triggeroffsetmean[orig_index]) *4096. /2000.;
178 }
179 }
180
181
182
183 //Set the event pointers
184// std::cout << "...Setting event pointers..." << std::endl;
185 datafile.SetRefAddress("EventNum", event_id);
186 datafile.SetRefAddress("TriggerType", event_triggertype);
187
188 event_data = new int16_t[data_ndata];
189 datafile.SetPtrAddress("Data", event_data, data_ndata);
190
191 event_offset = new int16_t[datafile.GetN("StartCellData")];
192 datafile.SetPtrAddress("StartCellData", event_offset, datafile.GetN("StartCellData"));
193
194 event_boardtimes = new int32_t[datafile.GetN("BoardTime")];
195 datafile.SetPtrAddress("BoardTime", event_boardtimes, datafile.GetN("BoardTime"));
196 }
197
198 ~CalFactFits() //Standard destructor
199 {
200 delete[] calib_baselinemean;
201 delete[] calib_gainmean;
202 delete[] calib_triggeroffsetmean;
203 delete[] event_data;
204 delete[] event_offset;
205 delete[] event_boardtimes;
206 delete[] baseline;
207 delete[] gain;
208 delete[] trigger_offset;
209 }
210
211 bool GetCalEvent() //Read calibrated event into the event variables
212 {
213 if(!npcaldata) {
214 ostringstream str;
215 str << "Pointer to the calibrated data not initialized!";
216#ifdef __EXCEPTIONS
217 throw runtime_error(str.str());
218#else
219 gLog << ___err___ << "ERROR - " << str.str() << endl;
220 return false;
221#endif
222 }
223 if(datafile.GetNextRow() == false) {
224// std::cout << "Last event reached..." << std::endl;
225 return false;
226 }
227 else
228 {
229 UInt_t drs_calib_offset;
230 double raw, raw_bsl, shifted;
231 for(UInt_t pixel = 0; pixel < data_npix; pixel++)
232 {
233 for(UInt_t slice = 0; slice < data_nroi; slice++)
234 {
235 drs_calib_offset = slice + event_offset[pixel];
236 raw = (double)event_data[pixel*data_nroi+slice];
237 raw_bsl = raw - baseline[pixel*calib_blm_size*2+drs_calib_offset];
238 shifted = raw_bsl - trigger_offset[pixel*data_nroi*2+slice];
239 npcaldata[pixel*data_nroi+slice] = shifted / gain[pixel*calib_blm_size*2+drs_calib_offset];
240 //Note: data_nroi=calib_nroi, calib_blm_size=calib_gm_size
241 }
242
243 // stolen from TBs MARS ... DrsCalib.h or something like that
244 double *p = npcaldata+(pixel*data_nroi);
245
246 // second loop for despiking
247 for (size_t i=1; i<data_nroi-2; i++)
248 {
249 if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
250 {
251 //std::cout << "single cand:\t" << pixel << "\t" << i << std::endl;
252 p[i] = (p[i-1]+p[i+1])/2;
253 }
254 if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
255 {
256 //std::cout << "double cand:\t" << pixel << "\t" << i << std::endl;
257 p[i] = (p[i-1]+p[i+2])/2;
258 p[i+1] = p[i];
259 }
260 }
261
262 }
263 }
264 return true;
265 }
266
267 void SetNpcaldataPtr(double *numpyptr) //Set the pointer for the calibrated data to the numpy array
268 {
269 npcaldata = numpyptr;
270 return;
271 }
272};
273#endif /* CALFACTFITS_H */
Note: See TracBrowser for help on using the repository browser.