source: trunk/Mars/hawc/quality.C@ 19999

Last change on this file since 19999 was 19921, checked in by tbretz, 6 years ago
quality macro for HAWC's Eye
File size: 10.8 KB
Line 
1#include <iostream>
2#include "TROOT.h"
3#include "TSystem.h"
4#include "TGraph.h"
5#include "TArrayD.h"
6#include "TCanvas.h"
7#include "TH1.h"
8#include "TAxis.h"
9#include "TText.h"
10#include "TLine.h"
11
12#include "fits.h"
13#include "MTime.h"
14
15#ifndef __CLING__
16#include <algorithm>
17#include <functional>
18#include <vector>
19#endif
20
21Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0)
22{
23 TArrayD *arr0 = vec[0];
24 TArrayD *arr1 = vec[1];
25 TArrayD *arr2 = vec[2];
26
27 for (int i=0; i<arr0->GetSize(); i++)
28 {
29 Double_t start = (*arr1)[i];
30 Double_t stop = (*arr2)[i];
31
32 if (stop>start+305./24/3600)
33 stop = start+305./24/3600;
34
35 if (t0>start-range && t0<stop+range)
36 //{
37 // if (fmod(t0, 1)>4./24 && fmod(t0,1)<4.1/24)
38 // cout << t0-start << " " <<t0 << " " << stop-t0 << " " << start-15779 << " " << stop-15779 << " " << (*arr0)[i] << endl;
39 return kTRUE;
40 //}
41 }
42
43 return arr0->GetSize()==0;
44}
45
46Int_t PlotThresholds(double start, TArrayD **vec, TString fname)
47{
48 fname += ".FTU_CONTROL_DAC.fits";
49
50 fits file(fname.Data());
51 if (!file)
52 {
53 cerr << fname << ": " << gSystem->GetError() << endl;
54 return -2;
55 }
56
57 //cout << fname << endl;
58
59 Double_t time;
60 UShort_t th[4];
61
62 if (!file.SetPtrAddress("Time", &time))
63 return -1;
64
65 if (!file.SetPtrAddress("dac_ch", th))
66 return -1;
67
68 TGraph g;
69 g.SetName("Threshold");
70 while (file.GetNextRow())
71 {
72 if (Contains(vec, time, 10./(24*3600)))
73 g.SetPoint(g.GetN(), time*24*3600, th[0]);
74 }
75
76
77 gROOT->SetSelectedPad(0);
78 gPad->GetCanvas()->cd(1);
79 gPad->SetGrid();
80 gPad->SetTopMargin(0);
81 gPad->SetRightMargin(0.001);
82 gPad->SetLeftMargin(0.04);
83 gPad->SetBottomMargin(0);
84
85 TH1C frame("Threshold_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
86 frame.SetMinimum(300);
87 frame.SetMaximum(1200);
88 frame.SetStats(kFALSE);
89 frame.GetXaxis()->SetTimeDisplay(true);
90 frame.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
91 frame.GetXaxis()->SetLabelSize(0.12);
92 frame.GetYaxis()->SetLabelSize(0.1);
93 frame.GetYaxis()->SetTitle("THRESHOLD");
94 frame.GetYaxis()->CenterTitle();
95 frame.GetYaxis()->SetTitleOffset(0.2);
96 frame.GetYaxis()->SetTitleSize(0.1);
97 frame.DrawCopy();
98
99 g.SetMarkerStyle(kFullDotMedium);
100 g.DrawClone("P");
101
102 return 0;
103}
104
105Int_t PlotCurrent(double start, TArrayD **vec, TString fname)
106{
107 fname += ".BIAS_CONTROL_AVERAGE_DATA.fits";
108
109 fits file(fname.Data());
110 if (!file)
111 {
112 cerr << fname << ": " << gSystem->GetError() << endl;
113 return -2;
114 }
115
116 //cout << fname << endl;
117
118 Double_t time;
119 Float_t Iavg[64];
120 Float_t Tavg[64];
121 Float_t Tpsu[2];
122 Float_t Vavg[64];
123
124 if (!file.SetPtrAddress("Time", &time))
125 return -1;
126
127 if (!file.SetPtrAddress("Iavg", Iavg))
128 return -1;
129
130 if (!file.SetPtrAddress("Tavg", Tavg))
131 return -1;
132
133 if (!file.SetPtrAddress("Tavg_psu", Tpsu))
134 return -1;
135
136 if (!file.SetPtrAddress("Vavg", Vavg))
137 return -1;
138
139 TGraph g1;
140 TGraph g2;
141 TGraph g3;
142 TGraph g4;
143 g1.SetName("CurrentAvg");
144 g2.SetName("TemperatureAvg");
145 g3.SetName("TempPSUavg");
146 g4.SetName("VoltageAvg");
147
148 while (file.GetNextRow())
149 if (Contains(vec, time))
150 {
151 sort(Iavg, Iavg+64);
152 sort(Tavg, Tavg+64);
153 sort(Vavg, Vavg+64);
154
155 g1.SetPoint(g1.GetN(), time*24*3600, Iavg[31]);
156 g2.SetPoint(g2.GetN(), time*24*3600, Tavg[31]);
157 g3.SetPoint(g3.GetN(), time*24*3600, (Tpsu[0]+Tpsu[1])/2);
158 g4.SetPoint(g4.GetN(), time*24*3600, Vavg[31]);
159 }
160
161
162 gROOT->SetSelectedPad(0);
163 gPad->GetCanvas()->cd(2);
164 gPad->SetGrid();
165 gPad->SetTopMargin(0);
166 gPad->SetRightMargin(0.001);
167 gPad->SetLeftMargin(0.04);
168 gPad->SetBottomMargin(0);
169
170 TH1C frame1("Current_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
171 frame1.SetMinimum(0);
172 frame1.SetMaximum(1000);
173 frame1.SetStats(kFALSE);
174 frame1.GetXaxis()->SetTimeDisplay(true);
175 frame1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
176 frame1.GetXaxis()->SetLabelSize(0.12);
177 frame1.GetYaxis()->SetLabelSize(0.1);
178 frame1.GetYaxis()->SetTitle("CURRENT");
179 frame1.GetYaxis()->CenterTitle();
180 frame1.GetYaxis()->SetTitleOffset(0.2);
181 frame1.GetYaxis()->SetTitleSize(0.1);
182 frame1.DrawCopy();
183
184 g1.SetMarkerStyle(kFullDotMedium);
185 g1.DrawClone("P");
186
187
188
189 gROOT->SetSelectedPad(0);
190 gPad->GetCanvas()->cd(5);
191 gPad->SetGrid();
192 gPad->SetTopMargin(0);
193 gPad->SetRightMargin(0.001);
194 gPad->SetLeftMargin(0.04);
195 gPad->SetBottomMargin(0);
196
197 TH1C frame4("Voltage_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
198 frame4.SetMinimum(29);
199 frame4.SetMaximum(30);
200 frame4.SetStats(kFALSE);
201 frame4.GetXaxis()->SetTimeDisplay(true);
202 frame4.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
203 frame4.GetXaxis()->SetLabelSize(0.12);
204 frame4.GetYaxis()->SetLabelSize(0.1);
205 frame4.GetYaxis()->SetTitle("VOLTAGE");
206 frame4.GetYaxis()->CenterTitle();
207 frame4.GetYaxis()->SetTitleOffset(0.2);
208 frame4.GetYaxis()->SetTitleSize(0.1);
209 frame4.DrawCopy();
210
211 g4.SetMarkerStyle(kFullDotMedium);
212 g4.DrawClone("P");
213
214
215
216 gROOT->SetSelectedPad(0);
217 gPad->GetCanvas()->cd(6);
218 gPad->SetGrid();
219 gPad->SetTopMargin(0);
220 gPad->SetRightMargin(0.001);
221 gPad->SetLeftMargin(0.04);
222 gPad->SetBottomMargin(0);
223
224 TH1C frame2("Temperature_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
225 frame2.SetMinimum(0);
226 frame2.SetMaximum(40);
227 frame2.SetStats(kFALSE);
228 frame2.GetXaxis()->SetTimeDisplay(true);
229 frame2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
230 frame2.GetXaxis()->SetLabelSize(0.12);
231 frame2.GetYaxis()->SetLabelSize(0.1);
232 frame2.GetYaxis()->SetTitle("TEMPERATURE");
233 frame2.GetYaxis()->CenterTitle();
234 frame2.GetYaxis()->SetTitleOffset(0.2);
235 frame2.GetYaxis()->SetTitleSize(0.1);
236 frame2.DrawCopy();
237
238 g2.SetMarkerStyle(kFullDotMedium);
239 g2.DrawClone("P");
240
241 g3.SetMarkerColor(kBlue);
242 g3.SetMarkerStyle(kFullDotMedium);
243 g3.DrawClone("P");
244
245 return 0;
246}
247
248Int_t PlotRate(double start, TArrayD **vec, TString fname)
249{
250 fname += ".FTM_CONTROL_DATA.fits";
251
252 fits file(fname.Data());
253 if (!file)
254 {
255 cerr << fname << ": " << gSystem->GetError() << endl;
256 return -2;
257 }
258
259 //cout << fname << endl;
260
261 Double_t time;
262 uint32_t trg_counter;
263 uint32_t dead_time, run_time;
264
265 if (!file.SetPtrAddress("Time", &time))
266 return -1;
267
268 if (!file.SetPtrAddress("trg_counter", &trg_counter))
269 return -1;
270 if (!file.SetPtrAddress("dead_time", &dead_time))
271 return -1;
272 if (!file.SetPtrAddress("run_time", &run_time))
273 return -1;
274
275 TGraph g1, g2;
276 g1.SetName("TriggerRate");
277 g2.SetName("RelOnTime");
278
279 uint32_t prev = 0;
280
281 while (file.GetNextRow())
282 if (Contains(vec, time))
283 {
284 g1.SetPoint(g1.GetN(), time*24*3600, 1e8*(trg_counter-prev)/(run_time-dead_time));
285 g2.SetPoint(g2.GetN(), time*24*3600, 1.-float(dead_time)/run_time);
286 prev = trg_counter;
287 }
288
289
290 gROOT->SetSelectedPad(0);
291 gPad->GetCanvas()->cd(3);
292 gPad->SetGrid();
293 gPad->SetTopMargin(0);
294 gPad->SetBottomMargin(0);
295 gPad->SetRightMargin(0.001);
296 gPad->SetLeftMargin(0.04);
297
298 TH1C frame1("Rate_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
299 frame1.SetMinimum(0);
300 frame1.SetMaximum(90);
301 frame1.SetStats(kFALSE);
302 frame1.GetXaxis()->SetTimeDisplay(true);
303 frame1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
304 frame1.GetXaxis()->SetLabelSize(0.12);
305 frame1.GetYaxis()->SetLabelSize(0.1);
306 frame1.GetYaxis()->SetTitle("TRIGGER RATE");
307 frame1.GetYaxis()->CenterTitle();
308 frame1.GetYaxis()->SetTitleOffset(0.2);
309 frame1.GetYaxis()->SetTitleSize(0.1);
310 frame1.DrawCopy();
311
312 g1.SetMarkerStyle(kFullDotMedium);
313 g1.DrawClone("P");
314
315
316
317 gROOT->SetSelectedPad(0);
318 gPad->GetCanvas()->cd(4);
319 gPad->SetGrid();
320 gPad->SetTopMargin(0);
321 gPad->SetBottomMargin(0);
322 gPad->SetRightMargin(0.001);
323 gPad->SetLeftMargin(0.04);
324
325 TH1C frame2("OnTime_frame", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
326 frame2.SetMinimum(0);
327 frame2.SetMaximum(1);
328 frame2.SetStats(kFALSE);
329 frame2.GetXaxis()->SetTimeDisplay(true);
330 frame2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
331 frame2.GetXaxis()->SetLabelSize(0.12);
332 frame2.GetYaxis()->SetLabelSize(0.1);
333 frame2.GetYaxis()->SetTitle("RELATIVE ON TIME");
334 frame2.GetYaxis()->CenterTitle();
335 frame2.GetYaxis()->SetTitleOffset(0.2);
336 frame2.GetYaxis()->SetTitleSize(0.1);
337 frame2.DrawCopy();
338
339 g2.SetMarkerStyle(kFullDotMedium);
340 g2.DrawClone("P");
341
342 return 0;
343}
344
345int quality(const char *basepath, UInt_t y=0, UInt_t m=0, UInt_t d=0, const char *outpath=0)
346{
347 // To get correct dates in the histogram you have to add
348 // the MJDREF offset (should be 40587) and 9131.
349
350 if (y==0)
351 {
352 UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt();
353 y = nt/10000;
354 m = (nt/100)%100;
355 d = nt%100;
356
357 cout << y << "/" << m << "/" << d << endl;
358 }
359
360 TString fname=Form("%s/%04d/%02d/%02d/%04d%02d%02d", basepath, y, m, d, y, m, d);
361
362 cout << "quality" << endl;
363 cout << "-------" << endl;
364 cout << endl;
365 cout << "Night: " << Form("%04d-%02d-%02d", y, m, d) << endl;
366 cout << endl;
367
368 double start = MTime(y, m, d).GetMjd()-40587;
369
370 TArrayD run, beg, end;
371
372 TArrayD *runs[3] = { &run, &beg, &end };
373
374 //check if the sqm was already installed on the telescope
375 TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 1120);
376 c->Divide(1, 7, 1e-5, 1e-5);
377
378 PlotThresholds(start, runs, fname);
379 PlotCurrent(start, runs, fname);
380 PlotRate(start, runs, fname);
381
382 gROOT->SetSelectedPad(0);
383 gPad->GetCanvas()->cd(7);
384 gPad->SetTopMargin(0);
385 gPad->SetBottomMargin(0.99);
386 gPad->SetRightMargin(0.001);
387 gPad->SetLeftMargin(0.04);
388
389 TH1C frame1("Axis", "", 1000, (start+0.45)*24*3600, (start+1.025)*24*3600);
390 frame1.SetMinimum(0);
391 frame1.SetMaximum(1);
392 frame1.SetStats(kFALSE);
393 frame1.GetXaxis()->SetTimeDisplay(true);
394 frame1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 17:00:00 GMT");
395 frame1.GetXaxis()->SetLabelSize(0.2);
396 frame1.GetXaxis()->SetTitleSize(0.2);
397 frame1.GetXaxis()->SetTitle("HAWC Local Time");
398 frame1.GetXaxis()->CenterTitle();
399 frame1.DrawCopy();
400
401 if (outpath)
402 c->SaveAs(outpath);
403
404 return 0;
405}
Note: See TracBrowser for help on using the repository browser.