Index: trunk/MagicSoft/Mars/mbadpixels/BadPixelsLinkDef.h
===================================================================
--- trunk/MagicSoft/Mars/mbadpixels/BadPixelsLinkDef.h	(revision 3068)
+++ trunk/MagicSoft/Mars/mbadpixels/BadPixelsLinkDef.h	(revision 3126)
@@ -9,4 +9,5 @@
 
 #pragma link C++ class MBadPixelsCalc+;
+#pragma link C++ class MBadPixelsTreat+;
 #pragma link C++ class MBadPixelsMerge+;
 
Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc
===================================================================
--- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc	(revision 3126)
+++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.cc	(revision 3126)
@@ -0,0 +1,320 @@
+/* ======================================================================== *\
+!
+! *
+! * 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
+//   MBadPixelsCam
+//
+//  Output Containers:
+//   MCerPhotEvt
+//
+/////////////////////////////////////////////////////////////////////////////
+#include "MBadPixelsTreat.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 "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 << "No camera geometry available... 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.
+//
+//  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 : Sqr(pix->GetErrorPhot());
+        pedrms[i] = nucp ? 0 : Sqr(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]   *= Sqr(ratio);
+	pedrms[i] *= Sqr(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]   += Sqr(nratio*evtpix->GetErrorPhot());
+	    pedrms[i] += Sqr(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, kUseInterpolation) && fGeomCam)
+        Interpolate();
+    else
+        Unmap();
+
+    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;
+}
Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.h
===================================================================
--- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.h	(revision 3126)
+++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsTreat.h	(revision 3126)
@@ -0,0 +1,55 @@
+#ifndef MARS_MBadPixelsTreat
+#define MARS_MBadPixelsTreat
+
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+class MGeomCam;
+class MCerPhotEvt;
+class MPedPhotCam;
+class MBadPixelsCam;
+
+class MBadPixelsTreat : public MTask
+{
+private:
+    MGeomCam      *fGeomCam;   //!
+    MPedPhotCam   *fPedPhot;   //!
+    MCerPhotEvt   *fEvt;       //!
+    MBadPixelsCam *fBadPixels; //!
+
+    Byte_t fFlags;       // flag for the method which is used
+
+    enum
+    {
+        kUseInterpolation = BIT(1),
+        kUseCentralPixel  = BIT(2),
+    };
+
+    static Double_t Sqr(Double_t x) { return x*x; }
+
+    void  Interpolate() const;
+    void  Unmap() const;
+    void  StreamPrimitive(ofstream &out) const;
+
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
+
+public:
+    MBadPixelsTreat(const char *name=NULL, const char *title=NULL);
+
+    void SetUseInterpolation(Bool_t b=kTRUE)
+    {
+        b ? SETBIT(fFlags, kUseInterpolation) : CLRBIT(fFlags, kUseInterpolation);
+    }
+    void SetUseCentralPixel(Bool_t b=kTRUE)
+    {
+        b ? SETBIT(fFlags, kUseCentralPixel) : CLRBIT(fFlags, kUseCentralPixel);
+    }
+
+
+    ClassDef(MBadPixelsTreat, 1) // Task to treat bad pixels (interpolation, unmapping)
+}; 
+
+#endif
+
Index: trunk/MagicSoft/Mars/mbadpixels/Makefile
===================================================================
--- trunk/MagicSoft/Mars/mbadpixels/Makefile	(revision 3068)
+++ trunk/MagicSoft/Mars/mbadpixels/Makefile	(revision 3126)
@@ -37,4 +37,5 @@
            MBadPixelsMerge.cc \
            MBadPixelsCalc.cc \
+           MBadPixelsTreat.cc \
            MMcBadPixelsSet.cc
 
