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

Last change on this file since 15515 was 15396, checked in by Daniela Dorner, 13 years ago
adapted query to allow for plots even when the ontime is not yet inserted in the DB
File size: 10.6 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 ontime=(*row2)[0];
150 bgevts=(*row2)[1];
151 sigevts=(*row2)[2];
152 excevts=(*row2)[3];
153 zd=(*row2)[6];
154 zd2=(*row2)[7];
155 zd3=(*row2)[8];
156 start.SetSqlDateTime((*row2)[4]);
157 stop.SetSqlDateTime((*row2)[5]);
158 if (counter2==1)
159 {
160 mjdstart=start.GetMjd();
161 cout << "mjdstart: " << (*row2)[4] << endl;
162 zdmin=zd3.Atof();
163 zdmax=zd2.Atof();
164 }
165
166 cout << "# " << counter2 << " " << (*row2)[0] << " " << (*row2)[1] << " " << (*row2)[2] << " " << (*row2)[3] << endl;
167
168 if (ontimesum+ontime.Atof()>minutes)
169 {
170 cout << "reached limit" << endl;
171 significance=MMath::SignificanceLiMa(sigevtssum, bgevtssum, 0.2);
172 cout << " significance: " << significance << endl;
173 cout << " excess: " << excevtssum << endl;
174 cout << " bg: " << bgevtssum << endl;
175 cout << " ontime: " << ontimesum/3600 << endl;
176 cout << " excessrate: " << excevtssum/ontimesum*3600 << endl;
177 cout << " bgrate: " << bgevtssum/ontimesum*3600 << endl;
178 signifrate=significance/ontimesum*3600;
179 bgrate=bgevtssum/ontimesum*3600;
180 excrate=excevtssum/ontimesum*3600;
181 excerr = MMath::ErrorExc(significance, bgevtssum, 0.2)/ontimesum*3600.;
182 //zdmean=zdsum/counter2;
183 zdmean=zdmin+(zdmax-zdmin)/2;
184 mjdstop=stop.GetMjd();
185 mjd=mjdstart+(mjdstop-mjdstart)/2;
186 cout << " mjd: " << mjd << " (" << mjdstart << " - " << mjdstop << ") " << endl;
187 if (mjd<50000)
188 cout << "==>" << (*row2)[0] << " " << (*row2)[1] << " " << (*row2)[2]
189 << " " << (*row2)[3] << " " << (*row2)[4] << " " << (*row2)[5]
190 << " " << (*row2)[6] << " " << (*row2)[7] << " " << (*row2)[8]
191 << endl;
192
193 excevtssum=0;
194 bgevtssum=0;
195 sigevtssum=0;
196 ontimesum=0;
197 zdsum=0;
198 exc.Fill(excrate);
199 bg.Fill(bgrate);
200 excessrate.SetPoint(counter, mjd, excrate);
201 excessrate.SetPointError(counter, ontimesum/3600./24/2, excerr);
202 backgroundrate.SetPoint(counter, mjd, bgrate);
203 backgroundrate.SetPointError(counter, ontimesum/3600./24/2, sqrt(bgevtssum)/ontimesum*3600);
204 excvszd.SetPoint(counter, zdmean, excrate);
205 excvszd.SetPointError(counter, (zdmax-zdmin)/2, excerr);
206 excvsbg.SetPoint(counter, bgrate, excrate);
207 excvsbg.SetPointError(counter, sqrt(bgevtssum)/ontimesum*3600, excerr);
208 bgvszd.SetPoint(counter, zdmean, bgrate);
209 bgvszd.SetPointError(counter, (zdmax-zdmin)/2, sqrt(bgevtssum)/ontimesum*3600);
210 signif.SetPoint(counter, mjd, significance);
211 signif.SetPointError(counter, ontimesum/3600./24/2, 0);
212 significancerate.SetPoint(counter, mjd, signifrate);
213 significancerate.SetPointError(counter, ontimesum/3600./24/2, 0);
214 counter++;
215 counter2=0;
216 }
217 ontimesum+=ontime.Atof();
218 excevtssum+=excevts.Atof();
219 bgevtssum+=bgevts.Atof();
220 sigevtssum+=sigevts.Atof();
221 zdsum+=zd.Atof();
222 if (zdmin>zd3.Atof())
223 zdmin=zd3.Atof();
224 if (zdmax<zd2.Atof())
225 zdmax=zd2.Atof();
226 //cout << "--> on " << ontimesum << " exc " << excevtssum << " bg " << bgevtssum << " sig " << sigevtssum << endl;
227
228 }
229
230 delete res2;
231 }
232 delete res3;
233
234 if (excessrate.GetN()<1)
235 return 1;
236
237 MStatusDisplay *d = new MStatusDisplay;
238
239 TCanvas &c1 = d->AddTab("Rates vs MJD", "Rates vs MJD");
240 excessrate.SetTitle("Rates vs MJD");
241 excessrate.SetMarkerStyle(7);
242 excessrate.SetMarkerSize(2);
243 excessrate.GetYaxis()->SetTitle("evts/h");
244 excessrate.GetXaxis()->SetTitle("MJD");
245 excessrate.GetXaxis()->CenterTitle();
246 excessrate.GetYaxis()->SetRangeUser(-50,200);
247 excessrate.DrawClone("AP");
248
249 backgroundrate.SetMarkerStyle(7);
250 backgroundrate.SetMarkerSize(2);
251 backgroundrate.SetMarkerColor(kBlue);
252 backgroundrate.SetLineColor(kBlue);
253 backgroundrate.DrawClone("Psame");
254
255 TCanvas &c2 = d->AddTab("Excess Rate vs MJD", "Excess Rate vs MJD");
256 excessrate.SetTitle("Excess Rate vs MJD");
257 excessrate.SetMarkerStyle(7);
258 excessrate.SetMarkerSize(2);
259 excessrate.GetYaxis()->SetTitle("exc evts/h");
260 excessrate.GetXaxis()->SetTitle("MJD");
261 excessrate.GetXaxis()->CenterTitle();
262 excessrate.DrawClone("AP");
263
264 TCanvas &c3 = d->AddTab("Exc vs Zd", "Exc vs Zd");
265 excvszd.SetTitle("Excess Rate vs Zd");
266 excvszd.SetMarkerStyle(7);
267 excvszd.SetMarkerSize(2);
268 excvszd.GetYaxis()->SetTitle("exc evts/h");
269 excvszd.GetXaxis()->SetTitle("zd");
270 excvszd.GetXaxis()->CenterTitle();
271 excvszd.DrawClone("AP");
272
273 TCanvas &c4 = d->AddTab("Exc vs Bg", "Exc vs Bg");
274 excvsbg.SetTitle("Excess Rate vs Background Rate");
275 excvsbg.SetMarkerStyle(7);
276 excvsbg.SetMarkerSize(2);
277 excvsbg.GetYaxis()->SetTitle("exc evts/h");
278 excvsbg.GetXaxis()->SetTitle("bg evts/h");
279 excvsbg.GetXaxis()->CenterTitle();
280 excvsbg.DrawClone("AP");
281
282 TCanvas &c5 = d->AddTab("Bg vs Zd", "Bg vs Zd");
283 bgvszd.SetTitle("Background Rate vs Zd");
284 bgvszd.SetMarkerStyle(7);
285 bgvszd.SetMarkerSize(2);
286 bgvszd.GetYaxis()->SetTitle("bg evts/h");
287 bgvszd.GetXaxis()->SetTitle("zd");
288 bgvszd.GetXaxis()->CenterTitle();
289 bgvszd.DrawClone("AP");
290
291 TCanvas &c6 = d->AddTab("Signifrate vs MJD", "Signifrate vs MJD");
292 significancerate.SetTitle("Significance Rate vs MJD");
293 significancerate.SetMarkerStyle(7);
294 significancerate.SetMarkerSize(2);
295 significancerate.GetYaxis()->SetTitle("sigma/sqrt(h)");
296 significancerate.GetXaxis()->SetTitle("MJD");
297 significancerate.GetXaxis()->CenterTitle();
298 significancerate.DrawClone("AP");
299
300 TCanvas &c7 = d->AddTab("Significance vs MJD", "Significance vs MJD");
301 signif.SetTitle("Significance vs MJD");
302 signif.SetMarkerStyle(7);
303 signif.SetMarkerSize(2);
304 signif.GetYaxis()->SetTitle("sigma");
305 signif.GetXaxis()->SetTitle("MJD");
306 signif.GetXaxis()->CenterTitle();
307 signif.DrawClone("AP");
308
309 TCanvas &c8 = d->AddTab("Rate Hist", "Rate Hist");
310 bg.SetLineColor(kBlue);
311 bg.DrawCopy("hist");
312 exc.DrawCopy("histsame");
313
314 d->SaveAs(outfile);
315
316 return 0;
317
318}
319
320
321
Note: See TracBrowser for help on using the repository browser.