diff --git a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt index 9f69999..a9ce83f 100644 --- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt +++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/NumassTransmission.kt @@ -10,9 +10,7 @@ import space.kscience.kmath.expressions.symbol import space.kscience.kmath.functions.PiecewisePolynomial import space.kscience.kmath.functions.UnivariateFunction import space.kscience.kmath.functions.asFunction -import space.kscience.kmath.integration.GaussIntegrator -import space.kscience.kmath.integration.integrate -import space.kscience.kmath.integration.value +import space.kscience.kmath.integration.* import space.kscience.kmath.operations.DoubleField import kotlin.jvm.Synchronized import kotlin.math.* @@ -251,7 +249,7 @@ public class NumassTransmission( * @return */ private fun getMargin(order: Int): Double { - return 80 + order * 50.0 + return 50 + order * 50.0 } /** @@ -264,13 +262,12 @@ public class NumassTransmission( @Synchronized private fun getNextLoss(margin: Double, loss: UnivariateFunction): PiecewisePolynomial { val res = { x: Double -> - integrator.integrate(5.0..margin) { y -> + DoubleField.simpsonIntegrator.integrate(5.0..margin, IntegrandMaxCalls(200)) { y -> loss(x - y) * singleScatterFunction(y) }.value } - return res.cache(0.0..margin, 200) - + return res.cache(0.0..margin, 100) } /** @@ -316,7 +313,6 @@ public class NumassTransmission( * порог по вероятности, до которого вычисляются компоненты функции потерь */ private const val SCATTERING_PROBABILITY_THRESHOLD = 1e-3 - private val integrator = GaussIntegrator(DoubleField) private val lossProbCache = HashMap>(100) @@ -343,48 +339,48 @@ public class NumassTransmission( } - /** - * A generic loss function for numass experiment in "Lobashev" - * parameterization - * - * @param exPos - * @param ionPos - * @param exW - * @param ionW - * @param exIonRatio - * @return - */ - public fun getSingleScatterFunction( - exPos: Double, - ionPos: Double, - exW: Double, - ionW: Double, - exIonRatio: Double, - ): UnivariateFunction { - val func: UnivariateFunction = { eps: Double -> - if (eps <= 0) { - 0.0 - } else { - val z1 = eps - exPos - val ex = exIonRatio * exp(-2.0 * z1 * z1 / exW / exW) - - val z = 4.0 * (eps - ionPos) * (eps - ionPos) - val ion = 1 / (1 + z / ionW / ionW) - - if (eps < exPos) { - ex - } else { - max(ex, ion) - } - } - } - - val cutoff = 25.0 - //caclulating lorentz integral analythically - val tailNorm = (atan((ionPos - cutoff) * 2.0 / ionW) + 0.5 * PI) * ionW / 2.0 - val norm: Double = integrator.integrate(range = 0.0..cutoff, function = func).value + tailNorm - return { e -> func(e) / norm } - } +// /** +// * A generic loss function for numass experiment in "Lobashev" +// * parameterization +// * +// * @param exPos +// * @param ionPos +// * @param exW +// * @param ionW +// * @param exIonRatio +// * @return +// */ +// public fun getSingleScatterFunction( +// exPos: Double, +// ionPos: Double, +// exW: Double, +// ionW: Double, +// exIonRatio: Double, +// ): UnivariateFunction { +// val func: UnivariateFunction = { eps: Double -> +// if (eps <= 0) { +// 0.0 +// } else { +// val z1 = eps - exPos +// val ex = exIonRatio * exp(-2.0 * z1 * z1 / exW / exW) +// +// val z = 4.0 * (eps - ionPos) * (eps - ionPos) +// val ion = 1 / (1 + z / ionW / ionW) +// +// if (eps < exPos) { +// ex +// } else { +// max(ex, ion) +// } +// } +// } +// +// val cutoff = 25.0 +// //calculating lorentz integral analytically +// val tailNorm = (atan((ionPos - cutoff) * 2.0 / ionW) + 0.5 * PI) * ionW / 2.0 +// val norm: Double = integrator.integrate(range = 0.0..cutoff, function = func).value + tailNorm +// return { e -> func(e) / norm } +// } public val exPos: Symbol by symbol @@ -393,15 +389,15 @@ public class NumassTransmission( public val ionW: Symbol by symbol public val exIonRatio: Symbol by symbol - public fun getSingleScatterFunction(set: Map): UnivariateFunction { - val exPos = set.getValue(exPos) - val ionPos = set.getValue(ionPos) - val exW = set.getValue(exW) - val ionW = set.getValue(ionW) - val exIonRatio = set.getValue(exIonRatio) - - return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio) - } +// public fun getSingleScatterFunction(set: Map): UnivariateFunction { +// val exPos = set.getValue(exPos) +// val ionPos = set.getValue(ionPos) +// val exW = set.getValue(exW) +// val ionW = set.getValue(ionW) +// val exIonRatio = set.getValue(exIonRatio) +// +// return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio) +// } public val trapFunction: (Double, Double) -> Double = { Ei: Double, Ef: Double -> val eps = Ei - Ef diff --git a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt index 1358c1a..fda66d1 100644 --- a/numass-model/src/commonMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt +++ b/numass-model/src/commonMain/kotlin/ru/inr/mass/models/SterileNeutrinoSpectrum.kt @@ -137,7 +137,11 @@ public class SterileNeutrinoSpectrum( arguments: Map, ): Double = DoubleField.gaussIntegrator.integrate(u..eIn, generateRanges( u..eIn, - *((u + 25)..(u + 6000) step 25.0).toDoubleArray() + u + 2.0, + u + 7.0, + u + 15.0, + u + 30.0, + *((u + 50)..(u + 6000) step 30.0).toDoubleArray() )) { eOut: Double -> transFunc(eIn, eOut, arguments) * resolution(eOut, u, arguments) }.value 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 index afd21c2..683c279 100644 --- a/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt +++ b/numass-workspace/src/main/kotlin/ru/inr/mass/scripts/sterileSpectrum.kt @@ -35,34 +35,35 @@ fun main() { ) Plotly.page { + + plot { + scatter { + name = "Computed spectrum" + mode = ScatterMode.lines + x.buffer = 14000.0..18600.0 step 10.0 + y.numbers = x.doubles.map { spectrum(it, args) } + } + layout { + title = "Sterile neutrino spectrum" + yaxis.type = AxisType.log + } + } + val spectrumTime = measureTimeMillis { plot { scatter { - name = "Computed spectrum" - mode = ScatterMode.lines - x.buffer = 14000.0..18600.0 step 10.0 - y.numbers = x.doubles.map { spectrum(it, args) } + mode = ScatterMode.markers + javaClass.getResource("/old-spectrum.dat").readText().lines().map { + val (u, w) = it.split("\t").map { it.toDouble() } + appendXY(u, w / spectrum(u, args) - 1.0) + } } layout { - title = "Sterile neutrino spectrum" - yaxis.type = AxisType.log + title = "Sterile neutrino old/new ratio" } } } println("Spectrum with 460 points computed in $spectrumTime millis") - - plot { - scatter { - mode = ScatterMode.markers - javaClass.getResource("/old-spectrum.dat").readText().lines().map { - val (u, w) = it.split("\t").map { it.toDouble() } - appendXY(u, w / spectrum(u, args) - 1.0) - } - } - layout { - title = "Sterile neutrino old/new ratio" - } - } plot { val resolution = NumassResolution() scatter {