source: trunk/MagicSoft/Mars/manalysis/MImgCleanStd.cc@ 761

/* ======================================================================== *\
!
! *
! * 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): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!   Author(s): Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

#include "MImgCleanStd.h"

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

#include "MParList.h"
#include "MGeomCam.h"
#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"

ClassImp(MImgCleanStd)

MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
                           const char *name, const char *title)
    : fCleanLvl1(lvl1), fCleanLvl2(lvl2)
{
    //
    //   the default constructor
    //
    *fName  = name  ? name  : "MImgCleanStd";
    *fTitle = title ? title : "Task which does a standard image cleaning";

    *fLog << "Cleaning initialized. Using noise level " << lvl1 << " and " << lvl2 << endl;
}

void MImgCleanStd::CleanStep1()
{
    //
    //  This method looks for all pixels with an entry (photons)
    //  that is three times bigger than the noise of the pixel
    //  (std: 3 sigma, clean level 1)
    //

    const Int_t entries = fEvt->GetNumPixels();

    //
    // check the number of all pixels against the noise level and
    // set them to 'unused' state if necessary
    //
    for (Int_t i=0; i<entries; i++ )
    {
        MCerPhotPix &pix = (*fEvt)[i];

        const Float_t entry = pix.GetNumPhotons();
        const Float_t noise = pix.GetErrorPhot();

        // COBB: '<=' to skip entry=noise=0
        if (entry <= fCleanLvl1 * noise )
            pix.SetPixelUnused();
    }
}

void MImgCleanStd::CleanStep2()
{
    //
    //  check if the survived pixel have a neighbor, that also
    //  survived
    //

    const Int_t entries = fEvt->GetNumPixels();

    for (Int_t i=0; i<entries; i++)
    {
        //
        // get entry i from list
        //
        MCerPhotPix &pix = (*fEvt)[i];

        //
        // check if pixel is in use, if not goto next pixel in list
        //
        if (!pix.IsPixelUsed())
            continue;

        //
        // get pixel id of this entry
        //
        const Int_t id = pix.GetPixId();

        //
        // count number of next neighbors of this pixel which
        // state is 'used'
        //
        MGeomPix   &gpix  = (*fCam)[id];
        const Int_t nnmax = gpix.GetNumNeighbors();

        Int_t cnt = 0;
        for (Int_t j=0; j<nnmax; j++)
        {
            const Int_t id2 = gpix.GetNeighbor(j); //GetNN(id, in) ;

            if (id2 < 0)
                continue;

            if (fEvt->IsPixelUsed(id2))
                cnt++;
        }

        //
        // check if no next neighbor has the state 'used'
        // set this pixel to 'unused', too.
        //
        if (cnt==0)
            pix.SetPixelUnused();
    }

    //
    // now we declare all pixels that survive as CorePixels
    //
    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];

        if (pix.IsPixelUsed())
            pix.SetCorePixel();
    }
} 

void MImgCleanStd::CleanStep3()
{
    //
    //   Look for the boundary pixels around the core pixels
    //   if a pixel has more than 2.5 (clean level 2) sigma, and
    //   a core neigbor it is declared as used.
    //
    const Int_t entries = fEvt->GetNumPixels();

    for (Int_t i=0; i<entries; i++)
    {
        //
        // get pixel as entry il from list
        //
        MCerPhotPix &pix = (*fEvt)[i];

        //
        // if pixel is a core pixel go to the next pixel
        //
        if (pix.IsCorePixel())
            continue;

        //
        // check the num of photons against the noise level
        //
        const Float_t entry = pix.GetNumPhotons();
        const Float_t noise = pix.GetErrorPhot();

        if (entry <= fCleanLvl2 * noise )
            continue;

        //
        // get pixel id of this entry
        //
        const Int_t id = pix.GetPixId();

        //
        // check if the pixel's next neighbor is a core pixel.
        // if it is a core pixel set pixel state to: used.
        //
        MGeomPix   &gpix  = (*fCam)[id];
        const Int_t nnmax = gpix.GetNumNeighbors();

        for (Int_t j=0; j<nnmax; j++)
        {
            const Int_t id2 = gpix.GetNeighbor(j);

            if (id2 <0)
                continue;

            if (!fEvt->IsPixelCore(id2))
                continue;

            pix.SetPixelUsed();

            break ;
        }
    }
}

Bool_t MImgCleanStd::PreProcess (MParList *pList)
{
    //
    //  check if MEvtHeader exists in the Parameter list already.
    //  if not create one and add them to the list
    //
    fCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fCam)
    {
        *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
        return kFALSE;
    }

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

    return kTRUE;
}
    
Bool_t MImgCleanStd::Process()
{
    CleanStep1();
    CleanStep2();
    CleanStep3();

    return kTRUE;
}