From 21c1edb3234ebbe0b7533539eb1b638aa1991cf8 Mon Sep 17 00:00:00 2001 From: darksnake Date: Mon, 25 Jul 2016 17:02:21 +0300 Subject: [PATCH] [no commit message] --- .../inr/numass/scripts/FindExIonRatio.groovy | 2 +- .../groovy/inr/numass/scripts/Loss2014.groovy | 2 +- .../groovy/inr/numass/scripts/OldTest.groovy | 12 +- .../inr/numass/scripts/ResolutionTest.groovy | 10 +- .../numass/scripts/SignificanceTest.groovy | 2 +- .../groovy/inr/numass/scripts/Simulate.groovy | 12 +- .../inr/numass/scripts/SimulateGun.groovy | 6 +- .../inr/numass/scripts/SterileDemo.groovy | 19 +- .../numass/scripts/SystTransmission.groovy | 12 +- .../inr/numass/scripts/Systematics.groovy | 12 +- ...estExperimentalVariableLossSpectrum.groovy | 2 +- .../inr/numass/scripts/TritiumTest.groovy | 10 +- .../java/inr/numass/models/BetaSpectrum.java | 5 +- .../inr/numass/models/sterile/NumassBeta.java | 201 ++++++++++++++++++ .../models/sterile/NumassResolution.java | 78 +++++++ .../models/sterile/NumassTransmission.java | 63 ++++++ .../sterile/SterileNeutrinoSpectrum.java | 49 ++++- 17 files changed, 435 insertions(+), 62 deletions(-) create mode 100644 numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java create mode 100644 numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java create mode 100644 numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java diff --git a/numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy b/numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy index 2a92c27a..5bb4e2fe 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy +++ b/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 diff --git a/numass-main/src/main/groovy/inr/numass/scripts/Loss2014.groovy b/numass-main/src/main/groovy/inr/numass/scripts/Loss2014.groovy index 792f006c..1934f561 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Loss2014.groovy +++ b/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 diff --git a/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy b/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy index 61c50817..bce84c8d 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/OldTest.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.MINUITPlugin +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.MINUITPlugin -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; -import hep.dataforge.likelihood.BayesianManager +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; +import hep.dataforge.fitting.likelihood.BayesianManager import inr.numass.data.SpectrumDataAdapter import inr.numass.models.BetaSpectrum; import inr.numass.models.ModularSpectrum; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/ResolutionTest.groovy b/numass-main/src/main/groovy/inr/numass/scripts/ResolutionTest.groovy index ea7e9a00..5fc03816 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/ResolutionTest.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.FitTask; -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.FitTask; +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; import hep.dataforge.exceptions.NamingException; import inr.numass.data.SpectrumDataAdapter; import inr.numass.data.SpectrumGenerator; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SignificanceTest.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SignificanceTest.groovy index 38c67c85..904a2e87 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SignificanceTest.groovy +++ b/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; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy b/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy index 26eb9b86..eaa947ff 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.FitTask; -import hep.dataforge.datafitter.MINUITPlugin +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.FitTask; +import hep.dataforge.fitting.MINUITPlugin -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; import hep.dataforge.exceptions.NamingException; import hep.dataforge.exceptions.PackFormatException; import inr.numass.data.SpectrumDataAdapter; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy index e494b2a3..07965b93 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SimulateGun.groovy +++ b/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.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; import hep.dataforge.exceptions.NamingException; import hep.dataforge.io.FittingIOUtils; import inr.numass.data.SpectrumDataAdapter; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy index fe2464f6..40fcd573 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.FitTask; -import hep.dataforge.datafitter.MINUITPlugin +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.FitTask; +import hep.dataforge.fitting.MINUITPlugin -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.fitting.ParamSet; +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); -beta.setCaching(false); +//ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d); +//beta.setCaching(false); +ParametricFunction beta = new SterileNeutrinoSpectrum(); NBkgSpectrum spectrum = new NBkgSpectrum(beta); XYModel model = new XYModel(spectrum, new SpectrumDataAdapter()); diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy index 7818ef51..9877144e 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SystTransmission.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.FitTask; -import hep.dataforge.datafitter.MINUITPlugin +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.FitTask; +import hep.dataforge.fitting.MINUITPlugin -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; import hep.dataforge.exceptions.NamingException; import hep.dataforge.exceptions.PackFormatException; import inr.numass.data.SpectrumDataAdapter; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy b/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy index f7747bca..8c7edd99 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.FitTask; -import hep.dataforge.datafitter.MINUITPlugin +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.FitTask; +import hep.dataforge.fitting.MINUITPlugin -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; import hep.dataforge.exceptions.NamingException; import hep.dataforge.exceptions.PackFormatException; import inr.numass.data.SpectrumDataAdapter; diff --git a/numass-main/src/main/groovy/inr/numass/scripts/TestExperimentalVariableLossSpectrum.groovy b/numass-main/src/main/groovy/inr/numass/scripts/TestExperimentalVariableLossSpectrum.groovy index 72c91a5c..add23256 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/TestExperimentalVariableLossSpectrum.groovy +++ b/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 diff --git a/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy b/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy index 2d6314b9..da811cb6 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy +++ b/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.datafitter.FitState; -import hep.dataforge.datafitter.ParamSet; -import hep.dataforge.datafitter.models.XYModel; -import hep.dataforge.likelihood.BayesianManager +import hep.dataforge.fitting.FitManager; +import hep.dataforge.fitting.FitState; +import hep.dataforge.fitting.ParamSet; +import hep.dataforge.fitting.models.XYModel; +import hep.dataforge.fitting.likelihood.BayesianManager import static hep.dataforge.maths.RandomUtils.setSeed; import inr.numass.data.SpectrumGenerator; import inr.numass.models.BetaSpectrum diff --git a/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java b/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java index b60adc6a..dd325ad3 100644 --- a/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java +++ b/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); } diff --git a/numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java b/numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java new file mode 100644 index 00000000..2e1643f3 --- /dev/null +++ b/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 + */ +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); + } + +} 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 new file mode 100644 index 00000000..343bcfe0 --- /dev/null +++ b/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 + */ +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)); + } + } + +} 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 new file mode 100644 index 00000000..7cf095a5 --- /dev/null +++ b/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 + */ +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); + } + +} 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 f0e9b8ff..cd8739a4 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 @@ -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,28 +31,37 @@ 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()); - fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs()); + 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 @@ -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; }