/* ======================================================================== *\
!
! *
! * 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): Abelardo Moralejo 10/2004 <mailto:moralejo@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MFSize                                                                //
//                                                                         //
//  A filter to select events as a function of Size. Depending on the Size //
//  value, a certain fraction of events is rejected. A histogram (fProb)   //
//  must be provided, containing the acceptance probability of events as a //
//  function of log10(Size) (which must be the x axis of the histogram).   //
//  There is no restriction on the histogram binning or ranges (events     //
//  with Size outside histogram ranges will simply be accepted). Histogram //
//  bin contents should be between 0. and 1. If they are not, negative     //
//  will be equivalent to 0 (events always rejected), and values above 1   //
//  will be equivalent to 1 (events always accepted). A warning will be    //
//  displayed by MFSize::SetProb in that case.                             //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MFSize.h"

#include <fstream>
#include <TRandom.h>

#include "MParList.h"

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

#include "MHillas.h"

ClassImp(MFSize);

using namespace std;

// --------------------------------------------------------------------------
//
//     Constructor
//
MFSize::MFSize(const char *name, const char *title): fNumSelectedEvts(0)
{
  fName  = name  ? name  : "MFSize";
  fTitle = title ? title : "Filter to select events according to Size";

  //
  // Default fProb histogram: 20 bins in the log10(Size) range [2,4]
  // Values are calculated to make the Size distribution of MC gammas
  // in the standard MC sample to look like the distribution for MC 
  // hadrons (p+He)
  //
  fProb = new TH1F("ProbVsLog10Size", "", 20, 2., 4.);

  Float_t frac[20] = {0.144928, 0.171911, 0.18763, 0.209461, 0.253219, 
 		      0.305425, 0.384593, 0.485204, 0.612452, 0.708138, 
 		      0.754522, 0.728028, 0.774046, 0.791422, 0.808775, 
 		      0.896842, 1., 1., 1., 1.};
  
  for (Int_t i = 0; i < 20; i++)
    fProb->SetBinContent(i+1, frac[i]);
}

// --------------------------------------------------------------------------
//
//   SetProb: Sets the histogram containing the probability acceptance for
//   events as a function of log10(Size). Checks that the values are between
//   0 and 1.
//

void MFSize::SetProb(TH1F* h)
{
  fProb = h;

  for (Int_t i = 1; i <= fProb->GetNbinsX(); i++)
    {
      if (fProb->GetBinContent(i) > 1.)
	*fLog << warn << dbginf << fName << 
	  "Probability value above 1 found in histogram. " <<
	  "Value will be interpreted as 1." << endl;

      else if(fProb->GetBinContent(i) < 0. )
	*fLog << warn << dbginf << fName << 
	  "Probability value below 0 found in histogram. " <<
	  "Value will be interpreted as 0." << endl;
    }
}

// --------------------------------------------------------------------------
//
//   Preprocess
//  
// Look for the MHillas container
//
Int_t MFSize::PreProcess(MParList *pList)
{
  fHillas = (MHillas*)pList->FindObject(AddSerialNumber("MHillas"));
    if (!fHillas)
      {
	*fLog << err << dbginf << fName << AddSerialNumber("MHillas") <<
	  " not found... aborting." << endl;
        return kFALSE;
      }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Process  
//
Int_t MFSize::Process()
{
  fResult = kTRUE;

  Int_t ibin = fProb->FindBin(log10(fHillas->GetSize()));

  //
  // If value is outside histogram range, accept event
  //
  if (ibin > fProb->GetNbinsX() || ibin < 1)
    {
      fNumSelectedEvts++;
      return kTRUE;
    }

  Float_t SelFrac = fProb->GetBinContent(ibin);

  const Float_t Nrnd = gRandom->Uniform();

  fResult = SelFrac > Nrnd;

  if (!fResult)
    return kTRUE;

  fNumSelectedEvts++;
  return kTRUE;
}

