From fbcc8eef9404bb4a4a7e6326aaf7c8646f0566cc Mon Sep 17 00:00:00 2001 From: darksnake Date: Tue, 26 Jul 2016 16:17:26 +0300 Subject: [PATCH] [no commit message] --- .../inr/numass/scripts/SterileDemo.groovy | 4 +- .../inr/numass/models/LossCalculator.java | 11 +-- .../models/sterile/NumassTransmission.java | 8 ++- .../sterile/SterileNeutrinoSpectrum.java | 67 +++++++++++++++++-- 4 files changed, 76 insertions(+), 14 deletions(-) diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy index b0ab6962..e7fd3dab 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/SterileDemo.groovy @@ -50,9 +50,9 @@ import hep.dataforge.io.FittingIOUtils setDefault(Locale.US); -//ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d); +ModularSpectrum beta = new ModularSpectrum(new BetaSpectrum(), 8.3e-5, 13990d, 18600d); -ParametricFunction beta = new SterileNeutrinoSpectrum(); +//ParametricFunction beta = new SterileNeutrinoSpectrum(); NBkgSpectrum spectrum = new NBkgSpectrum(beta); XYModel model = new XYModel(spectrum, new SpectrumDataAdapter()); diff --git a/numass-main/src/main/java/inr/numass/models/LossCalculator.java b/numass-main/src/main/java/inr/numass/models/LossCalculator.java index 223e98b7..dc1160bb 100644 --- a/numass-main/src/main/java/inr/numass/models/LossCalculator.java +++ b/numass-main/src/main/java/inr/numass/models/LossCalculator.java @@ -168,12 +168,12 @@ public class LossCalculator { final LossCalculator loss = LossCalculator.instance; final List probs = loss.getGunLossProbabilities(set.getDouble("X")); UnivariateFunction single = (double e) -> probs.get(1) * scatterFunction.value(e); - frame.add(PlottableXYFunction.plotFunction("Single scattering", x->single.value(x), 0, 100, 1000)); + frame.add(PlottableXYFunction.plotFunction("Single scattering", x -> single.value(x), 0, 100, 1000)); 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(PlottableXYFunction.plotFunction(j + " scattering", x->scatter.value(x), 0, 100, 1000)); + frame.add(PlottableXYFunction.plotFunction(j + " scattering", x -> scatter.value(x), 0, 100, 1000)); } UnivariateFunction total = (eps) -> { @@ -187,11 +187,11 @@ public class LossCalculator { return sum; }; - frame.add(PlottableXYFunction.plotFunction("Total loss", x->total.value(x), 0, 100, 1000)); + frame.add(PlottableXYFunction.plotFunction("Total loss", x -> total.value(x), 0, 100, 1000)); } else { - frame.add(PlottableXYFunction.plotFunction("Differential crosssection", x->scatterFunction.value(x), 0, 100, 2000)); + frame.add(PlottableXYFunction.plotFunction("Differential crosssection", x -> scatterFunction.value(x), 0, 100, 2000)); } } @@ -419,6 +419,9 @@ public class LossCalculator { * @return */ public double getTotalLossValue(double X, double Ei, double Ef) { + if (X == 0) { + return 0; + } List probs = getLossProbabilities(X); double sum = 0; diff --git a/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java b/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java index 7cf095a5..6778d128 100644 --- a/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java +++ b/numass-main/src/main/java/inr/numass/models/sterile/NumassTransmission.java @@ -50,9 +50,13 @@ public class NumassTransmission extends AbstractParametricBiFunction { return true; } - private double getX(double eIn, NamedValueSet set) { + public static double getX(double eIn, NamedValueSet set) { //From our article - return Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6; + return set.getDouble("X")*Math.log(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6; + } + + public static double p0(double eIn, NamedValueSet set) { + return LossCalculator.instance().getLossProbability(0, getX(eIn, set)); } @Override diff --git a/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java b/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java index 0dc0a6ba..ca085fed 100644 --- a/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/sterile/SterileNeutrinoSpectrum.java @@ -12,8 +12,11 @@ import hep.dataforge.fitting.parametric.AbstractParametricFunction; import hep.dataforge.fitting.parametric.ParametricBiFunction; import hep.dataforge.meta.Meta; import hep.dataforge.values.NamedValueSet; +import inr.numass.NumassContext; import inr.numass.models.FSS; import java.util.stream.DoubleStream; +import org.apache.commons.math3.analysis.BivariateFunction; +import org.apache.commons.math3.analysis.UnivariateFunction; import org.apache.commons.math3.distribution.EnumeratedRealDistribution; import org.apache.commons.math3.distribution.RealDistribution; import org.apache.commons.math3.random.JDKRandomGenerator; @@ -29,9 +32,9 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction { private static final String[] list = {"X", "trap", "E0", "mnu2", "msterile2", "U2"}; - private double eMax; private final RandomGenerator rnd; private RealDistribution fssDistribution; + private FSS fss; /** * variables:Eo offset,Ein; parameters: "mnu2", "msterile2", "U2" @@ -50,9 +53,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction { public SterileNeutrinoSpectrum(Context context, Meta configuration) { super(list); rnd = new SynchronizedRandomGenerator(new JDKRandomGenerator()); - this.eMax = 18600; if (configuration.hasValue("fssFile")) { - FSS fss = new FSS(context.io().getFile(configuration.getString("fssFile"))); + fss = new FSS(context.io().getFile(configuration.getString("fssFile"))); fssDistribution = new EnumeratedRealDistribution(rnd, fss.getEs(), fss.getPs()); } @@ -82,13 +84,21 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction { @Override public double value(double u, NamedValueSet set) { - return integrate(u, source, transmission, resolution, set); + if (useDirect()) { + return integrateDirect(u, source, transmission, resolution, set); + } else { + return integrate(u, source, transmission, resolution, set); + } } private int numCalls(double u) { return 100000; } + private boolean useDirect() { + return true; + } + @Override public boolean providesDeriv(String name) { return source.providesDeriv(name) && transmission.providesDeriv(name) && resolution.providesDeriv(name); @@ -114,7 +124,13 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction { ParametricBiFunction transmissionFunction, ParametricBiFunction resolutionFunction, NamedValueSet set) { + int num = numCalls(u); + double eMax = set.getDouble("E0") + 5d; + if (u > eMax) { + return 0; + } + double sum = DoubleStream.generate(() -> { // generate final state double fs; @@ -131,8 +147,8 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction { double res = sourceFunction.value(fs, eIn, set) * transmissionFunction.value(eIn, eOut, set) * resolutionFunction.value(eOut, u, set); - - if(Double.isNaN(res)){ + + if (Double.isNaN(res)) { throw new Error(); } return res; @@ -141,4 +157,43 @@ public class SterileNeutrinoSpectrum extends AbstractParametricFunction { return Math.pow(eMax - u, 2d) / 2d * sum / num; } + private double integrateDirect( + double u, + ParametricBiFunction sourceFunction, + ParametricBiFunction transmissionFunction, + ParametricBiFunction resolutionFunction, + NamedValueSet 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), u, set); + } + return res; + }; + } else { + fsSource = (eIn) -> sourceFunction.value(0, eIn, set); + } + + BivariateFunction transRes = (eIn, usp) -> { + UnivariateFunction integrand = (eOut) -> transmissionFunction.value(eIn, eOut, set) * resolutionFunction.value(eOut, usp, set); + return NumassContext.defaultIntegrator.integrate(integrand, usp, eIn); + }; + + UnivariateFunction integrand = (eIn) -> { + double p0 = NumassTransmission.p0(eIn, set); + return fsSource.value(eIn) * (p0 * resolutionFunction.value(eIn, u, set) + transRes.value(eIn, u)); + }; + + return NumassContext.defaultIntegrator.integrate(integrand, u, eMax); + } + }