[no commit message]

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

@@ -57,7 +57,7 @@ setDefault(Locale.US);
 //ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d);
 //beta.setCaching(false)
 
-Meta cfg = new GrindMetaBuilder().meta{
+Meta cfg = new GrindMetaBuilder().meta(){
     resolution(width: 8.3e-5)
 }.build();
 

numass-main/src/main/java/inr/numass/NumassPlugin.java

@@ -54,6 +54,7 @@ import inr.numass.models.RangedNamedSetSpectrum;
 import inr.numass.models.ResolutionFunction;
 import inr.numass.models.TransmissionInterpolator;
 import inr.numass.models.VariableLossSpectrum;
+import inr.numass.models.sterile.SterileNeutrinoSpectrum;
 import org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
 import org.apache.commons.math3.util.FastMath;
@@ -206,6 +207,13 @@ public class NumassPlugin extends BasicPlugin {
         });
 
         manager.addModel("sterile", (context, meta) -> {
+            SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(meta);
+            NBkgSpectrum spectrum = new NBkgSpectrum(sp);
+
+            return new XYModel(spectrum, getAdapter(meta));
+        });
+
+        manager.addModel("sterile-old", (context, 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);
@@ -227,7 +235,7 @@ public class NumassPlugin extends BasicPlugin {
             switch (meta.getString("trappingFunction", "default")) {
                 case "run2016":
                     sp.setTrappingFunction((Ei, Ef) -> {
-                        return 6.2e-5 * FastMath.exp(-(Ei - Ef) / 350d) + 1.97e-4 - 6.818e-9 * Ei;
+                        return 3.92e-5 * FastMath.exp(-(Ei - Ef) / 300d) + 1.97e-4 - 6.818e-9 * Ei;
                     });
                     break;
                 default:

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

@@ -12,10 +12,10 @@ import hep.dataforge.description.ValueDef;
 import hep.dataforge.exceptions.NotDefinedException;
 import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.fitting.parametric.ParametricBiFunction;
+import hep.dataforge.maths.integration.UnivariateIntegrator;
 import hep.dataforge.meta.Meta;
 import hep.dataforge.values.NamedValueSet;
 import inr.numass.NumassIntegrator;
-import inr.numass.NumassContext;
 import inr.numass.models.FSS;
 import java.util.stream.DoubleStream;
 import org.apache.commons.math3.analysis.BivariateFunction;
@@ -56,8 +56,9 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
      * variables:Eout,U; parameters: "X", "trap"
      */
     private final ParametricBiFunction resolution;
-    
+
     private boolean useMC;
+    private boolean fast;
 
     public SterileNeutrinoSpectrum(Context context, Meta configuration) {
         super(list);
@@ -69,13 +70,14 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
         transmission = new NumassTransmission(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);
     }
-    
+
     public SterileNeutrinoSpectrum(Meta configuration) {
-        this(GlobalContext.instance(),configuration);
+        this(GlobalContext.instance(), configuration);
     }
-    
+
     public SterileNeutrinoSpectrum() {
         this(GlobalContext.instance(), Meta.empty());
     }
@@ -120,9 +122,10 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     /**
      * Random E generator
+     *
      * @param a
      * @param b
-     * @return 
+     * @return
      */
     private double rndE(double a, double b) {
         return rnd.nextDouble() * (b - a) + a;
@@ -179,12 +182,13 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     /**
      * Direct Gauss-Legandre integration
+     *
      * @param u
      * @param sourceFunction
      * @param transmissionFunction
      * @param resolutionFunction
      * @param set
-     * @return 
+     * @return
      */
     private double integrateDirect(
             double u,
@@ -214,7 +218,20 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
         BivariateFunction transRes = (eIn, usp) -> {
             UnivariateFunction integrand = (eOut) -> transmissionFunction.value(eIn, eOut, set) * resolutionFunction.value(eOut, usp, set);
-            return NumassIntegrator.getDefaultIntegrator().integrate(integrand, usp, eIn);
+
+            double border = u + 30;
+            double firstPart = NumassIntegrator.getFastInterator().integrate(integrand, usp, Math.min(eIn, border));
+            double secondPart;
+            if (eIn > border) {
+                if(fast){
+                    secondPart = NumassIntegrator.getFastInterator().integrate(integrand, border, eIn);
+                } else {
+                    secondPart = NumassIntegrator.getDefaultIntegrator().integrate(integrand, border, eIn);
+                }
+            } else {
+                secondPart = 0;
+            }
+            return firstPart + secondPart;
         };
 
         UnivariateFunction integrand = (eIn) -> {
@@ -222,7 +239,14 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             return fsSource.value(eIn) * (p0 * resolutionFunction.value(eIn, u, set) + transRes.value(eIn, u));
         };
 
-        return NumassIntegrator.getDefaultIntegrator().integrate(integrand, u, eMax);
+        UnivariateIntegrator integrator;
+        if (fast && eMax - u < 500) {
+            integrator = NumassIntegrator.getFastInterator();
+        } else {
+            integrator = NumassIntegrator.getDefaultIntegrator();
+        }
+
+        return integrator.integrate(integrand, u, eMax);
     }
 
 }