source: trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc@ 2282

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

/////////////////////////////////////////////////////////////////////////////
//
//  MBlindPixelCalc
//
//  This is the specific image cleaning for a list of pixels. This task
//  sets the pixels listed in fPixelsIdx to unused so they should not be
//  used for analysis (eg calculation of hillas parameters).
//
//  You can use MBlindPixelCalc::SetUseInterpolation to replaced the
//  blind pixels by the average of its neighbors instead of unmapping
//  them. If you want to include the central pixel use
//  MBlindPixelCalc::SetUseCentralPixel.
//
//  You have three options:
//   1) Call SetUseBlindPixels():
//      This will take an existing MBlindPixels container filled from
//      elsewhere (eg. MCT1ReadPreProc) and use this pixels as blind
//      pixels.
//   2) Call SetPixels():
//      This will setup an array with pixel numbers. These pixels are used
//      as blind pixels.
//   3) Neither 1) nor 2)
//      This options tries to identify the starfield from the
//      MMcRunHeader container and tries to identifies it. If it is known
//      (eg. Crab) the fixed build in pixel numbers are used as blind
//      pixels.
//
//  If neither an array of pixels is given (or its size is 0) and
//  MMcRunHeader couldn't be found the task removes itself from the
//  tasklist.
//
//  Implemented star fields:
//   - Crab: 400, 401, 402, 437, 438, 439
//
//  Input Containers:
//   MCerPhotEvt[, MBlindPixels]
//
//  Output Containers:
//   MCerPhotEvt, MBlindPixels
//
/////////////////////////////////////////////////////////////////////////////
#include "MBlindPixelCalc.h"

#include <fstream>

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

#include "MParList.h"

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

#include "MMcRunHeader.hxx"

ClassImp(MBlindPixelCalc);

using namespace std;

static const TString gsDefName  = "MBlindPixelCalc";
static const TString gsDefTitle = "Find hot spots (star, broken pixels, etc)";

// --------------------------------------------------------------------------
//
// Default constructor.
//
MBlindPixelCalc::MBlindPixelCalc(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 MBlindPixelCalc::PreProcess (MParList *pList)
{
    if (TESTBIT(fFlags, kUseBlindPixels))
        fPixels = (MBlindPixels*)pList->FindObject("MBlindPixels");
    else
        fPixels = (MBlindPixels*)pList->FindCreateObj("MBlindPixels");
    if (!fPixels)
        return kFALSE; 
    
    fEvt = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
    if (!fEvt)
    {
        *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
        return kFALSE;
    }

    fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fGeomCam)
        *fLog << warn << dbginf << "No camera geometry available... can't use interpolation." << endl;

    if (TESTBIT(fFlags, kUseBlindPixels))
        return kTRUE;

    const UShort_t size = fPixelsIdx.GetSize();

    if (size == 0)
    {
        if (!pList->FindObject("MMcRunHeader"))
        {
            *fLog << warn << "Warning - Neither blind pixels are given nor a MMcRunHeader was found... removing MBlindPixelCalc from list." << endl;
            return kSKIP;
        }
        return kTRUE;
    }

    // Set as blind pixels the global blind pixels, which are given
    // through the macros

    UShort_t numids = fPixelsIdx.GetSize();

    for(Int_t i = 0; i<numids; i++)
        fPixels->SetPixelBlind(fPixelsIdx[i]);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  Replaces each pixel by the average of its surrounding pixels.
//  If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
//  included.
//
void MBlindPixelCalc::Interpolate() const
{
    const UShort_t entries = fEvt->GetNumPixels();

    Double_t *nphot = new Double_t[entries];
    Double_t *perr  = new Double_t[entries];

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

        const Int_t idx = pix.GetPixId();

        if (!fPixels->IsBlind(idx))
            continue;

        const MGeomPix &gpix = (*fGeomCam)[idx];

        Int_t num = TESTBIT(fFlags, kUseCentralPixel) ? 1 : 0;

        nphot[i] = TESTBIT(fFlags, kUseCentralPixel) ? pix.GetNumPhotons() : 0;
        perr[i]  = TESTBIT(fFlags, kUseCentralPixel) ? pix.GetErrorPhot()  : 0;

        nphot[i] *= fGeomCam->GetPixRatio(idx);
        // FIXME: perr[i] ???

        const Int_t n = gpix.GetNumNeighbors();
        for (int j=0; j<n; j++)
        {
            const UShort_t nidx = gpix.GetNeighbor(j);

            if (fPixels->IsBlind(nidx))
                continue;

            const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
            if (evtpix)
            {
                nphot[i] += evtpix->GetNumPhotons()*fGeomCam->GetPixRatio(nidx);
                perr[i]  += evtpix->GetErrorPhot();
                // FIXME: perr[i] ???
            }
            num++;
        }

        nphot[i] /= num*fGeomCam->GetPixRatio(idx);
        perr[i]  /= num/*FIXME:???*/;
    }

    if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
        for (UShort_t i=0; i<entries; i++)
        {
            MCerPhotPix &pix = (*fEvt)[i];

            if (fPixels->IsBlind(pix.GetPixId()))
                pix.Set(nphot[i], perr[i]);
        }

    delete nphot;
    delete perr;
}

// --------------------------------------------------------------------------
//
//  Removes all blind pixels from the analysis by setting their state
//  to unused.
//
void MBlindPixelCalc::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 (fPixels->IsBlind(pix.GetPixId()))
            pix.SetPixelUnused();
    }
}

// --------------------------------------------------------------------------
//
// Treat the blind pixels
//
Int_t MBlindPixelCalc::Process()
{
    if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
        Interpolate();
    else
        Unmap();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
//  - Check whether pixels to disable are available. If pixels are
//    given by the user nothing more is done.
//  - Otherwise try to determin the blind pixels from the starfield
//    given in MMcRunHeader.
//
Bool_t MBlindPixelCalc::ReInit(MParList *pList)
{
    if (TESTBIT(fFlags, kUseBlindPixels))
        return kTRUE;

    //
    // If pixels are given by the user, we are already done
    //
    if (fPixelsIdx.GetSize() > 0)
        return kTRUE;

    //
    // Delete the old array holding the blind pixels for the last file
    //
    fPixels->Clear();

    //
    // Set as blind some particular pixels because of a particular
    // Star Field of View.
    //
    MMcRunHeader *mcrun = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
    if (!mcrun)
    {
        *fLog << warn << "MBlindPixelCalc::ReInit: Warning - No run header available... no action." << endl;
        return kTRUE;
    }

    Int_t rah, ram, ras;
    Int_t ded, dem, des;
    mcrun->GetStarFieldRa(&rah, &ram, &ras);
    mcrun->GetStarFieldDec(&ded, &dem, &des);

    if (rah!=5 || ram!=34 || ras!=32 || ded!=22 || dem!=0 || des!=55)
    {
        *fLog << warn << "Warning - Starfield unknown..." << endl;
        return kTRUE;
    }

    //
    // Case for Crab Nebula FOV
    //
    fPixels->SetPixelBlind(400);
    fPixels->SetPixelBlind(401);
    fPixels->SetPixelBlind(402);
    fPixels->SetPixelBlind(437);
    fPixels->SetPixelBlind(438);
    fPixels->SetPixelBlind(439);

    *fLog << inf;
    *fLog << "FOV is centered at CRAB NEBULA: Setting 6 blind pixels" << endl;
    *fLog << "to avoid bias values of analysis due to CRAB NEBULA:" << endl;
    *fLog << "   Pixels: 400, 401, 402, 437, 438, 439" << endl;

    return kTRUE;
}

void MBlindPixelCalc::StreamPrimitive(ofstream &out) const
{
    out << "   MBlindPixelCalc " << 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;

    if (fPixelsIdx.GetSize()==0)
        return;

    out << "   {" << endl;
    out << "      TArrayS dummy;" << endl;
    for (int i=0; i<fPixelsIdx.GetSize(); i++)
        out << "      dummy[" << i << "]=" << ((TArrayS)fPixelsIdx)[i] << ";" << endl;
    out << "      " << GetUniqueName() << ".SetPixels(dummy);" << endl;
    out << "   }" << endl;
}