Sterile neutrino model

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

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<NumassBlock>) : 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<NumassEvent>
         get() = flow {

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

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<NumassFrame>
         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_

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

numass-model/src/jvmMain/kotlin/ru/inr/mass/maths/DerivativeCalculator.kt

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

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<Double> {
+public class FSS(public val es: DoubleBuffer, public val ps: DoubleBuffer) : ColumnarData<Double> {
 
     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)
+            }
+        }
     }
 }

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<Symbol, Double>): Double {
+    override fun invoke(fs: Double, eIn: Double, arguments: Map<Symbol, Double>): Double {
         val e0 = arguments.getValue(e0)
-        return rootsterile(y, e0 - x, arguments)
+        return rootsterile(eIn, e0 - fs, arguments)
     }
 
     override fun derivativeOrNull(symbols: List<Symbol>): Kernel? = when (symbols.size) {
@@ -180,41 +180,6 @@ 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

numass-model/src/jvmMain/kotlin/ru/inr/mass/models/NumassTransmission.kt

@@ -27,8 +27,8 @@ public fun PiecewisePolynomial<Double>.asFunction(defaultValue: Double = 0.0): U
  * @author [Alexander Nozik](mailto:altavir@gmail.com)
  */
 public class NumassTransmission(
-    public val trapFunc: Kernel,
+    public val trapFunc: Kernel = defaultTrapping,
-    private val adjustX: Boolean = false,
+//    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<Symbol, Double>, 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<Symbol, Double>): Double = getLossProbability(0, getX(set, eIn))
+        internal fun p0(eIn: Double, set: Map<Symbol, Double>): Double = getLossProbability(0, getX(set, eIn))
 
         private fun getGunLossProbabilities(X: Double): List<Double> {
             val res = ArrayList<Double>()
@@ -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<Double> = lossProbCache.getOrPut(x) { calculateLossProbabilities(x) }
+        private fun getLossProbabilities(x: Double): List<Double> =
+            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<Double>, Ei: Double, Ef: Double): Double {
+        private fun getLossValue(probs: List<Double>, 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)
+        }
     }
 
 }

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 inr.numass.models.sterile.NumassBeta.e0
-import ru.inr.mass.models.DifferentiableSpectrum
+import inr.numass.models.sterile.NumassBeta.mnu2
-import ru.inr.mass.models.FSS
+import inr.numass.models.sterile.NumassBeta.msterile2
-import ru.inr.mass.models.Spectrum
+import inr.numass.models.sterile.NumassBeta.u2
-import space.kscience.dataforge.context.Context
+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(
+public class SterileNeutrinoSpectrum(
-    context: Context,
     public val source: DifferentiableKernel = NumassBeta,
-    public val transmission: DifferentiableKernel,
+    public val transmission: DifferentiableKernel = NumassTransmission(),
-    public val resolution: DifferentiableKernel,
+    public val resolution: DifferentiableKernel = NumassResolution(),
-    public val fss: FSS?
+    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<Symbol, Double>): Double {
+    override fun invoke(u: Double, arguments: Map<Symbol, Double>): Double {
-        TODO("Not yet implemented")
+        return convolute(u, source, transRes, arguments)
     }
 
     override fun derivativeOrNull(symbols: List<Symbol>): Spectrum? {
-        TODO("Not yet implemented")
+        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 ->
-    override fun derivValue(parName: String, u: Double, set: Values): Double {
+                convolute(u, source, transRes.derivative(symbols), arguments)
-        return when (parName) {
+            }
-            "U2", "msterile2", "mnu2", "E0" -> integrate(u, source.derivative(parName), transRes, set)
+            else -> null
-            "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-Legendre integration
      *
      * @param u
      * @param sourceFunction
      * @param transResFunction
-     * @param set
+     * @param arguments
      * @return
      */
-    private fun integrate(
+    private fun convolute(
         u: Double,
-        sourceFunction: DifferentiableKernel,
+        sourceFunction: Kernel,
-        transResFunction: DifferentiableKernel,
+        transResFunction: Kernel,
-        set: Map<Symbol, Double>,
+        arguments: Map<Symbol, Double>,
     ): 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) {
+//        val integrator: UnivariateIntegrator<*> = if (fast) {
-            when {
+//            when {
-                eMax - u < 300 -> NumassIntegrator.getFastInterator()
+//                eMax - u < 300 -> getFastInterator()
-                eMax - u > 2000 -> NumassIntegrator.getHighDensityIntegrator()
+//                eMax - u > 2000 -> getHighDensityIntegrator()
-                else -> NumassIntegrator.getDefaultIntegrator()
+//                else -> getDefaultIntegrator()
+//            }
+//        } else {
+//            getHighDensityIntegrator()
+//        }
+
+        return DoubleField.integrate(u..eMax) { eIn ->
+            sumByFSS(eIn, sourceFunction, arguments) * transResFunction(eIn, u, arguments)
+        }.value ?: error("Integration failed")
     }
 
-        } else {
+    private fun sumByFSS(eIn: Double, sourceFunction: Kernel, arguments: Map<Symbol, Double>): Double {
-            NumassIntegrator.getHighDensityIntegrator()
-        }
-
-        return integrator.integrate(u, eMax) { eIn ->
-            sumByFSS(eIn, sourceFunction, set) * transResFunction.value(eIn,
-                u,
-                set)
-        }
-    }
-
-    private fun sumByFSS(eIn: Double, sourceFunction: ParametricBiFunction, set: Values): 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 {
+        override fun invoke(eIn: Double, u: Double, arguments: Map<Symbol, Double>): Double {
-            return true
+            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 {
+        override fun derivativeOrNull(symbols: List<Symbol>): Kernel? {
-            return when (parName) {
+            if (symbols.isEmpty()) return this
-                "X" -> throw NotDefinedException()//TODO implement p0 derivative
+            return when (symbols.singleOrNull() ?: TODO("First derivatives only")) {
-                "trap" -> lossRes(transmission.derivative(parName), eIn, u, set)
+                thickness -> null//TODO implement p0 derivative
-                else -> super.derivValue(parName, eIn, u, set)
+                trap -> Kernel { eIn, u, arguments -> lossRes(transmission.derivative(symbols), eIn, u, arguments) }
+                else -> null
             }
         }
 
-        override fun value(eIn: Double, u: Double, set: Values): Double {
+        private fun lossRes(transFunc: Kernel, eIn: Double, u: Double, arguments: Map<Symbol, Double>): Double {
-
+            val integrand = { eOut: Double -> transFunc(eIn, eOut, arguments) * resolution(eOut, u, arguments) }
-            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) }
 
             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) {
+                DoubleField.integrate(border..eIn, function = integrand).value ?: error("Integration failed")
-                    NumassIntegrator.getDefaultIntegrator().integrate(border, eIn, integrand)
-                } else {
-                    NumassIntegrator.getHighDensityIntegrator().integrate(border, eIn, integrand)
-                }
             } else {
                 0.0
             }
@@ -149,9 +137,5 @@ class SterileNeutrinoSpectrum(
 
     }
 
-    companion object {
+    public companion object
-
-        private val list = arrayOf("X", "trap", "E0", "mnu2", "msterile2", "U2")
-    }
-
 }

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

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<Symbol, Double> = 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()
+}

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<NumassPoint>.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,

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

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

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 {