source: trunk/Mars/fact/analysis/gain/extract_aux.C@ 17370

Last change on this file since 17370 was 17370, checked in by tbretz, 11 years ago
Renamed and added more auxiliary infomrations.
File size: 8.9 KB
Line 
1#include <iostream>
2
3#include <TString.h>
4#include <TFile.h>
5#include <TGraph.h>
6
7#include "zfits.h"
8#include "PixelMap.h"
9#include "Interpolator2D.h"
10
11#include "MMath.h"
12#include "MHCamera.h"
13#include "MGeomCamFACT.h"
14
15PixelMap pmap;
16
17bool extract_current(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
18{
19 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FEEDBACK_CALIBRATED_CURRENTS.fits",
20 day/10000, (day/100)%100, day%100, day);
21
22 // =====================================================================
23
24 // Now we read the temperatures of the sensors during the
25 // single pe run and average them
26 fits faux(naux.Data());
27 if (!faux)
28 {
29 cout << "Could not open " << naux << endl;
30 return false;
31 }
32
33 Double_t time;
34 Float_t I[416];
35
36 if (!faux.SetPtrAddress("Time", &time))
37 return false;
38 if (!faux.SetPtrAddress("I", I))
39 return false;
40
41 TArrayD avg(320);
42 TArrayD rms(320);
43 TArrayI cnt(320);
44
45 while (faux.GetNextRow())
46 {
47 if (time<tstart || time>tstop)
48 continue;
49
50 TArrayD med(320);
51
52 double mean = 0;
53
54 int j=0;
55 for (int i=0; i<320; i++)
56 {
57 if (i!=66 && i!=191 && i!=193 && i!=17 && i!=206)
58 med[j++] = I[i];
59
60 avg[i] += I[i];
61 rms[i] += I[i]*I[i];
62 cnt[i]++;
63 }
64
65 Double_t median;
66 MMath::MedianDev(315, med.GetArray(), median);
67 g.SetPoint(g.GetN(), time, median);
68 }
69
70 if (cnt[0]==0)
71 return 2;
72
73 for (int i=0; i<320; i++)
74 {
75 avg[i] /= cnt[i];
76 rms[i] /= cnt[i];
77 rms[i] -= avg[i]*avg[i];
78 rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
79 }
80
81 double mn = 0;
82 for (int i=0; i<1440; i++)
83 {
84 mn += avg[pmap.index(i).hv()];
85 havg.SetBinContent(i+1, avg[pmap.index(i).hv()]);
86 hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]);
87 }
88
89 return true;
90}
91
92bool extract_dac(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
93{
94 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.BIAS_CONTROL_DAC.fits",
95 day/10000, (day/100)%100, day%100, day);
96
97 // =====================================================================
98
99 // Now we read the temperatures of the sensors during the
100 // single pe run and average them
101 fits faux(naux.Data());
102 if (!faux)
103 {
104 cout << "Could not open " << naux << endl;
105 return false;
106 }
107
108 Double_t time;
109 UShort_t dac[416];
110
111 if (!faux.SetPtrAddress("Time", &time))
112 return false;
113 if (!faux.SetPtrAddress("U", dac))
114 return false;
115
116 TArrayD avg(320);
117 TArrayD rms(320);
118 TArrayI cnt(320);
119
120 while (faux.GetNextRow())
121 {
122 if (time<tstart || time>tstop)
123 continue;
124
125 TArrayD med(320);
126
127 for (int i=0; i<320; i++)
128 {
129 med[i] = dac[i];
130 avg[i] += dac[i];
131 rms[i] += dac[i]*dac[i];
132 cnt[i]++;
133 }
134
135 Double_t median;
136 MMath::MedianDev(320, med.GetArray(), median);
137 g.SetPoint(g.GetN(), time, median);
138 }
139
140 if (cnt[0]==0)
141 return 2;
142
143 for (int i=0; i<320; i++)
144 {
145 avg[i] /= cnt[i];
146 rms[i] /= cnt[i];
147 rms[i] -= avg[i]*avg[i];
148 rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
149 }
150
151 for (int i=0; i<1440; i++)
152 {
153 havg.SetBinContent(i+1, avg[pmap.index(i).hv()]);
154 hrms.SetBinContent(i+1, rms[pmap.index(i).hv()]);
155 }
156
157 return true;
158}
159
160bool extract_temp(UInt_t day, double tstart, double tstop, MHCamera &havg, MHCamera &hrms, TGraph &g)
161{
162 TString naux = Form("/fact/aux/%4d/%02d/%02d/%8d.FSC_CONTROL_TEMPERATURE.fits",
163 day/10000, (day/100)%100, day%100, day);
164
165 Interpolator2D interpol;
166
167 const auto sens = Interpolator2D::ReadGrid("operation/sensor-pos.txt");
168 if (sens.size()!=31)
169 {
170 cout << "Reading sensor-pos.txt failed: " << interpol.getInputGrid().size() << endl;
171 return false;
172 }
173
174 // =====================================================================
175
176 // Now we read the temperatures of the sensors during the
177 // single pe run and average them
178 fits faux(naux.Data());
179 if (!faux)
180 {
181 cout << "Could not open " << naux << endl;
182 return false;
183 }
184
185 Double_t time;
186 Float_t temp[31];
187
188 if (!faux.SetPtrAddress("Time", &time))
189 return false;
190 if (!faux.SetPtrAddress("T_sens", temp))
191 return false;
192
193 TArrayD avg(31);
194 TArrayD rms(31);
195 TArrayI cnt(31);
196
197 while (faux.GetNextRow())
198 {
199 if (time<tstart || time>tstop)
200 continue;
201
202 double mean = 0;
203 int count = 0;
204
205 for (int i=0; i<31; i++)
206 {
207 if (temp[i]==0)
208 continue;
209
210 double T = temp[i];
211
212 mean += T;
213 count++;
214
215 avg[i] += T;
216 rms[i] += T*T;
217 cnt[i]++;
218 }
219
220 g.SetPoint(g.GetN(), time, mean/count);
221
222 }
223
224 vector<double> zavg;
225 vector<double> zrms;
226
227 vector<Interpolator2D::vec> pos;
228 for (int i=0; i<31; i++)
229 {
230 if (cnt[i]>0)
231 {
232 avg[i] /= cnt[i];
233 rms[i] /= cnt[i];
234 rms[i] -= avg[i]*avg[i];
235 rms[i] = rms[i]<0 ? 0 : sqrt(rms[i]);
236
237 pos.push_back(sens[i]);
238 zavg.push_back(avg[i]);
239 zrms.push_back(rms[i]);
240 }
241 }
242
243 // ================================================================
244
245 interpol.SetInputGrid(pos);
246
247 if (!interpol.ReadOutputGrid("operation/bias-positions.txt"))
248 {
249 cout << "Error setting output grid." << endl;
250 return false;
251 }
252 if (interpol.getOutputGrid().size()!=320)
253 {
254 cout << "Reading bias-positions.txt failed: " << interpol.getOutputGrid().size() << endl;
255 return false;
256 }
257
258 // ================================================================
259
260 const auto &input = interpol.getInputGrid();
261
262 TGraph sensors(input.size());
263
264 for (unsigned int i=0; i<input.size(); i++)
265 sensors.SetPoint(i, input[i].x, input[i].y);
266
267 const vector<double> rc_avg = interpol.Interpolate(zavg);
268 const vector<double> rc_rms = interpol.Interpolate(zrms);
269
270 // ================================================================
271
272 for (int i=0; i<1440; i++)
273 {
274 havg.SetBinContent(i+1, rc_avg[pmap.index(i).hv()]);
275 hrms.SetBinContent(i+1, rc_rms[pmap.index(i).hv()]);
276 }
277
278 return true;
279
280}
281
282void extract_aux(UInt_t day=20131114, UInt_t run=136, const char *output=0)
283{
284 if (!pmap.Read("operation/FACTmap111030.txt"))
285 {
286 cout << "FACTmap111030.txt not found." << endl;
287 return;
288 }
289
290 // =====================================================================
291
292 TString nraw = Form("/fact/raw/%4d/%02d/%02d/%8d_%03d.fits.fz",
293 day/10000, (day/100)%100, day%100, day, run);
294
295
296 zfits fraw(nraw.Data());
297 if (!fraw)
298 {
299 cout << "Open raw file failed: " << nraw << endl;
300 return;
301 }
302
303 double tstart = fraw.GetUInt("TSTARTI") + fraw.GetFloat("TSTARTF");
304 double tstop = fraw.GetUInt("TSTOPI") + fraw.GetFloat("TSTOPF");
305
306 // =====================================================================
307
308 MGeomCamFACT geom;
309
310 TGraph gcur;
311 MHCamera hcur_m(geom);
312 MHCamera hcur_r(geom);
313 if (!extract_current(day, tstart, tstop, hcur_m, hcur_r, gcur))
314 return;
315
316 TGraph gtmp;
317 MHCamera htmp_m(geom);
318 MHCamera htmp_r(geom);
319 if (!extract_temp(day, tstart, tstop, htmp_m, htmp_r, gtmp))
320 return;
321
322 TGraph gdac;
323 MHCamera hdac_m(geom);
324 MHCamera hdac_r(geom);
325 if (!extract_dac(day, tstart, tstop, hdac_m, hdac_r, gdac))
326 return;
327
328 // ================================================================
329
330 TCanvas *c = new TCanvas;
331 c->Divide(3, 3);
332
333 c->cd(1);
334 hcur_m.SetName("CurrentAvg");
335 hcur_m.SetAllUsed();
336 hcur_m.DrawCopy();
337 c->cd(4);
338 hcur_r.SetName("CurrentRms");
339 hcur_r.SetAllUsed();
340 hcur_r.DrawCopy();
341 c->cd(7);
342 gcur.SetName("CurrentTime");
343 gcur.SetMarkerStyle(kFullDotMedium);
344 gcur.SetMinimum(0);
345 gcur.DrawClone("AP");
346
347 c->cd(2);
348 htmp_m.SetName("TempAvg");
349 htmp_m.SetAllUsed();
350 htmp_m.DrawCopy();
351 c->cd(5);
352 htmp_r.SetName("TempRms");
353 htmp_r.SetAllUsed();
354 htmp_r.DrawCopy();
355 c->cd(8);
356 gtmp.SetName("TempTime");
357 gtmp.SetMarkerStyle(kFullDotMedium);
358 gtmp.SetMinimum(0);
359 gtmp.DrawClone("AP");
360
361 c->cd(3);
362 hdac_m.SetName("DacAvg");
363 hdac_m.SetAllUsed();
364 hdac_m.DrawCopy();
365 c->cd(6);
366 hdac_r.SetName("DacRms");
367 hdac_r.SetAllUsed();
368 hdac_r.DrawCopy();
369 c->cd(9);
370 gdac.SetName("DacTime");
371 gdac.SetMarkerStyle(kFullDotMedium);
372 gdac.SetMinimum(0);
373 gdac.DrawClone("AP");
374
375 if (output)
376 c->SaveAs(output);
377}
Note: See TracBrowser for help on using the repository browser.