source: trunk/MagicSoft/Mars/datacenter/macros/plotoptical.C@ 8167

Last change on this file since 8167 was 8140, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 16.9 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: plotoptical.C,v 1.2 2006-10-20 18:26:00 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Thomas Bretz, 05/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
21! Author(s): Daniela Dorner, 05/2005 <mailto:dorner@astro.uni-wuerzburg.de>
22!
23! Copyright: MAGIC Software Development, 2000-2006
24!
25!
26\* ======================================================================== */
27
28/////////////////////////////////////////////////////////////////////////////
29//
30// plotdb.C
31// ========
32//
33// This macro is used to read quality parameters from the DB and plot them.
34//
35// The parameters are from the following files:
36// calib*.root:mean conversion factor, mean arrival time, rms arrival time
37// (each parameter for inner and outer camera)
38// signal*.root: mean pedestal rms (for inner and outer camera)
39// star*.root: PSF, # of Muons, Effective OnTime, Muon rate,
40// Ratio MC/Data(MuonSize) and mean number of islands
41//
42// In the DB these values are stored in the tables Calibration and Star.
43//
44// Usage:
45// .x plotdb.C --> all values in the DB are plotted
46// You can chose are certain period:
47// .x plotdb.C(25) --> all values from period 25 are plotted
48// or a time period from a certain date to a certain date
49// .x plotdb.C("2004-11-14 00:00:00", "2005-02-28 00:00:00")
50// --> all values from 14.11.2004 0h to 28.2.2005 0h are plotted
51// or all data, but with dataset data highlighted
52// .x plotdb.C("dataset.txt")
53// --> the sequences defined in dataset.txt are highlighted (blue:on, red:off)
54// --> You can also add a dataset-name as last argument to one of the
55// calls above
56//
57// Make sure, that database and password are corretly set in a resource
58// file called sql.rc and the resource file is found.
59//
60/////////////////////////////////////////////////////////////////////////////
61#include <iostream>
62#include <iomanip>
63
64#include <TH1.h>
65#include <TEnv.h>
66#include <TPad.h>
67#include <TLine.h>
68#include <TText.h>
69#include <TFrame.h>
70#include <TStyle.h>
71#include <TGraph.h>
72#include <TCanvas.h>
73#include <TSQLRow.h>
74#include <TSQLResult.h>
75
76#include "MTime.h"
77#include "MAstro.h"
78#include "MDataSet.h"
79#include "MSQLMagic.h"
80#include "MStatusDisplay.h"
81
82class MPlot : public MParContainer
83{
84private:
85 MSQLMagic &fServer;
86
87 MDataSet *fDataSet;
88
89 TString fRequestFrom;
90 TString fRequestTo;
91 Int_t fRequestPeriod;
92
93 Float_t fPlotMin;
94 Float_t fPlotMax;
95
96 Float_t fHistMin;
97 Float_t fHistMax;
98
99 TString fDescription;
100 TString fNameTab;
101
102 TString fCondition;
103 Bool_t fGroupBy;
104
105 void PlotTable(TSQLResult &res, TString name, Float_t fmin, Float_t fmax, Float_t resolution)
106 {
107 gStyle->SetOptStat(111111);
108
109 TSQLRow *row;
110
111 TGraph gt;
112 gt.SetNameTitle(name, Form("%s vs Time", name.Data()));
113 gt.SetMarkerStyle(kFullDotMedium);
114
115 TGraph gz;
116 gz.SetNameTitle(name, Form("%s vs <Zd>", name.Data()));
117 gz.SetMarkerStyle(kFullDotMedium);
118
119 TGraph gt0, gt1;
120 gt0.SetMarkerColor(kRed);
121 gt1.SetMarkerColor(kBlue);
122 gt0.SetMarkerStyle(kFullDotLarge);
123 gt1.SetMarkerStyle(kFullDotLarge);
124
125 TGraph gz0, gz1;
126 gz0.SetMarkerColor(kRed);
127 gz1.SetMarkerColor(kBlue);
128 gz0.SetMarkerStyle(kFullDotLarge);
129 gz1.SetMarkerStyle(kFullDotLarge);
130
131 Int_t first = -1;
132 Int_t last = -1;
133
134 while ((row=res.Next()))
135 {
136 const char *date = (*row)[0];
137 const char *zd = (*row)[1];
138 const char *val = (*row)[2];
139 const char *snum = res.GetFieldCount()>3 ? (*row)[3] : 0;
140 if (!date || !val || !zd)
141 continue;
142
143 MTime t(date);
144 if (!t.SetSqlDateTime(date))
145 continue;
146
147 if (fRequestPeriod>0 && MAstro::GetMagicPeriod(t.GetMjd())!=fRequestPeriod)
148 continue;
149
150 if (first<0)
151 first = TMath::Nint(TMath::Floor(t.GetMjd()));
152 last = TMath::Nint(TMath::Ceil(t.GetMjd()));
153
154 UInt_t seq = snum ? atoi(snum) : 0;
155
156 Float_t value = atof(val);
157 Float_t zenith = atof(zd);
158
159 if (fDataSet)
160 {
161 if (fDataSet->HasOnSequence(seq))
162 {
163 gt1.SetPoint(gt1.GetN(), t.GetAxisTime(), value);
164 gz1.SetPoint(gz1.GetN(), zenith, value);
165 }
166
167 if (fDataSet->HasOffSequence(seq))
168 {
169 gt0.SetPoint(gt0.GetN(), t.GetAxisTime(), value);
170 gz0.SetPoint(gz0.GetN(), zenith, value);
171 }
172 }
173
174 gt.SetPoint(gt.GetN(), t.GetAxisTime(), value);
175 gz.SetPoint(gz.GetN(), zenith, value);
176 }
177
178 TString title = fNameTab.IsNull() ? name(name.First('.')+2, name.Length()) : fNameTab;
179 cerr << setprecision(4) << setw(10) << title << ": ";
180 if (gt.GetN()==0)
181 {
182 cerr << " <empty>" << endl;
183 return;
184 }
185 cerr << setw(8) << gt.GetMean(2) << "+-" << setw(8) << gt.GetRMS(2) << " ";
186 if (gt0.GetN()>0 || gt1.GetN()>0)
187 {
188 cerr << setw(8) << gt1.GetMean(2) << "+-" << setw(8) << gt1.GetRMS(2) << " ";
189 cerr << setw(8) << gt0.GetMean(2) << "+-" << setw(8) << gt0.GetRMS(2);
190 }
191 cerr << endl;
192
193 // If this is done earlier the plots remain empty since root 5.12/00
194 if (fmax>fmin)
195 {
196 gt.SetMinimum(fmin);
197 gt.SetMaximum(fmax);
198 gz.SetMinimum(fmin);
199 gz.SetMaximum(fmax);
200 }
201
202 gROOT->SetSelectedPad(0);
203
204 TCanvas &c = fDisplay ? fDisplay->AddTab(title) : *new TCanvas;
205 c.SetFillColor(kWhite);
206 c.SetBorderMode(0);
207 c.Divide(1,2);
208
209 TVirtualPad *pad = gPad;
210 pad->cd(2);
211 gPad->SetBorderMode(0);
212 gPad->SetFrameBorderMode(0);
213 gPad->SetGridy();
214
215 gPad->SetLeftMargin(0.06);
216 gPad->SetRightMargin(0.06);
217 gPad->SetBottomMargin(0.08);
218
219 TH1 *h = gt.GetHistogram();
220
221 h->SetXTitle("Time");
222 h->SetYTitle(name);
223 h->GetXaxis()->SetTimeDisplay(1);
224 h->GetYaxis()->SetTitleOffset(0.8);
225 h->GetXaxis()->SetTitleOffset(1.0);
226 h->GetXaxis()->SetLabelOffset(0.01);
227
228 gt.DrawClone("AP");
229 if (gt0.GetN()>0)
230 gt0.DrawClone("P");
231 if (gt1.GetN()>0)
232 gt1.DrawClone("P");
233
234 TLine l;
235 TText t;
236 Int_t num=0;
237 l.SetLineStyle(kDotted);
238 l.SetLineColor(kBlue);
239 t.SetTextColor(kBlue);
240 l.SetLineWidth(1);
241 t.SetTextSize(h->GetXaxis()->GetLabelSize());
242 t.SetTextAlign(21);
243 Int_t p0 = MAstro::GetMagicPeriod(first);
244 for (Int_t p = first; p<last; p++)
245 {
246 Int_t p1 = MAstro::GetMagicPeriod(p);
247 if (p1!=p0)
248 {
249 l.DrawLine(MTime(p).GetAxisTime(), h->GetMinimum(), MTime(p).GetAxisTime(), h->GetMaximum());
250 t.DrawText(MTime(p+15).GetAxisTime(), h->GetMaximum(), Form("%d", p1));
251 num++;
252 }
253 p0 = p1;
254 }
255 if (num<4)
256 gPad->SetGridx();
257
258 const Double_t min = fHistMin>fHistMax ? h->GetMinimum()-resolution/2 : fHistMin;
259 const Double_t max = fHistMin>fHistMax ? h->GetMaximum()+resolution/2 : fHistMax;
260
261 // Use this to save the pad with the time development to a file
262 //gPad->SaveAs(Form("plotdb-%s.eps", title.Data()));
263
264 pad->cd(1);
265 gPad->SetBorderMode(0);
266 gPad->SetFrameBorderMode(0);
267 gPad->Divide(2,1);
268
269 TVirtualPad *pad2 = gPad;
270 pad2->cd(1);
271 gPad->SetBorderMode(0);
272 gPad->SetFrameBorderMode(0);
273 gPad->SetGridx();
274 gPad->SetGridy();
275
276 const Int_t n = resolution>0 ? TMath::Nint((max-min)/resolution) : 50;
277
278 TH1F hist("Hist", Form("Distribution of %s", fDescription.IsNull() ? name.Data() : fDescription.Data()), n, min, max);
279 hist.SetDirectory(0);
280
281 for (int i=0; i<gt.GetN(); i++)
282 hist.Fill(gt.GetY()[i]);
283
284 if (fDescription.IsNull())
285 hist.SetXTitle(name);
286 hist.SetYTitle("Counts");
287
288 hist.DrawCopy("");
289
290 pad2->cd(2);
291 gPad->SetBorderMode(0);
292 gPad->SetFrameBorderMode(0);
293 gPad->SetGridy();
294
295 TH1 *h2 = gz.GetHistogram();
296
297 h2->SetXTitle("Zd");
298 h2->SetYTitle(name);
299
300 gz.DrawClone("AP");
301
302 if (gz0.GetN()>0)
303 gz0.DrawClone("P");
304 if (gz1.GetN()>0)
305 gz1.DrawClone("P");
306 }
307
308public:
309 MPlot(MSQLMagic &server) : fServer(server), fDataSet(NULL),
310 fRequestPeriod(-1), fPlotMin(0), fPlotMax(-1), fHistMin(0), fHistMax(-1), fGroupBy(kFALSE)
311 {
312 }
313 ~MPlot()
314 {
315 if (fDataSet)
316 delete fDataSet;
317 }
318 void SetDataSet(const TString filename)
319 {
320 if (fDataSet)
321 {
322 delete fDataSet;
323 fDataSet = NULL;
324 }
325 if (!filename.IsNull())
326 fDataSet = new MDataSet(filename);
327 }
328 void SetPlotRange(Float_t min, Float_t max, Int_t n=5)
329 { fPlotMin = min; fPlotMax = max; }
330 void SetHistRange(Float_t min, Float_t max)
331 { fHistMin = min; fHistMax = max; }
332 void SetRequestRange(const char *from="", const char *to="")
333 { fRequestFrom = from; fRequestTo = to; }
334 void SetRequestPeriod(Int_t n=-1)
335 { fRequestPeriod = n; }
336 void SetCondition(const char *cond="")
337 { fCondition = cond; }
338 void SetDescription(const char *d, const char *t=0) { fDescription = d; fNameTab = t; }
339 void EnableGroupBy(Bool_t b=kTRUE) { fGroupBy=b; }
340
341 Int_t QueryKeyOfSource(TString src)
342 {
343 return fServer.QueryKeyOfName("Object", src, kFALSE);
344 }
345
346 Bool_t Plot(const char *value, Float_t min=0, Float_t max=-1, Float_t resolution=0)
347 {
348 TString named = "OpticalData.fTimeStamp";
349 TString named2 = fGroupBy ? "AVG(fZenithDistance)" : "fZenithDistance";
350 TString namev = value;
351 TString join = "fSequenceFirst";
352
353 TString tablev = namev(0, namev.First('.'));
354 TString valuev = namev(namev.First('.')+1, namev.Length());
355
356 TString tabled = named(0, named.First('.'));
357 TString valued = named(named.First('.')+1, named.Length());
358
359 TString query;
360 query = Form("select %s, %s, %s ", valued.Data(), named2.Data(), value);
361 query += Form("from %s ", tabled.Data());
362
363 //const Bool_t interval = !fRequestFrom.IsNull() && !fRequestTo.IsNull();
364
365 if (TString(value).Contains("Object."))
366 {
367 query += "left join Object on Object.fObjectKEY=OpticalData.fObjectKEY ";
368 }
369
370 if (!fCondition.IsNull())
371 {
372 query += "where ";
373 query += fCondition;
374 query += " ";
375 }
376
377 if (fGroupBy)
378 query += " GROUP BY fTimeStamp ";
379
380 query += "order by fTimeStamp";
381
382 TSQLResult *res = fServer.Query(query);
383 if (!res)
384 return kFALSE;
385
386 if (max>min)
387 PlotTable(*res, namev, min, max, resolution);
388 else
389 PlotTable(*res, namev, fPlotMin, fPlotMax, resolution);
390
391 delete res;
392 return kTRUE;
393 }
394};
395
396void plotall(MPlot &plot, TString source)
397{
398 TString cond = "fStatusKEY=13";
399 if (!source.IsNull())
400 {
401 const Int_t key = plot.QueryKeyOfSource(source);
402 if (key<0)
403 return;
404 cond += Form(" and fObjectKEY=%d", key);
405
406 }
407 plot.SetCondition(cond);
408
409 plot.SetDescription("Exposure;T_{E} [s]", "Expo");
410 plot.Plot("OpticalData.fExposure", 0, 900, 60);
411
412 plot.SetDescription("Sky Level;B [s^{-1}]", "SkyLvl");
413 plot.Plot("OpticalData.fSkyLevel/OpticalData.fExposure", 0, 5.0, 0.01);
414
415 plot.SetDescription("Full Width Half Maximum;FWHM [s^{-1}]", "Fwhm");
416 plot.Plot("OpticalData.fFWHM/OpticalData.fExposure", 0, 0.05, 0.001);
417
418 plot.SetDescription("Aperture Radius;R_{A}", "ApRad");
419 plot.Plot("OpticalData.fApertureRadius", 0, 10, 1);
420/*
421 plot.SetDescription("Instrumental Magnitude;M_{I}\\cdot s^{-1}", "InstMag/s");
422 plot.Plot("OpticalData.fInstrumentalMag/OpticalData.fExposure", 0, 0.2, 0.005);
423
424 plot.SetDescription("Instrumental Magnitude Error;\\sigma_{M}\\cdot s^{-1}", "MagErr/s");
425 plot.Plot("OpticalData.fInstrumentalMagErr/OpticalData.fExposure", 0, 0.01, 0.0002);
426
427 plot.SetDescription("Instrumental Magnitude;M_{I}", "InstMag");
428 plot.Plot("OpticalData.fInstrumentalMag", 0, 30, 0.5);
429
430 plot.SetDescription("Instrumental Magnitude Error;\\sigma_{M}", "MagErr");
431 plot.Plot("OpticalData.fInstrumentalMagErr", 0, 1, 0.01);
432 */
433 plot.SetDescription("m_{1};m_{1}", "M1");
434 plot.Plot("OpticalData.fInstrumentalMag+2.5*log10(OpticalData.fExposure)", 10, 35, 0.2);
435
436 cond += " and Object.fObjectName not like '%/BL' and not IsNull(Object.fMagnitude) ";
437 plot.SetCondition(cond);
438
439 TString ext("3080/25.0*pow(10, (OpticalData.fInstrumentalMag+2.5*log10(OpticalData.fExposure)-Object.fMagnitude)/-2.5)");
440 ext += "+0.0028*fZenithDistance-0.08";
441/*
442 plot.SetDescription("m_{1}-m_{true} (Extinction per Object);m_{1}-m_{true}", "ExtObj");
443 plot.Plot(ext, 0.5, 1.2, 0.01);
444 */
445
446 plot.EnableGroupBy();
447 plot.SetDescription("m_{1}-m_{true} (Extinction per Image);m_{1}-m_{true}", "ExtImg");
448 plot.Plot(Form("AVG(%s)", ext.Data()), 0.05, 1.2, 0.01);
449
450 plot.SetDescription("Error m_{1}-m_{true} (Extinction per Image);ERR m_{1}-m_{true}", "ExtImgErr");
451 plot.Plot(Form("STD(%s)", ext.Data()), 0, 0.3, 0.005);
452
453 plot.SetDescription("m_{1}-m_{true} (Extinction per Hour);m_{1}-m_{true}", "ExtHour");
454 plot.Plot(Form("AVG(%s), date_format(fTimeStamp, '%%Y-%%m-%%d %%H') as fTimeStamp", ext.Data()),
455 0.5, 1.2, 0.01);
456
457 plot.SetDescription("m_{1}-m_{true} (Extinction per Night);m_{1}-m_{true}", "ExtNight");
458 plot.Plot(Form("AVG(%s), date_format(adddate(fTimeStamp,Interval 12 hour),'%%Y-%%m-%%d') as fTimeStamp", ext.Data()),
459 0.5, 1.2, 0.01);
460}
461
462int plotoptical(TString from, TString to, const char *source=0)
463{
464 TEnv env("sql.rc");
465
466 MSQLMagic serv(env);
467 if (!serv.IsConnected())
468 {
469 cout << "ERROR - Connection to database failed." << endl;
470 return 0;
471 }
472
473 cout << "plotoptical" << endl;
474 cout << "-----------" << endl;
475 cout << endl;
476 cout << "Connected to " << serv.GetName() << endl;
477 cout << endl;
478
479 MStatusDisplay *d = new MStatusDisplay;
480 d->SetWindowName(serv.GetName());
481 d->SetTitle(serv.GetName());
482
483 MPlot plot(serv);
484// plot.SetDataSet(dataset);
485 plot.SetDisplay(d);
486 plot.SetRequestRange(from, to);
487 plotall(plot, source);
488 d->SaveAsRoot("plotoptical.root");
489 d->SaveAsPS("plotoptical.ps");
490
491 return 1;
492}
493
494int plotoptical(const char *source)
495{
496 TEnv env("sql.rc");
497
498 MSQLMagic serv(env);
499 if (!serv.IsConnected())
500 {
501 cout << "ERROR - Connection to database failed." << endl;
502 return 0;
503 }
504
505 cout << "plotoptical" << endl;
506 cout << "-----------" << endl;
507 cout << endl;
508 cout << "Connected to " << serv.GetName() << endl;
509 cout << endl;
510
511 MStatusDisplay *d = new MStatusDisplay;
512 d->SetWindowName(serv.GetName());
513 d->SetTitle(serv.GetName());
514
515 MPlot plot(serv);
516// plot.SetDataSet(ds);
517 plot.SetDisplay(d);
518 plot.SetRequestRange("", "");
519 plotall(plot, source);
520 d->SaveAsRoot("plotoptical.root");
521 d->SaveAsPS("plotoptical.ps");
522
523 return 1;
524}
525
526int plotoptical(Int_t period, const char *source="")
527{
528 TEnv env("sql.rc");
529
530 MSQLMagic serv(env);
531 if (!serv.IsConnected())
532 {
533 cout << "ERROR - Connection to database failed." << endl;
534 return 0;
535 }
536
537 cout << "plotoptical" << endl;
538 cout << "-----------" << endl;
539 cout << endl;
540 cout << "Connected to " << serv.GetName() << endl;
541 cout << endl;
542
543 MStatusDisplay *d = new MStatusDisplay;
544 d->SetWindowName(serv.GetName());
545 d->SetTitle(serv.GetName());
546
547 MPlot plot(serv);
548// plot.SetDataSet(dataset);
549 plot.SetDisplay(d);
550 plot.SetRequestPeriod(period);
551 plotall(plot, source);
552 d->SaveAsRoot("plotoptical.root");
553 d->SaveAsPS("plotoptical.ps");
554
555 return 1;
556}
557
558int plotoptical()
559{
560 return plotoptical("", "");
561}
Note: See TracBrowser for help on using the repository browser.