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

/* ======================================================================== *\
!
! *
! * 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 of the pixels have values of lower than 3 Pedestal RMS. 
//   fMaxEmptyPixels is set to 0.4 by default which is slightly higher 
//   than the number of outer pixels in MAGIC (for the case that 
//   the outer pixels have some defect).
//
//  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;

const TString MFCosmics::fgNamePedestalCam = "MPedestalCam";

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MFCosmics::MFCosmics(const char *name, const char *title)
    : fPedestals(NULL), fSignals(NULL), fRawEvt(NULL),
    fNamePedestalCam(fgNamePedestalCam), fMaxEmptyPixels(0.2)
{
    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 << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
    }

    fPedestals = (MPedestalCam*)pList->FindObject(fNamePedestalCam, "MPedestalCam");
    if (!fPedestals)
    {
        *fLog << err << fNamePedestalCam << " [MPedestalCam] not found... aborting." << endl;
        return kFALSE;
    }

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


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


// --------------------------------------------------------------------------
//
// Initialize number of used FADC slices
//
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() const
{
    MRawEvtPixelIter pixel(fRawEvt);

    Int_t cosmicpix = 0;
    Int_t allpix    = 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();

        allpix++;

        //
        // 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.
    //
    return cosmicpix > fMaxEmptyPixels*allpix;
}

Int_t MFCosmics::PostProcess()
{
    if (GetNumExecutions()==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << dec << setfill(' ');

    *fLog << " " << setw(7) << fCut[0] << " (" << setw(3) ;
    *fLog << (int)(fCut[0]*100/GetNumExecutions()) ;
    *fLog << "%) Detected cosmics " ;
    *fLog << " (with fMaxEmptyPixels = " << fMaxEmptyPixels*100 << "%)" << endl;

    *fLog << " " << setw(7) << fCut[1] << " (" << setw(3) ;
    *fLog << (int)(fCut[1]*100/GetNumExecutions()) ;
    *fLog << "%) No cosmics!" << endl;
    *fLog << endl;

    return kTRUE;
}