source: trunk/MagicSoft/Mars/mmontecarlo/MMcThresholdCalc.cc@ 855

/* ======================================================================== *\
!
! *
! * 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): Javier Lopez 05/2001 (jlopez@ifae.es)
!              Thomas Bretz 06/2001 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MMcThresholdCalc                                                       //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "MMcThresholdCalc.h"

#include <math.h>

#include "MParList.h"

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

#include "MMcEvt.hxx"
#include "MMcTrig.hxx"

#include "MHMcEnergy.h"

ClassImp(MMcThresholdCalc)

MMcThresholdCalc::MMcThresholdCalc(const UInt_t dim, const char* name,
                                   const char* title) : fDimension(dim)
{
    *fName  = name  ? name  : "MMcThresholdCalc";
    *fTitle = title ? title : "Task to calculate the energy threshold from Monte Carlo";

    fMcTrig     = new (MMcTrig*)[fDimension];
    fHMcEnergy  = new (MHMcEnergy*)[fDimension];
    fMustDelete = new Bool_t[fDimension];

    for (unsigned int i=0; i<fDimension; i++)
        fMustDelete[i]=kFALSE;
}

MMcThresholdCalc::~MMcThresholdCalc()
{
    for (unsigned int i=0; i<fDimension; i++)
        if (fMustDelete[i])
            delete fHMcEnergy[i];

    delete fMcTrig;
    delete fHMcEnergy;
    delete fMustDelete;
}

Bool_t MMcThresholdCalc::PreProcess(MParList* pList)
{
    // connect Monte Carlo data with this task

    char auxname[15]="MMcTrig"; // string to write container names

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

    for (unsigned int i=0; i<fDimension; i++)
    {
        if (fDimension>1)
            sprintf(auxname+7, ";%i", i+1);

        fMcTrig[i] = (MMcTrig*)pList->FindObject(auxname);
        if (fMcTrig[i])
            continue;

        *fLog << dbginf << "'MMcTrig";
        if (fDimension>1)
            *fLog << ";" << i+1;
        *fLog << "' not found... aborting." << endl;

        return kFALSE;
    }

    strcpy(auxname, "MHMcEnergy");
    for (unsigned int i=0; i<fDimension; i++)
    {
        if (fDimension>1)
            sprintf(auxname+10, ";%i", i+1);

        fHMcEnergy[i] = (MHMcEnergy*)pList->FindObject(auxname);
        if (fHMcEnergy[i])
            continue;

        *fLog << dbginf << "'" << auxname << "' not found in list... creating." << endl;

        fHMcEnergy[i] = new MHMcEnergy(fDimension>1 ? i+1 : 0);
        fMustDelete[i] = kTRUE;
        pList->AddToList(fHMcEnergy[i]);
    }

    return kTRUE;
}

Bool_t MMcThresholdCalc::Process()
{
    const Float_t energy   = fMcEvt->GetEnergy();
    const Float_t lg10     = log10(energy);
    const Float_t reciproc = 1./energy;

    for (unsigned int i=0; i<fDimension; i++)
    {
        if (fMcTrig[i]->GetFirstLevel()<=0)
            continue;

        fHMcEnergy[i]->Fill(lg10, reciproc);
    }

    return kTRUE;
}

Bool_t MMcThresholdCalc::PostProcess()
{
    // fit the energy distribution to get the threshold

    const Float_t sqrt2 = sqrt(2);

    for (unsigned int i=0; i<fDimension; i++)
    {
        MHMcEnergy &h = *fHMcEnergy[i];

        const Float_t peak  = h.GetGaussPeak();
        const Float_t sigma = h.GetGaussSigma();

        h.Fit(1, 3);
        h.Fit(peak - 2.   *sigma, peak + 2.   *sigma);
        h.Fit(peak - sqrt2*sigma, peak + sqrt2*sigma);
    }
    return kTRUE;
}