source: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc@ 3481

/* ======================================================================== *\
!
! *
! * 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): Oscar Blanch    12/2001 <mailto:blanch@ifae.es>
!   Author(s): Thomas Bretz    08/2002 <mailto:tbretz@astro.uni.wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MBadPixelsTreat
//
//  You can use MBadPixelsTreat::SetUseInterpolation to replaced the
//  bad pixels by the average of its neighbors instead of unmapping
//  them. If you want to include the central pixel use
//  MBadPixelsTreat::SetUseCentralPixel. The bad pixels are taken from
//  an existing MBadPixelsCam.
//
//  Input Containers:
//   MCerPhotEvt
//   MPedPhotCam
//   MBadPixelsCam
//   [MGeomCam]
//
//  Output Containers:
//   MCerPhotEvt
//
/////////////////////////////////////////////////////////////////////////////
#include "MBadPixelsTreat.h"

#include <fstream>

#include <TArrayD.h>

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

#include "MParList.h"

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

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

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

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

ClassImp(MBadPixelsTreat);

using namespace std;

static const TString gsDefName  = "MBadPixelsTreat";
static const TString gsDefTitle = "Task to treat bad pixels (interpolation, unmapping)";

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

// --------------------------------------------------------------------------
//
//  - Try to find or create MBlindPixels in parameter list.
//  - get the MCerPhotEvt from the parlist (abort if missing)
//  - if no pixels are given by the user try to determin the starfield
//    from the monte carlo run header.
//
Int_t MBadPixelsTreat::PreProcess (MParList *pList)
{
    fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
    if (!fBadPixels)
    {
        *fLog << err << "MBadPixelsCam not found... aborting." << endl;
        return kFALSE;
    }

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

    fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!fGeomCam && TESTBIT(fFlags, kUseInterpolation))
    {
        *fLog << err << "MGeomCam not found... can't use interpolation." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Replaces each pixel (signal, signal error, pedestal, pedestal rms)
//  by the average of its surrounding pixels.
//  If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
//  included.
//
void MBadPixelsTreat::InterpolateSignal() const
{
    const UShort_t entries = fGeomCam->GetNumPixels();

    //
    // Create arrays
    //
    TArrayD nphot(entries);
    TArrayD perr(entries);
 
    //
    // Loop over all pixels
    //
    for (UShort_t i=0; i<entries; i++)
    {
        //
        // Check whether pixel with idx i is blind
        //
        if (!(*fBadPixels)[i].IsBad())
            continue;

        //
        // Get a pointer to this pixel. If it is not yet existing
        // create a new entry for this pixel in MCerPhotEvt
        //
        MCerPhotPix *pix = fEvt->GetPixById(i);
        if (!pix)
            pix = fEvt->AddPixel(i, 0, 0);

        //
        // Get the corresponding geometry and pedestal
        //
        const MGeomPix &gpix = (*fGeomCam)[i];

        // Do Not-Use-Central-Pixel
        const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);

        Int_t num = nucp ? 0 : 1;

        nphot[i]  = nucp ? 0 : pix->GetNumPhotons();
        perr[i]   = nucp ? 0 : Pow2(pix->GetErrorPhot());

        //
        // The values are rescaled to the small pixels area for the right comparison
        //
        const Double_t ratio = fGeomCam->GetPixRatio(i);

        nphot[i] *= ratio;
        perr[i]  *= ratio;

        //
        // Loop over all its neighbors
        //
        const Int_t n = gpix.GetNumNeighbors();
        for (int j=0; j<n; j++)
        {
            const UShort_t nidx = gpix.GetNeighbor(j);

            //
            // Do not use blind neighbors
            //
            if (!(*fBadPixels)[i].IsBad())
                continue;

            //
            // Check whether the neighbor has a signal stored
            //
            const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
            if (!evtpix)
                continue;

            //
            // Get the geometry for the neighbor
            //
            const Double_t nratio = fGeomCam->GetPixRatio(nidx);

            //
            //The error is calculated as the quadratic sum of the errors
            //
            nphot[i] += nratio*evtpix->GetNumPhotons();
            perr[i]  += nratio* Pow2(evtpix->GetErrorPhot());

            num++;
        }

        //
        // Now the mean is calculated and the values rescaled back
        // to the pixel area
        //
        nphot[i] /= (num*ratio);
        perr[i]   = TMath::Sqrt(perr[i]/(num*ratio));
 
        // Check if there are enough neighbors to calculate the mean
        // If not, unmap the pixel. The maximum number of blind neighbors
        // should be 2
        if (num<3)
        {
            pix->SetPixelUnmapped();
            continue;
        }

        pix->Set(nphot[i], perr[i]);
    }
}

// --------------------------------------------------------------------------
//
void MBadPixelsTreat::InterpolatePedestals() const
{
    const Int_t entries = fPedPhot->GetSize();

    TArrayD ped(entries);
    TArrayD rms(entries);

    //
    // Loop over all pixels
    //
    for (UShort_t i=0; i<entries; i++)
    {
        //
        // Check whether pixel with idx i is blind
        //
        if (!(*fBadPixels)[i].IsBad())
            continue;

        //
        // Get the corresponding geometry and pedestal
        //
        const MGeomPix    &gpix = (*fGeomCam)[i];
        const MPedPhotPix &ppix = (*fPedPhot)[i];

        // Do Not-Use-Central-Pixel
        const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);

        Int_t num = nucp ? 0 : 1;

        ped[i] = nucp ? 0 : ppix.GetMean();
        rms[i] = nucp ? 0 : Pow2(ppix.GetRms());

        //
        // The values are rescaled to the small pixels area for the right comparison
        //
        const Double_t ratio = fGeomCam->GetPixRatio(i);

        ped[i] *= ratio;
        rms[i] *= ratio;

        //
        // Loop over all its neighbors
        //
        const Int_t n = gpix.GetNumNeighbors();
        for (int j=0; j<n; j++)
        {
            const UShort_t nidx = gpix.GetNeighbor(j);

            //
            // Do not use blind neighbors
            //
            if (!(*fBadPixels)[i].IsBad())
                continue;

            //
            // Get the geometry for the neighbor
            //
            const Double_t    nratio = fGeomCam->GetPixRatio(nidx);
            const MPedPhotPix &nppix = (*fPedPhot)[nidx];

            //
            //The error is calculated as the quadratic sum of the errors
            //
            ped[i] += (nratio*nppix.GetMean());
            rms[i] += nratio*Pow2(nppix.GetRms());

            num++;
        }

        // Check if there are enough neighbors to calculate the mean
        // If not, unmap the pixel. The minimum number of good neighbors
        // should be 3
        if (num < 3)
        {
            MCerPhotPix *pix =fEvt->GetPixById(i);
            if (!pix)
                pix = fEvt->AddPixel(i, 0, 0);
            pix->SetPixelUnmapped();
            continue;
        }

        //
        // Now the mean is calculated and the values rescaled back
        // to the pixel area
        //
        ped[i] /=  (num*ratio);
        rms[i]  = TMath::Sqrt(rms[i]/(num*ratio));

        (*fPedPhot)[i].Set(ped[i], rms[i]);
    }
}

// --------------------------------------------------------------------------
//
//  Replaces each pixel (signal, signal error, pedestal, pedestal rms)
//  by the average of its surrounding pixels.
//  If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
//  included.
//
//  NT: Informations about the interpolated pedestals are added. 
//      When the option Interpolated is called, the blind pixel with the new
//      values of signal and fluttuation is included in the calculation of
//      the Image Parameters.
//
void MBadPixelsTreat::Interpolate() const
{
    const UShort_t entries = fGeomCam->GetNumPixels();

    //
    // Create arrays
    //
    Double_t *nphot  = new Double_t[entries];
    Double_t *perr   = new Double_t[entries];
    Double_t *ped    = new Double_t[entries];
    Double_t *pedrms = new Double_t[entries];
 
    //
    // Loop over all pixels
    //
    for (UShort_t i=0; i<entries; i++)
    {
        //
        // Check whether pixel with idx i is blind
        //
        if ((*fBadPixels)[i].IsBad())
            continue;

        //
        // Get a pointer to this pixel. If it is not yet existing
        // create a new entry for this pixel in MCerPhotEvt
        //
        const MCerPhotPix *pix = fEvt->GetPixById(i);
        if (!pix)
            pix = fEvt->AddPixel(i, 0, 0);

        //
        // Get the corresponding geometry and pedestal
        //
        const MGeomPix &gpix    = (*fGeomCam)[i];
        const MPedPhotPix &ppix = (*fPedPhot)[i];

        // Do Not-Use-Central-Pixel
        const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);

        Int_t num = nucp ? 0 : 1;

        nphot[i]  = nucp ? 0 : pix->GetNumPhotons();
        ped[i]    = nucp ? 0 : ppix.GetMean();
        perr[i]   = nucp ? 0 : Pow2(pix->GetErrorPhot());
        pedrms[i] = nucp ? 0 : Pow2(ppix.GetRms());

        //
        // The values are rescaled to the small pixels area for the right comparison
        //
        const Double_t ratio = fGeomCam->GetPixRatio(i);

        nphot[i]  *= ratio;
        ped[i]    *= ratio;
        perr[i]   *= Pow2(ratio);
        pedrms[i] *= Pow2(pedrms[i]);

        //
        // Loop over all its neighbors
        //
        const Int_t n = gpix.GetNumNeighbors();
        for (int j=0; j<n; j++)
        {
            const UShort_t nidx = gpix.GetNeighbor(j);

            //
            // Do not use blind neighbors
            //
            if ((*fBadPixels)[i].IsBad())
                continue;

            //
            // Check whether the neighbor has a signal stored
            //
            const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
            if (!evtpix)
                continue;

            //
            // Get the geometry for the neighbor
            //
            const Double_t nratio = fGeomCam->GetPixRatio(nidx);
            MPedPhotPix &nppix    = (*fPedPhot)[nidx];

            //
            //The error is calculated as the quadratic sum of the errors
            //
            ped[i]    += nratio*nppix.GetMean();
            nphot[i]  += nratio*evtpix->GetNumPhotons();
            perr[i]   += Pow2(nratio*evtpix->GetErrorPhot());
            pedrms[i] += Pow2(nratio*nppix.GetRms());

            num++;
        }

        //
        // Now the mean is calculated and the values rescaled back to the pixel area
        //
        nphot[i] /= num*ratio;
        ped[i]   /= num*ratio;
        perr[i]   = TMath::Sqrt(perr[i]/num)*ratio;
        pedrms[i] = TMath::Sqrt(pedrms[i]/num)*ratio;

    }

    //
    // Now the new pixel values are calculated and can be replaced in
    // the corresponding containers
    //
    for (UShort_t i=0; i<entries; i++)
    {
        //
        // Do not use blind neighbors
        //
        if ((*fBadPixels)[i].IsBad())
            continue;

        //
        // It must exist, we have created it in the loop before.
        //
        fEvt->GetPixById(i)->Set(nphot[i], perr[i]);
        (*fPedPhot)[i].Set(ped[i], pedrms[i]);
    }

    //
    // Delete the intermediat arrays
    //
    delete nphot;
    delete perr;
    delete ped;
    delete pedrms;
}

// --------------------------------------------------------------------------
//
//  Removes all blind pixels from the analysis by setting their state
//  to unused.
//
void MBadPixelsTreat::Unmap() const
{
    const UShort_t entries = fEvt->GetNumPixels();

    //
    // remove the pixels in fPixelsIdx if they are set to be used,
    // (set them to 'unused' state)
    //
    for (UShort_t i=0; i<entries; i++)
    {
        MCerPhotPix &pix = (*fEvt)[i];

        if ((*fBadPixels)[pix.GetPixId()].IsBad())
            pix.SetPixelUnused();
    }
}

// --------------------------------------------------------------------------
//
// Treat the blind pixels
//
Int_t MBadPixelsTreat::Process()
{
    /*
    if (TESTBIT(fFlags, kCheckPedestalRms))
    {
        // if the number of blind pixels is too high, do not interpolate
       if (CheckPedestalRms()==kFALSE)
           return kTRUE;

       if (TESTBIT(fFlags, kUseInterpolation))
           InterpolatePedestals();
    }
    */

    if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
        Interpolate();
    else
        Unmap();


    //fErrors[0]++;

    return kTRUE;
}

void MBadPixelsTreat::StreamPrimitive(ofstream &out) const
{
    out << "   MBadPixelsTreat " << GetUniqueName();
    if (fName!=gsDefName || fTitle!=gsDefTitle)
    {
        out << "(\"" << fName << "\"";
        if (fTitle!=gsDefTitle)
            out << ", \"" << fTitle << "\"";
        out <<")";
    }
    out << ";" << endl;

    if (TESTBIT(fFlags, kUseInterpolation))
        out << "   " << GetUniqueName() << ".SetUseInterpolation();" << endl;
    if (TESTBIT(fFlags, kUseCentralPixel))
        out << "   " << GetUniqueName() << ".SetUseCentralPixel();" << endl;
}