/* ======================================================================== *\
!
! *
! * 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, 03/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2007
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//   MFMagicCuts
//
//  Predefinitions:
//  ---------------
//      width/length = MHillas.fWidth/MHillas.fLength
//      area         = MHillas.GetArea*fMm2Deg*fMm2Deg
//      lgsize       = log10(MHillas.fSize)
//
//      leakage1     = log10(MNewImagePar.fLeakage1+1)
//
//      alpha        = MHillasSrc.fAlpha
//      dist         = MHillasSrc.fDist*fMm2Deg
//      m3long       = MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*fMm2Deg
//      slope        = MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/fMm2Deg
//
//      antialpha    = MHillasSrcAnti.fAlpha
//      antidist     = MHillasSrcAnti.fDist*fMm2Deg
//      antim3long   = MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*fMm2Deg
//      antislope    = MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/fMm2Deg
//
//      had          = Hadronness.fVal
//
//  Coefficients/Cuts:
//  ------------------
//
//    c[0],  c[5], c[6], c[7],  c[8], c[9]:
//          xi          = c[0]+c[6]*pow(leakage1, c[7]);
//          p           = xi*(1-width/length);
//          sign1       = m3long-c[5]
//          sign2       = (dist-c[9])*c[8]-slope
//          disp        = sign1<0 ||sign2<0 ? -disp : disp
//          antisign1   = antim3long-c[5]
//          antisign2   = (antidist-c[9])*c[8]-antislope
//          antidisp    = antisign1<0 || antisign2<0 ? -antidisp : antidisp
//          thetasq     = disp^2 + dist^2 - 2*disp*dist*alpha
//          antithetasq = antidisp^2 + antidist^2 - 2*antidisp*antidist*antialpha
//
//    c[1]:
//          thetasq < c[1]*c[1]
//
//    AreaCut:
//    c[2], c[3], c[4]:
//          A = c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]))
//          area < A
//
//    HadronnessCut:
//    c[10], c[11]:
//          had       <= c[10]
//          10^lgsize >= c[11]
//
//
//  Options:
//  --------
//
//   HadronnessCut:
//    - none/nocut: No area cut
//    - area:       Area cut
//    - all:        same as area
//
//   ThetaCut:
//    - none/nocut: no theta cut
//    - on:         cut in theta only
//    - off:        anti-cut in anti-theta only
//    - wobble:     same as on|off (both)
//
//
//  Input Container:
//  ------
//   MGeomCam
//   MHillas
//   MHillasExt
//   MNewImagePar
//   MHillasSrc
//   [MHillasSrcAnti [MHillasSrc]]
//
//  Output:
//  -------
//   ThetaSquared [MParameterD]    // stores thetasq
//   Disp [MParameterD]            // stores -p*sign(MHillasSrc.fCosDeltaAlpha)
//   ThetaSquaredCut [MParameterD] // stores c[1]*c[1]
//
/////////////////////////////////////////////////////////////////////////////
#include "MFMagicCuts.h"

#include <fstream>
#include <errno.h> 

#include "MHMatrix.h"

#include "MLog.h"
#include "MLogManip.h"

#include "MMath.h"

#include "MParList.h"

#include "MGeomCam.h"
#include "MHillas.h"
#include "MHillasSrc.h"
#include "MHillasExt.h"
#include "MNewImagePar.h"
#include "MParameters.h"

ClassImp(MFMagicCuts);

using namespace std;

// --------------------------------------------------------------------------
//
// constructor
//
MFMagicCuts::MFMagicCuts(const char *name, const char *title)
    : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fNewImgPar(0),
    fThetaSq(0), fDisp(0), fHadronness(0), fMatrix(0), fVariables(30),
    fThetaCut(kNone), fHadronnessCut(kArea)
{
    fName  = name  ? name  : "MFMagicCuts";
    fTitle = title ? title : "Class to evaluate the MagicCuts";

    fMatrix = NULL;

    AddToBranchList("MHillas.fWidth");
    AddToBranchList("MHillas.fLength");
    AddToBranchList("MHillas.fSize");
    AddToBranchList("MHillasSrc.fAlpha");
    AddToBranchList("MHillasSrc.fDist");
    AddToBranchList("MHillasSrc.fCosDeltaAlpha");
    AddToBranchList("MHillasSrcAnti.fAlpha");
    AddToBranchList("MHillasSrcAnti.fDist");
    AddToBranchList("MHillasSrcAnti.fCosDeltaAlpha");
    AddToBranchList("MHillasExt.fM3Long");
    AddToBranchList("MHillasExt.fSlopeLong");
    AddToBranchList("MNewImagePar.fLeakage1");
    AddToBranchList("Hadronness.fVal");
/*
    fVariables[0] =  1.42547;      // Xi
    fVariables[1] =  0.233773;     // Theta^2
    fVariables[2] =  0.2539460;    // Area[0]
    fVariables[3] =  5.2149800;    // Area[1]
    fVariables[4] =  0.1139130;    // Area[2]
    fVariables[5] = -0.0889;       // M3long
    fVariables[6] =  0.18;         // Xi(theta)
    fVariables[7] =  0.18;         // Xi(theta)
    */
}

// --------------------------------------------------------------------------
//
// The theta cut value GetThetaSqCut() is propageted to the parameter list
// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else.
//
Int_t MFMagicCuts::PreProcess(MParList *pList)
{
    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!cam)
    {
        *fLog << err << "MGeomCam not found... aborting." << endl;
        return kFALSE;
    }

    fMm2Deg = cam->GetConvMm2Deg();

    fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
    if (!fThetaSq)
        return kFALSE;
    fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
    if (!fDisp)
        return kFALSE;

    // propagate Theta cut to the parameter list
    // so that it can be used somewhere else
    MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
    if (!thetacut)
        return kFALSE;
    thetacut->SetVal(GetThetaSqCut());

    Print();

    if (fMatrix)
        return kTRUE;

    //-----------------------------------------------------------

    fHil = (MHillas*)pList->FindObject("MHillas");
    if (!fHil)
    {
        *fLog << err << "MHillas not found... aborting." << endl;
        return kFALSE;
    }

    fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
    if (!fHilExt)
    {
        *fLog << err << "MHillasExt not found... aborting." << endl;
        return kFALSE;
    }

    fNewImgPar = (MNewImagePar*)pList->FindObject("MNewImagePar");
    if (!fNewImgPar)
    {
        *fLog << err << "MNewImagePar not found... aborting." << endl;
        return kFALSE;
    }

    fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
    if (!fHilSrc)
    {
        *fLog << err << "MHillasSrc not found... aborting." << endl;
        return kFALSE;
    }

    if (fThetaCut&kOff)
    {
        fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
        if (!fHilAnti)
        {
            *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
            return kFALSE;
        }
    }

    if (fHadronnessCut&kHadronness)
    {
        fHadronness = (MParameterD*)pList->FindObject("Hadronness", "MParameterD");
        if (!fHadronness)
        {
            *fLog << warn << "Hadronness [MParameterD] not found... aborting." << endl;
            return kFALSE;
        }
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Returns the mapped value from the Matrix
//
Double_t MFMagicCuts::GetVal(Int_t i) const
{
    return (*fMatrix)[fMap[i]];
}

// --------------------------------------------------------------------------
//
// You can use this function if you want to use a MHMatrix instead of the
// given containers. This function adds all necessary columns to the
// given matrix. Afterward you should fill the matrix with the corresponding
// data (eg from a file by using MHMatrix::Fill). If you now loop
// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
// will take the values from the matrix instead of the containers.
//
void MFMagicCuts::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
      return;

    fMatrix = mat;

//    fMap[kELength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
//    fMap[kEWidth]  = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");

    fMap[kEWdivL]   = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
    fMap[kESize]    = fMatrix->AddColumn("log10(MHillas.fSize)");
    fMap[kEArea]    = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");

    fMap[kELeakage] = fMatrix->AddColumn("log10(MNewImagePar.fLeakage1+1)");

    fMap[kEAlpha]   = fMatrix->AddColumn("MHillasSrc.fAlpha");
    fMap[kEDist]    = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
    fMap[kEM3Long]  = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");

    fMap[kESrcSign] = fMatrix->AddColumn("sign(MHillasSrc.fCosDeltaAlpha)");

    fMap[kESlope]   = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");

    if (fThetaCut&kOff)
    {
        fMap[kEAlphaAnti]  = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
        fMap[kEDistAnti]   = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
        fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
        fMap[kESlopeAnti]  = fMatrix->AddColumn("MHillasExt.fSlopeLong*sign(MHillasSrcAnti.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg");
    }

    if (fHadronnessCut&kHadronness)
        fMap[kEHadronness] = fMatrix->AddColumn("Hadronness.fVal");
}

// --------------------------------------------------------------------------
//
// Returns theta squared based on the formula:
//    p := c*(1-width/length)
//    return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
//
// If par!=NULL p is stored in this MParameterD
//
Double_t MFMagicCuts::GetThetaSq(Double_t p, Double_t wl, Double_t d, Double_t a) const
{
    // wl = width/l [1]
    // d  = dist    [deg]
    // a  = alpha   [deg]
    return d*d + p*p - 2*d*p*TMath::Cos(a*TMath::DegToRad());
}

// --------------------------------------------------------------------------
//
// Returns squared cut value in theta: fVariables[1]^2
//
Double_t MFMagicCuts::GetThetaSqCut() const
{
    return fVariables[1]*fVariables[1];
}

// ---------------------------------------------------------------------------
//
// Evaluate dynamical magic-cuts
//
Int_t MFMagicCuts::Process()
{
    // Get some variables
    const Double_t wdivl  = fMatrix ? GetVal(kEWdivL)   : fHil->GetWidth()/fHil->GetLength();
    const Double_t lgsize = fMatrix ? GetVal(kESize)    : TMath::Log10(fHil->GetSize());
    const Double_t leak   = fMatrix ? GetVal(kELeakage) : TMath::Log10(fNewImgPar->GetLeakage1()+1);

    // For simplicity
    const Double_t *c = fVariables.GetArray();

    // Value for Theta cut (Disp parametrization)
    const Double_t cut  = GetThetaSqCut();
    const Double_t xi   = c[0]+c[6]*pow(leak, c[7]);
    const Double_t disp = xi*(1-wdivl);

    // Default if we return before the end
    fResult = kFALSE;

    // ------------------- Most efficient cut -----------------------
    // ---------------------- Theta cut ON --------------------------
    const Double_t dist   = fMatrix ? GetVal(kEDist)    : fHilSrc->GetDist()*fMm2Deg;
    const Double_t alpha  = fMatrix ? GetVal(kEAlpha)   : fHilSrc->GetAlpha();
    const Double_t sign   = fMatrix ? GetVal(kESrcSign) : MMath::Sgn(fHilSrc->GetCosDeltaAlpha());
    const Double_t m3long = fMatrix ? GetVal(kEM3Long)  : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
    const Double_t slope  = fMatrix ? GetVal(kESlope)   : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());

        // Do Ghostbusting with slope and m3l
    const Double_t sign1 = (dist-c[9])*c[8]-slope;
    const Double_t sign2 = m3long-c[5];
    const Double_t p = sign1<0 || sign2<0 ? -disp : disp;

    // Align disp along source direction depending on third moment
    //const Double_t p = TMath::Sign(disp, m3long-c[5]);

    // Align sign of disp along shower axis like m3long
    fDisp->SetVal(-p*sign);

    // Calculate corresponding Theta^2
    const Double_t thetasq = GetThetaSq(p, wdivl, dist, alpha);
    fThetaSq->SetVal(thetasq);

    // Perform theta cut
    if (fThetaCut&kOn && thetasq >= cut)
        return kTRUE;

    // --------------------- Efficient  cut -------------------------
    // ---------------------- Area/M3l cut --------------------------
    if (fHadronnessCut&kArea)
    {
        const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
        const Double_t A    = lgsize>c[3] ? c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
        if (area>=A)
            return kTRUE;
    }

    // ---------------------------------------------------------------
    // ---------------------------------------------------------------

    if (fHadronnessCut&kHadronness)
    {
        const Double_t had = fMatrix ? GetVal(kEHadronness) : fHadronness->GetVal();
        if (had>c[10])
            return kTRUE;

        if (TMath::Power(10, lgsize)<c[11])
            return kTRUE;
    }

    // ------------------- Least efficient cut -----------------------
    // ---------------------- Theta cut OFF --------------------------

    if (fThetaCut&kOff)
    {
        const Double_t distanti    = fMatrix ? GetVal(kEDistAnti)   : fHilAnti->GetDist()*fMm2Deg;
        const Double_t alphaanti   = fMatrix ? GetVal(kEAlphaAnti)  : fHilAnti->GetAlpha();
        const Double_t m3lanti     = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
        const Double_t slopeanti   = fMatrix ? GetVal(kESlopeAnti)  : fHilExt->GetSlopeLong()/TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());

        // Do Ghostbusting with slope and m3l
        const Double_t sign3 = (distanti-c[9])*c[8]-slopeanti;
        const Double_t sign4 = m3lanti-c[5];
        const Double_t panti = sign3<0 || sign4<0 ? -disp : disp;

        // Align disp along source direction depending on third moment
        //const Double_t panti       = TMath::Sign(disp, m3lanti-c[5]);
        const Double_t thetasqanti = GetThetaSq(panti, wdivl, distanti, alphaanti);

        if (thetasqanti<cut)
            return kTRUE;

        // This cut throws away gamma-like events which would be assigned to the
        // anti-source. It is not necessary but it offers the possibility to
        // fill the MHDisp plot in one run (Having a hole instead of eg. Crab
        // the data can be rotated and subtracted)
    }

    fResult = kTRUE;

    return kTRUE;
}

void MFMagicCuts::SetVariables(const TArrayD &arr)
{
    fVariables = arr;
}

TString MFMagicCuts::GetParam(Int_t i) const
{
    if (i>=fVariables.GetSize())
        return "";

    return Form("Param[%2d]=%+f", i, fVariables[i]);
}

void MFMagicCuts::Print(Option_t *o) const
{
    *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;

    *fLog << "Using Theta cut: ";
    switch (fThetaCut)
    {
    case kNone:
        *fLog << "none" << endl;
        break;
    case kOn:
        *fLog << "on" << endl;
        break;
    case kOff:
        *fLog << "off" << endl;
        break;
    case kWobble:
        *fLog << "on and off (wobble)" << endl;
        break;
    }
    *fLog << "Using hadronness cut: ";
    switch (fHadronnessCut)
    {
    case kNoCut:
        *fLog << "none" << endl;
        break;
    case kArea:
        *fLog << "area" << endl;
        break;
    case kHadronness:
        *fLog << "hadronness" << endl;
        break;
    case kAll:
        *fLog << "all" << endl;
        break;
    }
    *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;

    const Int_t n = fVariables.GetSize();
    for (int i=0; i<n; i+=3)
    {
        *fLog << setw(25) << GetParam(i);
        *fLog << setw(25) << GetParam(i+1);
        *fLog << setw(25) << GetParam(i+2);
        *fLog << endl;
    }
    *fLog << setiosflags(ios::right);
}

Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
{
    ifstream fin(fname);
    if (!fin)
    {
        *fLog << err << "Cannot open file " << fname << ": ";
        *fLog << strerror(errno) << endl;
        return kFALSE;
    }

    for (int i=0; i<fVariables.GetSize(); i++)
    {
        fin >> fVariables[i];
        if (!fin)
        {
            *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
            return kERROR;
        }
    }
    return kTRUE;
}

Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
{
    ofstream fout(fname);
    if (!fout)
    {
        *fLog << err << "Cannot open file " << fname << ": ";
        *fLog << strerror(errno) << endl;
        return kFALSE;
    }

    for (int i=0; i<fVariables.GetSize(); i++)
        fout << setw(11) << fVariables[i] << endl;

    return kTRUE;
}

Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    if (MFilter::ReadEnv(env, prefix, print)==kERROR)
        return kERROR;


    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "ThetaCut", print))
    {
        TString str = GetEnvValue(env, prefix, "ThetaCut", "");
        str.ToLower();
        str=str.Strip(TString::kBoth);

        if (str==(TString)"nocut" || str==(TString)"none")
            fThetaCut = kNone;
        if (str==(TString)"on")
            fThetaCut = kOn;
        if (str==(TString)"off")
            fThetaCut = kOff;
        if (str==(TString)"wobble")
            fThetaCut = kWobble;
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "HadronnessCut", print))
    {
        TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
        str.ToLower();
        str=str.Strip(TString::kBoth);

        if (str==(TString)"nocut" || str==(TString)"none")
            fHadronnessCut = kNoCut;
        if (str==(TString)"area")
            fHadronnessCut = kArea;
        if (str==(TString)"hadronness")
            fHadronnessCut = kHadronness;
        if (str==(TString)"all")
            fHadronnessCut = kAll;

        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "File", print))
    {
        const TString fname = GetEnvValue(env, prefix, "File", "");
        if (!CoefficentsRead(fname))
            return kERROR;

        return kTRUE;
    }

    for (int i=0; i<fVariables.GetSize(); i++)
    {
        if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
        {
            // It is important that the last argument is a floating point
            fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
            rc = kTRUE;
        }
    }
    return rc;
    //return kTRUE; // means: can use default values
    //return rc;  // means: require something in resource file
}
