NPM/numass-framework
4
0
Fork 1
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Fixed fit internals

Browse Source
This commit is contained in:
Alexander Nozik 2018-07-06 22:07:16 +03:00
parent 60bddaa50c
commit 53aad7f6a4
14 changed files with 137 additions and 576 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
numass-main/src/main/groovy/inr/numass/scripts/FindExIonRatio.groovy
View File

@ -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<Double> 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";

1
numass-main/src/main/groovy/inr/numass/scripts/LossTailCalculation.groovy
View File

@ -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);

5
numass-main/src/main/groovy/inr/numass/scripts/PrintLossFunctions.groovy
View File

@ -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

33
numass-main/src/main/groovy/inr/numass/scripts/Systematics.groovy
View File

@ -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())

7
numass-main/src/main/java/inr/numass/models/EmpiricalLossSpectrum.java
View File

@ -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<Double> probs = LossCalculator.instance().getGunLossProbabilities(X);
final List<Double> 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);
}

442
numass-main/src/main/java/inr/numass/models/LossCalculator.java
View File

@ -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<Double, List<Double>> lossProbCache = Misc.getLRUCache(100);
private final Map<Integer, UnivariateFunction> 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<Double> 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<Double> getGunLossProbabilities(double X) {
List<Double> 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<Double> getLossProbDerivs(double X) {
List<Double> res = new ArrayList<>();
List<Double> 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;
}
/**
* рекурсивно вычисляем все вероятности, котрорые выше порога
* <p>
* дисер, стр.48
* </p>
* @param X
* @return
*/
public List<Double> getLossProbabilities(double X) {
return lossProbCache.computeIfAbsent(X, x -> {
List<Double> 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<Double> 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<Double> 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<Double> 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<Double> probs = getLossProbabilities(X);
double sum = 0;
for (int i = 1; i < probs.size(); i++) {
sum += probs.get(i) * getLossValue(i, Ei, Ef);
}
return sum;
}
}

6
numass-main/src/main/java/inr/numass/models/ModularSpectrum.java
View File

@ -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<NamedSpectrumCaching> 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);

7
numass-main/src/main/java/inr/numass/models/VariableLossSpectrum.java
View File

@ -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<Double> 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

22
numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt
View File

@ -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)

5
numass-main/src/main/kotlin/inr/numass/models/NumassModels.kt
View File

@ -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<Model>): ModelFactory {
return ModelFactory.build(name, descriptor, factory);
fun model(name: String,factory: ContextMetaFactory<Model>): ModelFactory {
return ModelFactory.build(name, factory);
}
// spectra operations

109
numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt
View File

@ -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<Int, UnivariateFunction>().also {
it[1] = singleScatterFunction
}
private val cache = HashMap<Int, Deferred<UnivariateFunction>>()
fun getGunLossProbabilities(X: Double): List<Double> {
@ -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<UnivariateFunction> {
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<Double> {
return (lossProbCache).getOrPut(x) {
val res = ArrayList<Double>()
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<Double>()
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)
}
}

9
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt
View File

@ -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<Double, UnivariateFunction>()
@ -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){

25
numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt
View File

@ -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<String>) {
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))

4
numass-main/src/main/kotlin/inr/numass/scripts/tristan/AnalyzeTristan.kt
View File

@ -12,14 +12,14 @@ fun main(args: Array<String>) {
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