source: trunk/MagicSoft/Mars/msim/MSimAbsorption.cc@ 9384

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular 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 appears 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,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MSimAbsorption
//
//  Task to calculate wavelength or incident angle dependent absorption
//
//  Input Containers:
//   MPhotonEvent
//   MCorsikaEvtHeader
//   MCorsikaRunHeader
//
//  Output Containers:
//   MPhotonEvent
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimAbsorption.h"

#include <fstream>

#include <TRandom.h>
#include <TSpline.h>

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

#include "MParList.h"

#include "MBinning.h"
#include "MArrayD.h"

#include "MSpline3.h"

#include "MCorsikaEvtHeader.h"
#include "MCorsikaRunHeader.h"
#include "MPhotonEvent.h"
#include "MPhotonData.h"

ClassImp(MSimAbsorption);

using namespace std;

// --------------------------------------------------------------------------
//
//  Default Constructor.
//
MSimAbsorption::MSimAbsorption(const char* name, const char *title)
    : fEvt(0), fHeader(0), fSpline(0), fUseTheta(kFALSE)
{
    fName  = name  ? name  : "MSimAbsorption";
    fTitle = title ? title : "Task to calculate wavelength dependent absorption";
}

// --------------------------------------------------------------------------
//
//  Calls Clear()
//
MSimAbsorption::~MSimAbsorption()
{
    Clear();
}

// --------------------------------------------------------------------------
//
// Delete fSpline and set to 0.
//
void MSimAbsorption::Clear(Option_t *)
{
    if (fSpline)
        delete fSpline;
    fSpline=0;
}

// --------------------------------------------------------------------------
//
// Read a TGraph from a file and return a newly allocated MSpline3.
//
MSpline3 *MSimAbsorption::ReadSpline(const char *fname)
{
    *fLog << inf << "Reading spline from " << fname << endl;

    TGraph g(fname);
    if (g.GetN()==0)
    {
        *fLog << err << "ERROR - No data points from " << fname << "." << endl;
        return 0;
    }

    // option: b1/e1 b2/e2   (first second derivative?)
    // option: valbeg/valend (first second derivative?)

    return new MSpline3(g);
}

// --------------------------------------------------------------------------
//
// Initializes a spline from min to max with n points with 1.
//
void MSimAbsorption::InitUnity(UInt_t n, Float_t min, Float_t max)
{
    // Delete the existing spline
    Clear();

    // We need at least two points (the edges)
    if (n<2)
        return;

    // Initialize an array with the x-values
    const MBinning bins(n-1, min, max);

    // Initialize an array with all one
    MArrayD y(n);
    y.Reset(1);

    // Set the spline
    fSpline = new MSpline3(bins.GetEdges(), y.GetArray(), n);
}

// --------------------------------------------------------------------------
//
// Takes the existing fSpline and multiplies it with the given spline.
// As x-points the values from fSpline are used.
//
void MSimAbsorption::Multiply(const TSpline3 &spline)
{
    if (!fSpline)
    {
        Clear();
        // WARNING WARNING WARNING: This is a very dangerous cast!
        fSpline = (MSpline3*)spline.Clone();
        return;
    }

    // Initialize a TGraph with the number of knots from a TSpline
    TGraph g(fSpline->GetNp());

    // Loop over all knots
    for (int i=0; i<fSpline->GetNp(); i++)
    {
        // Get th i-th knot
        Double_t x, y;
        fSpline->GetKnot(i, x, y);

        // Multiply y by the value from the given spline
        y *= spline.Eval(x);

        // push the point "back"
        g.SetPoint(i, x, y);
    }

    // Clear the spline and initialize a new one from the updated TGraph
    Clear();
    fSpline = new MSpline3(g);
}

// --------------------------------------------------------------------------
//
// Initializes a TSpline3 with n knots x and y. Call Multiply for it.
//
void MSimAbsorption::Multiply(UInt_t n, const Double_t *x, const Double_t *y)
{
    const TSpline3 spline("Spline", const_cast<Double_t*>(x), const_cast<Double_t*>(y), n);
    Multiply(spline);
}

// --------------------------------------------------------------------------
//
// Read a Spline from a file using ReadSpline.
// Call Multiply for it.
//
void MSimAbsorption::Multiply(const char *fname)
{
    const MSpline3 *spline = ReadSpline(fname);
    if (!spline)
        return;

    fFileName = "";

    Multiply(*spline);

    delete spline;
}

// --------------------------------------------------------------------------
//
// Read a spline from a file and set fSpline accfordingly.
//
Bool_t MSimAbsorption::ReadFile(const char *fname)
{
    if (fname)
        fFileName = fname;

    MSpline3 *spline = ReadSpline(fFileName);
    if (!spline)
        return kFALSE;


    // option: b1/e1 b2/e2   (first second derivative?)
    // option: valbeg/valend (first second derivative?)

    Clear();
    fSpline = spline;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Search for the needed parameter containers. Read spline from file
// calling ReadFile();
//
Int_t MSimAbsorption::PreProcess(MParList *pList)
{
    if (fFileName.IsNull())
    {
        *fLog << inf << "No file given... skipping." << endl;
        return kSKIP;
    }

    fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
    if (!fEvt)
    {
        *fLog << err << "MPhotonEvent not found... aborting." << endl;
        return kFALSE;
    }

    fHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader");
    if (!fHeader)
    {
        *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl;
        return kFALSE;
    }

    *fLog << inf << "Using " << (fUseTheta?"Theta":"Wavelength") << " for absorption." << endl;

    return ReadFile();
}

// --------------------------------------------------------------------------
//
Bool_t MSimAbsorption::ReInit(MParList *pList)
{
    if (fUseTheta)
        return kTRUE;

    MCorsikaRunHeader *h = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
    if (!h)
    {
        *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
        return kFALSE;
    }

    if (h->GetWavelengthMin()<fSpline->GetXmin())
        *fLog << warn << "WARNING - Lower bound of wavelength bandwidth exceeds lower bound of spline." << endl;

    if (h->GetWavelengthMax()>fSpline->GetXmax())
        *fLog << warn << "WARNING - Upper bound of wavelength bandwidth exceeds upper bound of spline." << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Throw all events out of the MPhotonEvent which don't survive.
// Depending on fUseTheta either the wavelength or the incident angle
// is used.
//
Int_t MSimAbsorption::Process()
{
    // Get the number of photons in the list
    const Int_t num = fEvt->GetNumPhotons();

    // Counter for number of total and final events
    Int_t cnt = 0;
    for (Int_t i=0; i<num; i++)
    {
        // Get i-th photon from the list
        const MPhotonData &ph = (*fEvt)[i];

        // Depending on fUseTheta get the incident angle of the wavelength
        const Double_t wl = fUseTheta ? ph.GetTheta()*TMath::RadToDeg() : ph.GetWavelength();

        // Get the efficiency (the probability that this photon will survive)
        // from the spline.
        const Double_t eff = fSpline->Eval(wl);

        // Get a random value between 0 and 1 to determine whether the photn will survive
        // gRandom->Rndm() = [0;1[
        if (gRandom->Rndm()>=eff)
            continue;

        // Copy the surviving events bakc in the list
        (*fEvt)[cnt++] = ph;
    }

    // Now we shrink the array to the number of new entries.
    fEvt->Shrink(cnt);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// FileName: reflectivity.txt
// UseTheta: No
//
Int_t MSimAbsorption::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "FileName", print))
    {
        rc = kTRUE;
        SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
    }

    if (IsEnvDefined(env, prefix, "UseTheta", print))
    {
        rc = kTRUE;
        SetUseTheta(GetEnvValue(env, prefix, "UseTheta", fUseTheta));
    }

    return rc;
}