source: trunk/Mars/fact/plots/quality.C@ 17227

Last change on this file since 17227 was 15244, checked in by Daniela Dorner, 13 years ago
added (macro to plot quality parameter)
File size: 11.9 KB
Line 
1Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0)
2{
3 TArrayD *arr0 = vec[0];
4 TArrayD *arr1 = vec[1];
5 TArrayD *arr2 = vec[2];
6
7 for (int i=0; i<arr0->GetSize(); i++)
8 {
9 Double_t start = (*arr1)[i];
10 Double_t stop = (*arr2)[i];
11
12 if (stop>start+305./24/3600)
13 stop = start+305./24/3600;
14
15 if (t0>start-range && t0<stop+range)
16 //{
17 // if (fmod(t0, 1)>4./24 && fmod(t0,1)<4.1/24)
18 // cout << t0-start << " " <<t0 << " " << stop-t0 << " " << start-15779 << " " << stop-15779 << " " << (*arr0)[i] << endl;
19 return kTRUE;
20 //}
21 }
22
23 return arr0->GetSize()==0;
24}
25
26Int_t PlotThresholds(TArrayD **vec, TString fname)
27{
28 fname += ".RATE_CONTROL_THRESHOLD.fits";
29
30 fits file(fname.Data());
31 if (!file)
32 {
33 cerr << fname << ": " << gSystem->GetError() << endl;
34 return -2;
35 }
36
37 //cout << fname << endl;
38
39 Double_t time;
40 UShort_t th;
41
42 if (!file.SetPtrAddress("Time", &time))
43 return -1;
44
45 if (!file.SetPtrAddress("threshold", &th))
46 return -1;
47
48 TGraph g;
49
50 while (file.GetNextRow())
51 {
52 if (Contains(vec, time, 10./(24*3600)))
53 g.SetPoint(g.GetN(), time*24*3600, th);
54 }
55
56 g.SetMarkerStyle(kFullDotMedium);
57 g.GetXaxis()->SetTimeDisplay(true);
58 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
59 g.GetXaxis()->SetLabelSize(0.12);
60 g.GetYaxis()->SetLabelSize(0.1);
61 g.GetYaxis()->SetTitle("THRESHOLD");
62 g.GetYaxis()->SetTitleOffset(0.2);
63 g.GetYaxis()->SetTitleSize(0.1);
64 g.DrawClone("AP");
65
66 return 0;
67}
68
69#include <vector>
70#include <pair>
71
72vector<pair<double, ln_equ_posn>> vec;
73
74ln_equ_posn FindPointing(Double_t time)
75{
76 for (int i=0; i<vec.size(); i++)
77 if (time<vec[i].first)
78 return i==0 ? ln_equ_posn() : vec[i-1].second;
79
80 return vec[vec.size()-1].second;
81}
82
83Float_t Prediction(Double_t time)
84{
85 Double_t jd = time + 40587 + 2400000.5;
86
87 Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
88
89 // Nova::EquPosn pos = FindPointing(time);
90
91 // get local position of moon
92 Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
93
94 // Distance between source and moon
95 //const double angle =
96 // MAstro::AngularDistance(TMath::DegToRad()*(90-moon.dec),
97 // TMath::DegToRad()*moon.ra,
98 // TMath::DegToRad()*(90-pos.dec),
99 // TMath::DegToRad()*pos.ra);
100
101 // Current prediction
102 //double cang = sin(angle);
103 double calt = sin(hrzm.alt*TMath::DegToRad());
104
105 double disk = Nova::GetLunarDisk(jd);
106
107 double lc = /*cang==0 ? 0 :*/ calt*pow(disk, 2.5);///sqrt(cang);
108
109 return 5+103*lc>4.5 ? 5+103*lc : 4.5;
110}
111
112Int_t ReadSources(TString fname)
113{
114 fname += ".DRIVE_CONTROL_SOURCE_POSITION.fits";
115
116 fits file(fname.Data());
117 if (!file)
118 {
119 cerr << fname << ": " << gSystem->GetError() << endl;
120 return -2;
121 }
122
123 Double_t time, ra, dec;
124 if (!file.SetPtrAddress("Time", &time))
125 return -1;
126 if (!file.SetPtrAddress("Ra_cmd", &ra))
127 return -1;
128 if (!file.SetPtrAddress("Dec_cmd", &dec))
129 return -1;
130
131 while (file.GetNextRow())
132 {
133 ln_equ_posn p;
134 p.ra = ra*15;
135 p.dec = dec;
136
137 vec.push_back(make_pair(time, p));
138 }
139
140 return 0;
141}
142
143Int_t PlotCurrent(TArrayD **vec, TString fname)
144{
145 Int_t rc = ReadSources(fname);
146 if (rc<0)
147 return rc;
148
149 fname += ".FEEDBACK_CALIBRATED_CURRENTS.fits";
150
151 fits file(fname.Data());
152 if (!file)
153 {
154 cerr << fname << ": " << gSystem->GetError() << endl;
155 return -2;
156 }
157
158 //cout << fname << endl;
159
160 Double_t time;
161 Float_t Imed;
162
163 if (!file.SetPtrAddress("Time", &time))
164 return -1;
165
166 if (!file.SetPtrAddress("I_med", &Imed))
167 return -1;
168
169 TGraph g1;
170 TGraph g2;
171
172 while (file.GetNextRow())
173 if (Contains(vec, time))
174 {
175 g1.SetPoint(g1.GetN(), time*24*3600, Imed);
176 g2.SetPoint(g2.GetN(), time*24*3600, Prediction(time));
177 }
178
179 g1.SetMinimum(0);
180 g1.SetMaximum(99);
181
182 g1.SetMarkerStyle(kFullDotMedium);
183 g1.GetXaxis()->SetTimeDisplay(true);
184 g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
185 g1.GetXaxis()->SetLabelSize(0.12);
186 g1.GetYaxis()->SetLabelSize(0.1);
187 g1.GetYaxis()->SetTitle("CURRENT");
188 g1.GetYaxis()->SetTitleOffset(0.2);
189 g1.GetYaxis()->SetTitleSize(0.1);
190 g1.DrawClone("AP");
191
192 g2.SetMarkerColor(kBlue);
193 g2.SetMarkerStyle(kFullDotMedium);
194 g2.GetXaxis()->SetTimeDisplay(true);
195 g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
196 g2.GetXaxis()->SetLabelSize(0.12);
197 g2.GetYaxis()->SetLabelSize(0.1);
198 g2.GetYaxis()->SetTitle("CURRENT");
199 g2.GetYaxis()->SetTitleOffset(0.2);
200 g2.GetYaxis()->SetTitleSize(0.1);
201 g2.DrawClone("P");
202
203 g1.DrawClone("P");
204
205 return 0;
206}
207
208Int_t PlotRate(TArrayD **vec, TString fname)
209{
210 fname += ".FTM_CONTROL_TRIGGER_RATES.fits";
211
212 fits file(fname.Data());
213 if (!file)
214 {
215 cerr << fname << ": " << gSystem->GetError() << endl;
216 return -2;
217 }
218
219 //cout << fname << endl;
220
221 Double_t time;
222 Float_t rate;
223 Float_t ontime, elapsed;
224
225 if (!file.SetPtrAddress("Time", &time))
226 return -1;
227
228 if (!file.SetPtrAddress("TriggerRate", &rate))
229 return -1;
230 if (!file.SetPtrAddress("OnTime", &ontime))
231 return -1;
232 if (!file.SetPtrAddress("ElapsedTime", &elapsed))
233 return -1;
234
235 TGraph g1, g2;
236
237 while (file.GetNextRow())
238 if (Contains(vec, time))
239 {
240 if (rate>=0)
241 {
242 g1.SetPoint(g1.GetN(), time*24*3600, rate);
243 g2.SetPoint(g2.GetN(), time*24*3600, ontime/elapsed);
244 }
245 }
246
247 g1.SetMinimum(0);
248 g1.SetMaximum(200);
249 g1.SetMarkerStyle(kFullDotMedium);
250 g1.GetXaxis()->SetTimeDisplay(true);
251 g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
252 g1.GetXaxis()->SetLabelSize(0.12);
253 g1.GetYaxis()->SetLabelSize(0.1);
254 g1.GetYaxis()->SetTitle("TRIGGER RATE");
255 g1.GetYaxis()->SetTitleOffset(0.2);
256 g1.GetYaxis()->SetTitleSize(0.1);
257 g1.DrawClone("AP");
258
259 gROOT->SetSelectedPad(0);
260 gPad->GetCanvas()->cd(4);
261
262 gPad->SetGridx();
263 gPad->SetTopMargin(0);
264 gPad->SetBottomMargin(0);
265 gPad->SetRightMargin(0.001);
266 gPad->SetLeftMargin(0.04);
267
268 g2.SetMinimum(0);
269 g2.SetMaximum(1);
270 g2.SetMarkerStyle(kFullDotMedium);
271 g2.GetXaxis()->SetTimeDisplay(true);
272 g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
273 g2.GetXaxis()->SetLabelSize(0.12);
274 g2.GetYaxis()->SetLabelSize(0.1);
275 g2.GetYaxis()->SetTitle("RELATIVE ON TIME");
276 g2.GetYaxis()->SetTitleOffset(0.2);
277 g2.GetYaxis()->SetTitleSize(0.1);
278 g2.DrawClone("AP");
279
280
281 return 0;
282}
283
284Int_t PlotPointing(TArrayD **vec, TString fname)
285{
286 fname += ".DRIVE_CONTROL_POINTING_POSITION.fits";
287
288 fits file(fname.Data());
289 if (!file)
290 {
291 cerr << fname << ": " << gSystem->GetError() << endl;
292 return -2;
293 }
294
295 //cout << fname << endl;
296
297 Double_t time;
298 Double_t zd;
299
300 if (!file.SetPtrAddress("Time", &time))
301 return -1;
302 if (!file.SetPtrAddress("Zd", &zd))
303 return -1;
304
305 TGraph g;
306
307 while (file.GetNextRow())
308 if (Contains(vec, time))
309 g.SetPoint(g.GetN(), time*24*3600, 90-zd);
310
311 g.SetMinimum(0);
312 g.SetMaximum(90);
313 g.SetMarkerStyle(kFullDotMedium);
314 g.GetXaxis()->SetTimeDisplay(true);
315 g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
316 g.GetXaxis()->SetLabelSize(0.1);
317 g.GetYaxis()->SetLabelSize(0.1);
318 g.GetYaxis()->SetTitle("ELEVATION");
319 g.GetYaxis()->SetTitleOffset(0.2);
320 g.GetYaxis()->SetTitleSize(0.1);
321 g.DrawClone("AP");
322
323 return 0;
324}
325/*
326Int_t PlotMem(TArrayD **vec, TString fname)
327{
328 fname += ".FAD_CONTROL_STATISTICS1.fits";
329
330 fits file(fname.Data());
331 if (!file)
332 {
333 cerr << fname << ": " << gSystem->GetError() << endl;
334 return -2;
335 }
336
337 //cout << fname << endl;
338
339 Double_t time;
340 UInt_t buf[5];
341 ULong64_t mem[4];
342 UInt_t rate[2];
343
344 if (!file.SetPtrAddress("Time", &time))
345 return -1;
346 if (!file.SetPtrAddress("bufferInfo", buf) ||
347 !file.SetPtrAddress("memInfo", mem) ||
348 !file.SetPtrAddress("rateNew", rate))
349 return -1;
350
351 TGraph g;
352
353 while (file.GetNextRow())
354 if (Contains(vec, time))
355 g.SetPoint(g.GetN(), time*24*3600, mem[1]);
356
357 g.SetMarkerStyle(kFullDotMedium);
358 g.GetXaxis()->SetTimeDisplay(true);
359 g.GetXaxis()->SetTimeFormat("%H:%M");
360 g.GetXaxis()->SetTimeOffset(0, "gmt");
361 g.GetXaxis()->SetLabelSize(0.12);
362 g.GetYaxis()->SetLabelSize(0.1);
363 //g.GetYaxis()->SetTitle("ELEVATION");
364 g.GetYaxis()->SetTitleOffset(0.2);
365 g.GetYaxis()->SetTitleSize(0.1);
366 g.DrawClone("AP");
367
368 return 0;
369}
370*/
371void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0)
372{
373 // To get correct dates in the histogram you have to add
374 // the MJDREF offset (should be 40587) and 9131.
375
376 if (y==0)
377 {
378 UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt();
379 y = nt/10000;
380 m = (nt/100)%100;
381 d = nt%100;
382
383 cout << y << "/" << m << "/" << d << endl;
384 }
385
386 TString fname=Form("/fact/aux/%04d/%02d/%02d/%04d%02d%02d", y, m, d, y, m, d);
387
388 UInt_t night = MTime(y, m, d, 0).GetNightAsInt();
389
390 MSQLMagic serv("sql.rc");
391 Bool_t con = serv.IsConnected();
392
393 cout << "quality" << endl;
394 cout << "-------" << endl;
395 cout << endl;
396 if (con)
397 {
398 cout << "Connected to " << serv.GetName() << endl;
399 cout << endl;
400 }
401 cout << "Night: " << Form("%04d-%02d-%02d", y, m, d) << endl;
402 cout << endl;
403
404 TArrayD run, beg, end;
405
406 TArrayD *runs[3] = { &run, &beg, &end };
407
408 if (con)
409 {
410 TString query;
411 query += "SELECT fRunID, fRunStart, fRunStop FROM runinfo";
412 query += " WHERE fNight=";
413 query += night;
414 query += " AND fRunTypeKEY=1 ORDER BY fRunStart";
415
416 TSQLResult *res = serv.Query(query);
417 if (!res)
418 return kFALSE;
419
420 run.Set(res->GetRowCount());
421 beg.Set(res->GetRowCount());
422 end.Set(res->GetRowCount());
423
424 Int_t n = 0;
425
426 TSQLRow *row = 0;
427 while ((row=res->Next()))
428 {
429 run[n] = atoi((*row)[0]);
430 beg[n] = MTime((*row)[1]).GetMjd()-40587;
431 end[n] = MTime((*row)[2]).GetMjd()-40587;
432 n++;
433 delete row;
434 }
435
436 delete res;
437
438 if (n==0)
439 cout << "WARNING - No data runs in db, displaying all data." << endl;
440 }
441
442 TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960);
443 c->Divide(1, 5, 1e-5, 1e-5);
444
445 gROOT->SetSelectedPad(0);
446 c->cd(1);
447 gPad->SetGridx();
448 gPad->SetTopMargin(0);
449 gPad->SetRightMargin(0.001);
450 gPad->SetLeftMargin(0.04);
451 gPad->SetBottomMargin(0);
452 cout << PlotThresholds(runs, fname) << endl;
453
454 gROOT->SetSelectedPad(0);
455 c->cd(2);
456 gPad->SetGridx();
457 gPad->SetTopMargin(0);
458 gPad->SetRightMargin(0.001);
459 gPad->SetLeftMargin(0.04);
460 gPad->SetBottomMargin(0);
461 cout << PlotCurrent(runs, fname) << endl;
462
463 gROOT->SetSelectedPad(0);
464 c->cd(3);
465 gPad->SetGridx();
466 gPad->SetTopMargin(0);
467 gPad->SetBottomMargin(0);
468 gPad->SetRightMargin(0.001);
469 gPad->SetLeftMargin(0.04);
470 cout << PlotRate(runs, fname) << endl;
471
472 gROOT->SetSelectedPad(0);
473 c->cd(5);
474 gPad->SetGridx();
475 gPad->SetTopMargin(0);
476 gPad->SetRightMargin(0.001);
477 gPad->SetLeftMargin(0.04);
478 cout << PlotPointing(runs, fname) << endl;
479/*
480 gROOT->SetSelectedPad(0);
481 c->cd(6);
482 gPad->SetGridx();
483 gPad->SetTopMargin(0);
484 gPad->SetRightMargin(0.001);
485 gPad->SetLeftMargin(0.04);
486 cout << PlotMem(runs, fname) << endl;
487*/
488 c->SaveAs(Form("%04d%02d%02d-quality.pdf", y, m, d));
489}
490
491//20130314_141
Note: See TracBrowser for help on using the repository browser.