numass model update

This commit is contained in:
darksnake 2017-08-15 14:38:17 +03:00
parent 21109003b8
commit 7655237a2c
20 changed files with 104 additions and 76 deletions

View File

@ -43,7 +43,7 @@ UnivariateIntegrator integrator = NumassContext.defaultIntegrator;
double border = 13.6; double border = 13.6;
UnivariateFunction ratioFunction = {e->integrator.integrate(scatterFunction, 0 , e) / integrator.integrate(scatterFunction, e, 100)} UnivariateFunction ratioFunction = {e->integrator.integrate(0, e, scatterFunction) / integrator.integrate(e, 100, scatterFunction)}
double ratio = ratioFunction.value(border); double ratio = ratioFunction.value(border);
println "The true excitation to ionization ratio with border energy $border is $ratio"; println "The true excitation to ionization ratio with border energy $border is $ratio";
@ -83,7 +83,7 @@ UnivariateFunction integral = {double u ->
return 0d; return 0d;
} else { } 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(integrand, 0d, u) return integrator.integrate(0d, u, integrand)
} }
} }

View File

@ -6,9 +6,8 @@
package inr.numass.scripts package inr.numass.scripts
import hep.dataforge.maths.integration.GaussRuleIntegrator; import hep.dataforge.maths.integration.GaussRuleIntegrator
import hep.dataforge.maths.integration.UnivariateIntegrator; import hep.dataforge.maths.integration.UnivariateIntegrator
import inr.numass.models.LossCalculator;
import org.apache.commons.math3.analysis.UnivariateFunction import org.apache.commons.math3.analysis.UnivariateFunction
UnivariateIntegrator integrator = new GaussRuleIntegrator(400); UnivariateIntegrator integrator = new GaussRuleIntegrator(400);
@ -42,6 +41,6 @@ UnivariateFunction func = {double eps ->
//caclulating lorentz integral analythically //caclulating lorentz integral analythically
double tailNorm = (Math.atan((ionPos - cutoff) * 2d / ionW) + 0.5 * Math.PI) * ionW / 2d; double tailNorm = (Math.atan((ionPos - cutoff) * 2d / ionW) + 0.5 * Math.PI) * ionW / 2d;
final double norm = integrator.integrate(func, 0d, cutoff) + tailNorm; final double norm = integrator.integrate(0d, cutoff, func) + tailNorm;
println 1/norm; println 1/norm;

View File

@ -24,11 +24,11 @@ def cutoff = 20d
UnivariateFunction loss = LossCalculator.getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio); UnivariateFunction loss = LossCalculator.getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio);
println integrator.integrate(loss,0,600); println integrator.integrate(0, 600, loss);
println integrator.integrate(loss,0, cutoff); println integrator.integrate(0, cutoff, loss);
println integrator.integrate(loss,cutoff,600d); println integrator.integrate(cutoff, 600d, loss);
println (integrator.integrate(loss,0,cutoff) + integrator.integrate(loss,cutoff,3000d)); println (integrator.integrate(0, cutoff, loss) + integrator.integrate(cutoff, 3000d, loss));
//double tailValue = (Math.atan((ionPos-cutoff)*2d/ionW) + 0.5*Math.PI)*ionW/2; //double tailValue = (Math.atan((ionPos-cutoff)*2d/ionW) + 0.5*Math.PI)*ionW/2;
//println tailValue //println tailValue
//println integrator.integrate(loss,0,100); //println integrator.integrate(loss,0,100);

View File

@ -15,17 +15,14 @@
*/ */
package inr.numass.scripts package inr.numass.scripts
import hep.dataforge.io.PrintFunction
import hep.dataforge.maths.integration.UnivariateIntegrator import hep.dataforge.maths.integration.UnivariateIntegrator
import hep.dataforge.plots.data.PlottableXYFunction import hep.dataforge.plots.data.PlottableXYFunction
import hep.dataforge.plots.jfreechart.JFreeChartFrame import hep.dataforge.plots.jfreechart.JFreeChartFrame
import hep.dataforge.stat.fit.ParamSet
import inr.numass.models.ExperimentalVariableLossSpectrum import inr.numass.models.ExperimentalVariableLossSpectrum
import org.apache.commons.math3.analysis.UnivariateFunction import org.apache.commons.math3.analysis.UnivariateFunction
import inr.numass.Numass
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.io.PrintFunction
//double exPos = 12.94 //double exPos = 12.94
//double exW = 1.31 //double exW = 1.31
//double ionPos = 14.13 //double ionPos = 14.13
@ -50,17 +47,17 @@ UnivariateIntegrator integrator = NumassContext.defaultIntegrator
UnivariateFunction exFunc = lsp.excitation(params.getValue("exPos"), params.getValue("exW")); UnivariateFunction exFunc = lsp.excitation(params.getValue("exPos"), params.getValue("exW"));
frame.add(PlottableXYFunction.plotFunction("ex", exFunc, 0d, 50d, 500)); frame.add(PlottableXYFunction.plotFunction("ex", exFunc, 0d, 50d, 500));
println "excitation norm factor " + integrator.integrate(exFunc, 0, 50) println "excitation norm factor " + integrator.integrate(0, 50, exFunc)
UnivariateFunction ionFunc = lsp.ionization(params.getValue("ionPos"), params.getValue("ionW")); UnivariateFunction ionFunc = lsp.ionization(params.getValue("ionPos"), params.getValue("ionW"));
frame.add(PlottableXYFunction.plotFunction("ion", ionFunc, 0d, 50d, 500)); frame.add(PlottableXYFunction.plotFunction("ion", ionFunc, 0d, 50d, 500));
println "ionization norm factor " + integrator.integrate(ionFunc, 0, 200) println "ionization norm factor " + integrator.integrate(0, 200, ionFunc)
UnivariateFunction sumFunc = lsp.singleScatterFunction(params); UnivariateFunction sumFunc = lsp.singleScatterFunction(params);
frame.add(PlottableXYFunction.plotFunction("sum", sumFunc, 0d, 50d, 500)); frame.add(PlottableXYFunction.plotFunction("sum", sumFunc, 0d, 50d, 500));
println "sum norm factor " + integrator.integrate(sumFunc, 0, 100) println "sum norm factor " + integrator.integrate(0, 100, sumFunc)
PrintFunction.printFunctionSimple(new PrintWriter(System.out), sumFunc, 0d, 50d, 100) PrintFunction.printFunctionSimple(new PrintWriter(System.out), sumFunc, 0d, 50d, 100)

View File

@ -0,0 +1,46 @@
package inr.numass.scripts.temp
import hep.dataforge.context.Context
import hep.dataforge.context.Global
import hep.dataforge.grind.Grind
import hep.dataforge.grind.GrindShell
import hep.dataforge.grind.helpers.PlotHelper
import hep.dataforge.meta.Meta
import hep.dataforge.plots.fx.FXPlotManager
import hep.dataforge.stat.fit.ParamSet
import hep.dataforge.utils.MetaMorph
import inr.numass.NumassPlugin
import inr.numass.models.sterile.SterileNeutrinoSpectrum
Context ctx = Global.instance()
ctx.pluginManager().load(FXPlotManager)
ctx.pluginManager().load(NumassPlugin.class)
new GrindShell(ctx).eval {
SterileNeutrinoSpectrum sp1 = new SterileNeutrinoSpectrum(context, Meta.empty());
SterileNeutrinoSpectrum sp2 = new SterileNeutrinoSpectrum(context, Grind.buildMeta(useFSS: false));
def params = MetaMorph.morph(ParamSet,
Grind.buildMeta("params") {
N(value: 6e5, err: 1e5, lower: 0)
bkg(value: 2, err: 0.1)
E0(value: 18575, err: 0.1)
mnu2(value: 0, err: 0.01)
msterile2(value: 1000**2, err: 1)
U2(value: 0, err: 1e-3)
X(value: 0, err: 0.01, lower: 0)
trap(value: 0, err: 0.05)
}
)
def xs = (1..400).collect{18000 + it*2}
def sp1Points = xs.collect { sp1.value(it, params) }
def sp2Points = xs.collect { sp2.value(it, params) }
(plots as PlotHelper).plot(xs,sp1Points,"FSS")
(plots as PlotHelper).plot(xs,sp2Points,"noFSS")
}

View File

@ -22,7 +22,7 @@ public class CustomNBkgSpectrum extends NBkgSpectrum {
UnivariateFunction differentialBkgFunction = NumassUtils.tritiumBackgroundFunction(amplitude); UnivariateFunction differentialBkgFunction = NumassUtils.tritiumBackgroundFunction(amplitude);
UnivariateFunction integralBkgFunction = UnivariateFunction integralBkgFunction =
(x) -> NumassIntegrator.getDefaultIntegrator() (x) -> NumassIntegrator.getDefaultIntegrator()
.integrate(differentialBkgFunction, x, 18580d); .integrate(x, 18580d, differentialBkgFunction);
return new CustomNBkgSpectrum(source, integralBkgFunction); return new CustomNBkgSpectrum(source, integralBkgFunction);
} }

View File

@ -65,7 +65,7 @@ public class EmpiricalLossSpectrum extends AbstractParametricFunction {
return LossCalculator.instance().getLossValue(probs, Ei, Ef); return LossCalculator.instance().getLossValue(probs, Ei, Ef);
}; };
UnivariateFunction integrand = (double x) -> transmission.value(x) * lossFunction.value(x, U - shift); UnivariateFunction integrand = (double x) -> transmission.value(x) * lossFunction.value(x, U - shift);
return noLossProb * transmission.value(U - shift) + integrator.integrate(integrand, U, eMax); return noLossProb * transmission.value(U - shift) + integrator.integrate(U, eMax, integrand);
} }
@Override @Override

View File

@ -70,7 +70,7 @@ public class GunSpectrum extends AbstractParametricFunction {
} else if (pos - cutoff * sigma > U * (1 + resA)) { } else if (pos - cutoff * sigma > U * (1 + resA)) {
return 0; return 0;
} else { } else {
return integrator.integrate(integrand, pos - cutoff * sigma, pos + cutoff * sigma); return integrator.integrate(pos - cutoff * sigma, pos + cutoff * sigma, integrand);
} }
} }
@ -148,7 +148,7 @@ public class GunSpectrum extends AbstractParametricFunction {
} else if (pos - cutoff * sigma > U * (1 + resA)) { } else if (pos - cutoff * sigma > U * (1 + resA)) {
return 1; return 1;
} else { } else {
return integrator.integrate(integrand, pos - cutoff * sigma, pos + cutoff * sigma); return integrator.integrate(pos - cutoff * sigma, pos + cutoff * sigma, integrand);
} }
} }

View File

@ -123,7 +123,7 @@ public class LossCalculator {
double cutoff = 25d; double cutoff = 25d;
//caclulating lorentz integral analythically //caclulating lorentz integral analythically
double tailNorm = (Math.atan((ionPos - cutoff) * 2d / ionW) + 0.5 * Math.PI) * ionW / 2d; double tailNorm = (Math.atan((ionPos - cutoff) * 2d / ionW) + 0.5 * Math.PI) * ionW / 2d;
final double norm = integrator.integrate(func, 0d, cutoff) + tailNorm; final double norm = integrator.integrate(0d, cutoff, func) + tailNorm;
return (e) -> func.value(e) / norm; return (e) -> func.value(e) / norm;
} }
@ -385,7 +385,7 @@ public class LossCalculator {
return 0; return 0;
} }
}; };
return integrator.integrate(integrand, 5d, margin); return integrator.integrate(5d, margin, integrand);
}; };
return FunctionCaching.cacheUnivariateFunction(0, margin, 200, res); return FunctionCaching.cacheUnivariateFunction(0, margin, 200, res);

View File

@ -38,7 +38,7 @@ class LossResConvolution implements BivariateFunction {
public double value(final double Ein, final double U) { public double value(final double Ein, final double U) {
UnivariateFunction integrand = (double Eout) -> loss.value(Ein, Eout) * res.value(Eout, U); UnivariateFunction integrand = (double Eout) -> loss.value(Ein, Eout) * res.value(Eout, U);
//Энергия в принципе не может быть больше начальной и меньше напряжения //Энергия в принципе не может быть больше начальной и меньше напряжения
return NumassIntegrator.getDefaultIntegrator().integrate(integrand, U, Ein); return NumassIntegrator.getDefaultIntegrator().integrate(U, Ein, integrand);
} }
} }

View File

@ -58,7 +58,7 @@ class TransmissionConvolution extends AbstractParametricFunction {
} }
return trans.value(E, U) * spectrum.derivValue(parName, E, set); return trans.value(E, U) * spectrum.derivValue(parName, E, set);
}; };
return NumassIntegrator.getDefaultIntegrator().integrate(integrand, Math.max(U, min), max + 1d); return NumassIntegrator.getDefaultIntegrator().integrate(Math.max(U, min), max + 1d, integrand);
} }
@Override @Override
@ -80,6 +80,6 @@ class TransmissionConvolution extends AbstractParametricFunction {
} }
return trans.value(E, U) * spectrum.value(E, set); return trans.value(E, U) * spectrum.value(E, set);
}; };
return NumassIntegrator.getDefaultIntegrator().integrate(integrand, Math.max(U, min), max + 1d); return NumassIntegrator.getDefaultIntegrator().integrate(Math.max(U, min), max + 1d, integrand);
} }
} }

View File

@ -106,7 +106,7 @@ public class VariableLossSpectrum extends AbstractParametricFunction {
} else { } else {
integrator = NumassIntegrator.getDefaultIntegrator(); integrator = NumassIntegrator.getDefaultIntegrator();
} }
return noLossProb * transmission.value(U - shift, set) + integrator.integrate(integrand, U, eMax); return noLossProb * transmission.value(U - shift, set) + integrator.integrate(U, eMax, integrand);
} }
public UnivariateFunction singleScatterFunction(ValueProvider set) { public UnivariateFunction singleScatterFunction(ValueProvider set) {

View File

@ -132,7 +132,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
double eMax = set.getDouble("E0") + 5d; double eMax = set.getDouble("E0") + 5d;
if (u > eMax) { if (u >= eMax) {
return 0; return 0;
} }
@ -163,7 +163,7 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
integrator = NumassIntegrator.getHighDensityIntegrator(); integrator = NumassIntegrator.getHighDensityIntegrator();
} }
return integrator.integrate(eIn -> fsSource.value(eIn) * transResFunction.value(eIn, u, set), u, eMax); return integrator.integrate(u, eMax, eIn -> fsSource.value(eIn) * transResFunction.value(eIn, u, set));
} }
private class TransRes extends AbstractParametricBiFunction { private class TransRes extends AbstractParametricBiFunction {
@ -201,13 +201,13 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction {
UnivariateFunction integrand = (eOut) -> transFunc.value(eIn, eOut, set) * resolution.value(eOut, u, set); UnivariateFunction integrand = (eOut) -> transFunc.value(eIn, eOut, set) * resolution.value(eOut, u, set);
double border = u + 30; double border = u + 30;
double firstPart = NumassIntegrator.getFastInterator().integrate(integrand, u, Math.min(eIn, border)); double firstPart = NumassIntegrator.getFastInterator().integrate(u, Math.min(eIn, border), integrand);
double secondPart; double secondPart;
if (eIn > border) { if (eIn > border) {
if (fast) { if (fast) {
secondPart = NumassIntegrator.getDefaultIntegrator().integrate(integrand, border, eIn); secondPart = NumassIntegrator.getDefaultIntegrator().integrate(border, eIn, integrand);
} else { } else {
secondPart = NumassIntegrator.getHighDensityIntegrator().integrate(integrand, border, eIn); secondPart = NumassIntegrator.getHighDensityIntegrator().integrate(border, eIn, integrand);
} }
} else { } else {
secondPart = 0; secondPart = 0;

View File

@ -12,6 +12,7 @@ import hep.dataforge.data.DataNode;
import hep.dataforge.data.DataSet; import hep.dataforge.data.DataSet;
import hep.dataforge.description.TypedActionDef; import hep.dataforge.description.TypedActionDef;
import hep.dataforge.meta.Laminate; import hep.dataforge.meta.Laminate;
import hep.dataforge.meta.Meta;
import hep.dataforge.stat.fit.FitResult; import hep.dataforge.stat.fit.FitResult;
import hep.dataforge.stat.fit.ParamSet; import hep.dataforge.stat.fit.ParamSet;
import hep.dataforge.stat.fit.UpperLimitGenerator; import hep.dataforge.stat.fit.UpperLimitGenerator;
@ -41,12 +42,11 @@ public class NumassFitScanSummaryTask extends AbstractTask<Table> {
} }
@Override @Override
protected TaskModel transformModel(TaskModel model) { protected void updateModel(TaskModel.Builder model, Meta meta) {
//Transmit meta as-is model.dependsOn("fitscan", meta, "fitscan");
model.dependsOn("fitscan", model.meta(), "fitscan");
return model;
} }
@Override @Override
public String getName() { public String getName() {
return "scansum"; return "scansum";

View File

@ -78,17 +78,14 @@ public class NumassFitScanTask extends AbstractTask<FitResult> {
} }
@Override @Override
protected TaskModel transformModel(TaskModel model) { protected void updateModel(TaskModel.Builder model, Meta meta) {
//Transmit meta as-is if (meta.hasMeta("filter")) {
MetaBuilder metaBuilder = new MetaBuilder(model.meta()).removeNode("fit").removeNode("scan"); model.dependsOn("filter", meta, "prepare");
if (model.meta().hasMeta("filter")) { } else if (meta.hasMeta("empty")) {
model.dependsOn("filter", metaBuilder.build(), "prepare"); model.dependsOn("subtractEmpty", meta, "prepare");
} else if (model.meta().hasMeta("empty")) {
model.dependsOn("subtractEmpty", metaBuilder.build(), "prepare");
} else { } else {
model.dependsOn("prepare", metaBuilder.build(), "prepare"); model.dependsOn("prepare", meta, "prepare");
} }
return model;
} }
@Override @Override

View File

@ -19,7 +19,6 @@ package inr.numass.tasks;
import hep.dataforge.actions.Action; import hep.dataforge.actions.Action;
import hep.dataforge.data.DataNode; import hep.dataforge.data.DataNode;
import hep.dataforge.meta.Meta; import hep.dataforge.meta.Meta;
import hep.dataforge.meta.MetaBuilder;
import hep.dataforge.stat.fit.FitState; import hep.dataforge.stat.fit.FitState;
import hep.dataforge.tables.Table; import hep.dataforge.tables.Table;
import hep.dataforge.workspace.SingleActionTask; import hep.dataforge.workspace.SingleActionTask;
@ -51,12 +50,7 @@ public class NumassFitSummaryTask extends SingleActionTask<FitState, Table> {
} }
@Override @Override
protected TaskModel transformModel(TaskModel model) { protected void updateModel(TaskModel.Builder model, Meta meta) {
//Transmit meta as-is model.dependsOn("fit", meta, "fit");
MetaBuilder meta = model.meta().getBuilder().removeNode("summary");
model.dependsOn("fit", meta.build(), "fit");
return model;
} }
} }

View File

@ -20,7 +20,6 @@ import hep.dataforge.actions.Action;
import hep.dataforge.actions.ActionUtils; import hep.dataforge.actions.ActionUtils;
import hep.dataforge.data.DataNode; import hep.dataforge.data.DataNode;
import hep.dataforge.meta.Meta; import hep.dataforge.meta.Meta;
import hep.dataforge.meta.MetaBuilder;
import hep.dataforge.plotfit.PlotFitResultAction; import hep.dataforge.plotfit.PlotFitResultAction;
import hep.dataforge.stat.fit.FitAction; import hep.dataforge.stat.fit.FitAction;
import hep.dataforge.stat.fit.FitResult; import hep.dataforge.stat.fit.FitResult;
@ -65,17 +64,15 @@ public class NumassFitTask extends SingleActionTask<Table, FitResult> {
return model.meta().getMeta("fit"); return model.meta().getMeta("fit");
} }
@Override @Override
protected TaskModel transformModel(TaskModel model) { protected void updateModel(TaskModel.Builder model, Meta meta) {
//Transmit meta as-is if (meta.hasMeta("filter")) {
MetaBuilder metaBuilder = new MetaBuilder(model.meta()).removeNode("fit"); model.dependsOn("filter", meta, "prepare");
if (model.meta().hasMeta("filter")) { } else if (meta.hasMeta("empty")) {
model.dependsOn("filter", metaBuilder.build(), "prepare"); model.dependsOn("subtractEmpty", meta, "prepare");
} else if (model.meta().hasMeta("empty")) {
model.dependsOn("subtractEmpty", metaBuilder.build(), "prepare");
} else { } else {
model.dependsOn("prepare", metaBuilder.build(), "prepare"); model.dependsOn("prepare", meta, "prepare");
} }
return model;
} }
} }

View File

@ -84,13 +84,12 @@ public class NumassPrepareTask extends AbstractTask<Table> {
} }
@Override @Override
protected TaskModel transformModel(TaskModel model) { protected void updateModel(TaskModel.Builder model, Meta meta) {
if (model.hasValue("data.from")) { if (meta.hasValue("data.from")) {
model.data(model.getString("data.from.*")); model.data(meta.getString("data.from.*"));
} else { } else {
model.data("*"); model.data("*");
} }
return model;
} }
// private DataSet.Builder<NumassData> readData(Work callback, Context context, URI numassRoot, Meta meta) { // private DataSet.Builder<NumassData> readData(Work callback, Context context, URI numassRoot, Meta meta) {

View File

@ -6,7 +6,7 @@ import hep.dataforge.context.Context;
import hep.dataforge.data.DataNode; import hep.dataforge.data.DataNode;
import hep.dataforge.description.TypedActionDef; import hep.dataforge.description.TypedActionDef;
import hep.dataforge.meta.Laminate; import hep.dataforge.meta.Laminate;
import hep.dataforge.meta.MetaBuilder; import hep.dataforge.meta.Meta;
import hep.dataforge.tables.Table; import hep.dataforge.tables.Table;
import hep.dataforge.tables.TableTransform; import hep.dataforge.tables.TableTransform;
import hep.dataforge.values.Value; import hep.dataforge.values.Value;
@ -34,15 +34,14 @@ public class NumassTableFilterTask extends SingleActionTask<Table, Table> {
return data.getCheckedNode("prepare", Table.class); return data.getCheckedNode("prepare", Table.class);
} }
@Override @Override
protected TaskModel transformModel(TaskModel model) { protected void updateModel(TaskModel.Builder model, Meta meta) {
MetaBuilder metaBuilder = new MetaBuilder(model.meta()).removeNode("filter"); if (meta.hasMeta("empty")) {
if (model.meta().hasMeta("empty")) { model.dependsOn("subtractEmpty", meta, "prepare");
model.dependsOn("subtractEmpty", metaBuilder.build(), "prepare");
} else { } else {
model.dependsOn("prepare", metaBuilder.build(), "prepare"); model.dependsOn("prepare", meta, "prepare");
} }
return model;
} }
@Override @Override

View File

@ -37,7 +37,7 @@ public class TestNeLossParametrisation {
UnivariateFunction oldFunction = LossCalculator.getSingleScatterFunction(); UnivariateFunction oldFunction = LossCalculator.getSingleScatterFunction();
UnivariateFunction newFunction = getSingleScatterFunction(12.86, 16.78, 1.65, 12.38, 4.79); UnivariateFunction newFunction = getSingleScatterFunction(12.86, 16.78, 1.65, 12.38, 4.79);
Double norm = new GaussRuleIntegrator(200).integrate(newFunction, 0d, 100d); Double norm = new GaussRuleIntegrator(200).integrate(0d, 100d, newFunction);
System.out.println(norm); System.out.println(norm);