source: trunk/MagicSoft/Mars/mimage/MHillasCalc.cc@ 2640

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

/////////////////////////////////////////////////////////////////////////////
//
//  MHillasCalc
//
//  This is a task to calculate the Hillas parameters from each event
//
//  By default MHillas, MHillasExt and MNewImagePar are calculated
//  with the information from MCerPhotEvt and MGeomCam.
//
//  To switch of the calculation you may use:
//   - Disable(MHillasCalc::kCalcHillas)
//   - Disable(MHillasCalc::kCalcHillasExt)
//   - Disable(MHillasCalc::kCalcNewImagePar)
//
//  If the calculation of MHillas is switched off a container MHillas
//  in the parameter list is nevertheless necessary for the calculation
//  of MHillasExt and MNewImagePar.
//
//  The names of the containers to be used can be set with:
//   - SetNameHillas("NewName")
//   - SetNameHillasExt("NewName")
//   - SetNameNewImgPar("NewName")
//
//  Input Containers:
//   MCerPhotEvt
//   MGeomCam
//   [MHillas]
//
//  Output Containers:
//   [MHillas]
//   MHillasExt
//   MNewImagePar
//
/////////////////////////////////////////////////////////////////////////////
#include "MHillasCalc.h"

#include "MParList.h"

#include "MHillas.h"
#include "MHillasExt.h"
#include "MNewImagePar.h"

#include "MCerPhotEvt.h"

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

ClassImp(MHillasCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MHillasCalc::MHillasCalc(const char *name, const char *title)
    : fHilName("MHillas"), fHilExtName("MHillasExt"),
    fImgParName("MNewImagePar"), fFlags(0xff), fErrors(5)
{
    fName  = name  ? name  : "MHillasCalc";
    fTitle = title ? title : "Calculate Hillas and other image parameters";
}

// --------------------------------------------------------------------------
//
// Check for a MCerPhotEvt object from which the Hillas are calculated.
// Try to find the Geometry conatiner. Depending on the flags
// try to find (and maybe create) the containers MHillas, MHillasExt,
// MNewImagePar, too.
//
Int_t MHillasCalc::PreProcess(MParList *pList)
{
    // necessary
    fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
    if (!fCerPhotEvt)
    {
        *fLog << err << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

    // necessary
    fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!fGeomCam)
    {
        *fLog << err << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
        return kFALSE;
    }

    // depend on whether MHillas is an in- or output container
    if (TestFlag(kCalcHillas))
        fHillas = (MHillas*)pList->FindCreateObj("MHillas", AddSerialNumber(fHilName));
    else
    {
        fHillas = (MHillas*)pList->FindObject(AddSerialNumber(fHilName), "MHillas");
        *fLog << err << fHilName << " [MHillas] not found... aborting." << endl;
    }
    if (!fHillas)
        return kFALSE;

    // if enabled
    if (TestFlag(kCalcHillasExt))
    {
        fHillasExt = (MHillasExt*)pList->FindCreateObj("MHillasExt", AddSerialNumber(fHilExtName));
        if (!fHillasExt)
            return kFALSE;
    }

    // if enabled
    if (TestFlag(kCalcNewImagePar))
    {
        fNewImgPar = (MNewImagePar*)pList->FindCreateObj("MNewImagePar", AddSerialNumber(fImgParName));
        if (!fNewImgPar)
            return kFALSE;
    }

    memset(fErrors.GetArray(), 0, sizeof(Char_t)*fErrors.GetSize());

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// If you want do complex descisions inside the calculations
// we must move the calculation code inside this function
//
// If the calculation wasn't sucessfull skip this event
//
Int_t MHillasCalc::Process()
{
    if (TestFlag(kCalcHillas))
    {
        Int_t rc = fHillas->Calc(*fGeomCam, *fCerPhotEvt);
        if (rc<0 || rc>4)
        {
            *fLog << err << dbginf << "MHillas::Calc returned unknown error code!" << endl;
            return kFALSE;
        }
        fErrors[rc]++;
        if (rc>0)
            return kCONTINUE;
    }

    if (TestFlag(kCalcHillasExt))
        fHillasExt->Calc(*fGeomCam, *fCerPhotEvt, *fHillas);

    if (TestFlag(kCalcNewImagePar))
        fNewImgPar->Calc(*fGeomCam, *fCerPhotEvt, *fHillas);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  This is used to print the output in the PostProcess. Later (having
//  millions of events) we may want to improve the output.
//
void MHillasCalc::PrintSkipped(int i, const char *str) const
{
    *fLog << " " << setw(7) << fErrors[i] << " (";
    *fLog << setw(3) << (int)(100.*fErrors[i]/GetNumExecutions());
    *fLog << "%) Evts skipped due to: " << str << endl;
}

// --------------------------------------------------------------------------
//
//  Prints some statistics about the hillas calculation. The percentage
//  is calculated with respect to the number of executions of this task.
//
Int_t MHillasCalc::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << dec << setfill(' ');
    PrintSkipped(1, "Event has less than 3 pixels");
    PrintSkipped(2, "Calculated Size == 0");
    PrintSkipped(3, "Number of used pixels < 3");
    PrintSkipped(4, "CorrXY==0");
    *fLog << " " << (int)fErrors[0] << " (" << (int)(100.*fErrors[0]/GetNumExecutions()) << "%) Evts survived Hillas calculation!" << endl;
    *fLog << endl;

    return kTRUE;
}