/* ======================================================================== *\
!
! *
! * 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 "TFile.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), fEmin(10.), fEmax(10000.), fEbins(10), fEsubbins(20),
    fCoeff(NULL), fEffA(NULL), fFile(NULL)
{
  // set the default function
  SetFunction("4.e9*pow(x,-2.6+1)");  

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

  // define some useful aliases
  fCini->SetAlias("logenergy","log10(MMcEvtBasic.fEnergy)");
  fCini->SetAlias("theta","MMcEvtBasic.fTelescopeTheta*180./3.14159");

  fCcut->SetAlias("logenergy","log10(MMcEvt.fEnergy)");
  fCcut->SetAlias("theta","MMcEvt.fTelescopeTheta*180./3.14159");

  fCcut->SetBranchAddress("MHillas.",&fHillas);
  fCcut->SetBranchAddress("MMcEvt.",&fMcEvt);
  
  // borra
  fFile = new TFile("test.root","RECREATE");
}

// -------------------------------------------------------------------------
//
// 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;
    }
  else
    {
      fEmin = emin;
      fEmax = emax;
    }
  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");
  // borra
  fHorig->Write();
}

// -------------------------------------------------------------------------
//
// 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();
  
  // borra
  TH1F hw("hw","hw",esbins,fEmin,fEmax);
  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;
	}
      // borra
      hw.SetBinContent(i+1,fWeight[i]);
    }
  // borra
  hw.Write();
}

// --------------------------------------------------------------
//
// 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);

  // borra
  TH1F* hest1  = new TH1F("hest1","Estimated energy",fEbins,logemin,logemax);
  TH1F* hmc1   = new TH1F("hmc1","MC energy",fEbins,logemin,logemax);
  TH1F* coef1  = new TH1F("coef1","coefficients",fEbins,logemin,logemax);

  // 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()/0.18/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]);
	  // borra
	  if(theta<fTbin[1])
	    {
	      hest1->Fill(TMath::Log10(estim),fWeight[wbin]); 
	      hmc1->Fill(TMath::Log10(emc),fWeight[wbin]);
	    }
	}
      else
	cout << "MEffAreaAndCoeffCalc::ComputeCoefficients Warning: event " << i << " with energy " << emc << " has energy out of bounds, skipping" << 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);
	//borra
	const Float_t num1 = hmc1->GetBinContent(i+1);
	const Float_t den1 = hest1->GetBinContent(i+1);
	if(den)
	  {
	    fCoeff->SetBinContent(i+1,j+1,num/den);
	    //borra
	    if(j==0 && den1)
	      coef1->SetBinContent(i+1,num1/den1);
	  }
	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.);
	  }
      }


  //borra
  hmc1->Write();
  hest1->Write();
  coef1->Write();

  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 Atot = 1.; //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
  const Double_t de = (logemax-logemin)/fEbins; // bin size (in log)
  const Double_t desub = de/fEsubbins; // subbin size (in log)

  // reset the effective areas
  if(fEffA)
    delete fEffA;
  fEffA = new TH2F("feffa","Effective area",fEbins,logemin,logemax,fNTbins,fTbin);
  //borra
  TH1F eff("eff","Effective area",fEbins,logemin,logemax);
  
  // 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<fEbins;i++)
      {
	Float_t effarea =0;
	Float_t wtot = 0;
	for(Int_t k=0;k<fEsubbins;k++)
	  {
	    Float_t numerator = hpass->GetBinContent(i*fEsubbins+k+1,j+1);
	    Float_t denominator = horig->GetBinContent(i*fEsubbins+k+1,j+1);
	    
	    if(denominator<=0)
	      cout << "MEffAreaAndCoeffCalc::ComputeEffectiveAreas Warning: energy subbin " << i*fEsubbins+k <<", theta bin " << j << " contains no events" << endl;
	    else
	      {
		const Double_t ew = TMath::Power(10,logemin+(i*fEsubbins+k+0.5)*desub);
		const Double_t ww = fSpec->Eval(ew);
		effarea+=Atot*numerator/denominator*ww;
		wtot   += ww;
	      }
	  }
	if(!wtot)
	  {
	    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,effarea/wtot);	
	    if(j==0)
	      {
		cout << "*****" << i << ": " << effarea/wtot << endl;
		eff.SetBinContent(i+1,effarea/wtot);
	      }
	  }
      }

  // borra
  eff.Write();
  
  delete hpass;
  delete horig;
}

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