/* ======================================================================== *\
!
! *
! * 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): Javier Rico     02/2005 <mailto:jrico@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// Computes the Effective areas and coefficients for unfolding for a given
// spectrum that can be parametrized by a function
//
//////////////////////////////////////////////////////////////////////////////

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

#include "MEffAreaAndCoeffCalc.h"

#include "TF1.h"
#include "MHillas.h"
#include "MMcEvt.hxx"
#include "TH2F.h"

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

ClassImp(MEffAreaAndCoeffCalc);

using namespace std;

const Int_t fNTbins = 2;
const Double_t fTbin[fNTbins+1] = {0,10,20};


// -------------------------------------------------------------------------
//
// Constructor
//
MEffAreaAndCoeffCalc::MEffAreaAndCoeffCalc()
  : fSpec(NULL), fHorig(NULL), fEbins(10), fEsubbins(20), fEmin(10.), fEmax(10000.),
    fCoeff(NULL), fEffA(NULL)
{
  // set the default function
  // OJO!! this should be changed
  SetFunction("4.e9*pow(x,-2.6+1)");  

  // create the TChains 
  // OJO!! make sure they read the appropriate branch
  fCini = new TChain("Events");
  fCcut = new TChain("Events");

  // define some usefull aliases
  fCini->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
  fCini->SetAlias("theta","MMcEvt.fTheta*180./3.14159");

  // OJO!! the estimated energy must be changed for a real one
  fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
  fCcut->SetAlias("theta","MMcEvt.fTheta*180./3.14159");

  fCcut->SetBranchStatus("*",0);
  fCcut->SetBranchStatus("MHillas.*",1);
  fCcut->SetBranchStatus("MMcEvt.*",1);
  fCcut->SetBranchAddress("MHillas.",&fHillas);
  fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
}

// -------------------------------------------------------------------------
//
// Destructor
//
MEffAreaAndCoeffCalc::~MEffAreaAndCoeffCalc()
{
  if(fSpec)
    delete fSpec;

  if(fWeight)
    delete [] fWeight;

  if(fCoeff)
    delete fCoeff;

  if(fEffA)
    delete fEffA;

  if(fHorig) 
    delete fHorig;

  delete fCini;
  delete fCcut;
}

// -------------------------------------------------------------------------
//
// Set the function by expression and minimum and maximum energies
//
void MEffAreaAndCoeffCalc::SetFunction(const Char_t* chfunc, Float_t emin, Float_t emax)
{
  if(fSpec)
    delete fSpec;
  if(emin<=0 || emax<=0 || emax<emin)
    {
      emin = fEmin;
      emax = fEmax;
    }
  fSpec = new TF1("fspec",chfunc,emin,emax);
}

// -------------------------------------------------------------------------
//
// Set the function by function pointer
//
void MEffAreaAndCoeffCalc::SetFunction(TF1* newf)
{
  if(fSpec)
    delete fSpec;
  fSpec = newf;
}

// -------------------------------------------------------------------------
//
// fill the histogram containing the original sample energy spectrum
//
void MEffAreaAndCoeffCalc::FillOriginalSpectrum()
{  
  const Double_t logemin = TMath::Log10(fEmin);
  const Double_t logemax = TMath::Log10(fEmax);
  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
  
  if(fHorig) delete fHorig;
  fHorig   = new TH1F("fhorig","Original energy spectrum",esbins,logemin,logemax);
  fCini->Draw("logenergy>>fhorig","","goff");
}

// -------------------------------------------------------------------------
//
// compute the weights for a particular input spectrum
//
void MEffAreaAndCoeffCalc::ComputeWeights()
{  
  const Double_t logemin = TMath::Log10(fEmin);
  const Double_t logemax = TMath::Log10(fEmax);
  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
  const Double_t desub = de/fEsubbins; // subbin size (in log)
  const Int_t esbins = fEbins*fEsubbins; // total number of subbins

  // reset the weights array
  if(fWeight)
    delete [] fWeight;
  fWeight = new Double_t[esbins];
    
  if(!fHorig)
    FillOriginalSpectrum();
  
  for(Int_t i=0;i<esbins;i++)
    {
      const Double_t denom = fHorig->GetBinContent(i+1);              // number of events in initial spectrum
      const Double_t ew    = TMath::Power(10,logemin+(i+0.5)*desub); // real energy
      const Double_t numer = fSpec->Eval(ew);                        // number of events for the required spectrum
      if(denom)
	fWeight[i]=numer/denom;
      else
	{
	  cout << "MEffAreaAndCoeffCalc::ComputeWeights Warning: no statistic to compute weight for energy " << ew << ", setting it to -1 "  << endl;
	  fWeight[i]=-1;
	}
    }
}

// --------------------------------------------------------------
//
// compute the coefficients used for the (iterative) unfolding 
//
void MEffAreaAndCoeffCalc::ComputeCoefficients()
{ 
  if(!fWeight)
    {
      cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: No weights computed! nothing done" << endl;
      return;
    }

  const Double_t logemin = TMath::Log10(fEmin);
  const Double_t logemax = TMath::Log10(fEmax);
  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
  const Double_t desub = de/fEsubbins; // subbin size (in log)
  const Int_t esbins = fEbins*fEsubbins; // total number of subbins
  const Int_t nentries = Int_t(fCcut->GetEntries());
 
  // declare needed histos
  TH2F* hest = new TH2F("hest","Estimated energy",fEbins,logemin,logemax,fNTbins,fTbin);
  TH2F* hmc  = new TH2F("hmc","MC energy",fEbins,logemin,logemax,fNTbins,fTbin);

  // reset the coefficients
  if(fCoeff)
    delete fCoeff;
  fCoeff = new TH2F("fcoeff","Coefficients for unfolding",fEbins,logemin,logemax,fNTbins,fTbin);

  // fill the histograms of estimated and measured energy for weighted events
  for(Int_t i=0;i<nentries;i++)
    {
      fCcut->GetEntry(i);  

      // mc and estimated energy
      // OJO!! Estimated energy will be taken directly from some input container
      Float_t emc   = fMcEvt->GetEnergy();
      Float_t estim = fHillas->GetSize()/15.;
      Float_t theta = fMcEvt->GetTheta()*180./3.14159; // zenith angle (in degrees)

      // compute the bin where the weight is taken from
      Int_t wbin = Int_t((TMath::Log10(emc)-logemin)/desub);

      // fill the histograms
      if(wbin<esbins)
	{
	  hest->Fill(TMath::Log10(estim),theta,fWeight[wbin]); 
	  hmc->Fill(TMath::Log10(emc),theta,fWeight[wbin]);
	}
      else
	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skpping" << endl;
    }

  // divide the previous histos to get the coefficients for unfolding
  for(Int_t j=0;j<fNTbins;j++)
    for(Int_t i=0;i<fEbins;i++)
      {
	const Float_t num = hmc->GetBinContent(i+1,j+1);
	const Float_t den = hest->GetBinContent(i+1,j+1);
	if(den)
	  fCoeff->SetBinContent(i+1,j+1,num/den);
	else
	  {
	    cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: energy bin " << i << ", thetabin " << j << " has no survivors, setting coefficient to 0" << endl;
	    fCoeff->SetBinContent(i+1,j+1,0.);
	  }
      }


  delete hmc;
  delete hest;
}

// --------------------------------------------------------------
//
// compute the coefficients used for the (iterative) unfolding 
//
void MEffAreaAndCoeffCalc::ComputeEffectiveAreas()
{
  if(!fWeight)
    {
      cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: No weights computed! nothing done" << endl;
      return;
    }

  //OJO!! radius should be read from somewhere!
  const Double_t radius = 30000.; // generation radius (in cm)
  const Double_t Atot = 3.14159*radius*radius; //total area (in cm^2)
  const Double_t logemin = TMath::Log10(fEmin);
  const Double_t logemax = TMath::Log10(fEmax);
  const Int_t esbins = fEbins*fEsubbins; // total number of subbins

  // reset the effective areas
  if(fEffA)
    delete fEffA;
  fEffA = new TH2F("feffa","Effective area",esbins,logemin,logemax,fNTbins,fTbin);

  // fill the spectrum of the survivors
  TH2F* hpass= new TH2F("hpass","Survivors",esbins,logemin,logemax,fNTbins,fTbin);
  TH2F* horig= new TH2F("horig","Original events",esbins,logemin,logemax,fNTbins,fTbin);

  fCcut->Draw("theta:logenergy>>hpass","","goff");
  fCini->Draw("theta:logenergy>>horig","","goff");
  
  // compute the effective areas
  for(Int_t j=0;j<fNTbins;j++)
    for(Int_t i=0;i<esbins;i++)
      {
	if(horig->GetBinContent(i+1,j+1)<=0)
	  {
	    cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy bin " << i <<", theta bin " << j << " contains no events setting effective area to 0" << endl;
	    fEffA->SetBinContent(i+1,j+1,0);
	  }
	else
	  fEffA->SetBinContent(i+1,j+1,Atot*hpass->GetBinContent(i+1,j+1)/horig->GetBinContent(i+1,j+1));
      }
  

  
  delete hpass;
  delete horig;
}

// --------------------------------------------------------------
//
// Call the internal functions to compute all the factors
//
void MEffAreaAndCoeffCalc::ComputeAllFactors()
{
  ComputeWeights();
  ComputeCoefficients();
  ComputeEffectiveAreas();
}
