source: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc@ 4494

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 02/2004 <mailto:tbretz@astro.uni.wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MBadPixelsCalc
//
//
// The job of the task is to determin bad pixels event-wise. This must be
// redone for each event. This particular task is for what is explained
// below.
// New algorithms may enter new tasks.
//
//
// Check the pedestal RMS of every pixel with respect to the mean
// pedestal RMS of the camera (Sigmabar).
//
// The pixels can be set as blind if the pedestalRMS is too big or 0.
//
// If you don't want to use this option set the PedestalLevel<=0;
//
//     MBadPixelsCalc calc;
//     calc.SetPedestalLevel(-1);
//
//
//  Input Containers:
//   [MPedPhotCam]
//   [MGeomCam]
//   [MSigmabar]
//
//  Output Containers:
//   MBadPixels
//
/////////////////////////////////////////////////////////////////////////////
#include "MBadPixelsCalc.h"

#include <TArrayD.h>

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MSigmabar.h"

#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

#include "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

ClassImp(MBadPixelsCalc);

using namespace std;

static const TString gsDefName  = "MBadPixelsCalc";
static const TString gsDefTitle = "Find hot spots (star, broken pixels, etc)";

// --------------------------------------------------------------------------
//
// Default constructor.
//
MBadPixelsCalc::MBadPixelsCalc(const char *name, const char *title)
    : fPedestalLevel(3)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();
}

// --------------------------------------------------------------------------
//
//
Int_t MBadPixelsCalc::PreProcess (MParList *pList)
{
    fBadPixels = (MBadPixelsCam*)pList->FindCreateObj(AddSerialNumber("MBadPixelsCam"));
    if (!fBadPixels)
        return kFALSE;

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

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

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

// --------------------------------------------------------------------------
//
// Check the pedestal RMS of every pixel with respect to the mean pedestal RMS 
// of the camera (Sigmabar).
//
// The pixels can be set as blind if the pedestalRMS is too big or 0.
//
// If you don't want to use this option set the PedestalLevel<=0;
//
//     MBadPixelsCalc calc;
//     calc.SetPedestalLevel(-1);
//
void MBadPixelsCalc::CheckPedestalRMS() const
{
    const Int_t entries = fPedPhotCam->GetSize();

    const Float_t meanPedRMS = fSigmabar->GetSigmabar();

    for (Int_t i=0; i<entries; i++)
    {
        //
        // get pixel as entry from list
        //
        const Double_t nratio    = fGeomCam->GetPixRatio(i);
        const Double_t pixPedRms = (*fPedPhotCam)[i].GetRms();

        if (pixPedRms*nratio > fPedestalLevel * meanPedRMS || pixPedRms == 0)
            (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
    }
}

// --------------------------------------------------------------------------
// Check the pedestal Rms of the pixels: compute with 2 iterations the mean 
// for inner and outer pixels. Set as blind the pixels with too small or 
// too high pedestal Rms with respect to the mean.
// 
Bool_t MBadPixelsCalc::CheckPedestalRms() const
{
    const Int_t entries = fPedPhotCam->GetSize();

    const Int_t na = fGeomCam->GetNumAreas();

    TArrayD meanrms(na);
    TArrayI npix(na);

    for (Int_t i=0; i<entries; i++)
    {
        const Double_t rms = (*fPedPhotCam)[i].GetRms();

        if (rms<=0 || rms>=200*fGeomCam->GetPixRatioSqrt(i))
            continue;

        const Byte_t aidx = (*fGeomCam)[i].GetAidx();

        meanrms[aidx] += rms;
        npix[aidx]++;
    }

    //if no pixel has a minimum signal, return
    for (int i=0; i<na; i++)
    {
        if (npix[i]==0 || meanrms[i]==0)
        {
            //fErrors[1]++;          //no valid Pedestals Rms
            return kFALSE;
        }

        meanrms[i] /= npix[i];
        npix[i]=0;
    }

    TArrayD meanrms2(na);
    for (Int_t i=0; i<entries; i++)
    {
        const Double_t rms = (*fPedPhotCam)[i].GetRms();
        const Byte_t  aidx = (*fGeomCam)[i].GetAidx();

        //Calculate the corrected means:

        if (rms<=0.5*meanrms[aidx] || rms>=1.5*meanrms[aidx])
            continue;

        meanrms2[aidx] += rms;
        npix[aidx]++;
    }

    //if no pixel has a minimum signal, return
    for (int i=0; i<na; i++)
    {
        if (npix[i]==0 || meanrms2[i]==0)
        {
            //fErrors[1]++;          //no valid Pedestals Rms
            return kFALSE;
        }

        meanrms2[i] /= npix[i];
    }

    Int_t bads = 0;

    //Blind the Bad Pixels
    for (Int_t i=0; i<entries; i++)
    {
        const Double_t rms = (*fPedPhotCam)[i].GetRms();
        const Byte_t  aidx = (*fGeomCam)[i].GetAidx();

        if (rms>meanrms2[aidx]/3 && rms<=meanrms2[aidx]*3)
            continue;

        (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt);
        bads++;
    }

    // Check if the number of pixels to blind is > 60% of total number of pixels
    //
    if (bads>0.6*entries)
    {
        //fErrors[2]++;
        return kFALSE;
    }

    return kTRUE;
}


// --------------------------------------------------------------------------
//
//
Int_t MBadPixelsCalc::Process()
{
    if (fPedestalLevel>0)
        CheckPedestalRms();

    return kTRUE;
}