source: trunk/MagicSoft/Mars/mfilter/MFCosmics.cc@ 3084

/* ======================================================================== *\
!
! *
! * 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): Markus Gaug  02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MFCosmics
//
//   Filter to reject cosmics by the criterion that at least 
//   fMaxEmptyPixels pixels have values of lower than 3 Pedestal RMS. 
//   fMaxEmptyPixels is set to 230 by default which is slightly higher 
//   than the number of outer pixels in MAGIC (for the case that 
//   the outer pixels have some defect).
//
//   ProProcess:  Search for MPedestalCam, MExtractedSignalCam
//               
//   ReInit:      Initialize number of used FADC slices
//
//   Process:     if fMaxEmptyPixels pixels lower than 3 pedRMS, 
//                the event is supposed to 
//                be a cosmic and kContinue is returned
//
//
//  Input Containers:
//   MRawEvtData
//   MPedestalCam
//   MExtractedSignalCam
//
//  Output Containers:
//
//////////////////////////////////////////////////////////////////////////////
#include "MFCosmics.h"

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MRawEvtPixelIter.h"

#include "MPedestalCam.h"
#include "MPedestalPix.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

ClassImp(MFCosmics);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
MFCosmics::MFCosmics(const char *name, const char *title)
    : fPedestals(NULL), fSignals(NULL),
      fRawEvt(NULL), fMaxEmptyPixels(230)
{

    fName  = name  ? name  : "MFCosmics";
    fTitle = title ? title : "Filter to reject cosmics";
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//  - MExtractedSignalCam
//
Int_t MFCosmics::PreProcess(MParList *pList)
{

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

    fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
    if (!fPedestals)
      {
        *fLog << err << dbginf << "Cannot find MPedestalCam ... aborting" << endl;
        return kFALSE;
      }

    fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
    if (!fSignals)
      {
        *fLog << err << dbginf << "Cannot find MExtractedSignalCam ... aborting" << endl;
        return kFALSE;
      }

    memset(fCut, 0, sizeof(fCut));
    
    return kTRUE;
}


// --------------------------------------------------------------------------
//
// The ReInit searches the following input containers for information:
//  - MExtractedSignalCam
//
Bool_t MFCosmics::ReInit(MParList *pList )
{
 
    fSqrtHiGainSamples = TMath::Sqrt((Float_t) fSignals->GetNumUsedHiGainFADCSlices());

    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Retrieve the integral of the FADC time slices and compare them to the 
// pedestal values.
//
Int_t MFCosmics::Process()
{

  fResult = CosmicsRejection();
  
  fCut[fResult ? 0 : 1]++;
  return kTRUE;
}

// ---------------------------------------------------------
//
// Cosmics rejection: 
// 
// Requiring less than fMaxEmptyPixels pixels to have values 
// lower than 3 Pedestal RMS. 
// 
// fMaxEmptyPixels is set to 230 by default which is slightly higher 
// than the number of outer pixels in MAGIC (for the case that 
// the outer pixels have some defect).
//
Bool_t MFCosmics::CosmicsRejection()
{
  
  MRawEvtPixelIter pixel(fRawEvt);
  
  Int_t cosmicpix = 0;
      
  //
  // Create a first loop to sort out the cosmics ...
  // 
  while (pixel.Next())
    {
      
      const UInt_t idx = pixel.GetPixelId();
      
      MExtractedSignalPix &sig =  (*fSignals)[idx];
      MPedestalPix        &ped =  (*fPedestals)[idx];
      const Float_t pedrms      = ped.GetPedestalRms()*fSqrtHiGainSamples;
      const Float_t sumhi       = sig.GetExtractedSignalHiGain();
          
      //
      // We consider a pixel as presumably due to cosmics 
      // if its sum of FADC slices is lower than 3 pedestal RMS
      //
      if (sumhi < 3.*pedrms )  
        cosmicpix++;
    }
  
  //
  // If the camera contains more than fMaxEmptyPixels
  // presumed pixels due to cosmics, then the event is discarted. 
  //
  if (cosmicpix > fMaxEmptyPixels)
    return kTRUE;

  return kFALSE;
}

Int_t MFCosmics::PostProcess()
{

  if (GetNumExecutions()==0)
    return kTRUE;
  
  *fLog << inf << endl;
  *fLog << GetDescriptor() << " execution statistics:" << endl;
  *fLog << dec << setfill(' ');
  
  *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ;
  *fLog << (int)(fCut[1]*100/GetNumExecutions()) ;
  *fLog << "%) Evts skipped due to: Cosmics Rejection applied " ;
  *fLog << " (with fMaxEmptyPixels = " << fMaxEmptyPixels << ")" << endl;

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

  return kTRUE;
}