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

/* ======================================================================== *\
!
! *
! * 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;
    fTheta                    = -1.;

} 

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

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

    return kTRUE;
}

Bool_t MMcCollectionAreaCalc::ReInit(MParList *plist)
{
    MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader");
    if (!runheader)
    {
        *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl;
        return kFALSE;
    }

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

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


    if (fTheta>=0 && fTheta!=runheader->GetTelesTheta())
    {
        *fLog << warn << dbginf << "Warning - Read files have different TelesTheta (";
        *fLog << fTheta << ", " << runheader->GetTelesTheta() << ")..." << endl;
    }
    fTheta = runheader->GetTelesTheta();


    if (fCorsikaVersion!=0 && fCorsikaVersion!=runheader->GetCorsikaVersion())
        *fLog << warn << dbginf << "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)
      {
        //
        // 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;
      }

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

    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()
{ 
    //
    //   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;
}