From b6af0f612cd920f8cfbb4273eee0194ec7359301 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 26 Nov 2017 22:15:16 +0300 Subject: [PATCH] Updating models --- build.gradle | 2 +- .../inr/numass/control/readvac/ConsoleVac.kt | 13 +- .../numass/scripts/ShowTransmission.groovy | 2 +- .../groovy/inr/numass/scripts/Simulate.groovy | 12 +- .../inr/numass/scripts/SterileDemo.groovy | 10 +- .../numass/scripts/models/TristanModel.groovy | 69 ++++++ .../scripts/underflow/UnderflowFitter.groovy | 4 +- .../inr/numass/models/GaussResolution.java | 199 ++++++++------- .../inr/numass/models/sterile/NumassBeta.java | 219 ----------------- .../models/sterile/NumassResolution.java | 99 -------- .../models/sterile/NumassTransmission.java | 114 --------- .../sterile/SterileNeutrinoSpectrum.java | 222 ----------------- .../inr/numass/models/sterile/TestModels.java | 99 -------- .../java/inr/numass/utils/DataModelUtils.java | 4 +- .../inr/numass/utils/NumassIntegrator.java | 7 +- .../src/main/kotlin/inr/numass/NumassUtils.kt | 24 +- .../kotlin/inr/numass/models/NumassModels.kt | 114 +++++++++ .../inr/numass/models/sterile/NumassBeta.kt | 230 ++++++++++++++++++ .../numass/models/sterile/NumassResolution.kt | 91 +++++++ .../models/sterile/NumassTransmission.kt | 105 ++++++++ .../models/sterile/SterileNeutrinoSpectrum.kt | 164 +++++++++++++ .../inr/numass/models/sterile/TestModels.kt | 104 ++++++++ .../kotlin/inr/numass/tasks/NumassTasks.kt | 4 +- 23 files changed, 1024 insertions(+), 887 deletions(-) create mode 100644 numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy delete mode 100644 numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java delete mode 100644 numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java delete mode 100644 numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java delete mode 100644 numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java delete mode 100644 numass-main/src/main/java/inr/numass/models/sterile/TestModels.java create mode 100644 numass-main/src/main/kotlin/inr/numass/models/sterile/NumassBeta.kt create mode 100644 numass-main/src/main/kotlin/inr/numass/models/sterile/NumassResolution.kt create mode 100644 numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt create mode 100644 numass-main/src/main/kotlin/inr/numass/models/sterile/SterileNeutrinoSpectrum.kt create mode 100644 numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt diff --git a/build.gradle b/build.gradle index bd715132..481c5884 100644 --- a/build.gradle +++ b/build.gradle @@ -1,5 +1,5 @@ plugins { - id "org.jetbrains.kotlin.jvm" version "1.1.60" apply false + id "org.jetbrains.kotlin.jvm" version "1.1.61" apply false } allprojects{ diff --git a/numass-control/vac/src/main/kotlin/inr/numass/control/readvac/ConsoleVac.kt b/numass-control/vac/src/main/kotlin/inr/numass/control/readvac/ConsoleVac.kt index df1e3522..fb2522f0 100644 --- a/numass-control/vac/src/main/kotlin/inr/numass/control/readvac/ConsoleVac.kt +++ b/numass-control/vac/src/main/kotlin/inr/numass/control/readvac/ConsoleVac.kt @@ -22,7 +22,7 @@ object ConsoleVac { options.addOption("p", "port", true, "Port name in dataforge-control notation") options.addOption("d", "delay", true, "A delay between measurements in Duration notation") - if (args.size == 0) { + if (args.isEmpty()) { HelpFormatter().printHelp("vac-console", options) return } @@ -30,19 +30,18 @@ object ConsoleVac { val parser = DefaultParser() val cli = parser.parse(options, args) - var className: String? = cli.getOptionValue("c") - if (!className!!.contains(".")) { + var className: String = cli.getOptionValue("c") ?: throw RuntimeException("Vacuumeter class not defined") + + if (!className.contains(".")) { className = "inr.numass.readvac.devices." + className } + val name = cli.getOptionValue("n", "sensor") val port = cli.getOptionValue("p", "com::/dev/ttyUSB0") val delay = Duration.parse(cli.getOptionValue("d", "PT1M")) - if (className == null) { - throw RuntimeException("Vacuumeter class not defined") - } val sensor = Class.forName(className) - .getConstructor(String::class.java).newInstance(port) as Sensor + .getConstructor(String::class.java).newInstance(port) as Sensor<*> try { sensor.init() while (true) { diff --git a/numass-main/src/main/groovy/inr/numass/scripts/ShowTransmission.groovy b/numass-main/src/main/groovy/inr/numass/scripts/ShowTransmission.groovy index 84373104..562e9f41 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/ShowTransmission.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/ShowTransmission.groovy @@ -22,7 +22,7 @@ shell.eval { PlotHelper plot = plots - ParamSet params = morph(ParamSet, "params") { + ParamSet params = morph(ParamSet,[:], "params") { N(value: 2.7e+06, err: 30, lower: 0) bkg(value: 5.0, err: 0.1) E0(value: 18575.0, err: 0.1) 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 26eec6b6..1005dbfd 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy @@ -47,11 +47,11 @@ SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(Global.instance(), Meta //beta.setCaching(false); NBkgSpectrum spectrum = new NBkgSpectrum(sp); -XYModel model = new XYModel(spectrum, new SpectrumDataAdapter()); +XYModel model = new XYModel(Meta.empty(), new SpectrumDataAdapter(), spectrum); ParamSet allPars = new ParamSet(); -allPars.setParValue("N", 2e6/100); +allPars.setParValue("N", 2e6 / 100); //значение 6е-6 соответствует полной интенстивности 6е7 распадов в секунду //Проблема была в переполнении счетчика событий в генераторе. Заменил на long. Возможно стоит поставить туда число с плавающей точкой allPars.setParError("N", 6); @@ -86,14 +86,14 @@ def data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration // allPars.setParValue("X", 0.4); -ColumnedDataWriter.writeTable(System.out,data,"--- DATA ---"); +ColumnedDataWriter.writeTable(System.out, data, "--- DATA ---"); FitState state = new FitState(data, model, allPars); //new PlotFitResultAction().eval(state); - - + + def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2"); - + res.print(out()); 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 7824ef43..99d57b5b 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy @@ -40,25 +40,25 @@ 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(); -ParametricFunction beta = new SterileNeutrinoSpectrum(cfg); +ParametricFunction beta = new SterileNeutrinoSpectrum(Global.instance(), cfg); NBkgSpectrum spectrum = new NBkgSpectrum(beta); XYModel model = new XYModel(spectrum, new SpectrumDataAdapter()); ParamSet allPars = new ParamSet(); -allPars.setPar("N", 6.6579e+05, 1.8e+03, 0d, Double.POSITIVE_INFINITY); +allPars.setPar("N", 6.6579e+05, 1.8e+03, 0d, Double.POSITIVE_INFINITY); allPars.setPar("bkg", 0.5387, 0.050); allPars.setPar("E0", 18574.94, 1.4); allPars.setPar("mnu2", 0d, 1d); -allPars.setPar("msterile2", 1000d * 1000d,0); +allPars.setPar("msterile2", 1000d * 1000d, 0); allPars.setPar("U2", 0.0, 1e-4, -1d, 1d); allPars.setPar("X", 0.04, 0.01, 0d, Double.POSITIVE_INFINITY); -allPars.setPar("trap", 1.634, 0.01,0d, Double.POSITIVE_INFINITY); +allPars.setPar("trap", 1.634, 0.01, 0d, Double.POSITIVE_INFINITY); FittingIOUtils.printSpectrum(Global.out(), spectrum, allPars, 14000, 18600.0, 400); diff --git a/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy b/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy new file mode 100644 index 00000000..941a3f53 --- /dev/null +++ b/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy @@ -0,0 +1,69 @@ +package inr.numass.scripts.models + +import hep.dataforge.context.Context +import hep.dataforge.context.Global +import hep.dataforge.fx.plots.PlotManager +import hep.dataforge.grind.Grind +import hep.dataforge.grind.GrindShell +import hep.dataforge.grind.helpers.PlotHelper +import hep.dataforge.meta.Meta +import hep.dataforge.stat.fit.ParamSet +import hep.dataforge.stat.models.XYModel +import hep.dataforge.stat.parametric.ParametricFunction +import inr.numass.NumassPlugin +import inr.numass.data.SpectrumGenerator +import inr.numass.models.GaussResolution +import inr.numass.models.NBkgSpectrum +import inr.numass.models.sterile.NumassBeta +import inr.numass.models.sterile.NumassTransmission +import inr.numass.models.sterile.SterileNeutrinoSpectrum + +import static hep.dataforge.grind.Grind.morph + +Context ctx = Global.instance() +ctx.getPluginManager().load(PlotManager) +ctx.getPluginManager().load(NumassPlugin) + +new GrindShell(ctx).eval { + ParametricFunction spectrum = new SterileNeutrinoSpectrum( + ctx, + Grind.buildMeta(useFSS: false), + new NumassBeta(), + new NumassTransmission(ctx, Meta.empty()), + new GaussResolution(5d) + ) + + def model = new XYModel(Meta.empty(), new NBkgSpectrum(spectrum)); + + ParamSet params = morph(ParamSet, [:], "params") { + N(value: 1e+04, err: 30, lower: 0) + bkg(value: 1.0, err: 0.1) + E0(value: 18575.0, err: 0.1) + mnu2(value: 0, err: 0.01) + msterile2(value: 1000**2, err: 1) + U2(value: 0.0, err: 1e-3) + X(value: 0.0, err: 0.01, lower: 0) + trap(value: 1.0, err: 0.05) + w(300, err: 5) + } + + SpectrumGenerator generator = new SpectrumGenerator(model, params, 12316); + + PlotHelper ph = plots + + ph.plot((10000..18500).step(100).collectEntries { + [it,model.value(it,params)] + }) + +// +// def data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(10000, 18500, 10*24*3600, 850)); +// +// FitState state = new FitState(data, model, params); +// +// def fm = ctx.getFeature(FitManager) +// +// def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2"); +// +// +// res.print(ctx.io.out()); +} \ No newline at end of file diff --git a/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowFitter.groovy b/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowFitter.groovy index e03374be..5843d699 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowFitter.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowFitter.groovy @@ -66,8 +66,8 @@ class UnderflowFitter { .map { new WeightedObservedPoint( 1d,//1d / p.getValue() , //weight - it.getDouble(CHANNEL_KEY), // x - it.getDouble(COUNT_RATE_KEY) / binning) //y + it.getDouble(CHANNEL_KEY) as double, // x + it.getDouble(COUNT_RATE_KEY) / binning as double) //y } .collect(Collectors.toList()); SimpleCurveFitter fitter = SimpleCurveFitter.create(new ExponentFunction(), [1d, 200d] as double[]) diff --git a/numass-main/src/main/java/inr/numass/models/GaussResolution.java b/numass-main/src/main/java/inr/numass/models/GaussResolution.java index d772c9fd..1fc3616b 100644 --- a/numass-main/src/main/java/inr/numass/models/GaussResolution.java +++ b/numass-main/src/main/java/inr/numass/models/GaussResolution.java @@ -15,135 +15,107 @@ */ package inr.numass.models; -import hep.dataforge.context.Global; -import hep.dataforge.exceptions.NameNotFoundException; -import hep.dataforge.stat.parametric.AbstractParametricFunction; -import hep.dataforge.stat.parametric.ParametricFunction; +import hep.dataforge.stat.parametric.AbstractParametricBiFunction; import hep.dataforge.values.ValueProvider; import hep.dataforge.values.Values; -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.analysis.integration.SimpsonIntegrator; -import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; -import static hep.dataforge.names.NamesUtils.combineNamesWithEquals; -import static java.lang.Double.isNaN; import static java.lang.Math.*; /** - * * @author Darksnake */ -public class GaussResolution extends AbstractParametricFunction implements Transmission { +public class GaussResolution extends AbstractParametricBiFunction { - private static final String[] list = {"w"}; + private static final String[] list = {"w", "shift"}; private double cutoff = 4d; - private final UnivariateIntegrator integrator = new SimpsonIntegrator(1e-3, Double.MAX_VALUE, 3, 64); +// private final UnivariateIntegrator integrator = new SimpsonIntegrator(1e-3, Double.MAX_VALUE, 3, 64); /** - * * @param cutoff - расстояние в сигмах от среднего значения, на котором - * функция считается обращающейся в ноль + * функция считается обращающейся в ноль */ public GaussResolution(double cutoff) { super(list); this.cutoff = cutoff; } - @Override - public double derivValue(String name, double X, Values pars) { - if (abs(X - getPos(pars)) > cutoff * getW(pars)) { - return 0; - } - switch (name) { - case "pos": - return this.value(X, pars) * (X - getPos(pars)) / getW(pars) / getW(pars); - case "w": - return this.value(X, pars) * ((X - getPos(pars)) * (X - getPos(pars)) / getW(pars) / getW(pars) / getW(pars) - 1 / getW(pars)); - default: - return 0; - } - } +// @Override +// public ParametricFunction getConvolutedSpectrum(final RangedNamedSetSpectrum bare) { +// return new AbstractParametricFunction(combineNamesWithEquals(this.namesAsArray(), bare.namesAsArray())) { +// int maxEval = Global.instance().getInt("INTEGR_POINTS", 500); +// +// @Override +// public double derivValue(String parName, double x, Values set) { +// double a = getLowerBound(set); +// double b = getUpperBound(set); +// assert b > a; +// return integrator.integrate(maxEval, getDerivProduct(parName, bare, set, x), a, b); +// } +// +// @Override +// public boolean providesDeriv(String name) { +// return "w".equals(name) || bare.providesDeriv(name); +// } +// +// @Override +// public double value(double x, Values set) { +// double a = getLowerBound(set); +// double b = getUpperBound(set); +// assert b > a; +// return integrator.integrate(maxEval, getProduct(bare, set, x), a, b); +// } +// }; +// } - @Override - public ParametricFunction getConvolutedSpectrum(final RangedNamedSetSpectrum bare) { - return new AbstractParametricFunction(combineNamesWithEquals(this.namesAsArray(), bare.namesAsArray())) { - int maxEval = Global.instance().getInt("INTEGR_POINTS", 500); +// @Override +// public double getDeriv(String name, Values set, double input, double output) { +// return this.derivValue(name, output - input, set); +// } - @Override - public double derivValue(String parName, double x, Values set) { - double a = getLowerBound(set); - double b = getUpperBound(set); - assert b > a; - return integrator.integrate(maxEval, getDerivProduct(parName, bare, set, x), a, b); - } +// private UnivariateFunction getDerivProduct(final String name, final ParametricFunction bare, final Values pars, final double x0) { +// return (double x) -> { +// double res1; +// double res2; +// try { +// res1 = bare.derivValue(name, x0 - x, pars) * GaussResolution.this.value(x, pars); +// } catch (NameNotFoundException ex1) { +// res1 = 0; +// } +// try { +// res2 = bare.value(x0 - x, pars) * GaussResolution.this.derivValue(name, x, pars); +// } catch (NameNotFoundException ex2) { +// res2 = 0; +// } +// return res1 + res2; +// }; +// } - @Override - public boolean providesDeriv(String name) { - if ("w".equals(name)) { - return true; - } - return bare.providesDeriv(name); - } +// private double getLowerBound(final Values pars) { +// return getPos(pars) - cutoff * getW(pars); +// } - @Override - public double value(double x, Values set) { - double a = getLowerBound(set); - double b = getUpperBound(set); - assert b > a; - return integrator.integrate(maxEval, getProduct(bare, set, x), a, b); - } - }; - } - @Override - public double getDeriv(String name, Values set, double input, double output) { - return this.derivValue(name, output - input, set); - } - private UnivariateFunction getDerivProduct(final String name, final ParametricFunction bare, final Values pars, final double x0) { - return (double x) -> { - double res1; - double res2; - try { - res1 = bare.derivValue(name, x0 - x, pars) * GaussResolution.this.value(x, pars); - } catch (NameNotFoundException ex1) { - res1 = 0; - } - try { - res2 = bare.value(x0 - x, pars) * GaussResolution.this.derivValue(name, x, pars); - } catch (NameNotFoundException ex2) { - res2 = 0; - } - return res1 + res2; - }; - } +// private UnivariateFunction getProduct(final ParametricFunction bare, final Values pars, final double x0) { +// return (double x) -> { +// double res = bare.value(x0 - x, pars) * GaussResolution.this.value(x, pars); +// assert !isNaN(res); +// return res; +// }; +// } - private double getLowerBound(final Values pars) { - return getPos(pars) - cutoff * getW(pars); - } +// private double getUpperBound(final Values pars) { +// return getPos(pars) + cutoff * getW(pars); +// } + +// @Override +// public double getValue(Values set, double input, double output) { +// return this.value(output - input, set); +// } private double getPos(ValueProvider pars) { -// return pars.getDouble("pos"); - // вряд ли стоит ожидать, что разрешение будет сдвигать среднее, поэтому оставляем один параметр - return 0; - } - - private UnivariateFunction getProduct(final ParametricFunction bare, final Values pars, final double x0) { - return (double x) -> { - double res = bare.value(x0 - x, pars) * GaussResolution.this.value(x, pars); - assert !isNaN(res); - return res; - }; - } - - private double getUpperBound(final Values pars) { - return getPos(pars) + cutoff * getW(pars); - } - - @Override - public double getValue(Values set, double input, double output) { - return this.value(output - input, set); + return pars.getDouble("shift", 0); } private double getW(ValueProvider pars) { @@ -155,12 +127,33 @@ public class GaussResolution extends AbstractParametricFunction implements Trans return true; } + @Override - public double value(double x, Values pars) { - if (abs(x - getPos(pars)) > cutoff * getW(pars)) { + public double value(double x, double y, Values pars) { + double d = x - y; + if (abs(d - getPos(pars)) > cutoff * getW(pars)) { return 0; } - double aux = (x - getPos(pars)) / getW(pars); + double aux = (d - getPos(pars)) / getW(pars); return exp(-aux * aux / 2) / getW(pars) / sqrt(2 * Math.PI); } + + @Override + public double derivValue(String parName, double x, double y, Values pars) { + double d = x - y; + if (abs(d - getPos(pars)) > cutoff * getW(pars)) { + return 0; + } + double pos = getPos(pars); + double w = getW(pars); + + switch (parName) { + case "shift": + return this.value(x, y, pars) * (d - pos) / w / w; + case "w": + return this.value(x, y, pars) * ((d - pos) * (d - pos) / w / w / w - 1 / w); + default: + return super.derivValue(parName, x, y, pars); + } + } } 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 deleted file mode 100644 index 1a68672c..00000000 --- a/numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java +++ /dev/null @@ -1,219 +0,0 @@ -/* - * 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.stat.parametric.AbstractParametricBiFunction; -import hep.dataforge.values.Values; - -import static java.lang.Math.*; - -/** - * A bi-function for beta-spectrum calculation taking energy and final state as - * input. - *

- * dissertation p.33 - *

- * - * @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); - } - - /** - * 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; - 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, Values 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; - } - - } - - /** - * The part independent of neutrino mass. Includes global normalization - * constant, momentum and Fermi correction - * - * @param E - * @return - */ - private double factor(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 K * res; - } - - @Override - public boolean providesDeriv(String name) { - return true; - } - - /** - * Bare beta spectrum with Mainz 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, Values 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": - case "U2": - case "msterile2": - return 0; - default: - return super.getDefaultParameter(name); - } - } - - @Override - public double derivValue(String parName, double fs, double eIn, Values pars) { - double E0 = getParameter("E0", pars); - return derivRootsterile(parName, eIn, E0 - fs, pars); - } - - @Override - public double value(double fs, double eIn, Values 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 deleted file mode 100644 index c93091a6..00000000 --- a/numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java +++ /dev/null @@ -1,99 +0,0 @@ -/* - * 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.context.Context; -import hep.dataforge.maths.MathPlugin; -import hep.dataforge.meta.Meta; -import hep.dataforge.stat.parametric.AbstractParametricBiFunction; -import hep.dataforge.values.Values; -import inr.numass.models.ResolutionFunction; -import inr.numass.utils.ExpressionUtils; -import org.apache.commons.math3.analysis.BivariateFunction; - -import java.util.HashMap; -import java.util.Map; - -import static java.lang.Math.sqrt; - -/** - * @author Alexander Nozik - */ -public class NumassResolution extends AbstractParametricBiFunction { - - private static final String[] list = {}; //leaving - - private final double resA; - private double resB = 0; - private BivariateFunction tailFunction = ResolutionFunction.getConstantTail(); - - public NumassResolution(Context context, Meta meta) { - super(list); - this.resA = meta.getDouble("A", 8.3e-5); - this.resB = meta.getDouble("B", 0); - if (meta.hasValue("tail")) { - String tailFunctionStr = meta.getString("tail"); - if (tailFunctionStr.startsWith("function::")) { - tailFunction = MathPlugin.buildFrom(context).buildBivariateFunction(tailFunctionStr.substring(10)); - } else { - tailFunction = (E, U) -> { - Map binding = new HashMap<>(); - binding.put("E", E); - binding.put("U", U); - binding.put("D", E - U); - return ExpressionUtils.function(tailFunctionStr, binding); - }; - } - } else if (meta.hasValue("tailAlpha")) { - //add polynomial function here - double alpha = meta.getDouble("tailAlpha"); - double beta = meta.getDouble("tailBeta", 0); - tailFunction = (double E, double U) -> 1 - (E - U) * (alpha + E / 1000d * beta) / 1000d; - - } else { - tailFunction = ResolutionFunction.getConstantTail(); - } - } - - @Override - public double derivValue(String parName, double x, double y, Values 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, Values set) { - assert resA > 0; - if (resB <= 0) { - 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 deleted file mode 100644 index d0992160..00000000 --- a/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java +++ /dev/null @@ -1,114 +0,0 @@ -/* - * 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.context.Context; -import hep.dataforge.maths.MathPlugin; -import hep.dataforge.meta.Meta; -import hep.dataforge.stat.parametric.AbstractParametricBiFunction; -import hep.dataforge.values.Values; -import inr.numass.models.LossCalculator; -import inr.numass.utils.ExpressionUtils; -import org.apache.commons.math3.analysis.BivariateFunction; -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.slf4j.LoggerFactory; - -import java.util.HashMap; -import java.util.Map; - -/** - * @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; - private Map lossCache = new HashMap<>(); - - private final boolean adjustX; - - public NumassTransmission(Context context, Meta meta) { - super(list); - this.calculator = LossCalculator.instance(); - adjustX = meta.getBoolean("adjustX", false); - if (meta.hasValue("trapping")) { - String trapFuncStr = meta.getString("trapping"); - if (trapFuncStr.startsWith("function::")) { - trapFunc = MathPlugin.buildFrom(context).buildBivariateFunction(trapFuncStr.substring(10)); - } else { - trapFunc = (Ei, Ef) -> { - Map binding = new HashMap<>(); - binding.put("Ei", Ei); - binding.put("Ef", Ef); - return ExpressionUtils.function(trapFuncStr, binding); - }; - } - } else { - LoggerFactory.getLogger(getClass()).warn("Trapping function not defined. Using default"); - trapFunc = MathPlugin.buildFrom(context).buildBivariateFunction("numass.trap.nominal"); - } - } - - public double getX(double eIn, Values set) { - if (adjustX) { - //From our article - return set.getDouble("X") * Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6; - } else { - return set.getDouble("X"); - } - } - - public double p0(double eIn, Values set) { - return LossCalculator.instance().getLossProbability(0, getX(eIn, set)); - } - - private BivariateFunction getTrapFunction(Context context, Meta meta) { - return MathPlugin.buildFrom(context).buildBivariateFunction(meta); - } - - @Override - public double derivValue(String parName, double eIn, double eOut, Values 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; - } - - @Override - public double value(double eIn, double eOut, Values set) { - //calculate X taking into account its energy dependence - double X = getX(eIn, set); - // loss part - double loss = calculator.getTotalLossValue(X, eIn, eOut); -// double loss; -// -// if(eIn-eOut >= 300){ -// loss = 0; -// } else { -// UnivariateFunction lossFunction = this.lossCache.computeIfAbsent(X, theX -> -// FunctionCaching.cacheUnivariateFunction(0, 300, 400, x -> calculator.getTotalLossValue(theX, eIn, eIn - x)) -// ); -// -// loss = lossFunction.value(eIn - eOut); -// } - - //trapping part - double trap = getParameter("trap", set) * trapFunc.value(eIn, eOut); - return loss + trap; - } - -} 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 deleted file mode 100644 index f2494b36..00000000 --- a/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java +++ /dev/null @@ -1,222 +0,0 @@ -/* - * 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.context.Context; -import hep.dataforge.context.Global; -import hep.dataforge.description.NodeDef; -import hep.dataforge.description.ValueDef; -import hep.dataforge.exceptions.NotDefinedException; -import hep.dataforge.maths.integration.UnivariateIntegrator; -import hep.dataforge.meta.Meta; -import hep.dataforge.stat.parametric.AbstractParametricBiFunction; -import hep.dataforge.stat.parametric.AbstractParametricFunction; -import hep.dataforge.stat.parametric.ParametricBiFunction; -import hep.dataforge.values.Values; -import inr.numass.models.FSS; -import inr.numass.utils.NumassIntegrator; -import org.apache.commons.math3.analysis.UnivariateFunction; - -import java.io.IOException; -import java.io.InputStream; - -import static hep.dataforge.values.ValueType.BOOLEAN; - -/** - * Compact all-in-one model for sterile neutrino spectrum - * - * @author Alexander Nozik - */ -@NodeDef(name = "resolution") -@NodeDef(name = "transmission") -@ValueDef(name = "fssFile", info = "The name for external FSS file. By default internal FSS file is used") -@ValueDef(name = "useFSS", type = {BOOLEAN}) -public class SterileNeutrinoSpectrum extends AbstractParametricFunction { - - private static final String[] list = {"X", "trap", "E0", "mnu2", "msterile2", "U2"}; - /** - * variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2" - */ - private final ParametricBiFunction source = new NumassBeta(); - /** - * variables:Ein,Eout; parameters: "A" - */ - private final NumassTransmission transmission; - /** - * variables:Eout,U; parameters: "X", "trap" - */ - private final ParametricBiFunction resolution; - /** - * auxiliary function for trans-res convolution - */ - private final ParametricBiFunction transRes; - private FSS fss; - // private boolean useMC; - private boolean fast; - - public SterileNeutrinoSpectrum(Context context, Meta configuration) { - super(list); - if (configuration.getBoolean("useFSS", true)) { - InputStream fssStream = configuration.optString("fssFile") - .map(fssFile -> { - try { - return context.getIo().optBinary(fssFile) - .orElseThrow(() -> new RuntimeException("Could not locate FSS file")) - .getStream(); - } catch (IOException e) { - throw new RuntimeException("Could not load FSS file", e); - } - }) - .orElse(getClass().getResourceAsStream("/data/FS.txt")); - - - fss = new FSS(fssStream); - } - - transmission = new NumassTransmission(context, configuration.getMetaOrEmpty("transmission")); - resolution = new NumassResolution(context, configuration.getMeta("resolution", Meta.empty())); - this.fast = configuration.getBoolean("fast", true); - transRes = new TransRes(); - } - - public SterileNeutrinoSpectrum(Meta configuration) { - this(Global.instance(), configuration); - } - - public SterileNeutrinoSpectrum() { - this(Global.instance(), Meta.empty()); - } - - @Override - public double derivValue(String parName, double u, Values set) { - switch (parName) { - case "U2": - case "msterile2": - case "mnu2": - case "E0": - return integrate(u, source.derivative(parName), transRes, set); - case "X": - case "trap": - return integrate(u, source, transRes.derivative(parName), set); - default: - throw new NotDefinedException(); - } - } - - @Override - public double value(double u, Values set) { - return integrate(u, source, transRes, set); - } - - @Override - public boolean providesDeriv(String name) { - return source.providesDeriv(name) && transmission.providesDeriv(name) && resolution.providesDeriv(name); - } - - - /** - * Direct Gauss-Legandre integration - * - * @param u - * @param sourceFunction - * @param transResFunction - * @param set - * @return - */ - private double integrate( - double u, - ParametricBiFunction sourceFunction, - ParametricBiFunction transResFunction, - Values 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), eIn, set); - } - return res; - }; - } else { - fsSource = (eIn) -> sourceFunction.value(0, eIn, set); - } - - UnivariateIntegrator integrator; - if (fast) { - if (eMax - u < 300) { - integrator = NumassIntegrator.getFastInterator(); - } else if (eMax - u > 2000) { - integrator = NumassIntegrator.getHighDensityIntegrator(); - } else { - integrator = NumassIntegrator.getDefaultIntegrator(); - } - - } else { - integrator = NumassIntegrator.getHighDensityIntegrator(); - } - - return integrator.integrate(u, eMax, eIn -> fsSource.value(eIn) * transResFunction.value(eIn, u, set)); - } - - private class TransRes extends AbstractParametricBiFunction { - - public TransRes() { - super(new String[]{"X", "trap"}); - } - - @Override - public boolean providesDeriv(String name) { - return true; - } - - @Override - public double derivValue(String parName, double eIn, double u, Values set) { - switch (parName) { - case "X": - //TODO implement p0 derivative - throw new NotDefinedException(); - case "trap": - return lossRes(transmission.derivative(parName), eIn, u, set); - default: - return super.derivValue(parName, eIn, u, set); - } - } - - @Override - public double value(double eIn, double u, Values set) { - - double p0 = transmission.p0(eIn, set); - return p0 * resolution.value(eIn, u, set) + lossRes(transmission, eIn, u, set); - } - - private double lossRes(ParametricBiFunction transFunc, double eIn, double u, Values set) { - UnivariateFunction integrand = (eOut) -> transFunc.value(eIn, eOut, set) * resolution.value(eOut, u, set); - - double border = u + 30; - double firstPart = NumassIntegrator.getFastInterator().integrate(u, Math.min(eIn, border), integrand); - double secondPart; - if (eIn > border) { - if (fast) { - secondPart = NumassIntegrator.getDefaultIntegrator().integrate(border, eIn, integrand); - } else { - secondPart = NumassIntegrator.getHighDensityIntegrator().integrate(border, eIn, integrand); - } - } else { - secondPart = 0; - } - return firstPart + secondPart; - } - - } - -} diff --git a/numass-main/src/main/java/inr/numass/models/sterile/TestModels.java b/numass-main/src/main/java/inr/numass/models/sterile/TestModels.java deleted file mode 100644 index 32c3de9a..00000000 --- a/numass-main/src/main/java/inr/numass/models/sterile/TestModels.java +++ /dev/null @@ -1,99 +0,0 @@ -/* - * 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.context.Context; -import hep.dataforge.maths.MathPlugin; -import hep.dataforge.meta.Meta; -import hep.dataforge.meta.MetaBuilder; -import hep.dataforge.stat.fit.ParamSet; -import hep.dataforge.stat.parametric.ParametricFunction; -import inr.numass.Numass; -import inr.numass.models.*; -import org.apache.commons.math3.analysis.BivariateFunction; - -/** - * - * @author Alexander Nozik - */ -public class TestModels { - - /** - * @param args the command line arguments - */ - public static void main(String[] args) { - Context context = Numass.buildContext(); - /* - - - - - */ - Meta meta = new MetaBuilder("model") - .putNode(new MetaBuilder("resolution") - .putValue("width", 1.22e-4) - .putValue("tailAlpha", 1e-2) - ) - .putNode(new MetaBuilder("transmission") - .putValue("trapping", "numass.trap2016") - ); - ParametricFunction oldFunc = oldModel(context, meta); - ParametricFunction newFunc = newModel(context, meta); - - ParamSet allPars = new ParamSet() - .setPar("N", 7e+05, 1.8e+03, 0d, Double.POSITIVE_INFINITY) - .setPar("bkg", 1d, 0.050) - .setPar("E0", 18575d, 1.4) - .setPar("mnu2", 0d, 1d) - .setPar("msterile2", 1000d * 1000d, 0) - .setPar("U2", 0.0, 1e-4, -1d, 1d) - .setPar("X", 0.04, 0.01, 0d, Double.POSITIVE_INFINITY) - .setPar("trap", 1, 0.01, 0d, Double.POSITIVE_INFINITY); - - for (double u = 14000; u < 18600; u += 100) { -// double oldVal = oldFunc.value(u, allPars); -// double newVal = newFunc.value(u, allPars); - double oldVal = oldFunc.derivValue("trap", u, allPars); - double newVal = newFunc.derivValue("trap", u, allPars); - System.out.printf("%f\t%g\t%g\t%g%n", u, oldVal, newVal, 1d - oldVal / newVal); - } - } - - @SuppressWarnings("deprecation") - private static ParametricFunction oldModel(Context context, Meta 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); - context.getChronicle().report("Setting up tritium model with real transmission function"); - BivariateFunction resolutionTail; - if (meta.hasValue("resolution.tailAlpha")) { - resolutionTail = ResolutionFunction.getAngledTail(meta.getDouble("resolution.tailAlpha"), meta.getDouble("resolution.tailBeta", 0)); - } else { - resolutionTail = ResolutionFunction.getRealTail(); - } - //RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); - RangedNamedSetSpectrum beta = new BetaSpectrum(); - ModularSpectrum sp = new ModularSpectrum(beta, new ResolutionFunction(A, resolutionTail), from, to); - if (meta.getBoolean("caching", false)) { - context.getChronicle().report("Caching turned on"); - sp.setCaching(true); - } - //Adding trapping energy dependence - - if (meta.hasValue("transmission.trapping")) { - BivariateFunction trap = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("transmission.trapping")); - sp.setTrappingFunction(trap); - } - - return new NBkgSpectrum(sp); - } - - private static ParametricFunction newModel(Context context, Meta meta) { - SterileNeutrinoSpectrum sp = new SterileNeutrinoSpectrum(context, meta); - return new NBkgSpectrum(sp); - } - -} diff --git a/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java b/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java index 73c64e56..096d49d8 100644 --- a/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java +++ b/numass-main/src/main/java/inr/numass/utils/DataModelUtils.java @@ -31,7 +31,7 @@ public class DataModelUtils { public static Table getUniformSpectrumConfiguration(double from, double to, double time, int numpoints) { assert to != from; - final String[] list = {SpectrumDataAdapter.X_VALUE_KEY, "time"}; + final String[] list = {SpectrumDataAdapter.X_AXIS, "time"}; ListTable.Builder res = new ListTable.Builder(list); for (int i = 0; i < numpoints; i++) { @@ -45,7 +45,7 @@ public class DataModelUtils { } public static Table getSpectrumConfigurationFromResource(String resource) { - final String[] list = {"x", "time"}; + final String[] list = {SpectrumDataAdapter.X_AXIS, "time"}; ListTable.Builder res = new ListTable.Builder(list); Scanner scan = new Scanner(DataModelUtils.class.getResourceAsStream(resource)); while (scan.hasNextLine()) { diff --git a/numass-main/src/main/java/inr/numass/utils/NumassIntegrator.java b/numass-main/src/main/java/inr/numass/utils/NumassIntegrator.java index 529f3804..d3bf36d0 100644 --- a/numass-main/src/main/java/inr/numass/utils/NumassIntegrator.java +++ b/numass-main/src/main/java/inr/numass/utils/NumassIntegrator.java @@ -20,7 +20,7 @@ public class NumassIntegrator { private static UnivariateIntegrator defaultIntegrator; private static UnivariateIntegrator highDensityIntegrator; - public synchronized static UnivariateIntegrator getFastInterator() { + public synchronized static UnivariateIntegrator getFastInterator() { if (fastInterator == null) { LoggerFactory.getLogger(NumassIntegrator.class).debug("Creating fast integrator"); fastInterator = new GaussRuleIntegrator((int) (mult * 100)); @@ -28,7 +28,7 @@ public class NumassIntegrator { return fastInterator; } - public synchronized static UnivariateIntegrator getDefaultIntegrator() { + public synchronized static UnivariateIntegrator getDefaultIntegrator() { if (defaultIntegrator == null) { LoggerFactory.getLogger(NumassIntegrator.class).debug("Creating default integrator"); defaultIntegrator = new GaussRuleIntegrator((int) (mult * 300)); @@ -36,12 +36,11 @@ public class NumassIntegrator { return defaultIntegrator; } - public synchronized static UnivariateIntegrator getHighDensityIntegrator() { + public synchronized static UnivariateIntegrator getHighDensityIntegrator() { if (highDensityIntegrator == null) { LoggerFactory.getLogger(NumassIntegrator.class).debug("Creating high precision integrator"); highDensityIntegrator = new GaussRuleIntegrator((int) (mult * 500)); } return highDensityIntegrator; } - } diff --git a/numass-main/src/main/kotlin/inr/numass/NumassUtils.kt b/numass-main/src/main/kotlin/inr/numass/NumassUtils.kt index 9eac8e6e..5643d787 100644 --- a/numass-main/src/main/kotlin/inr/numass/NumassUtils.kt +++ b/numass-main/src/main/kotlin/inr/numass/NumassUtils.kt @@ -35,6 +35,7 @@ import hep.dataforge.values.Values import inr.numass.data.api.NumassAnalyzer import inr.numass.data.api.NumassPoint import inr.numass.data.api.NumassSet +import inr.numass.models.FSS import inr.numass.utils.ExpressionUtils import org.apache.commons.math3.analysis.UnivariateFunction import org.jfree.chart.plot.IntervalMarker @@ -43,6 +44,7 @@ import tornadofx.* import java.awt.Color import java.awt.Font import java.io.IOException +import java.io.InputStream import java.io.OutputStream import java.lang.Math.* import java.util.* @@ -160,6 +162,26 @@ object NumassUtils { } +fun getFSS(context: Context, meta: Meta): FSS? { + return if (meta.getBoolean("useFSS", true)) { + val fssStream = meta.optString("fssFile") + .map { fssFile -> + try { + context.io.optBinary(fssFile) + .orElseThrow({ RuntimeException("Could not locate FSS file") }) + .stream + } catch (e: IOException) { + throw RuntimeException("Could not load FSS file", e) + } + } + .orElse(context.io.optResource("data/FS.txt").get().stream) + FSS(fssStream) + } else { + null + } +} + + /** * Evaluate groovy expression using numass point as parameter * @@ -200,7 +222,7 @@ fun addSetMarkers(frame: JFreeChartFrame, sets: Collection) { /** * Subtract one energy spectrum from the other one */ -fun subtract(context: Context, merge: Table, empty: Table): Table { +fun subtractAmplitudeSpectrum(context: Context, merge: Table, empty: Table): Table { val builder = ListTable.Builder(merge.format) merge.rows.forEach { point -> val pointBuilder = ValueMap.Builder(point) diff --git a/numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt b/numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt index 02eef9a7..2bd4fc6f 100644 --- a/numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt +++ b/numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt @@ -1,11 +1,125 @@ package inr.numass.models +import hep.dataforge.maths.integration.UnivariateIntegrator +import hep.dataforge.names.Names import hep.dataforge.stat.models.Model import hep.dataforge.stat.models.ModelDescriptor import hep.dataforge.stat.models.ModelFactory +import hep.dataforge.stat.parametric.AbstractParametricFunction +import hep.dataforge.stat.parametric.ParametricFunction import hep.dataforge.utils.ContextMetaFactory +import hep.dataforge.values.Values +import inr.numass.utils.NumassIntegrator +import java.util.stream.Stream fun model(name: String, descriptor: ModelDescriptor? = null, factory: ContextMetaFactory): ModelFactory { return ModelFactory.build(name, descriptor, factory); +} + +// spectra operations +//TODO move to stat + +/** + * Calculate sum of two parametric functions + */ +operator fun ParametricFunction.plus(func: ParametricFunction): ParametricFunction { + val mergedNames = Names.of(Stream.concat(names.stream(), func.names.stream()).distinct()) + return object : AbstractParametricFunction(mergedNames) { + override fun derivValue(parName: String, x: Double, set: Values): Double { + return func.derivValue(parName, x, set) + this@plus.derivValue(parName, x, set) + } + + override fun value(x: Double, set: Values): Double { + return this@plus.value(x, set) + func.value(x, set) + } + + override fun providesDeriv(name: String): Boolean { + return this@plus.providesDeriv(name) && func.providesDeriv(name) + } + + } +} + +operator fun ParametricFunction.minus(func: ParametricFunction): ParametricFunction { + val mergedNames = Names.of(Stream.concat(names.stream(), func.names.stream()).distinct()) + return object : AbstractParametricFunction(mergedNames) { + override fun derivValue(parName: String, x: Double, set: Values): Double { + return func.derivValue(parName, x, set) - this@minus.derivValue(parName, x, set) + } + + override fun value(x: Double, set: Values): Double { + return this@minus.value(x, set) - func.value(x, set) + } + + override fun providesDeriv(name: String): Boolean { + return this@minus.providesDeriv(name) && func.providesDeriv(name) + } + + } +} + +/** + * Calculate product of two parametric functions + * + */ +operator fun ParametricFunction.times(func: ParametricFunction): ParametricFunction { + val mergedNames = Names.of(Stream.concat(names.stream(), func.names.stream()).distinct()) + return object : AbstractParametricFunction(mergedNames) { + override fun derivValue(parName: String, x: Double, set: Values): Double { + return this@times.value(x, set) * func.derivValue(parName, x, set) + this@times.derivValue(parName, x, set) * func.value(x, set) + } + + override fun value(x: Double, set: Values): Double { + return this@times.value(x, set) * func.value(x, set) + } + + override fun providesDeriv(name: String): Boolean { + return this@times.providesDeriv(name) && func.providesDeriv(name) + } + + } +} + +/** + * Multiply parametric function by fixed value + */ +operator fun ParametricFunction.times(num: Number): ParametricFunction { + return object : AbstractParametricFunction(names) { + override fun derivValue(parName: String, x: Double, set: Values): Double { + return this@times.value(x, set) * num.toDouble() + } + + override fun value(x: Double, set: Values): Double { + return this@times.value(x, set) * num.toDouble() + } + + override fun providesDeriv(name: String): Boolean { + return this@times.providesDeriv(name) + } + + } +} + +/** + * Calculate convolution of two parametric functions + */ +fun ParametricFunction.convolute(func: ParametricFunction, + integrator: UnivariateIntegrator<*> = NumassIntegrator.getDefaultIntegrator()): ParametricFunction { + val mergedNames = Names.of(Stream.concat(names.stream(), func.names.stream()).distinct()) + return object : AbstractParametricFunction(mergedNames){ + override fun derivValue(parName: String, x: Double, set: Values): Double { + TODO("not implemented") //To change body of created functions use File | Settings | File Templates. + } + + override fun value(x: Double, set: Values): Double { + TODO("not implemented") //To change body of created functions use File | Settings | File Templates. + } + + override fun providesDeriv(name: String?): Boolean { + return this@convolute.providesDeriv(name) && func.providesDeriv(name) + } + + } + } \ No newline at end of file diff --git a/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassBeta.kt b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassBeta.kt new file mode 100644 index 00000000..0b7ead19 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassBeta.kt @@ -0,0 +1,230 @@ +/* + * 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.stat.parametric.AbstractParametricBiFunction +import hep.dataforge.stat.parametric.AbstractParametricFunction +import hep.dataforge.stat.parametric.ParametricFunction +import hep.dataforge.values.Values + +import java.lang.Math.* + +/** + * A bi-function for beta-spectrum calculation taking energy and final state as + * input. + * + * + * dissertation p.33 + * + * + * @author [Alexander Nozik](mailto:altavir@gmail.com) + */ +class NumassBeta : AbstractParametricBiFunction(list) { + + /** + * Beta spectrum derivative + * + * @param n parameter number + * @param E0 + * @param mnu2 + * @param E + * @return + * @throws NotDefinedException + */ + @Throws(NotDefinedException::class) + private fun derivRoot(n: Int, E0: Double, mnu2: Double, E: Double): Double { + val D = E0 - E//E0-E + if (D == 0.0) { + return 0.0 + } + + return if (mnu2 >= 0) { + if (E >= E0 - sqrt(mnu2)) { + 0.0 + } else { + val bare = sqrt(D * D - mnu2) + when (n) { + 0 -> factor(E) * (2.0 * D * D - mnu2) / bare + 1 -> -factor(E) * 0.5 * D / bare + else -> 0.0 + } + } + } else { + val mu = sqrt(-0.66 * mnu2) + if (E >= E0 + mu) { + 0.0 + } else { + val root = sqrt(Math.max(D * D - mnu2, 0.0)) + val exp = exp(-1 - D / mu) + when (n) { + 0 -> factor(E) * (D * (D + mu * exp) / root + root * (1 - exp)) + 1 -> factor(E) * (-(D + mu * exp) / root * 0.5 - root * exp * (1 + D / mu) / 3.0 / mu) + else -> 0.0 + } + } + } + } + + /** + * Derivative of spectrum with sterile neutrinos + * + * @param name + * @param E + * @param E0 + * @param pars + * @return + * @throws NotDefinedException + */ + @Throws(NotDefinedException::class) + private fun derivRootsterile(name: String, E: Double, E0: Double, pars: Values): Double { + val mnu2 = getParameter("mnu2", pars) + val mst2 = getParameter("msterile2", pars) + val u2 = getParameter("U2", pars) + + return when (name) { + "E0" -> { + if (u2 == 0.0) { + derivRoot(0, E0, mnu2, E) + } else { + u2 * derivRoot(0, E0, mst2, E) + (1 - u2) * derivRoot(0, E0, mnu2, E) + } + } + "mnu2" -> (1 - u2) * derivRoot(1, E0, mnu2, E) + "msterile2" -> { + if (u2 == 0.0) { + 0.0 + } else { + u2 * derivRoot(1, E0, mst2, E) + } + } + "U2" -> root(E0, mst2, E) - root(E0, mnu2, E) + else -> 0.0 + } + + } + + /** + * The part independent of neutrino mass. Includes global normalization + * constant, momentum and Fermi correction + * + * @param E + * @return + */ + private fun factor(E: Double): Double { + val me = 0.511006E6 + val eTot = E + me + val pe = sqrt(E * (E + 2.0 * me)) + val ve = pe / eTot + val yfactor = 2.0 * 2.0 * 1.0 / 137.039 * Math.PI + val y = yfactor / ve + val Fn = y / abs(1.0 - exp(-y)) + val Fermi = Fn * (1.002037 - 0.001427 * ve) + val res = Fermi * pe * eTot + return K * res + } + + override fun providesDeriv(name: String): Boolean { + return true + } + + /** + * Bare beta spectrum with Mainz negative mass correction + * + * @param E0 + * @param mnu2 + * @param E + * @return + */ + private fun root(E0: Double, mnu2: Double, E: Double): Double { + //bare beta-spectrum + val delta = E0 - E//E0-E + val res: Double + val bare = factor(E) * delta * sqrt(Math.max(delta * delta - mnu2, 0.0)) + if (delta == 0.0) { + return 0.0 + } + if (mnu2 >= 0) { + res = Math.max(bare, 0.0) + } else { + if (delta + 0.812 * sqrt(-mnu2) <= 0) { + return 0.0 //sqrt(0.66) + } + val aux = sqrt(-mnu2 * 0.66) / delta + res = Math.max(bare * (1 + aux * exp(-1 - 1 / aux)), 0.0) + } + return res + } + + /** + * beta-spectrum with sterile neutrinos + * + * @param E + * @param E0 + * @param pars + * @return + */ + private fun rootsterile(E: Double, E0: Double, pars: Values): Double { + val mnu2 = getParameter("mnu2", pars) + val mst2 = getParameter("msterile2", pars) + val u2 = getParameter("U2", pars) + + return if (u2 == 0.0) { + root(E0, mnu2, E) + } else { + u2 * root(E0, mst2, E) + (1 - u2) * root(E0, mnu2, E) + } +// P(rootsterile)+ (1-P)root + } + + override fun getDefaultParameter(name: String): Double { + when (name) { + "mnu2", "U2", "msterile2" -> return 0.0 + else -> return super.getDefaultParameter(name) + } + } + + override fun derivValue(parName: String, fs: Double, eIn: Double, pars: Values): Double { + val e0 = getParameter("E0", pars) + return derivRootsterile(parName, eIn, e0 - fs, pars) + } + + override fun value(fs: Double, eIn: Double, pars: Values): Double { + val e0 = getParameter("E0", pars) + return rootsterile(eIn, e0 - fs, pars) + } + + /** + * Get univariate spectrum with given final state + */ + fun getSpectrum(fs: Double = 0.0): ParametricFunction { + return BetaSpectrum(fs); + } + + inner class BetaSpectrum(val fs: Double) : AbstractParametricFunction(list) { + + override fun providesDeriv(name: String): Boolean { + return this@NumassBeta.providesDeriv(name) + } + + override fun derivValue(parName: String, x: Double, set: Values): Double { + return this@NumassBeta.derivValue(parName, fs, x, set) + } + + override fun value(x: Double, set: Values): Double { + return this@NumassBeta.value(fs, x, set) + } + + } + + + companion object { + + private val K = 1E-23 + private val list = arrayOf("E0", "mnu2", "msterile2", "U2") + } + +} diff --git a/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassResolution.kt b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassResolution.kt new file mode 100644 index 00000000..30a89e04 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassResolution.kt @@ -0,0 +1,91 @@ +/* + * 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.context.Context +import hep.dataforge.maths.MathPlugin +import hep.dataforge.meta.Meta +import hep.dataforge.stat.parametric.AbstractParametricBiFunction +import hep.dataforge.values.Values +import inr.numass.models.ResolutionFunction +import inr.numass.utils.ExpressionUtils +import org.apache.commons.math3.analysis.BivariateFunction +import java.lang.Math.sqrt +import java.util.* + +/** + * @author [Alexander Nozik](mailto:altavir@gmail.com) + */ +class NumassResolution(context: Context, meta: Meta) : AbstractParametricBiFunction(list) { + + private val resA: Double = meta.getDouble("A", 8.3e-5) + private val resB = meta.getDouble("B", 0.0) + private val tailFunction: BivariateFunction = when { + meta.hasValue("tail") -> { + val tailFunctionStr = meta.getString("tail") + if (tailFunctionStr.startsWith("function::")) { + MathPlugin.buildFrom(context).buildBivariateFunction(tailFunctionStr.substring(10)) + } else { + BivariateFunction { E, U -> + val binding = HashMap() + binding.put("E", E) + binding.put("U", U) + binding.put("D", E - U) + ExpressionUtils.function(tailFunctionStr, binding) + } + } + } + meta.hasValue("tailAlpha") -> { + //add polynomial function here + val alpha = meta.getDouble("tailAlpha")!! + val beta = meta.getDouble("tailBeta", 0.0)!! + BivariateFunction { E: Double, U: Double -> 1 - (E - U) * (alpha + E / 1000.0 * beta) / 1000.0 } + + } + else -> ResolutionFunction.getConstantTail() + } + + override fun derivValue(parName: String, x: Double, y: Double, set: Values): Double { + return 0.0 + } + + private fun getValueFast(E: Double, U: Double): Double { + val delta = resA * E + return if (E - U < 0) { + 0.0 + } else if (E - U > delta) { + tailFunction.value(E, U) + } else { + (E - U) / delta + } + } + + override fun providesDeriv(name: String): Boolean { + return true + } + + override fun value(E: Double, U: Double, set: Values): Double { + assert(resA > 0) + if (resB <= 0) { + return this.getValueFast(E, U) + } + assert(resB > 0) + val delta = resA * E + return if (E - U < 0) { + 0.0 + } else if (E - U > delta) { + tailFunction.value(E, U) + } else { + (1 - sqrt(1 - (E - U) / E * resB)) / (1 - sqrt(1 - resA * resB)) + } + } + + companion object { + + private val list = arrayOf() //leaving + } + +} diff --git a/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt new file mode 100644 index 00000000..85a96f5a --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt @@ -0,0 +1,105 @@ +/* + * 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.context.Context +import hep.dataforge.maths.MathPlugin +import hep.dataforge.meta.Meta +import hep.dataforge.stat.parametric.AbstractParametricBiFunction +import hep.dataforge.values.Values +import inr.numass.models.LossCalculator +import inr.numass.utils.ExpressionUtils +import org.apache.commons.math3.analysis.BivariateFunction +import org.slf4j.LoggerFactory +import java.util.* + +/** + * @author [Alexander Nozik](mailto:altavir@gmail.com) + */ +class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFunction(list) { + private val calculator: LossCalculator = LossCalculator.instance() + private val trapFunc: BivariateFunction + //private val lossCache = HashMap() + + private val adjustX: Boolean = meta.getBoolean("adjustX", false)!! + + init { + if (meta.hasValue("trapping")) { + val trapFuncStr = meta.getString("trapping") + trapFunc = if (trapFuncStr.startsWith("function::")) { + MathPlugin.buildFrom(context).buildBivariateFunction(trapFuncStr.substring(10)) + } else { + BivariateFunction { Ei: Double, Ef: Double -> + val binding = HashMap() + binding.put("Ei", Ei) + binding.put("Ef", Ef) + ExpressionUtils.function(trapFuncStr, binding) + } + } + } else { + LoggerFactory.getLogger(javaClass).warn("Trapping function not defined. Using default") + trapFunc = MathPlugin.buildFrom(context).buildBivariateFunction("numass.trap.nominal") + } + } + + fun getX(eIn: Double, set: Values): Double { + return if (adjustX) { + //From our article + set.getDouble("X") * Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6 + } else { + set.getDouble("X") + } + } + + fun p0(eIn: Double, set: Values): Double { + return LossCalculator.instance().getLossProbability(0, getX(eIn, set)) + } + + private fun getTrapFunction(context: Context, meta: Meta): BivariateFunction { + return MathPlugin.buildFrom(context).buildBivariateFunction(meta) + } + + override fun derivValue(parName: String, eIn: Double, eOut: Double, set: Values): Double { + return when (parName) { + "trap" -> trapFunc.value(eIn, eOut) + "X" -> calculator.getTotalLossDeriv(getX(eIn, set), eIn, eOut) + else -> super.derivValue(parName, eIn, eOut, set) + } + } + + override fun providesDeriv(name: String): Boolean { + return true + } + + override fun value(eIn: Double, eOut: Double, set: Values): Double { + //calculate X taking into account its energy dependence + val X = getX(eIn, set) + // loss part + val loss = calculator.getTotalLossValue(X, eIn, eOut) + // double loss; + // + // if(eIn-eOut >= 300){ + // loss = 0; + // } else { + // UnivariateFunction lossFunction = this.lossCache.computeIfAbsent(X, theX -> + // FunctionCaching.cacheUnivariateFunction(0, 300, 400, x -> calculator.getTotalLossValue(theX, eIn, eIn - x)) + // ); + // + // loss = lossFunction.value(eIn - eOut); + // } + + //trapping part + val trap = getParameter("trap", set) * trapFunc.value(eIn, eOut) + return loss + trap + } + + companion object { + + private val list = arrayOf("trap", "X") + private val ION_POTENTIAL = 15.4//eV + } + +} diff --git a/numass-main/src/main/kotlin/inr/numass/models/sterile/SterileNeutrinoSpectrum.kt b/numass-main/src/main/kotlin/inr/numass/models/sterile/SterileNeutrinoSpectrum.kt new file mode 100644 index 00000000..d416b793 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/sterile/SterileNeutrinoSpectrum.kt @@ -0,0 +1,164 @@ +/* + * 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.context.Context +import hep.dataforge.description.NodeDef +import hep.dataforge.description.NodeDefs +import hep.dataforge.description.ValueDef +import hep.dataforge.description.ValueDefs +import hep.dataforge.exceptions.NotDefinedException +import hep.dataforge.maths.integration.UnivariateIntegrator +import hep.dataforge.meta.Meta +import hep.dataforge.stat.parametric.AbstractParametricBiFunction +import hep.dataforge.stat.parametric.AbstractParametricFunction +import hep.dataforge.stat.parametric.ParametricBiFunction +import hep.dataforge.values.ValueType.BOOLEAN +import hep.dataforge.values.Values +import inr.numass.models.FSS +import inr.numass.utils.NumassIntegrator + +/** + * Compact all-in-one model for sterile neutrino spectrum + * + * @author Alexander Nozik + */ +@NodeDefs( + NodeDef(name = "resolution"), + NodeDef(name = "transmission") +) +@ValueDefs( + ValueDef(name = "fssFile", info = "The name for external FSS file. By default internal FSS file is used"), + ValueDef(name = "useFSS", type = arrayOf(BOOLEAN)) +) + +/** + * @param source variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2" + * @param transmission variables:Ein,Eout; parameters: "A" + * @param resolution variables:Eout,U; parameters: "X", "trap" + */ + +class SterileNeutrinoSpectrum @JvmOverloads constructor( + context: Context, + configuration: Meta, + val source: ParametricBiFunction = NumassBeta(), + val transmission: NumassTransmission = NumassTransmission(context, configuration.getMetaOrEmpty("transmission")), + val resolution: ParametricBiFunction = NumassResolution(context, configuration.getMeta("resolution", Meta.empty()))) : AbstractParametricFunction(list) { + + + /** + * auxiliary function for trans-res convolution + */ + private val transRes: ParametricBiFunction = TransRes() + private val fss: FSS? = getFSS(context, configuration) + // private boolean useMC; + private val fast: Boolean = configuration.getBoolean("fast", true) + + override fun derivValue(parName: String, u: Double, set: Values): Double { + return when (parName) { + "U2", "msterile2", "mnu2", "E0" -> integrate(u, source.derivative(parName), transRes, set) + "X", "trap" -> integrate(u, source, transRes.derivative(parName), set) + else -> throw NotDefinedException() + } + } + + override fun value(u: Double, set: Values): Double { + return integrate(u, source, transRes, set) + } + + override fun providesDeriv(name: String): Boolean { + return source.providesDeriv(name) && transmission.providesDeriv(name) && resolution.providesDeriv(name) + } + + + /** + * Direct Gauss-Legandre integration + * + * @param u + * @param sourceFunction + * @param transResFunction + * @param set + * @return + */ + private fun integrate( + u: Double, + sourceFunction: ParametricBiFunction, + transResFunction: ParametricBiFunction, + set: Values): Double { + + val eMax = set.getDouble("E0") + 5.0 + + if (u >= eMax) { + return 0.0 + } + + + val fsSource: (Double) -> Double = fss?.let { fss -> + { eIn: Double -> + (0 until fss.size()).sumByDouble { fss.getP(it) * sourceFunction.value(fss.getE(it), eIn, set) } + } + } ?: { eIn: Double -> sourceFunction.value(0.0, eIn, set) } + + + val integrator: UnivariateIntegrator<*> = if (fast) { + when { + eMax - u < 300 -> NumassIntegrator.getFastInterator() + eMax - u > 2000 -> NumassIntegrator.getHighDensityIntegrator() + else -> NumassIntegrator.getDefaultIntegrator() + } + + } else { + NumassIntegrator.getHighDensityIntegrator() + } + + return integrator.integrate(u, eMax) { eIn -> fsSource(eIn) * transResFunction.value(eIn, u, set) } + } + + private inner class TransRes : AbstractParametricBiFunction(arrayOf("X", "trap")) { + + override fun providesDeriv(name: String): Boolean { + return true + } + + override fun derivValue(parName: String, eIn: Double, u: Double, set: Values): Double { + return when (parName) { + "X" -> throw NotDefinedException()//TODO implement p0 derivative + "trap" -> lossRes(transmission.derivative(parName), eIn, u, set) + else -> super.derivValue(parName, eIn, u, set) + } + } + + override fun value(eIn: Double, u: Double, set: Values): Double { + + val p0 = transmission.p0(eIn, set) + return p0 * resolution.value(eIn, u, set) + lossRes(transmission, eIn, u, set) + } + + private fun lossRes(transFunc: ParametricBiFunction, eIn: Double, u: Double, set: Values): Double { + val integrand = { eOut: Double -> transFunc.value(eIn, eOut, set) * resolution.value(eOut, u, set) } + + val border = u + 30 + val firstPart = NumassIntegrator.getFastInterator().integrate(u, Math.min(eIn, border), integrand) + val secondPart: Double = if (eIn > border) { + if (fast) { + NumassIntegrator.getDefaultIntegrator().integrate(border, eIn, integrand) + } else { + NumassIntegrator.getHighDensityIntegrator().integrate(border, eIn, integrand) + } + } else { + 0.0 + } + return firstPart + secondPart + } + + } + + companion object { + + private val list = arrayOf("X", "trap", "E0", "mnu2", "msterile2", "U2") + } + +} diff --git a/numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt b/numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt new file mode 100644 index 00000000..c522ae66 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt @@ -0,0 +1,104 @@ +/* + * 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.context.Context +import hep.dataforge.maths.MathPlugin +import hep.dataforge.meta.Meta +import hep.dataforge.meta.MetaBuilder +import hep.dataforge.stat.fit.ParamSet +import hep.dataforge.stat.parametric.ParametricFunction +import inr.numass.Numass +import inr.numass.models.BetaSpectrum +import inr.numass.models.ModularSpectrum +import inr.numass.models.NBkgSpectrum +import inr.numass.models.ResolutionFunction +import org.apache.commons.math3.analysis.BivariateFunction + +/** + * + * @author Alexander Nozik + */ +object TestModels { + + /** + * @param args the command line arguments + */ + @JvmStatic + fun main(args: Array) { + val context = Numass.buildContext() + /* + + + + + */ + val meta = MetaBuilder("model") + .putNode(MetaBuilder("resolution") + .putValue("width", 1.22e-4) + .putValue("tailAlpha", 1e-2) + ) + .putNode(MetaBuilder("transmission") + .putValue("trapping", "numass.trap2016") + ) + val oldFunc = oldModel(context, meta) + val newFunc = newModel(context, meta) + + val allPars = ParamSet() + .setPar("N", 7e+05, 1.8e+03, 0.0, java.lang.Double.POSITIVE_INFINITY) + .setPar("bkg", 1.0, 0.050) + .setPar("E0", 18575.0, 1.4) + .setPar("mnu2", 0.0, 1.0) + .setPar("msterile2", 1000.0 * 1000.0, 0.0) + .setPar("U2", 0.0, 1e-4, -1.0, 1.0) + .setPar("X", 0.04, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY) + .setPar("trap", 1.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY) + + var u = 14000.0 + while (u < 18600) { + // double oldVal = oldFunc.value(u, allPars); + // double newVal = newFunc.value(u, allPars); + val oldVal = oldFunc.derivValue("trap", u, allPars) + val newVal = newFunc.derivValue("trap", u, allPars) + System.out.printf("%f\t%g\t%g\t%g%n", u, oldVal, newVal, 1.0 - oldVal / newVal) + u += 100.0 + } + } + + private fun oldModel(context: Context, meta: Meta): ParametricFunction { + val A = meta.getDouble("resolution", meta.getDouble("resolution.width", 8.3e-5)!!)!!//8.3e-5 + val from = meta.getDouble("from", 13900.0)!! + val to = meta.getDouble("to", 18700.0)!! + context.chronicle.report("Setting up tritium model with real transmission function") + val resolutionTail: BivariateFunction + if (meta.hasValue("resolution.tailAlpha")) { + resolutionTail = ResolutionFunction.getAngledTail(meta.getDouble("resolution.tailAlpha")!!, meta.getDouble("resolution.tailBeta", 0.0)!!) + } else { + resolutionTail = ResolutionFunction.getRealTail() + } + //RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); + val beta = BetaSpectrum() + val sp = ModularSpectrum(beta, ResolutionFunction(A, resolutionTail), from, to) + if (meta.getBoolean("caching", false)!!) { + context.chronicle.report("Caching turned on") + sp.setCaching(true) + } + //Adding trapping energy dependence + + if (meta.hasValue("transmission.trapping")) { + val trap = MathPlugin.buildFrom(context).buildBivariateFunction(meta.getString("transmission.trapping")) + sp.setTrappingFunction(trap) + } + + return NBkgSpectrum(sp) + } + + private fun newModel(context: Context, meta: Meta): ParametricFunction { + val sp = SterileNeutrinoSpectrum(context, meta) + return NBkgSpectrum(sp) + } + +} diff --git a/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt b/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt index 079ebb3c..2f230643 100644 --- a/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt +++ b/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt @@ -34,7 +34,7 @@ import inr.numass.addSetMarkers import inr.numass.data.analyzers.SmartAnalyzer import inr.numass.data.api.NumassPoint import inr.numass.data.api.NumassSet -import inr.numass.subtract +import inr.numass.subtractAmplitudeSpectrum import inr.numass.unbox import inr.numass.utils.ExpressionUtils import java.io.PrintWriter @@ -169,7 +169,7 @@ val subtractEmptyTask = task("dif") { node("empty", empty.meta) } val res = DataUtils.combine(input, empty, Table::class.java, resMeta) { mergeData, emptyData -> - subtract(context, mergeData, emptyData) + subtractAmplitudeSpectrum(context, mergeData, emptyData) } res.goal.onComplete { r, _ ->