source: trunk/MagicSoft/Mars/mimage/MImgCleanTGB.cc@ 2360

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

/////////////////////////////////////////////////////////////////////////////
//
//  MImgCleanTGB
//
// 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();
//
// MImgCleanTGB 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 MImgCleanTGB.cc to see (or change) the default values.
//
// Example: To modify this setting, use the member functions
// SetMethod(MImgCleanTGB::kDemocratic) and SetCleanRings(UShort_t n).
//
// MImgCleanTGB: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/MImgCleanTGB-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.
//
// MImgCleanTGB: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.
//
// MImgCleanTGB: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/MImgCleanTGB-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
// MImgCleanTGB: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/MImgCleanTGB-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:
//
// MImgCleanTGB 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/MImgCleanTGB-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
// MImgCleanTGB 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
// MImgCleanTGB clean;
// ...
// tlist.AddToList(&sbcalc);
// tlist.AddToList(&clean);
//
// Member Function:  SetMethod()
// ============================
// When you call the MImgCleanTGB task, the default method is kStandard.
//
// If you want to switch to the kDemocratic method you have to
// call this member function.
//
// Example:
//
// MImgCleanTGB clean;
// //creates a default Cleaning object, with default setting
//
// clean.SetMethod(MImgCleanTGB::kDemocratic);
// //now the method of cleaning is changed to Democratic
//
// FIRST AND SECOND CLEANING LEVEL
// ===============================
// When you call the MImgCleanTGB task, the default cleaning levels are
// fCleanLvl1 = 3, fCleanLvl2 = 2.5. You can change them easily when you
// create the MImgCleanTGB object.
//
// Example:
//
// MImgCleanTGB 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:
//
// MImgCleanTGB 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 "MImgCleanTGB.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 "MGeomPix.h"
#include "MGeomCam.h"

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

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

#include "MGGroupFrame.h" // MGGroupFrame

ClassImp(MImgCleanTGB);

using namespace std;

enum {
    kImgCleanLvl1,
    kImgCleanLvl2
};

static const TString gsDefName  = "MImgCleanTGB";
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.
//
MImgCleanTGB::MImgCleanTGB(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();
}


Int_t MImgCleanTGB::CleanStep3b(MCerPhotPix &pix)
{
    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();

    Int_t rc = 0;

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

        if (fEvt->GetPixById(id2) && fEvt->IsPixelUsed(id2))
            rc++;
    }
    return rc;
}

// --------------------------------------------------------------------------
//
//   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 MImgCleanTGB::CleanStep3(Int_t num1, Int_t num2)
{
    const Int_t entries = fEvt->GetNumPixels();

    Int_t *u = new Int_t[entries];

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

    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];
        if (u[i]<num1)
            pix.SetPixelUnused();
        if (u[i]>num2)
            pix.SetPixelUsed();
    }

    delete u;
}

void MImgCleanTGB::CleanStep3(Byte_t *nb, Int_t num1, Int_t num2)
{
    const Int_t entries = fEvt->GetNumPixels();

    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];

        const Int_t idx = pix.GetPixId();

        if (nb[idx]<num1 && pix.IsPixelUsed())
        {
            const MGeomPix &gpix = (*fCam)[idx];
            const Int_t nnmax = gpix.GetNumNeighbors();
            for (Int_t j=0; j<nnmax; j++)
                nb[gpix.GetNeighbor(j)]--;

            pix.SetPixelUnused();
        }
    }

    for (Int_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];

        const Int_t idx = pix.GetPixId();
        if (nb[idx]>num2 && !pix.IsPixelUsed())
        {
            const MGeomPix &gpix = (*fCam)[idx];
            const Int_t nnmax = gpix.GetNumNeighbors();
            for (Int_t j=0; j<nnmax; j++)
                nb[gpix.GetNeighbor(j)]++;

            pix.SetPixelUsed();
        }
    }
}

// --------------------------------------------------------------------------
//
//  Check if MEvtHeader exists in the Parameter list already.
//  if not create one and add them to the list
//
Int_t MImgCleanTGB::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;
    }

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Cleans the image.
//
Int_t MImgCleanTGB::Process()
{
    const Int_t entries = fEvt->GetNumPixels();

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

        const Double_t entry  = pix.GetNumPhotons();
        const Double_t factor = fCam->GetPixRatio(pix.GetPixId());
        const Float_t  noise  = (*fPed)[pix.GetPixId()].GetPedestalRms();

        if (entry * TMath::Sqrt(factor) <= fCleanLvl2 * noise)
        {
            sum += entry*factor;
            sq  += entry*entry*factor*factor;
            w   += factor;
            w2  += factor*factor;
        }
    }

    Double_t mean = sum/w;
    Double_t sdev = sqrt(sq/w2 - mean*mean);

    Byte_t *nb = new Byte_t[1000];
    memset(nb, 0, 577);

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

        const Double_t entry  = pix.GetNumPhotons();
        const Double_t factor = fCam->GetPixRatio(idx);

        if (entry*factor > fCleanLvl1*sdev)
        {
            pix.SetPixelUsed();

            const MGeomPix &gpix = (*fCam)[idx];
            const Int_t nnmax = gpix.GetNumNeighbors();
            for (Int_t j=0; j<nnmax; j++)
                nb[gpix.GetNeighbor(j)]++;
        }
        else
            pix.SetPixelUnused();
    }

    CleanStep3(nb, 2, 3);
    //CleanStep3(nb, 2, 3);
    //CleanStep3(nb, 2, 3);

    delete nb;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Print descriptor and cleaning levels.
//
void MImgCleanTGB::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 MImgCleanTGB::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 MImgCleanTGB::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 MImgCleanTGB::StreamPrimitive(ofstream &out) const
{
    out << "   MImgCleanTGB " << 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(MImgCleanTGB::kDemocratic);" << endl;

    if (fCleanRings==1)
        return;

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