source: branches/Mars_McMismatchStudy/fact/analysis/gain/extract_aux.C@ 19366

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