From 53aad7f6a4bdb3f5b184fb3f4a348b4893393ed1 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 6 Jul 2018 22:07:16 +0300 Subject: [PATCH] Fixed fit internals --- .../inr/numass/scripts/FindExIonRatio.groovy | 38 +- .../numass/scripts/LossTailCalculation.groovy | 1 - .../numass/scripts/PrintLossFunctions.groovy | 5 +- .../inr/numass/scripts/Systematics.groovy | 33 +- .../numass/models/EmpiricalLossSpectrum.java | 7 +- .../inr/numass/models/LossCalculator.java | 442 ------------------ .../inr/numass/models/ModularSpectrum.java | 6 +- .../numass/models/VariableLossSpectrum.java | 7 +- .../main/kotlin/inr/numass/NumassPlugin.kt | 22 +- .../kotlin/inr/numass/models/NumassModels.kt | 5 +- .../inr/numass/models/misc/LossCalculator.kt | 109 +++-- .../models/sterile/NumassTransmission.kt | 9 +- .../numass/scripts/models/IntegralSpectrum.kt | 25 +- .../numass/scripts/tristan/AnalyzeTristan.kt | 4 +- 14 files changed, 137 insertions(+), 576 deletions(-) delete mode 100644 numass-main/src/main/java/inr/numass/models/LossCalculator.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 86d426d7..45fb8ab6 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy @@ -20,16 +20,16 @@ import hep.dataforge.plots.PlotFrame import hep.dataforge.plots.XYFunctionPlot import hep.dataforge.plots.jfreechart.JFreeChartFrame import hep.dataforge.stat.fit.ParamSet -import inr.numass.models.LossCalculator +import inr.numass.models.misc.LossCalculator import org.apache.commons.math3.analysis.UnivariateFunction import org.apache.commons.math3.analysis.solvers.BisectionSolver ParamSet params = new ParamSet() -.setParValue("exPos", 12.76) -.setParValue("ionPos", 13.95) -.setParValue("exW", 1.2) -.setParValue("ionW", 13.5) -.setParValue("exIonRatio", 4.55) + .setParValue("exPos", 12.76) + .setParValue("ionPos", 13.95) + .setParValue("exW", 1.2) + .setParValue("ionW", 13.5) + .setParValue("exIonRatio", 4.55) @@ -43,7 +43,7 @@ UnivariateIntegrator integrator = NumassContext.defaultIntegrator; double border = 13.6; -UnivariateFunction ratioFunction = {e->integrator.integrate(0, e, scatterFunction) / integrator.integrate(e, 100, scatterFunction)} +UnivariateFunction ratioFunction = { e -> integrator.integrate(0, e, scatterFunction) / integrator.integrate(e, 100, scatterFunction) } double ratio = ratioFunction.value(border); println "The true excitation to ionization ratio with border energy $border is $ratio"; @@ -54,35 +54,35 @@ double resolution = 1.5d; def X = 0.527; -LossCalculator calculator = new LossCalculator(); +LossCalculator calculator = LossCalculator.INSTANCE; List lossProbs = calculator.getGunLossProbabilities(X); -UnivariateFunction newScatterFunction = { double d -> +UnivariateFunction newScatterFunction = { double d -> double res = scatterFunction.value(d); - for(i = 1; i < lossProbs.size(); i++){ + for (i = 1; i < lossProbs.size(); i++) { res += lossProbs.get(i) * calculator.getLossValue(i, d, 0); } return res; } -UnivariateFunction resolutionValue = {double e -> +UnivariateFunction resolutionValue = { double e -> if (e <= 0d) { return 0d; } else if (e >= resolution) { return 1d; } else { - return e/resolution; + return e / resolution; } }; -UnivariateFunction integral = {double u -> - if(u <= 0d){ +UnivariateFunction integral = { double u -> + if (u <= 0d) { return 0d; } else { - UnivariateFunction integrand = {double e -> resolutionValue.value(u-e) * newScatterFunction.value(e)}; + UnivariateFunction integrand = { double e -> resolutionValue.value(u - e) * newScatterFunction.value(e) }; return integrator.integrate(0d, u, integrand) } } @@ -92,9 +92,9 @@ frame.add(XYFunctionPlot.plotFunction("integral", integral, 0, 100, 800)); BisectionSolver solver = new BisectionSolver(1e-3); -UnivariateFunction integralShifted = {u -> +UnivariateFunction integralShifted = { u -> def integr = integral.value(u); - return integr/(1-integr) - ratio; + return integr / (1 - integr) - ratio; } double integralBorder = solver.solve(400, integralShifted, 10d, 20d); @@ -104,6 +104,6 @@ println "The integral border is $integralBorder"; double newBorder = 14.43 double integralValue = integral.value(newBorder); -double err = Math.abs(integralValue/(1-integralValue)/ratio - 1d) - +double err = Math.abs(integralValue / (1 - integralValue) / ratio - 1d) + println "The relative error ic case of using $newBorder instead of real one is $err"; \ No newline at end of file diff --git a/numass-main/src/main/groovy/inr/numass/scripts/LossTailCalculation.groovy b/numass-main/src/main/groovy/inr/numass/scripts/LossTailCalculation.groovy index 766166ce..cc018f8a 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/LossTailCalculation.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/LossTailCalculation.groovy @@ -8,7 +8,6 @@ package inr.numass.scripts import hep.dataforge.maths.integration.GaussRuleIntegrator import hep.dataforge.maths.integration.UnivariateIntegrator -import inr.numass.models.LossCalculator import org.apache.commons.math3.analysis.UnivariateFunction UnivariateIntegrator integrator = new GaussRuleIntegrator(400); diff --git a/numass-main/src/main/groovy/inr/numass/scripts/PrintLossFunctions.groovy b/numass-main/src/main/groovy/inr/numass/scripts/PrintLossFunctions.groovy index a68b035e..d91c4460 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/PrintLossFunctions.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/PrintLossFunctions.groovy @@ -6,9 +6,10 @@ package inr.numass.scripts -import inr.numass.models.LossCalculator +import inr.numass.models.misc.LossCalculator -LossCalculator loss = LossCalculator.instance() + +LossCalculator loss = LossCalculator.INSTANCE def X = 0.36 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 030f6a58..5c835369 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy @@ -15,12 +15,11 @@ */ package inr.numass.scripts -import hep.dataforge.stat.fit.FitManager -import hep.dataforge.stat.fit.FitState -import hep.dataforge.stat.fit.MINUITPlugin -import hep.dataforge.stat.fit.ParamSet +import hep.dataforge.context.Global +import hep.dataforge.io.FittingIOUtils +import hep.dataforge.meta.Meta +import hep.dataforge.stat.fit.* import hep.dataforge.stat.models.XYModel -import hep.dataforge.tables.ListTable import inr.numass.data.SpectrumAdapter import inr.numass.data.SpectrumGenerator import inr.numass.models.BetaSpectrum @@ -30,7 +29,6 @@ import inr.numass.models.ResolutionFunction import inr.numass.utils.DataModelUtils import org.apache.commons.math3.analysis.BivariateFunction -import static hep.dataforge.context.Global.out import static java.util.Locale.setDefault /** @@ -41,7 +39,8 @@ import static java.util.Locale.setDefault setDefault(Locale.US); new MINUITPlugin().startGlobal(); -FitManager fm = new FitManager(); +FitManager fm = Global.INSTANCE.get(FitManager) + BivariateFunction resolution = new ResolutionFunction(8.3e-5); @@ -49,16 +48,16 @@ ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), resolution, 13490 beta.setCaching(false); NBkgSpectrum spectrum = new NBkgSpectrum(beta); -XYModel model = new XYModel("tritium", spectrum, new SpectrumAdapter()); +XYModel model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), spectrum); ParamSet allPars = new ParamSet(); allPars.setPar("N", 6e5, 10, 0, Double.POSITIVE_INFINITY); -allPars.setPar("bkg", 2d, 0.1 ); +allPars.setPar("bkg", 2d, 0.1); -allPars.setPar("E0", 18575.0, 0.05 ); +allPars.setPar("E0", 18575.0, 0.05); allPars.setPar("mnu2", 0, 1); @@ -74,27 +73,27 @@ allPars.setPar("trap", 0, 0.01, 0d, Double.POSITIVE_INFINITY); SpectrumGenerator generator = new SpectrumGenerator(model, allPars, 12316); -ListTable data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(14000d, 18200, 1e6, 60)); +def data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(14000d, 18200, 1e6, 60)); // data = data.filter("X", Value.of(15510.0), Value.of(18610.0)); allPars.setParValue("U2", 0); FitState state = new FitState(data, model, allPars); //new PlotFitResultAction(Global.instance(), null).runOne(state); - + //double delta = 4e-6; //resolution.setTailFunction{double E, double U -> // 1-delta*(E-U); //} -resolution.setTailFunction(ResolutionFunction.getRealTail()) - +//resolution.setTailFunction(ResolutionFunction.getRealTail()) + //PlotFrame frame = JFreeChartFrame.drawFrame("Transmission function", null); //frame.add(new PlottableFunction("transmission",null, {U -> resolution.value(18500,U)},13500,18505,500)); -FitState res = fm.runTask(state, "QOW", FitTask.TASK_RUN, "N", "bkg", "E0", "U2", "trap"); +def res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "bkg", "E0", "U2", "trap"); - -res.print(out()); +res.printState(new PrintWriter(System.out)) +FittingIOUtils.printResiduals(new PrintWriter(System.out), res.optState().get()) diff --git a/numass-main/src/main/java/inr/numass/models/EmpiricalLossSpectrum.java b/numass-main/src/main/java/inr/numass/models/EmpiricalLossSpectrum.java index 6ffe618d..326f2f3a 100644 --- a/numass-main/src/main/java/inr/numass/models/EmpiricalLossSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/EmpiricalLossSpectrum.java @@ -21,6 +21,7 @@ import hep.dataforge.maths.integration.GaussRuleIntegrator; import hep.dataforge.maths.integration.UnivariateIntegrator; import hep.dataforge.stat.parametric.AbstractParametricFunction; import hep.dataforge.values.Values; +import inr.numass.models.misc.LossCalculator; import org.apache.commons.math3.analysis.BivariateFunction; import org.apache.commons.math3.analysis.UnivariateFunction; @@ -59,11 +60,9 @@ public class EmpiricalLossSpectrum extends AbstractParametricFunction { final double shift = set.getDouble("shift"); //FIXME тут толщины усреднены по длине источника, а надо брать чистого Пуассона - final List probs = LossCalculator.instance().getGunLossProbabilities(X); + final List probs = LossCalculator.INSTANCE.getGunLossProbabilities(X); final double noLossProb = probs.get(0); - final BivariateFunction lossFunction = (Ei, Ef) -> { - return LossCalculator.instance().getLossValue(probs, Ei, Ef); - }; + final BivariateFunction lossFunction = (Ei, Ef) -> LossCalculator.INSTANCE.getLossValue(probs, Ei, Ef); UnivariateFunction integrand = (double x) -> transmission.value(x) * lossFunction.value(x, U - shift); return noLossProb * transmission.value(U - shift) + integrator.integrate(U, eMax, integrand); } diff --git a/numass-main/src/main/java/inr/numass/models/LossCalculator.java b/numass-main/src/main/java/inr/numass/models/LossCalculator.java deleted file mode 100644 index 174b2061..00000000 --- a/numass-main/src/main/java/inr/numass/models/LossCalculator.java +++ /dev/null @@ -1,442 +0,0 @@ -/* - * Copyright 2015 Alexander Nozik. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package inr.numass.models; - -import hep.dataforge.maths.functions.FunctionCaching; -import hep.dataforge.maths.integration.GaussRuleIntegrator; -import hep.dataforge.maths.integration.UnivariateIntegrator; -import hep.dataforge.plots.PlotFrame; -import hep.dataforge.plots.XYFunctionPlot; -import hep.dataforge.utils.Misc; -import hep.dataforge.values.Values; -import org.apache.commons.math3.analysis.BivariateFunction; -import org.apache.commons.math3.analysis.UnivariateFunction; -import org.apache.commons.math3.exception.OutOfRangeException; -import org.slf4j.LoggerFactory; - -import java.util.ArrayList; -import java.util.HashMap; -import java.util.List; -import java.util.Map; - -import static java.lang.Math.exp; - -/** - * Вычисление произвольного порядка функции рассеяния. Не учитывается - * зависимость сечения от энергии электрона - * - * @author Darksnake - */ -public class LossCalculator { - - /** - * порог по вероятности, до которого вычисляются компоненты функции потерь - */ - public final static double SCATTERING_PROBABILITY_THRESHOLD = 1e-3; - - private static final LossCalculator instance = new LossCalculator(); - private static final UnivariateIntegrator integrator = new GaussRuleIntegrator(100); - private static final Map> lossProbCache = Misc.getLRUCache(100); - private final Map cache = new HashMap<>(); - - - - - private LossCalculator() { - cache.put(1, getSingleScatterFunction()); -// immutable.put(2, getDoubleScatterFunction()); - } - - public static UnivariateFunction getSingleScatterFunction() { - final double A1 = 0.204; - final double A2 = 0.0556; - final double b = 14.0; - final double pos1 = 12.6; - final double pos2 = 14.3; - final double w1 = 1.85; - final double w2 = 12.5; - - return (double eps) -> { - if (eps <= 0) { - return 0; - } - if (eps <= b) { - double z = eps - pos1; - return A1 * exp(-2 * z * z / w1 / w1); - } else { - double z = 4 * (eps - pos2) * (eps - pos2); - return A2 / (1 + z / w2 / w2); - } - }; - } - - /** - * A generic loss function for numass experiment in "Lobashev" - * parameterization - * - * @param exPos - * @param ionPos - * @param exW - * @param ionW - * @param exIonRatio - * @return - */ - public static UnivariateFunction getSingleScatterFunction( - final double exPos, - final double ionPos, - final double exW, - final double ionW, - final double exIonRatio) { - UnivariateFunction func = (double eps) -> { - if (eps <= 0) { - return 0; - } - double z1 = eps - exPos; - double ex = exIonRatio * exp(-2 * z1 * z1 / exW / exW); - - double z = 4 * (eps - ionPos) * (eps - ionPos); - double ion = 1 / (1 + z / ionW / ionW); - - double res; - if (eps < exPos) { - res = ex; - } else { - res = Math.max(ex, ion); - } - - return res; - }; - - double cutoff = 25d; - //caclulating lorentz integral analythically - double tailNorm = (Math.atan((ionPos - cutoff) * 2d / ionW) + 0.5 * Math.PI) * ionW / 2d; - final double norm = integrator.integrate(0d, cutoff, func) + tailNorm; - return (e) -> func.value(e) / norm; - } - - public static UnivariateFunction getSingleScatterFunction(Values set) { - - final double exPos = set.getDouble("exPos"); - final double ionPos = set.getDouble("ionPos"); - final double exW = set.getDouble("exW"); - final double ionW = set.getDouble("ionW"); - final double exIonRatio = set.getDouble("exIonRatio"); - - return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio); - } - - public static BivariateFunction getTrapFunction() { - return (double Ei, double Ef) -> { - double eps = Ei - Ef; - if (eps > 10) { - return 1.86e-04 * exp(-eps / 25.0) + 5.5e-05; - } else { - return 0; - } - }; - } - - /** - * Синглетон, так как кэши функций потреь можно вычислять один раз - * - * @return - */ - public static LossCalculator instance() { - return instance; - } - - public static void plotScatter(PlotFrame frame, Values set) { - //"X", "shift", "exPos", "ionPos", "exW", "ionW", "exIonRatio" - -// JFreeChartFrame frame = JFreeChartFrame.drawFrame("Differential scattering crosssection", null); - double X = set.getDouble("X"); - - final double exPos = set.getDouble("exPos"); - - final double ionPos = set.getDouble("ionPos"); - - final double exW = set.getDouble("exW"); - - final double ionW = set.getDouble("ionW"); - - final double exIonRatio = set.getDouble("exIonRatio"); - - UnivariateFunction scatterFunction = getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio); - - if (set.getNames().contains("X")) { - final LossCalculator loss = LossCalculator.instance; - final List probs = loss.getGunLossProbabilities(set.getDouble("X")); - UnivariateFunction single = (double e) -> probs.get(1) * scatterFunction.value(e); - frame.add(XYFunctionPlot.Companion.plot("Single scattering", 0, 100, 1000, single::value)); - - for (int i = 2; i < probs.size(); i++) { - final int j = i; - UnivariateFunction scatter = (double e) -> probs.get(j) * loss.getLossValue(j, e, 0d); - frame.add(XYFunctionPlot.Companion.plot(j + " scattering", 0, 100, 1000, scatter::value)); - } - - UnivariateFunction total = (eps) -> { - if (probs.size() == 1) { - return 0; - } - double sum = probs.get(1) * scatterFunction.value(eps); - for (int i = 2; i < probs.size(); i++) { - sum += probs.get(i) * loss.getLossValue(i, eps, 0); - } - return sum; - }; - - frame.add(XYFunctionPlot.Companion.plot("Total loss", 0, 100, 1000, total::value)); - - } else { - - frame.add(XYFunctionPlot.Companion.plot("Differential crosssection", 0, 100, 2000, scatterFunction::value)); - } - - } - - public List getGunLossProbabilities(double X) { - List res = new ArrayList<>(); - double prob; - if (X > 0) { - prob = Math.exp(-X); - } else { - // если x ==0, то выживает только нулевой член, первый равен 1 - res.add(1d); - return res; - } - res.add(prob); - - int n = 0; - while (prob > SCATTERING_PROBABILITY_THRESHOLD) { - /* - * prob(n) = prob(n-1)*X/n; - */ - n++; - prob *= X / n; - res.add(prob); - } - - return res; - } - - public double getGunZeroLossProb(double X) { - return Math.exp(-X); - } - - /** - * Ленивое рекурсивное вычисление функции потерь через предыдущие - * - * @param order - * @return - */ - private UnivariateFunction getLoss(int order) { - if (order <= 0) { - throw new IllegalArgumentException(); - } - if (cache.containsKey(order)) { - return cache.get(order); - } else { - synchronized (this) { - cache.computeIfAbsent(order, (i) -> { - LoggerFactory.getLogger(getClass()) - .debug("Scatter immutable of order {} not found. Updating", i); - return getNextLoss(getMargin(i), getLoss(i - 1)); - }); - return cache.get(order); - } - } - } - - public BivariateFunction getLossFunction(int order) { - assert order > 0; - return (double Ei, double Ef) -> getLossValue(order, Ei, Ef); - } - - public List getLossProbDerivs(double X) { - List res = new ArrayList<>(); - List probs = getLossProbabilities(X); - - double delta = Math.exp(-X); - res.add((delta - probs.get(0)) / X); - for (int i = 1; i < probs.size(); i++) { - delta *= X / i; - res.add((delta - probs.get(i)) / X); - } - - return res; - } - - /** - * рекурсивно вычисляем все вероятности, котрорые выше порога - *

- * дисер, стр.48 - *

- * @param X - * @return - */ - public List getLossProbabilities(double X) { - return lossProbCache.computeIfAbsent(X, x -> { - List res = new ArrayList<>(); - double prob; - if (x > 0) { - prob = 1 / x * (1 - Math.exp(-x)); - } else { - // если x ==0, то выживает только нулевой член, первый равен нулю - res.add(1d); - return res; - } - res.add(prob); - - while (prob > SCATTERING_PROBABILITY_THRESHOLD) { - /* - * prob(n) = prob(n-1)-1/n! * X^n * exp(-X); - */ - double delta = Math.exp(-x); - for (int i = 1; i < res.size() + 1; i++) { - delta *= x / i; - } - prob -= delta / x; - res.add(prob); - } - - return res; - }); - } - - public double getLossProbability(int order, double X) { - if (order == 0) { - if (X > 0) { - return 1 / X * (1 - Math.exp(-X)); - } else { - return 1d; - } - } - List probs = getLossProbabilities(X); - if (order >= probs.size()) { - return 0; - } else { - return probs.get(order); - } - } - - public double getLossValue(int order, double Ei, double Ef) { - if (Ei - Ef < 5d) { - return 0; - } else if (Ei - Ef >= getMargin(order)) { - return 0; - } else { - return getLoss(order).value(Ei - Ef); - } - } - - /** - * функция потерь с произвольными вероятностями рассеяния - * - * @param probs - * @param Ei - * @param Ef - * @return - */ - public double getLossValue(List probs, double Ei, double Ef) { - double sum = 0; - for (int i = 1; i < probs.size(); i++) { - sum += probs.get(i) * getLossValue(i, Ei, Ef); - } - return sum; - } - - /** - * граница интегрирования - * - * @param order - * @return - */ - private double getMargin(int order) { - return 80 + order * 50d; - } - - /** - * генерирует кэшированную функцию свертки loss со спектром однократных - * потерь - * - * @param loss - * @return - */ - private UnivariateFunction getNextLoss(double margin, UnivariateFunction loss) { - UnivariateFunction res = (final double x) -> { - UnivariateFunction integrand = (double y) -> { - try { - return loss.value(x - y) * getSingleScatterFunction().value(y); - } catch (OutOfRangeException ex) { - return 0; - } - }; - return integrator.integrate(5d, margin, integrand); - }; - - return FunctionCaching.cacheUnivariateFunction(0, margin, 200, res); - - } - - public BivariateFunction getTotalLossBivariateFunction(final double X) { - return (double Ei, double Ef) -> getTotalLossValue(X, Ei, Ef); - } - - /** - * Значение полной производной функции потерь с учетом всех неисчезающих - * порядков - * - * @param X - * @param Ei - * @param Ef - * @return - */ - public double getTotalLossDeriv(double X, double Ei, double Ef) { - List probs = getLossProbDerivs(X); - - double sum = 0; - for (int i = 1; i < probs.size(); i++) { - sum += probs.get(i) * getLossValue(i, Ei, Ef); - } - return sum; - } - - public BivariateFunction getTotalLossDerivBivariateFunction(final double X) { - return (double Ei, double Ef) -> getTotalLossDeriv(X, Ei, Ef); - } - - /** - * Значение полной функции потерь с учетом всех неисчезающих порядков - * - * @param X - * @param Ei - * @param Ef - * @return - */ - public double getTotalLossValue(double X, double Ei, double Ef) { - if (X == 0) { - return 0; - } - List probs = getLossProbabilities(X); - - double sum = 0; - for (int i = 1; i < probs.size(); i++) { - sum += probs.get(i) * getLossValue(i, Ei, Ef); - } - return sum; - } -} diff --git a/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java b/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java index 5122c514..80d8f821 100644 --- a/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java @@ -20,6 +20,7 @@ import hep.dataforge.stat.parametric.AbstractParametricFunction; import hep.dataforge.stat.parametric.ParametricFunction; import hep.dataforge.values.ValueProvider; import hep.dataforge.values.Values; +import inr.numass.models.misc.LossCalculator; import org.apache.commons.math3.analysis.BivariateFunction; import org.slf4j.LoggerFactory; @@ -35,7 +36,7 @@ import java.util.List; public class ModularSpectrum extends AbstractParametricFunction { private static final String[] list = {"X", "trap"}; - private LossCalculator calculator; + private LossCalculator calculator = LossCalculator.INSTANCE; List cacheList; NamedSpectrumCaching trappingCache; BivariateFunction resolution; @@ -61,7 +62,6 @@ public class ModularSpectrum extends AbstractParametricFunction { this.cacheMin = cacheMin; this.cacheMax = cacheMax; this.resolution = resolution; - this.calculator = LossCalculator.instance(); this.sourceSpectrum = source; setupCache(); } @@ -102,7 +102,7 @@ public class ModularSpectrum extends AbstractParametricFunction { //обновляем кэши для трэппинга и упругого прохождения //Using external trappingCache function if provided - BivariateFunction trapFunc = trappingFunction != null ? trappingFunction : LossCalculator.getTrapFunction(); + BivariateFunction trapFunc = trappingFunction != null ? trappingFunction : calculator.getTrapFunction(); BivariateFunction trapRes = new LossResConvolution(trapFunc, resolution); ParametricFunction elasticSpectrum = new TransmissionConvolution(sourceSpectrum, resolution, sourceSpectrum); diff --git a/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java b/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java index 8f083e66..73e410c1 100644 --- a/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java @@ -21,6 +21,7 @@ import hep.dataforge.stat.parametric.AbstractParametricFunction; import hep.dataforge.stat.parametric.ParametricFunction; import hep.dataforge.values.ValueProvider; import hep.dataforge.values.Values; +import inr.numass.models.misc.LossCalculator; import inr.numass.utils.NumassIntegrator; import org.apache.commons.math3.analysis.BivariateFunction; import org.apache.commons.math3.analysis.UnivariateFunction; @@ -40,7 +41,7 @@ public class VariableLossSpectrum extends AbstractParametricFunction { } public static VariableLossSpectrum withData(final UnivariateFunction transmission, double eMax) { - return new VariableLossSpectrum(new AbstractParametricFunction(new String[0]) { + return new VariableLossSpectrum(new AbstractParametricFunction() { @Override public double derivValue(String parName, double x, Values set) { @@ -82,7 +83,7 @@ public class VariableLossSpectrum extends AbstractParametricFunction { double X = set.getDouble("X"); final double shift = set.getDouble("shift"); - final LossCalculator loss = LossCalculator.instance(); + final LossCalculator loss = LossCalculator.INSTANCE; final List probs = loss.getGunLossProbabilities(X); final double noLossProb = probs.get(0); @@ -126,7 +127,7 @@ public class VariableLossSpectrum extends AbstractParametricFunction { final double exW, final double ionW, final double exIonRatio) { - return LossCalculator.getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio); + return LossCalculator.INSTANCE.getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio); } @Override diff --git a/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt b/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt index 28b42a44..42a6ab78 100644 --- a/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt +++ b/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt @@ -23,7 +23,7 @@ import hep.dataforge.meta.Meta import hep.dataforge.plots.jfreechart.JFreeChartFrame import hep.dataforge.providers.Provides import hep.dataforge.providers.ProvidesNames -import hep.dataforge.stat.models.ModelManager +import hep.dataforge.stat.models.ModelLibrary import hep.dataforge.stat.models.WeightedXYModel import hep.dataforge.stat.models.XYModel import hep.dataforge.tables.Adapters @@ -53,7 +53,7 @@ class NumassPlugin : BasicPlugin() { // StorageManager.buildFrom(context); super.attach(context) //TODO Replace by local providers - loadModels(context[ModelManager::class.java]) + loadModels(context[ModelLibrary::class.java]) loadMath(FunctionLibrary.buildFrom(context)) } @@ -116,9 +116,9 @@ class NumassPlugin : BasicPlugin() { /** * Load all numass model factories * - * @param manager + * @param library */ - private fun loadModels(manager: ModelManager) { + private fun loadModels(library: ModelLibrary) { // manager.addModel("modularbeta", (context, an) -> { // double A = an.getDouble("resolution", 8.3e-5);//8.3e-5 @@ -131,7 +131,7 @@ class NumassPlugin : BasicPlugin() { // return new XYModel(spectrum, getAdapter(an)); // }); - manager.addModel("scatter") { context, meta -> + library.addModel("scatter") { context, meta -> val A = meta.getDouble("resolution", 8.3e-5)//8.3e-5 val from = meta.getDouble("from", 0.0) val to = meta.getDouble("to", 0.0) @@ -148,7 +148,7 @@ class NumassPlugin : BasicPlugin() { XYModel(meta, getAdapter(meta), spectrum) } - manager.addModel("scatter-empiric") { context, meta -> + library.addModel("scatter-empiric") { context, meta -> val eGun = meta.getDouble("eGun", 19005.0) val interpolator = buildInterpolator(context, meta, eGun) @@ -161,7 +161,7 @@ class NumassPlugin : BasicPlugin() { WeightedXYModel(meta, getAdapter(meta), spectrum) { dp -> weightReductionFactor } } - manager.addModel("scatter-empiric-variable") { context, meta -> + library.addModel("scatter-empiric-variable") { context, meta -> val eGun = meta.getDouble("eGun", 19005.0) //builder transmisssion with given data, annotation and smoothing @@ -183,7 +183,7 @@ class NumassPlugin : BasicPlugin() { WeightedXYModel(meta, getAdapter(meta), spectrum) { dp -> weightReductionFactor } } - manager.addModel("scatter-analytic-variable") { context, meta -> + library.addModel("scatter-analytic-variable") { context, meta -> val eGun = meta.getDouble("eGun", 19005.0) val loss = VariableLossSpectrum.withGun(eGun + 5) @@ -200,7 +200,7 @@ class NumassPlugin : BasicPlugin() { XYModel(meta, getAdapter(meta), spectrum) } - manager.addModel("scatter-empiric-experimental") { context, meta -> + library.addModel("scatter-empiric-experimental") { context, meta -> val eGun = meta.getDouble("eGun", 19005.0) //builder transmisssion with given data, annotation and smoothing @@ -217,14 +217,14 @@ class NumassPlugin : BasicPlugin() { WeightedXYModel(meta, getAdapter(meta), spectrum) { dp -> weightReductionFactor } } - manager.addModel("sterile") { context, meta -> + library.addModel("sterile") { context, meta -> val sp = SterileNeutrinoSpectrum(context, meta) val spectrum = NBkgSpectrum(sp) XYModel(meta, getAdapter(meta), spectrum) } - manager.addModel("gun") { context, meta -> + library.addModel("gun") { context, meta -> val gsp = GunSpectrum() val tritiumBackground = meta.getDouble("tritiumBkg", 0.0) 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 17368b70..de8c4882 100644 --- a/numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt +++ b/numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt @@ -3,7 +3,6 @@ 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 @@ -14,8 +13,8 @@ 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); +fun model(name: String,factory: ContextMetaFactory): ModelFactory { + return ModelFactory.build(name, factory); } // spectra operations diff --git a/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt b/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt index fd732819..d399cf48 100644 --- a/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt +++ b/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt @@ -21,6 +21,10 @@ import hep.dataforge.plots.PlotFrame import hep.dataforge.plots.XYFunctionPlot import hep.dataforge.utils.Misc import hep.dataforge.values.Values +import kotlinx.coroutines.experimental.CompletableDeferred +import kotlinx.coroutines.experimental.Deferred +import kotlinx.coroutines.experimental.async +import kotlinx.coroutines.experimental.runBlocking import org.apache.commons.math3.analysis.BivariateFunction import org.apache.commons.math3.analysis.UnivariateFunction import org.apache.commons.math3.exception.OutOfRangeException @@ -35,9 +39,7 @@ import java.util.* * @author Darksnake */ object LossCalculator { - private val cache = HashMap().also { - it[1] = singleScatterFunction - } + private val cache = HashMap>() fun getGunLossProbabilities(X: Double): List { @@ -65,8 +67,23 @@ object LossCalculator { return res } - fun getGunZeroLossProb(X: Double): Double { - return Math.exp(-X) + fun getGunZeroLossProb(x: Double): Double { + return Math.exp(-x) + } + + + private fun getCachedSpectrum(order: Int): Deferred { + return when { + order <= 0 -> error("Non-positive loss cache order") + order == 1 -> CompletableDeferred(singleScatterFunction) + else -> cache.getOrPut(order) { + async { + LoggerFactory.getLogger(javaClass) + .debug("Scatter cache of order {} not found. Updating", order) + getNextLoss(getMargin(order), getCachedSpectrum(order - 1).await()) + } + } + } } /** @@ -75,22 +92,8 @@ object LossCalculator { * @param order * @return */ - private fun getLoss(order: Int): UnivariateFunction? { - if (order <= 0) { - throw IllegalArgumentException() - } - return if (cache.containsKey(order)) { - cache[order] - } else { - synchronized(this) { - cache.getOrPut(order) { - LoggerFactory.getLogger(javaClass) - .debug("Scatter cache of order {} not found. Updating", order) - getNextLoss(getMargin(order), getLoss(order - 1)) - } - return cache[order] - } - } + private fun getLoss(order: Int): UnivariateFunction { + return runBlocking { getCachedSpectrum(order).await() } } fun getLossFunction(order: Int): BivariateFunction { @@ -122,32 +125,30 @@ object LossCalculator { * @return */ fun getLossProbabilities(x: Double): List { - return (lossProbCache).getOrPut(x) { - val res = ArrayList() - var prob: Double - if (x > 0) { - prob = 1 / x * (1 - Math.exp(-x)) - } else { - // если x ==0, то выживает только нулевой член, первый равен нулю - res.add(1.0) - return@getOrPut res - } - res.add(prob) - - while (prob > SCATTERING_PROBABILITY_THRESHOLD) { - /* - * prob(n) = prob(n-1)-1/n! * X^n * exp(-X); - */ - var delta = Math.exp(-x) - for (i in 1 until res.size + 1) { - delta *= x / i - } - prob -= delta / x - res.add(prob) - } - - res + val res = ArrayList() + var prob: Double + if (x > 0) { + prob = 1 / x * (1 - Math.exp(-x)) + } else { + // если x ==0, то выживает только нулевой член, первый равен нулю + res.add(1.0) + return res } + res.add(prob) + + while (prob > SCATTERING_PROBABILITY_THRESHOLD) { + /* + * prob(n) = prob(n-1)-1/n! * X^n * exp(-X); + */ + var delta = Math.exp(-x) + for (i in 1 until res.size + 1) { + delta *= x / i + } + prob -= delta / x + res.add(prob) + } + + return res } fun getLossProbability(order: Int, X: Double): Double { @@ -167,12 +168,10 @@ object LossCalculator { } fun getLossValue(order: Int, Ei: Double, Ef: Double): Double { - return if (Ei - Ef < 5.0) { - 0.0 - } else if (Ei - Ef >= getMargin(order)) { - 0.0 - } else { - getLoss(order)!!.value(Ei - Ef) + return when { + Ei - Ef < 5.0 -> 0.0 + Ei - Ef >= getMargin(order) -> 0.0 + else -> getLoss(order).value(Ei - Ef) } } @@ -261,11 +260,11 @@ object LossCalculator { * @return */ fun getTotalLossValue(x: Double, Ei: Double, Ef: Double): Double { - if (x == 0.0) { - return 0.0 + return if (x == 0.0) { + 0.0 } else { val probs = getLossProbabilities(x) - return (0 until probs.size).sumByDouble { i -> + (1 until probs.size).sumByDouble { i -> probs[i] * getLossValue(i, Ei, Ef) } } 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 index 590c6427..345b4714 100644 --- a/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt +++ b/numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt @@ -10,7 +10,7 @@ import hep.dataforge.maths.functions.FunctionLibrary import hep.dataforge.meta.Meta import hep.dataforge.stat.parametric.AbstractParametricBiFunction import hep.dataforge.values.Values -import inr.numass.models.LossCalculator +import inr.numass.models.misc.LossCalculator import inr.numass.utils.ExpressionUtils import org.apache.commons.math3.analysis.BivariateFunction import org.slf4j.LoggerFactory @@ -20,7 +20,6 @@ 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() @@ -55,13 +54,13 @@ class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFun } fun p0(eIn: Double, set: Values): Double { - return LossCalculator.instance().getLossProbability(0, getX(eIn, set)) + return LossCalculator.getLossProbability(0, getX(eIn, set)) } 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) + "X" -> LossCalculator.getTotalLossDeriv(getX(eIn, set), eIn, eOut) else -> super.derivValue(parName, eIn, eOut, set) } } @@ -74,7 +73,7 @@ class NumassTransmission(context: Context, meta: Meta) : AbstractParametricBiFun //calculate X taking into account its energy dependence val X = getX(eIn, set) // loss part - val loss = calculator.getTotalLossValue(X, eIn, eOut) + val loss = LossCalculator.getTotalLossValue(X, eIn, eOut) // double loss; // // if(eIn-eOut >= 300){ diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt b/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt index 425268cc..5be1bac5 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt @@ -31,10 +31,11 @@ import hep.dataforge.stat.fit.FitStage import hep.dataforge.stat.fit.FitState import hep.dataforge.stat.fit.ParamSet import hep.dataforge.stat.models.XYModel -import hep.dataforge.tables.Adapters -import hep.dataforge.tables.ListTable +import hep.dataforge.tables.Adapters.X_AXIS +import hep.dataforge.values.ValueMap import inr.numass.NumassPlugin import inr.numass.data.SpectrumAdapter +import inr.numass.data.SpectrumGenerator import inr.numass.models.NBkgSpectrum import inr.numass.models.sterile.SterileNeutrinoSpectrum import kotlinx.coroutines.experimental.launch @@ -98,14 +99,20 @@ fun main(args: Array) { val x = (14000.0..18400.0).step(100.0).toList() - val table = ListTable.Builder(Adapters.getFormat(adapter)).apply { - x.forEach { u -> - row(adapter.buildSpectrumDataPoint(u, t * spectrum.value(u, params).toLong(), t.toDouble())) - } - }.build() - val model = XYModel(Meta.empty(), adapter, spectrum) - val state = FitState(table, model, params) + + val generator = SpectrumGenerator(model, paramsMod, 12316); + val configuration = x.map { ValueMap.ofPairs(X_AXIS to it, "time" to t) } + val data = generator.generateData(configuration); + +// val table = ListTable.Builder(Adapters.getFormat(adapter)).apply { +// x.forEach { u -> +// row(adapter.buildSpectrumDataPoint(u, t * spectrum.value(u, paramsMod).toLong(), t.toDouble())) +// } +// }.build() + + + val state = FitState(data, model, params) val res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "E0","bkg") res.printState(PrintWriter(System.out)) diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/tristan/AnalyzeTristan.kt b/numass-main/src/main/kotlin/inr/numass/scripts/tristan/AnalyzeTristan.kt index 6fc1364c..d3cfa110 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/tristan/AnalyzeTristan.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/tristan/AnalyzeTristan.kt @@ -12,14 +12,14 @@ fun main(args: Array) { val point = ProtoNumassPoint.readFile(file) point.plotAmplitudeSpectrum() - point.blocks.filter { it.channel == 0 }.firstOrNull()?.let { + point.blocks.firstOrNull { it.channel == 0 }?.let { it.plotAmplitudeSpectrum(plotName = "0") { "title" to "pixel 0" "binning" to 50 } } - point.blocks.filter { it.channel == 4 }.firstOrNull()?.let { + point.blocks.firstOrNull { it.channel == 4 }?.let { it.plotAmplitudeSpectrum(plotName = "4") { "title" to "pixel 4" "binning" to 50