#include <iomanip>

Double_t GetResidual(Double_t fRawEl,  Double_t fRawAz,
                     Double_t fStarEl, Double_t fStarAz)
{
    fRawEl *= TMath::DegToRad();
    fRawAz *= TMath::DegToRad();
    fStarEl *= TMath::DegToRad();
    fStarAz *= TMath::DegToRad();
    Double_t del = fRawEl-fStarEl;
    Double_t daz = fRawAz-fStarAz;
    Double_t dphi2 = daz/2.;
    Double_t cos2  = cos(dphi2)*cos(dphi2);
    Double_t sin2  = sin(dphi2)*sin(dphi2);
    Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;

    Double_t dist = acos(d);

    return dist * TMath::RadToDeg();
}

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 from " << 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",
               &alt, &az, &dalt, &daz, &mjd);

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

        //cout << dalt << " " << daz << " " << GetResidual(alt, az, alt+dalt, az+daz) << endl;

        mjd -=  53140.097505;

        Double_t res = GetResidual(alt, az, alt-dalt, az-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 (" << g[0].GetN() << ")" << endl;

    return g[0].GetN();
}

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

const Int_t counts = 13-4;
Description_t desc[counts] =
{
    // No good pointing model has been applied yet
  //  {"0404",  "TPoints Residuals 4/2004" ,   "tpoint/tpoint0404.txt"},
  //  {"0405",  "TPoints Residuals 5/2004" ,   "tpoint/tpoint0405.txt"},
  //  {"04081", "TPoints Residuals 8/2004-1" , "tpoint/tpoint0408-1.txt"},
    {"04082", "TPoints Residuals 8/2004-2" , "tpoint/tpoint0408-2.txt"},
    {"0409",  "TPoints Residuals 9/2004" ,   "tpoint/tpoint0409.txt"},
    // Culmination tests
  //  {"0410",  "TPoints Residuals 10/2004" ,  "tpoint/tpoint0410.txt"},
    {"0411",  "TPoints Residuals 11/2004" ,  "tpoint/tpoint0411.txt"},
    {"0412",  "TPoints Residuals 12/2004" ,  "tpoint/tpoint0412.txt"},
    // Worse pointing due to realignment of the mirror
    {"0503",  "TPoints Residuals 3/2005" ,   "tpoint/tpoint0503.txt"},
    {"0504",  "TPoints Residuals 4/2005" ,   "tpoint/tpoint0504.txt"},
    {"05051", "TPoints Residuals 5/2004-1" , "tpoint/tpoint0505-1.txt"},
    // Mirror alignment has been fixed
    {"05052", "TPoints Residuals 5/2004-2" , "tpoint/tpoint0505-2.txt"},
    // Fixes to pointing model due to fixing a screw
    {"0506", "TPoints Residuals 6/2004" , "tpoint/tpoint0506.txt"}
};

void plot()
{
    TGraph g[3];

    MBinning bins(10, 0, 0.2);

    TArrayI n(counts);

    TH1F hx[counts];
    for (int i=0; i<counts; i++)
    {
        hx[i].SetNameTitle(desc[i].fName, desc[i].fTitle);
        hx[i].SetDirectory(0);
        bins.Apply(hx[i]);

        cout << setw(2) << i << ": " << flush;
        n[i] = fill(desc[i].fFile, g, &hx[i]);
    }

    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");

    return;

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