diff --git a/doc/data.rst b/doc/data.rst
index 053ed1e4..35ebd947 100644
--- a/doc/data.rst
+++ b/doc/data.rst
@@ -637,11 +637,24 @@ and functions:
     Calculates Shirley background
     (useful in X-ray photoelectron spectroscopy).
 
+``snip_bg(@n, window_width, direction, filter_order, estimate_compton)``
+    Calculates the spectrum background applying a Sensitive Nonlinear Iterative Peak clipping algorithm (SNIP) (useful in gamma-ray spectroscopy).
+    This SNIP algorithm is taken from the `TSpectrum <https://root.cern.ch/doc/v608/classTSpectrum.html>`_ class of `ROOT <http://root.cern.ch>`_.
+    It requires the parameters:
+    1. ``window_width``: maximal width of clipping window;
+    2. ``direction``:  direction of change of clipping window, if it is a positive number the window will be increasing;
+    3. ``filter_order``:  order of clipping filter, possible values are 2, 4, 6 or 8;
+    4. ``estimate_compton``:  if a positive number the algorithm will try to estimate the Compton edge.
+    Active points are modified according to the SNIP algorithm.
+    Inactive points are left as they are, assuming that they are already considered background.
+    Subtracting the calculated background from the initial spectrum gives zero in all inactive points.
+
 Examples::
 
   @+ = @0 # duplicate the dataset
   @+ = @0 and @1 # create a new dataset from @0 and @1
   @0 = @0 - shirley_bg(@0) # remove Shirley background 
+  @0 = @0 - snip_bg(@0, 30, 1, 2, 0) # remove gamma background 
   @0 = @0 - @1 # subtract @1 from @0
   @0 = @0 - 0.28*@1 # subtract scaled dataset @1 from @0
 
diff --git a/fityk/Makefile.am b/fityk/Makefile.am
index 1eb1ab43..b9ba04bf 100644
--- a/fityk/Makefile.am
+++ b/fityk/Makefile.am
@@ -13,6 +13,7 @@ libfityk_la_SOURCES = logic.cpp view.cpp lexer.cpp eparser.cpp cparser.cpp \
 		 runner.cpp info.cpp common.cpp data.cpp var.cpp mgr.cpp \
 		 tplate.cpp func.cpp udf.cpp bfunc.cpp f_fcjasym.cpp ast.cpp \
 		 vm.cpp transform.cpp settings.cpp ui.cpp ui_api.cpp \
+		 root/background.cpp \
 		 luabridge.cpp GAfit.cpp LMfit.cpp guess.cpp NMfit.cpp \
 		 model.cpp fit.cpp voigt.cpp numfuncs.cpp fityk.cpp \
 		 \
@@ -20,6 +21,7 @@ libfityk_la_SOURCES = logic.cpp view.cpp lexer.cpp eparser.cpp cparser.cpp \
 		 runner.h info.h common.h data.h var.h mgr.h \
 		 tplate.h func.h udf.h bfunc.h f_fcjasym.h ast.h \
 		 vm.h transform.h settings.h ui.h luabridge.h \
+		 root/background.hpp \
 		 GAfit.h LMfit.h guess.h NMfit.h \
 		 model.h fit.h voigt.h numfuncs.h \
 		 swig/fityk_lua.cpp swig/luarun.h \
diff --git a/fityk/eparser.cpp b/fityk/eparser.cpp
index b21a7db2..99d9786e 100644
--- a/fityk/eparser.cpp
+++ b/fityk/eparser.cpp
@@ -79,6 +79,7 @@ const char* function_name(int op)
         case OP_DT_SUM_SAME_X: return "sum_same_x";
         case OP_DT_AVG_SAME_X: return "avg_same_x";
         case OP_DT_SHIRLEY_BG: return "shirley_bg";
+        case OP_DT_SNIP_BG: return "snip_bg";
         // 2-args functions
         case OP_MOD: return "mod";
         case OP_MIN2: return "min2";
@@ -138,6 +139,8 @@ int get_function_narg(int op)
         case OP_RANDNORM:
         case OP_RANDU:
             return 2;
+        case OP_DT_SNIP_BG:
+            return 5;
         // Fityk functions
         case OP_FUNC:
         case OP_SUM_F:
@@ -819,6 +822,8 @@ void ExpressionParser::parse_expr(Lexer& lex, int default_ds,
                         put_function(OP_DT_AVG_SAME_X);
                     else if (mode == kDatasetTrMode && word == "shirley_bg")
                         put_function(OP_DT_SHIRLEY_BG);
+                    else if (mode == kDatasetTrMode && word == "snip_bg")
+                        put_function(OP_DT_SNIP_BG);
 
                     else
                         lex.throw_syntax_error("unknown function: " + word);
diff --git a/fityk/root/background.cpp b/fityk/root/background.cpp
index c39459f2..f18bb323 100644
--- a/fityk/root/background.cpp
+++ b/fityk/root/background.cpp
@@ -15,6 +15,7 @@
 #include <cmath>
 #include <vector>
 
+#include "fityk.h"
 #include "background.hpp"
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -441,7 +442,7 @@
 /// }
 /// ~~~
 
-std::vector<double> background::background(std::vector<double> spectrum,
+std::vector<fityk::Point> ROOT::background(const std::vector<fityk::Point> spectrum,
                                            int numberIterations,
                                            int direction,
                                            int filterOrder,
@@ -455,18 +456,21 @@ std::vector<double> background::background(std::vector<double> spectrum,
     const unsigned int ssize = spectrum.size();
 
     if (ssize <= 0)
-        throw "Wrong Parameters";
+        throw fityk::ExecuteError("Wrong vector size");
 
     if (numberIterations < 1)
-        throw "Width of Clipping Window Must Be Positive";
+        throw fityk::ExecuteError("Width of Clipping Window Must Be Positive");
 
     if (ssize < 2 * numberIterations + 1)
-        throw "Too Large Clipping Window";
+        throw fityk::ExecuteError("Too Large Clipping Window");
+
+    if (filterOrder != kBackOrder2 && filterOrder != kBackOrder4 && filterOrder != kBackOrder6 && filterOrder != kBackOrder8)
+        throw fityk::ExecuteError("Incorrect order of clipping filter, possible values are: 2, 4, 6 or 8");
 
     if (smoothing == true && smoothWindow != kBackSmoothing3 && smoothWindow != kBackSmoothing5 && smoothWindow != kBackSmoothing7 && smoothWindow != kBackSmoothing9 && smoothWindow != kBackSmoothing11 && smoothWindow != kBackSmoothing13 && smoothWindow != kBackSmoothing15)
-        throw "Incorrect width of smoothing window";
+        throw fityk::ExecuteError("Incorrect width of smoothing window");
 
-   std::vector<double> working_space(2 * ssize);
+   std::vector<fityk::Point> working_space(2 * ssize);
 
    for (i = 0; i < ssize; i++)
     {
@@ -485,20 +489,20 @@ std::vector<double> background::background(std::vector<double> spectrum,
       do{
          for (j = i; j < ssize - i; j++) {
             if (smoothing == false){
-               a = working_space[ssize + j];
-               b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
+               a = working_space[ssize + j].y;
+               b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0;
                if (b < a)
                   a = b;
-               working_space[j] = a;
+               working_space[j].y = a;
             }
 
             else if (smoothing == true){
-               a = working_space[ssize + j];
+               a = working_space[ssize + j].y;
                av = 0;
                men = 0;
                for (w = j - bw; w <= j + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     av += working_space[ssize + w];
+                     av += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -507,7 +511,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j - i - bw; w <= j - i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     b += working_space[ssize + w];
+                     b += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -516,7 +520,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j + i - bw; w <= j + i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     c += working_space[ssize + w];
+                     c += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -524,11 +528,11 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b = (b + c) / 2;
                if (b < a)
                   av = b;
-               working_space[j]=av;
+               working_space[j].y=av;
             }
          }
          for (j = i; j < ssize - i; j++)
-            working_space[ssize + j] = working_space[j];
+            working_space[ssize + j].y = working_space[j].y;
          if (direction == kBackIncreasingWindow)
             i+=1;
          else if(direction == kBackDecreasingWindow)
@@ -540,28 +544,28 @@ std::vector<double> background::background(std::vector<double> spectrum,
       do{
          for (j = i; j < ssize - i; j++) {
             if (smoothing == false){
-               a = working_space[ssize + j];
-               b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
+               a = working_space[ssize + j].y;
+               b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0;
                c = 0;
                ai = i / 2;
-               c -= working_space[ssize + j - (int) (2 * ai)] / 6;
-               c += 4 * working_space[ssize + j - (int) ai] / 6;
-               c += 4 * working_space[ssize + j + (int) ai] / 6;
-               c -= working_space[ssize + j + (int) (2 * ai)] / 6;
+               c -= working_space[ssize + j - (int) (2 * ai)].y / 6;
+               c += 4 * working_space[ssize + j - (int) ai].y / 6;
+               c += 4 * working_space[ssize + j + (int) ai].y / 6;
+               c -= working_space[ssize + j + (int) (2 * ai)].y / 6;
                if (b < c)
                   b = c;
                if (b < a)
                   a = b;
-               working_space[j] = a;
+               working_space[j].y = a;
             }
 
             else if (smoothing == true){
-               a = working_space[ssize + j];
+               a = working_space[ssize + j].y;
                av = 0;
                men = 0;
                for (w = j - bw; w <= j + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     av += working_space[ssize + w];
+                     av += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -570,7 +574,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j - i - bw; w <= j - i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     b += working_space[ssize + w];
+                     b += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -579,7 +583,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j + i - bw; w <= j + i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     c += working_space[ssize + w];
+                     c += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -589,7 +593,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b4 = 0, men = 0;
                for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     b4 += working_space[ssize + w];
+                     b4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -597,7 +601,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                c4 = 0, men = 0;
                for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     c4 += working_space[ssize + w];
+                     c4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -605,7 +609,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                d4 = 0, men = 0;
                for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     d4 += working_space[ssize + w];
+                     d4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -613,7 +617,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                e4 = 0, men = 0;
                for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     e4 += working_space[ssize + w];
+                     e4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -623,11 +627,11 @@ std::vector<double> background::background(std::vector<double> spectrum,
                   b = b4;
                if (b < a)
                   av = b;
-               working_space[j]=av;
+               working_space[j].y=av;
             }
          }
          for (j = i; j < ssize - i; j++)
-            working_space[ssize + j] = working_space[j];
+            working_space[ssize + j].y = working_space[j].y;
          if (direction == kBackIncreasingWindow)
             i+=1;
          else if(direction == kBackDecreasingWindow)
@@ -639,38 +643,38 @@ std::vector<double> background::background(std::vector<double> spectrum,
       do{
          for (j = i; j < ssize - i; j++) {
             if (smoothing == false){
-               a = working_space[ssize + j];
-               b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
+               a = working_space[ssize + j].y;
+               b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0;
                c = 0;
                ai = i / 2;
-               c -= working_space[ssize + j - (int) (2 * ai)] / 6;
-               c += 4 * working_space[ssize + j - (int) ai] / 6;
-               c += 4 * working_space[ssize + j + (int) ai] / 6;
-               c -= working_space[ssize + j + (int) (2 * ai)] / 6;
+               c -= working_space[ssize + j - (int) (2 * ai)].y / 6;
+               c += 4 * working_space[ssize + j - (int) ai].y / 6;
+               c += 4 * working_space[ssize + j + (int) ai].y / 6;
+               c -= working_space[ssize + j + (int) (2 * ai)].y / 6;
                d = 0;
                ai = i / 3;
-               d += working_space[ssize + j - (int) (3 * ai)] / 20;
-               d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20;
-               d += 15 * working_space[ssize + j - (int) ai] / 20;
-               d += 15 * working_space[ssize + j + (int) ai] / 20;
-               d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20;
-               d += working_space[ssize + j + (int) (3 * ai)] / 20;
+               d += working_space[ssize + j - (int) (3 * ai)].y / 20;
+               d -= 6 * working_space[ssize + j - (int) (2 * ai)].y / 20;
+               d += 15 * working_space[ssize + j - (int) ai].y / 20;
+               d += 15 * working_space[ssize + j + (int) ai].y / 20;
+               d -= 6 * working_space[ssize + j + (int) (2 * ai)].y / 20;
+               d += working_space[ssize + j + (int) (3 * ai)].y / 20;
                if (b < d)
                   b = d;
                if (b < c)
                   b = c;
                if (b < a)
                   a = b;
-               working_space[j] = a;
+               working_space[j].y = a;
             }
 
             else if (smoothing == true){
-               a = working_space[ssize + j];
+               a = working_space[ssize + j].y;
                av = 0;
                men = 0;
                for (w = j - bw; w <= j + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     av += working_space[ssize + w];
+                     av += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -679,7 +683,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j - i - bw; w <= j - i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     b += working_space[ssize + w];
+                     b += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -688,7 +692,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j + i - bw; w <= j + i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     c += working_space[ssize + w];
+                     c += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -698,7 +702,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b4 = 0, men = 0;
                for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     b4 += working_space[ssize + w];
+                     b4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -706,7 +710,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                c4 = 0, men = 0;
                for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     c4 += working_space[ssize + w];
+                     c4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -714,7 +718,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                d4 = 0, men = 0;
                for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     d4 += working_space[ssize + w];
+                     d4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -722,7 +726,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                e4 = 0, men = 0;
                for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     e4 += working_space[ssize + w];
+                     e4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -732,7 +736,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b6 = 0, men = 0;
                for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     b6 += working_space[ssize + w];
+                     b6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -740,7 +744,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                c6 = 0, men = 0;
                for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     c6 += working_space[ssize + w];
+                     c6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -748,7 +752,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                d6 = 0, men = 0;
                for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     d6 += working_space[ssize + w];
+                     d6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -756,7 +760,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                e6 = 0, men = 0;
                for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     e6 += working_space[ssize + w];
+                     e6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -764,7 +768,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                f6 = 0, men = 0;
                for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     f6 += working_space[ssize + w];
+                     f6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -772,7 +776,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                g6 = 0, men = 0;
                for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     g6 += working_space[ssize + w];
+                     g6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -784,11 +788,11 @@ std::vector<double> background::background(std::vector<double> spectrum,
                   b = b4;
                if (b < a)
                   av = b;
-               working_space[j]=av;
+               working_space[j].y=av;
             }
          }
          for (j = i; j < ssize - i; j++)
-            working_space[ssize + j] = working_space[j];
+            working_space[ssize + j].y = working_space[j].y;
          if (direction == kBackIncreasingWindow)
             i+=1;
          else if(direction == kBackDecreasingWindow)
@@ -800,32 +804,32 @@ std::vector<double> background::background(std::vector<double> spectrum,
       do{
          for (j = i; j < ssize - i; j++) {
             if (smoothing == false){
-               a = working_space[ssize + j];
-               b = (working_space[ssize + j - i] + working_space[ssize + j + i]) / 2.0;
+               a = working_space[ssize + j].y;
+               b = (working_space[ssize + j - i].y + working_space[ssize + j + i].y) / 2.0;
                c = 0;
                ai = i / 2;
-               c -= working_space[ssize + j - (int) (2 * ai)] / 6;
-               c += 4 * working_space[ssize + j - (int) ai] / 6;
-               c += 4 * working_space[ssize + j + (int) ai] / 6;
-               c -= working_space[ssize + j + (int) (2 * ai)] / 6;
+               c -= working_space[ssize + j - (int) (2 * ai)].y / 6;
+               c += 4 * working_space[ssize + j - (int) ai].y / 6;
+               c += 4 * working_space[ssize + j + (int) ai].y / 6;
+               c -= working_space[ssize + j + (int) (2 * ai)].y / 6;
                d = 0;
                ai = i / 3;
-               d += working_space[ssize + j - (int) (3 * ai)] / 20;
-               d -= 6 * working_space[ssize + j - (int) (2 * ai)] / 20;
-               d += 15 * working_space[ssize + j - (int) ai] / 20;
-               d += 15 * working_space[ssize + j + (int) ai] / 20;
-               d -= 6 * working_space[ssize + j + (int) (2 * ai)] / 20;
-               d += working_space[ssize + j + (int) (3 * ai)] / 20;
+               d += working_space[ssize + j - (int) (3 * ai)].y / 20;
+               d -= 6 * working_space[ssize + j - (int) (2 * ai)].y / 20;
+               d += 15 * working_space[ssize + j - (int) ai].y / 20;
+               d += 15 * working_space[ssize + j + (int) ai].y / 20;
+               d -= 6 * working_space[ssize + j + (int) (2 * ai)].y / 20;
+               d += working_space[ssize + j + (int) (3 * ai)].y / 20;
                e = 0;
                ai = i / 4;
-               e -= working_space[ssize + j - (int) (4 * ai)] / 70;
-               e += 8 * working_space[ssize + j - (int) (3 * ai)] / 70;
-               e -= 28 * working_space[ssize + j - (int) (2 * ai)] / 70;
-               e += 56 * working_space[ssize + j - (int) ai] / 70;
-               e += 56 * working_space[ssize + j + (int) ai] / 70;
-               e -= 28 * working_space[ssize + j + (int) (2 * ai)] / 70;
-               e += 8 * working_space[ssize + j + (int) (3 * ai)] / 70;
-               e -= working_space[ssize + j + (int) (4 * ai)] / 70;
+               e -= working_space[ssize + j - (int) (4 * ai)].y / 70;
+               e += 8 * working_space[ssize + j - (int) (3 * ai)].y / 70;
+               e -= 28 * working_space[ssize + j - (int) (2 * ai)].y / 70;
+               e += 56 * working_space[ssize + j - (int) ai].y / 70;
+               e += 56 * working_space[ssize + j + (int) ai].y / 70;
+               e -= 28 * working_space[ssize + j + (int) (2 * ai)].y / 70;
+               e += 8 * working_space[ssize + j + (int) (3 * ai)].y / 70;
+               e -= working_space[ssize + j + (int) (4 * ai)].y / 70;
                if (b < e)
                   b = e;
                if (b < d)
@@ -834,16 +838,16 @@ std::vector<double> background::background(std::vector<double> spectrum,
                   b = c;
                if (b < a)
                   a = b;
-               working_space[j] = a;
+               working_space[j].y = a;
             }
 
             else if (smoothing == true){
-               a = working_space[ssize + j];
+               a = working_space[ssize + j].y;
                av = 0;
                men = 0;
                for (w = j - bw; w <= j + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     av += working_space[ssize + w];
+                     av += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -852,7 +856,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j - i - bw; w <= j - i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     b += working_space[ssize + w];
+                     b += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -861,7 +865,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                men = 0;
                for (w = j + i - bw; w <= j + i + bw; w++){
                   if ( w >= 0 && w < ssize){
-                     c += working_space[ssize + w];
+                     c += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -871,7 +875,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b4 = 0, men = 0;
                for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     b4 += working_space[ssize + w];
+                     b4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -879,7 +883,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                c4 = 0, men = 0;
                for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     c4 += working_space[ssize + w];
+                     c4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -887,7 +891,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                d4 = 0, men = 0;
                for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     d4 += working_space[ssize + w];
+                     d4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -895,7 +899,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                e4 = 0, men = 0;
                for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     e4 += working_space[ssize + w];
+                     e4 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -905,7 +909,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b6 = 0, men = 0;
                for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     b6 += working_space[ssize + w];
+                     b6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -913,7 +917,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                c6 = 0, men = 0;
                for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     c6 += working_space[ssize + w];
+                     c6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -921,7 +925,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                d6 = 0, men = 0;
                for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     d6 += working_space[ssize + w];
+                     d6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -929,7 +933,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                e6 = 0, men = 0;
                for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     e6 += working_space[ssize + w];
+                     e6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -937,7 +941,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                f6 = 0, men = 0;
                for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     f6 += working_space[ssize + w];
+                     f6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -945,7 +949,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                g6 = 0, men = 0;
                for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     g6 += working_space[ssize + w];
+                     g6 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -955,7 +959,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                b8 = 0, men = 0;
                for (w = j - (int)(4 * ai) - bw; w <= j - (int)(4 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     b8 += working_space[ssize + w];
+                     b8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -963,7 +967,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                c8 = 0, men = 0;
                for (w = j - (int)(3 * ai) - bw; w <= j - (int)(3 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     c8 += working_space[ssize + w];
+                     c8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -971,7 +975,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                d8 = 0, men = 0;
                for (w = j - (int)(2 * ai) - bw; w <= j - (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     d8 += working_space[ssize + w];
+                     d8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -979,7 +983,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                e8 = 0, men = 0;
                for (w = j - (int)ai - bw; w <= j - (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     e8 += working_space[ssize + w];
+                     e8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -987,7 +991,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                f8 = 0, men = 0;
                for (w = j + (int)ai - bw; w <= j + (int)ai + bw; w++){
                   if (w >= 0 && w < ssize){
-                     f8 += working_space[ssize + w];
+                     f8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -995,7 +999,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                g8 = 0, men = 0;
                for (w = j + (int)(2 * ai) - bw; w <= j + (int)(2 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     g8 += working_space[ssize + w];
+                     g8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -1003,7 +1007,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                h8 = 0, men = 0;
                for (w = j + (int)(3 * ai) - bw; w <= j + (int)(3 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     h8 += working_space[ssize + w];
+                     h8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -1011,7 +1015,7 @@ std::vector<double> background::background(std::vector<double> spectrum,
                i8 = 0, men = 0;
                for (w = j + (int)(4 * ai) - bw; w <= j + (int)(4 * ai) + bw; w++){
                   if (w >= 0 && w < ssize){
-                     i8 += working_space[ssize + w];
+                     i8 += working_space[ssize + w].y;
                      men +=1;
                   }
                }
@@ -1025,11 +1029,11 @@ std::vector<double> background::background(std::vector<double> spectrum,
                   b = b4;
                if (b < a)
                   av = b;
-               working_space[j]=av;
+               working_space[j].y = av;
             }
          }
          for (j = i; j < ssize - i; j++)
-            working_space[ssize + j] = working_space[j];
+            working_space[ssize + j].y = working_space[j].y;
          if (direction == kBackIncreasingWindow)
             i += 1;
          else if(direction == kBackDecreasingWindow)
@@ -1040,15 +1044,15 @@ std::vector<double> background::background(std::vector<double> spectrum,
    if (compton == true) {
       for (i = 0, b2 = 0; i < ssize; i++){
          b1 = b2;
-         a = working_space[i], b = spectrum[i];
+         a = working_space[i].y, b = spectrum[i].y;
          j = i;
          if (fabs(a - b) >= 1) {
             b1 = i - 1;
             if (b1 < 0)
                b1 = 0;
-            yb1 = working_space[b1];
+            yb1 = working_space[b1].y;
             for (b2 = b1 + 1, c = 0, priz = 0; priz == 0 && b2 < ssize; b2++){
-               a = working_space[b2], b = spectrum[b2];
+               a = working_space[b2].y, b = spectrum[b2].y;
                c = c + b - yb1;
                if (fabs(a - b) < 1) {
                   priz = 1;
@@ -1057,35 +1061,35 @@ std::vector<double> background::background(std::vector<double> spectrum,
             }
             if (b2 == ssize)
                b2 -= 1;
-            yb2 = working_space[b2];
+            yb2 = working_space[b2].y;
             if (yb1 <= yb2){
                for (j = b1, c = 0; j <= b2; j++){
-                  b = spectrum[j];
+                  b = spectrum[j].y;
                   c = c + b - yb1;
                }
                if (c > 1){
                   c = (yb2 - yb1) / c;
                   for (j = b1, d = 0; j <= b2 && j < ssize; j++){
-                     b = spectrum[j];
+                     b = spectrum[j].y;
                      d = d + b - yb1;
                      a = c * d + yb1;
-                     working_space[ssize + j] = a;
+                     working_space[ssize + j].y = a;
                   }
                }
             }
 
             else{
                for (j = b2, c = 0; j >= b1; j--){
-                  b = spectrum[j];
+                  b = spectrum[j].y;
                   c = c + b - yb2;
                }
                if (c > 1){
                   c = (yb1 - yb2) / c;
                   for (j = b2, d = 0;j >= b1 && j >= 0; j--){
-                     b = spectrum[j];
+                     b = spectrum[j].y;
                      d = d + b - yb2;
                      a = c * d + yb2;
-                     working_space[ssize + j] = a;
+                     working_space[ssize + j].y = a;
                   }
                }
             }
diff --git a/fityk/root/background.hpp b/fityk/root/background.hpp
index efe016b8..44862f1a 100644
--- a/fityk/root/background.hpp
+++ b/fityk/root/background.hpp
@@ -15,31 +15,35 @@
 #ifndef __BACKGROUND_HPP__
 #define __BACKGROUND_HPP__
 
-namespace background {
+#include <vector>
+
+#include "fityk.h"
+
+namespace ROOT {
 
     enum {
-        kBackIncreasingWindow =0,
-        kBackDecreasingWindow =1,
-        kBackOrder2 =2,
-        kBackOrder4 =4,
-        kBackOrder6 =6,
-        kBackOrder8 =8,
-        kBackSmoothing3 =3,
-        kBackSmoothing5 =5,
-        kBackSmoothing7 =7,
-        kBackSmoothing9 =9,
-        kBackSmoothing11 =11,
-        kBackSmoothing13 =13,
-        kBackSmoothing15 =15
+        kBackIncreasingWindow = 1,
+        kBackDecreasingWindow = -1,
+        kBackOrder2 = 2,
+        kBackOrder4 = 4,
+        kBackOrder6 = 6,
+        kBackOrder8 = 8,
+        kBackSmoothing3 = 3,
+        kBackSmoothing5 = 5,
+        kBackSmoothing7 = 7,
+        kBackSmoothing9 = 9,
+        kBackSmoothing11 = 11,
+        kBackSmoothing13 = 13,
+        kBackSmoothing15 = 15
    };
 
-    std::vector<double> background(std::vector<double> spectrum,
-                                   int numberIterations,
-                                   int direction,
-                                   int filterOrder,
-                                   bool smoothing,
-                                   int smoothWindow,
-                                   bool compton);
+    std::vector<fityk::Point> background(const std::vector<fityk::Point> spectrum,
+                                         int numberIterations,
+                                         int direction,
+                                         int filterOrder,
+                                         bool smoothing,
+                                         int smoothWindow,
+                                         bool compton);
 }
 
 #endif
diff --git a/fityk/transform.cpp b/fityk/transform.cpp
index 0e675713..ce4f0ce7 100644
--- a/fityk/transform.cpp
+++ b/fityk/transform.cpp
@@ -6,6 +6,8 @@
 #include "transform.h"
 #include "logic.h"
 #include "data.h"
+#include "root/background.hpp"
+#include <stdio.h>
 
 using namespace std;
 
@@ -95,6 +97,94 @@ void shirley_bg(vector<Point> &pp)
         pp[i].y = B[i];
 }
 
+template<class T>
+struct vector_slice {
+    typename vector<T>::iterator begin;
+    typename vector<T>::iterator end;
+};
+
+// Determines the first active slice of the data points.
+// It reads from the vector slice determined by the given
+// iterators.
+// It returns the iterators to the begin of the active
+// slice (first active datapoint) and to the end of the
+// slice (datapoint *after* the last active one).
+struct vector_slice<Point> first_active_slice(const vector<Point> &pp,
+                                              vector<Point>::iterator begin,
+                                              vector<Point>::iterator end)
+{
+    vector<Point>::iterator active_start = begin;
+
+    for (vector<Point>::iterator it = begin; it != end; it++) {
+        if (it->is_active) {
+            active_start = it;
+            break;
+        } else {
+            active_start += 1;
+        }
+    }
+
+    vector<Point>::iterator active_end = active_start;
+
+    for (vector<Point>::iterator it = active_start; it != end; it++) {
+        if (it->is_active) {
+            active_end = it + 1;
+        } else {
+            break;
+        }
+    }
+
+    return {active_start, active_end};
+}
+
+// Calculates the SNIP background iteratively in the
+// given slice of the vector of points.
+// The active points are modified in-place, inactive
+// points are left as they are, assuming that they
+// are already considered background.
+void snip_bg_slice(vector<Point> &pp,
+                   vector<Point>::iterator begin,
+                   vector<Point>::iterator end,
+                   int window_width,
+                   int direction,
+                   int filter_order,
+                   bool smoothing,
+                   bool smooth_window,
+                   bool estimate_compton)
+{
+    vector_slice<Point> slice = first_active_slice(pp, begin, end);
+        
+    if (slice.begin == slice.end) {
+        return;
+    } else {
+        vector<Point> bg_input(slice.begin, slice.end);
+
+        vector<Point> bg_output = ROOT::background(bg_input,
+                                                   window_width,
+                                                   direction,
+                                                   filter_order,
+                                                   smoothing,
+                                                   smooth_window,
+                                                   estimate_compton);
+
+        if (bg_output.size() == bg_input.size()) {
+            for (vector<Point>::iterator it = slice.begin; it != slice.end; it++) {
+                const unsigned int index = it - slice.begin;
+
+                it->y = bg_output[index].y;
+            }
+        }
+        
+        snip_bg_slice(pp, slice.end, pp.end(),
+                      window_width,
+                      direction,
+                      filter_order,
+                      smoothing,
+                      smooth_window,
+                      estimate_compton);
+    }
+}
+
 } // anonymous namespace
 
 namespace fityk {
@@ -202,6 +292,31 @@ void run_data_transform(const DataKeeper& dk, const VMData& vm, Data* data_out)
                 shirley_bg(stackPtr->points);
                 break;
 
+            case OP_DT_SNIP_BG:
+                stackPtr -= 4;
+
+                if ((stackPtr)->is_num) {
+                    throw ExecuteError(op2str(*i) + " is defined only for @n");
+                } else {
+                    const int window_width = (stackPtr+1)->num;
+                    const int direction = (stackPtr+2)->num >= 0 ? ROOT::kBackIncreasingWindow : ROOT::kBackDecreasingWindow;
+                    const int filter_order = (stackPtr+3)->num;
+                    const bool smoothing = false;
+                    const int smooth_window = ROOT::kBackSmoothing3;
+                    const bool estimate_compton = (stackPtr+4)->num > 0 ? true : false;
+
+                    snip_bg_slice(stackPtr->points,
+                                  stackPtr->points.begin(),
+                                  stackPtr->points.end(),
+                                  window_width,
+                                  direction,
+                                  filter_order,
+                                  smoothing,
+                                  smooth_window,
+                                  estimate_compton);
+                }
+                break;
+
             case OP_AND:
                 // do nothing
                 break;
diff --git a/fityk/vm.cpp b/fityk/vm.cpp
index 6e19b734..a6e467da 100644
--- a/fityk/vm.cpp
+++ b/fityk/vm.cpp
@@ -109,7 +109,7 @@ string op2str(int op)
         OP_(FUNC) OP_(SUM_F) OP_(SUM_Z)
         OP_(NUMAREA) OP_(FINDX) OP_(FIND_EXTR)
         OP_(TILDE)
-        OP_(DATASET) OP_(DT_SUM_SAME_X) OP_(DT_AVG_SAME_X) OP_(DT_SHIRLEY_BG)
+        OP_(DATASET) OP_(DT_SUM_SAME_X) OP_(DT_AVG_SAME_X) OP_(DT_SHIRLEY_BG) OP_(DT_SNIP_BG)
         OP_(OPEN_ROUND)  OP_(OPEN_SQUARE)
     }
     return S(op); // unreachable (if all OPs are listed above)
diff --git a/fityk/vm.h b/fityk/vm.h
index 0aa3437c..7532b011 100644
--- a/fityk/vm.h
+++ b/fityk/vm.h
@@ -101,6 +101,7 @@ enum Op
     OP_DT_SUM_SAME_X,
     OP_DT_AVG_SAME_X,
     OP_DT_SHIRLEY_BG,
+    OP_DT_SNIP_BG,
 
     // these two are not VM operators, but are handy to have here,
     // they and are used in implementation of shunting yard algorithm