Context Navigation

source: trunk/FACT++/src/makeplots.cc@ 17651

Visit:

Last change on this file since 17651 was 17651, checked in by tbretz, 11 years ago
Updated with the current prediction from the calibration paper; make plots only for VHE sources.
File size: 11.5 KB

Line
1	#include "externals/nova.h"
2
3	#include "Database.h"
4
5	#include "Time.h"
6	#include "Configuration.h"
7
8	#include <TROOT.h>
9	#include <TH1.h>
10	#include <TGraph.h>
11	#include <TCanvas.h>
12	#include <TLegend.h>
13
14	using namespace std;
15	using namespace Nova;
16
17	// ------------------------------------------------------------------------
18
19	void CheckForGap(TCanvas &c, TGraph &g, double axis)
20	{
21	if (g.GetN()==0 \|\| axis-g.GetX()[g.GetN()-1]<450)
22	return;
23
24	c.cd();
25	((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
26	while (g.GetN())
27	g.RemovePoint(0);
28	}
29
30	void DrawClone(TCanvas &c, TGraph &g)
31	{
32	if (g.GetN()==0)
33	return;
34
35	c.cd();
36	((TGraph*)g.DrawClone("C"))->SetBit(kCanDelete);
37	}
38
39	// ========================================================================
40	// ========================================================================
41	// ========================================================================
42
43	void SetupConfiguration(Configuration &conf)
44	{
45	po::options_description control("Smart FACT");
46	control.add_options()
47	("ra", var<double>(), "Source right ascension")
48	("dec", var<double>(), "Source declination")
49	("date-time", var<string>(), "SQL time (UTC)")
50	("source-database", var<string>(""), "Database link as in\n\tuser:password@server[:port]/database.")
51	("max-current", var<double>(75), "Maximum current to display in other plots.")
52	("max-zd", var<double>(75), "Maximum zenith distance to display in other plots")
53	("no-limits", po_switch(), "Switch off limits in plots")
54	;
55
56	po::positional_options_description p;
57	p.add("date-time", 1); // The first positional options
58
59	conf.AddOptions(control);
60	conf.SetArgumentPositions(p);
61	}
62
63	void PrintUsage()
64	{
65	cout <<
66	"makeplots - The astronomy plotter\n"
67	"\n"
68	"Calculates several plots for the sources in the database\n"
69	"helpful or needed for scheduling. The Plot is always calculated\n"
70	"for the night which starts at the same so. So no matter if\n"
71	"you specify '1974-09-09 00:00:00' or '1974-09-09 21:56:00'\n"
72	"the plots will refer to the night 1974-09-09/1974-09-10.\n"
73	"The advantage is that specification of the date as in\n"
74	"1974-09-09 is enough. Time axis starts and ends at nautical\n"
75	"twilight which is 12deg below horizon.\n"
76	"\n"
77	"Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
78	cout << endl;
79	}
80
81	int main(int argc, const char* argv[])
82	{
83	gROOT->SetBatch();
84
85	Configuration conf(argv[0]);
86	conf.SetPrintUsage(PrintUsage);
87	SetupConfiguration(conf);
88
89	if (!conf.DoParse(argc, argv))
90	return 127;
91
92	if (conf.Has("ra")^conf.Has("dec"))
93	{
94	cout << "ERROR - Either --ra or --dec missing." << endl;
95	return 1;
96	}
97
98	// ------------------ Eval config ---------------------
99
100	Time time;
101	if (conf.Has("date-time"))
102	time.SetFromStr(conf.Get<string>("date-time"));
103
104	const double max_current = conf.Get<double>("max-current");
105	const double max_zd = conf.Get<double>("max-zd");
106	const double no_limits = conf.Get<bool>("no-limits");
107
108	// -12: nautical
109	// Sun set with the same date than th provided date
110	// Sun rise on the following day
111	const RstTime sun_set = GetSolarRst(time.JD()-0.5, -12);
112	const RstTime sun_rise = GetSolarRst(time.JD()+0.5, -12);
113
114	const double jd = floor(time.Mjd())+2400001;
115	const double mjd = floor(time.Mjd())+49718+0.5; // root time
116
117	cout << "Time: " << time << endl;
118	cout << "Base: " << Time(jd-0.5) << endl;
119	cout << "Set: " << Time(sun_set.set) << endl;
120	cout << "Rise: " << Time(sun_rise.rise) << endl;
121
122	const double jd0 = fmod(sun_set.set, 1); // ~0.3
123	const double jd1 = fmod(sun_rise.rise, 1); // ~0.8
124
125	const string fDatabase = conf.Get<string>("source-database");
126
127	// ------------------ Precalc moon coord ---------------------
128
129	vector<tuple<EquPosn, double, double>> fMoonCoords;
130	vector<EquPosn> fSunCoords;
131
132	for (double h=0; h<1; h+=1./(24*12))
133	{
134	const EquPosn sun = GetSolarEquCoords(jd+h);
135	const EquPosn moon = GetLunarEquCoords(jd+h, 0.01);
136	const double disk = GetLunarDisk(jd+h);
137	const double edist = GetLunarEarthDist(jd+h)/384400;
138
139	fMoonCoords.emplace_back(moon, disk, edist);
140	fSunCoords.emplace_back(sun);
141	}
142
143	// ------------- Get Sources from databasse ---------------------
144
145	const mysqlpp::StoreQueryResult res =
146	Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1").store();
147
148	// ------------- Create canvases and frames ---------------------
149
150	// It is important to use an offset which is larger than
151	// 1970-01-01 00:00:00. This one will not work if your
152	// local time zone is positive!
153	TH1S hframe("", "", 1, (mjd+jd0)243600, (mjd+jd1)243600);
154	hframe.SetStats(kFALSE);
155	hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
156	hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
157	hframe.GetXaxis()->CenterTitle();
158	hframe.GetYaxis()->CenterTitle();
159	hframe.GetXaxis()->SetTimeDisplay(true);
160	hframe.GetYaxis()->SetTitleSize(0.040);
161	hframe.GetXaxis()->SetTitleSize(0.040);
162	hframe.GetXaxis()->SetTitleOffset(1.1);
163	hframe.GetYaxis()->SetLabelSize(0.040);
164	hframe.GetXaxis()->SetLabelSize(0.040);
165
166	TCanvas c1;
167	c1.SetFillColor(kWhite);
168	c1.SetBorderMode(0);
169	c1.SetFrameBorderMode(0);
170	c1.SetLeftMargin(0.085);
171	c1.SetRightMargin(0.01);
172	c1.SetTopMargin(0.03);
173	c1.SetGrid();
174	hframe.GetYaxis()->SetTitle("Altitude [deg]");
175	hframe.SetMinimum(15);
176	hframe.SetMaximum(90);
177	hframe.DrawCopy();
178
179	TCanvas c2;
180	c2.SetFillColor(kWhite);
181	c2.SetBorderMode(0);
182	c2.SetFrameBorderMode(0);
183	c2.SetLeftMargin(0.085);
184	c2.SetRightMargin(0.01);
185	c2.SetTopMargin(0.03);
186	c2.SetGrid();
187	hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
188	hframe.SetMinimum(0);
189	hframe.SetMaximum(100);
190	hframe.DrawCopy();
191
192	TCanvas c3;
193	c3.SetFillColor(kWhite);
194	c3.SetBorderMode(0);
195	c3.SetFrameBorderMode(0);
196	c3.SetLeftMargin(0.085);
197	c3.SetRightMargin(0.01);
198	c3.SetTopMargin(0.03);
199	c3.SetGrid();
200	c3.SetLogy();
201	hframe.GetYaxis()->SetTitle("Estimated relative threshold");
202	hframe.SetMinimum(0.9);
203	hframe.SetMaximum(180);
204	hframe.DrawCopy();
205
206	TCanvas c4;
207	c4.SetFillColor(kWhite);
208	c4.SetBorderMode(0);
209	c4.SetFrameBorderMode(0);
210	c4.SetLeftMargin(0.085);
211	c4.SetRightMargin(0.01);
212	c4.SetTopMargin(0.03);
213	c4.SetGrid();
214	hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
215	hframe.SetMinimum(0);
216	hframe.SetMaximum(180);
217	hframe.DrawCopy();
218
219	Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
220	Int_t style[] = { kSolid, kDashed, kDotted };
221
222	TLegend leg(0, 0, 1, 1);
223
224	// ------------- Loop over sources ---------------------
225
226	Int_t cnt=0;
227	for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
228	{
229	// Eval row
230	const string name = (*v)[0].c_str();
231
232	EquPosn pos;
233	pos.ra = double((v)[1])15;
234	pos.dec = double((*v)[2]);
235
236	// Create graphs
237	TGraph g1, g2, g3, g4, gr, gm;
238	g1.SetName(name.data());
239	g2.SetName(name.data());
240	g3.SetName(name.data());
241	g4.SetName(name.data());
242	g1.SetLineWidth(2);
243	g2.SetLineWidth(2);
244	g3.SetLineWidth(2);
245	g4.SetLineWidth(2);
246	gm.SetLineWidth(1);
247	g1.SetLineStyle(style[cnt/6]);
248	g2.SetLineStyle(style[cnt/6]);
249	g3.SetLineStyle(style[cnt/6]);
250	g4.SetLineStyle(style[cnt/6]);
251	g1.SetLineColor(color[cnt%6]);
252	g2.SetLineColor(color[cnt%6]);
253	g3.SetLineColor(color[cnt%6]);
254	g4.SetLineColor(color[cnt%6]);
255	gm.SetLineColor(kYellow);
256
257	if (cnt==0)
258	leg.AddEntry(gm.Clone(), "Moon", "l");
259	leg.AddEntry(g1.Clone(), name.data(), "l");
260
261	// Loop over 24 hours
262	int i=0;
263	for (double h=0; h<1; h+=1./(24*12), i++)
264	{
265	// check if it is betwene sun-rise and sun-set
266	if (h<jd0 \|\| h>jd1)
267	continue;
268
269	// get local position of source
270	const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
271
272	// Get sun and moon properties and
273	const EquPosn sun = fSunCoords[i];
274	const EquPosn moon = get<0>(fMoonCoords[i]);
275	const double disk = get<1>(fMoonCoords[i]);
276	const double edist = get<2>(fMoonCoords[i]);
277	const HrzPosn hrzm = GetHrzFromEqu(moon, jd+h);
278	const ZdAzPosn hrzs = GetHrzFromEqu(sun, jd+h);
279
280	if (v==res.begin())
281	cout << Time(jd+h) <<" " << 90-hrzm.alt << endl;
282
283	// Distance between source and moon
284	const double angle = GetAngularSeparation(moon, pos);
285
286	// Current prediction
287	const double sin_malt = hrzm.alt<0 ? 0 : sin(hrzm.alt*M_PI/180);
288	const double cos_mdist = cos(angle*M_PI/180);
289	const double sin_szd = sin(hrzs.zd*M_PI/180);
290
291	const double k0 = pow(disk, 2.63);
292	const double k1 = pow(sin_malt, 0.60);
293	const double k2 = pow(edist, -2.00);
294	const double k3 = exp(0.67cos_mdistcos_mdistcos_mdistcos_mdist);
295	const double k4 = exp(-97.8+105.8sin_szdsin_szd);
296
297	const double cur = 6.2 + 95.7k0k1k2k3 + k4;
298
299	// Relative energy threshold prediction
300	//const double cs = cos((90+hrz.alt)*M_PI/180);
301	//const double ratio = (10.sqrt(409600.cscs+9009.) + 6400.cs - 60.)/10.;
302	const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
303
304	// Add points to curve
305	const double axis = (mjd+h)243600;
306
307	// If there is a gap of more than one bin, start a new curve
308	CheckForGap(c1, g1, axis);
309	CheckForGap(c1, gm, axis);
310	CheckForGap(c2, g2, axis);
311	CheckForGap(c3, g3, axis);
312	CheckForGap(c4, g4, axis);
313
314	// Add data
315	if (no_limits \|\| cur<max_current)
316	g1.SetPoint(g1.GetN(), axis, hrz.alt);
317
318	if (no_limits \|\| 90-hrz.alt<max_zd)
319	g2.SetPoint(g2.GetN(), axis, cur);
320
321	if (no_limits \|\| (cur<max_current && 90-hrz.alt<max_zd))
322	g3.SetPoint(g3.GetN(), axis, ratio*cur/8.1);
323
324	if (no_limits \|\| (cur<max_current && 90-hrz.alt<max_zd))
325	g4.SetPoint(g4.GetN(), axis, angle);
326
327	if (cnt==0)
328	gm.SetPoint(gm.GetN(), axis, hrzm.alt);
329	}
330
331	if (cnt==0)
332	DrawClone(c1, gm);
333
334	DrawClone(c1, g1);
335	DrawClone(c2, g2);
336	DrawClone(c3, g3);
337	DrawClone(c4, g4);
338	}
339
340
341	// Save three plots
342	TCanvas c5;
343	c5.SetFillColor(kWhite);
344	c5.SetBorderMode(0);
345	c5.SetFrameBorderMode(0);
346	leg.Draw();
347
348	const string t = Time(jd).GetAsStr("%Y%m%d");
349
350	c1.SaveAs((t+"-ZenithDistance.eps").c_str());
351	c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
352	c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
353	c4.SaveAs((t+"-MoonDist.eps").c_str());
354	c5.SaveAs((t+"-Legend.eps").c_str());
355
356	c1.SaveAs((t+"-ZenithDistance.root").c_str());
357	c2.SaveAs((t+"-PredictedCurrent.root").c_str());
358	c3.SaveAs((t+"-RelativeThreshold.root").c_str());
359	c4.SaveAs((t+"-MoonDist.root").c_str());
360
361	return 0;
362	}

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

Download in other formats: