Index: /trunk/Mars/fact/analysis/lightcurve.C =================================================================== --- /trunk/Mars/fact/analysis/lightcurve.C (revision 15395) +++ /trunk/Mars/fact/analysis/lightcurve.C (revision 15395) @@ -0,0 +1,321 @@ +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "MSQLMagic.h" +#include "MAlphaFitter.h" +#include "MHThetaSq.h" + +#include "MObservatory.h" +#include "MAstro.h" +#include "MAstroSky2Local.h" +#include "MAstroCatalog.h" +#include "MVector3.h" +#include "MMath.h" +#include "MDirIter.h" +#include "MStatusArray.h" +#include "MStatusDisplay.h" +#include "MHCamera.h" +#include "MSequence.h" + +using namespace std; + + +int lightcurve(Int_t sourcekey=1, Int_t nightmin=0, Int_t nightmax=0, Int_t minutes=20, TString table="AnalysisResultsRunLP", TString outfile="lc") +{ + + cout << minutes << " min bins..." << endl; + minutes*=60; + cout << minutes << " sec bins..." << endl; + + MSQLServer serv("sql.rc"); + if (!serv.IsConnected()) + { + cout << "ERROR - Connection to database failed." << endl; + return 0; + } + + TString query; + query=Form("SELECT Min(fNight), Max(fNight) FROM %s ", table.Data()); + query+=" LEFT JOIN RunInfo USING(fNight, fRunId) "; + query+=Form(" WHERE fSourceKey=%d", sourcekey); + if (nightmin) + query+=Form(" AND fNight >=%d",nightmin); + if (nightmax) + query+=Form(" AND fNight <=%d",nightmax); + + cout << "Q: " << query << endl; + + TSQLResult *res = serv.Query(query); + if (!res) + return 1; + + TSQLRow *row=res->Next(); + TString nightstart=(*row)[0]; + TString nightstop=(*row)[1]; + delete res; + + query =" SELECT fNight FROM "+table; + query+=" LEFT JOIN RunInfo USING(fNight, fRunId) "; + query+=Form(" WHERE fSourceKey=%d", sourcekey); + query+=Form(" AND fNight BETWEEN %s AND %s ", nightstart.Data(), nightstop.Data()); + query+=" AND NOT ISNULL(fOnTime) ";//to be checked why sometimes NULL + query+=" GROUP BY fNight "; + query+=" ORDER BY fNight "; + + cout << "Q: " << query << endl; + + + TSQLResult *res3 = serv.Query(query); + if (!res3) + return 1; + + Int_t counter=0; + Int_t counter2=0; + TGraphErrors excessrate; + TGraphErrors backgroundrate; + TGraphErrors signif; + TGraphErrors significancerate; + TGraphErrors excratevsbgrate; + TGraphErrors excvszd; + TGraphErrors excvsbg; + TGraphErrors bgvszd; + TH1F exc("exc","rates", 43, -12.5, 202.5); + TH1F bg("bg","rates", 43, -12.5, 202.5); + + Float_t significance; + TString excevts; + TString bgevts; + TString sigevts; + TString ontime; + TString zd; + TString zd2; + TString zd3; + Float_t excevtssum=0; + Float_t bgevtssum=0; + Float_t sigevtssum=0; + Float_t ontimesum=0; + Float_t excrate=0; + Float_t bgrate=0; + Float_t signifrate=0; + Float_t excerr=0; + Float_t zdsum=0; + Float_t zdmean=0; + Float_t zdmin=0; + Float_t zdmax=0; + MTime start; + MTime stop; + Double_t mjdstart; + Double_t mjdstop; + Double_t mjd; + TSQLRow *row3=0; + while ((row3=res3->Next())) + { + cout << "NIGHT " << (*row3)[0] << endl; + query =" SELECT IF(ISNULL(fOnTime), Time_to_sec(fRunStop)-Time_to_sec(fRunStart), fOnTime), "; + query+=" fNumBgEvts, fNumSigEvts, fNumExcEvts, "; + query+=" fRunStart, fRunStop, (fZenithDistanceMax-fZenithDistanceMin)/2, "; + query+=" fZenithDistanceMax, fZenithDistanceMin FROM "+table; + query+=" LEFT JOIN RunInfo USING(fNight, fRunId) "; + query+=Form(" WHERE fSourceKey=%d", sourcekey); + query+=Form(" AND fNight= %s ", (*row3)[0]); + query+=" AND NOT ISNULL(fOnTime) ";//to be checked why sometimes NULL + query+=" ORDER BY fRunId "; + + cout << "Q: " << query << endl; + + TSQLResult *res2 = serv.Query(query); + if (!res2) + return 1; + + counter2=0; + excevtssum=0; + bgevtssum=0; + sigevtssum=0; + ontimesum=0; + TSQLRow *row2=0; + while ((row2=res2->Next())) + { + counter2++; + ontime=(*row2)[0]; + bgevts=(*row2)[1]; + sigevts=(*row2)[2]; + excevts=(*row2)[3]; + zd=(*row2)[6]; + zd2=(*row2)[7]; + zd3=(*row2)[8]; + start.SetSqlDateTime((*row2)[4]); + stop.SetSqlDateTime((*row2)[5]); + if (counter2==1) + { + mjdstart=start.GetMjd(); + cout << "mjdstart: " << (*row2)[4] << endl; + zdmin=zd3.Atof(); + zdmax=zd2.Atof(); + } + + cout << "# " << counter2 << " " << (*row2)[0] << " " << (*row2)[1] << " " << (*row2)[2] << " " << (*row2)[3] << endl; + + if (ontimesum+ontime.Atof()>minutes) + { + cout << "reached limit" << endl; + significance=MMath::SignificanceLiMa(sigevtssum, bgevtssum, 0.2); + cout << " significance: " << significance << endl; + cout << " excess: " << excevtssum << endl; + cout << " bg: " << bgevtssum << endl; + cout << " ontime: " << ontimesum/3600 << endl; + cout << " excessrate: " << excevtssum/ontimesum*3600 << endl; + cout << " bgrate: " << bgevtssum/ontimesum*3600 << endl; + signifrate=significance/ontimesum*3600; + bgrate=bgevtssum/ontimesum*3600; + excrate=excevtssum/ontimesum*3600; + excerr = MMath::ErrorExc(significance, bgevtssum, 0.2)/ontimesum*3600.; + //zdmean=zdsum/counter2; + zdmean=zdmin+(zdmax-zdmin)/2; + mjdstop=stop.GetMjd(); + mjd=mjdstart+(mjdstop-mjdstart)/2; + cout << " mjd: " << mjd << " (" << mjdstart << " - " << mjdstop << ") " << endl; + if (mjd<50000) + cout << "==>" << (*row2)[0] << " " << (*row2)[1] << " " << (*row2)[2] + << " " << (*row2)[3] << " " << (*row2)[4] << " " << (*row2)[5] + << " " << (*row2)[6] << " " << (*row2)[7] << " " << (*row2)[8] + << endl; + + excevtssum=0; + bgevtssum=0; + sigevtssum=0; + ontimesum=0; + zdsum=0; + exc.Fill(excrate); + bg.Fill(bgrate); + excessrate.SetPoint(counter, mjd, excrate); + excessrate.SetPointError(counter, ontimesum/3600./24/2, excerr); + backgroundrate.SetPoint(counter, mjd, bgrate); + backgroundrate.SetPointError(counter, ontimesum/3600./24/2, sqrt(bgevtssum)/ontimesum*3600); + excvszd.SetPoint(counter, zdmean, excrate); + excvszd.SetPointError(counter, (zdmax-zdmin)/2, excerr); + excvsbg.SetPoint(counter, bgrate, excrate); + excvsbg.SetPointError(counter, sqrt(bgevtssum)/ontimesum*3600, excerr); + bgvszd.SetPoint(counter, zdmean, bgrate); + bgvszd.SetPointError(counter, (zdmax-zdmin)/2, sqrt(bgevtssum)/ontimesum*3600); + signif.SetPoint(counter, mjd, significance); + signif.SetPointError(counter, ontimesum/3600./24/2, 0); + significancerate.SetPoint(counter, mjd, signifrate); + significancerate.SetPointError(counter, ontimesum/3600./24/2, 0); + counter++; + counter2=0; + } + ontimesum+=ontime.Atof(); + excevtssum+=excevts.Atof(); + bgevtssum+=bgevts.Atof(); + sigevtssum+=sigevts.Atof(); + zdsum+=zd.Atof(); + if (zdmin>zd3.Atof()) + zdmin=zd3.Atof(); + if (zdmax on " << ontimesum << " exc " << excevtssum << " bg " << bgevtssum << " sig " << sigevtssum << endl; + + } + + delete res2; + } + delete res3; + + if (excessrate.GetN()<1) + return 1; + + MStatusDisplay *d = new MStatusDisplay; + + TCanvas &c1 = d->AddTab("Rates vs MJD", "Rates vs MJD"); + excessrate.SetTitle("Rates vs MJD"); + excessrate.SetMarkerStyle(7); + excessrate.SetMarkerSize(2); + excessrate.GetYaxis()->SetTitle("evts/h"); + excessrate.GetXaxis()->SetTitle("MJD"); + excessrate.GetXaxis()->CenterTitle(); + excessrate.GetYaxis()->SetRangeUser(-50,200); + excessrate.DrawClone("AP"); + + backgroundrate.SetMarkerStyle(7); + backgroundrate.SetMarkerSize(2); + backgroundrate.SetMarkerColor(kBlue); + backgroundrate.SetLineColor(kBlue); + backgroundrate.DrawClone("Psame"); + + TCanvas &c2 = d->AddTab("Excess Rate vs MJD", "Excess Rate vs MJD"); + excessrate.SetTitle("Excess Rate vs MJD"); + excessrate.SetMarkerStyle(7); + excessrate.SetMarkerSize(2); + excessrate.GetYaxis()->SetTitle("exc evts/h"); + excessrate.GetXaxis()->SetTitle("MJD"); + excessrate.GetXaxis()->CenterTitle(); + excessrate.DrawClone("AP"); + + TCanvas &c3 = d->AddTab("Exc vs Zd", "Exc vs Zd"); + excvszd.SetTitle("Excess Rate vs Zd"); + excvszd.SetMarkerStyle(7); + excvszd.SetMarkerSize(2); + excvszd.GetYaxis()->SetTitle("exc evts/h"); + excvszd.GetXaxis()->SetTitle("zd"); + excvszd.GetXaxis()->CenterTitle(); + excvszd.DrawClone("AP"); + + TCanvas &c4 = d->AddTab("Exc vs Bg", "Exc vs Bg"); + excvsbg.SetTitle("Excess Rate vs Background Rate"); + excvsbg.SetMarkerStyle(7); + excvsbg.SetMarkerSize(2); + excvsbg.GetYaxis()->SetTitle("exc evts/h"); + excvsbg.GetXaxis()->SetTitle("bg evts/h"); + excvsbg.GetXaxis()->CenterTitle(); + excvsbg.DrawClone("AP"); + + TCanvas &c5 = d->AddTab("Bg vs Zd", "Bg vs Zd"); + bgvszd.SetTitle("Background Rate vs Zd"); + bgvszd.SetMarkerStyle(7); + bgvszd.SetMarkerSize(2); + bgvszd.GetYaxis()->SetTitle("bg evts/h"); + bgvszd.GetXaxis()->SetTitle("zd"); + bgvszd.GetXaxis()->CenterTitle(); + bgvszd.DrawClone("AP"); + + TCanvas &c6 = d->AddTab("Signifrate vs MJD", "Signifrate vs MJD"); + significancerate.SetTitle("Significance Rate vs MJD"); + significancerate.SetMarkerStyle(7); + significancerate.SetMarkerSize(2); + significancerate.GetYaxis()->SetTitle("sigma/sqrt(h)"); + significancerate.GetXaxis()->SetTitle("MJD"); + significancerate.GetXaxis()->CenterTitle(); + significancerate.DrawClone("AP"); + + TCanvas &c7 = d->AddTab("Significance vs MJD", "Significance vs MJD"); + signif.SetTitle("Significance vs MJD"); + signif.SetMarkerStyle(7); + signif.SetMarkerSize(2); + signif.GetYaxis()->SetTitle("sigma"); + signif.GetXaxis()->SetTitle("MJD"); + signif.GetXaxis()->CenterTitle(); + signif.DrawClone("AP"); + + TCanvas &c8 = d->AddTab("Rate Hist", "Rate Hist"); + bg.SetLineColor(kBlue); + bg.DrawCopy("hist"); + exc.DrawCopy("histsame"); + + d->SaveAs(outfile); + + return 0; + +} + + +