source: branches/removing_cpp11_features/fact/plots/quality.C@ 19260

Last change on this file since 19260 was 17956, checked in by tbretz, 11 years ago
Replaced current prediction by the new central current prediction.
File size: 17.1 KB
Line 
1#include <algorithm>
2#include <functional>
3
4Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0)
5{
6 TArrayD *arr0 = vec[0];
7 TArrayD *arr1 = vec[1];
8 TArrayD *arr2 = vec[2];
9
10 for (int i=0; i<arr0->GetSize(); i++)
11 {
12 Double_t start = (*arr1)[i];
13 Double_t stop = (*arr2)[i];
14
15 if (stop>start+305./24/3600)
16 stop = start+305./24/3600;
17
18 if (t0>start-range && t0<stop+range)
19 //{
20 // if (fmod(t0, 1)>4./24 && fmod(t0,1)<4.1/24)
21 // cout << t0-start << " " <<t0 << " " << stop-t0 << " " << start-15779 << " " << stop-15779 << " " << (*arr0)[i] << endl;
22 return kTRUE;
23 //}
24 }
25
26 return arr0->GetSize()==0;
27}
28
29Int_t PlotThresholds(TArrayD **vec, TString fname)
30{
31 fname += ".RATE_CONTROL_THRESHOLD.fits";
32
33 fits file(fname.Data());
34 if (!file)
35 {
36 cerr << fname << ": " << gSystem->GetError() << endl;
37 return -2;
38 }
39
40 //cout << fname << endl;
41
42 Double_t time;
43 UShort_t th;
44
45 if (!file.SetPtrAddress("Time", &time))
46 return -1;
47
48 if (!file.SetPtrAddress("threshold", &th))
49 return -1;
50
51 TGraph g;
52 g.SetName("Threshold");
53
54 while (file.GetNextRow())
55 {
56 if (Contains(vec, time, 10./(24*3600)))
57 g.SetPoint(g.GetN(), time*24*3600, th);
58 }
59
60 g.SetMinimum(281);
61 g.SetMarkerStyle(kFullDotMedium);
62 g.GetXaxis()->SetTimeDisplay(true);
63 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
64 g.GetXaxis()->SetLabelSize(0.12);
65 g.GetYaxis()->SetLabelSize(0.1);
66 g.GetYaxis()->SetTitle("THRESHOLD");
67 g.GetYaxis()->SetTitleOffset(0.2);
68 g.GetYaxis()->SetTitleSize(0.1);
69 g.DrawClone("AP");
70
71 return 0;
72}
73
74#include <vector>
75#include <pair>
76
77vector<pair<double, Nova::EquPosn>> vecp;
78
79Nova::EquPosn FindPointing(Double_t time)
80{
81 for (int i=0; i<vecp.size(); i++)
82 if (time<vecp[i].first)
83 {
84 if (i==0)
85 return Nova::EquPosn();
86 else
87 return vecp[i-1].second;
88 }
89
90 return vecp[vecp.size()-1].second;
91}
92
93Float_t Prediction(Double_t time)
94{
95 Double_t jd = time + 40587 + 2400000.5;
96
97 // Get source position
98 Nova::EquPosn pos = FindPointing(time);
99
100 return FACT::PredictI(jd, pos);
101}
102
103Int_t ReadSources(TString fname)
104{
105 fname += ".DRIVE_CONTROL_SOURCE_POSITION.fits";
106
107 fits file(fname.Data());
108 if (!file)
109 {
110 cerr << fname << ": " << gSystem->GetError() << endl;
111 return -2;
112 }
113
114 Double_t time, ra, dec;
115 if (!file.SetPtrAddress("Time", &time))
116 return -1;
117 if (!file.SetPtrAddress("Ra_cmd", &ra))
118 return -1;
119 if (!file.SetPtrAddress("Dec_cmd", &dec))
120 return -1;
121
122 while (file.GetNextRow())
123 {
124 Nova::EquPosn p;
125 p.ra = ra*15;
126 p.dec = dec;
127
128 vecp.push_back(make_pair(time, p));
129 }
130
131 return 0;
132}
133
134Int_t PlotCurrent(TArrayD **vec, TString fname)
135{
136 Int_t rc = ReadSources(fname);
137 if (rc<0)
138 return rc;
139
140 fname += ".FEEDBACK_CALIBRATED_CURRENTS.fits";
141
142 fits file(fname.Data());
143 if (!file)
144 {
145 cerr << fname << ": " << gSystem->GetError() << endl;
146 return -2;
147 }
148
149 //cout << fname << endl;
150
151 Double_t time;
152 Float_t Imed;
153 Float_t Idev;
154 Float_t I[416];
155
156 if (!file.SetPtrAddress("Time", &time))
157 return -1;
158
159 if (!file.SetPtrAddress("I_med", &Imed))
160 return -1;
161
162 if (!file.SetPtrAddress("I_dev", &Idev))
163 return -1;
164
165 if (!file.SetPtrAddress("I", I))
166 return -1;
167
168 TGraph g1;
169 TGraph g2;
170 TGraph g3;
171 TGraph g4;
172 TGraph g5;
173 g1.SetName("CurrentsMed");
174 g2.SetName("Prediction");
175 g3.SetName("CurrentsDev");
176 g4.SetName("CurrentsMax-4");
177 g5.SetName("CurrentsMax");
178
179 while (file.GetNextRow())
180 if (Contains(vec, time))
181 {
182 // crazy pixels
183 I[66] = 0;
184 I[191] = 0;
185 I[193] = 0;
186
187 sort(I, I+320);
188
189 g1.SetPoint(g1.GetN(), time*24*3600, Imed);
190 g2.SetPoint(g2.GetN(), time*24*3600, Prediction(time));
191 g3.SetPoint(g3.GetN(), time*24*3600, Imed+Idev);
192 g4.SetPoint(g4.GetN(), time*24*3600, I[315]);
193 g5.SetPoint(g5.GetN(), time*24*3600, I[319]);
194 }
195
196 g1.SetMinimum(0);
197 g1.SetMaximum(149);
198
199 g1.SetMarkerStyle(kFullDotMedium);
200 g1.GetXaxis()->SetTimeDisplay(true);
201 g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
202 g1.GetXaxis()->SetLabelSize(0.12);
203 g1.GetYaxis()->SetLabelSize(0.1);
204 g1.GetYaxis()->SetTitle("CURRENT");
205 g1.GetYaxis()->SetTitleOffset(0.2);
206 g1.GetYaxis()->SetTitleSize(0.1);
207 g1.DrawClone("AP");
208
209 g5.SetMarkerColor(kGray);
210 g5.DrawClone("P");
211
212 g4.SetMarkerColor(kGray+1);
213 g4.DrawClone("P");
214
215 g3.SetMarkerColor(kGray+2);
216 g3.DrawClone("P");
217
218 g2.SetMarkerColor(kBlue);
219 g2.SetMarkerStyle(kFullDotMedium);
220 g2.DrawClone("P");
221
222 g1.DrawClone("P");
223
224 return 0;
225}
226
227Int_t PlotRate(TArrayD **vec, TString fname)
228{
229 fname += ".FTM_CONTROL_TRIGGER_RATES.fits";
230
231 fits file(fname.Data());
232 if (!file)
233 {
234 cerr << fname << ": " << gSystem->GetError() << endl;
235 return -2;
236 }
237
238 //cout << fname << endl;
239
240 Double_t time;
241 Float_t rate;
242 Float_t ontime, elapsed;
243
244 if (!file.SetPtrAddress("Time", &time))
245 return -1;
246
247 if (!file.SetPtrAddress("TriggerRate", &rate))
248 return -1;
249 if (!file.SetPtrAddress("OnTime", &ontime))
250 return -1;
251 if (!file.SetPtrAddress("ElapsedTime", &elapsed))
252 return -1;
253
254 TGraph g1, g2;
255 g1.SetName("TriggerRate");
256 g2.SetName("RelOnTime");
257
258 while (file.GetNextRow())
259 if (Contains(vec, time))
260 {
261 if (rate>=0)
262 {
263 g1.SetPoint(g1.GetN(), time*24*3600, rate);
264 g2.SetPoint(g2.GetN(), time*24*3600, ontime/elapsed);
265 }
266 }
267
268 g1.SetMinimum(0);
269 g1.SetMaximum(269);
270 g1.SetMarkerStyle(kFullDotMedium);
271 g1.GetXaxis()->SetTimeDisplay(true);
272 g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
273 g1.GetXaxis()->SetLabelSize(0.12);
274 g1.GetYaxis()->SetLabelSize(0.1);
275 g1.GetYaxis()->SetTitle("TRIGGER RATE");
276 g1.GetYaxis()->SetTitleOffset(0.2);
277 g1.GetYaxis()->SetTitleSize(0.1);
278 g1.DrawClone("AP");
279
280 gROOT->SetSelectedPad(0);
281 gPad->GetCanvas()->cd(4);
282
283 gPad->SetGrid();
284 gPad->SetTopMargin(0);
285 gPad->SetBottomMargin(0);
286 gPad->SetRightMargin(0.001);
287 gPad->SetLeftMargin(0.04);
288
289 g2.SetMinimum(0);
290 g2.SetMaximum(1);
291 g2.SetMarkerStyle(kFullDotMedium);
292 g2.GetXaxis()->SetTimeDisplay(true);
293 g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
294 g2.GetXaxis()->SetLabelSize(0.12);
295 g2.GetYaxis()->SetLabelSize(0.1);
296 g2.GetYaxis()->SetTitle("RELATIVE ON TIME");
297 g2.GetYaxis()->SetTitleOffset(0.2);
298 g2.GetYaxis()->SetTitleSize(0.1);
299 g2.DrawClone("AP");
300
301 return 0;
302}
303
304void PlotRateQC(UInt_t night, MSQLServer &serv)
305{
306 TString query =
307 "LEFT JOIN AnalysisResultsRunLP USING(fNight, fRunID) "
308 "WHERE fRunTypeKey=1 AND NOT ISNULL (AnalysisResultsRunLP.fNumEvtsAfterCleaning) AND fNight=";
309 query += night;
310
311 TTree *t = serv.GetTree("RunInfo", query);
312 if (!t)
313 return;
314
315 int save = gErrorIgnoreLevel;
316 gErrorIgnoreLevel = kFatal;
317
318 gROOT->SetSelectedPad(0);
319 gPad->GetCanvas()->cd(3);
320
321 t->Draw("AnalysisResultsRunLP.fNumEvtsAfterCleaning/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same");
322 TGraph *g = (TGraph*)gPad->GetPrimitive("Graph");
323 if (g)
324 {
325 g->SetName("CleaningRate");
326 g->SetMarkerColor(kRed);
327 g->SetMarkerStyle(29);//kFullDotMedium);
328 }
329
330 t->Draw("AnalysisResultsRunLP.fNumEvtsAfterQualCuts/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same");
331 g = (TGraph*)gPad->GetPrimitive("Graph");
332 if (g)
333 {
334 g->SetName("RateAfterQC");
335 g->SetMarkerColor(kBlue);
336 g->SetMarkerStyle(29);//kFullDotMedium);
337 }
338
339 gErrorIgnoreLevel = save;
340}
341
342
343Int_t PlotPointing(TArrayD **vec, TString fname)
344{
345 fname += ".DRIVE_CONTROL_POINTING_POSITION.fits";
346
347 fits file(fname.Data());
348 if (!file)
349 {
350 cerr << fname << ": " << gSystem->GetError() << endl;
351 return -2;
352 }
353
354 //cout << fname << endl;
355
356 Double_t time;
357 Double_t zd;
358
359 if (!file.SetPtrAddress("Time", &time))
360 return -1;
361 if (!file.SetPtrAddress("Zd", &zd))
362 return -1;
363
364 TGraph g;
365 g.SetName("Zd");
366
367 while (file.GetNextRow())
368 if (Contains(vec, time))
369 g.SetPoint(g.GetN(), time*24*3600, 90-zd);
370
371 g.SetMinimum(1);
372 g.SetMaximum(90);
373 g.SetMarkerStyle(kFullDotMedium);
374 g.GetXaxis()->SetTimeDisplay(true);
375 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
376 g.GetXaxis()->SetLabelSize(0.12);
377 g.GetYaxis()->SetLabelSize(0.1);
378 g.GetYaxis()->SetTitle("ELEVATION");
379 g.GetYaxis()->SetTitleOffset(0.2);
380 g.GetYaxis()->SetTitleSize(0.1);
381 g.DrawClone("AP");
382
383 return 0;
384}
385
386Int_t PlotTemperature1(TArrayD **vec, TString fname)
387{
388 fname += ".TEMPERATURE_DATA.fits";
389
390 fits file(fname.Data());
391 if (!file)
392 {
393 cerr << fname << ": " << gSystem->GetError() << endl;
394 return -2;
395 }
396
397 //cout << fname << endl;
398
399 Double_t time;
400 Float_t temp;
401
402 if (!file.SetPtrAddress("Time", &time))
403 return -1;
404 if (!file.SetPtrAddress("T", &temp))
405 return -1;
406
407 TGraph g;
408 g.SetName("ContainerTemp");
409
410 while (file.GetNextRow())
411 if (Contains(vec, time))
412 g.SetPoint(g.GetN(), time*24*3600, temp);
413
414 g.SetMinimum(-5);
415 g.SetMaximum(49);
416 g.SetMarkerStyle(kFullDotMedium);
417 g.SetMarkerColor(kRed);
418 g.GetXaxis()->SetTimeDisplay(true);
419 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
420 g.GetXaxis()->SetLabelSize(0.1);
421 g.GetYaxis()->SetLabelSize(0.1);
422 g.GetYaxis()->SetTitle("TEMP");
423 g.GetYaxis()->SetTitleOffset(0.2);
424 g.GetYaxis()->SetTitleSize(0.1);
425 g.DrawClone("AP");
426
427 return 0;
428}
429
430Int_t PlotTemperature2(TArrayD **vec, TString fname)
431{
432 fname += ".MAGIC_WEATHER_DATA.fits";
433
434 fits file(fname.Data());
435 if (!file)
436 {
437 cerr << fname << ": " << gSystem->GetError() << endl;
438 return -2;
439 }
440
441 //cout << fname << endl;
442
443 Double_t time;
444 Float_t temp;
445
446 if (!file.SetPtrAddress("Time", &time))
447 return -1;
448 if (!file.SetPtrAddress("T", &temp))
449 return -1;
450
451 TGraph g;
452 g.SetName("OutsideTemp");
453
454 while (file.GetNextRow())
455 if (Contains(vec, time))
456 g.SetPoint(g.GetN(), time*24*3600, temp);
457
458 g.SetMarkerStyle(kFullDotMedium);
459 g.SetMarkerColor(kBlue);
460 g.DrawClone("P");
461
462 return 0;
463}
464
465Int_t PlotTemperature3(TArrayD **vec, TString fname)
466{
467 fname += ".FSC_CONTROL_TEMPERATURE.fits";
468
469 fits file(fname.Data());
470 if (!file)
471 {
472 cerr << fname << ": " << gSystem->GetError() << endl;
473 return -2;
474 }
475
476 //cout << fname << endl;
477
478 Double_t time;
479 Float_t temp[31];
480
481 if (!file.SetPtrAddress("Time", &time))
482 return -1;
483 if (!file.SetPtrAddress("T_sens", temp))
484 return -1;
485
486 TGraph g, g1, g2;
487 g.SetName("SensorTempAvg");
488 g1.SetName("SensorTempMin");
489 g2.SetName("SensorTempMax");
490
491 while (file.GetNextRow())
492 if (Contains(vec, time))
493 {
494 float min = 100;
495 float max = -100;
496 double avg = 0;
497 int num = 0;
498 for (int i=0; i<31; i++)
499 if (temp[i]!=0)
500 {
501 avg += temp[i];
502 num++;
503
504 min = TMath::Min(min, temp[i]);
505 max = TMath::Max(max, temp[i]);
506 }
507
508 g.SetPoint(g.GetN(), time*24*3600, avg/num);
509 g1.SetPoint(g1.GetN(), time*24*3600, min);
510 g2.SetPoint(g2.GetN(), time*24*3600, max);
511 }
512
513 g.SetMarkerStyle(kFullDotMedium);
514 g.DrawClone("P");
515
516 /*
517 g1.SetLineWidth(1);
518 g1.DrawClone("L");
519
520 g2.SetLineWidth(1);
521 g2.DrawClone("L");
522 */
523 return 0;
524}
525
526Int_t PlotTemperature4(TArrayD **vec, TString fname)
527{
528 fname += ".FAD_CONTROL_TEMPERATURE.fits";
529
530 fits file(fname.Data());
531 if (!file)
532 {
533 cerr << fname << ": " << gSystem->GetError() << endl;
534 return -2;
535 }
536
537 //cout << fname << endl;
538
539 Double_t time;
540 Float_t temp[160];
541
542 if (!file.SetPtrAddress("Time", &time))
543 return -1;
544 if (!file.SetPtrAddress("Data1", temp) &&
545 !file.SetPtrAddress("temp", temp))
546 return -1;
547
548 Int_t num = file.GetN("temp")==0 ? file.GetN("Data1") : file.GetN("temp");
549 Int_t beg = num==82 ? 2 : 0;
550
551 TGraphErrors g1;
552 TGraph g2,g3;
553
554 g1.SetName("FadTempAvg");
555 g2.SetName("FadTempMin");
556 g3.SetName("FadTempMax");
557
558 while (file.GetNextRow())
559 if (Contains(vec, time))
560 {
561 double avg = 0;
562 double rms = 0;
563 float min = 100;
564 float max = -100;
565 for (int i=beg; i<num; i++)
566 {
567 avg += temp[i];
568 rms += temp[i]*temp[i];
569
570 min = TMath::Min(min, temp[i]);
571 max = TMath::Max(max, temp[i]);
572 }
573
574 avg /= num-beg;
575 rms /= num-beg;
576
577 g1.SetPoint(g1.GetN(), time*24*3600, avg);
578 g1.SetPointError(g1.GetN()-1, 0, sqrt(rms-avg*avg));
579
580 g2.SetPoint(g2.GetN(), time*24*3600, min);
581 g3.SetPoint(g3.GetN(), time*24*3600, max);
582 }
583
584 g1.SetLineColor(kGreen);
585 g1.DrawClone("[]");
586
587 g2.SetLineColor(kGreen);
588 g2.SetLineWidth(1);
589 g2.DrawClone("L");
590
591 g3.SetLineColor(kGreen);
592 g3.SetLineWidth(1);
593 g3.DrawClone("L");
594
595 return 0;
596}
597
598int quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, const char *outpath="quality")
599{
600 // To get correct dates in the histogram you have to add
601 // the MJDREF offset (should be 40587) and 9131.
602
603 if (y==0)
604 {
605 UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt();
606 y = nt/10000;
607 m = (nt/100)%100;
608 d = nt%100;
609
610 cout << y << "/" << m << "/" << d << endl;
611 }
612
613 TString fname=Form("/fact/aux/%04d/%02d/%02d/%04d%02d%02d", y, m, d, y, m, d);
614
615 UInt_t night = MTime(y, m, d, 0).GetNightAsInt();
616
617 MSQLMagic serv("sql.rc");
618 Bool_t con = serv.IsConnected();
619
620 cout << "quality" << endl;
621 cout << "-------" << endl;
622 cout << endl;
623 if (con)
624 {
625 cout << "Connected to " << serv.GetName() << endl;
626 cout << endl;
627 }
628 cout << "Night: " << Form("%04d-%02d-%02d", y, m, d) << endl;
629 cout << endl;
630
631 TArrayD run, beg, end;
632
633 TArrayD *runs[3] = { &run, &beg, &end };
634
635 if (con)
636 {
637 TString query;
638 query += "SELECT fRunID, fRunStart, fRunStop FROM RunInfo";
639 query += " WHERE fNight=";
640 query += night;
641 query += " AND fRunTypeKey=1 ORDER BY fRunStart";
642
643 TSQLResult *res = serv.Query(query);
644 if (!res)
645 return 1;
646
647 run.Set(res->GetRowCount());
648 beg.Set(res->GetRowCount());
649 end.Set(res->GetRowCount());
650
651 Int_t n = 0;
652
653 TSQLRow *row = 0;
654 while ((row=res->Next()))
655 {
656 run[n] = atoi((*row)[0]);
657 beg[n] = MTime((*row)[1]).GetMjd()-40587;
658 end[n] = MTime((*row)[2]).GetMjd()-40587;
659 n++;
660 delete row;
661 }
662
663 delete res;
664
665 if (n==0)
666 cout << "WARNING - No data runs in db, displaying all data." << endl;
667 else
668 cout << "Num: " << n << "\n" << endl;
669 }
670
671 TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960);
672 c->Divide(1, 6, 1e-5, 1e-5);
673
674 gROOT->SetSelectedPad(0);
675
676 c->cd(1);
677 gPad->SetGrid();
678 gPad->SetTopMargin(0);
679 gPad->SetRightMargin(0.001);
680 gPad->SetLeftMargin(0.04);
681 gPad->SetBottomMargin(0);
682 cout << PlotThresholds(runs, fname) << endl;
683
684 gROOT->SetSelectedPad(0);
685 c->cd(2);
686 gPad->SetGrid();
687 gPad->SetTopMargin(0);
688 gPad->SetRightMargin(0.001);
689 gPad->SetLeftMargin(0.04);
690 gPad->SetBottomMargin(0);
691 cout << PlotCurrent(runs, fname) << endl;
692
693 gROOT->SetSelectedPad(0);
694 c->cd(3);
695 gPad->SetGrid();
696 gPad->SetTopMargin(0);
697 gPad->SetBottomMargin(0);
698 gPad->SetRightMargin(0.001);
699 gPad->SetLeftMargin(0.04);
700 cout << PlotRate(runs, fname) << endl;
701 cout << PlotRateQC(night, serv) << endl;
702
703 gROOT->SetSelectedPad(0);
704 c->cd(5);
705 gPad->SetGrid();
706 gPad->SetTopMargin(0);
707 gPad->SetBottomMargin(0);
708 gPad->SetRightMargin(0.001);
709 gPad->SetLeftMargin(0.04);
710 cout << PlotPointing(runs, fname) << endl;
711
712 gROOT->SetSelectedPad(0);
713 c->cd(6);
714 gPad->SetGrid();
715 gPad->SetTopMargin(0);
716 gPad->SetRightMargin(0.001);
717 gPad->SetLeftMargin(0.04);
718 cout << PlotTemperature1(runs, fname) << endl;
719 cout << PlotTemperature2(runs, fname) << endl;
720 cout << PlotTemperature3(runs, fname) << endl;
721 cout << PlotTemperature4(runs, fname) << endl;
722
723 c->SaveAs(Form("%s/%04d%02d%02d.png", outpath, y, m, d));
724
725 return 0;
726}
Note: See TracBrowser for help on using the repository browser.