Context Navigation

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

Visit:

Last change on this file since 17667 was 17656, checked in by tbretz, 11 years ago
To calculate the relative energy threshold not the current must be used but the trigger threshold -- used the numbers from the paper; adapted y-scale of the plot accordingly
File size: 11.8 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>(100), "Maximum current to display in other plots.")
52	("max-zd", var<double>(50), "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\n";
78	// "Usage: makeplots sql-datetime [--ra={ra} --dec={dec}]\n";
79	cout << endl;
80	}
81
82	int main(int argc, const char* argv[])
83	{
84	gROOT->SetBatch();
85
86	Configuration conf(argv[0]);
87	conf.SetPrintUsage(PrintUsage);
88	SetupConfiguration(conf);
89
90	if (!conf.DoParse(argc, argv))
91	return 127;
92
93	// if (conf.Has("ra")^conf.Has("dec"))
94	// {
95	// cout << "ERROR - Either --ra or --dec missing." << endl;
96	// return 1;
97	// }
98
99	// ------------------ Eval config ---------------------
100
101	Time time;
102	if (conf.Has("date-time"))
103	time.SetFromStr(conf.Get<string>("date-time"));
104
105	const double max_current = conf.Get<double>("max-current");
106	const double max_zd = conf.Get<double>("max-zd");
107	const double no_limits = conf.Get<bool>("no-limits");
108
109	// -12: nautical
110	// Sun set with the same date than th provided date
111	// Sun rise on the following day
112	const RstTime sun_set = GetSolarRst(time.JD()-0.5, -12);
113	const RstTime sun_rise = GetSolarRst(time.JD()+0.5, -12);
114
115	const double jd = floor(time.Mjd())+2400001;
116	const double mjd = floor(time.Mjd())+49718+0.5; // root time
117
118	cout << "Time: " << time << endl;
119	cout << "Base: " << Time(jd-0.5) << endl;
120	cout << "Set: " << Time(sun_set.set) << endl;
121	cout << "Rise: " << Time(sun_rise.rise) << endl;
122
123	const double jd0 = fmod(sun_set.set, 1); // ~0.3
124	const double jd1 = fmod(sun_rise.rise, 1); // ~0.8
125
126	const string fDatabase = conf.Get<string>("source-database");
127
128	// ------------------ Precalc moon coord ---------------------
129
130	vector<tuple<EquPosn, double, double>> fMoonCoords;
131	vector<EquPosn> fSunCoords;
132
133	for (double h=0; h<1; h+=1./(24*12))
134	{
135	const EquPosn sun = GetSolarEquCoords(jd+h);
136	const EquPosn moon = GetLunarEquCoords(jd+h, 0.01);
137	const double disk = GetLunarDisk(jd+h);
138	const double edist = GetLunarEarthDist(jd+h)/384400;
139
140	fMoonCoords.emplace_back(moon, disk, edist);
141	fSunCoords.emplace_back(sun);
142	}
143
144	// ------------- Get Sources from databasse ---------------------
145
146	const mysqlpp::StoreQueryResult res =
147	Database(fDatabase).query("SELECT fSourceName, fRightAscension, fDeclination FROM Source WHERE fSourceTypeKEY=1").store();
148
149	// ------------- Create canvases and frames ---------------------
150
151	// It is important to use an offset which is larger than
152	// 1970-01-01 00:00:00. This one will not work if your
153	// local time zone is positive!
154	TH1S hframe("", "", 1, (mjd+jd0)243600, (mjd+jd1)243600);
155	hframe.SetStats(kFALSE);
156	hframe.GetXaxis()->SetTimeFormat("%Hh%M%F1995-01-01 00:00:00 GMT");
157	hframe.GetXaxis()->SetTitle((Time(jd).GetAsStr("%d/%m/%Y")+" - "+Time(jd+1).GetAsStr("%d/%m/%Y")+" [UTC]").c_str());
158	hframe.GetXaxis()->CenterTitle();
159	hframe.GetYaxis()->CenterTitle();
160	hframe.GetXaxis()->SetTimeDisplay(true);
161	hframe.GetYaxis()->SetTitleSize(0.040);
162	hframe.GetXaxis()->SetTitleSize(0.040);
163	hframe.GetXaxis()->SetTitleOffset(1.1);
164	hframe.GetYaxis()->SetLabelSize(0.040);
165	hframe.GetXaxis()->SetLabelSize(0.040);
166
167	TCanvas c1;
168	c1.SetFillColor(kWhite);
169	c1.SetBorderMode(0);
170	c1.SetFrameBorderMode(0);
171	c1.SetLeftMargin(0.085);
172	c1.SetRightMargin(0.01);
173	c1.SetTopMargin(0.03);
174	c1.SetGrid();
175	hframe.GetYaxis()->SetTitle("Altitude [deg]");
176	hframe.SetMinimum(15);
177	hframe.SetMaximum(90);
178	hframe.DrawCopy();
179
180	TCanvas c2;
181	c2.SetFillColor(kWhite);
182	c2.SetBorderMode(0);
183	c2.SetFrameBorderMode(0);
184	c2.SetLeftMargin(0.085);
185	c2.SetRightMargin(0.01);
186	c2.SetTopMargin(0.03);
187	c2.SetGrid();
188	hframe.GetYaxis()->SetTitle("Predicted Current [\\muA]");
189	hframe.SetMinimum(0);
190	hframe.SetMaximum(100);
191	hframe.DrawCopy();
192
193	TCanvas c3;
194	c3.SetFillColor(kWhite);
195	c3.SetBorderMode(0);
196	c3.SetFrameBorderMode(0);
197	c3.SetLeftMargin(0.085);
198	c3.SetRightMargin(0.01);
199	c3.SetTopMargin(0.03);
200	c3.SetGrid();
201	c3.SetLogy();
202	hframe.GetYaxis()->SetTitle("Estimated relative threshold");
203	hframe.GetYaxis()->SetMoreLogLabels();
204	hframe.SetMinimum(0.9);
205	hframe.SetMaximum(11);
206	hframe.DrawCopy();
207
208	TCanvas c4;
209	c4.SetFillColor(kWhite);
210	c4.SetBorderMode(0);
211	c4.SetFrameBorderMode(0);
212	c4.SetLeftMargin(0.085);
213	c4.SetRightMargin(0.01);
214	c4.SetTopMargin(0.03);
215	c4.SetGrid();
216	hframe.GetYaxis()->SetTitle("Distance to moon [deg]");
217	hframe.SetMinimum(0);
218	hframe.SetMaximum(180);
219	hframe.DrawCopy();
220
221	Int_t color[] = { kBlack, kRed, kBlue, kGreen, kCyan, kMagenta };
222	Int_t style[] = { kSolid, kDashed, kDotted };
223
224	TLegend leg(0, 0, 1, 1);
225
226	// ------------- Loop over sources ---------------------
227
228	Int_t cnt=0;
229	for (vector<mysqlpp::Row>::const_iterator v=res.begin(); v<res.end(); v++, cnt++)
230	{
231	// Eval row
232	const string name = (*v)[0].c_str();
233
234	EquPosn pos;
235	pos.ra = double((v)[1])15;
236	pos.dec = double((*v)[2]);
237
238	// Create graphs
239	TGraph g1, g2, g3, g4, gr, gm;
240	g1.SetName(name.data());
241	g2.SetName(name.data());
242	g3.SetName(name.data());
243	g4.SetName(name.data());
244	g1.SetLineWidth(2);
245	g2.SetLineWidth(2);
246	g3.SetLineWidth(2);
247	g4.SetLineWidth(2);
248	gm.SetLineWidth(1);
249	g1.SetLineStyle(style[cnt/6]);
250	g2.SetLineStyle(style[cnt/6]);
251	g3.SetLineStyle(style[cnt/6]);
252	g4.SetLineStyle(style[cnt/6]);
253	g1.SetLineColor(color[cnt%6]);
254	g2.SetLineColor(color[cnt%6]);
255	g3.SetLineColor(color[cnt%6]);
256	g4.SetLineColor(color[cnt%6]);
257	gm.SetLineColor(kYellow);
258
259	if (cnt==0)
260	leg.AddEntry(gm.Clone(), "Moon", "l");
261	leg.AddEntry(g1.Clone(), name.data(), "l");
262
263	// Loop over 24 hours
264	int i=0;
265	for (double h=0; h<1; h+=1./(24*12), i++)
266	{
267	// check if it is betwene sun-rise and sun-set
268	if (h<jd0 \|\| h>jd1)
269	continue;
270
271	// get local position of source
272	const HrzPosn hrz = GetHrzFromEqu(pos, jd+h);
273
274	// Get sun and moon properties and
275	const EquPosn sun = fSunCoords[i];
276	const EquPosn moon = get<0>(fMoonCoords[i]);
277	const double disk = get<1>(fMoonCoords[i]);
278	const double edist = get<2>(fMoonCoords[i]);
279	const HrzPosn hrzm = GetHrzFromEqu(moon, jd+h);
280	const ZdAzPosn hrzs = GetHrzFromEqu(sun, jd+h);
281
282	if (v==res.begin())
283	cout << Time(jd+h) <<" " << 90-hrzm.alt << endl;
284
285	// Distance between source and moon
286	const double angle = GetAngularSeparation(moon, pos);
287
288	// Current prediction
289	const double sin_malt = hrzm.alt<0 ? 0 : sin(hrzm.alt*M_PI/180);
290	const double cos_mdist = cos(angle*M_PI/180);
291	const double sin_szd = sin(hrzs.zd*M_PI/180);
292
293	const double k0 = pow(disk, 2.63);
294	const double k1 = pow(sin_malt, 0.60);
295	const double k2 = pow(edist, -2.00);
296	const double k3 = exp(0.67cos_mdistcos_mdistcos_mdistcos_mdist);
297	const double k4 = exp(-97.8+105.8sin_szdsin_szd);
298
299	const double cur = 6.2 + 95.7k0k1k2k3 + k4;
300
301	// Relative energy threshold prediction
302	//const double cs = cos((90+hrz.alt)*M_PI/180);
303	//const double ratio = (10.sqrt(409600.cscs+9009.) + 6400.cs - 60.)/10.;
304	const double ratio = pow(cos((90-hrz.alt)*M_PI/180), -2.664);
305
306	// Add points to curve
307	const double axis = (mjd+h)243600;
308
309	// If there is a gap of more than one bin, start a new curve
310	CheckForGap(c1, g1, axis);
311	CheckForGap(c1, gm, axis);
312	CheckForGap(c2, g2, axis);
313	CheckForGap(c3, g3, axis);
314	CheckForGap(c4, g4, axis);
315
316	// Add data
317	if (no_limits \|\| cur<max_current)
318	g1.SetPoint(g1.GetN(), axis, hrz.alt);
319
320	if (no_limits \|\| 90-hrz.alt<max_zd)
321	g2.SetPoint(g2.GetN(), axis, cur);
322
323	if (no_limits \|\| (cur<max_current && 90-hrz.alt<max_zd))
324	g3.SetPoint(g3.GetN(), axis, ratio*pow(cur/6.2, 0.394));
325
326	if (no_limits \|\| (cur<max_current && 90-hrz.alt<max_zd))
327	g4.SetPoint(g4.GetN(), axis, angle);
328
329	if (cnt==0)
330	gm.SetPoint(gm.GetN(), axis, hrzm.alt);
331	}
332
333	if (cnt==0)
334	DrawClone(c1, gm);
335
336	DrawClone(c1, g1);
337	DrawClone(c2, g2);
338	DrawClone(c3, g3);
339	DrawClone(c4, g4);
340	}
341
342
343	// Save three plots
344	TCanvas c5;
345	c5.SetFillColor(kWhite);
346	c5.SetBorderMode(0);
347	c5.SetFrameBorderMode(0);
348	leg.Draw();
349
350	const string t = Time(jd).GetAsStr("%Y%m%d");
351
352	c1.SaveAs((t+"-ZenithDistance.eps").c_str());
353	c2.SaveAs((t+"-PredictedCurrent.eps").c_str());
354	c3.SaveAs((t+"-RelativeThreshold.eps").c_str());
355	c4.SaveAs((t+"-MoonDist.eps").c_str());
356	c5.SaveAs((t+"-Legend.eps").c_str());
357
358	c1.SaveAs((t+"-ZenithDistance.root").c_str());
359	c2.SaveAs((t+"-PredictedCurrent.root").c_str());
360	c3.SaveAs((t+"-RelativeThreshold.root").c_str());
361	c4.SaveAs((t+"-MoonDist.root").c_str());
362
363	c1.Print((t+".pdf(").c_str(), "pdf");
364	c2.Print((t+".pdf" ).c_str(), "pdf");
365	c3.Print((t+".pdf" ).c_str(), "pdf");
366	c4.Print((t+".pdf" ).c_str(), "pdf");
367	c5.Print((t+".pdf)").c_str(), "pdf");
368
369	return 0;
370	}

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

Download in other formats: