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

Updating models

Browse Source
This commit is contained in:
Alexander Nozik 2017-11-26 22:15:16 +03:00
parent a1187722da
commit b6af0f612c
23 changed files with 1024 additions and 887 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
build.gradle
View File

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

13
numass-control/vac/src/main/kotlin/inr/numass/control/readvac/ConsoleVac.kt
View File

@ -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<Double>
.getConstructor(String::class.java).newInstance(port) as Sensor<*>
try {
sensor.init()
while (true) {

2
numass-main/src/main/groovy/inr/numass/scripts/ShowTransmission.groovy
View File

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

6
numass-main/src/main/groovy/inr/numass/scripts/Simulate.groovy
View File

@ -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,7 +86,7 @@ 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);

8
numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy
View File

@ -40,11 +40,11 @@ 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());
@ -55,10 +55,10 @@ 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);

69
numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy Normal file
View File

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

4
numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowFitter.groovy
View File

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

197
numass-main/src/main/java/inr/numass/models/GaussResolution.java
View File

@ -15,33 +15,23 @@
*/
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 - расстояние в сигмах от среднего значения, на котором
* функция считается обращающейся в ноль
*/
@ -50,100 +40,82 @@ public class GaussResolution extends AbstractParametricFunction implements Trans
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);
}
}
}

219
numass-main/src/main/java/inr/numass/models/sterile/NumassBeta.java
View File

@ -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.
* <p>
* dissertation p.33
* </p>
*
* @author <a href="mailto:altavir@gmail.com">Alexander Nozik</a>
*/
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);
}
}

99
numass-main/src/main/java/inr/numass/models/sterile/NumassResolution.java
View File

@ -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 <a href="mailto:altavir@gmail.com">Alexander Nozik</a>
*/
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<String, Object> 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));
}
}
}

114
numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java
View File

@ -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 <a href="mailto:altavir@gmail.com">Alexander Nozik</a>
*/
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<Double, UnivariateFunction> 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<String, Object> 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;
}
}

222
numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java
View File

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

99
numass-main/src/main/java/inr/numass/models/sterile/TestModels.java
View File

@ -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();
/*
<model modelName="sterile" fssFile="FS.txt">
<resolution width = "1.22e-4" tailAlpha="1e-2"/>
<transmission trapping="numass.trap2016"/>
</model>
*/
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);
}
}

4
numass-main/src/main/java/inr/numass/utils/DataModelUtils.java
View File

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

7
numass-main/src/main/java/inr/numass/utils/NumassIntegrator.java
View File

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

24
numass-main/src/main/kotlin/inr/numass/NumassUtils.kt
View File

@ -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<InputStream> { 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<NumassSet>) {
/**
* 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)

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

@ -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<Model>): 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)
}
}
}

230
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassBeta.kt Normal file
View File

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

91
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassResolution.kt Normal file
View File

@ -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<String, Any>()
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<String>() //leaving
}
}

105
numass-main/src/main/kotlin/inr/numass/models/sterile/NumassTransmission.kt Normal file
View File

@ -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<Double, UnivariateFunction>()
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<String, Any>()
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
}
}

164
numass-main/src/main/kotlin/inr/numass/models/sterile/SterileNeutrinoSpectrum.kt Normal file
View File

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

104
numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt Normal file
View File

@ -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<String>) {
val context = Numass.buildContext()
/*
<model modelName="sterile" fssFile="FS.txt">
<resolution width = "1.22e-4" tailAlpha="1e-2"/>
<transmission trapping="numass.trap2016"/>
</model>
*/
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)
}
}

4
numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt
View File

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