source: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc@ 5653

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MMcWeightEnergySlopeCalc
//               
//  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 spectrum. The new spectrum can be 
//  pass to this class in different ways:
//    1. Is the new spectrum will be a power law, just introduce the slope
//       of this power law.
//    2. Is the new spectrum will have a general shape, different options are
//       available:
//       a) The new spectrum is pass as a TF1 function 
//       b) The new spectrum is pass as a char*
//       c) The new spectrum is pass as a "interpreted function", i.e., a 
//          function defined inside a ROOT macro, which will be invoked by the
//          ROOT Cint itself. This is the case when we use ROOT macros.
//       d) The new spectrum is pass as a "real function", i.e., a 
//          function defined inside normal c++ file.
//
//  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) 
//
//  In the case the new spectrum is simply a power law: dN/dE = B * E^(b), we
//  have: 
//
//     W(E) = B/A * E^(b-a) 
//
//  (The factor B/A is used in order both the original and new spectrum have 
//   the same area (i.e. in order they represent the same number of showers))
//
//  Note:
//  ------
//   -If the the new spectrum is just a power law (i.e. the user only specify
//     the slope), the needed calculations (such as the integral of the 
//     spectrum) are done analytically. But if the new spectrum is given as a 
//     TF1 object, the calculations is done numerically.
//
//  ToDo:
//  ----- 
//   -Give to the user also the possibility to specify the integral of the 
//    spectrum as another TF1 object (or C++ function)
//
//
//  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;



void MMcWeightEnergySpecCalc::Init(const char *name, const char *title)
{

    fName  = name  ? name  : "MMcWeightEnergySpecCalc";
    fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 
    
    AddToBranchList("MMcEvt.fEnergy");

    fAllEvtsTriggered         =  kFALSE;
    fTotalNumSimulatedShowers =  0;
}


// ---------------------------------------------------------------------------
//
// Constructor. The new spectrum will be just a power law.
//
MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Float_t slope, 
                                                 const char *name, const char *title)
{
    fNewSpecIsPowLaw = kTRUE;
    fNewSlope = slope; 

    fNewSpectrum = NULL; 

    Init(name,title);
} 

// ---------------------------------------------------------------------------
//
// Constructor. The new spectrum will have a general shape, given by the user 
// as a TF1 function.
//
MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(const TF1& spectrum, 
                                                 const char *name, const char *title)
{
    fNewSpecIsPowLaw = kFALSE;
    fNewSpectrum = (TF1*)spectrum.Clone();

    Init(name,title);
} 

// ---------------------------------------------------------------------------
//
// As before, but the function which represent the new spectrum is given as
// a char* . Starting from it, we build a TF1 function
//
MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(const char* spectrum, 
                                                 const char *name, const char *title)
{
    fNewSpecIsPowLaw = kFALSE;
    fNewSpectrum = new TF1("NewSpectrum",spectrum);

    Init(name,title);
} 

// ---------------------------------------------------------------------------
//
// As before, but the new spectrum is given as a intrepreted C++ function. 
// Starting from it we build a TF1 function.
// This constructor is called for interpreted functions by CINT, i.e., when
// the functions are declared inside a ROOT macro.
//
// NOTE: you muss do a casting to (void*) of the function that you pass to this
//       constructor before invoking it in a macro, e.g.
//
//       Double_t myfunction(Double_t *x, Double_t *par)
//         {
//           ...
//         }
//
//        MMcWeightEnergySpecCalc wcalc((void*)myfunction); 
//
//        tasklist.AddToList(&wcalc);
//
//        otherwise ROOT will invoke the constructor McWeightEnergySpecCalc(
//         const char* spectrum, const char *name, const char *title)
//
MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(void* function, 
                                                 const char *name, const char *title)
{
    fNewSpecIsPowLaw = kFALSE;
    fNewSpectrum = new TF1("NewSpectrum",function,0,1,1);

    Init(name,title);
} 

// ---------------------------------------------------------------------------
//
// As before, but this is the constructor for real functions, i.e. it is called
// when invoked with the normal C++ compiler, i.e. not inside a ROOT macro.
//
MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Double_t (*function)(Double_t*x, Double_t* par), const Int_t npar,  const char *name, const char *title)
{
    fNewSpecIsPowLaw = kFALSE;
    fNewSpectrum = new TF1("NewSpectrum",function,0,1,1);

    Init(name,title);
} 



// ----------------------------------------------------------------------------
//
// Destructor. Deletes the cloned fNewSpectrum.
//
MMcWeightEnergySpecCalc::~MMcWeightEnergySpecCalc()
{
  if (fNewSpectrum)
    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->FindCreateObj("MWeight");
    if (!fWeight)
      {
        *fLog << err << dbginf << "MWeight not found... exit." << endl;
        return kFALSE;
      }
    
    return kTRUE;
}



// ----------------------------------------------------------------------------
//
// Executed each time a new root file is loaded
// We will need 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 = (Double_t)corrunheader->GetSlopeSpec();
    fELowLim      = (Double_t)corrunheader->GetELowLim();
    fEUppLim      = (Double_t)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;



    //
    // Sanity checks to be sure that we won't divide by zero later on
    //
    if(fCorsikaSlope == -1. || fNewSlope == -1.) 
      {
        *fLog << err << "The Slope of the power law must be different of -1... exit" << endl;
        return kFALSE;
      }


    //
    // 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:
    //
    fCorSpecInt = ( pow(fEUppLim,1+fCorsikaSlope) - pow(fELowLim,1+fCorsikaSlope) ) / ( 1+fCorsikaSlope );
 
 
    //
    // For the new spectrum:
    //
    if (fNewSpecIsPowLaw)     // just the integral of a power law
      {
        fNewSpecInt = ( pow(fEUppLim,1+fNewSlope) - pow(fELowLim,1+fNewSlope) )/ ( 1+fNewSlope );
      }
    else
      { 
        fNewSpectrum->SetRange(fELowLim, fEUppLim);

        // 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(1000);
        TGraph gr(fNewSpectrum,"i");
        Int_t Npx = gr.GetN();
        Double_t* y = gr.GetY();

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


    return kTRUE;
}



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

    const Double_t C = fCorSpecInt / fNewSpecInt;
    Double_t weight;


    if (fNewSpecIsPowLaw)   
      weight = C * pow(energy,fNewSlope-fCorsikaSlope);
    else
      weight = C * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope);
     

    fWeight->SetWeight( weight );

    return kTRUE;
}