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

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MImgCleanStd                                                           //
//                                                                         //
//  This is the standard image cleaning. If you want to know how it works  //
//  Please look at the three CleanSteps and Process                        //
//                                                                         //
//   FIXME: MImgCleanStd is not yet completely optimized for speed.        //
//          Maybe we don't have to loop over all pixels all the time...    //
//                                                                         //
//  Input Containers:                                                      //
//   MGeomCam, MCerPhotEvt                                                 //
//                                                                         //
//  Output Containers:                                                     //
//   -/-                                                                   //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MImgCleanStd.h"

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

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

ClassImp(MImgCleanStd);

// --------------------------------------------------------------------------
//
// Default constructor. Here you can specify the cleaning levels. If you
// don't specify them the 'common standard' values 2.5 and 3.0 (sigma
// above mean) are used
//
MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
                           const char *name, const char *title)
    : fCleanLvl1(lvl1), fCleanLvl2(lvl2)
{
    fName  = name  ? name  : "MImgCleanStd";
    fTitle = title ? title : "Task which does a standard image cleaning";

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

// --------------------------------------------------------------------------
//
//  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)
//
void MImgCleanStd::CleanStep1()
{
    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();
    }
}

// --------------------------------------------------------------------------
//
//  check if the survived pixel have a neighbor, that also
//  survived
//
void MImgCleanStd::CleanStep2()
{
    const Int_t entries = fEvt->GetNumPixels();

    //
    // In the worst case we have to loop 6 times 577 times, to
    // catch the behaviour of all next neighbors. Here we can gain
    // much by using an array instead of checking through all pixels
    // (MCerPhotEvt::IsPixelUsed) all the time.
    //
    Byte_t ispixused[577];

    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];
        ispixused[pix.GetPixId()] = pix.IsPixelUsed();
    }

    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'
        //
        const MGeomPix &gpix  = (*fCam)[id];
        const Int_t     nnmax = gpix.GetNumNeighbors();

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

            if (!ispixused[id2])
                continue;

            cnt = kTRUE;
            break;
        }
        if (cnt)
            continue;

        //
        // check if no next neighbor has the state 'used'
        // set this pixel to 'unused', too.
        //
        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();
    }
} 

// --------------------------------------------------------------------------
//
//   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.
//
void MImgCleanStd::CleanStep3()
{
    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 (!fEvt->IsPixelCore(id2))
                continue;

            pix.SetPixelUsed();
            break;
        }
    }
}

// --------------------------------------------------------------------------
//
//  check if MEvtHeader exists in the Parameter list already.
//  if not create one and add them to the list
//
Bool_t MImgCleanStd::PreProcess (MParList *pList)
{
    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;
}
    
// --------------------------------------------------------------------------
//
// Cleans the image.
//
Bool_t MImgCleanStd::Process()
{
    CleanStep1();
    CleanStep2();
    CleanStep3();

    return kTRUE;
}