source: trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc@ 2591

/* ======================================================================== *\
!
! *
! * 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, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer, 1/2001
!   Author(s): Nadia Tonello, 4/2003 <mailto:tonello@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MImgCleanStd
//
// The Image Cleaning task selects the pixels you use for the Hillas
// parameters calculation.
//
// There are two methods to make the selection: the standard one, as done
// in the analysis of CT1 data, and the democratic one, as suggested by
// W.Wittek. The number of photo-electrons of a pixel is compared with the
// pedestal RMS of the pixel itself (standard method) or with the average
// RMS of the inner pixels (democratic method).
// In both cases, the possibility to have a camera with pixels of
// different area is taken into account.
// The too noisy pixels can be recognized and eventally switched off
// (Unmap: set blind pixels to UNUSED) separately, using the
// MBlindPixelCalc Class. In the MBlindPixelCalc class there is also the
// function to replace the value of the noisy pixels with the interpolation
// of the content of the neighbors (SetUseInterpolation).
//
// Example:
// ...
// MBlindPixelCalc blind;
// blind.SetUseInterpolation();
// blind.SetUseBlindPixels();
//
// MImgCleanStd clean;
// ...
// tlist.AddToList(&blind);
// tlist.AddToList(&clean);
//
// Look at the MBlindPixelCalc Class for more details.
//
// Starting point: default values ----------------------------------------
//
// When an event is read, before the image cleaning, all the pixels that
// are in MCerPhotEvt are set as USED and NOT CORE. All the pixels belong
// to RING number 1 (like USED pixels).
// Look at MCerPhotPix.h to see how these informations of the pixel are
// stored.
// The default  cleaning METHOD is the STANDARD one and the number of the
// rings around the CORE pixel it analyzes is 1. Look at the Constructor
// of the class in MImgCleanStd.cc to see (or change) the default values.
//
// Example: To modify this setting, use the member functions
// SetMethod(MImgCleanStd::kDemocratic) and SetCleanRings(UShort_t n).
//
// MImgCleanStd:CleanStep1 -----------------------------------------------
//
// The first step of cleaning defines the CORE pixels. The CORE pixels are
// the ones which contain the informations about the core of the electro-
// magnetic shower.
// The ratio  (A_0/A_i) is calculated from fCam->GetPixRatio(i). A_0 is
// the area of the central pixel of the camera, A_i is the area of the
// examined pixel. In this way, if we have a MAGIC-like camera, with the
// outer pixels bigger than the inner ones, the level of cleaning in the
// two different regions is weighted.
// This avoids problems of deformations of the shower images.
// The signal S_i and the pedestal RMS Prms_i of the pixel are called from
// the object MCerPhotPix.
// If (default method = kStandard)
//Begin_Html
//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f1.png">
//End_Html
// the pixel is set as CORE pixel. L_1 (n=1) is called "first level of
// cleaning" (default: fCleanLvl1 = 3).
// All the other pixels are set as UNUSED and belong to RING 0.
// After this point, only the CORE pixels are set as USED, with RING
// number 1.
//
// MImgCleanStd:CleanStep2 ----------------------------------------------
//
// The second step of cleaning looks at the isolated CORE pixels and sets
// them to UNUSED. An isolated pixel is a pixel without CORE neighbors.
// At the end of this point, we have set as USED only CORE pixels with at
// least one CORE neighbor.
//
// MImgCleanStd:CleanStep3  ----------------------------------------------
//
// The third step of cleaning looks at all the pixels (USED or UNUSED) that
// surround the USED pixels.
// If the content of the analyzed pixel survives at the second level of
// cleaning, i.e. if
//Begin_Html
//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f1.png">
//End_Html
// the pixel is set as USED. L_2 (n=2) is called "second level of cleaning"
// (default:fCleanLvl2 = 2.5).
//
// When the number of RINGS to analyze is 1 (default value), only the
// pixels that have a neighbor CORE pixel are analyzed.
//
// There is the option to decide the number of times you want to repeat
// this procedure (number of RINGS analyzed around the core pixels = n).
// Every time the level of cleaning is the same (fCleanLvl2) and the pixel
// will belong to ring r+1, 1 < r < n+1. This is described in
// MImgCleanStd:CleanStep4 .
//
// Dictionary and member functions ---------------------------------------
//
// Here there is the detailed description of the member functions and of
// the terms commonly used in the class.
//
// STANDARD CLEANING:
// =================
// This is the method used for the CT1 data analysis. It is the default
// method of the class.
// The number of photo-electrons of a pixel (S_i) is compared to the
// pedestal RMS of the pixel itself (Prms_i). To have the comparison to
// the same photon density for all the pixels, taking into account they
// can have different areas, we have to keep in mind that the number of
// photons that hit each pixel, goes linearly with the area of the pixel.
// The fluctuations of the LONS are proportional to sqrt(A_i), so when we
// compare S_i with Prms_i, only a factor  sqrt(A_0/A_i) is missing to
// have the same (N.photons/Area) threshold for all the pixels.
//
//              !!WARNING: if noise independent from the
//         pixel size (example: electronic noise) is introduced,
//         then the noise fluctuations are no longer proportional
//         to sqrt(A_i), and then the cut value (for a camera with
//         pixels of different sizes) resulting from the above
//         procedure  would not be proportional to pixel size as we 
//         intend. In that case, democratic cleaning is preferred.
//
// If
//Begin_Html
//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f1.png">
//End_Html
// the pixel survives the cleaning and it is set as CORE (when L_n is the
// first level of cleaning, fCleanLvl1) or USED (when L_n is the second
// level of cleaning, fCleanLvl2).
//
// Example:
//
// MImgCleanStd clean;
// //creates a default Cleaning object, with default setting
// ...
// tlist.AddToList(&clean);
// // add the image cleaning to the main task list
//
// DEMOCRATIC CLEANING:
// ===================
// You use this cleaning method when you want to compare the number of
// photo-electons of each pixel with the average pedestal RMS
// (fInnerNoise = fSgb->GetSigmabarInner()) of the inner pixels (for the
// MAGIC camera they are the smaller ones):
//Begin_Html
//&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<img src="images/MImgCleanStd-f2.png">
//End_Html
// In this case, the simple ratio (A_0/A_i) is used to weight the level of
// cleaning, because both the inner and the outer pixels (that in MAGIC
// have a different area) are compared to the same pedestal RMS, coming
// from the inner pixels.
// To calculate the average pedestal RMS of the inner pixels, you have to
// add to the main task list an object of type MSigmabarCalc before the
// MImgCleanStd object. To know how the calculation of fInnerNoise is done
// look at the MSigmabarCalc Class.
//
// Example:
//
// MSigmabarCalc   sbcalc;
// //creates an object that calcutates the average pedestal RMS
// MImgCleanStd clean;
// ...
// tlist.AddToList(&sbcalc);
// tlist.AddToList(&clean);
//
// Member Function:  SetMethod()
// ============================
// When you call the MImgCleanStd task, the default method is kStandard.
//
// If you want to switch to the kDemocratic method you have to
// call this member function.
//
// Example:
//
// MImgCleanStd clean;
// //creates a default Cleaning object, with default setting
//
// clean.SetMethod(MImgCleanStd::kDemocratic);
// //now the method of cleaning is changed to Democratic
//
// FIRST AND SECOND CLEANING LEVEL
// ===============================
// When you call the MImgCleanStd task, the default cleaning levels are
// fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you
// create the MImgCleanStd object.
//
// Example:
//
// MImgCleanStd clean(Float_t lvl1,Float_t lvl2);
// //creates a default cleaning object, but the cleaning levels are now
// //lvl1 and lvl2.
//
// RING NUMBER
// ===========
// The standard cleaning procedure is such that it looks for the
// informations of the boundary part of the shower only on the first
// neighbors of the CORE pixels.
// There is the possibility now to look not only at the firs neighbors
// (first ring),but also further away, around the CORE pixels. All the new
// pixels you can find with this method, are tested with the second level
// of cleaning and have to have at least an USED neighbor.
//
// They will be also set as USED and will be taken into account during the
// calculation of the image parameters.
// The only way to distinguish them from the other USED pixels, is the
// Ring number, that is bigger than 1.
//
// Example: You can decide how many rings you want to analyze using:
//
// MImgCleanStd clean;
// //creates a default cleaning object (default number of rings =1)
// clean.SetCleanRings(UShort_t r);
// //now it looks r times around the CORE pixels to find new pixels with
// //signal.
//
//
//  Input Containers:
//   MGeomCam
//   MCerPhotEvt
//   [MSigmabar]
//
//  Output Containers:
//   MCerPhotEvt
//
/////////////////////////////////////////////////////////////////////////////
#include "MImgCleanStd.h"

#include <stdlib.h>       // atof                                         
#include <fstream>        // ofstream, SavePrimitive

#include <TGFrame.h>      // TGFrame
#include <TGLabel.h>      // TGLabel
#include <TGTextEntry.h>  // TGTextEntry

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

#include "MParList.h"
#include "MSigmabar.h"
#include "MCameraData.h"

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

#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"

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

#include "MGGroupFrame.h" // MGGroupFrame

ClassImp(MImgCleanStd);

using namespace std;

enum {
    kImgCleanLvl1,
    kImgCleanLvl2
};

static const TString gsDefName  = "MImgCleanStd";
static const TString gsDefTitle = "Task to perform image cleaning";

// --------------------------------------------------------------------------
//
// Default constructor. Here you can specify the cleaning method and levels. 
// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
// above mean) are used.
// Here you can also specify how many rings around the core pixels you want 
// to analyze (with the fixed lvl2). The default value for "rings" is 1.
//
MImgCleanStd::MImgCleanStd(const Float_t lvl1, const Float_t lvl2,
                           const char *name, const char *title)
    : fSgb(NULL), fCleaningMethod(kStandard), fCleanLvl1(lvl1),
    fCleanLvl2(lvl2), fCleanRings(1)

{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

    Print();
}

// --------------------------------------------------------------------------
//
// The first step of cleaning defines the CORE pixels. All the other pixels
// are set as UNUSED and belong to RING 0.
// After this point, only the CORE pixels are set as USED, with RING
// number 1.
//
//  NT 28/04/2003: now the option to use the standard method or the
//  democratic method is implemented:
//
//  kStandard:   This method looks for all pixels with an entry (photons)
//               that is three times bigger than the noise of the pixel
//               (default: 3 sigma, clean level 1)
//
//  kDemocratic: this method looks for all pixels with an entry (photons)
//               that is n times bigger than the noise of the mean of the
//               inner pixels (default: 3 sigmabar, clean level 1)
//
//
void MImgCleanStd::CleanStep1()
{
    const Int_t entries = fEvt->GetNumPixels();
    const TArrayD &data = fData->GetData();

    //
    // 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];

        if (data[pix.GetPixId()] <= fCleanLvl1)
            pix.SetPixelUnused();
    }
}

// --------------------------------------------------------------------------
//
//  Check if the survived pixel have a neighbor, that also
//  survived, otherwise set pixel to unused (removes pixels without
//  neighbors).
//
//  Takes the maximum pixel id from CleanStep1 as an argument
//
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.
    //
    // We allocate the array ourself because the TArrays always do
    // range check which slows down the access to the array
    // by 25-50%
    //
    Byte_t *ispixused = new Byte_t[fCam->GetNumPixels()];
    memset(ispixused, 0, sizeof(Byte_t)*fCam->GetNumPixels());

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

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

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

        // check if pixel is in use, if not goto next pixel in list
        if (ispixused[idx] == 0)
            continue;
 
        // check for 'used' neighbors of this pixel
        const MGeomPix &gpix  = (*fCam)[idx];
        const Int_t     nnmax = gpix.GetNumNeighbors();

        Bool_t hasNeighbor = kFALSE;

        //loop on the neighbors to check if they are used
        for (Int_t j=0; j<nnmax; j++)
        {
            const Int_t idx2 = gpix.GetNeighbor(j);

            // when you find an used neighbor, break the loop
            if (ispixused[idx2] == 1)
            {
                hasNeighbor = kTRUE;
                break;
            }
        }

        if (hasNeighbor == kFALSE)
            pix.SetPixelUnused();
    }

    delete ispixused;

    //
    // 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.SetPixelCore();
    }
} 

void MImgCleanStd::CleanStep3b(MCerPhotPix &pix)
{
    const Int_t idx = 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)[idx];
    const Int_t nnmax = gpix.GetNumNeighbors();

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

        if (!fEvt->GetPixById(idx2) || !fEvt->IsPixelCore(idx2))
            continue;

        pix.SetPixelUsed();
        break;
    }
}

// --------------------------------------------------------------------------
//
//   NT: Add option "rings": default value = 1. 
//   Look n (n>1) times for the boundary pixels around the used pixels.
//   If a pixel has more than 2.5 (clean level 2.5) sigma,
//   it is declared as used.
//
//   If a value<2 for fCleanRings is used, no CleanStep4 is done.
//
void MImgCleanStd::CleanStep4(UShort_t r, MCerPhotPix &pix)
{
    //
    // check if the pixel's next neighbor is a used pixel.
    // if it is a used pixel set pixel state to: used,
    // and tell to which ring it belongs to.
    //
    const Int_t idx = pix.GetPixId();
    MGeomPix  &gpix  = (*fCam)[idx];

    const Int_t nnmax = gpix.GetNumNeighbors();

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

        MCerPhotPix *npix = fEvt->GetPixById(idx2);

        if (!npix || !npix->IsPixelUsed() || npix->GetRing()>r-1 )
            continue;

        pix.SetRing(r);
        break;
    }
}

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

    for (UShort_t r=1; r<fCleanRings+1; r++)
    {
        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.IsPixelCore())
                continue;

            if (data[pix.GetPixId()] <= fCleanLvl2)
                continue;

            if (r==1)
                CleanStep3b(pix);
            else
                CleanStep4(r, pix);
        }
    }
}

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

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

    if (fCleaningMethod == kDemocratic)
    {
        fSgb = (MSigmabar*)pList->FindObject(AddSerialNumber("MSigmabar"));
        if (!fSgb)
        {
            *fLog << dbginf << "MSigmabar not found... aborting." << endl;
            return kFALSE;
        }
    }
    else
    {
        fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
        if (!fPed)
        {
            *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
            return kFALSE;
        }
    }

    fData = (MCameraData*)pList->FindCreateObj(AddSerialNumber("MCameraData"));
    if (!fData)
        return kFALSE;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Cleans the image.
//
Int_t MImgCleanStd::Process()
{
    if (fSgb)
        fData->Calc(*fEvt, *fSgb, *fCam);
    else
        fData->Calc(*fEvt, *fPed, *fCam);

#ifdef DEBUG
    *fLog << all << "CleanStep 1" << endl;
#endif
    CleanStep1();
#ifdef DEBUG
    *fLog << all << "CleanStep 2" << endl;
#endif
    CleanStep2();
#ifdef DEBUG
    *fLog << all << "CleanStep 3" << endl;
#endif
    CleanStep3();
#ifdef DEBUG
    *fLog << all << "Done." << endl;
#endif

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Print descriptor and cleaning levels.
//
void MImgCleanStd::Print(Option_t *o) const
{
    *fLog << all << GetDescriptor() << " using ";
    switch (fCleaningMethod)
    {
    case kDemocratic:
        *fLog << "democratic";
        break;
    case kStandard:
        *fLog << "standard";
        break;
    }
    *fLog << " cleaning initialized with noise level " << fCleanLvl1 << " and " << fCleanLvl2;
    *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
}

// --------------------------------------------------------------------------
//
//  Create two text entry fields, one for each cleaning level and a
//  describing text line.
//
void MImgCleanStd::CreateGuiElements(MGGroupFrame *f)
{
  //
  // Create a frame for line 3 and 4 to be able
  // to align entry field and label in one line
  //
  TGHorizontalFrame *f1 = new TGHorizontalFrame(f, 0, 0);
  TGHorizontalFrame *f2 = new TGHorizontalFrame(f, 0, 0);
  
  /*
   * --> use with root >=3.02 <--
   *
   
   TGNumberEntry *fNumEntry1 = new TGNumberEntry(frame, 3.0, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);
   TGNumberEntry *fNumEntry2 = new TGNumberEntry(frame, 2.5, 2, M_NENT_LVL1, kNESRealOne, kNEANonNegative);

  */
  TGTextEntry *entry1 = new TGTextEntry(f1, "****", kImgCleanLvl1);
  TGTextEntry *entry2 = new TGTextEntry(f2, "****", kImgCleanLvl2);
  
  // --- doesn't work like expected (until root 3.02?) --- fNumEntry1->SetAlignment(kTextRight);
  // --- doesn't work like expected (until root 3.02?) --- fNumEntry2->SetAlignment(kTextRight);
  
  entry1->SetText("3.0");
  entry2->SetText("2.5");
  
  entry1->Associate(f);
  entry2->Associate(f);
  
  TGLabel *l1 = new TGLabel(f1, "Cleaning Level 1");
  TGLabel *l2 = new TGLabel(f2, "Cleaning Level 2");
  
  l1->SetTextJustify(kTextLeft);
  l2->SetTextJustify(kTextLeft);
  
  //
  // Align the text of the label centered, left in the row
  // with a left padding of 10
  //
  TGLayoutHints *laylabel = new TGLayoutHints(kLHintsCenterY|kLHintsLeft, 10);
  TGLayoutHints *layframe = new TGLayoutHints(kLHintsCenterY|kLHintsLeft,  5, 0, 10);
  
  //
  // Add one entry field and the corresponding label to each line
  //
  f1->AddFrame(entry1);
  f2->AddFrame(entry2);
  
  f1->AddFrame(l1, laylabel);
  f2->AddFrame(l2, laylabel);
  
  f->AddFrame(f1, layframe);
  f->AddFrame(f2, layframe);
  
  f->AddToList(entry1);
  f->AddToList(entry2);
  f->AddToList(l1);
  f->AddToList(l2);
  f->AddToList(laylabel);
  f->AddToList(layframe);
}

// --------------------------------------------------------------------------
//
//  Process the GUI Events comming from the two text entry fields.
//
Bool_t MImgCleanStd::ProcessMessage(Int_t msg, Int_t submsg, Long_t param1, Long_t param2)
{
    if (msg!=kC_TEXTENTRY || submsg!=kTE_ENTER)
        return kTRUE;

    TGTextEntry *txt = (TGTextEntry*)FindWidget(param1);

    if (!txt)
        return kTRUE;

    Float_t lvl = atof(txt->GetText());

    switch (param1)
    {
    case kImgCleanLvl1:
        fCleanLvl1 = lvl;
        *fLog << "Cleaning level 1 set to " << lvl << " sigma." << endl;
        return kTRUE;

    case kImgCleanLvl2:
        fCleanLvl2 = lvl;
        *fLog << "Cleaning level 2 set to " << lvl << " sigma." << endl;
        return kTRUE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Implementation of SavePrimitive. Used to write the call to a constructor
// to a macro. In the original root implementation it is used to write
// gui elements to a macro-file.
//
void MImgCleanStd::StreamPrimitive(ofstream &out) const
{
    out << "   MImgCleanStd " << GetUniqueName() << "(";
    out << fCleanLvl1 << ", " << fCleanLvl2;

    if (fName!=gsDefName || fTitle!=gsDefTitle)
    {
        out << ", \"" << fName << "\"";
        if (fTitle!=gsDefTitle)
            out << ", \"" << fTitle << "\"";
    }
    out << ");" << endl;

    if (fCleaningMethod!=kDemocratic)
        return;

    out << "   " << GetUniqueName() << ".SetMethod(MImgCleanStd::kDemocratic);" << endl;

    if (fCleanRings==1)
        return;

    out << "   " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl;
}