From a1dea9cac976d74cd723753c796539a88a2e242e Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Sun, 16 May 2021 22:41:21 +0300 Subject: [PATCH] Sterile neutrino model --- build.gradle.kts | 2 +- .../kotlin/ru/inr/mass/data/api/MetaBlock.kt | 5 +- .../ru/inr/mass/data/api/NumassPoint.kt | 3 +- .../inr/mass/data/proto/ProtoNumassPoint.kt | 12 +- .../src/commonMain/resources/data/FS.txt | 193 ++++++++++++++++++ .../ru/inr/mass/maths/DerivativeCalculator.kt | 109 ---------- .../jvmMain/kotlin/ru/inr/mass/models/FSS.kt | 15 +- .../kotlin/ru/inr/mass/models/NumassBeta.kt | 49 +---- .../ru/inr/mass/models/NumassTransmission.kt | 31 +-- .../mass/models/SterileNeutrinoSpectrum.kt | 142 ++++++------- numass-workspace/build.gradle.kts | 3 +- .../ru/inr/mass/scripts/sterileSpectrum.kt | 31 +++ .../inr/mass/workspace/amplitudeSpectrum.kt | 2 + .../kotlin/ru/inr/mass/workspace/workspace.kt | 3 +- .../src/main/resources/data/FS.txt | 193 ++++++++++++++++++ settings.gradle.kts | 2 +- 16 files changed, 538 insertions(+), 257 deletions(-) create mode 100644 numass-model/src/commonMain/resources/data/FS.txt delete mode 100644 numass-model/src/jvmMain/kotlin/ru/inr/mass/maths/DerivativeCalculator.kt create mode 100644 numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt create mode 100644 numass-workspace/src/main/resources/data/FS.txt diff --git a/build.gradle.kts b/build.gradle.kts index 75f6298..dcc09bd 100644 --- a/build.gradle.kts +++ b/build.gradle.kts @@ -13,7 +13,7 @@ allprojects { } val dataforgeVersion by extra("0.4.0") -val kmathVersion by extra("0.3.0-dev-9") +val kmathVersion by extra("0.3.0-dev-10") apiValidation{ validationDisabled = true diff --git a/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/MetaBlock.kt b/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/MetaBlock.kt index 7006abf..d0d555f 100644 --- a/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/MetaBlock.kt +++ b/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/MetaBlock.kt @@ -3,7 +3,7 @@ package ru.inr.mass.data.api import kotlinx.coroutines.flow.* import kotlinx.datetime.Instant import kotlin.time.Duration -import kotlin.time.nanoseconds +import kotlin.time.DurationUnit public interface ParentBlock : NumassBlock { @@ -26,7 +26,8 @@ public class MetaBlock(private val blocks: List) : ParentBlock { override val startTime: Instant get() = blocks.first().startTime - override suspend fun getLength(): Duration = blocks.sumOf { it.getLength().inNanoseconds }.nanoseconds + override suspend fun getLength(): Duration = + Duration.nanoseconds(blocks.sumOf { it.getLength().toDouble(DurationUnit.NANOSECONDS) }) override val events: Flow get() = flow { diff --git a/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/NumassPoint.kt b/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/NumassPoint.kt index 4dc07f9..cfb155c 100644 --- a/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/NumassPoint.kt +++ b/numass-data-model/src/commonMain/kotlin/ru/inr/mass/data/api/NumassPoint.kt @@ -22,6 +22,7 @@ import space.kscience.dataforge.meta.double import space.kscience.dataforge.meta.get import space.kscience.dataforge.meta.int import kotlin.time.Duration +import kotlin.time.DurationUnit import kotlin.time.nanoseconds /** @@ -57,7 +58,7 @@ public interface NumassPoint : ParentBlock { * Get the length key of meta or calculate length as a sum of block lengths. The latter could be a bit slow */ override suspend fun getLength(): Duration = - flowBlocks().filter { it.channel == 0 }.toList().sumOf { it.getLength().inNanoseconds }.nanoseconds + flowBlocks().filter { it.channel == 0 }.toList().sumOf { it.getLength().toDouble(DurationUnit.NANOSECONDS) }.nanoseconds /** * Get all events it all blocks as a single sequence diff --git a/numass-data-proto/src/main/kotlin/ru/inr/mass/data/proto/ProtoNumassPoint.kt b/numass-data-proto/src/main/kotlin/ru/inr/mass/data/proto/ProtoNumassPoint.kt index 732631d..efd2278 100644 --- a/numass-data-proto/src/main/kotlin/ru/inr/mass/data/proto/ProtoNumassPoint.kt +++ b/numass-data-proto/src/main/kotlin/ru/inr/mass/data/proto/ProtoNumassPoint.kt @@ -32,8 +32,6 @@ import java.io.ByteArrayOutputStream import java.io.InputStream import java.util.zip.Inflater import kotlin.time.Duration -import kotlin.time.milliseconds -import kotlin.time.nanoseconds /** * Protobuf based numass point @@ -67,7 +65,7 @@ internal class ProtoNumassPoint( } ?: Instant.DISTANT_PAST override suspend fun getLength(): Duration = meta["acquisition_time"].double?.let { - (it * 1000).toLong().milliseconds + Duration.milliseconds(it * 1000) } ?: super.getLength() override fun toString(): String = "ProtoNumassPoint(index = ${index}, hv = $voltage)" @@ -130,15 +128,15 @@ public class ProtoNumassBlock( } override suspend fun getLength(): Duration = when { - block.length > 0 -> block.length.nanoseconds + block.length > 0 -> Duration.nanoseconds(block.length) parent?.meta["acquisition_time"] != null -> - (parent?.meta["acquisition_time"].double ?: 0.0 * 1000).milliseconds + Duration.milliseconds((parent?.meta["acquisition_time"].double ?: 0.0 * 1000)) else -> { LoggerFactory.getLogger(javaClass) .error("No length information on block. Trying to infer from first and last events") val times = runBlocking { events.map { it.timeOffset }.toList() } val nanos = (times.maxOrNull()!! - times.minOrNull()!!) - nanos.nanoseconds + Duration.nanoseconds(nanos) } } @@ -172,7 +170,7 @@ public class ProtoNumassBlock( override val frames: Flow get() { - val tickSize = block.bin_size.nanoseconds + val tickSize = Duration.nanoseconds(block.bin_size) return block.frames.asFlow().map { frame -> val time = startTime.plus(frame.time, DateTimeUnit.NANOSECOND) val frameData = frame.data_ diff --git a/numass-model/src/commonMain/resources/data/FS.txt b/numass-model/src/commonMain/resources/data/FS.txt new file mode 100644 index 0000000..93b3006 --- /dev/null +++ b/numass-model/src/commonMain/resources/data/FS.txt @@ -0,0 +1,193 @@ +0.000 0.008 +0.097 0.005 +0.197 0.028 +0.297 0.055 +0.397 0.056 +0.497 0.218 +0.597 0.191 +0.697 0.434 +0.797 0.429 +0.897 0.688 +0.997 1.300 +1.097 1.078 +1.197 2.793 +1.297 3.715 +1.397 4.480 +1.497 7.176 +1.597 6.825 +1.697 5.171 +1.797 6.187 +1.897 5.023 +1.997 3.334 +2.097 2.239 +2.197 1.493 +2.297 1.008 +2.397 1.562 +2.647 0.940 +2.897 0.518 +3.147 0.249 +3.397 0.116 +3.647 0.055 +3.897 0.036 +4.397 0.007 +4.897 0.001 +20.881 0.003 +21.881 0.021 +22.881 0.109 +23.881 0.385 +24.881 0.973 +25.881 1.833 +26.881 2.671 +27.881 3.093 +28.881 2.913 +29.881 2.276 +30.881 1.503 +31.881 0.882 +32.881 0.727 +33.881 1.389 +34.881 2.175 +35.881 2.086 +36.881 1.310 +37.881 0.676 +38.725 0.010 +38.881 0.416 +39.881 0.370 +40.881 0.350 +41.881 0.269 +42.732 0.965 +42.881 0.166 +43.405 0.029 +43.881 0.091 +43.963 0.372 +44.147 0.128 +44.881 0.043 +45.881 0.016 +46.881 0.004 +47.913 0.129 +50.599 1.216 +52.553 0.440 +55.109 0.065 +55.852 0.154 +57.004 0.159 +58.092 0.000 +58.592 0.001 +59.092 0.003 +59.592 0.010 +60.092 0.026 +60.592 0.058 +61.092 0.126 +61.592 0.206 +62.092 0.301 +62.592 0.377 +63.092 0.418 +63.592 0.377 +64.092 0.301 +64.386 0.003 +64.592 0.206 +64.886 0.007 +65.092 0.126 +65.386 0.023 +65.592 0.058 +65.886 0.060 +66.092 0.026 +66.386 0.133 +66.592 0.010 +66.886 0.288 +67.092 0.003 +67.386 0.471 +67.592 0.001 +67.886 0.688 +68.092 0.000 +68.386 0.863 +68.886 0.956 +69.386 0.863 +69.886 0.688 +70.386 0.471 +70.886 0.288 +71.386 0.133 +71.725 0.306 +71.886 0.060 +72.386 0.023 +72.886 0.007 +73.386 0.003 +74.820 0.245 +76.169 0.088 +76.868 0.100 +77.221 0.273 +79.427 0.020 +80.865 0.238 +81.965 0.137 +83.429 0.151 +84.170 0.212 +84.218 0.112 +86.123 0.014 +87.374 0.010 +88.259 0.009 +88.876 0.013 +89.871 0.026 +90.690 0.023 +91.784 0.052 +93.247 0.178 +94.333 0.133 +96.192 0.026 +96.701 0.054 +97.543 0.023 +98.514 0.005 +98.840 0.010 +100.263 0.014 +100.784 0.003 +101.620 0.003 +102.426 0.005 +102.842 0.001 +103.170 0.001 +103.594 0.006 +104.236 0.002 +105.008 0.001 +105.799 0.002 +106.990 0.006 +108.711 0.010 +109.189 0.008 +109.975 0.007 +111.148 0.005 +112.339 0.013 +113.145 0.010 +113.882 0.005 +114.892 0.002 +115.612 0.002 +116.455 0.001 +117.594 0.005 +118.481 0.023 +119.245 0.023 +120.360 0.009 +121.764 0.013 +123.594 0.009 +124.247 0.005 +125.709 0.012 +127.715 0.003 +129.373 0.002 +130.271 0.004 +132.887 0.060 +133.402 0.025 +134.813 0.082 +135.371 0.006 +136.379 0.005 +136.916 0.003 +138.243 0.008 +139.737 0.010 +141.093 0.006 +142.461 0.047 +144.001 0.004 +144.391 0.007 +147.073 0.021 +148.311 0.015 +148.895 0.001 +150.849 0.004 +151.442 0.001 +152.854 0.000 +154.169 0.002 +156.093 0.001 +157.003 0.003 +158.134 0.003 +159.271 0.002 +162.054 0.007 +164.173 0.002 \ No newline at end of file diff --git a/numass-model/src/jvmMain/kotlin/ru/inr/mass/maths/DerivativeCalculator.kt b/numass-model/src/jvmMain/kotlin/ru/inr/mass/maths/DerivativeCalculator.kt deleted file mode 100644 index 9c3532f..0000000 --- a/numass-model/src/jvmMain/kotlin/ru/inr/mass/maths/DerivativeCalculator.kt +++ /dev/null @@ -1,109 +0,0 @@ -/* - * Copyright 2015 Alexander Nozik. - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * - * http://www.apache.org/licenses/LICENSE-2.0 - * - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ -package ru.inr.mass.maths - -import hep.dataforge.stat.fit.Param -import hep.dataforge.stat.parametric.ParametricValue - -/** - * - * DerivativeCalculator class. - * - * @author Alexander Nozik - * @version $Id: $Id - */ -object DerivativeCalculator { - private const val numPoints = 3 - - /** - * Calculates finite differences derivative via 3 points differentiator. - * - * @param function a [hep.dataforge.stat.parametric.ParametricValue] object. - * @param point a [hep.dataforge.stat.fit.ParamSet] object. - * @param parName a [String] object. - * @return a double. - */ - fun calculateDerivative(function: ParametricValue?, point: ParamSet, parName: String?): Double { - val projection: UnivariateFunction = ParametricUtils.getNamedProjection(function, parName, point) - val par: Param = point.getByName(parName) - val diff = - FiniteDifferencesDifferentiator(numPoints, par.getErr() / 2.0, par.getLowerBound(), par.getUpperBound()) - val derivative: UnivariateDifferentiableFunction = diff.differentiate(projection) - val x = DerivativeStructure(1, 1, 0, point.getDouble(parName)) - val y: DerivativeStructure = derivative.value(x) - return y.getPartialDerivative(1) - } - - /** - * Calculates finite differences derivative via 3 points differentiator. - * - * @param function a [org.apache.commons.math3.analysis.UnivariateFunction] object. - * @param point a double. - * @param step a double. - * @return a double. - */ - fun calculateDerivative(function: UnivariateFunction?, point: Double, step: Double): Double { - val diff = FiniteDifferencesDifferentiator(numPoints, step) - val derivative: UnivariateDifferentiableFunction = diff.differentiate(function) - val x = DerivativeStructure(1, 1, 0, point) - val y: DerivativeStructure = derivative.value(x) - return y.getPartialDerivative(1) - } - - /** - * - * providesValidDerivative. - * - * @param function a [hep.dataforge.stat.parametric.ParametricValue] object. - * @param point a [hep.dataforge.stat.fit.ParamSet] object. - * @param tolerance a double. - * @param parName a [String] object. - * @return a boolean. - */ - fun providesValidDerivative( - function: ParametricValue, - point: ParamSet, - tolerance: Double, - parName: String? - ): Boolean { - if (!function.providesDeriv(parName)) { - return false - } - val calculatedDeriv = calculateDerivative(function, point, parName) - val providedDeriv: Double = function.derivValue(parName, point) - return safeRelativeDifference(calculatedDeriv, providedDeriv) <= tolerance - } - - /** - * Returns safe from (no devision by zero) relative difference between two - * input values - * - * @param val1 - * @param val2 - * @return - */ - private fun safeRelativeDifference(val1: Double, val2: Double): Double { - if (Precision.equals(val1, val2, Precision.EPSILON)) { - return 0 - } - val average: Double = Math.abs(val1 + val2) / 2 - return if (average > Precision.EPSILON) { - Math.abs(val1 - val2) / average - } else { - Double.POSITIVE_INFINITY - } - } -} \ No newline at end of file diff --git a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/FSS.kt b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/FSS.kt index a022a9d..9a64c78 100644 --- a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/FSS.kt +++ b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/FSS.kt @@ -24,7 +24,7 @@ import space.kscience.kmath.structures.DoubleBuffer @OptIn(UnstableKMathAPI::class) -public class FSS(public val ps: DoubleBuffer, public val es: DoubleBuffer) : ColumnarData { +public class FSS(public val es: DoubleBuffer, public val ps: DoubleBuffer) : ColumnarData { override val size: Int get() = ps.size @@ -37,5 +37,18 @@ public class FSS(public val ps: DoubleBuffer, public val es: DoubleBuffer) : Col public companion object { public val p: Symbol by symbol public val e: Symbol by symbol + + public val default: FSS by lazy { + val stream = FSS::class.java.getResourceAsStream("/data/FS.txt") ?: error("Default FS resource not found") + stream.use { inputStream -> + val data = inputStream.bufferedReader().lineSequence().map { + val (e, p) = it.split("\t") + e.toDouble() to p.toDouble() + }.toList() + val es = DoubleBuffer(data.size) { data[it].first } + val ps = DoubleBuffer(data.size) { data[it].second } + FSS(es, ps) + } + } } } \ No newline at end of file diff --git a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassBeta.kt b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassBeta.kt index 3a82931..2d8d258 100644 --- a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassBeta.kt +++ b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassBeta.kt @@ -167,9 +167,9 @@ public object NumassBeta : DifferentiableKernel { override val x: Symbol = StringSymbol("fs") override val y: Symbol = StringSymbol("eIn") - override fun invoke(x: Double, y: Double, arguments: Map): Double { + override fun invoke(fs: Double, eIn: Double, arguments: Map): Double { val e0 = arguments.getValue(e0) - return rootsterile(y, e0 - x, arguments) + return rootsterile(eIn, e0 - fs, arguments) } override fun derivativeOrNull(symbols: List): Kernel? = when (symbols.size) { @@ -180,47 +180,12 @@ public object NumassBeta : DifferentiableKernel { } else -> null } -// -// override fun getDefaultParameter(name: String): Double { -// return when (name) { -// "mnu2", "U2", "msterile2" -> 0.0 -// else -> 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) -// } -// -// /** -// * 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) -// } -// -// } - private const val K: Double = 1E-23 - public val e0: Symbol by symbol - public val mnu2: Symbol by symbol - public val msterile2: Symbol by symbol - public val u2: Symbol by symbol + private const val K: Double = 1E-23 + public val e0: Symbol by symbol + public val mnu2: Symbol by symbol + public val msterile2: Symbol by symbol + public val u2: Symbol by symbol } diff --git a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassTransmission.kt b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassTransmission.kt index 8b059a7..efedac0 100644 --- a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassTransmission.kt +++ b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassTransmission.kt @@ -27,8 +27,8 @@ public fun PiecewisePolynomial.asFunction(defaultValue: Double = 0.0): U * @author [Alexander Nozik](mailto:altavir@gmail.com) */ public class NumassTransmission( - public val trapFunc: Kernel, - private val adjustX: Boolean = false, + public val trapFunc: Kernel = defaultTrapping, +// private val adjustX: Boolean = false, ) : DifferentiableKernel { // private val trapFunc: Kernel = if (meta.hasValue("trapping")) { // val trapFuncStr = meta.getString("trapping") @@ -83,7 +83,7 @@ public class NumassTransmission( // } //trapping part - val trap = arguments.getOrElse(trap) { 1.0 } * trapFunc(x, y, arguments) + val trap = arguments[trap] ?: 1.0 * trapFunc(x, y, arguments) return loss + trap } @@ -97,15 +97,16 @@ public class NumassTransmission( private fun getX(arguments: Map, eIn: Double, adjustX: Boolean = false): Double { + val thickness = arguments[thickness] ?: 0.0 return if (adjustX) { //From our article - arguments.getValue(thickness) * ln(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6 + thickness * ln(eIn / ION_POTENTIAL) * eIn * ION_POTENTIAL / 1.9580741410115568e6 } else { - arguments.getValue(thickness) + thickness } } - private fun p0(eIn: Double, set: Map): Double = getLossProbability(0, getX(set, eIn)) + internal fun p0(eIn: Double, set: Map): Double = getLossProbability(0, getX(set, eIn)) private fun getGunLossProbabilities(X: Double): List { val res = ArrayList() @@ -132,7 +133,7 @@ public class NumassTransmission( return res } - fun getGunZeroLossProb(x: Double): Double { + private fun getGunZeroLossProb(x: Double): Double { return exp(-x) } @@ -207,9 +208,10 @@ public class NumassTransmission( return res } - fun getLossProbabilities(x: Double): List = lossProbCache.getOrPut(x) { calculateLossProbabilities(x) } + private fun getLossProbabilities(x: Double): List = + lossProbCache.getOrPut(x) { calculateLossProbabilities(x) } - fun getLossProbability(order: Int, X: Double): Double { + private fun getLossProbability(order: Int, X: Double): Double { if (order == 0) { return if (X > 0) { 1 / X * (1 - exp(-X)) @@ -225,7 +227,7 @@ public class NumassTransmission( } } - fun getLossValue(order: Int, Ei: Double, Ef: Double): Double { + private fun getLossValue(order: Int, Ei: Double, Ef: Double): Double { return when { Ei - Ef < 5.0 -> 0.0 Ei - Ef >= getMargin(order) -> 0.0 @@ -241,7 +243,7 @@ public class NumassTransmission( * @param Ef * @return */ - fun getLossValue(probs: List, Ei: Double, Ef: Double): Double { + private fun getLossValue(probs: List, Ei: Double, Ef: Double): Double { var sum = 0.0 for (i in 1 until probs.size) { sum += probs[i] * getLossValue(i, Ei, Ef) @@ -305,7 +307,7 @@ public class NumassTransmission( * @param Ef * @return */ - fun getTotalLossValue(x: Double, Ei: Double, Ef: Double): Double { + private fun getTotalLossValue(x: Double, Ei: Double, Ef: Double): Double { return if (x == 0.0) { 0.0 } else { @@ -464,6 +466,11 @@ public class NumassTransmission( // } // // } + + internal val defaultTrapping: Kernel = Kernel { e, u, arguments -> + val d = e - u + 0.99797 - 3.05346E-7 * d - 5.45738E-10 * d.pow(2.0) - 6.36105E-14 * d.pow(3.0) + } } } diff --git a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt index f51b6ed..6921cbe 100644 --- a/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt +++ b/numass-model/src/jvmMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt @@ -6,24 +6,30 @@ package inr.numass.models.sterile -import ru.inr.mass.models.DifferentiableKernel -import ru.inr.mass.models.DifferentiableSpectrum -import ru.inr.mass.models.FSS -import ru.inr.mass.models.Spectrum -import space.kscience.dataforge.context.Context +import inr.numass.models.sterile.NumassBeta.e0 +import inr.numass.models.sterile.NumassBeta.mnu2 +import inr.numass.models.sterile.NumassBeta.msterile2 +import inr.numass.models.sterile.NumassBeta.u2 +import inr.numass.models.sterile.NumassTransmission.Companion.thickness +import inr.numass.models.sterile.NumassTransmission.Companion.trap +import ru.inr.mass.models.* +import space.kscience.kmath.expressions.derivative +import space.kscience.kmath.integration.integrate +import space.kscience.kmath.integration.value import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.operations.DoubleField +import kotlin.math.min /** * @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( - context: Context, +public class SterileNeutrinoSpectrum( public val source: DifferentiableKernel = NumassBeta, - public val transmission: DifferentiableKernel, - public val resolution: DifferentiableKernel, - public val fss: FSS? + public val transmission: DifferentiableKernel = NumassTransmission(), + public val resolution: DifferentiableKernel = NumassResolution(), + public val fss: FSS? = FSS.default, ) : DifferentiableSpectrum { @@ -36,111 +42,93 @@ class SterileNeutrinoSpectrum( //private val fast: Boolean = configuration.getBoolean("fast", true) - override fun invoke(x: Double, arguments: Map): Double { - TODO("Not yet implemented") + override fun invoke(u: Double, arguments: Map): Double { + return convolute(u, source, transRes, arguments) } override fun derivativeOrNull(symbols: List): Spectrum? { - TODO("Not yet implemented") - } - - 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() + if (symbols.isEmpty()) return this + return when (symbols.singleOrNull() ?: TODO("First derivatives only")) { + u2, msterile2, mnu2, e0 -> Spectrum { u, arguments -> + convolute(u, source.derivative(symbols), transRes, arguments) + } + thickness, trap -> Spectrum { u, arguments -> + convolute(u, source, transRes.derivative(symbols), arguments) + } + else -> null } } - 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-Legendre integration * * @param u * @param sourceFunction * @param transResFunction - * @param set + * @param arguments * @return */ - private fun integrate( + private fun convolute( u: Double, - sourceFunction: DifferentiableKernel, - transResFunction: DifferentiableKernel, - set: Map, + sourceFunction: Kernel, + transResFunction: Kernel, + arguments: Map, ): Double { - val eMax = set.getDouble("E0") + 5.0 + val eMax = arguments.getValue(e0) + 5.0 if (u >= eMax) { return 0.0 } - val integrator: UnivariateIntegrator<*> = if (fast) { - when { - eMax - u < 300 -> NumassIntegrator.getFastInterator() - eMax - u > 2000 -> NumassIntegrator.getHighDensityIntegrator() - else -> NumassIntegrator.getDefaultIntegrator() - } +// val integrator: UnivariateIntegrator<*> = if (fast) { +// when { +// eMax - u < 300 -> getFastInterator() +// eMax - u > 2000 -> getHighDensityIntegrator() +// else -> getDefaultIntegrator() +// } +// } else { +// getHighDensityIntegrator() +// } - } else { - NumassIntegrator.getHighDensityIntegrator() - } - - return integrator.integrate(u, eMax) { eIn -> - sumByFSS(eIn, sourceFunction, set) * transResFunction.value(eIn, - u, - set) - } + return DoubleField.integrate(u..eMax) { eIn -> + sumByFSS(eIn, sourceFunction, arguments) * transResFunction(eIn, u, arguments) + }.value ?: error("Integration failed") } - private fun sumByFSS(eIn: Double, sourceFunction: ParametricBiFunction, set: Values): Double { + private fun sumByFSS(eIn: Double, sourceFunction: Kernel, arguments: Map): Double { return if (fss == null) { - sourceFunction.value(0.0, eIn, set) + sourceFunction(0.0, eIn, arguments) } else { - (0 until fss.size()).sumByDouble { fss.getP(it) * sourceFunction.value(fss.getE(it), eIn, set) } + (0 until fss.size).sumOf { fss.ps[it] * sourceFunction(fss.es[it], eIn, arguments) } } } private inner class TransRes : DifferentiableKernel { - override fun providesDeriv(name: String): Boolean { - return true + override fun invoke(eIn: Double, u: Double, arguments: Map): Double { + val p0 = NumassTransmission.p0(eIn, arguments) + return p0 * resolution(eIn, u, arguments) + lossRes(transmission, eIn, u, arguments) } - 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 derivativeOrNull(symbols: List): Kernel? { + if (symbols.isEmpty()) return this + return when (symbols.singleOrNull() ?: TODO("First derivatives only")) { + thickness -> null//TODO implement p0 derivative + trap -> Kernel { eIn, u, arguments -> lossRes(transmission.derivative(symbols), eIn, u, arguments) } + else -> null } } - override fun value(eIn: Double, u: Double, set: Values): Double { - - val p0 = LossCalculator.p0(set, eIn) - 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) } + private fun lossRes(transFunc: Kernel, eIn: Double, u: Double, arguments: Map): Double { + val integrand = { eOut: Double -> transFunc(eIn, eOut, arguments) * resolution(eOut, u, arguments) } val border = u + 30 - val firstPart = NumassIntegrator.getFastInterator().integrate(u, Math.min(eIn, border), integrand) + val firstPart = DoubleField.integrate(u..min(eIn, border), function = integrand).value + ?: error("Integration failed") val secondPart: Double = if (eIn > border) { - if (fast) { - NumassIntegrator.getDefaultIntegrator().integrate(border, eIn, integrand) - } else { - NumassIntegrator.getHighDensityIntegrator().integrate(border, eIn, integrand) - } + DoubleField.integrate(border..eIn, function = integrand).value ?: error("Integration failed") } else { 0.0 } @@ -149,9 +137,5 @@ class SterileNeutrinoSpectrum( } - companion object { - - private val list = arrayOf("X", "trap", "E0", "mnu2", "msterile2", "U2") - } - + public companion object } diff --git a/numass-workspace/build.gradle.kts b/numass-workspace/build.gradle.kts index 19d592d..10782a5 100644 --- a/numass-workspace/build.gradle.kts +++ b/numass-workspace/build.gradle.kts @@ -10,11 +10,12 @@ kotlin { } val dataforgeVersion: String by rootProject.extra -val plotlyVersion: String by rootProject.extra("0.4.0-dev-1") +val plotlyVersion: String by rootProject.extra("0.4.0") val kmathVersion: String by rootProject.extra dependencies { implementation(project(":numass-data-proto")) + implementation(project(":numass-model")) implementation("space.kscience:dataforge-workspace:$dataforgeVersion") implementation("space.kscience:plotlykt-core:$plotlyVersion") implementation("space.kscience:kmath-histograms:$kmathVersion") diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt new file mode 100644 index 0000000..da263a3 --- /dev/null +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt @@ -0,0 +1,31 @@ +package ru.inr.mass.scripts + +import inr.numass.models.sterile.NumassBeta.e0 +import inr.numass.models.sterile.NumassBeta.mnu2 +import inr.numass.models.sterile.NumassBeta.msterile2 +import inr.numass.models.sterile.NumassBeta.u2 +import inr.numass.models.sterile.SterileNeutrinoSpectrum +import space.kscience.kmath.misc.Symbol +import space.kscience.kmath.real.step +import space.kscience.kmath.structures.asIterable +import space.kscience.plotly.Plotly +import space.kscience.plotly.makeFile +import space.kscience.plotly.models.ScatterMode +import space.kscience.plotly.scatter + +fun main() { + val spectrum = SterileNeutrinoSpectrum() + val args: Map = mapOf( + mnu2 to 0.0, + e0 to 18575.0, + msterile2 to 1e6, + u2 to 1e-2 + ) + Plotly.plot { + scatter { + mode = ScatterMode.lines + x.numbers = (14000.0..18600.0 step 1.0).asIterable() + y.numbers = x.numbers.map { spectrum(it.toDouble(), args) } + } + }.makeFile() +} \ No newline at end of file diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/amplitudeSpectrum.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/amplitudeSpectrum.kt index e67e241..d52b17a 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/amplitudeSpectrum.kt +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/amplitudeSpectrum.kt @@ -8,6 +8,7 @@ import space.kscience.dataforge.context.warn import space.kscience.kmath.histogram.UnivariateHistogram import space.kscience.kmath.histogram.center import space.kscience.kmath.histogram.put +import space.kscience.kmath.misc.UnstableKMathAPI import space.kscience.kmath.structures.DoubleBuffer import space.kscience.kmath.structures.asBuffer @@ -41,6 +42,7 @@ fun Collection.spectrum(): UnivariateHistogram { /** * Re-shape the spectrum with the increased bin size and range. Throws a error if new bin is smaller than before. */ +@OptIn(UnstableKMathAPI::class) fun UnivariateHistogram.reShape( binSize: Int, channelRange: IntRange, diff --git a/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/workspace.kt b/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/workspace.kt index f249ae2..2da8f3e 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/workspace.kt +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/workspace/workspace.kt @@ -4,7 +4,8 @@ import ru.inr.mass.data.proto.NumassProtoPlugin import space.kscience.dataforge.workspace.Workspace val NUMASS = Workspace { - context("NUMASS") { + context{ + name("NUMASS") plugin(NumassProtoPlugin) } } \ No newline at end of file diff --git a/numass-workspace/src/main/resources/data/FS.txt b/numass-workspace/src/main/resources/data/FS.txt new file mode 100644 index 0000000..93b3006 --- /dev/null +++ b/numass-workspace/src/main/resources/data/FS.txt @@ -0,0 +1,193 @@ +0.000 0.008 +0.097 0.005 +0.197 0.028 +0.297 0.055 +0.397 0.056 +0.497 0.218 +0.597 0.191 +0.697 0.434 +0.797 0.429 +0.897 0.688 +0.997 1.300 +1.097 1.078 +1.197 2.793 +1.297 3.715 +1.397 4.480 +1.497 7.176 +1.597 6.825 +1.697 5.171 +1.797 6.187 +1.897 5.023 +1.997 3.334 +2.097 2.239 +2.197 1.493 +2.297 1.008 +2.397 1.562 +2.647 0.940 +2.897 0.518 +3.147 0.249 +3.397 0.116 +3.647 0.055 +3.897 0.036 +4.397 0.007 +4.897 0.001 +20.881 0.003 +21.881 0.021 +22.881 0.109 +23.881 0.385 +24.881 0.973 +25.881 1.833 +26.881 2.671 +27.881 3.093 +28.881 2.913 +29.881 2.276 +30.881 1.503 +31.881 0.882 +32.881 0.727 +33.881 1.389 +34.881 2.175 +35.881 2.086 +36.881 1.310 +37.881 0.676 +38.725 0.010 +38.881 0.416 +39.881 0.370 +40.881 0.350 +41.881 0.269 +42.732 0.965 +42.881 0.166 +43.405 0.029 +43.881 0.091 +43.963 0.372 +44.147 0.128 +44.881 0.043 +45.881 0.016 +46.881 0.004 +47.913 0.129 +50.599 1.216 +52.553 0.440 +55.109 0.065 +55.852 0.154 +57.004 0.159 +58.092 0.000 +58.592 0.001 +59.092 0.003 +59.592 0.010 +60.092 0.026 +60.592 0.058 +61.092 0.126 +61.592 0.206 +62.092 0.301 +62.592 0.377 +63.092 0.418 +63.592 0.377 +64.092 0.301 +64.386 0.003 +64.592 0.206 +64.886 0.007 +65.092 0.126 +65.386 0.023 +65.592 0.058 +65.886 0.060 +66.092 0.026 +66.386 0.133 +66.592 0.010 +66.886 0.288 +67.092 0.003 +67.386 0.471 +67.592 0.001 +67.886 0.688 +68.092 0.000 +68.386 0.863 +68.886 0.956 +69.386 0.863 +69.886 0.688 +70.386 0.471 +70.886 0.288 +71.386 0.133 +71.725 0.306 +71.886 0.060 +72.386 0.023 +72.886 0.007 +73.386 0.003 +74.820 0.245 +76.169 0.088 +76.868 0.100 +77.221 0.273 +79.427 0.020 +80.865 0.238 +81.965 0.137 +83.429 0.151 +84.170 0.212 +84.218 0.112 +86.123 0.014 +87.374 0.010 +88.259 0.009 +88.876 0.013 +89.871 0.026 +90.690 0.023 +91.784 0.052 +93.247 0.178 +94.333 0.133 +96.192 0.026 +96.701 0.054 +97.543 0.023 +98.514 0.005 +98.840 0.010 +100.263 0.014 +100.784 0.003 +101.620 0.003 +102.426 0.005 +102.842 0.001 +103.170 0.001 +103.594 0.006 +104.236 0.002 +105.008 0.001 +105.799 0.002 +106.990 0.006 +108.711 0.010 +109.189 0.008 +109.975 0.007 +111.148 0.005 +112.339 0.013 +113.145 0.010 +113.882 0.005 +114.892 0.002 +115.612 0.002 +116.455 0.001 +117.594 0.005 +118.481 0.023 +119.245 0.023 +120.360 0.009 +121.764 0.013 +123.594 0.009 +124.247 0.005 +125.709 0.012 +127.715 0.003 +129.373 0.002 +130.271 0.004 +132.887 0.060 +133.402 0.025 +134.813 0.082 +135.371 0.006 +136.379 0.005 +136.916 0.003 +138.243 0.008 +139.737 0.010 +141.093 0.006 +142.461 0.047 +144.001 0.004 +144.391 0.007 +147.073 0.021 +148.311 0.015 +148.895 0.001 +150.849 0.004 +151.442 0.001 +152.854 0.000 +154.169 0.002 +156.093 0.001 +157.003 0.003 +158.134 0.003 +159.271 0.002 +162.054 0.007 +164.173 0.002 \ No newline at end of file diff --git a/settings.gradle.kts b/settings.gradle.kts index 5a13cd7..f03c7b7 100644 --- a/settings.gradle.kts +++ b/settings.gradle.kts @@ -6,7 +6,7 @@ pluginManagement { gradlePluginPortal() } - val toolsVersion = "0.9.6" + val toolsVersion = "0.9.7" val kotlinVersion = "1.5.0" plugins {