#include <iostream>

#include <TH1.h>
#include <TH2.h>
#include <TGraph.h>
#include <TGraphErrors.h>

#include <TFile.h>
#include <TSystem.h>
#include <TSQLResult.h>
#include <TSQLRow.h>
#include <TPRegexp.h>

#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) ";//at the moment filled at ISDC
    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(fOnTimeAfterCuts), IF(ISNULL(fOnTime), Time_to_sec(Timediff(fRunStop,fRunStart)), fOnTime), fOnTimeAfterCuts), ";
        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) ";//at the moment filled at isdc
        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++;
            cout << "-----" << (*row2)[0] << " " << (*row2)[3] << endl;
            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 << " evts" << endl;
                cout << "  bg: " << bgevtssum << " evts" << endl;
                cout << "  ontime: " << ontimesum << " s " << endl;
                cout << "  excessrate: " << excevtssum/ontimesum*3600 << " evts/h" << endl;
                cout << "  bgrate: " << bgevtssum/ontimesum*3600 << " evts/h" << endl;
                signifrate=significance/sqrt(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 (excrate>300)
                    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<zd2.Atof())
                zdmax=zd2.Atof();
            //cout << "--> on " << ontimesum << " exc " << excevtssum << " bg " << bgevtssum << " sig " << sigevtssum << endl;

        }
        cout << "reached last run" << endl;
        // if ontime is larger than 90% of the timebin width, the last point is filled
        if (ontimesum>(minutes-0.1*minutes))
        {
            cout << "ontime > " << minutes-0.1*minutes << endl;
            significance=MMath::SignificanceLiMa(sigevtssum, bgevtssum, 0.2);
            cout << "  significance: " << significance << endl;
            cout << "  excess: " << excevtssum << " evts" << endl;
            cout << "  bg: " << bgevtssum << " evts" << endl;
            cout << "  ontime: " << ontimesum << " s " << endl;
            cout << "  excessrate: " << excevtssum/ontimesum*3600 << " evts/h" << endl;
            cout << "  bgrate: " << bgevtssum/ontimesum*3600 << " evts/h" << endl;
            signifrate=significance/sqrt(ontimesum*3600);
            bgrate=bgevtssum/ontimesum*3600;
            excrate=excevtssum/ontimesum*3600;
            excerr = MMath::ErrorExc(significance, bgevtssum, 0.2)/ontimesum*3600.;
            zdmean=zdmin+(zdmax-zdmin)/2;
            mjdstop=stop.GetMjd();
            mjd=mjdstart+(mjdstop-mjdstart)/2;
            cout << "  mjd: " << mjd <<  " (" << mjdstart << " - " << mjdstop << ") " << endl;

            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);
        }
        delete res2;
    }
    delete res3;

    if (excessrate.GetN()<1)
        return 1;

    MStatusDisplay *d = new MStatusDisplay;

    TCanvas &c1 = d->AddTab("Rates vs MJD", "Rates vs MJD");
    gPad->SetGridy();
    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,400);
    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");
    gPad->SetGridy();
    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");
    gPad->SetGridy();
    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;

}