/* ======================================================================== *\
!
! *
! * 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): Marcos Lopez  10/2003 <mailto:marcos@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MWeightEnergySlopeCalc
//               
//  Change the spectrum of the MC showers simulated with Corsika (a power law)
//  to a new one, which can be either, again a power law but with a different
//  spectral index, or a generalizeed spectra (given as a TF1 object)
//
//  Method:
//  ------ 
//                                 
//  -Corsika spectrun: dN/dE = A * E^(-a)
//    with a = fCorsikaSlope, and A = N/integral{E*de} from ELowLim to EUppLim
//
//  -New spectrum:     dN/dE = B * g(E)
//    where B = N/integral{g*dE} from ELowLim to EUppLim, and N=NumEvents 
//
//  For converting the spectrum simulated with Corsika to the new one, we apply
//  a weight to each event, given by:
//
//     W(E) = B/A * g(E)/E^(-a) 
//
//  (The factor B/A is introduced in order both the original and new spectrum 
//  have the same area (i.e. number of showers))
//
//  Note:
//  ------
//   -In the calculations, the new spectrum is treated internally as a TF1 
//    object, to give to the user the option of applying a general spectrum. 
//    Is the new spectrum is just a power law (i.e. the user only specify the 
//    slope),it is converted to a TF1 object for consistency. 
//   -With the original corsika spectrum, as it is always a power law, not TF1 
//    object is used to represent the spectrum.
//
// 
//  Input Containers:                         
//   MMcEvt, MMcRunHeader, MMcCorsikaRunHeader
//                                            
//  Output Container:                        
//   MWeight  
//
//////////////////////////////////////////////////////////////////////////////

#include "MMcWeightEnergySpecCalc.h"

#include "MParList.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MMcEvt.hxx"
#include "MMcRunHeader.hxx"
#include "MMcCorsikaRunHeader.h"
#include "MWeight.h"

#include "TF1.h"
#include "TGraph.h"

ClassImp(MMcWeightEnergySpecCalc);

using namespace std;



// ---------------------------------------------------------------------------
//
//
MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Float_t slope, 
					  const char *name, const char *title)
{
    fNewSpecIsPowLaw = kTRUE;

    fNewSlope = slope; 

    fNewSpectrum = new TF1("NewSpectrum","pow(x,[0])");
    fNewSpectrum->SetParameter(0,fNewSlope);


    fName  = name  ? name  : "MMcWeightEnergySpecCalc";
    fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 

    AddToBranchList("MMcEvt.fEnergy");

    fAllEvtsTriggered         =  kFALSE;
    fTotalNumSimulatedShowers =  0;
} 

MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(TF1* spectrum, 
					  const char *name, const char *title)
{
    fNewSpecIsPowLaw = kFALSE;

    fNewSpectrum = spectrum;

    fName  = name  ? name  : "MMcWeightEnergySpecCalc";
    fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 

    AddToBranchList("MMcEvt.fEnergy");

    fAllEvtsTriggered         =  kFALSE;
    fTotalNumSimulatedShowers =  0;
} 


MMcWeightEnergySpecCalc::~MMcWeightEnergySpecCalc()
{
  delete fNewSpectrum; 
}


// ---------------------------------------------------------------------------
//
//
Int_t MMcWeightEnergySpecCalc::PreProcess (MParList *pList)
{

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

    fWeight = (MWeight*)pList->FindObject("MWeight");
    if (!fWeight)
    {
      fWeight = (MWeight*)pList->FindCreateObj("MWeight");
        if (!fWeight)
            return kFALSE;
    }

    return kTRUE;
}


// ----------------------------------------------------------------------------
//
// Executed each time a new root file is loaded
// We will need the fCorsikaSlope, and fE{Upp,Low}Lim to calculate the weights
//
Bool_t MMcWeightEnergySpecCalc::ReInit(MParList *plist)
{
    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
    if (!runheader)
    {
      *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
      return kFALSE;
    }

    MMcCorsikaRunHeader *corrunheader  = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
    if (!corrunheader)
      {
	*fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl;
	return kFALSE;
      }


    fCorsikaSlope = corrunheader->GetSlopeSpec();
    fELowLim      = corrunheader->GetELowLim();
    fEUppLim      = corrunheader->GetEUppLim();

    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();    
    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();    



    *fLog << inf << "Slope of primaries' energy spectrum of Simulated showers: " << fCorsikaSlope << endl;
    *fLog << inf << "Limits of energy range of Simulated showers: " 
	  << fELowLim <<" - " << fEUppLim << endl;
    *fLog << inf << "New Slope for Simulated showers: " << fNewSlope << endl;
    *fLog << inf << "Total Number of Simulated showers: " 
	  << fTotalNumSimulatedShowers << endl;    
    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;

    

    //
    // Starting from fCorsikaSlope, and fE{Upp,Low}Lim, calculate the integrals
    // of both, the original Corsika spectrum and the new one.
    //
    //
    // For the Corsika simulated spectrum (just a power law), we have:
    //
    Double_t a = fCorsikaSlope;
    fCorSpecInt = ( pow(fEUppLim,1+a) - pow(fELowLim,1+a) ) / ( 1+a );
 
    //
    // For the new spectrum:
    //
    fNewSpectrum->SetRange(fELowLim, fEUppLim);
    
    if (fNewSpecIsPowLaw)     // just the integral of a power law
      {
	Double_t b = fNewSlope;
	fNewSpecInt = ( pow(fEUppLim,1+b) - pow(fELowLim,1+b) ) / ( 1+b );
      }
    else
      { // In this case we have to integrate the new spectrum numerically. We 
	// could do simply fNewSpectrum->Integral(fELowLim,fEUppLim), but ROOT
	// fails integrating up to fEUppLim for a sharp cutoff spectrum
	
	//
	// Trick to calculate the integral numerically (it works better than 
	// fNewSpectrum->Integral(fELowLim,fEUppLim) (although not perfectlly)
	//
	fNewSpectrum->SetNpx(1e4);
	TGraph gr(fNewSpectrum,"i");
	Int_t Npx = gr.GetN();
	Double_t* y = gr.GetY();

	Double_t integral = y[Npx-1];
	fNewSpecInt = integral; 
    }


    return kTRUE;
}


// ----------------------------------------------------------------------------
//
//
Int_t MMcWeightEnergySpecCalc::Process()
{
    const Double_t energy = fMcEvt->GetEnergy();

    Double_t weight = fCorSpecInt/fNewSpecInt * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope);
     
    fWeight->SetWeight( weight );

    return kTRUE;
}


















