source: trunk/MagicSoft/Mars/mfilter/MFCT1SelBasic.cc@ 1985

/* ======================================================================== *\
!
! *
! * 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): Wolfgang Wittek  04/2003 <mailto:wittek@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MFCT1SelBasic                                                          //
//                                                                         //
//  This is a class to evaluate basic cuts                                 //
//                                                                         //
//  to be called after the calibration (when the number of photons is      //
//               available for all pixels)                                  //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "MFCT1SelBasic.h"

#include "MParList.h"

#include "MHillas.h"
#include "MCerPhotEvt.h"
#include "MMcEvt.hxx"
#include "MRawRunHeader.h"
#include "MGeomCam.h"
#include "MPedestalCam.h"
#include "MPedestalPix.h"
#include "MGeomPix.h"

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

ClassImp(MFCT1SelBasic);

// --------------------------------------------------------------------------
//
// Default constructor.
//
MFCT1SelBasic::MFCT1SelBasic(const char *name, const char *title)
{
    fName  = name  ? name  : "MFCT1SelBasic";
    fTitle = title ? title : "Filter to evaluate basic cuts";

    // default values of cuts
    SetCuts(13.0, 0.0, 60.0);
}

// --------------------------------------------------------------------------
//
// Set the cut values
// 
//
void MFCT1SelBasic::SetCuts(Float_t minphotons, 
                            Float_t thetamin, Float_t thetamax)
{
  fMinPhotons = minphotons;
  fThetaMin   = thetamin;
  fThetaMax   = thetamax;

  *fLog << inf << "MFCT1SelBasic cut values : fMinPhotons, fThetaMin, fThetaMax = "
        << fMinPhotons <<",  " << fThetaMin << ",  " << fThetaMax << endl;
}

// --------------------------------------------------------------------------
//
// Set the pointers
// 
//
Bool_t MFCT1SelBasic::PreProcess(MParList *pList)
{
    fRawRun = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!fRawRun)
    {
        *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
        return kFALSE;
    }

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

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


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

    fPed = (MPedestalCam*)pList->FindObject("MPedestalCam");
    if (!fPed)
    {
        *fLog << dbginf << "MPedestalCam missing in Parameter List... aborting." << endl;
        return kFALSE;
    }

    memset(fCut, 0, sizeof(fCut));

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Evaluate basic cuts
// 
// if cuts are fulfilled      : fResult = kTRUE;
// if they are not fullfilled : fResult = kFALSE;
//
Bool_t MFCT1SelBasic::Process()
{
    Int_t rc = 0;
    fResult  = kFALSE;

    //if ( fRawRun->GetRunNumber() < 16279 ) 
    //{
    //   rc = 1;
    //   fResult = kTRUE;
    //}

    Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
    if ( theta < fThetaMin )
    {
      rc = 1;
      fResult = kTRUE;
    }    

    else if ( theta > fThetaMax )
    {
      rc = 2;
      fResult = kTRUE;
    }    

    else if ( !SwTrigger() )
    {
      //*fLog << "MFCT1SelBasic::Process; SwTrigger = " << SwTrigger << endl;
      rc = 3;
      fResult = kTRUE;
    }    

    fCut[rc]++;

    return kTRUE;
}
// --------------------------------------------------------------------------
//
// Software trigger
// 
// require 2 neighboring pixels (which are not in the outermost ring), 
//                       each having at least 'fMinPhotons' photons
// 
// 
Bool_t MFCT1SelBasic::SwTrigger()
{
    const Int_t entries = fEvt->GetNumPixels();
 
    for (Int_t i=0; i<entries; i++)
    {
      MCerPhotPix &pix = (*fEvt)[i];
      Int_t id = pix.GetPixId();
      if (!pix.IsPixelUsed()) continue;

      Double_t photons = pix.GetNumPhotons();
      if (photons < fMinPhotons) continue;

      // this pixel is used and has the required no.of photons
      // check whether this is also true for a neigboring pixel

      MGeomPix &gpix = (*fCam)[id];
      if ( gpix.IsInOutermostRing() ) continue;

      const Int_t nneighbors = gpix.GetNumNeighbors();
      for (Int_t n=0; n<nneighbors; n++)
      {
        const Int_t id1 =  gpix.GetNeighbor(n);
        if ( !fEvt->IsPixelUsed(id1) ) continue;

        MGeomPix &gpix1 = (*fCam)[id1];
        if ( gpix1.IsInOutermostRing() ) continue;

        MCerPhotPix &pix1 = *(fEvt->GetPixById(id1));      

        Double_t photons1 = pix1.GetNumPhotons();
        if (photons1 >= fMinPhotons) return kTRUE;
      }
    }
    return kFALSE;
}

// --------------------------------------------------------------------------
//
//  Prints some statistics about the Basic selections.
//
Bool_t MFCT1SelBasic::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << dec << setfill(' ');
    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) 
          << (int)(fCut[1]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Zenith angle < " << fThetaMin << endl;

    *fLog << " " << setw(7) << fCut[2] << " (" << setw(3) 
          << (int)(fCut[2]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Zenith angle > " << fThetaMax << endl;

    *fLog << " " << setw(7) << fCut[3] << " (" << setw(3) 
          << (int)(fCut[3]*100/GetNumExecutions()) 
          << "%) Evts skipped due to: Software trigger not fullfilled" 
          << " (with fMinPhotons = " << fMinPhotons << ")" << endl;

    *fLog << " " << fCut[0] << " (" << (int)(fCut[0]*100/GetNumExecutions()) 
          << "%) Evts survived Basic selections!" << endl;
    *fLog << endl;

    return kTRUE;
}