diff --git a/numass-main/src/main/java/inr/numass/Main.java b/numass-main/src/main/java/inr/numass/Main.java
index 59acfce3..480c2fd1 100644
--- a/numass-main/src/main/java/inr/numass/Main.java
+++ b/numass-main/src/main/java/inr/numass/Main.java
@@ -17,17 +17,17 @@ package inr.numass;
 
 import hep.dataforge.actions.ActionUtils;
 import hep.dataforge.context.Context;
+import hep.dataforge.context.GlobalContext;
 import static hep.dataforge.context.GlobalContext.out;
 import hep.dataforge.data.FileDataFactory;
 import hep.dataforge.fitting.MINUITPlugin;
 import hep.dataforge.io.IOManager;
 import hep.dataforge.io.MetaFileReader;
 import hep.dataforge.meta.Meta;
-import static inr.numass.NumassContext.printDescription;
+import static inr.numass.Numass.printDescription;
 import java.io.File;
 import java.io.FileNotFoundException;
 import java.util.Locale;
-import static java.util.Locale.setDefault;
 import javax.swing.JFileChooser;
 import javax.swing.JFrame;
 import javax.swing.filechooser.FileFilter;
@@ -41,20 +41,6 @@ import org.apache.commons.cli.ParseException;
 import org.slf4j.Logger;
 import org.slf4j.LoggerFactory;
 import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
-import static java.util.Locale.setDefault;
 
 /**
  *
@@ -64,13 +50,13 @@ public class Main {
     public static void main(String[] args) throws Exception {
         setDefault(Locale.US);
 
-        NumassContext context = new NumassContext();
+        Context context = Numass.buildContext();
         context.loadPlugin(new MINUITPlugin());
         run(context, args);
     }
 
     @SuppressWarnings("deprecation")
-    public static void run(NumassContext context, String[] args) throws Exception {
+    public static void run(Context context, String[] args) throws Exception {
         Logger logger = LoggerFactory.getLogger("numass-main");
 
         Options options = prepareOptions();
diff --git a/numass-main/src/main/java/inr/numass/NumassContext.java b/numass-main/src/main/java/inr/numass/Numass.java
similarity index 78%
rename from numass-main/src/main/java/inr/numass/NumassContext.java
rename to numass-main/src/main/java/inr/numass/Numass.java
index a5c9ba38..0ecfd250 100644
--- a/numass-main/src/main/java/inr/numass/NumassContext.java
+++ b/numass-main/src/main/java/inr/numass/Numass.java
@@ -30,27 +30,18 @@ import java.io.PrintWriter;
  *
  * @author Darksnake
  */
-public class NumassContext extends Context {
+public class Numass {
 
-    public NumassContext(Context parent, Meta config) {
-        super(parent, "numass", config);
-        init();
+    public static Context buildContext(Context parent, Meta meta) {
+        Context numassContext = new Context(parent, "numass", meta);
+        GlobalContext.registerContext(numassContext);
+        numassContext.loadPlugin("inr.numass:numass");
+        numassContext.setIO(new NumassIO());
+        return numassContext;
     }
 
-    public NumassContext(Context parent) {
-        super(parent, "numass");
-        init();
-    }
-
-    public NumassContext() {
-        super("numass");
-        init();
-    }
-
-    private void init() {
-        GlobalContext.registerContext(this);
-        loadPlugin("inr.numass:numass");
-        setIO(new NumassIO());
+    public static Context buildContext() {
+        return buildContext(GlobalContext.instance(), Meta.empty());
     }
 
     public static void printDescription(Context context, boolean allowANSI) throws DescriptorException {
diff --git a/numass-main/src/main/java/inr/numass/NumassPlugin.java b/numass-main/src/main/java/inr/numass/NumassPlugin.java
index 0553f1db..f1a281e6 100644
--- a/numass-main/src/main/java/inr/numass/NumassPlugin.java
+++ b/numass-main/src/main/java/inr/numass/NumassPlugin.java
@@ -48,7 +48,6 @@ import inr.numass.models.EmpiricalLossSpectrum;
 import inr.numass.models.ExperimentalVariableLossSpectrum;
 import inr.numass.models.GaussSourceSpectrum;
 import inr.numass.models.GunSpectrum;
-import inr.numass.models.LossCalculator;
 import inr.numass.models.ModularSpectrum;
 import inr.numass.models.NBkgSpectrum;
 import inr.numass.models.RangedNamedSetSpectrum;
@@ -97,11 +96,18 @@ public class NumassPlugin extends BasicPlugin {
     public void detach() {
 
     }
-    
-    private void loadMath(MathPlugin math){
-        
+
+    private void loadMath(MathPlugin math) {
+        math.registerBivariate("numass.trap2016", (Ei, Ef) -> {
+            return 3.92e-5 * FastMath.exp(-(Ei - Ef) / 300d) + 1.97e-4 - 6.818e-9 * Ei;
+        });
+
+        math.registerBivariate("numass.resolutionTail", meta -> {
+            double alpha = meta.getDouble("tailAlpha", 0);
+            double beta = meta.getDouble("tailBeta", 0);
+            return (double E, double U) -> 1 - (E - U) * (alpha + E / 1000d * beta) / 1000d;
+        });
     }
-    
 
     /**
      * Load all numass model factories
@@ -213,7 +219,7 @@ public class NumassPlugin extends BasicPlugin {
         });
 
         manager.addModel("sterile", (context, meta) -> {
-            SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(meta);
+            SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(context, meta);
             NBkgSpectrum spectrum = new NBkgSpectrum(sp);
 
             return new XYModel(spectrum, getAdapter(meta));
@@ -237,20 +243,12 @@ public class NumassPlugin extends BasicPlugin {
                 sp.setCaching(true);
             }
             //Adding trapping energy dependence
-
-            switch (meta.getString("trappingFunction", "default")) {
-                case "run2016":
-                    sp.setTrappingFunction((Ei, Ef) -> {
-                        return 3.92e-5 * FastMath.exp(-(Ei - Ef) / 300d) + 1.97e-4 - 6.818e-9 * Ei;
-                    });
-                    break;
-                default:
-                    //Intercept = 4.95745, B1 = -0.36879, B2 = 0.00827
-                    sp.setTrappingFunction((Ei, Ef) -> LossCalculator.getTrapFunction().value(Ei, Ef) * (4.95745 - 0.36879 * Ei + 0.00827 * Ei * Ei));
+            
+            if(meta.hasValue("transmission.trapping")){
+                BivariateFunction trap = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("transmussion.trapping"));
+                sp.setTrappingFunction(trap);
             }
 
-            context.getReport().report("Using folowing trapping formula: {}",
-                    "6.2e-5 * FastMath.exp(-(Ei - Ef) / 350d) + 1.97e-4 - 6.818e-9 * Ei");
             NBkgSpectrum spectrum = new NBkgSpectrum(sp);
 
             return new XYModel(spectrum, getAdapter(meta));
diff --git a/numass-main/src/main/java/inr/numass/WorkspaceTest.java b/numass-main/src/main/java/inr/numass/WorkspaceTest.java
index 5ad1c17c..49f4f991 100644
--- a/numass-main/src/main/java/inr/numass/WorkspaceTest.java
+++ b/numass-main/src/main/java/inr/numass/WorkspaceTest.java
@@ -23,7 +23,7 @@ public class WorkspaceTest {
         
         String storagepath = "D:\\Work\\Numass\\data\\";
         Workspace workspace = BasicWorkspace.builder()
-                .setContext(new NumassContext())
+                .setContext(Numass.buildContext())
                 .loadData("", new StorageDataFactory(), new MetaBuilder("storage").putValue("path", storagepath))
                 
                 .build();
diff --git a/numass-main/src/main/java/inr/numass/actions/ShowLossSpectrumAction.java b/numass-main/src/main/java/inr/numass/actions/ShowLossSpectrumAction.java
index 34353c80..1979a1c2 100644
--- a/numass-main/src/main/java/inr/numass/actions/ShowLossSpectrumAction.java
+++ b/numass-main/src/main/java/inr/numass/actions/ShowLossSpectrumAction.java
@@ -41,7 +41,7 @@ import hep.dataforge.tables.Table;
 import hep.dataforge.tables.XYAdapter;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.NumassIntegrator;
-import inr.numass.NumassContext;
+import inr.numass.Numass;
 import inr.numass.models.ExperimentalVariableLossSpectrum;
 import inr.numass.models.LossCalculator;
 import java.io.OutputStreamWriter;
diff --git a/numass-main/src/main/java/inr/numass/models/CustomNBkgSpectrum.java b/numass-main/src/main/java/inr/numass/models/CustomNBkgSpectrum.java
index 8aac6298..698e1be2 100644
--- a/numass-main/src/main/java/inr/numass/models/CustomNBkgSpectrum.java
+++ b/numass-main/src/main/java/inr/numass/models/CustomNBkgSpectrum.java
@@ -8,7 +8,7 @@ package inr.numass.models;
 import hep.dataforge.fitting.parametric.ParametricFunction;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.NumassIntegrator;
-import inr.numass.NumassContext;
+import inr.numass.Numass;
 import inr.numass.utils.TritiumUtils;
 import org.apache.commons.math3.analysis.UnivariateFunction;
 
diff --git a/numass-main/src/main/java/inr/numass/models/LossCalculator.java b/numass-main/src/main/java/inr/numass/models/LossCalculator.java
index 3dfc2154..5cf211b6 100644
--- a/numass-main/src/main/java/inr/numass/models/LossCalculator.java
+++ b/numass-main/src/main/java/inr/numass/models/LossCalculator.java
@@ -326,7 +326,7 @@ public class LossCalculator {
     }
 
     public synchronized double getLossValue(int order, double Ei, double Ef) {
-        if (Ei < Ef) {
+        if (Ei - Ef < 5d) {
             return 0;
         } else if (Ei - Ef >= getMargin(order)) {
             return 0;
@@ -377,7 +377,7 @@ public class LossCalculator {
                     return 0;
                 }
             };
-            return integrator.integrate(integrand, 0d, margin);
+            return integrator.integrate(integrand, 5d, margin);
         };
 
         return FunctionCaching.cacheUnivariateFunction(res, 0, margin, 200);
diff --git a/numass-main/src/main/java/inr/numass/models/LossResConvolution.java b/numass-main/src/main/java/inr/numass/models/LossResConvolution.java
index 5bdf0668..59ce44b0 100644
--- a/numass-main/src/main/java/inr/numass/models/LossResConvolution.java
+++ b/numass-main/src/main/java/inr/numass/models/LossResConvolution.java
@@ -16,7 +16,6 @@
 package inr.numass.models;
 
 import inr.numass.NumassIntegrator;
-import inr.numass.NumassContext;
 import org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
 
diff --git a/numass-main/src/main/java/inr/numass/models/TransmissionConvolution.java b/numass-main/src/main/java/inr/numass/models/TransmissionConvolution.java
index 4166c7f8..98f287f9 100644
--- a/numass-main/src/main/java/inr/numass/models/TransmissionConvolution.java
+++ b/numass-main/src/main/java/inr/numass/models/TransmissionConvolution.java
@@ -19,7 +19,7 @@ import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.fitting.parametric.ParametricFunction;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.NumassIntegrator;
-import inr.numass.NumassContext;
+import inr.numass.Numass;
 import org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
 
diff --git a/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java b/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java
index 46fbd5ab..2b92e73c 100644
--- a/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java
+++ b/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java
@@ -22,7 +22,7 @@ import hep.dataforge.maths.integration.UnivariateIntegrator;
 import hep.dataforge.values.NamedValueSet;
 import hep.dataforge.values.ValueProvider;
 import inr.numass.NumassIntegrator;
-import inr.numass.NumassContext;
+import inr.numass.Numass;
 import java.util.List;
 import org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
diff --git a/numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java b/numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java
index f5f5c83a..32e27934 100644
--- a/numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java
+++ b/numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java
@@ -9,7 +9,6 @@ import hep.dataforge.fitting.parametric.AbstractParametricBiFunction;
 import hep.dataforge.meta.Meta;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.models.ResolutionFunction;
-import static java.lang.Double.isNaN;
 import static java.lang.Math.sqrt;
 import org.apache.commons.math3.analysis.BivariateFunction;
 
diff --git a/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java b/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java
index d504b7c7..9e463521 100644
--- a/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java
+++ b/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java
@@ -5,7 +5,9 @@
  */
 package inr.numass.models.sterile;
 
+import hep.dataforge.context.Context;
 import hep.dataforge.fitting.parametric.AbstractParametricBiFunction;
+import hep.dataforge.maths.MathPlugin;
 import hep.dataforge.meta.Meta;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.models.LossCalculator;
@@ -22,21 +24,18 @@ public class NumassTransmission extends AbstractParametricBiFunction {
     private final LossCalculator calculator;
     private final BivariateFunction trapFunc;
 
-//    private BicubicInterpolatingFunction cache;
-//    private double cachedX;
-//    private Meta cacheMeta;
-
-    public NumassTransmission(Meta meta) {
+    public NumassTransmission(Context context, Meta meta) {
         super(list);
         this.calculator = LossCalculator.instance();
-        trapFunc = getTrapFunction(meta.getNodeOrEmpty("trapping"));
-//        if (meta.hasNode("cache")) {
-//            cacheMeta = meta.getNode("cache");
-//        }
+        if (meta.hasValue("trapping")) {
+            trapFunc = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("trapping"));
+        } else {
+            trapFunc = getTrapFunction(context, meta.getNodeOrEmpty("trapping"));
+        }
     }
 
-    private BivariateFunction getTrapFunction(Meta meta) {
-        return LossCalculator.getTrapFunction();
+    private BivariateFunction getTrapFunction(Context context, Meta meta) {
+        return MathPlugin.buildFrom(context).buildBivariateFunction(meta);
     }
 
     @Override
@@ -65,18 +64,6 @@ public class NumassTransmission extends AbstractParametricBiFunction {
         return LossCalculator.instance().getLossProbability(0, getX(eIn, set));
     }
 
-//    private synchronized void setupCache(double X) {
-//        if (this.cachedX != X) {
-//            double cacheLo = cacheMeta.getDouble("lo", 14000);
-//            double cacheHi = cacheMeta.getDouble("hi", 18575);
-//            int numPoints = cacheMeta.getInt("numPoints", 1000);
-//            double[] eIns = GridCalculator.getUniformUnivariateGrid(cacheLo, cacheHi, numPoints);
-//            double[] eOuts = GridCalculator.getUniformUnivariateGrid(cacheLo, cacheHi, numPoints);
-//            double[][] vals = MathUtils.calculateFunction(calculator.getTotalLossBivariateFunction(X), eOuts, eOuts);
-//            this.cachedX = X;
-//            this.cache = new BicubicInterpolator().interpolate(eIns, eOuts, vals);
-//        }
-//    }
     @Override
     public double value(double eIn, double eOut, NamedValueSet set) {
         //calculate X taking into account its energy dependence
diff --git a/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java b/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java
index cba62e23..c3c094b3 100644
--- a/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java
+++ b/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java
@@ -10,6 +10,7 @@ import hep.dataforge.context.GlobalContext;
 import hep.dataforge.description.NodeDef;
 import hep.dataforge.description.ValueDef;
 import hep.dataforge.exceptions.NotDefinedException;
+import hep.dataforge.fitting.parametric.AbstractParametricBiFunction;
 import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.fitting.parametric.ParametricBiFunction;
 import hep.dataforge.maths.integration.UnivariateIntegrator;
@@ -17,14 +18,7 @@ import hep.dataforge.meta.Meta;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.NumassIntegrator;
 import inr.numass.models.FSS;
-import java.util.stream.DoubleStream;
-import org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
-import org.apache.commons.math3.distribution.EnumeratedRealDistribution;
-import org.apache.commons.math3.distribution.RealDistribution;
-import org.apache.commons.math3.random.JDKRandomGenerator;
-import org.apache.commons.math3.random.RandomGenerator;
-import org.apache.commons.math3.random.SynchronizedRandomGenerator;
 
 /**
  * Compact all-in-one model for sterile neutrino spectrum
@@ -39,8 +33,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     private static final String[] list = {"X", "trap", "E0", "mnu2", "msterile2", "U2"};
 
-    private final RandomGenerator rnd;
-    private RealDistribution fssDistribution;
+//    private final RandomGenerator rnd;
+//    private RealDistribution fssDistribution;
     private FSS fss;
 
     /**
@@ -57,21 +51,27 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
      */
     private final ParametricBiFunction resolution;
 
-    private boolean useMC;
+    /**
+     * auxiliary function for trans-res convolution
+     */
+    private final ParametricBiFunction transRes;
+
+//    private boolean useMC;
     private boolean fast;
 
     public SterileNeutrinoSpectrum(Context context, Meta configuration) {
         super(list);
-        rnd = new SynchronizedRandomGenerator(new JDKRandomGenerator());
+//        rnd = new SynchronizedRandomGenerator(new JDKRandomGenerator());
         if (configuration.hasValue("fssFile")) {
             fss = new FSS(context.io().getFile(configuration.getString("fssFile")));
-            fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs());
+//            fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs());
         }
 
-        transmission = new NumassTransmission(configuration.getNodeOrEmpty("transmission"));
+        transmission = new NumassTransmission(context, configuration.getNodeOrEmpty("transmission"));
         resolution = new NumassResolution(configuration.getNode("resolution", Meta.empty()));
-        this.useMC = configuration.getBoolean("useMC", false);
+//        this.useMC = configuration.getBoolean("useMC", false);
         this.fast = configuration.getBoolean("fast", true);
+        transRes = new TransRes();
     }
 
     public SterileNeutrinoSpectrum(Meta configuration) {
@@ -89,10 +89,10 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             case "msterile2":
             case "mnu2":
             case "E0":
-                return integrate(u, source.derivative(parName), transmission, resolution, set);
+                return integrate(u, source.derivative(parName), transRes, set);
             case "X":
             case "trap":
-                return integrate(u, source, transmission, resolution.derivative(parName), set);
+                return integrate(u, source, transRes.derivative(parName), set);
             default:
                 throw new NotDefinedException();
         }
@@ -100,39 +100,94 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     @Override
     public double value(double u, NamedValueSet set) {
-        if (useDirect()) {
-            return integrateDirect(u, source, transmission, resolution, set);
-        } else {
-            return integrate(u, source, transmission, resolution, set);
-        }
-    }
-
-    private int numCalls(double u) {
-        return 100000;
-    }
-
-    private boolean useDirect() {
-        return !useMC;
+        return integrate(u, source, transRes, set);
     }
 
+//    private int numCalls(double u) {
+//        return 100000;
+//    }
+//
+//    private boolean useDirect() {
+//        return !useMC;
+//    }
     @Override
     public boolean providesDeriv(String name) {
         return source.providesDeriv(name) && transmission.providesDeriv(name) && resolution.providesDeriv(name);
     }
 
+//    /**
+//     * Random E generator
+//     *
+//     * @param a
+//     * @param b
+//     * @return
+//     */
+//    private double rndE(double a, double b) {
+//        return rnd.nextDouble() * (b - a) + a;
+//    }
+//
+//    private double integrate(
+//            double u,
+//            ParametricBiFunction sourceFunction,
+//            ParametricBiFunction transmissionFunction,
+//            ParametricBiFunction resolutionFunction,
+//            NamedValueSet set) {
+//        if (useDirect()) {
+//            return integrateDirect(u, sourceFunction, transmissionFunction, resolutionFunction, set);
+//        } else {
+//            return integrateRandom(u, sourceFunction, transmissionFunction, resolutionFunction, set);
+//        }
+//    }
+//    /**
+//     * Monte-Carlo integration of spectrum
+//     *
+//     * @param u
+//     * @param sourceFunction
+//     * @param transmissionFunction
+//     * @param resolutionFunction
+//     * @param set
+//     * @return
+//     */
+//    private double integrateRandom(
+//            double u,
+//            ParametricBiFunction sourceFunction,
+//            ParametricBiFunction transmissionFunction,
+//            ParametricBiFunction resolutionFunction,
+//            NamedValueSet set) {
+//
+//        int num = numCalls(u);
+//        double eMax = set.getDouble("E0") + 5d;
+//        if (u > eMax) {
+//            return 0;
+//        }
+//
+//        double sum = DoubleStream.generate(() -> {
+//            // generate final state
+//            double fs;
+//            if (fssDistribution != null) {
+//                fs = fssDistribution.sample();
+//            } else {
+//                fs = 0;
+//            }
+//
+//            double eIn = rndE(u, eMax);
+//
+//            double eOut = rndE(u, eIn);
+//
+//            double res = sourceFunction.value(fs, eIn, set)
+//                    * transmissionFunction.value(eIn, eOut, set)
+//                    * resolutionFunction.value(eOut, u, set);
+//
+//            if (Double.isNaN(res)) {
+//                throw new Error();
+//            }
+//            return res;
+//        }).parallel().limit(num).sum();
+//        //triangle surface
+//        return Math.pow(eMax - u, 2d) / 2d * sum / num;
+//    }
     /**
-     * Random E generator
-     *
-     * @param a
-     * @param b
-     * @return
-     */
-    private double rndE(double a, double b) {
-        return rnd.nextDouble() * (b - a) + a;
-    }
-
-    /**
-     * Monte-Carlo integration of spectrum
+     * Direct Gauss-Legandre integration
      *
      * @param u
      * @param sourceFunction
@@ -144,57 +199,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
     private double integrate(
             double u,
             ParametricBiFunction sourceFunction,
-            ParametricBiFunction transmissionFunction,
-            ParametricBiFunction resolutionFunction,
-            NamedValueSet set) {
-
-        int num = numCalls(u);
-        double eMax = set.getDouble("E0") + 5d;
-        if (u > eMax) {
-            return 0;
-        }
-
-        double sum = DoubleStream.generate(() -> {
-            // generate final state
-            double fs;
-            if (fssDistribution != null) {
-                fs = fssDistribution.sample();
-            } else {
-                fs = 0;
-            }
-
-            double eIn = rndE(u, eMax);
-
-            double eOut = rndE(u, eIn);
-
-            double res = sourceFunction.value(fs, eIn, set)
-                    * transmissionFunction.value(eIn, eOut, set)
-                    * resolutionFunction.value(eOut, u, set);
-
-            if (Double.isNaN(res)) {
-                throw new Error();
-            }
-            return res;
-        }).parallel().limit(num).sum();
-        //triangle surface
-        return Math.pow(eMax - u, 2d) / 2d * sum / num;
-    }
-
-    /**
-     * Direct Gauss-Legandre integration
-     *
-     * @param u
-     * @param sourceFunction
-     * @param transmissionFunction
-     * @param resolutionFunction
-     * @param set
-     * @return
-     */
-    private double integrateDirect(
-            double u,
-            ParametricBiFunction sourceFunction,
-            ParametricBiFunction transmissionFunction,
-            ParametricBiFunction resolutionFunction,
+            ParametricBiFunction transResFunction,
             NamedValueSet set) {
 
         double eMax = set.getDouble("E0") + 5d;
@@ -216,14 +221,55 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             fsSource = (eIn) -> sourceFunction.value(0, eIn, set);
         }
 
-        BivariateFunction transRes = (eIn, usp) -> {
-            UnivariateFunction integrand = (eOut) -> transmissionFunction.value(eIn, eOut, set) * resolutionFunction.value(eOut, usp, set);
+        UnivariateIntegrator integrator;
+        if (fast && eMax - u < 500) {
+            integrator = NumassIntegrator.getFastInterator();
+        } else {
+            integrator = NumassIntegrator.getDefaultIntegrator();
+        }
+
+        return integrator.integrate(eIn -> fsSource.value(eIn) * transResFunction.value(eIn, u, set), u, eMax);
+    }
+
+    private class TransRes extends AbstractParametricBiFunction {
+
+        public TransRes() {
+            super(new String[]{"X", "trap"});
+        }
+
+        @Override
+        public boolean providesDeriv(String name) {
+            return true;
+        }
+
+        @Override
+        public double derivValue(String parName, double eIn, double u, NamedValueSet set) {
+            switch (parName) {
+                case "X":
+                    //TODO implement p0 derivative
+                    throw new NotDefinedException();
+                case "trap":
+                    return lossRes(transmission.derivative(parName), eIn, u, set);
+                default:
+                    return super.derivValue(parName, eIn, u, set);
+            }
+        }
+
+        @Override
+        public double value(double eIn, double u, NamedValueSet set) {
+
+            double p0 = NumassTransmission.p0(eIn, set);
+            return p0 * resolution.value(eIn, u, set) + lossRes(transmission, eIn, u, set);
+        }
+
+        private double lossRes(ParametricBiFunction transFunc, double eIn, double u, NamedValueSet set) {
+            UnivariateFunction integrand = (eOut) -> transFunc.value(eIn, eOut, set) * resolution.value(eOut, u, set);
 
             double border = u + 30;
-            double firstPart = NumassIntegrator.getFastInterator().integrate(integrand, usp, Math.min(eIn, border));
+            double firstPart = NumassIntegrator.getFastInterator().integrate(integrand, u, Math.min(eIn, border));
             double secondPart;
             if (eIn > border) {
-                if(fast){
+                if (fast) {
                     secondPart = NumassIntegrator.getFastInterator().integrate(integrand, border, eIn);
                 } else {
                     secondPart = NumassIntegrator.getDefaultIntegrator().integrate(integrand, border, eIn);
@@ -232,21 +278,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
                 secondPart = 0;
             }
             return firstPart + secondPart;
-        };
-
-        UnivariateFunction integrand = (eIn) -> {
-            double p0 = NumassTransmission.p0(eIn, set);
-            return fsSource.value(eIn) * (p0 * resolutionFunction.value(eIn, u, set) + transRes.value(eIn, u));
-        };
-
-        UnivariateIntegrator integrator;
-        if (fast && eMax - u < 500) {
-            integrator = NumassIntegrator.getFastInterator();
-        } else {
-            integrator = NumassIntegrator.getDefaultIntegrator();
         }
 
-        return integrator.integrate(integrand, u, eMax);
     }
 
 }
diff --git a/numass-main/src/main/java/inr/numass/models/sterile/TestModels.java b/numass-main/src/main/java/inr/numass/models/sterile/TestModels.java
new file mode 100644
index 00000000..3f049862
--- /dev/null
+++ b/numass-main/src/main/java/inr/numass/models/sterile/TestModels.java
@@ -0,0 +1,102 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package inr.numass.models.sterile;
+
+import hep.dataforge.context.Context;
+import hep.dataforge.fitting.ParamSet;
+import hep.dataforge.fitting.parametric.ParametricFunction;
+import hep.dataforge.maths.MathPlugin;
+import hep.dataforge.meta.Meta;
+import hep.dataforge.meta.MetaBuilder;
+import inr.numass.Numass;
+import inr.numass.models.BetaSpectrum;
+import inr.numass.models.ModularSpectrum;
+import inr.numass.models.NBkgSpectrum;
+import inr.numass.models.RangedNamedSetSpectrum;
+import inr.numass.models.ResolutionFunction;
+import org.apache.commons.math3.analysis.BivariateFunction;
+
+/**
+ *
+ * @author Alexander Nozik
+ */
+public class TestModels {
+
+    /**
+     * @param args the command line arguments
+     */
+    public static void main(String[] args) {
+        Context context = Numass.buildContext();
+        /*
+        		<model modelName="sterile" fssFile="FS.txt">
+			<resolution width = "1.22e-4" tailAlpha="1e-2"/>
+			<transmission trapping="numass.trap2016"/>
+		</model>
+         */
+        Meta meta = new MetaBuilder("model")
+                .putNode(new MetaBuilder("resolution")
+                        .putValue("width", 1.22e-4)
+                        .putValue("tailAlpha", 1e-2)
+                )
+                .putNode(new MetaBuilder("transmission")
+                        .putValue("trapping", "numass.trap2016")
+                );
+        ParametricFunction oldFunc = oldModel(context, meta);
+        ParametricFunction newFunc = newModel(context, meta);
+
+        ParamSet allPars = new ParamSet()
+                .setPar("N", 7e+05, 1.8e+03, 0d, Double.POSITIVE_INFINITY)
+                .setPar("bkg", 1d, 0.050)
+                .setPar("E0", 18575d, 1.4)
+                .setPar("mnu2", 0d, 1d)
+                .setPar("msterile2", 1000d * 1000d, 0)
+                .setPar("U2", 0.0, 1e-4, -1d, 1d)
+                .setPar("X", 0.04, 0.01, 0d, Double.POSITIVE_INFINITY)
+                .setPar("trap", 1, 0.01, 0d, Double.POSITIVE_INFINITY);
+
+        for (double u = 14000; u < 18600; u += 100) {
+//            double oldVal = oldFunc.value(u, allPars);
+//            double newVal = newFunc.value(u, allPars);
+            double oldVal = oldFunc.derivValue("trap", u, allPars);
+            double newVal = newFunc.derivValue("trap", u, allPars);
+            System.out.printf("%f\t%g\t%g\t%g%n", u, oldVal, newVal, 1d - oldVal / newVal);
+        }
+    }
+
+    private static ParametricFunction oldModel(Context context, Meta meta) {
+        double A = meta.getDouble("resolution", meta.getDouble("resolution.width", 8.3e-5));//8.3e-5
+        double from = meta.getDouble("from", 13900d);
+        double to = meta.getDouble("to", 18700d);
+        context.getReport().report("Setting up tritium model with real transmission function");
+        BivariateFunction resolutionTail;
+        if (meta.hasValue("resolution.tailAlpha")) {
+            resolutionTail = ResolutionFunction.getAngledTail(meta.getDouble("resolution.tailAlpha"), meta.getDouble("resolution.tailBeta", 0));
+        } else {
+            resolutionTail = ResolutionFunction.getRealTail();
+        }
+        //RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt"));
+        RangedNamedSetSpectrum beta = new BetaSpectrum();
+        ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, resolutionTail), from, to);
+        if (meta.getBoolean("caching", false)) {
+            context.getReport().report("Caching turned on");
+            sp.setCaching(true);
+        }
+        //Adding trapping energy dependence
+
+        if (meta.hasValue("transmission.trapping")) {
+            BivariateFunction trap = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("transmission.trapping"));
+            sp.setTrappingFunction(trap);
+        }
+
+        return new NBkgSpectrum(sp);
+    }
+
+    private static ParametricFunction newModel(Context context, Meta meta) {
+        SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(context, meta);
+        return new NBkgSpectrum(sp);
+    }
+
+}
diff --git a/numass-main/src/main/java/inr/numass/workbench/NumassWorkbenchController.java b/numass-main/src/main/java/inr/numass/workbench/NumassWorkbenchController.java
index b2be8299..ef57a5c2 100644
--- a/numass-main/src/main/java/inr/numass/workbench/NumassWorkbenchController.java
+++ b/numass-main/src/main/java/inr/numass/workbench/NumassWorkbenchController.java
@@ -30,7 +30,6 @@ import hep.dataforge.meta.MetaBuilder;
 import hep.dataforge.plots.PlotFrame;
 import hep.dataforge.plots.PlotHolder;
 import hep.dataforge.plots.PlotsPlugin;
-import hep.dataforge.utils.MetaFactory;
 import hep.dataforge.values.Value;
 import inr.numass.NumassIO;
 import inr.numass.NumassProperties;
@@ -60,6 +59,7 @@ import javafx.scene.control.TitledPane;
 import javafx.scene.control.ToggleButton;
 import javafx.stage.FileChooser;
 import org.controlsfx.control.StatusBar;
+import hep.dataforge.utils.ContextMetaFactory;
 
 /**
  * FXML Controller class
@@ -69,7 +69,7 @@ import org.controlsfx.control.StatusBar;
 public class NumassWorkbenchController implements Initializable, StagePaneHolder, ActionStateListener, PlotHolder {
 
     Context parentContext;
-    MetaFactory<Context> contextFactory;
+    ContextMetaFactory<Context> contextFactory;
 
     List<MetaEditor> actionEditors = new ArrayList<>();
     MetaEditor dataEditor;
@@ -402,7 +402,7 @@ public class NumassWorkbenchController implements Initializable, StagePaneHolder
         this.parentContext = parentContext;
     }
 
-    public void setContextFactory(MetaFactory<Context> contextFactory) {
+    public void setContextFactory(ContextMetaFactory<Context> contextFactory) {
         this.contextFactory = contextFactory;
     }
 
diff --git a/numass-main/src/main/java/inr/numass/workbench/Workbench.java b/numass-main/src/main/java/inr/numass/workbench/Workbench.java
index fe776ddb..3ac34c62 100644
--- a/numass-main/src/main/java/inr/numass/workbench/Workbench.java
+++ b/numass-main/src/main/java/inr/numass/workbench/Workbench.java
@@ -6,7 +6,7 @@
 package inr.numass.workbench;
 
 import hep.dataforge.context.GlobalContext;
-import inr.numass.NumassContext;
+import inr.numass.Numass;
 import java.io.IOException;
 import java.text.ParseException;
 import javafx.application.Application;
@@ -30,7 +30,7 @@ public class Workbench extends Application {
         Scene scene = new Scene(parent, 800, 600);
 
         NumassWorkbenchController controller = loader.getController();
-        controller.setContextFactory(NumassContext::new);
+        controller.setContextFactory(Numass::buildContext);
 
         primaryStage.setTitle("Numass workbench");
         primaryStage.setScene(scene);
diff --git a/numass-main/src/test/java/inr/numass/MainTest.java b/numass-main/src/test/java/inr/numass/MainTest.java
new file mode 100644
index 00000000..a21b2e90
--- /dev/null
+++ b/numass-main/src/test/java/inr/numass/MainTest.java
@@ -0,0 +1,60 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package inr.numass;
+
+import hep.dataforge.context.Context;
+import org.apache.commons.cli.CommandLine;
+import org.junit.Test;
+import static org.junit.Assert.*;
+
+/**
+ *
+ * @author Alexander Nozik
+ */
+public class MainTest {
+    
+    public MainTest() {
+    }
+
+    /**
+     * Test of main method, of class Main.
+     */
+    @Test
+    public void testMain() throws Exception {
+        System.out.println("main");
+        String[] args = null;
+        Main.main(args);
+        // TODO review the generated test code and remove the default call to fail.
+        fail("The test case is a prototype.");
+    }
+
+    /**
+     * Test of run method, of class Main.
+     */
+    @Test
+    public void testRun() throws Exception {
+        System.out.println("run");
+        Context context = null;
+        String[] args = null;
+        Main.run(context, args);
+        // TODO review the generated test code and remove the default call to fail.
+        fail("The test case is a prototype.");
+    }
+
+    /**
+     * Test of applyCLItoContext method, of class Main.
+     */
+    @Test
+    public void testApplyCLItoContext() throws Exception {
+        System.out.println("applyCLItoContext");
+        CommandLine line = null;
+        Context context = null;
+        Main.applyCLItoContext(line, context);
+        // TODO review the generated test code and remove the default call to fail.
+        fail("The test case is a prototype.");
+    }
+    
+}
diff --git a/numass-main/src/test/java/inr/numass/NumassContextTest.java b/numass-main/src/test/java/inr/numass/NumassContextTest.java
deleted file mode 100644
index 3a801db9..00000000
--- a/numass-main/src/test/java/inr/numass/NumassContextTest.java
+++ /dev/null
@@ -1,32 +0,0 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-package inr.numass;
-
-import inr.numass.storage.SetDirectionUtility;
-import org.junit.Test;
-
-/**
- *
- * @author Alexander Nozik
- */
-public class NumassContextTest {
-    
-    public NumassContextTest() {
-    }
-
-    /**
-     * Test of close method, of class NumassContext.
-     */
-    @Test
-    public void testClose() throws Exception {
-        System.out.println("close");
-        NumassContext instance = new NumassContext();
-        instance.close();
-        assert SetDirectionUtility.cacheFile(instance).exists();
-        
-    }
-    
-}
diff --git a/numass-main/src/test/java/inr/numass/NumassTest.java b/numass-main/src/test/java/inr/numass/NumassTest.java
new file mode 100644
index 00000000..db12649d
--- /dev/null
+++ b/numass-main/src/test/java/inr/numass/NumassTest.java
@@ -0,0 +1,33 @@
+/*
+ * To change this license header, choose License Headers in Project Properties.
+ * To change this template file, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package inr.numass;
+
+import hep.dataforge.context.Context;
+import hep.dataforge.maths.MathPlugin;
+import hep.dataforge.plots.PlotsPlugin;
+import org.junit.Test;
+
+/**
+ *
+ * @author Alexander Nozik
+ */
+public class NumassTest {
+    
+    public NumassTest() {
+    }
+
+    /**
+     * Test of buildContext method, of class Numass.
+     */
+    @Test
+    public void testBuildContext() {
+        
+        Context context = Numass.buildContext();
+        MathPlugin.buildFrom(context);
+        PlotsPlugin.buildFrom(context);
+    }
+
+}