/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz, 2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MBending
//
//    Double_t fIe   ; // [rad] Index Error in Elevation
//    Double_t fIa   ; // [rad] Index Error in Azimuth
//
//    Double_t fFlop ; // [rad] Vertical Sag
//     * do not use if not data: Zd<0
//
//    Double_t fNpae ; // [rad] Az-El Nonperpendicularity
//
//    Double_t fCa   ; // [rad] Left-Right Collimation Error
//
//    Double_t fAn   ; // [rad] Azimuth Axis Misalignment (N-S)
//    Double_t fAw   ; // [rad] Azimuth Axis Misalignment (E-W)
//
//    Double_t fTf   ; // [rad] Tube fluxture (sin)
//     * same as ecec if no data: Zd<0
//    Double_t fTx   ; // [rad] Tube fluxture (tan)
//     * do not use with NPAE if no data: Zd<0
//
//    Double_t fNrx  ; // [rad] Nasmyth rotator displacement, horizontan
//    Double_t fNry  ; // [rad] Nasmyth rotator displacement, vertical
//
//    Double_t fCrx  ; // [rad] Alt/Az Coude Displacement (N-S)
//    Double_t fCry  ; // [rad] Alt/Az Coude Displacement (E-W)
//
//    Double_t fEces ; // [rad] Elevation Centering Error (sin)
//    Double_t fAces ; // [rad] Azimuth Centering Error (sin)
//    Double_t fEcec ; // [rad] Elevation Centering Error (cos)
//    Double_t fAcec ; // [rad] Azimuth Centering Error (cos)
//
////////////////////////////////////////////////////////////////////////////
#include "MBending.h"

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

#include <TMinuit.h>

#include "MTime.h"

ClassImp(MBending);

using namespace std;

#undef DEBUG
//#define DEBUG(txt) txt
#define DEBUG(txt)

const Int_t MBending::fNumPar=19;

void MBending::Init()
{
    fCoeff = new Double_t*[fNumPar];
    fName  = new TString[fNumPar];
    fDescr = new TString[fNumPar];

    fCoeff[ 0] = &fIa;      fName[ 0] = "IA";
    fCoeff[ 1] = &fIe;      fName[ 1] = "IE";
    fCoeff[ 2] = &fFlop;    fName[ 2] = "FLOP";
    fCoeff[ 3] = &fAn;      fName[ 3] = "AN";
    fCoeff[ 4] = &fAw;      fName[ 4] = "AW";
    fCoeff[ 5] = &fNpae;    fName[ 5] = "NPAE";
    fCoeff[ 6] = &fCa;      fName[ 6] = "CA";
    fCoeff[ 7] = &fTf;      fName[ 7] = "TF";
    fCoeff[ 8] = &fTx;      fName[ 8] = "TX";
    fCoeff[ 9] = &fEces;    fName[ 9] = "ECES";
    fCoeff[10] = &fAces;    fName[10] = "ACES";
    fCoeff[11] = &fEcec;    fName[11] = "ECEC";
    fCoeff[12] = &fAcec;    fName[12] = "ACEC";
    fCoeff[13] = &fNrx;     fName[13] = "NRX";
    fCoeff[14] = &fNry;     fName[14] = "NRY";
    fCoeff[15] = &fCrx;     fName[15] = "CRX";
    fCoeff[16] = &fCry;     fName[16] = "CRY";
    fCoeff[17] = &fMagic1;  fName[17] = "MAGIC1";
    fCoeff[18] = &fMagic2;  fName[18] = "MAGIC2";

    fDescr[ 0] =  "Index Error Azimuth";
    fDescr[ 1] =  "Index Error Zenith Distance";
    fDescr[ 2] =  "Vertical Sag";
    fDescr[ 3] =  "Azimuth Axis Misalignment (N-S)";
    fDescr[ 4] =  "Azimuth Axis Misalignment (E-W)";
    fDescr[ 5] =  "Az-El Nonperpendicularity";
    fDescr[ 6] =  "Left-Right Collimation Error";
    fDescr[ 7] =  "Tube fluxture (sin)";
    fDescr[ 8] =  "Tube fluxture (tan)";
    fDescr[ 9] =  "Elevation Centering Error (sin)";
    fDescr[10] =  "Azimuth Centering Error (sin)";
    fDescr[11] =  "Elevation Centering Error (cos)";
    fDescr[12] =  "Azimuth Centering Error (cos)";
    fDescr[13] =  "Nasmyth rotator displacement (horizontal)";
    fDescr[14] =  "Nasmyth rotator displacement (vertical)";
    fDescr[15] =  "Alt/Az Coude Displacement (N-S)";
    fDescr[16] =  "Alt/Az Coude Displacement (E-W)";
    fDescr[17] =  "n/a <ZA Hysteresis>";
    fDescr[18] =  "n/a";
}
void MBending::Reset()
{
    Clear();
}

void MBending::Load(const char *name)
{
    /*
     ! MMT 1987 July 8
     ! T   36   7.3622   41.448  -0.0481
     !   IA        -37.5465    20.80602
     !   IE        -13.9180     1.25217
     !   NPAE       +7.0751    26.44763
     !   CA         -6.9149    32.05358
     !   AN         +0.5053     1.40956
     !   AW         -2.2016     1.37480
     ! END
     */

    ifstream fin(name);
    if (!fin)
    {
        cout << "Error: Cannot open file '" << name << "'" << endl;
        return;
    }

    char c;
    while (fin && fin.get()!='\n');
    fin >> c;

    if (c!='S' && c!='s')
    {
        cout << "Error: This in not a model correcting the star position (" << c << ")" << endl;
        return;
    }

    Clear();

    cout << endl;

    Double_t val;
    fin >> val;
    cout << "Number of observed stars: " << val << endl;
    fin >> val;
    cout << "Sky RMS: " << val << "\"" << endl;
    fin >> val;
    cout << "Refraction Constant A: " << val << "\"" << endl;
    fin >> val;
    cout << "Refraction Constant B: " << val << "\"" << endl;

    cout << endl;

    cout << "  & = Name            Value                  Sigma" << endl;
    cout << "--------------------------------------------------" << endl;

    while (fin)
    {
        TString str;
        fin >> str;

        if (str=="END")
            break;

        if (str[0]=='&')
        {
            cout << " & ";
            str.Remove(0);
        }
        else
            cout << "   ";

        if (str[1]=='=')
        {
            cout << "=  ";
            str.Remove(0);
        }
        else
            cout << "   ";

        fin >> val;
        cout << str << "\t" << setw(11) << val << "°     \t";
        val *= kDeg2Rad;

        // Find parameter
        Int_t n = -1;
        for (int i=0; i<fNumPar; i++)
            if (str==fName[i])
            {
                n = i;
                *fCoeff[i] = val;
                break;
            }

        fin >> val;
        cout << setw(9) << val << "°" << endl;

        // corresponding error
        fError[n] = val*kDeg2Rad;
    }
    cout << endl;
}

void MBending::Save(const char *name)
{
    /*
     ! MMT 1987 July 8
     ! T   36   7.3622   41.448  -0.0481
     !   IA        -37.5465    20.80602
     !   IE        -13.9180     1.25217
     !   NPAE       +7.0751    26.44763
     !   CA         -6.9149    32.05358
     !   AN         +0.5053     1.40956
     !   AW         -2.2016     1.37480
     ! END
     */

    ofstream fout(name);
    if (!fout)
    {
        cout << "Error: Cannot open file '" << name << "'" << endl;
        return;
    }

    MTime t;
    t.Now();

    fout << "MAGIC1 " << t << endl;
    fout << "S   00   000000   000000  0000000" << endl;
    fout << setprecision(8);
    for (int i=0; i<fNumPar; i++)
    {
        fout << " " << setw(6) << GetName(i) << " ";
        fout << setw(13) << *fCoeff[i]*kRad2Deg << "   ";
        fout << setw(11) << fError[i]*kRad2Deg << endl;
    }
    fout << "END" << endl;
}

Double_t MBending::Sign(Double_t val, Double_t alt)
{
    // Some pointing corrections are defined as Delta ZA, which
    // is (P. Wallace) defined [0,90]deg while Alt is defined
    // [0,180]deg
    return (TMath::Pi()/2-alt < 0 ? -val : val);
}

AltAz MBending::AddOffsets(const AltAz &aa) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p = aa;

    const AltAz I(fIe, fIa);
    p += I;

    return p;
}

AltAz MBending::SubtractOffsets(const AltAz &aa) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p = aa;

    const AltAz I(fIe, fIa);
    p -= I;

    return p;
}

AltAz MBending::CalcAnAw(const AltAz &p, Int_t sign) const
{
    // Corrections for AN and AW without approximations
    // as done by Patrick Wallace. The approximation cannot
    // be used for MAGIC because the correctioon angle
    // AW (~1.5deg) is not small enough.

    // Vector in cartesian coordinates
    TVector3 v1;

    // Set Azimuth and Elevation
    v1.SetMagThetaPhi(1, TMath::Pi()/2-p.Alt(), p.Az());


    TVector3 v2(v1);
//    cout << sign << endl;

//    cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;

    // Rotate around the x- and y-axis
    v1.RotateY(sign*fAn);
    v1.RotateX(sign*fAw);

//    cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
//    cout << "v2: " << v2.Theta()*TMath::RadToDeg() << " " << v2.Theta()*TMath::RadToDeg() << endl;

   // cout << "dv: " << (v2.Theta()-v1.Theta())*TMath::RadToDeg() << " " << (v2.Phi()-v1.Phi())*TMath::RadToDeg() << endl;

    Double_t dalt = v1.Theta()-v2.Theta();
    Double_t daz  = v1.Phi()  -v2.Phi();

    //cout << dalt*TMath::RadToDeg() << " " << daz*TMath::RadToDeg() << endl;

    if (daz>TMath::Pi())
        daz -= TMath::TwoPi();
    if (daz<-TMath::Pi())
        daz += TMath::TwoPi();

//    if (daz>TMath::Pi()/2)
//    {
//    }

    AltAz d(dalt, daz);
    return d;

    // Calculate Delta Azimuth and Delta Elevation
    /*
    AltAz d(TMath::Pi()/2-v1.Theta(), v1.Phi());

    cout << "p :  " << p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
    cout << "d :  " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
    d -= p;
    cout << "d-p: " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
    d *= sign;
    cout << "d* : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;


    cout << "p2:  " << 90-p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
    cout << "d2:  " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;

    Int_t s1 = 90-d.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
    Int_t s2 = 90-p.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;


    if (s1 != s2)
    {
        //90-d.Alt() <-- -90+d.Alt()

        d.Alt(d.Alt()-TMath::Pi());
        cout << "Alt-" << endl;
    }
    cout << "d': " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
 /*
    // Fix 'direction' of output depending on input vector
    if (TMath::Pi()/2-sign*p.Alt()<0)
    {
        d.Alt(d.Alt()-TMath::Pi());
        cout << "Alt-" << endl;
    }
    //if (TMath::Pi()/2-sign*p.Alt()>TMath::Pi())
    //{
    //    d.Alt(TMath::Pi()-d.Alt());
    //    cout << "Alt+" << endl;
    //}

    // Align correction into [-180,180]
    while (d.Az()>TMath::Pi())
    {
        d.Az(d.Az()-TMath::Pi()*2);
        cout << "Az-" << endl;
    }
    while (d.Az()<-TMath::Pi())
    {
        d.Az(d.Az()+TMath::Pi()*2);
        cout << "Az+" << endl;
    }
   */
    return d;
}


AltAz MBending::Correct(const AltAz &aa) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p = aa;

    DEBUG(cout << setprecision(16));
    DEBUG(cout << "Bend7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()),  fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
    const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
    p += CRX;
    p += CRY;

    DEBUG(cout << "Bend6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
    const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
    p += NRX;
    p += NRY;

    DEBUG(cout << "Bend5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
    const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
    p += CES;
    p += CEC;

    DEBUG(cout << "Bend4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
    const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
    //p += TX;
    p += TF;


    DEBUG(cout << "Bend3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    /*
     //New Corrections for NPAE and CA:
     TVector3 v(1.,1.,1.); // Vector in cartesian coordinates

     //Set Azimuth and Elevation
     v.SetPhi(p.Az());
     v.SetTheta(TMath::Pi()/2-p.Alt());
     //Rotation Vectors:
     TVector3 vNpae(             cos(p.Az()),              sin(p.Az()),             0);
     TVector3   vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
     //Rotate around the vectors vNpae and vCa
     v.Rotate(fNpae, vNpae);
     v.Rotate(fCa,     vCa);

     p.Az(v.Phi());
     p.Alt(TMath::Pi()/2-v.Theta());
     */

    //Old correction terms for Npae and Ca:
    const AltAz CA(0, -fCa/cos(p.Alt()));
    p += CA;

    const AltAz NPAE(0, -fNpae*tan(p.Alt()));
    p += NPAE;

    DEBUG(cout << "Bend2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz ANAW(CalcAnAw(p, -1));
    p += ANAW;

    /* Old correction terms for An and Aw:
     const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
     const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
     p += AW;
     p += AN;
    */

    DEBUG(cout << "Bend1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
    p += FLOP;

    const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
    p += MAGIC1;

    const AltAz I(fIe, fIa);
    p += I;

    DEBUG(cout << "Bend0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    return p;
}

AltAz MBending::CorrectBack(const AltAz &aa) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p = aa;

    DEBUG(cout << setprecision(16));
    DEBUG(cout << "Back0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz I(fIe, fIa);
    p -= I;

    const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
    p -= MAGIC1;

    //const AltAz MAGIC1(fMagic1*sin(p.Az()), 0);
    //p -= MAGIC1;

    const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
    p -= FLOP;

    DEBUG(cout << "Back1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    /* Old correction terms for An and Aw:
     const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
     const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
     p -= AN;
     p -= AW;
     */

    const AltAz ANAW(CalcAnAw(p, -1));
    p -= ANAW;

    DEBUG(cout << "Back2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    //Old Correction terms for Npae and Ca:
    const AltAz NPAE(0, -fNpae*tan(p.Alt()));
    p -= NPAE;

    const AltAz CA(0, -fCa/cos(p.Alt()));
    p -= CA;

    /*
     //New Correction term for Npae and Ca:
     TVector3 v2(1.,1.,1.); // Vector in cartesian coordinates
     //Set Azimuth and Elevation
     v2.SetPhi(p.Az());
     v2.SetTheta(TMath::Pi()/2-p.Alt());
     //Rotation Vectors:
     TVector3 vNpae(             cos(p.Az()),              sin(p.Az()),             0);
     TVector3   vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
     //Rotate around the vectors vCa and vNpae
     v2.Rotate(-fCa,     vCa);
     v2.Rotate(-fNpae, vNpae);

     p.Az(v2.Phi());
     p.Alt(TMath::Pi()/2-v2.Theta());
    */

    DEBUG(cout << "Back3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
    const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
    p -= TF;
    //p -= TX;

    DEBUG(cout << "Back4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
    const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
    p -= CEC;
    p -= CES;

    DEBUG(cout << "Back5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
    const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
    p -= NRY;
    p -= NRX;

    DEBUG(cout << "Back6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
    const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()),  fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
    p -= CRY;
    p -= CRX;

    DEBUG(cout << "Back7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);

    return p;
}

ZdAz MBending::Correct(const ZdAz &zdaz) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
    AltAz c = Correct(p);
    return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
}

ZdAz MBending::CorrectBack(const ZdAz &zdaz) const
{
    // Correct [rad]
    // zdaz    [rad]
    AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
    AltAz c = CorrectBack(p);
    return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
}

void MBending::SetParameters(const Double_t *par, Int_t n)
{
    Clear();

    while (n--)
        *fCoeff[n] = par[n]/kRad2Deg;
}

void MBending::GetParameters(Double_t *par, Int_t n) const
{
    while (n--)
        par[n] = *fCoeff[n]*kRad2Deg;
}

void MBending::SetMinuitParameters(TMinuit &m, Int_t n) const
{
    if (n<0)
        n = fNumPar;

    Int_t ierflg = 0;

    while (n--)
        m.mnparm(n, fName[n], *fCoeff[n]*kRad2Deg,  1, -360, 360, ierflg);
}

void MBending::GetMinuitParameters(TMinuit &m, Int_t n)
{
    if (n<0 || n>m.GetNumPars())
        n = m.GetNumPars();

    while (n--)
    {
        m.GetParameter(n, *fCoeff[n], fError[n]);
        *fCoeff[n] /= kRad2Deg;
        fError[n]  /= kRad2Deg;
    }
}
/*
void FormatPar(TMinuit &m, Int_t n)
{
    Double_t par, err;
    m.GetParameter(n, par, err);

    int expp = (int)log10(par);
    int expe = (int)log10(err);

    if (err<2*pow(10, expe))
        expe--;

    Int_t exp = expe>expp ? expp : expe;

    par = (int)(par/pow(10, exp)) * pow(10, exp);
    err = (int)(err/pow(10, exp)) * pow(10, exp);

    cout << par << " +- " << err << flush;
}
*/
void MBending::PrintMinuitParameters(TMinuit &m, Int_t n) const
{
    if (n<0)
        n = m.GetNumPars();

    cout << setprecision(3);

    Double_t par, err;

    while (n--)
    {
        m.GetParameter(n, par, err);
        cout << Form(" %2d %6s: ", n, (const char*)fName[n]);
        cout << setw(8) << par << " \xb1 " << setw(6) <<  err << endl;
    }
}