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;
    g.SetName("Threshold");

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

    g.SetMinimum(281);
    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, Nova::EquPosn>> vecp;

Nova::EquPosn FindPointing(Double_t time)
{
    for (int i=0; i<vecp.size(); i++)
        if (time<vecp[i].first)
        {
            if (i==0)
                return Nova::EquPosn();
            else
                return vecp[i-1].second;
        }

    return vecp[vecp.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 = Nova::GetAngularSeparation(moon, hrzm)*TMath::DegToRad();

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

    double disk = Nova::GetLunarDisk(jd);

    //                                                                           semi-major axis
    double lc = calt*pow(disk, 2.2)*pow(Nova::GetLunarEarthDist(jd)/384400, -2);///sqrt(cang);

    return lc>0 ? 4+103*lc : 4;
}*/

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 = Nova::GetAngularSeparation(moon, pos);

    // Distance between earth and moon relative to major semi-axis
    const double dist  = Nova::GetLunarEarthDist(jd)/384400;

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

    double disk = Nova::GetLunarDisk(jd);

    // light condition
    double lc = sqrt(calt)*pow(disk, 2.3)*pow(cang+1, 1)*pow(dist, -2);

    return lc>0 ? 7.2 + 69*lc : 7.2;
    */

    // Sun properties
    Nova::EquPosn  sun  = Nova::GetSolarEquCoords(jd);
    Nova::ZdAzPosn hrzs = Nova::GetHrzFromEqu(sun, jd);

    // Get source position
    Nova::EquPosn pos = FindPointing(time);

    // Moon properties
    Nova::EquPosn moon = Nova::GetLunarEquCoords(jd, 0.01);
    Nova::HrzPosn hrzm = Nova::GetHrzFromEqu(moon, jd);
    double        disk = Nova::GetLunarDisk(jd);

    // Derived moon properties
    double angle = Nova::GetAngularSeparation(moon, pos);
    double edist = Nova::GetLunarEarthDist(jd)/384400;

    // Current prediction
    double sin_malt  = hrzm.alt<0 ? 0 : sin(hrzm.alt*TMath::DegToRad());
    double sin_mdist = sin(angle*TMath::DegToRad());
    double cos_salt  = cos(hrzs.zd*TMath::DegToRad());

    double c0 = pow(disk,     2.52);
    double c1 = pow(sin_malt, 0.72);
    double c2 = pow(edist,   -2.00);
    double c3 = exp(1.46*(1-sin_mdist)*(1-sin_mdist));
    double c4 = exp(33.0)*exp(-20.1*(1-cos_salt)*(1-cos_salt));

    double cur = 6.4 + 96.9*c0*c1*c2*c3 + c4;

   // cout << cur << " " << hrzm.alt << " " << c1 << endl;

    return cur;

}

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())
    {
        Nova::EquPosn p;
        p.ra  = ra*15;
        p.dec = dec;

        vecp.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;
    g1.SetName("Currents");
    g2.SetName("Prediction");

    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;
    g1.SetName("TriggerRate");
    g2.SetName("RelOnTime");

    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(269);
    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->SetGrid();
    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;
}

void PlotRateQC(UInt_t night, MSQLServer &serv)
{
    TString query = 
        "LEFT JOIN AnalysisResultsRunLP USING(fNight, fRunID) "
        "WHERE fRunTypeKey=1 AND NOT ISNULL (AnalysisResultsRunLP.fNumEvtsAfterCleaning) AND fNight=";
    query += night;

    TTree *t = serv.GetTree("RunInfo", query);
    if (!t)
        return;

    int save = gErrorIgnoreLevel;
    gErrorIgnoreLevel = kFatal;

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

    t->Draw("AnalysisResultsRunLP.fNumEvtsAfterCleaning/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same");
    TGraph *g = (TGraph*)gPad->GetPrimitive("Graph");
    if (g)
    {
        g->SetName("CleaningRate");
        g->SetMarkerColor(kRed);
        g->SetMarkerStyle(29);//kFullDotMedium);
    }

    t->Draw("AnalysisResultsRunLP.fNumEvtsAfterQualCuts/AnalysisResultsRunLP.fOnTimeAfterCuts:(RunInfo.fRunStart+RunInfo.fRunStop)/2+9131*24*3600", "", "same");
    g = (TGraph*)gPad->GetPrimitive("Graph");
    if (g)
    {
        g->SetName("RateAfterQC");
        g->SetMarkerColor(kBlue);
        g->SetMarkerStyle(29);//kFullDotMedium);
    }

    gErrorIgnoreLevel = save;
}


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;
    g.SetName("Zd");

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

    g.SetMinimum(1);
    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.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;
}

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

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

    //cout << fname << endl;

    Double_t time;
    Float_t  temp;

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

    TGraph g;
    g.SetName("ContainerTemp");

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

    g.SetMinimum(-5);
    g.SetMaximum(49);
    g.SetMarkerStyle(kFullDotMedium);
    g.SetMarkerColor(kRed);
    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("TEMP");
    g.GetYaxis()->SetTitleOffset(0.2);
    g.GetYaxis()->SetTitleSize(0.1);
    g.DrawClone("AP");

    return 0;
}

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

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

    //cout << fname << endl;

    Double_t time;
    Float_t  temp;

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

    TGraph g;
    g.SetName("OutsideTemp");

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

    g.SetMarkerStyle(kFullDotMedium);
    g.SetMarkerColor(kBlue);
    g.DrawClone("P");

    return 0;
}

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

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

    //cout << fname << endl;

    Double_t time;
    Float_t  temp[31];

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

    TGraph g, g1, g2;
    g.SetName("SensorTempAvg");
    g1.SetName("SensorTempMin");
    g2.SetName("SensorTempMax");

    while (file.GetNextRow())
        if (Contains(vec, time))
        {
            float min =  100;
            float max = -100;
            double avg = 0;
            int num = 0;
            for (int i=0; i<31; i++)
                if (temp[i]!=0)
                {
                    avg += temp[i];
                    num++;

                    min = TMath::Min(min, temp[i]);
                    max = TMath::Max(max, temp[i]);
                }

            g.SetPoint(g.GetN(), time*24*3600, avg/num);
            g1.SetPoint(g1.GetN(), time*24*3600, min);
            g2.SetPoint(g2.GetN(), time*24*3600, max);
        }

    g.SetMarkerStyle(kFullDotMedium);
    g.DrawClone("P");

    /*
     g1.SetLineWidth(1);
     g1.DrawClone("L");

     g2.SetLineWidth(1);
     g2.DrawClone("L");
    */
    return 0;
}

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

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

    //cout << fname << endl;

    Double_t time;
    Float_t  temp[160];

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

    Int_t num = file.GetN("temp")==0 ? file.GetN("Data1") : file.GetN("temp");
    Int_t beg = num==82 ? 2 : 0;

    TGraphErrors g1;
    TGraph g2,g3;

    g1.SetName("FadTempAvg");
    g2.SetName("FadTempMin");
    g3.SetName("FadTempMax");

    while (file.GetNextRow())
        if (Contains(vec, time))
        {
            double avg = 0;
            double rms = 0;
            float min =  100;
            float max = -100;
            for (int i=beg; i<num; i++)
            {
                avg += temp[i];
                rms += temp[i]*temp[i];

                min = TMath::Min(min, temp[i]);
                max = TMath::Max(max, temp[i]);
            }

            avg /= num-beg;
            rms /= num-beg;

            g1.SetPoint(g1.GetN(), time*24*3600, avg);
            g1.SetPointError(g1.GetN()-1, 0, sqrt(rms-avg*avg));

            g2.SetPoint(g2.GetN(), time*24*3600, min);
            g3.SetPoint(g3.GetN(), time*24*3600, max);
        }

    g1.SetLineColor(kGreen);
    g1.DrawClone("[]");

    g2.SetLineColor(kGreen);
    g2.SetLineWidth(1);
    g2.DrawClone("L");

    g3.SetLineColor(kGreen);
    g3.SetLineWidth(1);
    g3.DrawClone("L");

    return 0;
}

void quality(UInt_t y=0, UInt_t m=0, UInt_t d=0, TString outpath="./")
{
    // 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;
        else
            cout << "Num: " << n << "\n" << endl;
    }

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

    gROOT->SetSelectedPad(0);

    c->cd(1);
    gPad->SetGrid();
    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->SetGrid();
    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->SetGrid();
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    cout << PlotRate(runs, fname) << endl;
    cout << PlotRateQC(night, serv) << endl;

    gROOT->SetSelectedPad(0);
    c->cd(5);
    gPad->SetGrid();
    gPad->SetTopMargin(0);
    gPad->SetBottomMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    cout << PlotPointing(runs, fname) << endl;

    gROOT->SetSelectedPad(0);
    c->cd(6);
    gPad->SetGrid();
    gPad->SetTopMargin(0);
    gPad->SetRightMargin(0.001);
    gPad->SetLeftMargin(0.04);
    cout << PlotTemperature1(runs, fname) << endl;
    cout << PlotTemperature2(runs, fname) << endl;
    cout << PlotTemperature3(runs, fname) << endl;
    cout << PlotTemperature4(runs, fname) << endl;

    c->SaveAs(Form("%s/%04d%02d%02d-quality.png", outpath.Data(), y, m, d));
}

//20130314_141