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

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

//////////////////////////////////////////////////////////////////////////////
//
//  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:
//       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))
//
//
//  Input Containers:
//    MMcEvt
//    MMcCorsikaRunHeader
//
//  Output Container:
//    MWeight [MParameterD]
//
//////////////////////////////////////////////////////////////////////////////
#include "MMcSpectrumWeight.h"

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

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

#include "MParList.h"
#include "MParameters.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  = -1;
    fOldSlope  = -1;

    fEnergyMin = -1;
    fEnergyMax = -2;

    fNorm      =  1;

    fAllowChange = kFALSE;

    fFunc      = NULL;
    fMcEvt     = NULL;
    fWeight    = NULL;
    fZdWeights = 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;
}

// ---------------------------------------------------------------------------
//
// 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 (!fZdWeights)
        return kTRUE;

    fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
    if (!fPointing)
    {
        *fLog << err << "MPointingPos 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
{
    return fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula;
}

// ---------------------------------------------------------------------------
//
// Return the formula to calculate weights.
// Is is compiled by
//   o = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX());
//   n = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX());
//
//   result: fNorm*o/n*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 = GetSpecOldIntegral();
    const Double_t inew = GetSpecNewIntegral();

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

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

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

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

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

    TString form(GetFormulaWeightsX());

    if (fFunc)
        delete fFunc;

    fFunc = new TF1("", form);
    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 << " User normalization:       " << fNorm << 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 (fZdWeights)
    {
        const Int_t i = fZdWeights->GetXaxis()->FindFixBin(fPointing->GetZd());
        w = fZdWeights->GetBinContent(i);
    }

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read the setup from a TEnv, eg:
//   MMcSpectrumWeight.NewSlope: -2.6
//   MMcSpectrumWeight.Norm:      1.0
//   MMcSpectrumWeight.Formula:  pow(MMcEvt.fEnergy, -2.6)
//
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, "Formula", print))
    {
        rc = kTRUE;
        SetFormula(GetEnvValue(env, prefix, "Formula", fFormula));
    }

    return rc;
}