Bool_t Contains(TArrayD **vec, Double_t t0, Double_t range=0)
{
    TArrayD *arr0 = vec[0];
    TArrayD *arr1 = vec[1];
    TArrayD *arr2 = vec[2];

    for (int i=0; i<arr0->GetSize(); i++)
    {
        Double_t start = (*arr1)[i];
        Double_t stop  = (*arr2)[i];

        if (stop>start+305./24/3600)
            stop = start+305./24/3600;

        if (t0>start-range && t0<stop+range)
        //{
        //    if (fmod(t0, 1)>4./24 && fmod(t0,1)<4.1/24)
        //        cout << t0-start << " " <<t0 << " " << stop-t0 << " " << start-15779 << " " << stop-15779 << " " << (*arr0)[i] << endl;
            return kTRUE;
        //}
    }

    return arr0->GetSize()==0;
}

Int_t PlotThresholds(TArrayD **vec, TString fname)
{
    fname += ".RATE_CONTROL_THRESHOLD.fits";

    fits file(fname.Data());
    if (!file)
    {
        cerr << fname << ": " << gSystem->GetError() << endl;
        return -2;
    }

    //cout << fname << endl;

    Double_t time;
    UShort_t th;

    if (!file.SetPtrAddress("Time", &time))
        return -1;

    if (!file.SetPtrAddress("threshold", &th))
        return -1;

    TGraph g;

    while (file.GetNextRow())
    {
        if (Contains(vec, time, 10./(24*3600)))
            g.SetPoint(g.GetN(), time*24*3600, th);
    }

    g.SetMarkerStyle(kFullDotMedium);
    g.GetXaxis()->SetTimeDisplay(true);
    g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    g.GetXaxis()->SetLabelSize(0.12);
    g.GetYaxis()->SetLabelSize(0.1);
    g.GetYaxis()->SetTitle("THRESHOLD");
    g.GetYaxis()->SetTitleOffset(0.2);
    g.GetYaxis()->SetTitleSize(0.1);
    g.DrawClone("AP");

    return 0;
}

#include <vector>
#include <pair>

vector<pair<double, ln_equ_posn>> vec;

ln_equ_posn FindPointing(Double_t time)
{
    for (int i=0; i<vec.size(); i++)
        if (time<vec[i].first)
            return i==0 ? ln_equ_posn() : vec[i-1].second;

    return vec[vec.size()-1].second;
}

Float_t Prediction(Double_t time)
{
    Double_t jd = time + 40587 + 2400000.5;

    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);

    // Nova::EquPosn pos = FindPointing(time);

    // get local position of moon
    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);

    // Distance between source and moon
    //const double angle =
    //    MAstro::AngularDistance(TMath::DegToRad()*(90-moon.dec),
    //                            TMath::DegToRad()*moon.ra,
    //                            TMath::DegToRad()*(90-pos.dec),
    //                            TMath::DegToRad()*pos.ra);

    // Current prediction
    //double cang = sin(angle);
    double calt = sin(hrzm.alt*TMath::DegToRad());

    double disk = Nova::GetLunarDisk(jd);

    double lc = /*cang==0 ? 0 :*/ calt*pow(disk, 2.5);///sqrt(cang);

    return 5+103*lc>4.5 ? 5+103*lc : 4.5;
}

Int_t ReadSources(TString fname)
{
    fname += ".DRIVE_CONTROL_SOURCE_POSITION.fits";

    fits file(fname.Data());
    if (!file)
    {
        cerr << fname << ": " << gSystem->GetError() << endl;
        return -2;
    }

    Double_t time, ra, dec;
    if (!file.SetPtrAddress("Time",  &time))
        return -1;
    if (!file.SetPtrAddress("Ra_cmd", &ra))
        return -1;
    if (!file.SetPtrAddress("Dec_cmd", &dec))
        return -1;

    while (file.GetNextRow())
    {
        ln_equ_posn p;
        p.ra  = ra*15;
        p.dec = dec;

        vec.push_back(make_pair(time, p));
    }

    return 0;
}

Int_t PlotCurrent(TArrayD **vec, TString fname)
{
    Int_t rc = ReadSources(fname);
    if (rc<0)
        return rc;

    fname += ".FEEDBACK_CALIBRATED_CURRENTS.fits";

    fits file(fname.Data());
    if (!file)
    {
        cerr << fname << ": " << gSystem->GetError() << endl;
        return -2;
    }

    //cout << fname << endl;

    Double_t time;
    Float_t  Imed;

    if (!file.SetPtrAddress("Time",  &time))
        return -1;

    if (!file.SetPtrAddress("I_med", &Imed))
        return -1;

    TGraph g1;
    TGraph g2;

    while (file.GetNextRow())
        if (Contains(vec, time))
        {
            g1.SetPoint(g1.GetN(), time*24*3600, Imed);
            g2.SetPoint(g2.GetN(), time*24*3600, Prediction(time));
        }

    g1.SetMinimum(0);
    g1.SetMaximum(99);

    g1.SetMarkerStyle(kFullDotMedium);
    g1.GetXaxis()->SetTimeDisplay(true);
    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    g1.GetXaxis()->SetLabelSize(0.12);
    g1.GetYaxis()->SetLabelSize(0.1);
    g1.GetYaxis()->SetTitle("CURRENT");
    g1.GetYaxis()->SetTitleOffset(0.2);
    g1.GetYaxis()->SetTitleSize(0.1);
    g1.DrawClone("AP");

    g2.SetMarkerColor(kBlue);
    g2.SetMarkerStyle(kFullDotMedium);
    g2.GetXaxis()->SetTimeDisplay(true);
    g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    g2.GetXaxis()->SetLabelSize(0.12);
    g2.GetYaxis()->SetLabelSize(0.1);
    g2.GetYaxis()->SetTitle("CURRENT");
    g2.GetYaxis()->SetTitleOffset(0.2);
    g2.GetYaxis()->SetTitleSize(0.1);
    g2.DrawClone("P");

    g1.DrawClone("P");

    return 0;
}

Int_t PlotRate(TArrayD **vec, TString fname)
{
    fname += ".FTM_CONTROL_TRIGGER_RATES.fits";

    fits file(fname.Data());
    if (!file)
    {
        cerr << fname << ": " << gSystem->GetError() << endl;
        return -2;
    }

    //cout << fname << endl;

    Double_t time;
    Float_t  rate;
    Float_t  ontime, elapsed;

    if (!file.SetPtrAddress("Time",  &time))
        return -1;

    if (!file.SetPtrAddress("TriggerRate", &rate))
        return -1;
    if (!file.SetPtrAddress("OnTime", &ontime))
        return -1;
    if (!file.SetPtrAddress("ElapsedTime", &elapsed))
        return -1;

    TGraph g1, g2;

    while (file.GetNextRow())
        if (Contains(vec, time))
        {
            if (rate>=0)
            {
                g1.SetPoint(g1.GetN(), time*24*3600, rate);
                g2.SetPoint(g2.GetN(), time*24*3600, ontime/elapsed);
            }
        }

    g1.SetMinimum(0);
    g1.SetMaximum(200);
    g1.SetMarkerStyle(kFullDotMedium);
    g1.GetXaxis()->SetTimeDisplay(true);
    g1.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    g1.GetXaxis()->SetLabelSize(0.12);
    g1.GetYaxis()->SetLabelSize(0.1);
    g1.GetYaxis()->SetTitle("TRIGGER RATE");
    g1.GetYaxis()->SetTitleOffset(0.2);
    g1.GetYaxis()->SetTitleSize(0.1);
    g1.DrawClone("AP");

    gROOT->SetSelectedPad(0);
    gPad->GetCanvas()->cd(4);

    gPad->SetGridx();
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);

    g2.SetMinimum(0);
    g2.SetMaximum(1);
    g2.SetMarkerStyle(kFullDotMedium);
    g2.GetXaxis()->SetTimeDisplay(true);
    g2.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    g2.GetXaxis()->SetLabelSize(0.12);
    g2.GetYaxis()->SetLabelSize(0.1);
    g2.GetYaxis()->SetTitle("RELATIVE ON TIME");
    g2.GetYaxis()->SetTitleOffset(0.2);
    g2.GetYaxis()->SetTitleSize(0.1);
    g2.DrawClone("AP");


    return 0;
}

Int_t PlotPointing(TArrayD **vec, TString fname)
{
    fname += ".DRIVE_CONTROL_POINTING_POSITION.fits";

    fits file(fname.Data());
    if (!file)
    {
        cerr << fname << ": " << gSystem->GetError() << endl;
        return -2;
    }

    //cout << fname << endl;

    Double_t time;
    Double_t zd;

    if (!file.SetPtrAddress("Time",  &time))
        return -1;
    if (!file.SetPtrAddress("Zd", &zd))
        return -1;

    TGraph g;

    while (file.GetNextRow())
        if (Contains(vec, time))
            g.SetPoint(g.GetN(), time*24*3600, 90-zd);

    g.SetMinimum(0);
    g.SetMaximum(90);
    g.SetMarkerStyle(kFullDotMedium);
    g.GetXaxis()->SetTimeDisplay(true);
    g.GetXaxis()->SetTimeFormat("%H:%M %F1995-01-01 00:00:00 GMT");
    g.GetXaxis()->SetLabelSize(0.1);
    g.GetYaxis()->SetLabelSize(0.1);
    g.GetYaxis()->SetTitle("ELEVATION");
    g.GetYaxis()->SetTitleOffset(0.2);
    g.GetYaxis()->SetTitleSize(0.1);
    g.DrawClone("AP");

    return 0;
}
/*
Int_t PlotMem(TArrayD **vec, TString fname)
{
    fname += ".FAD_CONTROL_STATISTICS1.fits";

    fits file(fname.Data());
    if (!file)
    {
        cerr << fname << ": " << gSystem->GetError() << endl;
        return -2;
    }

    //cout << fname << endl;

    Double_t time;
    UInt_t buf[5];
    ULong64_t mem[4];
    UInt_t rate[2];

    if (!file.SetPtrAddress("Time",  &time))
        return -1;
    if (!file.SetPtrAddress("bufferInfo", buf) ||
        !file.SetPtrAddress("memInfo",    mem) ||
        !file.SetPtrAddress("rateNew",    rate))
        return -1;

    TGraph g;

    while (file.GetNextRow())
        if (Contains(vec, time))
            g.SetPoint(g.GetN(), time*24*3600, mem[1]);

    g.SetMarkerStyle(kFullDotMedium);
    g.GetXaxis()->SetTimeDisplay(true);
    g.GetXaxis()->SetTimeFormat("%H:%M");
    g.GetXaxis()->SetTimeOffset(0, "gmt");
    g.GetXaxis()->SetLabelSize(0.12);
    g.GetYaxis()->SetLabelSize(0.1);
    //g.GetYaxis()->SetTitle("ELEVATION");
    g.GetYaxis()->SetTitleOffset(0.2);
    g.GetYaxis()->SetTitleSize(0.1);
    g.DrawClone("AP");

    return 0;
}
*/
void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0)
{
    // To get correct dates in the histogram you have to add
    // the MJDREF offset (should be 40587) and 9131.

    if (y==0)
    {
        UInt_t nt = MTime(MTime(-1).GetMjd()-1.5).GetNightAsInt();
        y = nt/10000;
        m = (nt/100)%100;
        d = nt%100;

        cout << y << "/" << m << "/" << d << endl;
    }

    TString fname=Form("/fact/aux/%04d/%02d/%02d/%04d%02d%02d", y, m, d, y, m, d);

    UInt_t night = MTime(y, m, d, 0).GetNightAsInt();

    MSQLMagic serv("sql.rc");
    Bool_t con = serv.IsConnected();

    cout << "quality" << endl;
    cout << "-------" << endl;
    cout << endl;
    if (con)
    {
        cout << "Connected to " << serv.GetName() << endl;
        cout << endl;
    }
    cout << "Night: " << Form("%04d-%02d-%02d", y, m, d) << endl;
    cout << endl;

    TArrayD run, beg, end;

    TArrayD *runs[3] = { &run, &beg, &end };

    if (con)
    {
        TString query;
        query += "SELECT fRunID, fRunStart, fRunStop FROM runinfo";
        query += " WHERE fNight=";
        query += night;
        query += " AND fRunTypeKEY=1 ORDER BY fRunStart";

        TSQLResult *res = serv.Query(query);
        if (!res)
            return kFALSE;

        run.Set(res->GetRowCount());
        beg.Set(res->GetRowCount());
        end.Set(res->GetRowCount());

        Int_t n = 0;

        TSQLRow *row = 0;
        while ((row=res->Next()))
        {
            run[n] = atoi((*row)[0]);
            beg[n] = MTime((*row)[1]).GetMjd()-40587;
            end[n] = MTime((*row)[2]).GetMjd()-40587;
            n++;
            delete row;
        }

        delete res;

        if (n==0)
            cout << "WARNING - No data runs in db, displaying all data." << endl;
    }

    TCanvas *c = new TCanvas("quality", Form("Quality %04d/%02d/%02d", y, m, d), 1280, 960);
    c->Divide(1, 5, 1e-5, 1e-5);

    gROOT->SetSelectedPad(0);
    c->cd(1);
    gPad->SetGridx();
    gPad->SetTopMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    gPad->SetBottomMargin(0);
    cout << PlotThresholds(runs, fname) << endl;

    gROOT->SetSelectedPad(0);
    c->cd(2);
    gPad->SetGridx();
    gPad->SetTopMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    gPad->SetBottomMargin(0);
    cout << PlotCurrent(runs, fname) << endl;

    gROOT->SetSelectedPad(0);
    c->cd(3);
    gPad->SetGridx();
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    cout << PlotRate(runs, fname) << endl;

    gROOT->SetSelectedPad(0);
    c->cd(5);
    gPad->SetGridx();
    gPad->SetTopMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    cout << PlotPointing(runs, fname) << endl;
/*
    gROOT->SetSelectedPad(0);
    c->cd(6);
    gPad->SetGridx();
    gPad->SetTopMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    cout << PlotMem(runs, fname) << endl;
*/
    c->SaveAs(Form("%04d%02d%02d-quality.pdf", y, m, d));
}

//20130314_141