source: trunk/MagicSoft/Mars/manalysis/MSigmabarCalc.cc@ 2946

/* ======================================================================== *\
!
! *
! * 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): Robert Wagner, 10/2002 <rwagner@mppmu.mpg.de>
!   Author(s): Thomas Bretz, 4/2003 <tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MSigmabarCalc
//
//  This task calculates Sigmabar using the MSigmabar container and stores
//  its extremal values together with the corresponding Theta values       
//  in MSigmabarParam. For the time being, Theta is taken from a Monte
//  Carlo container. This is preliminary and to be changed ASAP.
//
//  Input Containers:
//   MGeomCam
//   MPedPhotCam
//   MRawRunHeader
//   MMcEvt  (FIXME: Must be replaced by a 'real-data' container)
//   MCerPhotEvt
//
//  Output Containers:
//   MSigmabar
//   MSigmabarParam
//
/////////////////////////////////////////////////////////////////////////////
#include "MSigmabarCalc.h"

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MPedPhotCam.h"

#include "MSigmabar.h"
#include "MSigmabarParam.h"

#include "MMcEvt.hxx"

ClassImp(MSigmabarCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MSigmabarCalc::MSigmabarCalc(const char *name, const char *title)
{
    fName  = name  ? name  : "MSigmabarCalc";
    fTitle = title ? title : "Task to calculate Sigmabar";

    Reset();
}

// --------------------------------------------------------------------------
//
//  check if necessary containers exists in the Parameter list already.
//  if not create one and add them to the list
//
Int_t MSigmabarCalc::PreProcess(MParList *pList)
{
    fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!fCam)
    {
        *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }

    fPed = (MPedPhotCam*)pList->FindObject(AddSerialNumber("MPedPhotCam"));
    if (!fPed)
    {
        *fLog << err << "MPedPhotCam not found... aborting." << endl;
        return kFALSE;
    }

    fRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRun)
    {
        *fLog << err << "MRawRunHeader not found... aborting." << endl;
        return kFALSE;
    }

    // This is needed for determining min/max Theta
    fMcEvt = (MMcEvt*)pList->FindObject(AddSerialNumber("MMcEvt"));
    if (!fMcEvt)
    {
        *fLog << warn << "MMcEvt not found... aborting." << endl;
        fThetaMin =  0;
        fThetaMax = 90;
    }

    fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
    if (!fEvt)
      {
        *fLog << err << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
      }

    fSig = (MSigmabar*)pList->FindCreateObj(AddSerialNumber("MSigmabar"));
    if (!fSig)
        return kFALSE;

    // FIXME: This might have to change in a
    //        'more-than-one-telescope-environment'
    fSigParam = (MSigmabarParam*)pList->FindCreateObj(AddSerialNumber("MSigmabarParam"));
    if (!fSigParam)
        return kFALSE;


    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculating a new Sigmabar is not necessary on event basis as long as
// we deal with CT1 data. Therefore, the real calculation is done in
// the ReInit function. Process just takes care for finding the extremal
// values. Preliminary.
//
Int_t MSigmabarCalc::Process()
{
    const Double_t rc = fSig->Calc(*fCam, *fPed, *fEvt);

    if (rc>fSigmabarMax) fSigmabarMax=rc;
    if (rc<fSigmabarMin) fSigmabarMin=rc;

    if (!fMcEvt)
        return kTRUE;

    const Double_t theta = fMcEvt->GetTelescopeTheta()*kRad2Deg;

    if (theta>fSigmabarMax) fSigmabarMax=theta;
    if (theta<fSigmabarMax) fSigmabarMin=theta;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculates Sigmabar (for CT1 only needed at Reinit, i.e. when reading
// a new file)
//
Bool_t MSigmabarCalc::ReInit(MParList *pList)
{
    fSigParam->SetParams(1, fSigmabarMin, fSigmabarMax, fThetaMin, fThetaMax);
    fSigParam->SetRunNumber(fRun->GetRunNumber());

    Reset();

    return kTRUE;
}

void MSigmabarCalc::Reset()
{
    if (!fMcEvt)
        return;

    fThetaMin    = FLT_MAX;
    fThetaMax    = 0;
    fSigmabarMin = FLT_MAX;
    fSigmabarMax = 0;
}