Index: /trunk/Mars/mimage/MImgCleanSimple.cc
===================================================================
--- /trunk/Mars/mimage/MImgCleanSimple.cc	(revision 14202)
+++ /trunk/Mars/mimage/MImgCleanSimple.cc	(revision 14202)
@@ -0,0 +1,411 @@
+/* ======================================================================== *\
+!
+! *
+! * 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, 06/2012 <mailto:thomas.bretz@epfl.ch>
+!
+!   Copyright: MAGIC Software Development, 2000-2012
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//
+//  MImgCleanSimple
+//
+/////////////////////////////////////////////////////////////////////////////
+#include "MImgCleanSimple.h"
+
+#include <TEnv.h>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MParList.h"
+
+#include "MGeomPix.h"
+#include "MGeomCam.h"
+
+#include "MSignalPix.h"
+#include "MSignalCam.h"
+
+ClassImp(MImgCleanSimple);
+
+using namespace std;
+
+static const TString gsDefName  = "MImgCleanSimple";
+static const TString gsDefTitle = "Task to perform image cleaning";
+
+const TString MImgCleanSimple::gsNameSignalCam ="MSignalCam";  // default name of the 'MSignalCam' container
+const TString MImgCleanSimple::gsNameGeomCam   ="MGeomCam";    // default name of the 'MGeomCam' container
+
+// --------------------------------------------------------------------------
+//
+// 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.
+//
+MImgCleanSimple::MImgCleanSimple(const Float_t lvl1, const Float_t lvl2,
+                                 const char *name, const char *title) : fCleaningMethod(kStandard),
+    fCleanLvl1(lvl1), fCleanLvl2(-1), fTimeLvl1(lvl2), fTimeLvl2(-1),
+    fNameGeomCam(gsNameGeomCam), fNameSignalCam(gsNameSignalCam)
+{
+    fName  = name  ? name  : gsDefName.Data();
+    fTitle = title ? title : gsDefTitle.Data();
+}
+
+void MImgCleanSimple::ResetCleaning() const
+{
+    //
+    // check the number of all pixels against the noise level and
+    // set them to 'unused' state if necessary
+    //
+    const UInt_t npixevt = fEvt->GetNumPixels();
+    for (UInt_t idx=0; idx<npixevt; idx++)
+    {
+        MSignalPix &pix = (*fEvt)[idx];
+        if (pix.IsPixelUnmapped())
+            continue;
+
+        pix.SetPixelUnused();
+        pix.SetPixelCore(kFALSE);
+    }
+}
+
+void MImgCleanSimple::DoCleanStd() const
+{
+    const UInt_t npixevt = fEvt->GetNumPixels();
+    for (UInt_t idx=0; idx<npixevt; idx++)
+    {
+        MSignalPix &pix = (*fEvt)[idx];
+
+        // Check if this is above the cleaning level
+        if (pix.GetNumPhotons() <= fCleanLvl1)
+            continue;
+
+        // Ignore unmapped pixels
+        if (pix.IsPixelUnmapped())
+            continue;
+
+        const MGeom &gpix = (*fCam)[idx];
+
+        const Int_t n = gpix.GetNumNeighbors();
+        for (Int_t i=0; i<n; i++)
+        {
+            const MSignalPix &npix = (*fEvt)[gpix.GetNeighbor(i)];
+
+            // Check if one neighbor pixel is also above the cleaning level
+            if (npix.GetNumPhotons() <= fCleanLvl1)
+                continue;
+
+            // Check if both pixels are within the given time window
+            if (TMath::Abs(npix.GetArrivalTime()-pix.GetArrivalTime())>=fTimeLvl1)
+                continue;
+
+            // Ignore unmapped pixels
+            if (npix.IsPixelUnmapped())
+                continue;
+
+            // The pixel is now marked to survive the cleaning. It's neighbor will
+            // be marked later, because the condition is symmetric and the neighbor
+            // will fullfil the condition, too.
+            pix.SetPixelUsed();
+            pix.SetPixelCore();
+
+            break;
+        }
+    }
+}
+
+void MImgCleanSimple::DoCleanSum2() const
+{
+    const UInt_t npixevt = fEvt->GetNumPixels();
+    for (UInt_t idx=0; idx<npixevt; idx++)
+    {
+        MSignalPix &pix = (*fEvt)[idx];
+
+        // Ignore unmapped pixels
+        if (pix.IsPixelUnmapped())
+            continue;
+
+        const MGeom &gpix = (*fCam)[idx];
+
+        const Int_t n = gpix.GetNumNeighbors();
+        for (Int_t i=0; i<n; i++)
+        {
+            const MSignalPix &npix = (*fEvt)[gpix.GetNeighbor(i)];
+
+            // Check if the sum (or average) charge of the two neighbors is below threshold
+            if (pix.GetNumPhotons()+npix.GetNumPhotons()<=fCleanLvl1*2)
+                continue;
+
+            // Check if both pixels are within the given time window
+            if (TMath::Abs(npix.GetArrivalTime()-pix.GetArrivalTime())>=fTimeLvl1)
+                continue;
+
+            // Ignore unmapped pixels
+            if (npix.IsPixelUnmapped())
+                continue;
+
+            // The pixel is now marked to survive the cleaning. It's neighbor will
+            // be marked later, because the condition is symmetric and the neighbor
+            // will fullfil the condition, too.
+            pix.SetPixelUsed();
+            pix.SetPixelCore();
+            break;
+        }
+    }
+}
+
+void MImgCleanSimple::DoCleanSum3() const
+{
+    const Double_t clvl = fCleanLvl2>0 ? fCleanLvl2*3 : fCleanLvl1*3;
+    const Double_t tlvl = fTimeLvl2>0  ? fTimeLvl2    : fTimeLvl1;
+
+    const UInt_t npixevt = fEvt->GetNumPixels();
+    for (UInt_t idx=0; idx<npixevt; idx++)
+    {
+        MSignalPix &pix = (*fEvt)[idx];
+
+        // Ignore unmapped pixels
+        if (pix.IsPixelUnmapped())
+            continue;
+
+        const MGeom &gpix = (*fCam)[idx];
+
+        const Int_t n = gpix.GetNumNeighbors();
+        for (Int_t i=0; i<n; i++)
+        {
+            const MSignalPix &npix1 = (*fEvt)[gpix.GetNeighbor(i)];
+            const MSignalPix &npix2 = (*fEvt)[gpix.GetNeighbor((i+1)%n)];
+
+            // Check if the sum (or average) charge of the three neighbors (closed package)
+            // is below threshold
+            if (pix.GetNumPhotons()+npix1.GetNumPhotons()+npix2.GetNumPhotons()<=clvl)
+                continue;
+
+            // Check if all thre pixels fullfil the time window condition
+            if (TMath::Abs(pix.GetArrivalTime()-npix1.GetArrivalTime())>=tlvl)
+                continue;
+            if (TMath::Abs(npix1.GetArrivalTime()-npix2.GetArrivalTime())>=tlvl)
+                continue;
+            if (TMath::Abs(npix2.GetArrivalTime()-pix.GetArrivalTime())>=tlvl)
+                continue;
+
+            // Ignore unmapped pixels
+            if (npix1.IsPixelUnmapped() || npix2.IsPixelUnmapped())
+                continue;
+
+            // The pixel is now marked to survive the cleaning. It's neighbors will
+            // be marked later, because the condition is symmetric and the neighbors
+            // will fullfil the condition, too.
+            pix.SetPixelUsed();
+            pix.SetPixelCore();
+            break;
+        }
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+//  Check if MEvtHeader exists in the Parameter list already.
+//  if not create one and add them to the list
+//
+Int_t MImgCleanSimple::PreProcess(MParList *pList)
+{
+    fCam = (MGeomCam*)pList->FindObject(AddSerialNumber(fNameGeomCam), "MGeomCam");
+    if (!fCam)
+    {
+        *fLog << err << fNameGeomCam << " [MGeomCam] not found (no geometry information available)... aborting." << endl;
+        return kFALSE;
+    }
+    fEvt = (MSignalCam*)pList->FindObject(AddSerialNumber(fNameSignalCam), "MSignalCam");
+    if (!fEvt)
+    {
+        *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    Print();
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Cleans the image.
+//
+Int_t MImgCleanSimple::Process()
+{
+    ResetCleaning();
+
+    switch (fCleaningMethod)
+    {
+    case kStandard:
+        DoCleanStd();   // Check each of the neighbors charge individually
+        break;
+    case kSum2:
+        DoCleanSum2();  // Check the sum of two neighbors
+    case kSum3:
+        DoCleanSum2();  // Check the sum of two   neighbors
+        DoCleanSum3();  // Check the sum of three neighbors
+        break;
+    }
+
+    // Takes roughly 10% of the time
+    fEvt->CalcIslands(*fCam);
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+//  Print descriptor and cleaning levels.
+//
+void MImgCleanSimple::Print(Option_t *o) const
+{
+    *fLog << all << GetDescriptor() << " using ";
+    switch (fCleaningMethod)
+    {
+    case kStandard:
+        *fLog << "standard";
+        break;
+    case kSum2:
+        *fLog << "sum2";
+        break;
+    case kSum3:
+        *fLog << "sum3";
+        break;
+
+    }
+    /*
+    *fLog << " cleaning" << endl;
+    *fLog << "initialized with level " << fCleanLvl1 << " and " << fCleanLvl2;
+    *fLog << " (CleanRings=" << fCleanRings << ")" << endl;
+
+    if (fPostCleanType)
+    {
+        *fLog << "Time post cleaning is switched on with:" << endl;
+        if (fPostCleanType&2)
+            *fLog << " - Two pixel coincidence window: " << fTimeLvl2 << "ns" << endl;
+        if (fPostCleanType&1)
+            *fLog << " - One pixel coincidence window: " << fTimeLvl1 << "ns" << endl;
+    }
+    */
+}
+
+/*
+// --------------------------------------------------------------------------
+//
+// 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 MImgCleanSimple::StreamPrimitive(ostream &out) const
+{
+    out << "   MImgCleanSimple " << GetUniqueName() << "(";
+    out << fCleanLvl1 << ", " << fCleanLvl2;
+
+    if (fName!=gsDefName || fTitle!=gsDefTitle)
+    {
+        out << ", \"" << fName << "\"";
+        if (fTitle!=gsDefTitle)
+            out << ", \"" << fTitle << "\"";
+    }
+    out << ");" << endl;
+
+    if (fCleaningMethod!=kStandard)
+    {
+        out << "   " << GetUniqueName() << ".SetMethod(MImgCleanSimple::k";
+        switch (fCleaningMethod)
+        {
+        case kScaled:      out << "Scaled";      break;
+        case kDemocratic:  out << "Democratic";  break;
+        case kProbability: out << "Probability"; break;
+        case kAbsolute:    out << "Absolute";    break;
+        case kTime    :    out << "Time";        break;
+        default:
+            break;
+        }
+        out << ");" << endl;
+    }
+    if (fCleanRings!=1)
+        out << "   " << GetUniqueName() << ".SetCleanRings(" << fCleanRings << ");" << endl;
+
+    if (gsNamePedPhotCam!=fNamePedPhotCam)
+        out << "   " << GetUniqueName() << ".SetNamePedPhotCam(\"" << fNamePedPhotCam << "\");" << endl;
+    if (gsNameGeomCam!=fNameGeomCam)
+        out << "   " << GetUniqueName() << ".SetNameGeomCam(\"" << fNameGeomCam << "\");" << endl;
+    if (gsNameSignalCam!=fNameSignalCam)
+        out << "   " << GetUniqueName() << ".SetNameSignalCam(\"" << fNameSignalCam << "\");" << endl;
+    if (fKeepIsolatedPixels)
+        out << "   " << GetUniqueName() << ".SetKeepIsolatedPixels();" << endl;
+    if (fRecoverIsolatedPixels)
+        out << "   " << GetUniqueName() << ".SetRecoverIsolatedPixels(" << fRecoverIsolatedPixels << ");" << endl;
+
+}
+*/
+
+// --------------------------------------------------------------------------
+//
+// Read the setup from a TEnv, eg:
+//   MImgCleanSimple.CleanLevel1: 3.0
+//   MImgCleanSimple.CleanLevel2: 2.5
+//   MImgCleanSimple.CleanMethod: Standard, Scaled, Democratic, Probability, Absolute
+//   MImgCleanSimple.CleanRings:  1
+//   MImgCleanSimple.KeepIsolatedPixels: yes, no
+//
+Int_t MImgCleanSimple::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "CleanLevel1", print))
+    {
+        rc = kTRUE;
+        fCleanLvl1 = GetEnvValue(env, prefix, "CleanLevel1", fCleanLvl1);
+    }
+    if (IsEnvDefined(env, prefix, "CleanLevel2", print))
+    {
+        rc = kTRUE;
+        fCleanLvl2 = GetEnvValue(env, prefix, "CleanLevel2", fCleanLvl2);
+    }
+    if (IsEnvDefined(env, prefix, "TimeLevel1", print))
+    {
+        rc = kTRUE;
+        fTimeLvl1 = GetEnvValue(env, prefix, "TimeLevel1", fTimeLvl1);
+    }
+    if (IsEnvDefined(env, prefix, "TimeLevel2", print))
+    {
+        rc = kTRUE;
+        fTimeLvl2 = GetEnvValue(env, prefix, "TimeLevel2", fTimeLvl2);
+    }
+
+    if (IsEnvDefined(env, prefix, "CleanMethod", print))
+    {
+        rc = kTRUE;
+        TString s = GetEnvValue(env, prefix, "CleanMethod", "");
+        s.ToLower();
+        if (s.BeginsWith("standard"))
+            SetMethod(kStandard);
+        if (s.BeginsWith("sum2"))
+            SetMethod(kSum2);
+        if (s.BeginsWith("sum3"))
+            SetMethod(kSum3);
+    }
+
+    return rc;
+}
Index: /trunk/Mars/mimage/MImgCleanSimple.h
===================================================================
--- /trunk/Mars/mimage/MImgCleanSimple.h	(revision 14202)
+++ /trunk/Mars/mimage/MImgCleanSimple.h	(revision 14202)
@@ -0,0 +1,77 @@
+#ifndef MARS_MImgCleanSimple
+#define MARS_MImgCleanSimple
+
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+class MGeomCam;
+class MGeom;
+class MSignalCam;
+class MArrivalTime;
+
+class MImgCleanSimple : public MTask
+{
+public:
+    typedef enum {
+        kStandard,
+        kSum2,
+        kSum3,
+    } CleaningMethod_t;
+
+private:
+    static const TString gsNameGeomCam;    // default name of the 'MGeomCam' container
+    static const TString gsNameSignalCam;  // default name of the 'MSignalCam' container
+
+    const MGeomCam   *fCam;  //!
+          MSignalCam *fEvt;  //!
+
+    CleaningMethod_t fCleaningMethod;
+
+    Float_t  fCleanLvl1;
+    Float_t  fCleanLvl2;
+
+    Float_t  fTimeLvl1;
+    Float_t  fTimeLvl2;
+
+    TString  fNameGeomCam;    // name of the 'MGeomCam' container
+    TString  fNameSignalCam;  // name of the 'MSignalCam' container
+
+    // MImgCleanSimple
+    void ResetCleaning() const;
+    void DoCleanStd() const;
+    void DoCleanSum2() const;
+    void DoCleanSum3() const;
+
+    // MGTask, MTask, MParContainer
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+    //void  StreamPrimitive(ostream &out) const;
+
+    Int_t PreProcess(MParList *pList);
+    Int_t Process();
+
+public:
+    MImgCleanSimple(const Float_t lvl1=2.0, const Float_t lvl2=1.75,
+                 const char *name=NULL, const char *title=NULL);
+
+    void Print(Option_t *o="") const;
+
+    Float_t GetCleanLvl1() const { return fCleanLvl1; }
+    Float_t GetCleanLvl2() const { return fCleanLvl2; }
+    Float_t GetTimeLvl1() const  { return fTimeLvl1; }
+    Float_t GetTimeLvl2() const  { return fTimeLvl2; }
+
+    void SetCleanLvl1(Float_t lvl) { fCleanLvl1=lvl; }
+    void SetCleanLvl2(Float_t lvl) { fCleanLvl2=lvl; }
+    void SetTimeLvl1(Float_t lvl)  { fTimeLvl1=lvl; }
+    void SetTimeLvl2(Float_t lvl)  { fTimeLvl2=lvl; }
+
+    void SetMethod(CleaningMethod_t m) { fCleaningMethod = m; }
+
+    void SetNameSignalCam(const char *name)   { fNameSignalCam = name; }
+    void SetNameGeomCam(const char *name)     { fNameGeomCam = name; }
+
+    ClassDef(MImgCleanSimple, 1)    // task doing the image cleaning
+}; 
+
+#endif