[no commit message]

numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy

@@ -16,7 +16,7 @@
 package inr.numass.scripts
 
 
-import hep.dataforge.datafitter.ParamSet
+import hep.dataforge.fitting.ParamSet
 import hep.dataforge.maths.integration.RiemanIntegrator
 import hep.dataforge.maths.integration.UnivariateIntegrator
 import hep.dataforge.plots.PlotFrame

numass-main/src/main/groovy/inr/numass/scripts/Loss2014.groovy

@@ -16,7 +16,7 @@
 package inr.numass.scripts
 
 import hep.dataforge.data.DataNode
-import hep.dataforge.datafitter.FitTaskResult
+import hep.dataforge.fitting.FitTaskResult
 import inr.numass.Main
 import inr.numass.NumassContext
 

numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy

@@ -17,13 +17,13 @@ package inr.numass.scripts;
 
 import hep.dataforge.context.GlobalContext;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.MINUITPlugin
+import hep.dataforge.fitting.MINUITPlugin
 
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
-import hep.dataforge.likelihood.BayesianManager
+import hep.dataforge.fitting.likelihood.BayesianManager
 import inr.numass.data.SpectrumDataAdapter
 import inr.numass.models.BetaSpectrum;
 import inr.numass.models.ModularSpectrum;

numass-main/src/main/groovy/inr/numass/scripts/ResolutionTest.groovy

@@ -18,11 +18,11 @@ package inr.numass.scripts;
 import hep.dataforge.context.GlobalContext;
 import static hep.dataforge.context.GlobalContext.out;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.FitTask;
+import hep.dataforge.fitting.FitTask;
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
 import hep.dataforge.exceptions.NamingException;
 import inr.numass.data.SpectrumDataAdapter;
 import inr.numass.data.SpectrumGenerator;

numass-main/src/main/groovy/inr/numass/scripts/SignificanceTest.groovy

@@ -17,7 +17,7 @@ package inr.numass.scripts;
 
 import hep.dataforge.meta.MetaBuilder;
 import hep.dataforge.context.GlobalContext;
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
 import inr.numass.data.SpectrumInformation;
 import inr.numass.models.ModularSpectrum;
 import inr.numass.models.BetaSpectrum;

numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy

@@ -18,13 +18,13 @@ package inr.numass.scripts;
 import hep.dataforge.context.GlobalContext;
 import static hep.dataforge.context.GlobalContext.out;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.FitTask;
+import hep.dataforge.fitting.FitTask;
-import hep.dataforge.datafitter.MINUITPlugin
+import hep.dataforge.fitting.MINUITPlugin
 
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
 import hep.dataforge.exceptions.NamingException;
 import hep.dataforge.exceptions.PackFormatException;
 import inr.numass.data.SpectrumDataAdapter;

numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy

@@ -16,9 +16,9 @@
 package inr.numass.scripts;
 
 import hep.dataforge.context.GlobalContext;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
 import hep.dataforge.exceptions.NamingException;
 import hep.dataforge.io.FittingIOUtils;
 import inr.numass.data.SpectrumDataAdapter;

numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy

@@ -18,13 +18,14 @@ package inr.numass.scripts;
 import hep.dataforge.context.GlobalContext;
 import static hep.dataforge.context.GlobalContext.out;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.FitTask;
+import hep.dataforge.fitting.FitTask;
-import hep.dataforge.datafitter.MINUITPlugin
+import hep.dataforge.fitting.MINUITPlugin
 
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
+import hep.dataforge.fitting.parametric.ParametricFunction
 import hep.dataforge.exceptions.NamingException;
 import hep.dataforge.exceptions.PackFormatException;
 import inr.numass.data.SpectrumDataAdapter;
@@ -32,6 +33,7 @@ import inr.numass.data.SpectrumGenerator;
 import inr.numass.models.BetaSpectrum
 import inr.numass.models.ModularSpectrum;
 import inr.numass.models.NBkgSpectrum;
+import inr.numass.models.sterile.SterileNeutrinoSpectrum
 import inr.numass.utils.DataModelUtils;
 import hep.dataforge.plotfit.PlotFitResultAction;
 import java.io.FileNotFoundException;
@@ -48,8 +50,9 @@ import hep.dataforge.io.FittingIOUtils
 
 setDefault(Locale.US);
 
-ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d);
+//ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d);
-beta.setCaching(false);
+//beta.setCaching(false);
+ParametricFunction beta = new SterileNeutrinoSpectrum();
 
 NBkgSpectrum spectrum = new NBkgSpectrum(beta);
 XYModel model = new XYModel(spectrum, new SpectrumDataAdapter());

numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy

@@ -18,13 +18,13 @@ package inr.numass.scripts;
 import hep.dataforge.context.GlobalContext;
 import static hep.dataforge.context.GlobalContext.out;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.FitTask;
+import hep.dataforge.fitting.FitTask;
-import hep.dataforge.datafitter.MINUITPlugin
+import hep.dataforge.fitting.MINUITPlugin
 
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
 import hep.dataforge.exceptions.NamingException;
 import hep.dataforge.exceptions.PackFormatException;
 import inr.numass.data.SpectrumDataAdapter;

numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy

@@ -18,13 +18,13 @@ package inr.numass.scripts;
 import hep.dataforge.context.GlobalContext;
 import static hep.dataforge.context.GlobalContext.out;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.FitTask;
+import hep.dataforge.fitting.FitTask;
-import hep.dataforge.datafitter.MINUITPlugin
+import hep.dataforge.fitting.MINUITPlugin
 
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
 import hep.dataforge.exceptions.NamingException;
 import hep.dataforge.exceptions.PackFormatException;
 import inr.numass.data.SpectrumDataAdapter;

numass-main/src/main/groovy/inr/numass/scripts/TestExperimentalVariableLossSpectrum.groovy

@@ -22,7 +22,7 @@ import inr.numass.models.ExperimentalVariableLossSpectrum
 import org.apache.commons.math3.analysis.UnivariateFunction
 
 import inr.numass.NumassContext
-import hep.dataforge.datafitter.ParamSet
+import hep.dataforge.fitting.ParamSet
 import hep.dataforge.io.PrintFunction
 
 

numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy

@@ -18,11 +18,11 @@ package inr.numass.scripts;
 import hep.dataforge.context.GlobalContext;
 import hep.dataforge.data.DataSet;
 import hep.dataforge.tables.ListTable;
-import hep.dataforge.datafitter.FitManager;
+import hep.dataforge.fitting.FitManager;
-import hep.dataforge.datafitter.FitState;
+import hep.dataforge.fitting.FitState;
-import hep.dataforge.datafitter.ParamSet;
+import hep.dataforge.fitting.ParamSet;
-import hep.dataforge.datafitter.models.XYModel;
+import hep.dataforge.fitting.models.XYModel;
-import hep.dataforge.likelihood.BayesianManager
+import hep.dataforge.fitting.likelihood.BayesianManager
 import static hep.dataforge.maths.RandomUtils.setSeed;
 import inr.numass.data.SpectrumGenerator;
 import inr.numass.models.BetaSpectrum

numass-main/src/main/java/inr/numass/models/BetaSpectrum.java

@@ -20,9 +20,9 @@ import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.values.NamedValueSet;
 import hep.dataforge.values.ValueProvider;
 import java.io.File;
+import static java.lang.Math.abs;
 import static java.lang.Math.exp;
 import static java.lang.Math.sqrt;
-import static java.lang.Math.abs;
 
 /**
  *
@@ -130,12 +130,9 @@ public class BetaSpectrum extends AbstractParametricFunction implements RangedNa
             res += fss.getP(i) * this.derivRootsterile(name, E + fss.getE(i), pars);
         }
         return res;
-
-//        return this.derivRootsterile(name, E, pars);
     }
 
     double factor(double E) {
-//        return K*(4.632e10+2.21e6*E);
         return K * pfactor(E);
     }
 

numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java

@@ -0,0 +1,201 @@
+/*
+ * 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.exceptions.NotDefinedException;
+import hep.dataforge.fitting.parametric.AbstractParametricBiFunction;
+import hep.dataforge.values.NamedValueSet;
+import static java.lang.Math.abs;
+import static java.lang.Math.exp;
+import static java.lang.Math.sqrt;
+
+/**
+ * A bi-function for beta-spectrum calculation taking energy and final state as input.
+ * @author Alexander Nozik <altavir@gmail.com>
+ */
+public class NumassBeta extends AbstractParametricBiFunction {
+
+    private static final double K = 1E-23;
+    private static final String[] list = {"E0", "mnu2", "msterile2", "U2"};
+
+    public NumassBeta() {
+        super(list);
+    }
+
+    private double derivRoot(int n, double E0, double mnu2, double E) throws NotDefinedException {
+        double D = E0 - E;//E0-E
+        double res;
+        if (D == 0) {
+            return 0;
+        }
+
+        if (mnu2 >= 0) {
+            if (E >= (E0 - sqrt(mnu2))) {
+                return 0;
+            }
+            double bare = sqrt(D * D - mnu2);
+            switch (n) {
+                case 0:
+                    res = factor(E) * (2 * D * D - mnu2) / bare;
+                    break;
+                case 1:
+                    res = -factor(E) * 0.5 * D / bare;
+                    break;
+                default:
+                    return 0;
+            }
+        } else {
+            double mu = sqrt(-0.66 * mnu2);
+            if (E >= (E0 + mu)) {
+                return 0;
+            }
+            double root = sqrt(Math.max(D * D - mnu2, 0));
+            double exp = exp(-1 - D / mu);
+            switch (n) {
+                case 0:
+                    res = factor(E) * (D * (D + mu * exp) / root + root * (1 - exp));
+                    break;
+                case 1:
+                    res = factor(E) * (-(D + mu * exp) / root * 0.5 - root * exp * (1 + D / mu) / 3 / mu);
+                    break;
+                default:
+                    return 0;
+            }
+        }
+
+        return res;
+    }
+
+    /**
+     * Derivative of spectrum with sterile neutrinos
+     * @param name
+     * @param E
+     * @param E0
+     * @param pars
+     * @return
+     * @throws NotDefinedException 
+     */
+    private double derivRootsterile(String name, double E, double E0, NamedValueSet pars) throws NotDefinedException {
+        double mnu2 = getParameter("mnu2", pars);
+        double mst2 = getParameter("msterile2", pars);
+        double u2 = getParameter("U2", pars);
+
+        switch (name) {
+            case "E0":
+                if (u2 == 0) {
+                    return derivRoot(0, E0, mnu2, E);
+                }
+                return u2 * derivRoot(0, E0, mst2, E) + (1 - u2) * derivRoot(0, E0, mnu2, E);
+            case "mnu2":
+                return (1 - u2) * derivRoot(1, E0, mnu2, E);
+            case "msterile2":
+                if (u2 == 0) {
+                    return 0;
+                }
+                return u2 * derivRoot(1, E0, mst2, E);
+            case "U2":
+                return root(E0, mst2, E) - root(E0, mnu2, E);
+            default:
+                return 0;
+        }
+
+    }
+
+    private double factor(double E) {
+        return K * pfactor(E);
+    }
+
+    private double pfactor(double E) {
+        double me = 0.511006E6;
+        double Etot = E + me;
+        double pe = sqrt(E * (E + 2d * me));
+        double ve = pe / Etot;
+        double yfactor = 2d * 2d * 1d / 137.039 * Math.PI;
+        double y = yfactor / ve;
+        double Fn = y / abs(1d - exp(-y));
+        double Fermi = Fn * (1.002037 - 0.001427 * ve);
+        double res = Fermi * pe * Etot;
+        return res;
+    }
+
+    @Override
+    public boolean providesDeriv(String name) {
+        return true;
+    }
+
+    /**
+     * Bare beta spectrum with Mintz negative mass correction
+     * @param E0
+     * @param mnu2
+     * @param E
+     * @return 
+     */
+    private double root(double E0, double mnu2, double E) {
+        /*чистый бета-спектр*/
+        double delta = E0 - E;//E0-E
+        double res;
+        double bare = factor(E) * delta * sqrt(Math.max(delta * delta - mnu2, 0));
+        if (delta == 0) {
+            return 0;
+        }
+        if (mnu2 >= 0) {
+            res = Math.max(bare, 0);
+        } else {
+            if (delta + 0.812 * sqrt(-mnu2) <= 0) {
+                return 0;              //sqrt(0.66)
+            }
+            double aux = sqrt(-mnu2 * 0.66) / delta;
+            res = Math.max(bare * (1 + aux * exp(-1 - 1 / aux)), 0);
+        }
+        return res;
+    }
+
+    /**
+     * beta-spectrum with sterile neutrinos
+     * @param E
+     * @param E0
+     * @param pars
+     * @return 
+     */
+    private double rootsterile(double E, double E0, NamedValueSet pars) {
+        double mnu2 = getParameter("mnu2", pars);
+        double mst2 = getParameter("msterile2", pars);
+        double u2 = getParameter("U2", pars);
+
+        if (u2 == 0) {
+            return root(E0, mnu2, E);
+        }
+        return u2 * root(E0, mst2, E) + (1 - u2) * root(E0, mnu2, E);
+        // P(rootsterile)+ (1-P)root
+    }
+
+    @Override
+    protected double getDefaultParameter(String name) {
+        switch (name) {
+            case "mnu2":
+                return 0;
+            case "U2":
+                return 0;
+            case "msterile2":
+                return 0;
+            default:
+                return super.getDefaultParameter(name);
+        }
+    }
+
+    @Override
+    public double derivValue(String parName, double fs, double eIn, NamedValueSet pars) {
+        double E0 = getParameter("E0", pars);
+        return derivRootsterile(parName,eIn, E0 - fs, pars);
+    }
+
+    @Override
+    public double value(double fs, double eIn, NamedValueSet pars) {
+        double E0 = getParameter("E0", pars);
+        return rootsterile(eIn, E0 - fs, pars);
+    }
+
+}

numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java

@@ -0,0 +1,78 @@
+/*
+ * 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.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;
+
+/**
+ *
+ * @author Alexander Nozik <altavir@gmail.com>
+ */
+public class NumassResolution extends AbstractParametricBiFunction {
+
+    private static final String[] list = {}; //leaving 
+
+    private final double resA;
+    private double resB = Double.NaN;
+    private BivariateFunction tailFunction = ResolutionFunction.getConstantTail();
+
+    public NumassResolution(Meta meta) {
+        super(list);
+        this.resA = meta.getDouble("A", 8.3e-5);
+        this.resB = meta.getDouble("B", 0);
+        if (meta.hasValue("tailAlpha")) {
+            //add polinomial function here
+            tailFunction = ResolutionFunction.getAngledTail(meta.getDouble("tailAlpha"), meta.getDouble("tailBeta", 0));
+        } else {
+            tailFunction = ResolutionFunction.getConstantTail();
+        }
+    }
+
+    @Override
+    public double derivValue(String parName, double x, double y, NamedValueSet set) {
+        return 0;
+    }
+
+    private double getValueFast(double E, double U) {
+        double delta = resA * E;
+        if (E - U < 0) {
+            return 0;
+        } else if (E - U > delta) {
+            return tailFunction.value(E, U);
+        } else {
+            return (E - U) / delta;
+        }
+    }
+
+    @Override
+    public boolean providesDeriv(String name) {
+        return true;
+    }
+
+    @Override
+    public double value(double E, double U, NamedValueSet set) {
+        assert resA > 0;
+        if (isNaN(resB)) {
+            return this.getValueFast(E, U);
+        }
+        assert resB > 0;
+        double delta = resA * E;
+        if (E - U < 0) {
+            return 0;
+        } else if (E - U > delta) {
+            return tailFunction.value(E, U);
+        } else {
+            return (1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB));
+        }
+    }
+
+}

numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java

@@ -0,0 +1,63 @@
+/*
+ * 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.fitting.parametric.AbstractParametricBiFunction;
+import hep.dataforge.values.NamedValueSet;
+import inr.numass.models.LossCalculator;
+import org.apache.commons.math3.analysis.BivariateFunction;
+
+/**
+ *
+ * @author Alexander Nozik <altavir@gmail.com>
+ */
+public class NumassTransmission extends AbstractParametricBiFunction {
+
+    private static final String[] list = {"trap", "X"};
+    private static final double ION_POTENTIAL = 15.4;//eV
+    private final LossCalculator calculator;
+    private final BivariateFunction trapFunc;
+
+    public NumassTransmission() {
+        super(list);
+        this.calculator = LossCalculator.instance();
+        trapFunc = LossCalculator.getTrapFunction();
+    }
+
+    public NumassTransmission(BivariateFunction trapFunc) {
+        super(list);
+        this.calculator = LossCalculator.instance();
+        this.trapFunc = trapFunc;
+    }
+
+    @Override
+    public double derivValue(String parName, double eIn, double eOut, NamedValueSet set) {
+        switch (parName) {
+            case "trap":
+                return trapFunc.value(eIn, eOut);
+            case "X":
+                return calculator.getTotalLossDeriv(getX(eIn, set), eIn, eOut);
+            default:
+                return super.derivValue(parName, eIn, eOut, set);
+        }
+    }
+
+    @Override
+    public boolean providesDeriv(String name) {
+        return true;
+    }
+
+    private double getX(double eIn, NamedValueSet set) {
+        //From our article
+        return Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6;
+    }
+
+    @Override
+    public double value(double eIn, double eOut, NamedValueSet set) {
+        return calculator.getTotalLossValue(getX(eIn, set), eIn, eOut) + getParameter("trap", set) * trapFunc.value(eIn, eOut);
+    }
+
+}

numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java

@@ -6,6 +6,7 @@
 package inr.numass.models.sterile;
 
 import hep.dataforge.context.Context;
+import hep.dataforge.context.GlobalContext;
 import hep.dataforge.exceptions.NotDefinedException;
 import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.fitting.parametric.ParametricBiFunction;
@@ -30,30 +31,39 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     private double eMax;
     private final RandomGenerator rnd;
-    private final RealDistribution fssDistribution;
+    private RealDistribution fssDistribution;
 
     /**
      * variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2"
      */
-    private ParametricBiFunction source;
+    private final ParametricBiFunction source = new NumassBeta();
 
     /**
      * variables:Ein,Eout; parameters: "A"
      */
-    private ParametricBiFunction transmission;
+    private final ParametricBiFunction transmission;
     /**
      * variables:Eout,U; parameters: "X", "trap"
      */
-    private ParametricBiFunction resolution;
+    private final ParametricBiFunction resolution;
 
     public SterileNeutrinoSpectrum(Context context, Meta configuration) {
         super(list);
-        this.eMax = 18600;
-        FSS fss = new FSS(context.io().getFile(configuration.getString("fssFile", "FS.txt")));
         rnd = new SynchronizedRandomGenerator(new JDKRandomGenerator());
+        this.eMax = 18600;
+        if (configuration.hasValue("fssFile")) {
+            FSS fss = new FSS(context.io().getFile(configuration.getString("fssFile")));
             fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs());
         }
 
+        transmission = new NumassTransmission();
+        resolution = new NumassResolution(configuration.getNode("resolution", Meta.empty()));
+    }
+
+    public SterileNeutrinoSpectrum() {
+        this(GlobalContext.instance(), Meta.empty());
+    }
+
     @Override
     public double derivValue(String parName, double u, NamedValueSet set) {
         switch (parName) {
@@ -88,6 +98,16 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
         return rnd.nextDouble() * (b - a) + a;
     }
 
+    /**
+     * Monte-Carlo integration of spectrum
+     *
+     * @param u
+     * @param sourceFunction
+     * @param transmissionFunction
+     * @param resolutionFunction
+     * @param set
+     * @return
+     */
     private double integrate(
             double u,
             ParametricBiFunction sourceFunction,
@@ -97,16 +117,27 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
         int num = numCalls(u);
         double sum = DoubleStream.generate(() -> {
             // generate final state
-            double fs = fssDistribution.sample();
+            double fs;
+            if (fssDistribution != null) {
+                fs = fssDistribution.sample();
+            } else {
+                fs = 0;
+            }
 
             double eIn = rndE(u, eMax);
 
             double eOut = rndE(u, eIn);
 
-            return sourceFunction.value(fs, eIn, set)
+            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;
         }).limit(num).parallel().sum();
+        //triangle surface
         return Math.pow(eMax - u, 2d) / 2d * sum / num;
     }