[no commit message]

numass-main/build.gradle

@@ -14,6 +14,7 @@ dependencies {
     compile project(':dataforge-fx')
     compile project(':dataforge-plots')
     compile project(':numass-storage')
+    compile project(':dataforge-grind')
 }
 
 task workbench(dependsOn: classes, type: JavaExec){

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

@@ -42,6 +42,8 @@ import static java.util.Locale.setDefault;
 import inr.numass.utils.TritiumUtils;
 import inr.numass.data.SpectrumDataAdapter;
 import hep.dataforge.io.FittingIOUtils
+import hep.dataforge.meta.Meta
+import hep.dataforge.grind.GrindMetaBuilder
 
 /**
  *
@@ -50,9 +52,16 @@ import hep.dataforge.io.FittingIOUtils
 
 setDefault(Locale.US);
 
-ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d);
+//ParametricFunction beta = new BetaSpectrum();
 
-//ParametricFunction beta = new SterileNeutrinoSpectrum();
+//ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d);
+//beta.setCaching(false)
+
+Meta cfg = new GrindMetaBuilder().meta{
+    resolution(width: 8.3e-5)
+}.build();
+
+ParametricFunction beta = new SterileNeutrinoSpectrum(cfg);
 
 NBkgSpectrum spectrum = new NBkgSpectrum(beta);
 XYModel model = new XYModel(spectrum, new SpectrumDataAdapter());
@@ -65,10 +74,10 @@ allPars.setPar("E0", 18574.94, 1.4);
 allPars.setPar("mnu2", 0d, 1d);
 allPars.setPar("msterile2", 1000d * 1000d,0);
 allPars.setPar("U2", 0.0, 1e-4, -1d, 1d);
-allPars.setPar("X", 0.04000, 0.01, 0d, Double.POSITIVE_INFINITY);
+allPars.setPar("X", 0.04, 0.01, 0d, Double.POSITIVE_INFINITY);
 allPars.setPar("trap", 1.634, 0.01,0d, Double.POSITIVE_INFINITY);
 
-FittingIOUtils.printSpectrum(GlobalContext.out(), spectrum, allPars, 14000.0, 18600.0, 400);
+FittingIOUtils.printSpectrum(GlobalContext.out(), spectrum, allPars, 14000, 18600.0, 400);
 
 //SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316);
 //

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

@@ -34,10 +34,6 @@ import java.io.PrintWriter;
  */
 public class NumassContext extends Context {
 
-    public static UnivariateIntegrator defaultIntegrator = new GaussRuleIntegrator(300);
-    public static UnivariateIntegrator highDensityIntegrator = new GaussRuleIntegrator(500);
-    
-
     public NumassContext(Context parent, Meta config) {
         super(parent, "numass", config);
         init();

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

@@ -0,0 +1,42 @@
+/*
+ * 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.maths.integration.GaussRuleIntegrator;
+import hep.dataforge.maths.integration.UnivariateIntegrator;
+
+/**
+ *
+ * @author Alexander Nozik
+ */
+public class NumassIntegrator {
+
+    private static UnivariateIntegrator fastInterator;
+    private static UnivariateIntegrator defaultIntegrator;
+    private static UnivariateIntegrator highDensityIntegrator;
+
+    public static UnivariateIntegrator getFastInterator() {
+        if (fastInterator == null) {
+            fastInterator = new GaussRuleIntegrator(100);
+        }
+        return fastInterator;
+    }
+
+    public static UnivariateIntegrator getDefaultIntegrator() {
+        if (defaultIntegrator == null) {
+            defaultIntegrator = new GaussRuleIntegrator(300);
+        }
+        return defaultIntegrator;
+    }
+
+    public static UnivariateIntegrator getHighDensityIntegrator() {
+        if (highDensityIntegrator == null) {
+            highDensityIntegrator = new GaussRuleIntegrator(500);
+        }
+        return highDensityIntegrator;
+    }
+
+}

numass-main/src/main/java/inr/numass/actions/ShowLossSpectrumAction.java

@@ -40,6 +40,7 @@ import hep.dataforge.tables.MapPoint;
 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.models.ExperimentalVariableLossSpectrum;
 import inr.numass.models.LossCalculator;
@@ -79,7 +80,7 @@ public class ShowLossSpectrumAction extends OneToOneAction<FitState, FitState> {
                         new MetaBuilder("plot")
                         .setValue("plotTitle", "Differential scattering crossection for " + name)
                 );
-        switch (input.getModel().meta().getString("name","")) {
+        switch (input.getModel().meta().getString("name", "")) {
             case "scatter-variable":
                 scatterFunction = LossCalculator.getSingleScatterFunction(pars);
                 calculateRatio = true;
@@ -178,7 +179,7 @@ public class ShowLossSpectrumAction extends OneToOneAction<FitState, FitState> {
     }
 
     public static double calcultateIonRatio(NamedValueSet set, double threshold) {
-        UnivariateIntegrator integrator = NumassContext.highDensityIntegrator;
+        UnivariateIntegrator integrator = NumassIntegrator.getHighDensityIntegrator();
         UnivariateFunction integrand = LossCalculator.getSingleScatterFunction(set);
         return 1d - integrator.integrate(integrand, 5d, threshold);
     }

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

@@ -7,6 +7,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.utils.TritiumUtils;
 import org.apache.commons.math3.analysis.UnivariateFunction;
@@ -21,7 +22,7 @@ public class CustomNBkgSpectrum extends NBkgSpectrum {
     public static CustomNBkgSpectrum tritiumBkgSpectrum(ParametricFunction source, double amplitude){
         UnivariateFunction differentialBkgFunction = TritiumUtils.tritiumBackgroundFunction(amplitude);
         UnivariateFunction integralBkgFunction = 
-                (x) -> NumassContext.defaultIntegrator
+                (x) -> NumassIntegrator.getDefaultIntegrator()
                         .integrate(differentialBkgFunction, x, 18580d);
         return new CustomNBkgSpectrum(source, integralBkgFunction);
     }

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

@@ -19,14 +19,12 @@ import hep.dataforge.exceptions.NotDefinedException;
 import hep.dataforge.fitting.parametric.AbstractParametricFunction;
 import hep.dataforge.maths.integration.UnivariateIntegrator;
 import hep.dataforge.values.NamedValueSet;
-import inr.numass.NumassContext;
+import inr.numass.NumassIntegrator;
-import static java.lang.Math.abs;
 import static java.lang.Math.exp;
 import static java.lang.Math.sqrt;
 import org.apache.commons.math3.analysis.UnivariateFunction;
 import static java.lang.Math.abs;
 import static java.lang.Math.abs;
-import static java.lang.Math.abs;
 
 /**
  *
@@ -38,10 +36,9 @@ public class GunSpectrum extends AbstractParametricFunction {
     private final double cutoff = 4d;
     protected final UnivariateIntegrator integrator;
 
-
     public GunSpectrum() {
         super(list);
-        integrator = NumassContext.defaultIntegrator;
+        integrator = NumassIntegrator.getDefaultIntegrator();
     }
 
     @Override
@@ -50,7 +47,9 @@ public class GunSpectrum extends AbstractParametricFunction {
         final double sigma = set.getDouble("sigma");
         final double resA = set.getDouble("resA");
 
-        if(sigma == 0) throw new NotDefinedException();
+        if (sigma == 0) {
+            throw new NotDefinedException();
+        }
 
         UnivariateFunction integrand;
         switch (parName) {
@@ -140,7 +139,7 @@ public class GunSpectrum extends AbstractParametricFunction {
         final double sigma = set.getDouble("sigma");
         final double resA = set.getDouble("resA");
 
-        if (sigma <1e-5 ) {
+        if (sigma < 1e-5) {
             return transmissionValueFast(U, pos, resA);
         }
 

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

@@ -46,7 +46,7 @@ public class LossCalculator {
 
     private static final LossCalculator instance = new LossCalculator();
     private static final UnivariateIntegrator integrator = new GaussRuleIntegrator(100);
-//    private static final UnivariateIntegrator tailIntegrator = new GaussRuleIntegrator(50);
+
 
     public static UnivariateFunction getSingleScatterFunction() {
         final double A1 = 0.204;
@@ -282,6 +282,7 @@ public class LossCalculator {
      * @return
      */
     public List<Double> getLossProbabilities(double X) {
+
         List<Double> res = new ArrayList<>();
         double prob;
         if (X > 0) {

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

@@ -15,6 +15,7 @@
  */
 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;
@@ -38,7 +39,7 @@ class LossResConvolution implements BivariateFunction {
     public double value(final double Ein, final double U) {
         UnivariateFunction integrand = (double Eout) -> loss.value(Ein, Eout) * res.value(Eout, U);
         //Энергия в принципе не может быть больше начальной и меньше напряжения
-        return NumassContext.defaultIntegrator.integrate(integrand, U, Ein);
+        return NumassIntegrator.getDefaultIntegrator().integrate(integrand, U, Ein);
 
     }
 }

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

@@ -18,6 +18,7 @@ package inr.numass.models;
 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 org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
@@ -58,7 +59,7 @@ class TransmissionConvolution extends AbstractParametricFunction {
             }
             return trans.value(E, U) * spectrum.derivValue(parName, E, set);
         };
-        return NumassContext.defaultIntegrator.integrate(integrand, Math.max(U, min), max + 1d);
+        return NumassIntegrator.getDefaultIntegrator().integrate(integrand, Math.max(U, min), max + 1d);
     }
 
     @Override
@@ -80,6 +81,6 @@ class TransmissionConvolution extends AbstractParametricFunction {
             }
             return trans.value(E, U) * spectrum.value(E, set);
         };
-        return NumassContext.defaultIntegrator.integrate(integrand, Math.max(U, min), max + 1d);
+        return NumassIntegrator.getDefaultIntegrator().integrate(integrand, Math.max(U, min), max + 1d);
     }
 }

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

@@ -21,6 +21,7 @@ import hep.dataforge.fitting.parametric.ParametricFunction;
 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 java.util.List;
 import org.apache.commons.math3.analysis.BivariateFunction;
@@ -101,9 +102,9 @@ public class VariableLossSpectrum extends AbstractParametricFunction {
         UnivariateFunction integrand = (double x) -> transmission.value(x, set) * lossFunction.value(x, U - shift);
         UnivariateIntegrator integrator;
         if (eMax - U > 150) {
-            integrator = NumassContext.highDensityIntegrator;
+            integrator = NumassIntegrator.getHighDensityIntegrator();
         } else {
-            integrator = NumassContext.defaultIntegrator;
+            integrator = NumassIntegrator.getDefaultIntegrator();
         }
         return noLossProb * transmission.value(U - shift, set) + integrator.integrate(integrand, U, eMax);
     }

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

@@ -13,7 +13,12 @@ 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.
+ * A bi-function for beta-spectrum calculation taking energy and final state as
+ * input.
+ * <p>
+ * dissertation p.33
+ * </p>
+ *
  * @author Alexander Nozik <altavir@gmail.com>
  */
 public class NumassBeta extends AbstractParametricBiFunction {
@@ -25,6 +30,16 @@ public class NumassBeta extends AbstractParametricBiFunction {
         super(list);
     }
 
+    /**
+     * Beta spectrum derivative
+     *
+     * @param n parameter number
+     * @param E0
+     * @param mnu2
+     * @param E
+     * @return
+     * @throws NotDefinedException
+     */
     private double derivRoot(int n, double E0, double mnu2, double E) throws NotDefinedException {
         double D = E0 - E;//E0-E
         double res;
@@ -71,6 +86,7 @@ public class NumassBeta extends AbstractParametricBiFunction {
 
     /**
      * Derivative of spectrum with sterile neutrinos
+     *
      * @param name
      * @param E
      * @param E0
@@ -104,11 +120,14 @@ public class NumassBeta extends AbstractParametricBiFunction {
 
     }
 
+    /**
+     * The part independent of neutrino mass. Includes global normalization
+     * constant, momentum and Fermi correction
+     *
+     * @param E
+     * @return
+     */
     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));
@@ -118,7 +137,7 @@ public class NumassBeta extends AbstractParametricBiFunction {
         double Fn = y / abs(1d - exp(-y));
         double Fermi = Fn * (1.002037 - 0.001427 * ve);
         double res = Fermi * pe * Etot;
-        return res;
+        return K * res;
     }
 
     @Override
@@ -127,7 +146,8 @@ public class NumassBeta extends AbstractParametricBiFunction {
     }
 
     /**
-     * Bare beta spectrum with Mintz negative mass correction
+     * Bare beta spectrum with Mainz negative mass correction
+     *
      * @param E0
      * @param mnu2
      * @param E
@@ -155,6 +175,7 @@ public class NumassBeta extends AbstractParametricBiFunction {
 
     /**
      * beta-spectrum with sterile neutrinos
+     *
      * @param E
      * @param E0
      * @param pars
@@ -189,7 +210,7 @@ public class NumassBeta extends AbstractParametricBiFunction {
     @Override
     public double derivValue(String parName, double fs, double eIn, NamedValueSet pars) {
         double E0 = getParameter("E0", pars);
-        return derivRootsterile(parName,eIn, E0 - fs, pars);
+        return derivRootsterile(parName, eIn, E0 - fs, pars);
     }
 
     @Override

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

@@ -6,6 +6,7 @@
 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.LossCalculator;
 import org.apache.commons.math3.analysis.BivariateFunction;
@@ -21,16 +22,21 @@ public class NumassTransmission extends AbstractParametricBiFunction {
     private final LossCalculator calculator;
     private final BivariateFunction trapFunc;
 
-    public NumassTransmission() {
+//    private BicubicInterpolatingFunction cache;
+//    private double cachedX;
+//    private Meta cacheMeta;
+
+    public NumassTransmission(Meta meta) {
         super(list);
         this.calculator = LossCalculator.instance();
-        trapFunc = LossCalculator.getTrapFunction();
+        trapFunc = getTrapFunction(meta.getNodeOrEmpty("trapping"));
+//        if (meta.hasNode("cache")) {
+//            cacheMeta = meta.getNode("cache");
+//        }
     }
 
-    public NumassTransmission(BivariateFunction trapFunc) {
+    private BivariateFunction getTrapFunction(Meta meta) {
-        super(list);
+        return LossCalculator.getTrapFunction();
-        this.calculator = LossCalculator.instance();
-        this.trapFunc = trapFunc;
     }
 
     @Override
@@ -52,16 +58,34 @@ public class NumassTransmission extends AbstractParametricBiFunction {
 
     public static double getX(double eIn, NamedValueSet set) {
         //From our article
-        return set.getDouble("X")*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));
     }
 
+//    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) {
-        return calculator.getTotalLossValue(getX(eIn, set), eIn, eOut) + getParameter("trap", set) * trapFunc.value(eIn, eOut);
+        //calculate X taking into account its energy dependence
+        double X = getX(eIn, set);
+        // loss part
+        double loss = calculator.getTotalLossValue(X, eIn, eOut);
+        //trapping part
+        double trap = getParameter("trap", set) * trapFunc.value(eIn, eOut);
+        return loss + trap;
     }
 
 }

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

@@ -7,11 +7,14 @@ package inr.numass.models.sterile;
 
 import hep.dataforge.context.Context;
 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.AbstractParametricFunction;
 import hep.dataforge.fitting.parametric.ParametricBiFunction;
 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;
@@ -28,6 +31,10 @@ import org.apache.commons.math3.random.SynchronizedRandomGenerator;
  *
  * @author Alexander Nozik
  */
+@NodeDef(name = "resolution")
+@NodeDef(name = "transmission")
+@ValueDef(name = "fssFile")
+
 public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
     private static final String[] list = {"X", "trap", "E0", "mnu2", "msterile2", "U2"};
@@ -50,6 +57,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
      */
     private final ParametricBiFunction resolution;
     
+    private boolean useMC;
+
     public SterileNeutrinoSpectrum(Context context, Meta configuration) {
         super(list);
         rnd = new SynchronizedRandomGenerator(new JDKRandomGenerator());
@@ -58,8 +67,13 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs());
         }
 
-        transmission = new NumassTransmission();
+        transmission = new NumassTransmission(configuration.getNodeOrEmpty("transmission"));
         resolution = new NumassResolution(configuration.getNode("resolution", Meta.empty()));
+        this.useMC = configuration.getBoolean("useMC",false);
+    }
+    
+    public SterileNeutrinoSpectrum(Meta configuration) {
+        this(GlobalContext.instance(),configuration);
     }
     
     public SterileNeutrinoSpectrum() {
@@ -96,7 +110,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
     }
 
     private boolean useDirect() {
-        return true;
+        return !useMC;
     }
 
     @Override
@@ -104,6 +118,12 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
         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;
     }
@@ -157,6 +177,15 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
         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,
@@ -175,7 +204,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             fsSource = (eIn) -> {
                 double res = 0;
                 for (int i = 0; i < fss.size(); i++) {
-                    res += fss.getP(i) * sourceFunction.value(fss.getE(i), u, set);
+                    res += fss.getP(i) * sourceFunction.value(fss.getE(i), eIn, set);
                 }
                 return res;
             };
@@ -185,7 +214,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
 
         BivariateFunction transRes = (eIn, usp) -> {
             UnivariateFunction integrand = (eOut) -> transmissionFunction.value(eIn, eOut, set) * resolutionFunction.value(eOut, usp, set);
-            return NumassContext.defaultIntegrator.integrate(integrand, usp, eIn);
+            return NumassIntegrator.getDefaultIntegrator().integrate(integrand, usp, eIn);
         };
 
         UnivariateFunction integrand = (eIn) -> {
@@ -193,7 +222,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
             return fsSource.value(eIn) * (p0 * resolutionFunction.value(eIn, u, set) + transRes.value(eIn, u));
         };
 
-        return NumassContext.defaultIntegrator.integrate(integrand, u, eMax);
+        return NumassIntegrator.getDefaultIntegrator().integrate(integrand, u, eMax);
     }
 
 }