source: branches/Mars_use_drstimefiles/fact/processing/calibrate.C@ 19478

Last change on this file since 19478 was 17101, checked in by Daniela Dorner, 12 years ago
updated paths (new cluster)
File size: 7.9 KB
Line 
1/*
2#ifndef HAVE_ZLIB
3#define HAVE_ZLIB
4#endif
5
6#include <TMath.h>
7
8#include "fits.h"
9#include "ofits.h"
10
11#include "MArrayF.h"
12#include "MMath.h"
13*/
14
15using namespace std;
16
17bool Read(const char *fname, int *npix)
18{
19 std::ifstream fin(fname);
20
21 int l = 0;
22
23 while (1)
24 {
25 if (l>1439)
26 break;
27
28 TString line;
29 line.ReadLine(fin);
30 if (!fin)
31 break;
32
33 line = line.Strip(TString::kBoth);
34 if (line[0]=='#')
35 continue;
36
37 TObjArray *arr = line.Tokenize(' ');
38 int cbpx = atoi(((*arr)[1])->GetName());
39 int board = atoi(((*arr)[6])->GetName());
40 int ch = atoi(((*arr)[7])->GetName());
41 delete arr;
42
43 int hv = board*32+ch;
44 if (hv>=416)
45 return false;
46
47 int n = cbpx%10>3 ? 5 : 4;
48 if (npix[hv]>0 && npix[hv]!=n)
49 return false;
50
51 npix[hv] = n;
52 l++;
53 }
54
55 return l==1440;
56}
57
58int calibrate(UInt_t night, TString path="/fact/aux")
59{
60 int npix[416];
61 memset(npix, 0, sizeof(int)*416);
62 if (!Read("/gpfs/scratch/fact/FACTmap111030.txt", npix))
63 return 42;
64
65 UInt_t d = night%100;
66 UInt_t m = (night/100)%100;
67 UInt_t y = night/10000;
68
69 path += Form("/%04d/%02d/%02d/%d", y, m, d, night);
70 TString path2 = Form("/gpfs/scratch/fact/calibrated_currents/%d", night);
71
72 TString filec = path + ".FEEDBACK_CALIBRATION.fits";
73 TString fileu = path + ".BIAS_CONTROL_VOLTAGE.fits";
74 TString filei = path + ".BIAS_CONTROL_CURRENT.fits";
75 TString fileo = path2 + ".CALIBRATED_CURRENTS.fits";
76
77 cout << filec << endl;
78 cout << fileu << endl;
79 cout << filei << endl;
80 cout << fileo << endl;
81
82 fits fitsc1(filec.Data());
83 fits fitsc2(filec.Data());
84 if (!fitsc2)
85 {
86 cout << "ERROR - Open " << filec << " failed." << endl;
87 return 1;
88 }
89
90 fits fitsu1(fileu.Data());
91 fits fitsu2(fileu.Data());
92 if (!fitsu2)
93 {
94 cout << "ERROR - Open " << fileu << " failed." << endl;
95 return 2;
96 }
97
98 fits fitsi1(filei.Data());
99 fits fitsi2(filei.Data());
100 if (!fitsi2)
101 {
102 cout << "ERROR - Open " << filei << " failed." << endl;
103 return 3;
104 }
105
106 ofits fout(fileo.Data());
107 if (!fout)
108 {
109 cout << "ERROR - Open " << fileo << " failed." << endl;
110 return 4;
111 }
112
113 //const MTime now;
114 fout.SetStr( "TELESCOP", "FACT", "Telescope that acquired this data");
115 fout.SetStr( "PACKAGE", "MARS", "Package name");
116 fout.SetStr( "VERSION", "1.0", "Package description");
117 fout.SetFloat("EXTREL", 1.0, "Release Number");
118 //fout.SetStr( "COMPILED", __DATE__" "__TIME__, "Compile time");
119 fout.SetStr( "ORIGIN", "FACT", "Institution that wrote the file");
120 //fout.SetStr( "DATE", now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), "File creation date");
121 //fout.SetInt( "NIGHT", now.GetNightAsInt(), "Night as int");
122 fout.SetStr( "TIMESYS", "UTC", "Time system");
123 fout.SetStr( "TIMEUNIT", "d", "Time given in days w.r.t. to MJDREF");
124 fout.SetInt( "MJDREF", 40587, "MJD to UNIX time (seconds since 1970/1/1)");
125
126 fout.AddColumnDouble("Time", "MJD", "Event time: MJD - MJDREF (= Unix Time)");
127 fout.AddColumnDouble("T_cal", "MJD", "Calibration time: MJD - MJDREF (= Unix Time)");
128 fout.AddColumnInt( "N", "", "Number of valid values");
129 fout.AddColumnFloat( "I_avg", "uA", "Average calibrated current (320 channels)");
130 fout.AddColumnFloat( "I_rms", "uA", "RMS of calibrated current (320 channels)");
131 fout.AddColumnFloat( "I_med", "uA", "Median of calibrated currents (320 channels)");
132 fout.AddColumnFloat( "I_dev", "uA", "Deviation of calibrated currents (320 channels)");
133 fout.AddColumnFloat(416, "I", "uA", "My comment");
134
135 fout.WriteTableHeader("CALIBRATED");
136
137 Double_t timec, timeu, timei;
138 Double_t nextc, nextu, nexti;
139
140 if (!fitsc1.SetPtrAddress("Time", &timec))
141 {
142 cout << "ERROR couldn't get pointer to 1st Time for R" << endl;
143 return 2;
144 }
145 if (!fitsu1.SetPtrAddress("Time", &timeu))
146 {
147 cout << "ERROR couldn't get pointer to 1st Time for U" << endl;
148 return 2;
149 }
150 if (!fitsi1.SetPtrAddress("Time", &timei))
151 {
152 cout << "ERROR couldn't get pointer to 1st Time for I" << endl;
153 return 2;
154 }
155
156 if (!fitsc2.SetPtrAddress("Time", &nextc))
157 {
158 cout << "ERROR couldn't get pointer to 2nd Time for R" << endl;
159 return 2;
160 }
161 if (!fitsu2.SetPtrAddress("Time", &nextu))
162 {
163 cout << "ERROR couldn't get pointer to 2nd Time for U" << endl;
164 return 2;
165 }
166 if (!fitsi2.SetPtrAddress("Time", &nexti))
167 {
168 cout << "ERROR couldn't get pointer to 2nd Time for I" << endl;
169 return 2;
170 }
171
172 Float_t R[416];
173 Float_t U[416];
174 UShort_t I[416];
175
176 if (!fitsc1.SetPtrAddress("R", R))
177 {
178 cout << "ERROR couldn't get pointer to R" << endl;
179 return 2;
180 }
181 if (!fitsu1.SetPtrAddress("Uout", U))
182 {
183 cout << "ERROR couldn't get pointer to U" << endl;
184 return 2;
185 }
186 if (!fitsi1.SetPtrAddress("I", I))
187 {
188 cout << "ERROR couldn't get pointer to I" << endl;
189 return 2;
190 }
191
192 fitsc2.GetNextRow();
193 fitsu2.GetNextRow();
194 fitsi2.GetNextRow();
195
196 if (!fitsc1.GetNextRow())
197 return 0;
198 if (!fitsc2.GetNextRow())
199 nextc = 1e50;
200
201 while (1)
202 {
203 // First the next two rows (voltage, current)
204 // with identical time stamp
205 if (nextu<nextc && nextu<nexti)
206 {
207 if (!fitsu1.GetNextRow())
208 break;
209 fitsu2.GetNextRow();
210 continue;
211 }
212
213 if (nexti<nextc && nexti<nextu)
214 {
215 if (!fitsi1.GetNextRow())
216 break;
217 fitsi2.GetNextRow();
218 continue;
219 }
220
221 // At this point we have found currents and voltages
222 // with the same time-stamp
223 if (nexti!=nextu)
224 break;
225
226 // Now read the next row, voltages and currenst, which
227 // should be evaluated
228 if (!fitsu1.GetNextRow() || !fitsi1.GetNextRow())
229 break;
230
231 // Read also one event in advance to get the next time-stamps
232 fitsu2.GetNextRow();
233 fitsi2.GetNextRow();
234
235 // Read new calibration events as long as the difference
236 // in time to the next calibration evenet is smaller
237 // than to the current one. In this way always the
238 // closest one is used.
239 while (fabs(nextc-nextu)<fabs(timec-timeu))
240 {
241 // If there is no next one, make sure that
242 // the the condition is never met again.
243 if (!fitsc1.GetNextRow())
244 {
245 nextc=1e50;
246 break;
247 }
248 fitsc2.GetNextRow();
249 }
250
251 // Calibrate the data (subtract offset)
252 MArrayF rc(9+416);
253 MArrayF stat(416);
254
255 Float_t *cur = rc.GetArray()+9;
256
257 UInt_t cnt = 0;
258 for (int i=0; i<416; i++)
259 {
260 if (R[i]<=0)
261 continue;
262
263 // Convert dac counts to uA
264 // Measued current minus leakage current (bias crate calibration)
265 cur[i] = I[i] * 5000./4096 - U[i]/R[i]*1e6;
266
267 if (npix[i]>0)
268 {
269 cur[i] /= npix[i];
270 stat[cnt++] = cur[i];
271 }
272 }
273
274 memcpy(rc.GetArray(), &timeu, sizeof(Double_t));
275 memcpy(rc.GetArray()+2, &timec, sizeof(Double_t));
276 memcpy(rc.GetArray()+4, &cnt, sizeof(UInt_t));
277
278 Double_t med;
279 rc[5] = TMath::Mean(cnt, stat.GetArray());
280 rc[6] = TMath::RMS(cnt, stat.GetArray());
281 rc[8] = MMath::MedianDev(cnt, stat.GetArray(), med);
282 rc[7] = med;
283
284 fout.WriteRow(rc.GetArray(), rc.GetSize()*sizeof(Float_t));
285 }
286
287 return 0;
288}
Note: See TracBrowser for help on using the repository browser.