source: trunk/MagicSoft/Mars/manalysis/MSelBasic.cc@ 1895

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

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

#include "MSelBasic.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(MSelBasic);

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

    fThetaMin =  0.0;
    fThetaMax = 60.0;
}

// --------------------------------------------------------------------------
//
// 
// 
// 
//
Bool_t MSelBasic::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      : return 0
// if they are not fullfilled : skip remaining tasks for this event
//
Bool_t MSelBasic::Process()
{
    Int_t rc = 0;

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

    Double_t theta = kRad2Deg*fMcEvt->GetTelescopeTheta();
    if ( theta < fThetaMin )
    {
      *fLog << "MSelBasic::Process; Run, Event, Theta = " 
            << fRawRun->GetRunNumber()<< ",  " 
            << fMcEvt->GetEvtNumber() << ",  " << theta << endl;
      rc = 1;
    }    

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

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

    fCut[rc]++;

    return rc==0 ? kTRUE : kCONTINUE;
}
// --------------------------------------------------------------------------
//
// Software trigger
// 
// require 2 neighboring pixels (which are not in the outermost ring), 
//                       each having at least 'minphotons' photons
// 
// 
Bool_t MSelBasic::SwTrigger()
{
    // minimum number of photons required
    Double_t minphotons = 13.0;

    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 < minphotons) 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)[id1];      
        MCerPhotPix &pix1 = *(fEvt->GetPixById(id1));      
        if ( &pix1 == NULL ) 
        {
          *fLog << "MSelBasic::SwTrigger; &pix1 is NULL" << endl;
          continue;
        }

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

// --------------------------------------------------------------------------
//
//  Prints some statistics about the Basic selections.
//
Bool_t MSelBasic::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" << endl;

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

    return kTRUE;
}