source: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc@ 8045

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz 5/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Marcos Lopez 10/2003 <mailto:marcos@gae.ucm.es>
!
!   Copyright: MAGIC Software Development, 2000-2006
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MMcSpectrumWeight
//               
//  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:
//       The new spectrum is passed as a char* (SetFormula())
//
//  Method:
//  -------
//                                 
//   - Corsika spectrun: dN/dE = A * E^(a)
//     with a = fOldSlope, 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))
//
//
//  If using SetFormula you can specify formulas accepted by TF1, eg:
//      pow(X, -2.6)
//  (Rem: all capital (!) 'X' are replaced by the corresponding %s.fEnergy
//        automatically)
//
//  For more details of the setup see MMcSpectrumWeight::ReadEnv
//
//
//  Input Containers:
//    MMcEvt
//    MMcCorsikaRunHeader
//
//  Output Container:
//    MWeight [MParameterD]
//
//////////////////////////////////////////////////////////////////////////////
#include "MMcSpectrumWeight.h"

#include <TF1.h>
#include <TH1.h>
#include <TSpline.h>

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

#include "MParList.h"
#include "MParameters.h"

#include "MHillas.h"
#include "MPointingPos.h"

#include "MMcEvt.hxx"
#include "MMcCorsikaRunHeader.h"

ClassImp(MMcSpectrumWeight);

using namespace std;

void MMcSpectrumWeight::Init(const char *name, const char *title)
{
    fName  = name  ? name  : "MMcSpectrumWeight";
    fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; 

    AddToBranchList("MMcEvt.fEnergy");

    fNameWeight  = "MWeight";
    fNameMcEvt   = "MMcEvt";

    fNewSlope    = -9;
    fOldSlope    = -9;

    fEnergyMin   = -1;
    fEnergyMax   = -2;

    fNorm        =  1;
    fNormEnergy  = -1;

    fAllowChange = kFALSE;

    fFunc        = NULL;
    fMcEvt       = NULL;
    fHillas      = NULL;
    fWeight      = NULL;
    fWeightsZd   = NULL;
    fWeightsSize = NULL;
    fPointing    = NULL;
}

// ---------------------------------------------------------------------------
//
// Default Constructor.
//
MMcSpectrumWeight::MMcSpectrumWeight(const char *name, const char *title)
{
    Init(name,title);
} 

// ---------------------------------------------------------------------------
//
// Destructor. If necessary delete fFunc
//
MMcSpectrumWeight::~MMcSpectrumWeight()
{
    if (fFunc)
        delete fFunc;
//    if (fWeightsSize)
//        delete fWeightsSize;
}

// ---------------------------------------------------------------------------
//
// Search for
//   - fNameMcEvt [MMcEvtBasic]
//
// Find/Create:
//   - fNameWeight [MWeight]
//
Int_t MMcSpectrumWeight::PreProcess(MParList *pList)
{
    fMcEvt = (MMcEvt*)pList->FindObject(fNameMcEvt, "MMcEvtBasic");
    if (!fMcEvt)
    {
        *fLog << err << fNameMcEvt << " [MMcEvtBasic] not found... abort." << endl;
        return kFALSE;
    }

    fWeight = (MParameterD*)pList->FindCreateObj("MParameterD", fNameWeight);
    if (!fWeight)
        return kFALSE;

    if (fWeightsZd)
    {
        fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
        if (!fPointing)
        {
            *fLog << err << "MPointingPos not found... abort." << endl;
            return kFALSE;
        }
    }

    if (fWeightsSize)
    {
        fHillas = (MHillas*)pList->FindObject("MHillas");
        if (!fHillas)
        {
            *fLog << err << "MHillas not found... abort." << endl;
            return kFALSE;
        }
    }

    return kTRUE;
}

// ---------------------------------------------------------------------------
//
// Replace {fNameMcEvt}.fEnergy by "(x)" and return the result.
//
TString MMcSpectrumWeight::ReplaceX(TString str) const
{
    return str.ReplaceAll(Form("%s.fEnergy", fNameMcEvt.Data()), "(x)");
}

// ---------------------------------------------------------------------------
//
// Return the function corresponding to the mc spectrum with
// slope fOldSlope: pow({fNameMcEvt}.fEnergy, fOldSlope)
//
// The slope is returned as %.3f
//
TString MMcSpectrumWeight::GetFormulaSpecOld() const
{
    return Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fOldSlope);
}

// ---------------------------------------------------------------------------
//
// Return the function corresponding to the new spectrum with
// slope fNewSlope: pow({fNameMcEvt}.fEnergy, fNewSlope)
//
// The slope is returned as %.3f
//
// If a different formula is set (SetFormula()) this formula is returned
// unchanged.
//
TString MMcSpectrumWeight::GetFormulaSpecNew() const
{
    TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula.Data();
    if (!fFormula.IsNull())
        str.ReplaceAll("X", Form("(%s.fEnergy)", fNameMcEvt.Data()));

    return str;
}

// ---------------------------------------------------------------------------
//
// Return the formula to calculate weights.
// Is is compiled by
//   o1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX());
//   n1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX());
//   o2 = CalcSpecOld(fNormEnergy);
//   n2 = CalcSpecNew(fNormEnergy);
//
//   result (fNormEnergy<0):
//      fNorm*o1/n1*GetFormulaNewSpec()/GetFormulaOldSpec()
//
//   result (fNormEnergy>=0):
//      fNorm*o2/n2*GetFormulaNewSpec()/GetFormulaOldSpec()
//
// fNorm is 1 by default but can be overwritten using SetNorm()
//
// If the formulas GetFormulaSpecOldX() and GetFormulaSpecNewX()
// are equal only fNorm is returned.
//
// The normalization constant is returned as %.16f
//
// Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600))
//
TString MMcSpectrumWeight::GetFormulaWeights() const
{
    if (GetFormulaSpecOld()==GetFormulaSpecNew())
        return Form("%.16f", fNorm);

    const Double_t iold = fNormEnergy<0 ? GetSpecOldIntegral() : CalcSpecOld(fNormEnergy);
    const Double_t inew = fNormEnergy<0 ? GetSpecNewIntegral() : CalcSpecNew(fNormEnergy);

    const Double_t norm = fNorm*iold/inew;

    return Form("%.16f*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data());
}

// ---------------------------------------------------------------------------
//
// Returns the integral between fEnergyMin and fEnergyMax of
// GetFormulaSpecNewX() describing the destination spectrum
//
Double_t MMcSpectrumWeight::GetSpecNewIntegral() const
{
    TF1 funcnew("Dummy", GetFormulaSpecNewX());
    return funcnew.Integral(fEnergyMin, fEnergyMax);
}

// ---------------------------------------------------------------------------
//
// Returns the integral between fEnergyMin and fEnergyMax of
// GetFormulaSpecOldX() describing the simulated spectrum
//
Double_t MMcSpectrumWeight::GetSpecOldIntegral() const
{
    TF1 funcold("Dummy", GetFormulaSpecOldX());
    return funcold.Integral(fEnergyMin, fEnergyMax);
}

// ---------------------------------------------------------------------------
//
// Returns the value of GetFormulaSpecNewX() at the energy e describing
// the destination spectrum
//
Double_t MMcSpectrumWeight::CalcSpecNew(Double_t e) const
{
    TF1 funcnew("Dummy", GetFormulaSpecNewX());
    return funcnew.Eval(e);
}

// ---------------------------------------------------------------------------
//
// Returns the value of GetFormulaSpecOldX() at the energy e describing
// the simulated spectrum
//
Double_t MMcSpectrumWeight::CalcSpecOld(Double_t e) const
{
    TF1 funcnew("Dummy", GetFormulaSpecOldX());
    return funcnew.Eval(e);
}

void MMcSpectrumWeight::SetWeightsSize(TH1D *h)
{
    fWeightsSize=h;
    /*
    if (h==0)
    {
        fWeightsSize=0;
        return;
    }

    if (fWeightsSize)
        delete fWeightsSize;

    const Double_t xmin = TMath::Log10(h->GetXaxis()->GetXmin());
    const Double_t xmax = TMath::Log10(h->GetXaxis()->GetXmax());
    const Double_t xnum = h->GetNbinsX()+1;

    fWeightsSize = new TSpline3("WeightsSize", xmin, xmax,
                                h->GetArray()+1, xnum);*/
}

// ---------------------------------------------------------------------------
//
// Initialize fEnergyMin, fEnergymax and fOldSlope from MMcCorsikaRunHeader
// by GetELowLim(), GetEUppLim() and GetSlopeSpec().
//
// If fEnergyMax>fEnergyMin (means: the values have already been
// initialized) and !fAllowChange the consistency of the new values
// with the present values is checked with a numerical precision of 1e-10.
// If one doesn't match kFALSE is returned.
//
// If the mc slope is -1 kFALSE is returned.
//
// If the new slope for the spectrum is -1 it is set to the original MC
// slope.
//
// fFunc is set to the formula returned by GetFormulaWeightsX()
//
Bool_t MMcSpectrumWeight::Set(const MMcCorsikaRunHeader &rh)
{
    if (fEnergyMax>fEnergyMin && !fAllowChange)
    {
        if (TMath::Abs(fOldSlope-rh.GetSlopeSpec())>1e-10)
        {
            *fLog << err;
            *fLog << "ERROR - The slope of the Monte Carlo is not allowed to change ";
            *fLog << "(" << fOldSlope << " --> " << rh.GetSlopeSpec() << ")... abort." << endl;
            return kFALSE;
        }
        if (TMath::Abs(fEnergyMin-rh.GetELowLim())>1e-10)
        {
            *fLog << err;
            *fLog << "ERROR - The minimum simulated Monte Carlo energy is not allowed to change ";
            *fLog << "(" << fEnergyMin << " --> " << rh.GetELowLim() << ")... abort." << endl;
            return kFALSE;
        }
        if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10)
        {
            *fLog << err;
            *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change (";
            *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl;
            return kFALSE;
        }
        return kTRUE;
    }

    fOldSlope  = rh.GetSlopeSpec();
    fEnergyMin = rh.GetELowLim();
    fEnergyMax = rh.GetEUppLim();

    if (fNewSlope==-9)
    {
        *fLog << inf << "The new slope of the power law is undefined (-9)... no weighting applied." << endl;
        fNewSlope = fOldSlope;
    }

    if (fFunc)
        delete fFunc;

    fFunc = new TF1("", GetFormulaWeightsX());
    gROOT->GetListOfFunctions()->Remove(fFunc);
    fFunc->SetName("SpectralWeighs");

    return kTRUE;
}

// ---------------------------------------------------------------------------
//
// The current contants are printed
//
void MMcSpectrumWeight::Print(Option_t *o) const
{
    *fLog << all << GetDescriptor() << endl;
    *fLog << " Simulated energy range:   " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
    *fLog << " Simulated spectral slope: " << fOldSlope << endl;
    *fLog << " New spectral slope:       " << fNewSlope << endl;
    *fLog << " Additional user norm.:    " << fNorm << endl;
    *fLog << " Spectra are normalized:   " << (fNormEnergy<0?"by integral":Form("at %.1fGeV", fNormEnergy)) << endl;
    *fLog << " Old Spectrum:     " << GetFormulaSpecOldX() << "   (I=" << GetSpecOldIntegral() << ")" << endl;
    *fLog << " New Spectrum:     " << GetFormulaSpecNewX() << "   (I=" << GetSpecNewIntegral() << ")" << endl;
    if (fFunc)
        *fLog << " Weight function:  " << fFunc->GetTitle()   << endl;
}

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

    return Set(*rh);
}

// ----------------------------------------------------------------------------
//
// Fill the result of the evaluation of fFunc at fEvEvt->GetEnergy
// into the weights container.
//
Int_t MMcSpectrumWeight::Process()
{
    const Double_t e = fMcEvt->GetEnergy();

    Double_t w = 1;

    if (fWeightsZd)
    {
        const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
        w = fWeightsZd->GetBinContent(i);
    }
    if (fWeightsSize)
    {
        const Int_t i = fWeightsSize->GetXaxis()->FindFixBin(fHillas->GetSize());
        w *= fWeightsSize->GetBinContent(i);
       // w *= fWeightsSize->Eval(TMath::Log10(fHillas->GetSize()));
    }

    fWeight->SetVal(fFunc->Eval(e)*w);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//
//   MMcSpectrumWeight.NewSlope:   -2.6
//    The new slope of the spectrum
//
//   MMcSpectrumWeight.Norm:        1.0
//    An additional artificial scale factor
//
//   MMcSpectrumWeight.NormEnergy:  200
//    To normalize at a given energy instead of the integral
//
//   MMcSpectrumWeight.Formula:     pow(X, -2.6)
//    A formula to which the spectrum is weighted (use a capital X for
//    the energy)
//
Int_t MMcSpectrumWeight::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "NewSlope", print))
    {
        rc = kTRUE;
        SetNewSlope(GetEnvValue(env, prefix, "NewSlope", fNewSlope));
    }
    if (IsEnvDefined(env, prefix, "Norm", print))
    {
        rc = kTRUE;
        SetNorm(GetEnvValue(env, prefix, "Norm", fNorm));
    }
    if (IsEnvDefined(env, prefix, "NormEnergy", print))
    {
        rc = kTRUE;
        SetNormEnergy(GetEnvValue(env, prefix, "NormEnergy", fNormEnergy));
    }
    if (IsEnvDefined(env, prefix, "Formula", print))
    {
        rc = kTRUE;
        SetFormula(GetEnvValue(env, prefix, "Formula", fFormula));
    }

    return rc;
}