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

/* ======================================================================== *\
!
! *
! * 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);

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");
} 

Bool_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;
    }

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

    fTheta                    =  -1;
    fEmin                     =  -1;
    fEmax                     =  -1;
    fSlope                    =   0;
    fTotalNumSimulatedShowers =  0;
    fCorsikaVersion           =  0;
    fAllEvtsTriggered         = 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();

    if ( fEmin > 0 &&
         (fEmin  != corrunheader->GetELowLim() ||
          fEmax  != corrunheader->GetEUppLim() ||
          fSlope != corrunheader->GetSlopeSpec()) )
      *fLog << warn << dbginf << "Warning - Read files have different energy distribution..." << endl;

    fEmin  = corrunheader->GetELowLim();
    fEmax  = corrunheader->GetEUppLim();
    fSlope = corrunheader->GetSlopeSpec();

    fAllEvtsTriggered |= runheader->GetAllEvtsTriggered();

    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;

    if (fAllEvtsTriggered)
        return kTRUE;

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

    return kTRUE;
}

Bool_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;
}

Bool_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->CalcEfficiency(fTotalNumSimulatedShowers, fEmin, fEmax, fSlope);
    }

    return kTRUE;
}