Context Navigation

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

Visit:

Last change on this file since 8198 was 8186, checked in by tbretz, 19 years ago
* empty log message *
File size: 19.2 KB

Line
1	/* ======================================================================== *\
2	! $Name: not supported by cvs2svn $:$Id: plotoptical.C,v 1.4 2006-11-01 08:53:26 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 <TCanvas.h>
72	#include <TPRegexp.h>
73	#include <TSQLRow.h>
74	#include <TSQLResult.h>
75	#include <TGraphErrors.h>
76
77	#include "MTime.h"
78	#include "MAstro.h"
79	#include "MDataSet.h"
80	#include "MSQLMagic.h"
81	#include "MStatusDisplay.h"
82
83	class MPlot : public MParContainer
84	{
85	public:
86	enum GroupBy_t
87	{
88	kNone,
89	kGroupByPrimary,
90	kGroupByHour,
91	kGroupByNight,
92	kGroupByWeek,
93	kGroupByMonth,
94	kGroupByYear
95	};
96	private:
97	MSQLMagic &fServer;
98
99	MDataSet *fDataSet;
100
101	TString fPrimaryDate;
102	TString fPrimaryNumber;
103	TString fSecondary;
104
105	TString fRequestFrom;
106	TString fRequestTo;
107	Int_t fRequestPeriod;
108
109	Float_t fPlotMin;
110	Float_t fPlotMax;
111
112	Float_t fHistMin;
113	Float_t fHistMax;
114
115	TString fDescription;
116	TString fNameTab;
117
118	TString fCondition;
119	GroupBy_t fGroupBy;
120
121	void PlotTable(TSQLResult &res, TString name, Float_t fmin, Float_t fmax, Float_t resolution)
122	{
123	gStyle->SetOptStat(111111);
124
125	TSQLRow *row;
126
127	TGraph &gt = res.GetFieldCount()>4 ? new TGraphErrors : new TGraph;
128	gt.SetNameTitle(name, Form("%s vs Time", name.Data()));
129	gt.SetMarkerStyle(kFullDotMedium);
130
131	TGraph gz;
132	gz.SetNameTitle(name, Form("%s vs <Zd>", name.Data()));
133	gz.SetMarkerStyle(kFullDotMedium);
134
135	TGraph gt0, gt1;
136	gt0.SetMarkerColor(kRed);
137	gt1.SetMarkerColor(kBlue);
138	gt0.SetMarkerStyle(kFullDotLarge);
139	gt1.SetMarkerStyle(kFullDotLarge);
140
141	TGraph gz0, gz1;
142	gz0.SetMarkerColor(kRed);
143	gz1.SetMarkerColor(kBlue);
144	gz0.SetMarkerStyle(kFullDotLarge);
145	gz1.SetMarkerStyle(kFullDotLarge);
146
147	Int_t first = -1;
148	Int_t last = -1;
149
150	while ((row=res.Next()))
151	{
152	const char date = (row)[0];
153	const char zd = (row)[1];
154	const char val = (row)[2];
155	const char snum = res.GetFieldCount()>3 ? (row)[3] : 0;
156	const char verr = res.GetFieldCount()>4 ? (row)[5] : 0;
157	if (!date \|\| !val \|\| !zd)
158	continue;
159
160	MTime t(date);
161	if (!t.SetSqlDateTime(date))
162	continue;
163
164	if (fRequestPeriod>0 && MAstro::GetMagicPeriod(t.GetMjd())!=fRequestPeriod)
165	continue;
166
167	if (first<0)
168	first = TMath::Nint(TMath::Floor(t.GetMjd()));
169	last = TMath::Nint(TMath::Ceil(t.GetMjd()));
170
171	UInt_t seq = snum ? atoi(snum) : 0;
172
173	Float_t value = atof(val);
174	Float_t zenith = atof(zd);
175
176	if (fDataSet)
177	{
178	if (fDataSet->HasOnSequence(seq))
179	{
180	gt1.SetPoint(gt1.GetN(), t.GetAxisTime(), value);
181	gz1.SetPoint(gz1.GetN(), zenith, value);
182	}
183
184	if (fDataSet->HasOffSequence(seq))
185	{
186	gt0.SetPoint(gt0.GetN(), t.GetAxisTime(), value);
187	gz0.SetPoint(gz0.GetN(), zenith, value);
188	}
189	}
190
191	gt.SetPoint(gt.GetN(), t.GetAxisTime(), value);
192	gz.SetPoint(gz.GetN(), zenith, value);
193
194	if (verr)
195	static_cast<TGraphErrors&>(gt).SetPointError(gt.GetN()-1, 0, atof(verr));
196	}
197
198	// If this is done earlier the plots remain empty since root 5.12/00
199	if (fmax>fmin)
200	{
201	gt.SetMinimum(fmin);
202	gt.SetMaximum(fmax);
203	gz.SetMinimum(fmin);
204	gz.SetMaximum(fmax);
205	}
206
207	gROOT->SetSelectedPad(0);
208
209	TString title = fNameTab.IsNull() ? name(name.First('.')+2, name.Length()) : fNameTab;
210	TCanvas &c = fDisplay ? fDisplay->AddTab(title) : *new TCanvas;
211	c.SetFillColor(kWhite);
212	c.SetBorderMode(0);
213	c.Divide(1,2);
214
215	cerr << setprecision(4) << setw(10) << title << ": ";
216	if (gt.GetN()==0)
217	{
218	cerr << " <empty>" << endl;
219	return;
220	}
221	cerr << setw(8) << gt.GetMean(2) << "+-" << setw(8) << gt.GetRMS(2) << " ";
222	if (gt0.GetN()>0 \|\| gt1.GetN()>0)
223	{
224	cerr << setw(8) << gt1.GetMean(2) << "+-" << setw(8) << gt1.GetRMS(2) << " ";
225	cerr << setw(8) << gt0.GetMean(2) << "+-" << setw(8) << gt0.GetRMS(2);
226	}
227	cerr << endl;
228
229	TVirtualPad *pad = gPad;
230	pad->cd(2);
231	gPad->SetBorderMode(0);
232	gPad->SetFrameBorderMode(0);
233	gPad->SetGridy();
234
235	gPad->SetLeftMargin(0.06);
236	gPad->SetRightMargin(0.06);
237	gPad->SetBottomMargin(0.08);
238
239	TH1 *h = gt.GetHistogram();
240
241	h->SetXTitle("Time");
242	h->SetYTitle(name);
243	h->GetXaxis()->SetTimeDisplay(1);
244	h->GetYaxis()->SetTitleOffset(0.8);
245	h->GetXaxis()->SetTitleOffset(1.0);
246	h->GetXaxis()->SetLabelOffset(0.01);
247
248	gt.DrawClone("AP");
249	if (gt0.GetN()>0)
250	gt0.DrawClone("P");
251	if (gt1.GetN()>0)
252	gt1.DrawClone("P");
253
254	TLine l;
255	TText t;
256	Int_t num=0;
257	l.SetLineStyle(kDotted);
258	l.SetLineColor(kBlue);
259	t.SetTextColor(kBlue);
260	l.SetLineWidth(1);
261	t.SetTextSize(h->GetXaxis()->GetLabelSize());
262	t.SetTextAlign(21);
263	Int_t p0 = MAstro::GetMagicPeriod(first);
264	for (Int_t p = first; p<last; p++)
265	{
266	Int_t p1 = MAstro::GetMagicPeriod(p);
267	if (p1!=p0)
268	{
269	l.DrawLine(MTime(p).GetAxisTime(), h->GetMinimum(), MTime(p).GetAxisTime(), h->GetMaximum());
270	t.DrawText(MTime(p+15).GetAxisTime(), h->GetMaximum(), Form("%d", p1));
271	num++;
272	}
273	p0 = p1;
274	}
275	if (num<4)
276	gPad->SetGridx();
277
278	const Double_t min = fHistMin>fHistMax ? h->GetMinimum()-resolution/2 : fHistMin;
279	const Double_t max = fHistMin>fHistMax ? h->GetMaximum()+resolution/2 : fHistMax;
280
281	// Use this to save the pad with the time development to a file
282	//gPad->SaveAs(Form("plotdb-%s.eps", title.Data()));
283
284	pad->cd(1);
285	gPad->SetBorderMode(0);
286	gPad->SetFrameBorderMode(0);
287	gPad->Divide(2,1);
288
289	TVirtualPad *pad2 = gPad;
290	pad2->cd(1);
291	gPad->SetBorderMode(0);
292	gPad->SetFrameBorderMode(0);
293	gPad->SetGridx();
294	gPad->SetGridy();
295
296	const Int_t n = resolution>0 ? TMath::Nint((max-min)/resolution) : 50;
297
298	TH1F hist("Hist", Form("Distribution of %s", fDescription.IsNull() ? name.Data() : fDescription.Data()), n, min, max);
299	hist.SetDirectory(0);
300
301	for (int i=0; i<gt.GetN(); i++)
302	hist.Fill(gt.GetY()[i]);
303
304	if (fDescription.IsNull())
305	hist.SetXTitle(name);
306	hist.SetYTitle("Counts");
307
308	hist.DrawCopy("");
309
310	pad2->cd(2);
311	gPad->SetBorderMode(0);
312	gPad->SetFrameBorderMode(0);
313	gPad->SetGridy();
314
315	TH1 *h2 = gz.GetHistogram();
316
317	h2->SetXTitle("Zd");
318	h2->SetYTitle(name);
319
320	gz.DrawClone("AP");
321
322	if (gz0.GetN()>0)
323	gz0.DrawClone("P");
324	if (gz1.GetN()>0)
325	gz1.DrawClone("P");
326	}
327
328	public:
329	MPlot(MSQLMagic &server) : fServer(server), fDataSet(NULL),
330	fRequestPeriod(-1), fPlotMin(0), fPlotMax(-1), fHistMin(0), fHistMax(-1), fGroupBy(kNone)
331	{
332	}
333	~MPlot()
334	{
335	if (fDataSet)
336	delete fDataSet;
337	}
338	void SetDataSet(const TString filename)
339	{
340	if (fDataSet)
341	{
342	delete fDataSet;
343	fDataSet = NULL;
344	}
345	if (!filename.IsNull())
346	fDataSet = new MDataSet(filename);
347	}
348	void SetPlotRange(Float_t min, Float_t max, Int_t n=5) { fPlotMin = min; fPlotMax = max; }
349	void SetHistRange(Float_t min, Float_t max) { fHistMin = min; fHistMax = max; }
350	void SetRequestRange(const char from="", const char to="") { fRequestFrom = from; fRequestTo = to; }
351	void SetRequestPeriod(Int_t n=-1) { fRequestPeriod = n; }
352	void SetCondition(const char *cond="") { fCondition = cond; }
353	void SetDescription(const char d, const char t=0) { fDescription = d; fNameTab = t; }
354	void SetGroupBy(GroupBy_t b=kGroupByWeek) { fGroupBy=b; }
355	void SetPrimaryDate(const char *ts) { fPrimaryDate=ts; }
356	void SetPrimaryNumber(const char *ts) { fPrimaryNumber=ts; }
357	void SetSecondary(const char *ts) { fSecondary=ts; }
358
359	Int_t QueryKeyOfSource(TString src)
360	{
361	return fServer.QueryKeyOfName("Object", src, kFALSE);
362	}
363
364	Bool_t Plot(const char *value, Float_t min=0, Float_t max=-1, Float_t resolution=0)
365	{
366	TString named = fPrimaryDate;
367	TString named2 = fSecondary;
368	TString namev = value;
369
370	TString tablev = namev(0, namev.First('.'));
371	TString valuev = namev(namev.First('.')+1, namev.Length());
372
373	TString tabled = named(0, named.First('.'));
374	TString valued = named(named.First('.')+1, named.Length());
375
376	TString query="SELECT ";
377	query += valued;
378	if (fGroupBy==kNone)
379	{
380	query += ", ";
381	query += fSecondary;
382	query += ", ";
383	query += value;
384	query += ", ";
385	query += fPrimaryNumber;
386	query += " ";
387	}
388	else
389	{
390	query += ", AVG(";
391	query += fSecondary;
392	query += "), AVG(";
393	query += value;
394	query += "), ";
395	query += fPrimaryNumber;
396	query += ", STD(";
397	query += fSecondary;
398	query += "), STD(";
399	query += value;
400	query += ") ";
401	}
402
403	switch (fGroupBy)
404	{
405	case kNone:
406	case kGroupByPrimary:
407	break;
408	case kGroupByHour:
409	query += Form(", DATE_FORMAT(%s, '%%Y-%%m-%%d %%H') AS %s ", named.Data(), valued.Data());
410	break;
411	case kGroupByNight:
412	query += Form(", DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y-%%m-%%d') AS %s ", named.Data(), valued.Data());
413	break;
414	case kGroupByWeek:
415	query += Form(", DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%x%%v') AS %s ", named.Data(), valued.Data());
416	break;
417	case kGroupByMonth:
418	query += Form(", DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y-%%m') AS %s ", named.Data(), valued.Data());
419	break;
420	case kGroupByYear:
421	query += Form(", DATE_FORMAT(ADDDATE(%s,Interval 12 hour), '%%Y') AS %s ", named.Data(), valued.Data());
422	break;
423	}
424
425	query += Form("FROM %s ", tabled.Data());
426
427	TString where(fCondition);
428
429	if (!where.IsNull())
430	where += " AND ";
431
432	where += fCondition;
433	where += " ";
434
435	// ------------------------------
436
437	query += fServer.GetJoins(tabled, query+" "+where);
438
439	if (!where.IsNull())
440	{
441	query += "WHERE ";
442	query += where;
443	}
444
445	if (fGroupBy!=kNone)
446	{
447	query += Form("GROUP BY %s ", valued.Data());
448	//query += Form(" HAVING COUNT(%s)=(COUNT(*)+1)/2 ", valuev.Data());
449	}
450	query += Form("ORDER BY %s ", fPrimaryDate.Data());
451
452
453	// ------------------------------
454
455	TSQLResult *res = fServer.Query(query);
456	if (!res)
457	return kFALSE;
458
459	if (max>min)
460	PlotTable(*res, namev, min, max, resolution);
461	else
462	PlotTable(*res, namev, fPlotMin, fPlotMax, resolution);
463
464	delete res;
465	return kTRUE;
466	}
467	};
468
469	void plotall(MPlot &plot, TString source)
470	{
471	plot.SetPrimaryDate("OpticalData.fTimestamp");
472	plot.SetPrimaryNumber("OpticalData.fTimestamp");
473	plot.SetSecondary("OpticalData.fZenithDistance");
474
475	TString cond = "fStatusKEY=13";
476	if (!source.IsNull())
477	{
478	const Int_t key = plot.QueryKeyOfSource(source);
479	if (key<0)
480	return;
481	cond += Form(" AND Object.fObjectKEY=%d", key);
482
483	}
484	plot.SetCondition(cond);
485
486	plot.SetDescription("Exposure;T_{E} [s]", "Expo");
487	plot.Plot("OpticalData.fExposure", 0, 900, 60);
488
489	plot.SetDescription("Sky Level;B [s^{-1}]", "SkyLvl");
490	plot.Plot("OpticalData.fSkyLevel/OpticalData.fExposure", 0, 5.0, 0.01);
491
492	plot.SetDescription("Full Width Half Maximum;FWHM [s^{-1}]", "Fwhm");
493	plot.Plot("OpticalData.fFWHM/OpticalData.fExposure", 0, 0.05, 0.001);
494
495	plot.SetDescription("Aperture Radius;R_{A}", "ApRad");
496	plot.Plot("OpticalData.fApertureRadius", 0, 10, 1);
497
498	/*
499	plot.SetDescription("Instrumental Magnitude;M_{I}\\cdot s^{-1}", "InstMag/s");
500	plot.Plot("OpticalData.fInstrumentalMag/OpticalData.fExposure", 0, 0.2, 0.005);
501
502	plot.SetDescription("Instrumental Magnitude Error;\\sigma_{M}\\cdot s^{-1}", "MagErr/s");
503	plot.Plot("OpticalData.fInstrumentalMagErr/OpticalData.fExposure", 0, 0.01, 0.0002);
504
505	plot.SetDescription("Instrumental Magnitude;M_{I}", "InstMag");
506	plot.Plot("OpticalData.fInstrumentalMag", 0, 30, 0.5);
507
508	plot.SetDescription("Instrumental Magnitude Error;\\sigma_{M}", "MagErr");
509	plot.Plot("OpticalData.fInstrumentalMagErr", 0, 1, 0.01);
510	*/
511
512	plot.SetDescription("m_{1};m_{1}", "M1");
513	plot.Plot("OpticalData.fInstrumentalMag+2.5*log10(OpticalData.fExposure)", 10, 35, 0.2);
514
515	cond += " AND Object.fObjectName NOT LIKE '%/BL' AND NOT ISNULL(Object.fMagnitude) ";
516	plot.SetCondition(cond);
517
518	TString ext("3080/25.0pow(10, (OpticalData.fInstrumentalMag+2.5log10(OpticalData.fExposure)-Object.fMagnitude)/-2.5)");
519	ext += "+0.0028*fZenithDistance-0.08";
520
521	plot.SetGroupBy(MPlot::kGroupByPrimary);
522	plot.SetDescription("m_{1}-m_{true} (Extinction per Image);m_{1}-m_{true}", "ExtImg");
523	plot.Plot(ext/Form("AVG(%s)", ext.Data())/, 0.05, 1.2, 0.01);
524
525	plot.SetGroupBy(MPlot::kGroupByHour);
526	plot.SetDescription("m_{1}-m_{true} (Extinction per Hour);m_{1}-m_{true}", "ExtHour");
527	plot.Plot(ext/Form("AVG(%s), date_format(fTimeStamp, '%%Y-%%m-%%d %%H') as fTimeStamp", ext.Data())/,
528	0.5, 1.2, 0.01);
529
530	plot.SetGroupBy(MPlot::kGroupByNight);
531	plot.SetDescription("m_{1}-m_{true} (Extinction per Night);m_{1}-m_{true}", "ExtNight");
532	plot.Plot(ext/Form("AVG(%s), date_format(adddate(fTimeStamp,Interval 12 hour),'%%Y-%%m-%%d') as fTimeStamp", ext.Data())/,
533	0.5, 1.2, 0.01);
534	}
535
536	int plotoptical(TString from, TString to, const char *source=0)
537	{
538	TEnv env("sql.rc");
539
540	MSQLMagic serv(env);
541	if (!serv.IsConnected())
542	{
543	cout << "ERROR - Connection to database failed." << endl;
544	return 0;
545	}
546
547	cout << "plotoptical" << endl;
548	cout << "-----------" << endl;
549	cout << endl;
550	cout << "Connected to " << serv.GetName() << endl;
551	cout << endl;
552
553	MStatusDisplay *d = new MStatusDisplay;
554	d->SetWindowName(serv.GetName());
555	d->SetTitle(serv.GetName());
556
557	MPlot plot(serv);
558	// plot.SetDataSet(dataset);
559	plot.SetDisplay(d);
560	plot.SetRequestRange(from, to);
561	plotall(plot, source);
562	d->SaveAsRoot("plotoptical.root");
563	d->SaveAsPS("plotoptical.ps");
564
565	return 1;
566	}
567
568	int plotoptical(const char *source)
569	{
570	TEnv env("sql.rc");
571
572	MSQLMagic serv(env);
573	if (!serv.IsConnected())
574	{
575	cout << "ERROR - Connection to database failed." << endl;
576	return 0;
577	}
578
579	cout << "plotoptical" << endl;
580	cout << "-----------" << endl;
581	cout << endl;
582	cout << "Connected to " << serv.GetName() << endl;
583	cout << endl;
584
585	MStatusDisplay *d = new MStatusDisplay;
586	d->SetWindowName(serv.GetName());
587	d->SetTitle(serv.GetName());
588
589	MPlot plot(serv);
590	// plot.SetDataSet(ds);
591	plot.SetDisplay(d);
592	plot.SetRequestRange("", "");
593	plotall(plot, source);
594	d->SaveAsRoot("plotoptical.root");
595	d->SaveAsPS("plotoptical.ps");
596
597	return 1;
598	}
599
600	int plotoptical(Int_t period, const char *source="")
601	{
602	TEnv env("sql.rc");
603
604	MSQLMagic serv(env);
605	if (!serv.IsConnected())
606	{
607	cout << "ERROR - Connection to database failed." << endl;
608	return 0;
609	}
610
611	cout << "plotoptical" << endl;
612	cout << "-----------" << endl;
613	cout << endl;
614	cout << "Connected to " << serv.GetName() << endl;
615	cout << endl;
616
617	MStatusDisplay *d = new MStatusDisplay;
618	d->SetWindowName(serv.GetName());
619	d->SetTitle(serv.GetName());
620
621	MPlot plot(serv);
622	// plot.SetDataSet(dataset);
623	plot.SetDisplay(d);
624	plot.SetRequestPeriod(period);
625	plotall(plot, source);
626	d->SaveAsRoot("plotoptical.root");
627	d->SaveAsPS("plotoptical.ps");
628
629	return 1;
630	}
631
632	int plotoptical()
633	{
634	return plotoptical("", "");
635	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: