diff --git a/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy b/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy index b0f79c7e..80927d7e 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/TritiumTest.groovy @@ -28,7 +28,6 @@ import inr.numass.models.BetaSpectrum import inr.numass.models.ModularSpectrum import inr.numass.models.NBkgSpectrum -import static hep.dataforge.maths.RandomUtils.setSeed import static inr.numass.utils.DataModelUtils.getUniformSpectrumConfiguration PrintWriter out = Global.out(); diff --git a/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy b/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy index af5fd8eb..88af28dd 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/models/TristanModel.groovy @@ -36,26 +36,32 @@ new GrindShell(ctx).eval { def response = new ModGauss(5.0) ParametricFunction spectrum = NumassModelsKt.convolute(beta, response) - def model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), new NBkgSpectrum(spectrum)); + def model = new XYModel(Meta.empty(), new SpectrumAdapter(Meta.empty()), new NBkgSpectrum(spectrum)) ParamSet params = morph(ParamSet, [:], "params") { N(value: 1e+14, err: 30, lower: 0) bkg(value: 5.0, err: 0.1) E0(value: 18575.0, err: 0.1) mnu2(value: 0, err: 0.01) - msterile2(value: 1000**2, err: 1) + msterile2(value: 7000**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(value: 150, err: 5) //shift(value: 1, err: 1e-2) - tailAmp(value: 0.01, err: 1e-2) + //tailAmp(value: 0.01, err: 1e-2) tailW(value: 300, err: 1) } - ph.plotFunction(-2000d, 500d, 400, "actual", "response") { double x -> - response.value(x, params) - } +// double norm = NumassIntegrator.defaultIntegrator.integrate(1000d, 18500d) { +// model.value(it, params) +// } + +// println("The total number of events is $norm") +// +// ph.plotFunction(-2000d, 500d, 400, "actual", "response") { double x -> +// response.value(x, params) +// } SpectrumGenerator generator = new SpectrumGenerator(model, params, 12316); @@ -65,7 +71,7 @@ new GrindShell(ctx).eval { .configure(showLine: true, showSymbol: false, showErrors: false, thickness: 2, connectionType: "spline", color: "red") - Table data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(10000, 19500, 1, 950)); + Table data = generator.generateData(DataModelUtils.getUniformSpectrumConfiguration(7000, 19500, 1, 1000)); //params.setParValue("w", 151) //params.setParValue("tailAmp", 0.011) diff --git a/numass-main/src/main/java/inr/numass/Numass.java b/numass-main/src/main/java/inr/numass/Numass.java index 7a06de38..9451f045 100644 --- a/numass-main/src/main/java/inr/numass/Numass.java +++ b/numass-main/src/main/java/inr/numass/Numass.java @@ -48,7 +48,7 @@ public class Numass { .ln() .text("\t") .content( - MarkupUtils.markupDescriptor(DescriptorUtils.buildDescriptor("method::hep.dataforge.data.DataManager.read")) + MarkupUtils.INSTANCE.markupDescriptor(DescriptorUtils.buildDescriptor("method::hep.dataforge.data.DataManager.read")) ) .ln() .text("***Allowed actions***", "red") @@ -60,7 +60,7 @@ public class Numass { am.getAllActions() .map(name -> am.optAction(name).get()) .map(ActionDescriptor::build).forEach(descriptor -> - builder.text("\t").content(MarkupUtils.markupDescriptor(descriptor)) + builder.text("\t").content(MarkupUtils.INSTANCE.markupDescriptor(descriptor)) ); builder.text("***End of actions list***", "red"); diff --git a/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java b/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java index 8554b912..2ddfcd26 100644 --- a/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java +++ b/numass-main/src/main/java/inr/numass/data/SpectrumGenerator.java @@ -16,6 +16,7 @@ package inr.numass.data; import hep.dataforge.meta.Meta; +import hep.dataforge.stat.RandomKt; import hep.dataforge.stat.fit.ParamSet; import hep.dataforge.stat.models.Generator; import hep.dataforge.stat.models.XYModel; @@ -27,7 +28,6 @@ import org.apache.commons.math3.random.JDKRandomGenerator; import org.apache.commons.math3.random.RandomDataGenerator; import org.apache.commons.math3.random.RandomGenerator; -import static hep.dataforge.maths.RandomUtils.getDefaultRandomGenerator; import static java.lang.Double.isNaN; import static java.lang.Math.sqrt; @@ -62,7 +62,7 @@ public class SpectrumGenerator implements Generator { } public SpectrumGenerator(XYModel source, ParamSet params) { - this(source, params, INSTANCE.getDefaultRandomGenerator()); + this(source, params, RandomKt.getDefaultGenerator()); } @Override 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 fa62d412..a9e0ba0a 100644 --- a/numass-main/src/main/java/inr/numass/models/LossCalculator.java +++ b/numass-main/src/main/java/inr/numass/models/LossCalculator.java @@ -180,12 +180,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(XYFunctionPlot.plotFunction("Single scattering", single::value, 0, 100, 1000)); + frame.add(XYFunctionPlot.plotFunction("Single scattering", 0, 100, 1000, single::value)); for (int i = 2; i < probs.size(); i++) { final int j = i; UnivariateFunction scatter = (double e) -> probs.get(j) * loss.getLossValue(j, e, 0d); - frame.add(XYFunctionPlot.plotFunction(j + " scattering", scatter::value, 0, 100, 1000)); + frame.add(XYFunctionPlot.plotFunction(j + " scattering", 0, 100, 1000, scatter::value)); } UnivariateFunction total = (eps) -> { @@ -199,11 +199,11 @@ public class LossCalculator { return sum; }; - frame.add(XYFunctionPlot.plotFunction("Total loss", total::value, 0, 100, 1000)); + frame.add(XYFunctionPlot.plotFunction("Total loss", 0, 100, 1000, total::value)); } else { - frame.add(XYFunctionPlot.plotFunction("Differential crosssection", scatterFunction::value, 0, 100, 2000)); + frame.add(XYFunctionPlot.plotFunction("Differential crosssection", 0, 100, 2000, scatterFunction::value)); } } diff --git a/numass-main/src/main/kotlin/inr/numass/data/Generator.kt b/numass-main/src/main/kotlin/inr/numass/data/Generator.kt index 16b9e3cb..fd69a995 100644 --- a/numass-main/src/main/kotlin/inr/numass/data/Generator.kt +++ b/numass-main/src/main/kotlin/inr/numass/data/Generator.kt @@ -3,13 +3,13 @@ package inr.numass.data import hep.dataforge.maths.chain.Chain import hep.dataforge.maths.chain.MarkovChain import hep.dataforge.maths.chain.StatefulChain +import hep.dataforge.stat.defaultGenerator import inr.numass.data.api.NumassBlock import inr.numass.data.api.NumassEvent import inr.numass.data.api.SimpleBlock import kotlinx.coroutines.experimental.channels.takeWhile import kotlinx.coroutines.experimental.channels.toList import kotlinx.coroutines.experimental.runBlocking -import org.apache.commons.math3.random.JDKRandomGenerator import org.apache.commons.math3.random.RandomGenerator import java.time.Duration import java.time.Instant @@ -24,12 +24,10 @@ private fun RandomGenerator.nextDeltaTime(cr: Double): Long { fun generateBlock(start: Instant, length: Long, chain: Chain): NumassBlock { - val events = runBlocking { chain.asChannel().takeWhile { it.timeOffset < length }.toList()} + val events = runBlocking { chain.channel.takeWhile { it.timeOffset < length }.toList()} return SimpleBlock(start, Duration.ofNanos(length), events) } -internal val defaultGenerator = JDKRandomGenerator() - internal val defaultAmplitudeGenerator: RandomGenerator.(NumassEvent?, Long) -> Short = { _, _ -> ((nextDouble() + 2.0) * 100).toShort() } fun buildSimpleEventChain( diff --git a/numass-main/src/main/kotlin/inr/numass/models/mc/BetaSampler.kt b/numass-main/src/main/kotlin/inr/numass/models/mc/BetaSampler.kt new file mode 100644 index 00000000..8d1ca7de --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/mc/BetaSampler.kt @@ -0,0 +1,55 @@ +package inr.numass.models.mc + +import hep.dataforge.fx.plots.PlotManager +import hep.dataforge.kodex.buildMeta +import hep.dataforge.kodex.global +import hep.dataforge.maths.chain.Chain +import hep.dataforge.plots.data.XYFunctionPlot +import hep.dataforge.stat.PolynomialDistribution +import hep.dataforge.stat.fit.ParamSet +import inr.numass.NumassPlugin +import inr.numass.models.sterile.SterileNeutrinoSpectrum + +fun sampleBeta(params: ParamSet): Chain { + TODO() +} + + +fun main(args: Array) { + NumassPlugin().startGlobal() + val pm = PlotManager().apply { startGlobal() } + val meta = buildMeta("model") { + "fast" to true + node("resolution") { + "width" to 1.22e-4 + "tailAlpha" to 1e-2 + } + } + 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.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY) + .setPar("trap", 1.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY) + + val sp = SterileNeutrinoSpectrum(global, meta) + + val spectrumPlot = XYFunctionPlot.plotFunction("spectrum", 14000.0, 18600.0, 500) { + sp.value(it, allPars) + } + + val distribution = PolynomialDistribution(0.0, 5000.0, 3.0); + + val distributionPlot = XYFunctionPlot.plotFunction("distribution", 14000.0, 18500.0, 500) { + 50*distribution.density(18575.0 - it) + } + + pm.getPlotFrame("beta").apply { + add(spectrumPlot) + add(distributionPlot) + } + +} \ No newline at end of file diff --git a/numass-main/src/main/kotlin/inr/numass/models/misc/FunctionSupport.kt b/numass-main/src/main/kotlin/inr/numass/models/misc/FunctionSupport.kt index eec46c5b..d0c3a6c2 100644 --- a/numass-main/src/main/kotlin/inr/numass/models/misc/FunctionSupport.kt +++ b/numass-main/src/main/kotlin/inr/numass/models/misc/FunctionSupport.kt @@ -1,6 +1,6 @@ package inr.numass.models.misc -import hep.dataforge.maths.HyperSquareDomain +import hep.dataforge.maths.domains.HyperSquareDomain import hep.dataforge.values.Values interface FunctionSupport { diff --git a/numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt b/numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt deleted file mode 100644 index c522ae66..00000000 --- a/numass-main/src/main/kotlin/inr/numass/models/sterile/TestModels.kt +++ /dev/null @@ -1,104 +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.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) { - val context = Numass.buildContext() - /* - - - - - */ - 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) - } - -} diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/ScanTree.kt b/numass-main/src/main/kotlin/inr/numass/scripts/utils/ScanTree.kt similarity index 95% rename from numass-main/src/main/kotlin/inr/numass/scripts/ScanTree.kt rename to numass-main/src/main/kotlin/inr/numass/scripts/utils/ScanTree.kt index f043c353..9cd0e2af 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/ScanTree.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/utils/ScanTree.kt @@ -1,4 +1,4 @@ -package inr.numass.scripts +package inr.numass.scripts.utils import hep.dataforge.io.XMLMetaWriter import hep.dataforge.kodex.buildMeta @@ -74,7 +74,7 @@ fun main(args: Array) { val directory = if (args.isNotEmpty()) { args.first() } else { - "." + "" } val path = Paths.get(directory) diff --git a/numass-main/src/test/java/inr/numass/run/NumassSpectrumTest.java b/numass-main/src/test/java/inr/numass/models/NumassSpectrumTest.java similarity index 93% rename from numass-main/src/test/java/inr/numass/run/NumassSpectrumTest.java rename to numass-main/src/test/java/inr/numass/models/NumassSpectrumTest.java index c01e7210..5a000b89 100644 --- a/numass-main/src/test/java/inr/numass/run/NumassSpectrumTest.java +++ b/numass-main/src/test/java/inr/numass/models/NumassSpectrumTest.java @@ -13,13 +13,11 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -package inr.numass.run; +package inr.numass.models; import hep.dataforge.exceptions.NamingException; import hep.dataforge.stat.fit.MINUITPlugin; import hep.dataforge.stat.fit.ParamSet; -import inr.numass.models.BetaSpectrum; -import inr.numass.models.ModularSpectrum; import java.io.File; import java.io.FileNotFoundException; diff --git a/numass-main/src/test/java/inr/numass/models/TestModels.kt b/numass-main/src/test/java/inr/numass/models/TestModels.kt new file mode 100644 index 00000000..8378a20d --- /dev/null +++ b/numass-main/src/test/java/inr/numass/models/TestModels.kt @@ -0,0 +1,93 @@ +/* + * 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 + +import hep.dataforge.context.Context +import hep.dataforge.kodex.step +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.sterile.SterileNeutrinoSpectrum +import org.apache.commons.math3.analysis.BivariateFunction + + +/** + * @param args the command line arguments + */ + +fun main(args: Array) { + val context = Numass.buildContext() + /* + + + + + */ + val meta = MetaBuilder("model") + .putNode(MetaBuilder("resolution") + .putValue("width", 1.22e-4) + .putValue("tailAlpha", 1e-2) + ) + .putNode(MetaBuilder("transmission") +// .putValue("trapping", "function::numass.trap.nominal") + ) + 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.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY) + .setPar("trap", 1.0, 0.01, 0.0, java.lang.Double.POSITIVE_INFINITY) + + for (u in 14000.0..18600.0 step 100.0) { + val oldVal = oldFunc.value(u, allPars); + val 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) + } +} + +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")) { + ResolutionFunction.getAngledTail(meta.getDouble("resolution.tailAlpha"), meta.getDouble("resolution.tailBeta", 0.0)) + } else { + ResolutionFunction.getRealTail() + } + //RangedNamedSetSpectrum beta = new BetaSpectrum(context.io().getFile("FS.txt")); + val sp = ModularSpectrum(BetaSpectrum(), 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) +} + diff --git a/numass-main/src/test/java/inr/numass/models/TestNeLossParametrisation.java b/numass-main/src/test/java/inr/numass/models/TestNeLossParametrisation.java index f658a5e9..44cd5bef 100644 --- a/numass-main/src/test/java/inr/numass/models/TestNeLossParametrisation.java +++ b/numass-main/src/test/java/inr/numass/models/TestNeLossParametrisation.java @@ -41,8 +41,8 @@ public class TestNeLossParametrisation { System.out.println(norm); - frame.add(XYFunctionPlot.plotFunction("old", oldFunction::value, 0, 30, 300)); - frame.add(XYFunctionPlot.plotFunction("new", newFunction::value, 0, 30, 300)); + frame.add(XYFunctionPlot.plotFunction("old", 0, 30, 300, oldFunction::value)); + frame.add(XYFunctionPlot.plotFunction("new", 0, 30, 300, newFunction::value)); } public static UnivariateFunction getSingleScatterFunction( diff --git a/numass-main/src/test/java/inr/numass/models/TransmissionInterpolatorTest.java b/numass-main/src/test/java/inr/numass/models/TransmissionInterpolatorTest.java index b2f4f3a6..dea75d03 100644 --- a/numass-main/src/test/java/inr/numass/models/TransmissionInterpolatorTest.java +++ b/numass-main/src/test/java/inr/numass/models/TransmissionInterpolatorTest.java @@ -33,7 +33,7 @@ public class TransmissionInterpolatorTest { TransmissionInterpolator interpolator = TransmissionInterpolator.fromFile(Global.instance(), "d:\\sterile-new\\loss2014-11\\.dataforge\\merge\\empty_sum.onComplete", "Uset", "CR", 15, 0.8, 19002d); frame.add(DataPlot.plot("data", interpolator.getX(), interpolator.getY())); - frame.add(XYFunctionPlot.plotFunction("interpolated", interpolator::value, interpolator.getXmin(), interpolator.getXmax(), 2000)); + frame.add(XYFunctionPlot.plotFunction("interpolated", interpolator.getXmin(), interpolator.getXmax(), 2000, interpolator::value)); // PrintFunction.printFuntionSimple(new PrintWriter(System.onComplete), interpolator, interpolator.getXmin(), interpolator.getXmax(), 500); }