source: trunk/MagicSoft/Mars/mmontecarlo/MMcCollectionAreaCalc.cc@ 5618

/* ======================================================================== *\
!
! *
! * 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  12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MHMcCollectionAreaCalc
//
//  Remark: The initialization is mainly done in the ReInit function.
//          Please make sure, that you don't use MReadTree when processing
//          a file. Use a 'ReInit'-calling task like MReadMarsFile
//
//////////////////////////////////////////////////////////////////////////////

#include "MMcCollectionAreaCalc.h"

#include "MParList.h"

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

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

#include "MHMcCollectionArea.h"

ClassImp(MMcCollectionAreaCalc);

using namespace std;

MMcCollectionAreaCalc::MMcCollectionAreaCalc(const char *input,
                                             const char *name, const char *title)
{
    fName  = name  ? name  : "MMcCollectionAreaCalc";
    fTitle = title ? title : "Task to calculate the collection area";

    fObjName = input ? input : "MMcTrig";
    AddToBranchList(Form("%s.fNumFirstLevel", input));

    AddToBranchList("MMcEvt.fEnergy");
    AddToBranchList("MMcEvt.fImpact");

    fCorsikaVersion           =  0;
    fAllEvtsTriggered         =  kFALSE;
    fTotalNumSimulatedShowers =  0;
    fThetaMin                 = -1.;
    fThetaMax                 = -1.;

} 

Int_t MMcCollectionAreaCalc::PreProcess (MParList *pList)
{
    // connect the raw data with this task

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

    return kTRUE;
}

Bool_t MMcCollectionAreaCalc::ReInit(MParList *plist)
{
    fCollArea=0;

    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
    if (!runheader)
    {
        *fLog << err << "Error - MMcRunHeader not found... abort." << endl;
        return kFALSE;
    }

    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();
    *fLog << inf << "Total Number of Simulated showers: " << fTotalNumSimulatedShowers << endl;


    if ( (fThetaMin >= 0 && fThetaMin != runheader->GetShowerThetaMin()) ||
         (fThetaMax >= 0 && fThetaMax != runheader->GetShowerThetaMax()) )
    {
        *fLog << warn << "Warning - Read files have different Theta ranges (";
        *fLog << "(" << fThetaMin << ", " << fThetaMax << ")   vs   (" << 
          runheader->GetShowerThetaMin() << ", " << 
          runheader->GetShowerThetaMax() << ")..." << endl;
    }
    fThetaMin = runheader->GetShowerThetaMin();
    fThetaMax = runheader->GetShowerThetaMax();


    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
        *fLog << warn << "Warning - Read files have different Corsika versions..." << endl;
    fCorsikaVersion = runheader->GetCorsikaVersion();


    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();
    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;


    fCollArea = (MHMcCollectionArea*)plist->FindCreateObj("MHMcCollectionArea");
    if (!fCollArea)
        return kFALSE;

    if (!fAllEvtsTriggered)
    {
        fMcTrig = (MMcTrig*)plist->FindObject(fObjName);
        if (!fMcTrig)
        {
            *fLog << err << fObjName << " [MMcTrig] not found... anort." << endl;
            return kFALSE;
        }
        return kTRUE;
    }

    MMcCorsikaRunHeader *corrunheader  = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
    if (!corrunheader)
    {
        *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
        return kFALSE;
    }

    //
    // Calculate approximately the original number of events in each
    // energy bin:
    //

    const Float_t emin = corrunheader->GetELowLim();
    const Float_t emax = corrunheader->GetEUppLim();
    const Float_t expo = 1 + corrunheader->GetSlopeSpec();
    const Float_t k = runheader->GetNumSimulatedShowers() /
        (pow(emax,expo) - pow(emin,expo)) ;

    TH2 &hall = *fCollArea->GetHistAll();

    const Int_t nbinx = hall.GetNbinsX();

    TAxis &axe = *hall.GetXaxis();
    for (Int_t i = 1; i <= nbinx; i++)
    {
        const Float_t e1 = axe.GetBinLowEdge(i);
        const Float_t e2 = axe.GetBinLowEdge(i+1);

        if (e1 < emin || e2 > emax)
            continue;

        const Float_t events = k * (pow(e2, expo) - pow(e1, expo));
        //
        // We fill the i-th energy bin, with the total number of events
        // Second argument of Fill would be impact parameter of each
        // event, but we don't really need it for the collection area,
        // so we just put a dummy value (1.)
        //

        const Float_t energy = (e1+e2)/2.;
        hall.Fill(energy, 1., events);
    }

    return kTRUE;
}

Int_t MMcCollectionAreaCalc::Process()
{
    const Double_t energy = fMcEvt->GetEnergy();
    const Double_t impact = fMcEvt->GetImpact()/100.;

    //
    // This happens for camera files created with Camera 0.5
    //
    if (TMath::IsNaN(impact))
    {
        *fLog << err << dbginf << "ERROR - Impact=NaN (Not a number)... abort." << endl;
        return kERROR;
    }

    if (!fAllEvtsTriggered)
    {
        fCollArea->FillAll(energy, impact);

        if (fMcTrig->GetFirstLevel() <= 0)
            return kTRUE;
    }

    fCollArea->FillSel(energy, impact);

    return kTRUE;
}

Int_t MMcCollectionAreaCalc::PostProcess()
{
    if (!fCollArea)
        return kTRUE;

    //
    //   do the calculation of the effectiv area
    //
    *fLog << inf << "Calculation Collection Area..." << endl;

    if (!fAllEvtsTriggered)
        fCollArea->CalcEfficiency();
    else
    {
        *fLog << inf << "Total number of showers: " << fTotalNumSimulatedShowers << endl;
        fCollArea->CalcEfficiency2();
    }

    return kTRUE;
}