Index: trunk/Mars/mcore/DrsCalib.h
===================================================================
--- trunk/Mars/mcore/DrsCalib.h	(revision 14864)
+++ trunk/Mars/mcore/DrsCalib.h	(revision 14865)
@@ -278,8 +278,8 @@
     };
 
-    static Step AverageSteps(const std::vector<Step> &list, uint32_t n)
+    static Step AverageSteps(const std::vector<Step>::iterator beg, const std::vector<Step>::iterator end)
     {
         Step rc;
-        for (auto it=list.begin(); it!=list.begin()+n; it++)
+        for (auto it=beg; it!=end; it++)
         {
             rc.pos += it->pos;
@@ -288,5 +288,5 @@
         }
 
-        rc.cnt = n;
+        rc.cnt = end-beg;
 
         rc.pos /= rc.cnt;
@@ -329,15 +329,21 @@
             return Step();
 
-        Step rc = AverageSteps(list, list.size());;
+        Step rc = AverageSteps(list.begin(), list.begin()+list.size());;
 
         if (rc.avg==0)
             return Step();
 
+        // std::cout << "   A0=" << rc.avg << "    rms=" << rc.rms << std::endl;
         if (rc.rms>5)
         {
-            partial_sort(list.begin(), list.begin()+list.size()/2,
-                         list.end(), Step::sort);
-
-            rc = AverageSteps(list, list.size()/2);
+            sort(list.begin(), list.end(), Step::sort);
+
+            //for (auto it=list.begin(); it!=list.end(); it++)
+            //    std::cout << "     " << it->avg << std::endl;
+
+            const Int_t skip = list.size()/10;
+            rc = AverageSteps(list.begin()+skip, list.begin()+list.size()-skip);
+
+            // std::cout << "   A1=" << rc.avg << "    rms=" << rc.rms << std::endl;
         }
 
@@ -387,26 +393,23 @@
             std::vector<float> Ameas(p, p+roi);
 
-            std::vector<float> N1mean(roi);
-            for (size_t i=2; i<roi-2; i++)
-            {
-                N1mean[i] = (p[i-1] + p[i+1])/2.;
-            }
+            std::vector<float> diff(roi);
+            for (size_t i=1; i<roi-1; i++)
+                diff[i] = (p[i-1] + p[i+1])/2 - p[i];
+
+            //std::vector<float> N1mean(roi);
+            //for (size_t i=1; i<roi-1; i++)
+            //    N1mean[i] = (p[i-1] + p[i+1])/2;
 
             const float fract = 0.8;
-            float x, xp, xpp;
 
             for (size_t i=0; i<roi-3; i++)
             {
-                x = Ameas[i] - N1mean[i];
-                if ( x > -5.)
+                if (diff[i]<5)
                     continue;
 
-                xp  = Ameas[i+1] - N1mean[i+1];
-                xpp = Ameas[i+2] - N1mean[i+2];
-
-                if(Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2. > 10.)
+                if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2 > 10)
                 {
-                    p[i+1]=(Ameas[i] + Ameas[i+3])/2.;
-                    p[i+2]=(Ameas[i] + Ameas[i+3])/2.;
+                    p[i+1]=   (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+                    p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
 
                     i += 3;
@@ -414,13 +417,100 @@
                     continue;
                 }
-                if ( (xp > -2.*x*fract) && (xpp < -10.) )
+
+                if ( (diff[i+1]<-diff[i]*fract*2) && (diff[i+2]>10) )
                 {
-                    p[i+1] = N1mean[i+1];
-                    N1mean[i+2] = (Ameas[i+1] - Ameas[i+3] / 2.);
+                    p[i+1]    = (Ameas[i]+Ameas[i+2])/2;
+                    diff[i+2] = (p[i+1] + Ameas[i+3])/2 - Ameas[i+2];
 
                     i += 2;
                 }
-            }
-        }
+
+                // const float x = Ameas[i] - N1mean[i];
+                // if (x > -5.)
+                //     continue;
+
+                // if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2. > 10.)
+                // {
+                //     p[i+1]=   (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+                //     p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
+                //     i += 3;
+                //     continue;
+                // }
+
+                // const float xp  = Ameas[i+1] - N1mean[i+1];
+                // const float xpp = Ameas[i+2] - N1mean[i+2];
+
+                // if ( (xp > -2.*x*fract) && (xpp < -10.) )
+                // {
+                //     p[i+1] = N1mean[i+1];
+                //     N1mean[i+2] = Ameas[i+1] - Ameas[i+3]/2;
+                //
+                //     i += 2;
+                // }
+            }
+        }
+    }
+
+    static void RemoveSpikes3(float *vec, uint32_t roi)//from Werner
+    {
+        const float SingleCandidateTHR = -10.;
+        const float DoubleCandidateTHR =  -5.;
+
+        const std::vector<float> src(vec, vec+roi);
+ 
+        std::vector<float> diff(roi);
+        for (size_t i=1; i<roi-1; i++)
+            diff[i] = src[i] - (src[i-1] + src[i+1])/2;
+
+        // find the spike and replace it by mean value of neighbors
+        for (unsigned int i=1; i<roi-3; i++)
+        {
+            // Speed up (no leading edge)
+            if (diff[i]>=DoubleCandidateTHR)
+                continue;
+
+            //bool checkDouble = false;
+
+            // a single spike candidate
+            if (diff[i]<SingleCandidateTHR)
+            {
+                // check consistency with a single channel spike
+                if (diff[i+1] > -1.6*diff[i])
+                {
+                    vec[i+1] = (src[i] + src[i+2]) / 2;
+
+                    i += 2;
+
+                    /*** NEW ***/
+                    continue;
+                    /*** NEW ***/
+                }
+                /*
+                else
+                {
+                    // do nothing - not really a single spike,
+                    // but check if it is a double
+                    checkDouble = true;
+                }*/
+            }
+
+            // a double spike candidate
+            //if (diff[i]>DoubleCandidateTHR || checkDouble == 1)
+            { 
+                // check the consistency with a double spike
+                if ((diff[i+1] > -DoubleCandidateTHR) &&
+                    (diff[i+2] > -DoubleCandidateTHR))
+                {
+                    vec[i+1] =   (src[i+3] - src[i])/3 + src[i];
+                    vec[i+2] = 2*(src[i+3] - src[i])/3 + src[i];
+
+                    //vec[i]   = (src[i-1] + src[i+2]) / 2.;
+                    //vec[i+1] = (src[i-1] + src[i+2]) / 2.;
+
+                    //do not care about the next sample it was the spike
+                    i += 3;
+                }
+            }
+        } 
     }