/* ======================================================================== *\
!
! *
! * 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): Eva Domingo,     12/2004 <mailto:domingo@ifae.es> 
!              Wolfgang Wittek, 12/2004 <mailto:wittek@mppmu.mpg.de>
!
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MFDisp                                                                  //
//                                                                          //
//  This is a class to set cuts to select an event sample                   //
//  for the Disp optimization                                               //
//                                                                          //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////

#include "MFDisp.h"

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

#include "MParList.h"

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


ClassImp(MFDisp);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MFDisp::MFDisp(const char *name, const char *title)
{
    fName  = name  ? name  : "MFDisp";
    fTitle = title ? title : "Cuts for Disp optimization";

    // default values of cuts
    SetCuts(0, 1000,   0, 1000,   0, 1000, 
            0.0, 1.e10,  0.0, 1.0,  0.0, 1.0,   0.0, 0.0);
}


// --------------------------------------------------------------------------
//
// Set the values for the cuts 
// 
void MFDisp::SetCuts(Int_t islandsmin,      Int_t islandsmax,
                     Int_t usedpixelsmin,   Int_t usedpixelsmax, 
                     Int_t corepixelsmin,   Int_t corepixelsmax,
                     Float_t sizemin,       Float_t sizemax,
                     Float_t leakage1min,   Float_t leakage1max,
                     Float_t leakage2min,   Float_t leakage2max,
                     Float_t lengthmin,     Float_t widthmin)
{ 
  fIslandsMin    = islandsmin;
  fIslandsMax    = islandsmax;

  fUsedPixelsMin = usedpixelsmin;
  fUsedPixelsMax = usedpixelsmax;

  fCorePixelsMin = corepixelsmin;
  fCorePixelsMax = corepixelsmax;

  fSizeMin       = sizemin;
  fSizeMax       = sizemax;

  fLeakage1Min   = leakage1min;
  fLeakage1Max   = leakage1max;

  fLeakage2Min   = leakage2min;
  fLeakage2Max   = leakage2max;

  fLengthMin     = lengthmin;
  fWidthMin      = widthmin;
}


// --------------------------------------------------------------------------
//
Int_t MFDisp::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();


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

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

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


    memset(fCut, 0, sizeof(fCut));


    //--------------------
    *fLog << inf << "MFDisp cut values :" << endl; 

    *fLog << inf << "     fIslandsMin, fIslandsMax = "
          << fIslandsMin << ",  " << fIslandsMax << endl;

    *fLog << inf << "     fUsedPixelsMin, fUsedPixelsMax = "
          << fUsedPixelsMin << ",  " << fUsedPixelsMax << endl;

    *fLog << inf << "     fCorePixelsMin, fCorePixelsMax = "
          << fCorePixelsMin << ",  " << fCorePixelsMax << endl;

    *fLog << inf << "     fSizeMin, fSizeMax = " 
          << fSizeMin << ",  " << fSizeMax << endl;

    *fLog << inf << "     fLeakage1Min, fLeakage1Max = " 
          << fLeakage1Min << ",  " << fLeakage1Max << endl;

    *fLog << inf << "     fLeakage2Min, fLeakage2Max = " 
          << fLeakage2Min << ",  " << fLeakage2Max << endl;

    *fLog << inf << "     fLengthMin, fWidthMin = " 
          << fLengthMin << ",  " << fWidthMin << endl;
    //--------------------

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Sum up the number of events passing each cut
//
Bool_t MFDisp::Set(Int_t rc)
{
    fResult = kTRUE;
    fCut[rc]++;
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Evaluate cuts for Disp optimization
// if selections are fulfilled set fResult = kTRUE;
// 
Int_t MFDisp::Process()
{
    const Double_t length     = fHil->GetLength() * fMm2Deg;
    const Double_t width      = fHil->GetWidth()  * fMm2Deg;
    //const Double_t dist       = fHilSrc->GetDist()* fMm2Deg;
    //const Double_t delta      = fHil->GetDelta()  * kRad2Deg;
    const Double_t size       = fHil->GetSize();

    const Double_t leakage1   = fNewImgPar->GetLeakage1();
    const Double_t leakage2   = fNewImgPar->GetLeakage2();

    const Int_t numusedpixels = fNewImgPar->GetNumUsedPixels();
    const Int_t numcorepixels = fNewImgPar->GetNumCorePixels();
    //const Int_t numislands    = fNewImgPar->GetNumIslands();

    fResult = kFALSE;

    //if  (numislands<fIslandsMin  || numislands>fIslandsMax )
    //    return Set(1);

    if  (numusedpixels<fUsedPixelsMin  || numusedpixels>fUsedPixelsMax )
        return Set(2);

    if  (numcorepixels<fCorePixelsMin  || numcorepixels>fCorePixelsMax )
        return Set(3);

    if (size<fSizeMin  ||  size>fSizeMax)
        return Set(4);

    if (leakage1<fLeakage1Min || leakage1>fLeakage1Max)
        return Set(5);

    if (leakage2<fLeakage2Min || leakage2>fLeakage2Max)
        return Set(6);

    if (length<fLengthMin || width<fWidthMin)
        return Set(7);

    fCut[0]++;
    return kTRUE;
}


// --------------------------------------------------------------------------
//
//  Prints some statistics about the cuts selections
//
Int_t MFDisp::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;

    *fLog << dec << setfill(' ');
    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3);
    *fLog << (int)(fCut[1]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: No.of islands < " << fIslandsMin ;
    *fLog << " or > " << fIslandsMax << endl;

    *fLog << " " << setw(7) << fCut[2] << " (" << setw(3);
    *fLog << (int)(fCut[2]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: No.of used pixels < " << fUsedPixelsMin ;
    *fLog << " or > " << fUsedPixelsMax << endl;

    *fLog << " " << setw(7) << fCut[3] << " (" << setw(3);
    *fLog << (int)(fCut[3]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: No.of core pixels < " << fCorePixelsMin ;
    *fLog << " or > " << fCorePixelsMax << endl;

    *fLog << " " << setw(7) << fCut[4] << " (" << setw(3) ;
    *fLog << (int)(fCut[4]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: SIZE < " << fSizeMin;
    *fLog << " or > " << fSizeMax << endl;

    *fLog << " " << setw(7) << fCut[5] << " (" << setw(3) ;
    *fLog << (int)(fCut[5]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: Leakage1 < " << fLeakage1Min;
    *fLog << " or > " << fLeakage1Max << endl;

    *fLog << " " << setw(7) << fCut[6] << " (" << setw(3) ;
    *fLog << (int)(fCut[6]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: Leakage2 < " << fLeakage2Min;
    *fLog << " or > " << fLeakage2Max << endl;

    *fLog << " " << setw(7) << fCut[7] << " (" << setw(3) ;
    *fLog << (int)(fCut[7]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts skipped due to: LENGTH <= " << fLengthMin;
    *fLog << " or WIDTH <= " << fWidthMin << endl;

    *fLog << " " << fCut[0] << " (" ;
    *fLog << (int)(fCut[0]*100/GetNumExecutions()+0.5) ;
    *fLog << "%) Evts survived the Disp cuts!" << endl;
    *fLog << endl;

    return kTRUE;
}