Index: trunk/Mars/mmuon/MHSingleMuon.cc
===================================================================
--- trunk/Mars/mmuon/MHSingleMuon.cc	(revision 10166)
+++ trunk/Mars/mmuon/MHSingleMuon.cc	(revision 14896)
@@ -211,6 +211,9 @@
     for (Int_t i=0; i<entries; i++)
     {
-        const MSignalPix &pix  = (*fSignalCam)[i];
-        const MGeom      &gpix = (*fGeomCam)[i];
+        const MSignalPix &pix = (*fSignalCam)[i];
+        if (fUseCleanedSignal && !pix.IsPixelUsed())
+            continue;
+
+        const MGeom  &gpix = (*fGeomCam)[i];
 
         const Float_t dx = gpix.GetX() - cenx;
@@ -227,5 +230,5 @@
         }
 
-        // use only the inner pixles. FIXME: This is geometry dependent
+        // use only the inner pixels. FIXME: This is geometry dependent
         if (gpix.GetAidx()>0)
             continue;
@@ -234,26 +237,34 @@
     }
 
-    // Setup the function and perform the fit
-    TF1 g1("g1", "gaus");//, -fHistTime.GetXmin(), fHistTime.GetXmax());
-
-    // Choose starting values as accurate as possible
-    g1.SetParameter(0, fHistTime.GetMaximum());
-    g1.SetParameter(1, 0);
-    g1.SetParameter(2, 0.7); // FIXME! GetRMS instead???
-
-    // According to fMuonSearchPar->GetTimeRMS() identified muons
-    // do not have an arrival time rms>3
-    g1.SetParLimits(1, -1.7, 1.7);
-    g1.SetParLimits(2,  0,   3.4);
-
-    // options : N  do not store the function, do not draw
-    //           I  use integral of function in bin rather than value at bin center
-    //           R  use the range specified in the function range
-    //           Q  quiet mode
-    if (fHistTime.Fit(&g1, "QNB"))
-        return kTRUE;
-
-    fRelTimeMean  = g1.GetParameter(1);
-    fRelTimeSigma = g1.GetParameter(2);
+    if (!fUseCleanedSignal)
+    {
+        // Setup the function and perform the fit
+        TF1 g1("g1", "gaus");//, -fHistTime.GetXmin(), fHistTime.GetXmax());
+
+        // Choose starting values as accurate as possible
+        g1.SetParameter(0, fHistTime.GetMaximum());
+        g1.SetParameter(1, 0);
+        g1.SetParameter(2, 0.7); // FIXME! GetRMS instead???
+
+        // According to fMuonSearchPar->GetTimeRMS() identified muons
+        // do not have an arrival time rms>3
+        g1.SetParLimits(1, -1.7, 1.7);
+        g1.SetParLimits(2,  0,   3.4);
+
+        // options : N  do not store the function, do not draw
+        //           I  use integral of function in bin rather than value at bin center
+        //           R  use the range specified in the function range
+        //           Q  quiet mode
+        if (fHistTime.Fit(&g1, "QNB"))
+            return kTRUE;
+
+        fRelTimeMean  = g1.GetParameter(1);
+        fRelTimeSigma = g1.GetParameter(2);
+    }
+    else
+    {
+        fRelTimeMean  = fMuonSearchPar->GetTime();
+        fRelTimeSigma = fMuonSearchPar->GetTimeRms();
+    }
 
     // The mean arrival time which was subtracted before will
@@ -264,5 +275,8 @@
     {
         const MSignalPix &pix  = (*fSignalCam)[i];
-        const MGeom      &gpix = (*fGeomCam)[i];
+        if (fUseCleanedSignal && !pix.IsPixelUsed())
+            continue;
+
+        const MGeom &gpix = (*fGeomCam)[i];
 
         const Float_t dx = gpix.GetX() - cenx;
@@ -273,5 +287,5 @@
         // if the signal is not near the estimated circle, it is ignored.
         if (TMath::Abs(dist-fMuonSearchPar->GetRadius())<fMargin &&
-            TMath::Abs(pix.GetArrivalTime()-tm0) < 2*fRelTimeSigma)
+            (fUseCleanedSignal || TMath::Abs(pix.GetArrivalTime()-tm0) < 2*fRelTimeSigma))
         {
             fHistPhi.Fill(TMath::ATan2(dx, dy)*TMath::RadToDeg(), pix.GetNumPhotons());
@@ -393,5 +407,4 @@
 {
     Int_t first, last;
-
     if (!FindRangeAboveThreshold(fHistWidth, thres, first, last))
         return kFALSE;
Index: trunk/Mars/mmuon/MHSingleMuon.h
===================================================================
--- trunk/Mars/mmuon/MHSingleMuon.h	(revision 10166)
+++ trunk/Mars/mmuon/MHSingleMuon.h	(revision 14896)
@@ -23,4 +23,5 @@
 
     Double_t fMargin;               //!
+    Bool_t   fUseCleanedSignal;     //!
 
     TProfile fHistPhi;    // Histogram of photon distribution along the arc.
@@ -49,4 +50,6 @@
     Double_t GetRelTimeSigma() const { return fRelTimeSigma; }
 
+    void SetUseCleanedSignal(Bool_t b=kTRUE) { fUseCleanedSignal = b; }
+
     Float_t CalcSize() const;
 
Index: trunk/Mars/mmuon/MMuonSearchPar.cc
===================================================================
--- trunk/Mars/mmuon/MMuonSearchPar.cc	(revision 10166)
+++ trunk/Mars/mmuon/MMuonSearchPar.cc	(revision 14896)
@@ -295,5 +295,97 @@
 
     //SetReadyToSave();
-} 
+}
+
+Bool_t MMuonSearchPar::CalcFact(const MGeomCam &geom, const MSignalCam &evt)
+{
+    // ===================== Reset cleaning ======================
+
+    for (UInt_t idx=0; idx<evt.GetNumPixels(); idx++)
+    {
+        MSignalPix &pix = evt[idx];
+        if (pix.IsPixelUnmapped())
+            continue;
+
+        pix.SetPixelUnused();
+        pix.SetPixelCore(kFALSE);
+    }
+
+    // ============ Do muon cleaning / calculate COG =============
+
+    // The window must be large enough that the earliest and latest
+    // events do not get a biased timerms
+    double time_min   =  -19.5;
+    double time_max   =  -4.5;
+    double signal_min =  2.30;
+    double delta_t    =  1.75;
+
+    Double_t sumx = 0.;
+    Double_t sumy = 0.;
+    Double_t sumw = 0.;
+    for (UInt_t i=0; i<evt.GetNumPixels(); i++)
+    {
+        MSignalPix &pix = evt[i];
+
+        if (pix.GetNumPhotons()<signal_min)
+            continue;
+        if (pix.GetArrivalTime()>=time_max || pix.GetArrivalTime()<time_min)
+            continue;
+        if (pix.IsPixelUnmapped())
+            continue;
+
+        const MGeom &gpix = geom[i];
+
+        const double x = gpix.GetX();
+        const double y = gpix.GetY();
+
+        int counter = 0;
+        for (int j=0; j<gpix.GetNumNeighbors(); j++)
+        {
+            const int idx = gpix.GetNeighbor(j);
+
+            const MSignalPix &spix = evt[idx];
+
+            if (spix.GetNumPhotons()<signal_min)
+                continue;
+            if (spix.GetArrivalTime()>=time_max || spix.GetArrivalTime()<time_min)
+                continue;
+            if (spix.IsPixelUnmapped())
+                continue;
+
+            if (TMath::Abs(pix.GetArrivalTime()-spix.GetArrivalTime())>delta_t)
+                continue;
+
+            counter++;
+        }
+
+        if (counter==0)
+            continue;
+
+        sumx += pix.GetNumPhotons()*x;
+        sumy += pix.GetNumPhotons()*y;
+        sumw += pix.GetNumPhotons();
+
+        pix.SetPixelUsed();
+        pix.SetPixelCore();
+    }
+
+    if (sumw==0)
+        return kFALSE;
+
+    sumx /= sumw;
+    sumy /= sumw;
+
+    // ==============  Fit circle to resulting data ===============
+
+    Double_t sigma, rad;
+    CalcMinimumDeviation(geom, evt, sumx, sumy, sigma, rad);
+
+    fCenterX   = sumx;
+    fCenterY   = sumy;
+    fRadius    = rad;
+    fDeviation = sigma;
+
+    return kTRUE;
+}
 
 void MMuonSearchPar::Print(Option_t *) const
Index: trunk/Mars/mmuon/MMuonSearchPar.h
===================================================================
--- trunk/Mars/mmuon/MMuonSearchPar.h	(revision 10166)
+++ trunk/Mars/mmuon/MMuonSearchPar.h	(revision 14896)
@@ -52,4 +52,5 @@
     void   Calc(const MGeomCam &geom, const MSignalCam &evt,
 		const MHillas &hillas);
+    Bool_t CalcFact(const MGeomCam &geom, const MSignalCam &evt);
 
     // TObject
Index: trunk/Mars/mmuon/MMuonSearchParCalc.cc
===================================================================
--- trunk/Mars/mmuon/MMuonSearchParCalc.cc	(revision 10166)
+++ trunk/Mars/mmuon/MMuonSearchParCalc.cc	(revision 14896)
@@ -76,9 +76,13 @@
 Int_t MMuonSearchParCalc::PreProcess(MParList *pList)
 {
-    fHillas = (MHillas*)pList->FindObject("MHillas");
-    if (!fHillas)
+    fHillas = 0;
+    if (!fUseFactAlgorithm)
     {
-        *fLog << err << "MHillas not found... aborting." << endl;
-        return kFALSE;
+        fHillas = (MHillas*)pList->FindObject("MHillas");
+        if (!fHillas)
+        {
+            *fLog << err << "MHillas not found... aborting." << endl;
+            return kFALSE;
+        }
     }
 
@@ -108,5 +112,9 @@
 Int_t MMuonSearchParCalc::Process()
 {
-    fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas);
+    if (fUseFactAlgorithm)
+        fMuonPar->CalcFact(*fGeomCam, *fSignalCam);
+    else
+        fMuonPar->Calc(*fGeomCam, *fSignalCam, *fHillas);
+
     return kTRUE;
 }
Index: trunk/Mars/mmuon/MMuonSearchParCalc.h
===================================================================
--- trunk/Mars/mmuon/MMuonSearchParCalc.h	(revision 10166)
+++ trunk/Mars/mmuon/MMuonSearchParCalc.h	(revision 14896)
@@ -19,4 +19,6 @@
     MMuonSearchPar *fMuonPar;     //! Pointer to the output container for the new image parameters
 
+    Bool_t fUseFactAlgorithm;
+
     Int_t PreProcess(MParList *plist);
     Int_t Process();
@@ -25,4 +27,6 @@
     MMuonSearchParCalc(const char *name=NULL, const char *title=NULL);
 
+    void SetUseFactAlgorithm(Bool_t b=kTRUE) { fUseFactAlgorithm = kTRUE; }
+
     ClassDef(MMuonSearchParCalc, 0) // task to calculate muon parameters
 };