source: trunk/Mars/fact/analysis/lightcurve.C@ 16621

Last change on this file since 16621 was 15527, checked in by Daniela Dorner, 12 years ago
bugfix: sometimes last point was missing, adapted y-range for rates
File size: 13.0 KB
Line 
1#include <iostream>
2
3#include <TH1.h>
4#include <TH2.h>
5#include <TGraph.h>
6#include <TGraphErrors.h>
7
8#include <TFile.h>
9#include <TSystem.h>
10#include <TSQLResult.h>
11#include <TSQLRow.h>
12#include <TPRegexp.h>
13
14#include "MSQLMagic.h"
15#include "MAlphaFitter.h"
16#include "MHThetaSq.h"
17
18#include "MObservatory.h"
19#include "MAstro.h"
20#include "MAstroSky2Local.h"
21#include "MAstroCatalog.h"
22#include "MVector3.h"
23#include "MMath.h"
24#include "MDirIter.h"
25#include "MStatusArray.h"
26#include "MStatusDisplay.h"
27#include "MHCamera.h"
28#include "MSequence.h"
29
30using namespace std;
31
32
33int lightcurve(Int_t sourcekey=1, Int_t nightmin=0, Int_t nightmax=0, Int_t minutes=20, TString table="AnalysisResultsRunLP", TString outfile="lc")
34{
35
36 cout << minutes << " min bins..." << endl;
37 minutes*=60;
38 cout << minutes << " sec bins..." << endl;
39
40 MSQLServer serv("sql.rc");
41 if (!serv.IsConnected())
42 {
43 cout << "ERROR - Connection to database failed." << endl;
44 return 0;
45 }
46
47 TString query;
48 query=Form("SELECT Min(fNight), Max(fNight) FROM %s ", table.Data());
49 query+=" LEFT JOIN RunInfo USING(fNight, fRunId) ";
50 query+=Form(" WHERE fSourceKey=%d", sourcekey);
51 if (nightmin)
52 query+=Form(" AND fNight >=%d",nightmin);
53 if (nightmax)
54 query+=Form(" AND fNight <=%d",nightmax);
55
56 cout << "Q: " << query << endl;
57
58 TSQLResult *res = serv.Query(query);
59 if (!res)
60 return 1;
61
62 TSQLRow *row=res->Next();
63 TString nightstart=(*row)[0];
64 TString nightstop=(*row)[1];
65 delete res;
66
67 query =" SELECT fNight FROM "+table;
68 query+=" LEFT JOIN RunInfo USING(fNight, fRunId) ";
69 query+=Form(" WHERE fSourceKey=%d", sourcekey);
70 query+=Form(" AND fNight BETWEEN %s AND %s ", nightstart.Data(), nightstop.Data());
71 //query+=" AND NOT ISNULL(fOnTime) ";//at the moment filled at ISDC
72 query+=" GROUP BY fNight ";
73 query+=" ORDER BY fNight ";
74
75 cout << "Q: " << query << endl;
76
77
78 TSQLResult *res3 = serv.Query(query);
79 if (!res3)
80 return 1;
81
82 Int_t counter=0;
83 Int_t counter2=0;
84 TGraphErrors excessrate;
85 TGraphErrors backgroundrate;
86 TGraphErrors signif;
87 TGraphErrors significancerate;
88 TGraphErrors excratevsbgrate;
89 TGraphErrors excvszd;
90 TGraphErrors excvsbg;
91 TGraphErrors bgvszd;
92 TH1F exc("exc","rates", 43, -12.5, 202.5);
93 TH1F bg("bg","rates", 43, -12.5, 202.5);
94
95 Float_t significance;
96 TString excevts;
97 TString bgevts;
98 TString sigevts;
99 TString ontime;
100 TString zd;
101 TString zd2;
102 TString zd3;
103 Float_t excevtssum=0;
104 Float_t bgevtssum=0;
105 Float_t sigevtssum=0;
106 Float_t ontimesum=0;
107 Float_t excrate=0;
108 Float_t bgrate=0;
109 Float_t signifrate=0;
110 Float_t excerr=0;
111 Float_t zdsum=0;
112 Float_t zdmean=0;
113 Float_t zdmin=0;
114 Float_t zdmax=0;
115 MTime start;
116 MTime stop;
117 Double_t mjdstart;
118 Double_t mjdstop;
119 Double_t mjd;
120 TSQLRow *row3=0;
121 while ((row3=res3->Next()))
122 {
123 cout << "NIGHT " << (*row3)[0] << endl;
124 query =" SELECT IF(ISNULL(fOnTime), Time_to_sec(fRunStop)-Time_to_sec(fRunStart), fOnTime), ";
125 query+=" fNumBgEvts, fNumSigEvts, fNumExcEvts, ";
126 query+=" fRunStart, fRunStop, (fZenithDistanceMax-fZenithDistanceMin)/2, ";
127 query+=" fZenithDistanceMax, fZenithDistanceMin FROM "+table;
128 query+=" LEFT JOIN RunInfo USING(fNight, fRunId) ";
129 query+=Form(" WHERE fSourceKey=%d", sourcekey);
130 query+=Form(" AND fNight= %s ", (*row3)[0]);
131 //query+=" AND NOT ISNULL(fOnTime) ";//at the moment filled at isdc
132 query+=" ORDER BY fRunId ";
133
134 cout << "Q: " << query << endl;
135
136 TSQLResult *res2 = serv.Query(query);
137 if (!res2)
138 return 1;
139
140 counter2=0;
141 excevtssum=0;
142 bgevtssum=0;
143 sigevtssum=0;
144 ontimesum=0;
145 TSQLRow *row2=0;
146 while ((row2=res2->Next()))
147 {
148 counter2++;
149 cout << "-----" << (*row2)[0] << " " << (*row2)[3] << endl;
150 ontime=(*row2)[0];
151 bgevts=(*row2)[1];
152 sigevts=(*row2)[2];
153 excevts=(*row2)[3];
154 zd=(*row2)[6];
155 zd2=(*row2)[7];
156 zd3=(*row2)[8];
157 start.SetSqlDateTime((*row2)[4]);
158 stop.SetSqlDateTime((*row2)[5]);
159 if (counter2==1)
160 {
161 mjdstart=start.GetMjd();
162 cout << "mjdstart: " << (*row2)[4] << endl;
163 zdmin=zd3.Atof();
164 zdmax=zd2.Atof();
165 }
166
167 cout << "# " << counter2 << " " << (*row2)[0] << " " << (*row2)[1] << " " << (*row2)[2] << " " << (*row2)[3] << endl;
168
169 if (ontimesum+ontime.Atof()>minutes)
170 {
171 cout << "reached limit" << endl;
172 significance=MMath::SignificanceLiMa(sigevtssum, bgevtssum, 0.2);
173 cout << " significance: " << significance << endl;
174 cout << " excess: " << excevtssum << " evts" << endl;
175 cout << " bg: " << bgevtssum << " evts" << endl;
176 cout << " ontime: " << ontimesum << " s " << endl;
177 cout << " excessrate: " << excevtssum/ontimesum*3600 << " evts/h" << endl;
178 cout << " bgrate: " << bgevtssum/ontimesum*3600 << " evts/h" << endl;
179 signifrate=significance/sqrt(ontimesum*3600);
180 bgrate=bgevtssum/ontimesum*3600;
181 excrate=excevtssum/ontimesum*3600;
182 excerr = MMath::ErrorExc(significance, bgevtssum, 0.2)/ontimesum*3600.;
183 //zdmean=zdsum/counter2;
184 zdmean=zdmin+(zdmax-zdmin)/2;
185 mjdstop=stop.GetMjd();
186 mjd=mjdstart+(mjdstop-mjdstart)/2;
187 cout << " mjd: " << mjd << " (" << mjdstart << " - " << mjdstop << ") " << endl;
188 if (excrate>300)
189 cout << "==>" << (*row2)[0] << " " << (*row2)[1] << " " << (*row2)[2]
190 << " " << (*row2)[3] << " " << (*row2)[4] << " " << (*row2)[5]
191 << " " << (*row2)[6] << " " << (*row2)[7] << " " << (*row2)[8]
192 << endl;
193
194 excevtssum=0;
195 bgevtssum=0;
196 sigevtssum=0;
197 ontimesum=0;
198 zdsum=0;
199 exc.Fill(excrate);
200 bg.Fill(bgrate);
201 excessrate.SetPoint(counter, mjd, excrate);
202 excessrate.SetPointError(counter, ontimesum/3600./24/2, excerr);
203 backgroundrate.SetPoint(counter, mjd, bgrate);
204 backgroundrate.SetPointError(counter, ontimesum/3600./24/2, sqrt(bgevtssum)/ontimesum*3600);
205 excvszd.SetPoint(counter, zdmean, excrate);
206 excvszd.SetPointError(counter, (zdmax-zdmin)/2, excerr);
207 excvsbg.SetPoint(counter, bgrate, excrate);
208 excvsbg.SetPointError(counter, sqrt(bgevtssum)/ontimesum*3600, excerr);
209 bgvszd.SetPoint(counter, zdmean, bgrate);
210 bgvszd.SetPointError(counter, (zdmax-zdmin)/2, sqrt(bgevtssum)/ontimesum*3600);
211 signif.SetPoint(counter, mjd, significance);
212 signif.SetPointError(counter, ontimesum/3600./24/2, 0);
213 significancerate.SetPoint(counter, mjd, signifrate);
214 significancerate.SetPointError(counter, ontimesum/3600./24/2, 0);
215 counter++;
216 counter2=0;
217 }
218 ontimesum+=ontime.Atof();
219 excevtssum+=excevts.Atof();
220 bgevtssum+=bgevts.Atof();
221 sigevtssum+=sigevts.Atof();
222 zdsum+=zd.Atof();
223 if (zdmin>zd3.Atof())
224 zdmin=zd3.Atof();
225 if (zdmax<zd2.Atof())
226 zdmax=zd2.Atof();
227 //cout << "--> on " << ontimesum << " exc " << excevtssum << " bg " << bgevtssum << " sig " << sigevtssum << endl;
228
229 }
230 cout << "reached last run" << endl;
231 // if ontime is larger than 90% of the timebin width, the last point is filled
232 if (ontimesum>(minutes-0.1*minutes))
233 {
234 cout << "ontime > " << minutes-0.1*minutes << endl;
235 significance=MMath::SignificanceLiMa(sigevtssum, bgevtssum, 0.2);
236 cout << " significance: " << significance << endl;
237 cout << " excess: " << excevtssum << " evts" << endl;
238 cout << " bg: " << bgevtssum << " evts" << endl;
239 cout << " ontime: " << ontimesum << " s " << endl;
240 cout << " excessrate: " << excevtssum/ontimesum*3600 << " evts/h" << endl;
241 cout << " bgrate: " << bgevtssum/ontimesum*3600 << " evts/h" << endl;
242 signifrate=significance/sqrt(ontimesum*3600);
243 bgrate=bgevtssum/ontimesum*3600;
244 excrate=excevtssum/ontimesum*3600;
245 excerr = MMath::ErrorExc(significance, bgevtssum, 0.2)/ontimesum*3600.;
246 zdmean=zdmin+(zdmax-zdmin)/2;
247 mjdstop=stop.GetMjd();
248 mjd=mjdstart+(mjdstop-mjdstart)/2;
249 cout << " mjd: " << mjd << " (" << mjdstart << " - " << mjdstop << ") " << endl;
250
251 exc.Fill(excrate);
252 bg.Fill(bgrate);
253 excessrate.SetPoint(counter, mjd, excrate);
254 excessrate.SetPointError(counter, ontimesum/3600./24/2, excerr);
255 backgroundrate.SetPoint(counter, mjd, bgrate);
256 backgroundrate.SetPointError(counter, ontimesum/3600./24/2, sqrt(bgevtssum)/ontimesum*3600);
257 excvszd.SetPoint(counter, zdmean, excrate);
258 excvszd.SetPointError(counter, (zdmax-zdmin)/2, excerr);
259 excvsbg.SetPoint(counter, bgrate, excrate);
260 excvsbg.SetPointError(counter, sqrt(bgevtssum)/ontimesum*3600, excerr);
261 bgvszd.SetPoint(counter, zdmean, bgrate);
262 bgvszd.SetPointError(counter, (zdmax-zdmin)/2, sqrt(bgevtssum)/ontimesum*3600);
263 signif.SetPoint(counter, mjd, significance);
264 signif.SetPointError(counter, ontimesum/3600./24/2, 0);
265 significancerate.SetPoint(counter, mjd, signifrate);
266 significancerate.SetPointError(counter, ontimesum/3600./24/2, 0);
267 }
268 delete res2;
269 }
270 delete res3;
271
272 if (excessrate.GetN()<1)
273 return 1;
274
275 MStatusDisplay *d = new MStatusDisplay;
276
277 TCanvas &c1 = d->AddTab("Rates vs MJD", "Rates vs MJD");
278 gPad->SetGridy();
279 excessrate.SetTitle("Rates vs MJD");
280 excessrate.SetMarkerStyle(7);
281 excessrate.SetMarkerSize(2);
282 excessrate.GetYaxis()->SetTitle("evts/h");
283 excessrate.GetXaxis()->SetTitle("MJD");
284 excessrate.GetXaxis()->CenterTitle();
285 excessrate.GetYaxis()->SetRangeUser(-50,400);
286 excessrate.DrawClone("AP");
287
288 backgroundrate.SetMarkerStyle(7);
289 backgroundrate.SetMarkerSize(2);
290 backgroundrate.SetMarkerColor(kBlue);
291 backgroundrate.SetLineColor(kBlue);
292 backgroundrate.DrawClone("Psame");
293
294 TCanvas &c2 = d->AddTab("Excess Rate vs MJD", "Excess Rate vs MJD");
295 gPad->SetGridy();
296 excessrate.SetTitle("Excess Rate vs MJD");
297 excessrate.SetMarkerStyle(7);
298 excessrate.SetMarkerSize(2);
299 excessrate.GetYaxis()->SetTitle("exc evts/h");
300 excessrate.GetXaxis()->SetTitle("MJD");
301 excessrate.GetXaxis()->CenterTitle();
302 excessrate.DrawClone("AP");
303
304 TCanvas &c3 = d->AddTab("Exc vs Zd", "Exc vs Zd");
305 excvszd.SetTitle("Excess Rate vs Zd");
306 excvszd.SetMarkerStyle(7);
307 excvszd.SetMarkerSize(2);
308 excvszd.GetYaxis()->SetTitle("exc evts/h");
309 excvszd.GetXaxis()->SetTitle("zd");
310 excvszd.GetXaxis()->CenterTitle();
311 excvszd.DrawClone("AP");
312
313 TCanvas &c4 = d->AddTab("Exc vs Bg", "Exc vs Bg");
314 excvsbg.SetTitle("Excess Rate vs Background Rate");
315 excvsbg.SetMarkerStyle(7);
316 excvsbg.SetMarkerSize(2);
317 excvsbg.GetYaxis()->SetTitle("exc evts/h");
318 excvsbg.GetXaxis()->SetTitle("bg evts/h");
319 excvsbg.GetXaxis()->CenterTitle();
320 excvsbg.DrawClone("AP");
321
322 TCanvas &c5 = d->AddTab("Bg vs Zd", "Bg vs Zd");
323 bgvszd.SetTitle("Background Rate vs Zd");
324 bgvszd.SetMarkerStyle(7);
325 bgvszd.SetMarkerSize(2);
326 bgvszd.GetYaxis()->SetTitle("bg evts/h");
327 bgvszd.GetXaxis()->SetTitle("zd");
328 bgvszd.GetXaxis()->CenterTitle();
329 bgvszd.DrawClone("AP");
330
331 TCanvas &c6 = d->AddTab("Signifrate vs MJD", "Signifrate vs MJD");
332 gPad->SetGridy();
333 significancerate.SetTitle("Significance Rate vs MJD");
334 significancerate.SetMarkerStyle(7);
335 significancerate.SetMarkerSize(2);
336 significancerate.GetYaxis()->SetTitle("sigma/sqrt(h)");
337 significancerate.GetXaxis()->SetTitle("MJD");
338 significancerate.GetXaxis()->CenterTitle();
339 significancerate.DrawClone("AP");
340
341 TCanvas &c7 = d->AddTab("Significance vs MJD", "Significance vs MJD");
342 signif.SetTitle("Significance vs MJD");
343 signif.SetMarkerStyle(7);
344 signif.SetMarkerSize(2);
345 signif.GetYaxis()->SetTitle("sigma");
346 signif.GetXaxis()->SetTitle("MJD");
347 signif.GetXaxis()->CenterTitle();
348 signif.DrawClone("AP");
349
350 TCanvas &c8 = d->AddTab("Rate Hist", "Rate Hist");
351 bg.SetLineColor(kBlue);
352 bg.DrawCopy("hist");
353 exc.DrawCopy("histsame");
354
355 d->SaveAs(outfile);
356
357 return 0;
358
359}
360
361
362
Note: See TracBrowser for help on using the repository browser.