/* ======================================================================== *\
!
! *
! * 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.h>
#include <iomanip.h>
#include <TVector3.h>

#include <TMinuit.h>

#include "timer.h"

ClassImp(MBending);

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

        Double_t *dest=NULL;

        if (str=="IA")   dest = &fIa;
        if (str=="IE")   dest = &fIe;
        if (str=="FLOP") dest = &fFlop;
        if (str=="NPAE") dest = &fNpae;
        if (str=="CA")   dest = &fCa;
        if (str=="AN")   dest = &fAn;
        if (str=="AW")   dest = &fAw;
        if (str=="TF")   dest = &fTf;
        if (str=="TX")   dest = &fTx;
        if (str=="NRX")  dest = &fNrx;
        if (str=="NRY")  dest = &fNry;
        if (str=="CRX")  dest = &fCrx;
        if (str=="CRY")  dest = &fCry;
        if (str=="ECES") dest = &fEces;
        if (str=="ACES") dest = &fAces;
        if (str=="ECEC") dest = &fEcec;
        if (str=="ACEC") dest = &fAcec;

        if (dest)
            *dest = val;

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

        // Find corresponding error
        for (int i=0; i<MBending::GetNumPar(); i++)
            if (dest==fCoeff[i])
            {
                fError[i] = val*kDeg2Rad;
                break;
            }
    }
    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;
    }

    Timer t;
    t.Now();

    fout << "MAGIC1 " << t.GetTimeStr() << 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::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);

    // New Corrections for AW and AN:
    TVector3 v2(1.,1.,1.); // Vector in cartesian coordinates
    // Set Azimuth and Elevation
    v2.SetPhi(p.Az());
    v2.SetTheta(TMath::Pi()/2-p.Alt());
    // Rotate around the x- and y-axis
    v2.RotateX(fAw);
    v2.RotateY(fAn);
    // Get Azimuth and Elevation
    p.Az(v2.Phi());
    p.Alt(TMath::Pi()/2-v2.Theta());

    /*
    // Old correction terms for An and Aw:
    cout << fAw << " " << cos(p.Az()) << " " << tan(p.Alt()) << endl;
    cout << fAw* cos(p.Az())*tan(p.Alt())*180/3.1415 << endl;

    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 I(fIe, fIa);
    p += I;

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

    return p;
}

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::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 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;
     */

    // New Corrections for AN and AW:
    TVector3 v(1.,1.,1.); // Vector in cartesian coordinates
    // Set Azimuth and Elevation
    v.SetPhi(p.Az());
    v.SetTheta(TMath::Pi()/2-p.Alt());
    // Rotate around the x- and y-axis
    v.RotateY(-fAn);
    v.RotateX(-fAw);
    // Get Azimuth and Elevation
    p.Az(v.Phi());
    p.Alt(TMath::Pi()/2-v.Theta());


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

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

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

    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=-1) 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=-1)
{
    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=-1) 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;
    }
}