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

/* ======================================================================== *\
!
! *
! * 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
//                                                                         
//  Input Containers:
//   MMcEvt, MMcTrig;*
//
//  Output Containers:
//   MHMcEnergies
//
/////////////////////////////////////////////////////////////////////////////

#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)

// --------------------------------------------------------------------------
//
// Default Constructor.
//
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";

    // Arrays of MMcTrig* and MHMcEnergy* are created in order to be 
    // able to work with root files with several trigger conditions.
    fMcTrig     = new MMcTrig*[fDimension];
    fHMcEnergy  = new MHMcEnergy*[fDimension];
    fMustDelete = new Bool_t[fDimension];

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

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

    delete fMcTrig;
    delete fHMcEnergy;
    delete fMustDelete;
}

// --------------------------------------------------------------------------
//
// connect Monte Carlo data with this task
//
Bool_t MMcThresholdCalc::PreProcess(MParList* pList)
{
    // This task takes into accout if the root file has one trigger
    // condition (MMcTrig) or severl of them (MMcTrig;#.)

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

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

    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&&i!=0)
            sprintf(auxname+10, ";%i", i);

        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!=0 ? i : 0);
        fMustDelete[i] = kTRUE;
        pList->AddToList(fHMcEnergy[i]);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// The histograms are filled with log10 of the energy for triggered
// events and weighted with 1/E because it is needed the dN/dE vs. logE
// distribution to get the energy threshold.
//
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)
            fHMcEnergy[i]->Fill(lg10, reciproc);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// fit the energy distribution to get the threshold
// Some iterations are done to be sure the fit parameters converge.
//
Bool_t MMcThresholdCalc::PostProcess()
{
    const Float_t sqrt2 = sqrt(2);

    for (unsigned int i=0; i<fDimension; i++)
    {
      Float_t peak;
      Float_t sigma;

      fHMcEnergy[i]->Fit(1, 3);
      
      peak  = fHMcEnergy[i]->GetGaussPeak();
      sigma = fHMcEnergy[i]->GetGaussSigma();
      fHMcEnergy[i]->Fit(peak - 2.   *sigma, peak + 2.   *sigma);
        
      peak  = fHMcEnergy[i]->GetGaussPeak();
      sigma = fHMcEnergy[i]->GetGaussSigma();
      fHMcEnergy[i]->Fit(peak - sqrt2*sigma, peak + sqrt2*sigma);
    }
    return kTRUE;
}