#include <iomanip>

void DrawMarker(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1)
    {
        TView *view = pad->GetView();

        if (!view)
        {
            cout << "No View!" << endl;
            return;
        }

        TMarker mark0;
        TMarker mark1;
        mark0.SetMarkerStyle(kStar);
        mark1.SetMarkerStyle(kStar);
        mark1.SetMarkerColor(kRed);
    
        r0 /= 90;
        r1 /= 90;
        phi0 *= TMath::DegToRad();
        phi1 *= TMath::DegToRad();
    
        Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0};
        Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0};

        mark0.DrawMarker(x0[0], x0[1]);
        mark1.DrawMarker(x1[0], x1[1]);

        return;
        Double_t y0[3], y1[3];
    
        view->WCtoNDC(x0, y0);
        view->WCtoNDC(x1, y1);
    
        mark0.DrawMarker(y0[0], y0[1]);
        mark1.DrawMarker(y1[0], y1[1]);
    }

int fill(const char *fname, TGraph *g, TH1 *h)
{
/*
    TH2F h2res1("Res2D1", " Dataset positions on the sky ", 32, 0, 360,  10, 0, 90);
    h2res1.SetBit(TH1::kNoStats);
    h2res1.DrawCopy("surf1pol");
    gPad->Modified();
    gPad->Update();
    gPad->SetTheta(90);
    gPad->SetPhi(-90);

    DrawMarker(gPad, 45, 0, 0, 0);
    gPad->Modified();
    gPad->Update();

    return;

  */
    ifstream fin(fname);

    cout << "Reading " << setw(23) << fname << "..." << flush;

    while (1)
    {
        TString str;
        str.ReadLine(fin);
        if (!fin)
            break;

        if (str.Contains("#"))
            continue;

        Float_t alt, az, dalt, daz, mjd;
        sscanf(str.Data(), "%f %f %*f %*f %*f %*f %f %f %f",
               &az, &alt, &dalt, &daz, &mjd);

        if (dalt==0/* ||  GetResidual(alt, az, alt+dalt, az+daz)>0.1*/)
            continue;

        mjd -=  53140.097505;

        Double_t res = MAstro::GetDevAbs(90-(alt-dalt), -dalt, -daz);

        g[0].SetPoint(g[0].GetN(), g[0].GetN(), fabs(dalt));
        g[1].SetPoint(g[1].GetN(), g[1].GetN(), fabs(daz));
        g[2].SetPoint(g[2].GetN(), g[2].GetN(), res);

        h->Fill(res);
    }

    cout << "done (" << setw(3) << (Int_t)h->GetEntries() << "/";
    cout << setw(4) << g[0].GetN() << ") " << flush;

    return g[0].GetN();
}

struct Description_t
{
    const char *fName;
    const char *fTitle;
    const char *fFile;
};

const Int_t counts = 27;
Description_t desc[counts] =
{
    {"090401",  "TPoints Residuals 8/2004-2" , "tpoint/m2/first/tpoints_m2_1.txt"},
    {"+090401",  "TPoints Residuals 9/2004" ,  "tpoint/m2/first/tpoints_m2_2.txt"},
    {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090403_025841.txt"},
    {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090403_214619.txt"},
    {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090404_235955.txt"},
    {"+090401",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090405_015920.txt"},
    // New pointing model
    {"090402",   "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090411_231731.txt"},
    {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090412_232919.txt"},
    {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090414_001746.txt"},
    {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090418_012742.txt"},
    {"+090402",  "TPoints Residuals 11/2004" , "tpoint/m2/first/tpoint_20090418_024412.txt"},
    // New pointing model
    {"090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_12/tpoint_20090512_011933.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_12/tpoint_20090512_013112.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090512_210644.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_002758.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_025123.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_13/tpoint_20090513_033504.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_05_14/tpoint_20090514_013332.txt"},
    // Just start of a new observation period
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_01/tpoint_20090531_215139.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_01/tpoint_20090531_222548.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_02/tpoint_20090601_223001.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090602_213508.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090602_230944.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090603_000616.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_03/tpoint_20090603_011935.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_04/tpoint_20090603_215841.txt"},
    {"+090512",  "TPoints Residuals 11/2004" , "tpoint/m2/2009_06_05/tpoint_20090604_215941.txt"},
};

void plot_m2()
{
    TGraph g[3];

    MBinning bins(100, 0, 0.2);

    TArrayI n(counts);

    TH1F hx[counts];
    Int_t num = -1;
    for (int i=0; i<counts; i++)
    {
        if (desc[i].fName[0]!='+')
        {
            num++;

            hx[num].SetNameTitle(desc[i].fName, desc[i].fTitle);
            hx[num].SetDirectory(0);
            bins.Apply(hx[num]);
        }

        cout << setw(2) << num << ": " << flush;
        n[num] = fill(desc[i].fFile, g, &hx[num]);
        cout << " Mean: " << setw(5) << setprecision(2) << hx[num].GetMean() << " deg +/- " <<  hx[num].GetRMS()<< endl;
    }

    n.Set(++num);

    g[0].SetMarkerColor(kGreen);
    g[1].SetMarkerColor(kMagenta);
    g[2].SetMarkerColor(kBlack);
    //g[2].SetLineColor(kBlack);
    g[0].SetMarkerStyle(kFullDotMedium);
    g[1].SetMarkerStyle(kFullDotMedium);
    g[2].SetMarkerStyle(kFullDotLarge);

    // --------- First Canvas ----------

    new TCanvas("Vs Time", "");

    TObject *obj[4];

    for (int i=0; i<3; i++)
    {
        g[i].SetFillColor(kWhite);
        g[i].SetLineColor(kWhite);
        obj[i] = g[i].Clone();
        obj[i]->SetBit(kCanDelete);
        obj[i]->Draw(i==0?"AP":"P");
    }

    TLegend leg(0.905, 0.86, 0.99, 0.99);
    leg.AddEntry(obj[0], "  \\Delta\\theta");
    leg.AddEntry(obj[1], "  \\Delta\\phi");
    leg.AddEntry(obj[2], "  \\Delta");
    leg.DrawClone()->SetBit(kCanDelete);

    TLine l;
    l.SetLineColor(kBlue);
    for (int i=0; i<n.GetSize(); i++)
        l.DrawLine(n[i], 0, n[i], 0.2);

    // --------- Second Canvas ----------

    new TCanvas("Distrib", "");

    Double_t max=0;
    for (int i=0; i<n.GetSize(); i++)
    {
        if (hx[i].GetEntries()==0)
        {
            cout << "Skip #" << i << endl;
            continue;
        }

        hx[i].Scale(1./hx[i].GetEntries());
        max = TMath::Max(max, hx[i].GetMaximum());
    }
    for (int i=0; i<n.GetSize(); i++)
    {
        hx[i].SetMaximum(max*1.05);
        if (i<6)
            hx[i].SetLineColor(kRed+i);
        else
            hx[i].SetMarkerStyle(kPlus+i-6);


    }

    for (int i=0; i<counts; i++)
        hx[i].DrawCopy(i==0?"LP":"LPsame");

    Double_t time[] = {
        MTime(2009,  3, 30).GetAxisTime(),
        MTime(2009,  4, 11).GetAxisTime(),
        MTime(2009,  5, 11).GetAxisTime(),
        MTime(2009,  5, 31).GetAxisTime(),
        MTime(2009,  6,  6).GetAxisTime(),
    };

    TH1D histres[4];

    MBinning bins;
    bins.SetEdges(TArrayD(5, time));

    bins.Apply(histres[0]);
    bins.Apply(histres[1]);
    bins.Apply(histres[2]);
    bins.Apply(histres[3]);

    TGraphAsymmErrors result[4];
    TGraph resultm;
    for (int i=0; i<n.GetSize(); i++)
    {
        cout << i+1 << " - Mean:  " << Form("%.4f +- %.4f", hx[i].GetMean(), hx[i].GetRMS());
        cout << "   (Overflows=" << hx[i].GetBinContent(hx[i].GetNbinsX()+1)*hx[i].GetEntries() << ")   " << (int)hx[i].GetEntries() << endl;

        /*
         TF1 fg("fg", "gaus", 0, 1);
         hx[i].Fit(&fg);
         result2.SetPoint(result.GetN(), result.GetN()+7, hx[i].GetMean());
         result2.SetPointError(result.GetN()-1, 0, hx[i].GetRMS()/sqrt(hx[i].GetEntries()));
         result.SetPoint(result.GetN(), result.GetN()+7, fg.GetParameter(1));
         result.SetPointError(result.GetN()-1, 0, fg.GetParameter(2));
         */

        //Double_t q[4] = { MMath::GaussProb(0.5), MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
        Double_t q[4] = { 0.5, MMath::GaussProb(1), MMath::GaussProb(2), MMath::GaussProb(3) };
        Double_t rc[4];
        hx[i].GetQuantiles(4, rc, q);

        for (int j=0; j<4; j++)
        {
            histres[j].SetBinContent(i+1, 60*rc[j]);

            result[j].SetPoint(i, time[i], 60*rc[j]);
            result[j].SetPointError(i, 0, 0, 60*rc[j], 0);

            //result[j].SetPoint(result[j].GetN(), result[j].GetN()+1, 60*rc[j]);
            //result[j].SetPointError(result[j].GetN()-1, 0, 0, 60*rc[j], 0);
        }

//        result.SetPoint(result.GetN(), result.GetN()+1, rc[1]);
//        result2.SetPoint(result2.GetN(), result2.GetN()+1, rc[2]);

        //result.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);
        //result2.SetPointError(result.GetN()-1, 0.5, 0);//rc[2]-rc[1]);

//        result.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);
//        result2.SetPointError(result.GetN()-1, 0, 0, 0/*rc[1]-rc[0]*/, rc[2]-rc[1]);

//        result2.SetPointError(result.GetN()-1, 0, 0, rc[2], 0);

//        resultm.SetPoint(resultm.GetN(), resultm.GetN()+1, 60*hx[i].GetMean());
        resultm.SetPoint(resultm.GetN(), (time[i]+time[i+1])/2, 60*hx[i].GetMean());

    }

    new TCanvas;

    gPad->SetBorderMode(0);
    gPad->SetFrameBorderMode(0);
    gPad->SetFillColor(kWhite);
    gPad->SetRightMargin(0.01);
    gPad->SetTopMargin(0.01);
    gPad->SetLeftMargin(0.06);
    gPad->SetGridy();

    //Int_t col[] = { 12, 15, 17, 19 };
    //Int_t col[] = { 12, 16, 18, 0 };
    Int_t col[] = { 13, 16, 19, 0 };

    TH1 *h = &histres[2];//result[3].GetHistogram();
    h->SetXTitle("");
    h->SetYTitle("Residual / arcmin");
    h->SetBit(TH1::kNoStats);
    h->GetXaxis()->CenterTitle();
    h->GetYaxis()->CenterTitle();
    h->GetYaxis()->SetTitleOffset(0.75);
//    h->GetXaxis()->SetTimeFormat("%m/%y %F1995-01-01 00:00:00 GMT");
//    h->GetXaxis()->SetTimeDisplay(1);

    h->GetXaxis()->SetLabelColor(kWhite);

    TLine line;

    for (int j=2; j>=0; j--)
    {
        histres[j].SetMinimum(0);
        histres[j].SetFillColor(col[j]);//12+2*j);

        histres[j].DrawCopy(j==2?"":"same");

        /*
        //result[j].SetLineColor(kBlue);
        //result[j].SetLineWidth(2);
        //result[j].SetMarkerColor(kBlue);
        result[j].SetMinimum(0);
        //result2.SetMarkerStyle(kFullDotMedium);
        //result[j].SetMarkerStyle(23);
        result[j].SetFillColor(col[j]);//12+2*j);
        result[j].DrawClone(j==3 ? "ABX" : "B"); // E3 B
        */
    }


    resultm.SetMarkerStyle(20);
    resultm.SetMarkerSize(0.8);
    resultm.DrawClone("P");

    line.DrawLine(time[0], 0, time[0], histres[2].GetMaximum());
    for (int i=0; i<bins.GetNumBins(); i++)
        line.DrawLine(time[i+1], 0, time[i+1], histres[2].GetBinContent(i+1));

    TText txt;
    txt.SetTextSize(0.037);
//    txt.SetTextAngle(-45);
    txt.SetTextAlign(23);
    /*
    for (int m=4; m<13; m++)
    {
        Double_t monl = MTime(2005, m,   1,0).GetAxisTime();
        Double_t mont = MTime(2005, m,   15,0).GetAxisTime();
//        txt.DrawText(mon, -0.12, Form("%02d/05", m));
        txt.DrawText(mont, -0.10, Form("%d", m));
        line.DrawLine(monl, -0.12, monl, 0.12);
    }
    for (int m=1; m<13; m++)
    {
        Double_t monl = MTime(2006, m,   1,0).GetAxisTime();
        Double_t mont = MTime(2006, m,   15,0).GetAxisTime();
//        txt.DrawText(mon, -0.12, Form("%02d/06", m));
        txt.DrawText(mont, -0.10, Form("%d", m));
        line.DrawLine(monl, -0.12, monl, 0.12);
    }
    for (int m=1; m<13; m++)
    {
        Double_t monl = MTime(2007, m,   1,0).GetAxisTime();
        Double_t mont = MTime(2007, m,   15,0).GetAxisTime();
//        txt.DrawText(mon, -0.12, Form("%02d/07", m));
        txt.DrawText(mont, -0.10, Form("%d", m));
        line.DrawLine(monl, -0.12, monl, 0.12);
    }*/
    for (int m=1; m<13; m++)
    {
        Double_t monl = MTime(2009, m,   1,0).GetAxisTime();
        Double_t mont = MTime(2009, m,   15,0).GetAxisTime();
//        txt.DrawText(mon, -0.12, Form("%02d/07", m));
        txt.DrawText(mont, -0.10, Form("%d", m));
        line.DrawLine(monl, -0.12, monl, 0.12);
    }

//    Double_t y6 = MTime(2006,1,1,0).GetAxisTime();
//    Double_t y7 = MTime(2007,1,1,0).GetAxisTime();
//    Double_t y8 = MTime(2008,1,1,0).GetAxisTime();
    Double_t y9 = MTime(2009,1,1,0).GetAxisTime();
    Double_t y0 = MTime(2010,1,1,0).GetAxisTime();

    txt.SetTextSize(0.042);
//    txt.DrawText(y6-(y7-y6)/2, -0.6, "2005");
//    txt.DrawText((y6+y7)/2,    -0.6, "2006");
//    txt.DrawText((y7+y8)/2,    -0.6, "2007");
    txt.DrawText((y9+y0)/2,    -0.6, "2009");

//    line.DrawLine(y6, -0.7, y6, 0.26);
//    line.DrawLine(y7, -0.7, y7, 0.26);
//    line.DrawLine(y8, -0.7, y8, 0.26);
    line.DrawLine(y9, -0.7, y9, 0.26);
    line.SetLineStyle(3);
//    line.DrawLine(y6, 0.26, y6, 1.05*histres[2].GetMaximum());
//    line.DrawLine(y7, 0.26, y7, 1.05*histres[2].GetMaximum());
//    line.DrawLine(y8, 0.26, y8, 1.05*histres[2].GetMaximum());
    line.DrawLine(y9, 0.26, y9, 1.05*histres[2].GetMaximum());

    line.SetLineColor(kBlue);
    line.SetLineWidth(2);
    line.SetLineStyle(kSolid);
    line.DrawLine(time[0], 1*360/16384.*60, time[bins.GetNumBins()], 1*360/16384.*60);
    line.SetLineStyle(9);
    line.DrawLine(time[0], 2*360/16384.*60, time[bins.GetNumBins()], 2*360/16384.*60);
    line.SetLineStyle(7);
    line.DrawLine(time[0], 3*360/16384.*60, time[bins.GetNumBins()], 3*360/16384.*60);

    /*    result.SetMinimum(-0.06);
    result.SetMarkerStyle(kFullDotMedium);
    result.DrawClone("LP");*/
}