#include <fstream.h>
#include <fstream.h>
#include <fstream.h>
#include <iostream.h>
#include <iomanip.h>

#include <TError.h>

#include <TGFrame.h>
#include <TGLabel.h>
#include <TGButton.h>
#include <TGFileDialog.h>

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

#include <TList.h>
#include <TStyle.h>
#include <TMinuit.h>

#include <TView.h>
#include <TLine.h>
#include <TMarker.h>
#include <TCanvas.h>

#include "coord.h"

#include "MGList.h"
#include "MBending.h"

class Set : public TObject
{
    friend ifstream &operator>>(ifstream &fin, Set &set);
private:
    Double_t fStarAz;
    Double_t fStarEl;

    Double_t fRawAz;
    Double_t fRawEl;

    Double_t fMag;
public:
    Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) :
        fStarAz(saz*TMath::DegToRad()),
        fStarEl(sel*TMath::DegToRad()),
        fRawAz(raz*TMath::DegToRad()),
        fRawEl(rel*TMath::DegToRad()), fMag(-25)
    {
    }

    Double_t GetMag() const { return fMag; }
    Double_t GetResidual(Double_t *err=0) const
    {
      /*
       TVector3 v1, v2;
       v1.SetMagThetaPhi(1, TMath::Pi()/2-fRawEl, fRawAz);
       v2.SetMagThetaPhi(1, TMath::Pi()/2-fStarEl, fStarAz);

       return v1.Angle(v2)*TMath::RadToDeg();
       */

        const Double_t del = fRawEl-fStarEl;
        const Double_t daz = fRawAz-fStarAz;
        /*
        const Double_t dphi2 = daz/2.;
        const Double_t cos2  = cos(dphi2)*cos(dphi2);
        const Double_t sin2  = sin(dphi2)*sin(dphi2);
        const Double_t d = cos(del)*cos2 - cos(fRawEl+fStarEl)*sin2;
        */

        const Double_t d  = cos(del) - cos(fRawEl)*cos(fStarEl)*(1.-cos(daz));

        if (err)
        {
            // Error of one pixel in the CCD
            const Double_t e1 = 32./3600*TMath::DegToRad()   * 0.5;

            // Error of one SE unit
            const Double_t e2 = 360./16384*TMath::DegToRad() * 0.5;

            const Double_t e11 =  sin(del)+cos(fRawEl)*sin(fStarEl)*(1-cos(daz));
            const Double_t e12 =  cos(fRawEl)*cos(fStarEl)*sin(daz);

            const Double_t e21 = -sin(del)+sin(fRawEl)*cos(fStarEl)*(1-cos(daz));
            const Double_t e22 = -cos(fRawEl)*cos(fStarEl)*sin(daz);

            const Double_t err1  = sqrt(1-d*d);
            const Double_t err2  = (e11*e11 + e12*e12)*e1*e1;
            const Double_t err3  = (e21*e21 + e22*e22)*e2*e2;

            *err = sqrt(err2+err3)/err1 * TMath::RadToDeg();
        }

        const Double_t dist = acos(d);
        return dist * TMath::RadToDeg();
    }

    void operator=(Set &set)
    {
        fStarAz = set.fStarAz;
        fStarEl = set.fStarEl;
        fRawAz  = set.fRawAz;
        fRawEl  = set.fRawEl;
        fMag    = set.fMag;
    }

    Double_t GetDEl() const     { return (fRawEl-fStarEl)*TMath::RadToDeg(); }
    Double_t GetDZd() const     { return -GetDEl(); }
    Double_t GetDAz() const     { return (fRawAz-fStarAz)*TMath::RadToDeg(); }
    Double_t GetStarEl() const  { return fStarEl*TMath::RadToDeg(); }
    Double_t GetStarZd() const  { return 90.-fStarEl*TMath::RadToDeg(); }
    Double_t GetStarAz() const  { return fStarAz*TMath::RadToDeg(); }
    Double_t GetRawEl() const   { return fRawEl*TMath::RadToDeg(); }
    Double_t GetRawAz() const   { return fRawAz*TMath::RadToDeg(); }
    Double_t GetRawZd() const   { return 90.-fRawEl*TMath::RadToDeg(); }

    ZdAz  GetStarZdAz() const   { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); }
    AltAz GetStarAltAz() const  { return AltAz(fStarEl, fStarAz); }

    ZdAz  GetRawZdAz() const    { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); }
    AltAz GetRawAltAz() const   { return AltAz(fRawEl, fRawAz); }

    void AdjustEl(Double_t del) { fStarEl += del*TMath::DegToRad(); }
    void AdjustAz(Double_t daz) { fStarAz += daz*TMath::DegToRad(); }

    void Adjust(const MBending &bend)
    {
        AltAz p = bend(GetStarAltAz());
        fStarEl = p.Alt();
        fStarAz = p.Az();
    }
    void AdjustBack(const MBending &bend)
    {
        AltAz p = bend.CorrectBack(GetRawAltAz());
        fRawEl = p.Alt();
        fRawAz = p.Az();
    }
    ClassDef(Set, 0)
};

ClassImp(Set);

ifstream &operator>>(ifstream &fin, Set &set)
{
    TString str;
    do
    {
        str.ReadLine(fin);
        if (!fin)
            return fin;
    } while (str[0]=='#');

    Float_t v[4], mag;
    Int_t n = sscanf(str.Data(), "%f %f %f %f %*f %*f %*f %*f %*f %*f %f", v, v+1, v+2, v+3, &mag);
    if (n<4)
    {
        cout << "Read: ERROR - Not enough numbers" << endl;
        return fin;
    }
    set.fMag = n<5 ? -25 : mag;

    set.fStarAz = v[0]*TMath::DegToRad();
    set.fStarEl = v[1]*TMath::DegToRad();

    set.fRawAz  = v[2]*TMath::DegToRad();
    set.fRawEl  = v[3]*TMath::DegToRad();

    if (fin)
    {
        Double_t res, err;
        res = set.GetResidual(&err);
        cout << "Read: " << v[0] << " " << v[1] << "  :  " << v[2] << " " << v[3] << "  :  " << v[2]-v[0] << " " << v[3]-v[1] << "  :  " << res << " " << err << " " << err/res << endl;
    }

    return fin;
}

class MFit : public TGMainFrame
{
private:
    enum {
        kTbFit = 19, //MBending::GetNumPar(), // FIXME!!!
        kTbLoad,
        kTbSave,
        kTbLoadStars,
        kTbReset,
        kTbResetStars
    };

    MGList *fList;

    TList fOriginal;
    TList fCoordinates;
    TList fLabel;

    MBending fBending;

    FontStruct_t fFont;

    void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
    {
        f = 0;

        MBending bend;
        bend.SetParameters(par); // Set Parameters [deg] to MBending

        for (int i=0; i<fCoordinates.GetSize(); i++)
        {
            Set set = *(Set*)fCoordinates.At(i);

            set.Adjust(bend);

            Double_t err;// = 0.005; // [deg] = 0.25SE
            Double_t res = set.GetResidual(&err);
            res /= err;

            f += res*res;
        }

        //f /= (fCoordinates.GetSize()-7)*(fCoordinates.GetSize()-7);
        //f /= fCoordinates.GetSize()*fCoordinates.GetSize();
        //cout << f << ": " << fCoordinates.GetSize() << endl;
        f /= fCoordinates.GetSize();
    }

    static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
    {
        ((MFit*)gMinuit->GetObjectFit())->Fcn(npar, gin, f, par, iflag);
    }

    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};

        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]);
    }
    
    void DrawPolLine(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;
        }
        /*
        if (r0<0)
        {
            r0 = -r0;
            phi0 += 180;
        }
        if (r1<0)
        {
            r1 = -r1;
            phi1 += 180;
        }
    
        phi0 = fmod(phi0+360, 360);
        phi1 = fmod(phi1+360, 360);
    
        if (phi1-phi0<-180)
            phi1+=360;
        */
        TLine line;
        line.SetLineWidth(2);
        line.SetLineColor(kBlue);
    
        Double_t p0 = phi0<phi1?phi0:phi1;
        Double_t p1 = phi0<phi1?phi1:phi0;
    
        if (phi0>phi1)
        {
            Double_t d = r1;
            r1 = r0;
            r0 = d;
        }
    
        r0 /= 90;
        r1 /= 90;
    
        Double_t dr = r1-r0;
        Double_t dp = p1-p0;
    
        Double_t x0[3] = { r0*cos(p0*TMath::DegToRad()), r0*sin(p0*TMath::DegToRad()), 0};
    
        for (double i=p0+10; i<p1+10; i+=10)
        {
            if (i>p1)
                i=p1;
    
            Double_t r = dr/dp*(i-p0)+r0;
            Double_t p = TMath::DegToRad()*i;
    
            Double_t x1[3] = { r*cos(p), r*sin(p), 0};
    
            Double_t y0[3], y1[3];
    
            view->WCtoNDC(x0, y0);
            view->WCtoNDC(x1, y1);
    
            line.DrawLine(y0[0], y0[1], y1[0], y1[1]);
    
            x0[0] = x1[0];
            x0[1] = x1[1];
        }
    }
    
    void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1, Float_t angle=0)
    {
        Double_t r0   = set.GetRawZd();
        Double_t phi0 = set.GetRawAz()-angle;
        Double_t r1   = set.GetStarZd();
        Double_t phi1 = set.GetStarAz()-angle;
    
        if (r0<0)
        {
            r0 = -r0;
            phi0 += 180;
        }
        if (r1<0)
        {
            r1 = -r1;
            phi1 += 180;
        }
    
        phi0 = fmod(phi0+360, 360);
        phi1 = fmod(phi1+360, 360);
    
        if (phi1-phi0<-180)
            phi1+=360;
    
        if (scale<0 || scale>1000)
            scale = -1;
    
        if (scale>0)
        {
            Double_t d = r1-r0;
            r0 += scale*d;
            r1 -= scale*d;
            d = phi1-phi0;
            phi0 += scale*d;
            phi1 -= scale*d;
    
            DrawPolLine(pad, r0, phi0, r1, phi1);
            DrawMarker(pad,  r0, phi0, r1, phi1);
        }
        else
            DrawMarker(pad,  r1, phi1, 0 ,0);
    }

    void Fit(Double_t &before, Double_t &after, Double_t &backw)
    {
        if (fOriginal.GetSize()==0)
        {
            cout << "Sorry, no input data loaded..." << endl;
            return;
        }

        fCoordinates.Delete();
        for (int i=0; i<fOriginal.GetSize(); i++)
            fCoordinates.Add(new Set(*(Set*)fOriginal.At(i)));

        cout << "-----------------------------------------------------------------------" << endl;

        gStyle->SetOptStat("emro");

        TH1F hres1("Res1", " Residuals before correction ", fOriginal.GetSize()/3, 0, 0.3);
        TH1F hres2("Res2", " Residuals after correction ",  fOriginal.GetSize()/3, 0, 0.3);
        TH1F hres3("Res3", " Residuals after backward correction ",  fOriginal.GetSize()/3, 0, 0.3);
    
        hres1.SetXTitle("\\Delta [\\circ]");
        hres1.SetYTitle("Counts");
    
        hres2.SetXTitle("\\Delta [\\circ]");
        hres2.SetYTitle("Counts");
    
        hres3.SetXTitle("\\Delta [\\circ]");
        hres3.SetYTitle("Counts");
    
        TGraph gdaz;
        TGraph gdzd;
        TGraph gaz;
        TGraph gzd;
        TGraph graz;
        TGraph grzd;
        TGraph grmag;
        TGraph gmaz;
        TGraph gmzd;
    
        gdaz.SetTitle(" \\Delta Az vs. Zd ");
        gdzd.SetTitle(" \\Delta Zd vs. Az ");
    
        gaz.SetTitle(" \\Delta Az vs. Az ");
        gzd.SetTitle(" \\Delta Zd vs. Zd ");

        gmaz.SetTitle(" \\Delta Az vs. Mag ");
        gmzd.SetTitle(" \\Delta Zd vs. Mag ");

        graz.SetTitle(" \\Delta vs. Az ");
        grzd.SetTitle(" \\Delta vs. Zd ");
        grmag.SetTitle(" \\Delta vs. Mag ");
    
        TMinuit minuit(MBending::GetNumPar());  //initialize TMinuit with a maximum of 5 params
        minuit.SetObjectFit(this);
        minuit.SetPrintLevel(-1);
        minuit.SetFCN(fcn);
    
        fBending.SetMinuitParameters(minuit, MBending::GetNumPar()); // Init Parameters [deg]
    
        for (int i=0; i<MBending::GetNumPar(); i++)
        {
            TGButton *l = (TGButton*)fList->FindWidget(i);
            minuit.FixParameter(i);
            if (l->GetState()==kButtonDown)
                minuit.Release(i);
        }
    
        //minuit.Command("SHOW PARAMETERS");
        //minuit.Command("SHOW LIMITS");
    
        cout << endl;
        cout << "Starting fit..." << endl;
        cout << "For the fit an measurement error in the residual of ";
        cout << "0.02deg (=1SE) is assumed." << endl;
        cout << endl;
    
        Int_t ierflg = 0;
        ierflg = minuit.Migrad();
        cout << "Migrad returns " << ierflg << endl;
        // minuit.Release(2);
        ierflg = minuit.Migrad();
        cout << "Migrad returns " << ierflg << endl << endl;
    
        //
        // Get Fit Results
        //
        fBending.GetMinuitParameters(minuit);
        fBending.PrintMinuitParameters(minuit);
        //fBending.Save("bending_magic.txt");
    
    
        //
        // Make a copy of all list entries
        //
        TList list;
        list.SetOwner();
        for (int i=0; i<fCoordinates.GetSize(); i++)
            list.Add(new Set(*(Set*)fCoordinates.At(i)));
    
        //
        // Correct for Offsets only
        //
        TArrayD par;
        fBending.GetParameters(par);
        for (int i=2; i<MBending::GetNumPar(); i++)
            par[i]=0;
    
        MBending b2;
        b2.SetParameters(par);
    
        //
        // Calculate correction and residuals
        //
        for (int i=0; i<fCoordinates.GetSize(); i++)
        {
            Set &set0 = *(Set*)fCoordinates.At(i);
    
            ZdAz za = set0.GetStarZdAz()*kRad2Deg;
    
            //
            // Correct for offsets only
            //
            Set set1(set0);
            set1.Adjust(b2);
    
            hres1.Fill(set1.GetResidual());
    
            set0.Adjust(fBending);
            hres2.Fill(set0.GetResidual());
    
            Double_t dz = fmod(set0.GetDAz()+720, 360);
            if (dz>180)
                dz -= 360;
    
            gdzd.SetPoint(i, za.Az(), set0.GetDZd());
            gdaz.SetPoint(i, za.Zd(), dz);
            graz.SetPoint(i, za.Az(), set0.GetResidual());
            grzd.SetPoint(i, za.Zd(), set0.GetResidual());
            gaz.SetPoint( i, za.Az(), dz);
            gzd.SetPoint( i, za.Zd(), set0.GetDZd());
            if (set0.GetMag()>=-20)
            {
                grmag.SetPoint(i, set0.GetMag(), set0.GetResidual());
                gmaz.SetPoint( i, set0.GetMag(), dz);
                gmzd.SetPoint( i, set0.GetMag(), set0.GetDZd());
            }
        }
    
        //
        // Check for overflows
        //
        const Stat_t ov = hres2.GetBinContent(hres2.GetNbinsX()+1);
        if (ov>0)
            cout << "WARNING: " << ov << " overflows in residuals." << endl;
    
    
    
        cout << dec << endl;
        cout << "              Number of calls to FCN: " << minuit.fNfcn << endl;
        cout << "Minimum value found for FCN (Chi^2?): " << minuit.fAmin << endl;
        cout << "                  Fit-Probability(?): " << TMath::Prob(minuit.fAmin/*fOriginal.GetSize()*/, fOriginal.GetSize()-minuit.GetNumFreePars())*100 << "%" << endl;
        cout << "                           Chi^2/NDF: " << minuit.fAmin/(fOriginal.GetSize()-minuit.GetNumFreePars()) << endl;
        //cout << "Prob(?): " << TMath::Prob(fChisquare,ndf);
    
    
    
        //
        // Print all data sets for which the backward correction is
        // twice times worse than the residual gotten from the
        // bending correction itself
        //
        cout << endl;
        cout << "Checking backward correction (raw-->star):" << endl;
        for (int i=0; i<fCoordinates.GetSize(); i++)
        {
            Set set0(*(Set*)list.At(i));
            Set &set1 = *(Set*)list.At(i);
    
            set0.AdjustBack(fBending);
            set1.Adjust(fBending);
    
            const Double_t res0 = set0.GetResidual();
            const Double_t res1 = set1.GetResidual();
            const Double_t diff = TMath::Abs(res0-res1);
    
            hres3.Fill(res0);
    
            if (diff<hres2.GetMean()*0.66)
                continue;
    
            cout << "DBack: " << setw(6) << set0.GetStarZd() << " " << setw(7) << set0.GetStarAz() << ":  ";
            cout << "ResB="<< setw(7) << res0*60 << "  ResF=" << setw(7) << res1*60 << "  |ResB-ResF|=" << setw(7) << diff*60 << " arcmin" << endl;
        }
        cout << "OK." << endl;
        cout << endl;
    

        TCanvas *c1;

        if ((c1 = (TCanvas*)gROOT->FindObject("CanvGraphs")))
            delete c1;
        if ((c1 = (TCanvas*)gROOT->FindObject("CanvResiduals")))
            delete c1;


        const Double_t max1 = TMath::Max(gaz.GetHistogram()->GetMaximum(), gdaz.GetHistogram()->GetMaximum());
        const Double_t max2 = TMath::Max(gzd.GetHistogram()->GetMaximum(), gdzd.GetHistogram()->GetMaximum());
        const Double_t max3 = TMath::Max(grzd.GetHistogram()->GetMaximum(), graz.GetHistogram()->GetMaximum());

        const Double_t min1 = TMath::Min(gaz.GetHistogram()->GetMinimum(), gdaz.GetHistogram()->GetMinimum());
        const Double_t min2 = TMath::Min(gzd.GetHistogram()->GetMinimum(), gdzd.GetHistogram()->GetMinimum());
        const Double_t min3 = TMath::Min(grzd.GetHistogram()->GetMinimum(), graz.GetHistogram()->GetMinimum());

        const Double_t absmax1 = TMath::Max(max1, TMath::Abs(min1));
        const Double_t absmax2 = TMath::Max(max2, TMath::Abs(min2));
        const Double_t absmax3 = TMath::Max(max3, TMath::Abs(min3));

        gaz.SetMaximum(absmax1);
        gzd.SetMaximum(absmax2);
        gdaz.SetMaximum(absmax1);
        gdzd.SetMaximum(absmax2);
        gmaz.SetMaximum(absmax1);
        gmzd.SetMaximum(absmax2);
        graz.SetMaximum(absmax3);
        grzd.SetMaximum(absmax3);
        grmag.SetMaximum(absmax3);
        gaz.SetMinimum(-absmax1);
        gzd.SetMinimum(-absmax2);
        gdaz.SetMinimum(-absmax1);
        gdzd.SetMinimum(-absmax2);
        gmaz.SetMinimum(-absmax1);
        gmzd.SetMinimum(-absmax2);
        graz.SetMinimum(0);
        grzd.SetMinimum(0);
        grmag.SetMinimum(0);

        c1=new TCanvas("CanvGraphs", "Graphs");
        c1->Divide(3,3,0,0);

        TLine line;
        line.SetLineColor(kGreen);
        line.SetLineWidth(2);

        c1->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        TGraph *g=(TGraph*)gaz.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Az [\\circ]");
        g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);

        c1->cd(2);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)gdaz.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Zd [\\circ]");
        g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
        cout << "Mean dAz: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) <<  endl;

        c1->cd(3);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)gmaz.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Mag");
        g->GetHistogram()->SetYTitle("\\Delta Az [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
    
        c1->cd(4);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)gdzd.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Az [\\circ]");
        g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);
        cout << "Mean dZd: " << g->GetMean(2) << " \xb1 " << g->GetRMS(2) <<  endl;
        cout << endl;
    
        c1->cd(5);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)gzd.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Zd [\\circ]");
        g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);

        c1->cd(6);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)gmzd.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Mag");
        g->GetHistogram()->SetYTitle("\\Delta Zd [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
        line.DrawLine(g->GetXaxis()->GetXmin(), -360./16384, g->GetXaxis()->GetXmax(), -360./16384);

        c1->cd(7);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)graz.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Az [\\circ]");
        g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);
    
        c1->cd(8);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)grzd.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Zd [\\circ]");
        g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);

        c1->cd(9);
        gPad->SetBorderMode(0);
        gPad->SetGridx();
        gPad->SetGridy();
        g=(TGraph*)grmag.DrawClone("A*");
        g->SetBit(kCanDelete);
        g->GetHistogram()->SetXTitle("Mag");
        g->GetHistogram()->SetYTitle("\\Delta [\\circ]");
        line.DrawLine(g->GetXaxis()->GetXmin(),  360./16384, g->GetXaxis()->GetXmax(),  360./16384);


        //
        // Print out the residual before and after correction in several
        // units
        //
        cout << fCoordinates.GetSize() << " data sets." << endl << endl;
        cout << "Total Spread of Residual:" << endl;
        cout << "-------------------------" << endl;
        cout << "before: " << Form("%6.4f", hres1.GetMean()) << " \xb1 " << Form("%6.4f", hres1.GetRMS()) << " deg \t";
        cout << "before: " << Form("%4.1f", hres1.GetMean()*60) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60) << " arcmin" << endl;
        cout << "after:  " << Form("%6.4f", hres2.GetMean()) << " \xb1 " << Form("%6.4f", hres2.GetRMS()) << " deg \t";
        cout << "after:  " << Form("%4.1f", hres2.GetMean()*60) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60) << " arcmin" << endl;
        cout << "backw:  " << Form("%6.4f", hres3.GetMean()) << " \xb1 " << Form("%6.4f", hres3.GetRMS()) << " deg \t";
        cout << "backw:  " << Form("%4.1f", hres3.GetMean()*60) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60) << " arcmin" << endl;
        cout << endl;
        cout << "before: " << Form("%4.1f", hres1.GetMean()*16348/360) << " \xb1 " << Form("%.1f", hres1.GetRMS()*16384/360) << " SE \t\t";
        cout << "before: " << Form("%4.1f", hres1.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres1.GetRMS()*60*60/23.4) << " pix" << endl;
        cout << "after:  " << Form("%4.1f", hres2.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres2.GetRMS()*16384/360) << " SE \t\t";
        cout << "after:  " << Form("%4.1f", hres2.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres2.GetRMS()*60*60/23.4) << " pix" << endl;
        cout << "backw:  " << Form("%4.1f", hres3.GetMean()*16384/360) << " \xb1 " << Form("%.1f", hres3.GetRMS()*16384/360) << " SE \t\t";
        cout << "backw:  " << Form("%4.1f", hres3.GetMean()*60*60/23.4) << " \xb1 " << Form("%.1f", hres3.GetRMS()*60*60/23.4) << " pix" << endl;
        cout << endl;
        cout << endl; // ±
    

        before = hres1.GetMean()*16384/360;
        after  = hres2.GetMean()*16384/360;
        backw  = hres3.GetMean()*16384/360;

    
        gStyle->SetOptStat(1110);
        gStyle->SetStatFormat("6.2g");


        c1=new TCanvas("CanvResiduals", "Residuals", 800, 800);
        c1->Divide(2, 2, 0, 0);
    
        c1->cd(2);
        gPad->SetBorderMode(0);
        hres1.SetLineColor(kRed);
        hres1.DrawCopy();

        line.DrawLine(360./16384, 0, 360./16384, hres1.GetMaximum());
    
        c1->cd(4);
        gPad->SetBorderMode(0);
        hres2.SetLineColor(kBlue);
        TH1 *h=hres2.DrawCopy();
        TF1 f("mygaus", "(gaus)", 0, 1);
        f.SetLineColor(kMagenta/*6*/);
        f.SetLineWidth(1);
        f.SetParameter(0, h->GetBinContent(1));
        f.FixParameter(1, 0);
        f.SetParameter(2, h->GetRMS());
        h->Fit("mygaus", "QR");
        hres3.SetLineColor(kCyan);
        hres3.SetLineStyle(kDashed);
        hres3.DrawCopy("same");
        cout << "Gaus-Fit  Sigma: " << f.GetParameter(2) << "\xb0" << endl;
        cout << "Fit-Probability: " << f.GetProb()*100 << "%" << endl;
        cout << "      Chi^2/NDF: " << f.GetChisquare() << "/" << f.GetNDF() << " = " << f.GetChisquare()/f.GetNDF() << endl;
        line.DrawLine(360./16384, 0, 360./16384, h->GetMaximum());

        c1->cd(1);
        gPad->SetBorderMode(0);
        gPad->SetTheta(90);
        gPad->SetPhi(90);
        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();
        for (int i=0; i<fOriginal.GetSize(); i++)
            DrawSet(gPad, *(Set*)fOriginal.At(i));//, 10./hres1.GetMean());
    
        c1->cd(3);
        gPad->SetBorderMode(0);
        gPad->SetTheta(90);
        gPad->SetPhi(90);
        h2res1.SetTitle(" Arb. Residuals after correction (scaled) ");
        h2res1.DrawCopy("surf1pol");
        gPad->Modified();
        gPad->Update();
//        for (int i=0; i<fCoordinates.GetSize(); i++)
//            DrawSet(gPad, *(Set*)fCoordinates.At(i), 10./hres2.GetMean(), par[0]);

        RaiseWindow();
    }

    void LoadStars(TString fname="tpoint.txt")
    {
        const Int_t size = fOriginal.GetSize();

        ifstream fin(fname);

        while (fin && fin.get()!='\n');
        while (fin && fin.get()!='\n');
        while (fin && fin.get()!='\n');
        if (!fin)
        {
            cout << "File '" << fname << "' not found!" << endl;
            return;
        }

        while (1)
        {
            Set set;

            fin >> set;  // Read data from file [deg], it is stored in [rad]
            if (!fin)
                break;

            fOriginal.Add(new Set(set));
        }

        cout << "Found " << fOriginal.GetSize()-size << " sets of coordinates ";
        cout << "(Total=" << fOriginal.GetSize() << ")" << endl;
    }

    // --------------------------------------------------------------------------
    //
    //  Opens an open dialog
    //
    TString OpenDialog(EFileDialogMode mode=kFDOpen)
    {
        static const char *gOpenTypes[] =
        {
            "TPoint files", "*.txt",
            "All files",    "*",
            NULL,           NULL
        };

        static TString dir(".");

        TGFileInfo fi; // fFileName and fIniDir deleted in ~TGFileInfo

        fi.fFileTypes = (const char**)gOpenTypes;
        fi.fIniDir    = StrDup(dir);

        new TGFileDialog(fClient->GetRoot(), this, mode, &fi);

        if (!fi.fFilename)
            return 0;

        dir = fi.fIniDir;

        return fi.fFilename;
    }

    Bool_t ProcessMessage(Long_t msg, Long_t mp1, Long_t mp2)
    {
        // cout << "Msg: " << hex << GET_MSG(msg) << endl;
        // cout << "SubMsg: " << hex << GET_SUBMSG(msg) << dec << endl;
        switch (GET_MSG(msg))
        {
        case kC_COMMAND:
            switch (GET_SUBMSG(msg))
            {
            case kCM_BUTTON:
                switch (mp1)
                {
                case kTbFit:
                    {
                        Double_t before=0;
                        Double_t after=0;
                        Double_t backw=0;
                        Fit(before, after, backw);
                        DisplayBending();
                        DisplayResult(before, after, backw);
                    }
                    return kTRUE;
                case kTbLoad:
                    fBending.Load(OpenDialog());
                    DisplayBending();
                    return kTRUE;
                case kTbSave:
                    fBending.Save(OpenDialog(kFDSave));
                    return kTRUE;
                case kTbLoadStars:
                    LoadStars(OpenDialog());
                    DisplayData();
                    return kTRUE;
                case kTbReset:
                    fBending.Clear();
                    DisplayBending();
                    return kTRUE;
                case kTbResetStars:
                    fOriginal.Delete();
                    DisplayData();
                    return kTRUE;
                }
                return kTRUE;
            }
            return kTRUE;
        }
        return kTRUE;
    }

    void AddTextButton(TGCompositeFrame *f, TString txt, Int_t id=-1, TGLayoutHints *h=0)
    {
        TGButton *but = new TGTextButton(f, txt, id);
        but->Associate(this);
        f->AddFrame(but, h);
        fList->Add(but);

    }

    void AddCheckButton(TGCompositeFrame *f, TString txt, Int_t id=-1, TGLayoutHints *h=0)
    {
        TGButton *but = new TGCheckButton(f, txt, id);
        but->Associate(this);
        f->AddFrame(but, h);
        fList->Add(but);
    }

    TGLabel *AddLabel(TGCompositeFrame *f, TString txt, TGLayoutHints *h=0)
    {
        TGLabel *l = new TGLabel(f, txt/*, TGLabel::GetDefaultGC()(), fFont*/);
        f->AddFrame(l, h);
        fList->Add(l);
        fLabel.Add(l);
        return l;
    }

    void DisplayBending()
    {
        TArrayD par, err;
        fBending.GetParameters(par);
        fBending.GetError(err);

        TGLabel *l;

        for (int i=0; i<MBending::GetNumPar(); i++)
        {
            l = (TGLabel*)fLabel.At(i);
            l->SetText(Form("%.4f\xb0", par[i]));

            l = (TGLabel*)fLabel.At(MBending::GetNumPar()+i);
            l->SetText(Form("\xb1 %8.4f\xb0", err[i]));
        }
    }

    void DisplayData()
    {
        TGLabel *l = (TGLabel*)fLabel.At(3*MBending::GetNumPar());
        l->SetText(Form("%d data sets loaded.", fOriginal.GetSize()));
    }

    void DisplayResult(Double_t before, Double_t after, Double_t backw)
    {
        TGLabel *l1 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+1);
        l1->SetText(Form("Before: %.1f +- %.1f SE", before, 0));

        TGLabel *l2 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+2);
        l2->SetText(Form("After:  %.1f +- %.1f SE", after, 0));

        TGLabel *l3 = (TGLabel*)fLabel.At(3*MBending::GetNumPar()+3);
        l3->SetText(Form("Backw:  %.1f +- %.1f SE", backw, 0));
    }

public:
    ~MFit()
    {
        if (fFont)
            gVirtualX->DeleteFont(fFont);
    }
    MFit() : TGMainFrame(gClient->GetRoot(), 750, 370, kHorizontalFrame)
    {
        fCoordinates.SetOwner();
        fOriginal.SetOwner();

        fList = new MGList;
        fList->SetOwner();

        fFont = gVirtualX->LoadQueryFont("7x13bold");

        TGLayoutHints *hints0 = new TGLayoutHints(kLHintsExpandY, 7, 5, 5, 6);
        TGLayoutHints *hints1 = new TGLayoutHints(kLHintsExpandX|kLHintsExpandY, 5, 7, 5, 6);
        TGLayoutHints *hints2 = new TGLayoutHints(kLHintsCenterX|kLHintsCenterY, 5, 5, 5, 5);
        fList->Add(hints0);
        fList->Add(hints1);
        fList->Add(hints2);

        TGGroupFrame *grp1 = new TGGroupFrame(this, "Control", kVerticalFrame);
        AddFrame(grp1, hints0);
        fList->Add(grp1);

        TGGroupFrame *grp2 = new TGGroupFrame(this, "Parameters", kHorizontalFrame);
        AddFrame(grp2, hints1);
        fList->Add(grp2);



        TGLayoutHints *hints4 = new TGLayoutHints(kLHintsExpandX, 5, 5,  5);
        TGLayoutHints *hints5 = new TGLayoutHints(kLHintsExpandX, 5, 5, 15);
        AddTextButton(grp1, "Load Pointing Model", kTbLoad,       hints5);
        AddTextButton(grp1, "Save Pointing Model", kTbSave,       hints4);
        AddTextButton(grp1, "Fit Parameters",      kTbFit,        hints5);
        AddTextButton(grp1, "Reset Parameters",    kTbReset,      hints4);
        AddTextButton(grp1, "Load Stars",          kTbLoadStars,  hints5);
        AddTextButton(grp1, "Reset Stars",         kTbResetStars, hints4);
        fList->Add(hints4);
        fList->Add(hints5);


        TGHorizontalFrame *comp = new TGHorizontalFrame(grp2, 1, 1);
        grp2->AddFrame(comp);
        fList->Add(comp);

        TGLayoutHints *hints3 = new TGLayoutHints(kLHintsLeft|kLHintsTop, 0, 20, 5, 0);
        fList->Add(hints3);

        TGLabel *l;

        TGVerticalFrame *vframe = new TGVerticalFrame(comp, 1, 1);
        comp->AddFrame(vframe, hints3);
        fList->Add(vframe);

        for (int i=0; i<MBending::GetNumPar(); i++)
            AddCheckButton(vframe, fBending.GetName(i), i);

        vframe = new TGVerticalFrame(comp, 1, 1);
        comp->AddFrame(vframe, hints3);
        fList->Add(vframe);

        l = new TGLabel(vframe, "+000.0000");
        l->SetTextJustify(kTextRight);
        fList->Add(l);
        fLabel.Add(l);

        TGButton *but = (TGButton*)fList->FindWidget(0);

        TGLayoutHints *h = new TGLayoutHints(kLHintsCenterY, 0, 0, but->GetHeight()-l->GetHeight());
        fList->Add(h);

        vframe->AddFrame(l,h);

        for (int i=1; i<MBending::GetNumPar(); i++)
            AddLabel(vframe, "+000.0000", h)->SetTextJustify(kTextRight);

        vframe = new TGVerticalFrame(comp, 1, 1);
        comp->AddFrame(vframe, hints3);
        fList->Add(vframe);

        for (int i=0; i<MBending::GetNumPar(); i++)
            AddLabel(vframe, "\xb1 00.0000\xb0", h)->SetTextJustify(kTextRight);

        vframe = new TGVerticalFrame(comp, 1, 1);
        comp->AddFrame(vframe, hints3);
        fList->Add(vframe);

        for (int i=0; i<MBending::GetNumPar(); i++)
            AddLabel(vframe, fBending.GetDescription(i), h);

        l = new TGLabel(grp1, "0000000 Data Sets loaded.");
        grp1->AddFrame(l, hints5);
        fList->Add(l);
        fLabel.Add(l);

        l = new TGLabel(grp1, "");
        l->SetTextJustify(kTextLeft);
        grp1->AddFrame(l, hints5);
        fList->Add(l);
        fLabel.Add(l);

        l = new TGLabel(grp1, "");
        l->SetTextJustify(kTextLeft);
        grp1->AddFrame(l, hints5);
        fList->Add(l);
        fLabel.Add(l);

        l = new TGLabel(grp1, "");
        l->SetTextJustify(kTextLeft);
        grp1->AddFrame(l, hints5);
        fList->Add(l);
        fLabel.Add(l);


        ((TGCheckButton*)fList->FindWidget(0))->SetState(kButtonDown);
        ((TGCheckButton*)fList->FindWidget(1))->SetState(kButtonDown);
        ((TGCheckButton*)fList->FindWidget(6))->SetState(kButtonDown);
        ((TGCheckButton*)fList->FindWidget(10))->SetState(kButtonDown);
        ((TGCheckButton*)fList->FindWidget(12))->SetState(kButtonDown);
        ((TGCheckButton*)fList->FindWidget(17))->SetState(kButtonDown);

        SetWindowName("TPoint Fitting Window");
        SetIconName("TPoint++");

        Layout();

        MapSubwindows();
        MapWindow();

        DisplayBending();
        DisplayData();
    }
    ClassDef(MFit, 0)
};

ClassImp(MFit);

void gui()
{
    gErrorIgnoreLevel = kError;
    new MFit;
    // TF1 f1("f1", "[0]/cos((90-x)*3.1415/180)", 0, 90)
}