[no commit message]

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

@@ -50,9 +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);
 
-ParametricFunction beta = new SterileNeutrinoSpectrum();
+//ParametricFunction beta = new SterileNeutrinoSpectrum();
 
 NBkgSpectrum spectrum = new NBkgSpectrum(beta);
 XYModel model = new XYModel(spectrum, new SpectrumDataAdapter());

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

@@ -168,12 +168,12 @@ public class LossCalculator {
             final LossCalculator loss = LossCalculator.instance;
             final List<Double> probs = loss.getGunLossProbabilities(set.getDouble("X"));
             UnivariateFunction single = (double e) -> probs.get(1) * scatterFunction.value(e);
-            frame.add(PlottableXYFunction.plotFunction("Single scattering",  x->single.value(x), 0, 100, 1000));
+            frame.add(PlottableXYFunction.plotFunction("Single scattering", x -> single.value(x), 0, 100, 1000));
 
             for (int i = 2; i < probs.size(); i++) {
                 final int j = i;
                 UnivariateFunction scatter = (double e) -> probs.get(j) * loss.getLossValue(j, e, 0d);
-                frame.add(PlottableXYFunction.plotFunction(j + " scattering", x->scatter.value(x), 0, 100, 1000));
+                frame.add(PlottableXYFunction.plotFunction(j + " scattering", x -> scatter.value(x), 0, 100, 1000));
             }
 
             UnivariateFunction total = (eps) -> {
@@ -187,11 +187,11 @@ public class LossCalculator {
                 return sum;
             };
 
-            frame.add(PlottableXYFunction.plotFunction("Total loss", x->total.value(x), 0, 100, 1000));
+            frame.add(PlottableXYFunction.plotFunction("Total loss", x -> total.value(x), 0, 100, 1000));
 
         } else {
 
-            frame.add(PlottableXYFunction.plotFunction("Differential crosssection", x->scatterFunction.value(x), 0, 100, 2000));
+            frame.add(PlottableXYFunction.plotFunction("Differential crosssection", x -> scatterFunction.value(x), 0, 100, 2000));
         }
 
     }
@@ -419,6 +419,9 @@ public class LossCalculator {
      * @return
      */
     public double getTotalLossValue(double X, double Ei, double Ef) {
+        if (X == 0) {
+            return 0;
+        }
         List<Double> probs = getLossProbabilities(X);
 
         double sum = 0;

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

@@ -50,9 +50,13 @@ public class NumassTransmission extends AbstractParametricBiFunction {
         return true;
     }
 
-    private double getX(double eIn, NamedValueSet set) {
+    public static double getX(double eIn, NamedValueSet set) {
         //From our article
-        return Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6;
+        return set.getDouble("X")*Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6;
+    }
+    
+    public static double p0(double eIn, NamedValueSet set) {
+        return LossCalculator.instance().getLossProbability(0, getX(eIn, set));
     }
 
     @Override

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

@@ -12,8 +12,11 @@ import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.fitting.parametric.ParametricBiFunction;
 import hep.dataforge.meta.Meta;
 import hep.dataforge.values.NamedValueSet;
+import inr.numass.NumassContext;
 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;
@@ -29,9 +32,9 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     private static final String[] list = {"X", "trap", "E0", "mnu2", "msterile2", "U2"};
 
-    private double eMax;
     private final RandomGenerator rnd;
     private RealDistribution fssDistribution;
+    private FSS fss;
 
     /**
      * variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2"
@@ -50,9 +53,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
     public SterileNeutrinoSpectrum(Context context, Meta configuration) {
         super(list);
         rnd = new SynchronizedRandomGenerator(new JDKRandomGenerator());
-        this.eMax = 18600;
         if (configuration.hasValue("fssFile")) {
-            FSS fss = new FSS(context.io().getFile(configuration.getString("fssFile")));
+            fss = new FSS(context.io().getFile(configuration.getString("fssFile")));
             fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs());
         }
 
@@ -82,13 +84,21 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     @Override
     public double value(double u, NamedValueSet set) {
-        return integrate(u, source, transmission, resolution, 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 true;
+    }
+
     @Override
     public boolean providesDeriv(String name) {
         return source.providesDeriv(name) && transmission.providesDeriv(name) && resolution.providesDeriv(name);
@@ -114,7 +124,13 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             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;
@@ -131,8 +147,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             double res = sourceFunction.value(fs, eIn, set)
                     * transmissionFunction.value(eIn, eOut, set)
                     * resolutionFunction.value(eOut, u, set);
-            
+
-            if(Double.isNaN(res)){
+            if (Double.isNaN(res)) {
                 throw new Error();
             }
             return res;
@@ -141,4 +157,43 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
         return Math.pow(eMax - u, 2d) / 2d * sum / num;
     }
 
+    private double integrateDirect(
+            double u,
+            ParametricBiFunction sourceFunction,
+            ParametricBiFunction transmissionFunction,
+            ParametricBiFunction resolutionFunction,
+            NamedValueSet set) {
+
+        double eMax = set.getDouble("E0") + 5d;
+
+        if (u > eMax) {
+            return 0;
+        }
+
+        UnivariateFunction fsSource;
+        if (fss != null) {
+            fsSource = (eIn) -> {
+                double res = 0;
+                for (int i = 0; i < fss.size(); i++) {
+                    res += fss.getP(i) * sourceFunction.value(fss.getE(i), u, set);
+                }
+                return res;
+            };
+        } else {
+            fsSource = (eIn) -> sourceFunction.value(0, eIn, set);
+        }
+
+        BivariateFunction transRes = (eIn, usp) -> {
+            UnivariateFunction integrand = (eOut) -> transmissionFunction.value(eIn, eOut, set) * resolutionFunction.value(eOut, usp, set);
+            return NumassContext.defaultIntegrator.integrate(integrand, usp, eIn);
+        };
+
+        UnivariateFunction integrand = (eIn) -> {
+            double p0 = NumassTransmission.p0(eIn, set);
+            return fsSource.value(eIn) * (p0 * resolutionFunction.value(eIn, u, set) + transRes.value(eIn, u));
+        };
+
+        return NumassContext.defaultIntegrator.integrate(integrand, u, eMax);
+    }
+
 }