/* ======================================================================== *\
!
! *
! * 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): Wolfgang Wittek, 04/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MCT1SupercutsCalc                                                     //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MCT1SupercutsCalc.h"

#include <math.h>
#include <fstream>

#include "MParList.h"
#include "MHillasExt.h"
#include "MHillasSrc.h"
#include "MMcEvt.hxx"
#include "MCerPhotEvt.h"
#include "MGeomCam.h"
#include "MHadronness.h"
#include "MHMatrix.h"

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

ClassImp(MCT1SupercutsCalc);

using namespace std;

void MCT1SupercutsCalc::InitParams()
{
    fLengthUp.Set(8);
    fLengthLo.Set(8);
    fWidthUp.Set(8);
    fWidthLo.Set(8);
    fDistUp.Set(8);
    fDistLo.Set(8);
    fAsymUp.Set(8);
    fAsymLo.Set(8);
    fAlphaUp.Set(8);

    //---------------------------------
    // cut parameters
    fLengthUp[0] = 0.315585; 
    fLengthUp[1] = 0.001455; 
    fLengthUp[2] = 0.203198; 
    fLengthUp[3] = 0.005532; 
    fLengthUp[4] =-0.001670;
    fLengthUp[5] =-0.020362; 
    fLengthUp[6] = 0.007388; 
    fLengthUp[7] =-0.013463;

    fWidthUp[0] = 0.145412; 
    fWidthUp[1] =-0.001771; 
    fWidthUp[2] = 0.054462; 
    fWidthUp[3] = 0.022280; 
    fWidthUp[4] =-0.009893;
    fWidthUp[5] = 0.056353; 
    fWidthUp[6] = 0.020711; 
    fWidthUp[7] =-0.016703;

    fDistUp[0] = 1.787943; 
    fDistUp[1] = 0.; 
    fDistUp[2] = 2.942310; 
    fDistUp[3] = 0.199815; 
    fDistUp[4] = 0.; 
    fDistUp[5] = 0.249909;
    fDistUp[6] = 0.189697; 
    fDistUp[7] = 0.;

    fLengthLo[0] = 0.151530; 
    fLengthLo[1] = 0.028323; 
    fLengthLo[2] = 0.510707; 
    fLengthLo[3] = 0.053089; 
    fLengthLo[4] = 0.013708;
    fLengthLo[5] = 2.357993; 
    fLengthLo[6] = 0.000080; 
    fLengthLo[7] =-0.007157;

    fWidthLo[0] = 0.089187; 
    fWidthLo[1] =-0.006430; 
    fWidthLo[2] = 0.074442; 
    fWidthLo[3] = 0.003738;
    fWidthLo[4] =-0.004256; 
    fWidthLo[5] =-0.014101; 
    fWidthLo[6] = 0.006126; 
    fWidthLo[7] =-0.002849;

    fDistLo[0] = 0.589406;
    fDistLo[1] = 0.;
    fDistLo[2] =-0.083964;
    fDistLo[3] =-0.007975;
    fDistLo[4] = 0.;
    fDistLo[5] = 0.045374;
    fDistLo[6] =-0.001750;
    fDistLo[7] = 0.;

    fAsymUp[0] = 0.061267; 
    fAsymUp[1] = 0.014462; 
    fAsymUp[2] = 0.014327; 
    fAsymUp[3] = 0.014540; 
    fAsymUp[4] = 0.013391;
    fAsymUp[5] = 0.012319; 
    fAsymUp[6] = 0.010444; 
    fAsymUp[7] = 0.008328;

    fAsymLo[0] =-0.012055; 
    fAsymLo[1] = 0.009157; 
    fAsymLo[2] = 0.005441; 
    fAsymLo[3] = 0.000399; 
    fAsymLo[4] = 0.001433;
    fAsymLo[5] =-0.002050; 
    fAsymLo[6] =-0.000104; 
    fAsymLo[7] =-0.001188;

    fAlphaUp[0] = 13.123440; 
    fAlphaUp[1] = 0.; 
    fAlphaUp[2] = 0.; 
    fAlphaUp[3] = 0.; 
    fAlphaUp[4] = 0.; 
    fAlphaUp[5] = 0.; 
    fAlphaUp[6] = 0.; 
    fAlphaUp[7] = 0.;
    //---------------------------------
}

// --------------------------------------------------------------------------
//
// Set the parameter values from vector 'par'
//
// Attention : it is assumed that there are (9*ncutpar) values
//
void MCT1SupercutsCalc::SetParams(Double_t *par)
{
    UInt_t ncutpar = fLengthUp.GetSize();
    UInt_t k0 = 0;

    TArrayD lup(ncutpar,     par + k0);
    SetLengthUp(lup);
    k0 += ncutpar;

    TArrayD wup(ncutpar,     par + k0);
    SetWidthUp(wup);
    k0 += ncutpar;

    TArrayD dup(ncutpar,     par + k0);
    SetDistUp(dup);
    k0 += ncutpar;

    TArrayD llo(ncutpar,     par + k0);
    SetLengthLo(llo);
    k0 += ncutpar;

    TArrayD wlo(ncutpar,     par + k0);
    SetWidthLo(wlo);
    k0 += ncutpar;

    TArrayD dlo(ncutpar,     par + k0);
    SetDistLo(dlo);
    k0 += ncutpar;

    TArrayD aup(ncutpar,     par + k0);
    SetAsymUp(aup);
    k0 += ncutpar;

    TArrayD alo(ncutpar,     par + k0);
    SetAsymLo(alo);
    k0 += ncutpar;

    TArrayD alphaup(ncutpar, par + k0);
    SetAlphaUp(alphaup);
}

// --------------------------------------------------------------------------
//
// Get the parameter values 
//
// Attention : it is assumed that there are (9*ncutpar) values
//
void MCT1SupercutsCalc::GetParams(Double_t *par)
{
    UInt_t ncutpar = fLengthUp.GetSize();
    UInt_t k0 = 0;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fLengthUp[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fWidthUp[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fDistUp[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fLengthLo[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fWidthLo[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fDistLo[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fAsymUp[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fAsymLo[j];
    }
    k0 += ncutpar;

    for (UInt_t j=0; j<ncutpar; j++)
    {
      UInt_t k = k0 + j;
      par[k] = fAlphaUp[j];
    }
}

// --------------------------------------------------------------------------
//
MCT1SupercutsCalc::MCT1SupercutsCalc(const char *hilname, 
                                     const char *hilsrcname, const char *name, const char *title)
: fHadronnessName("MHadronness"), fHilName(hilname), fHilSrcName(hilsrcname)
{
    fName  = name  ? name  : "MCT1SupercutsCalc";
    fTitle = title ? title : "Class to evaluate the Supercuts";

    InitParams();
}

// --------------------------------------------------------------------------
//
Int_t MCT1SupercutsCalc::PreProcess(MParList *pList)
{
    MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!cam)
    {
        *fLog << err << "MGeomCam (Camera Geometry) not found... aborting." << endl;
        return kFALSE;
    }

    fMm2Deg = cam->GetConvMm2Deg();

    fHadronness = (MHadronness*)pList->FindCreateObj("MHadronness", fHadronnessName);
    if (!fHadronness)
    {
        *fLog << err << fHadronnessName << " [MHadronness] not found... aborting." << endl;
        return kFALSE;
    }

    if (fMatrix)
        return kTRUE;

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

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

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculation of upper and lower limits
//
Double_t MCT1SupercutsCalc::CtsMCut(
#if ROOT_VERSION_CODE > ROOT_VERSION(3,05,00)
const
#endif
                                    TArrayD &a,  Double_t ls, Double_t ct,
                                    Double_t ls2, Double_t dd2)
{
    // define cut-function
    //
    //    dNOMLOGSIZE = 4.1 (=log(60.0)
    //    dNOMCOSZA   = 1.0
    //
    //      a: array of cut parameters
    //     ls: log(SIZE) - dNOMLOGSIZE
    //    ls2: ls^2
    //     ct: Cos(ZA.) - dNOMCOSZA
    //    dd2: DIST^2
    const Double_t limit =
        a[0] + a[1] * dd2 + a[2] * ct  +
        ls  * (a[3] + a[4] * dd2 + a[5] * ct) +
        ls2 * (a[6] + a[7] * dd2);

    //*fLog << "MCT1SupercutsCalc::CtsMCut; *a = "
    //      << *a     << ",  " << *(a+1) << ",  " << *(a+2) << ",  "
    //      << *(a+3) << ",  " << *(a+4) << ",  " << *(a+5) << ",  "
    //      << *(a+6) << ",  " << *(a+7) << endl;

    //*fLog << "MCT1SupercutsCalc::CtsMCut; ls, ls2, ct, dd2, limit = " << ls
    //      << ",  " << ls2 << ",  " << ct << ",  " << dd2 << ",  "
    //      << limit << endl;

    return limit;
}

// --------------------------------------------------------------------------
//
// Returns the mapped value from the Matrix
//
Double_t MCT1SupercutsCalc::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 MCT1SupercutsCalc::InitMapping(MHMatrix *mat)
{
    if (fMatrix)
        return;

    fMatrix = mat;

    fMap[0] = fMatrix->AddColumn("MMcEvt.fTelescopeTheta");
    fMap[1] = fMatrix->AddColumn("MHillas.fWidth");
    fMap[2] = fMatrix->AddColumn("MHillas.fLength");
    fMap[3] = fMatrix->AddColumn("MHillas.fSize");
    fMap[4] = fMatrix->AddColumn("MHillas.fMeanX");
    fMap[5] = fMatrix->AddColumn("MHillas.fMeanY");
    fMap[6] = fMatrix->AddColumn("MHillasSrc.fDist");
    fMap[7] = fMatrix->AddColumn("MHillasSrc.fAlpha");
}

// ---------------------------------------------------------------------------
//
// Evaluate dynamical supercuts for CT1 Mkn421 2001 (Daniel Kranich)
// optimized for mkn 421 2001 data
// 
//          set hadronness to 0.25 if cuts are fullfilled
//                            0.75 otherwise
//
Int_t MCT1SupercutsCalc::Process()
{
    const Double_t kNomLogSize = 4.1;
    const Double_t kNomCosZA   = 1.0;

    const Double_t theta   = fMatrix ? GetVal(0) : fMcEvt->GetTelescopeTheta();
    const Double_t width0  = fMatrix ? GetVal(1) : fHil->GetWidth();
    const Double_t length0 = fMatrix ? GetVal(2) : fHil->GetLength();
    const Double_t size    = fMatrix ? GetVal(3) : fHil->GetSize();
    const Double_t meanx   = fMatrix ? GetVal(4) : fHil->GetMeanX();
    const Double_t meany   = fMatrix ? GetVal(5) : fHil->GetMeanY();
    const Double_t dist0   = fMatrix ? GetVal(6) : fHilSrc->GetDist();

    const Double_t newdist = dist0 * fMm2Deg;

    const Double_t dist2   = meanx*meanx + meany*meany;
    const Double_t dd2     = dist2*fMm2Deg;
    const Double_t dist    = sqrt(dist2) * fMm2Deg;

    const Double_t dmls    = log(size) - kNomLogSize;
    const Double_t dmls2   = dmls * dmls;

    const Double_t dmcza   = cos(theta) - kNomCosZA;

    const Double_t length  = length0 * fMm2Deg;
    const Double_t width   = width0  * fMm2Deg;

    if (newdist < 1.05                                         &&
        newdist < CtsMCut (fDistUp,   dmls, dmcza, dmls2, dd2) &&
        newdist > CtsMCut (fDistLo,   dmls, dmcza, dmls2, dd2) &&
        dist    < 1.05                                         &&
        length  < CtsMCut (fLengthUp, dmls, dmcza, dmls2, dd2) &&
        length  > CtsMCut (fLengthLo, dmls, dmcza, dmls2, dd2) &&
        width   < CtsMCut (fWidthUp,  dmls, dmcza, dmls2, dd2) &&
        width   > CtsMCut (fWidthLo,  dmls, dmcza, dmls2, dd2) &&
        //asym  < CtsMCut (asymup,    dmls, dmcza, dmls2, dd2) &&
        //asym  > CtsMCut (asymlow,   dmls, dmcza, dmls2, dd2) &&
        dist    < CtsMCut (fDistUp,   dmls, dmcza, dmls2, dd2) &&
        dist    > CtsMCut (fDistLo,   dmls, dmcza, dmls2, dd2)  )
        fHadronness->SetHadronness(0.25);
    else
        fHadronness->SetHadronness(0.75);

    fHadronness->SetReadyToSave();

    return kTRUE;
}
//==========================================================================
