diff --git a/build.gradle b/build.gradle index 82c3959c..47732cab 100644 --- a/build.gradle +++ b/build.gradle @@ -28,6 +28,7 @@ allprojects { //maven { url = "https://jitpack.io" } maven { url = "http://dl.bintray.com/kotlin/ktor" } maven { url = "https://dl.bintray.com/kotlin/kotlinx" } + maven { url = "https://dl.bintray.com/kotlin/kotlin-eap"} } dependencies { diff --git a/numass-control/cryotemp/src/main/kotlin/inr/numass/control/cryotemp/PKT8Display.kt b/numass-control/cryotemp/src/main/kotlin/inr/numass/control/cryotemp/PKT8Display.kt index 3dac7178..95ec7c73 100644 --- a/numass-control/cryotemp/src/main/kotlin/inr/numass/control/cryotemp/PKT8Display.kt +++ b/numass-control/cryotemp/src/main/kotlin/inr/numass/control/cryotemp/PKT8Display.kt @@ -146,7 +146,7 @@ class PKT8Display : DeviceDisplayFX(), PKT8ValueListener { private val plotFrame: PlotFrame by lazy { JFreeChartFrame(plotFrameMeta).apply { - plots.descriptor = Descriptors.forElement(TimePlot::class.java) + plots.descriptor = Descriptors.forType(TimePlot::class.java) PlotUtils.setXAxis(this, "timestamp", "", "time") } } diff --git a/numass-core/src/main/kotlin/inr/numass/data/SpectrumAdapter.kt b/numass-core/src/main/kotlin/inr/numass/data/SpectrumAdapter.kt index c1f7a92a..416ecdf2 100644 --- a/numass-core/src/main/kotlin/inr/numass/data/SpectrumAdapter.kt +++ b/numass-core/src/main/kotlin/inr/numass/data/SpectrumAdapter.kt @@ -61,8 +61,10 @@ class SpectrumAdapter : BasicAdapter { } fun buildSpectrumDataPoint(x: Double, count: Long, countErr: Double, t: Double): Values { - return ValueMap.of(arrayOf(getComponentName(X_VALUE_KEY), getComponentName(Y_VALUE_KEY), getComponentName(Y_ERROR_KEY), getComponentName(POINT_LENGTH_NAME)), - x, count, countErr, t) + return ValueMap.of( + arrayOf(getComponentName(X_VALUE_KEY), getComponentName(Y_VALUE_KEY), getComponentName(Y_ERROR_KEY), getComponentName(POINT_LENGTH_NAME)), + x, count, countErr, t + ) } diff --git a/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt b/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt new file mode 100644 index 00000000..fd732819 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/models/misc/LossCalculator.kt @@ -0,0 +1,417 @@ +/* + * Copyright 2018 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 inr.numass.models.misc + +import hep.dataforge.maths.functions.FunctionCaching +import hep.dataforge.maths.integration.GaussRuleIntegrator +import hep.dataforge.plots.PlotFrame +import hep.dataforge.plots.XYFunctionPlot +import hep.dataforge.utils.Misc +import hep.dataforge.values.Values +import org.apache.commons.math3.analysis.BivariateFunction +import org.apache.commons.math3.analysis.UnivariateFunction +import org.apache.commons.math3.exception.OutOfRangeException +import org.slf4j.LoggerFactory +import java.lang.Math.exp +import java.util.* + +/** + * Вычисление произвольного порядка функции рассеяния. Не учитывается + * зависимость сечения от энергии электрона + * + * @author Darksnake + */ +object LossCalculator { + private val cache = HashMap().also { + it[1] = singleScatterFunction + } + + + fun getGunLossProbabilities(X: Double): List { + val res = ArrayList() + var prob: Double + if (X > 0) { + prob = Math.exp(-X) + } else { + // если x ==0, то выживает только нулевой член, первый равен 1 + res.add(1.0) + return res + } + res.add(prob) + + var n = 0 + while (prob > SCATTERING_PROBABILITY_THRESHOLD) { + /* + * prob(n) = prob(n-1)*X/n; + */ + n++ + prob *= X / n + res.add(prob) + } + + return res + } + + fun getGunZeroLossProb(X: Double): Double { + return Math.exp(-X) + } + + /** + * Ленивое рекурсивное вычисление функции потерь через предыдущие + * + * @param order + * @return + */ + private fun getLoss(order: Int): UnivariateFunction? { + if (order <= 0) { + throw IllegalArgumentException() + } + return if (cache.containsKey(order)) { + cache[order] + } else { + synchronized(this) { + cache.getOrPut(order) { + LoggerFactory.getLogger(javaClass) + .debug("Scatter cache of order {} not found. Updating", order) + getNextLoss(getMargin(order), getLoss(order - 1)) + } + return cache[order] + } + } + } + + fun getLossFunction(order: Int): BivariateFunction { + assert(order > 0) + return BivariateFunction { Ei: Double, Ef: Double -> getLossValue(order, Ei, Ef) } + } + + fun getLossProbDerivs(x: Double): List { + val res = ArrayList() + val probs = getLossProbabilities(x) + + var delta = Math.exp(-x) + res.add((delta - probs[0]) / x) + for (i in 1 until probs.size) { + delta *= x / i + res.add((delta - probs[i]) / x) + } + + return res + } + + /** + * рекурсивно вычисляем все вероятности, котрорые выше порога + * + * + * дисер, стр.48 + * + * @param X + * @return + */ + fun getLossProbabilities(x: Double): List { + return (lossProbCache).getOrPut(x) { + val res = ArrayList() + var prob: Double + if (x > 0) { + prob = 1 / x * (1 - Math.exp(-x)) + } else { + // если x ==0, то выживает только нулевой член, первый равен нулю + res.add(1.0) + return@getOrPut res + } + res.add(prob) + + while (prob > SCATTERING_PROBABILITY_THRESHOLD) { + /* + * prob(n) = prob(n-1)-1/n! * X^n * exp(-X); + */ + var delta = Math.exp(-x) + for (i in 1 until res.size + 1) { + delta *= x / i + } + prob -= delta / x + res.add(prob) + } + + res + } + } + + fun getLossProbability(order: Int, X: Double): Double { + if (order == 0) { + return if (X > 0) { + 1 / X * (1 - Math.exp(-X)) + } else { + 1.0 + } + } + val probs = getLossProbabilities(X) + return if (order >= probs.size) { + 0.0 + } else { + probs[order] + } + } + + fun getLossValue(order: Int, Ei: Double, Ef: Double): Double { + return if (Ei - Ef < 5.0) { + 0.0 + } else if (Ei - Ef >= getMargin(order)) { + 0.0 + } else { + getLoss(order)!!.value(Ei - Ef) + } + } + + /** + * функция потерь с произвольными вероятностями рассеяния + * + * @param probs + * @param Ei + * @param Ef + * @return + */ + 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) + } + return sum + } + + /** + * граница интегрирования + * + * @param order + * @return + */ + private fun getMargin(order: Int): Double { + return 80 + order * 50.0 + } + + /** + * генерирует кэшированную функцию свертки loss со спектром однократных + * потерь + * + * @param loss + * @return + */ + private fun getNextLoss(margin: Double, loss: UnivariateFunction?): UnivariateFunction { + val res = { x: Double -> + val integrand = UnivariateFunction { y: Double -> + try { + loss!!.value(x - y) * singleScatterFunction.value(y) + } catch (ex: OutOfRangeException) { + 0.0 + } + } + integrator.integrate(5.0, margin, integrand) + } + + return FunctionCaching.cacheUnivariateFunction(0.0, margin, 200, res) + + } + + fun getTotalLossBivariateFunction(X: Double): BivariateFunction { + return BivariateFunction { Ei: Double, Ef: Double -> getTotalLossValue(X, Ei, Ef) } + } + + /** + * Значение полной производной функции потерь с учетом всех неисчезающих + * порядков + * + * @param X + * @param Ei + * @param Ef + * @return + */ + fun getTotalLossDeriv(X: Double, Ei: Double, Ef: Double): Double { + val probs = getLossProbDerivs(X) + + var sum = 0.0 + for (i in 1 until probs.size) { + sum += probs[i] * getLossValue(i, Ei, Ef) + } + return sum + } + + fun getTotalLossDerivBivariateFunction(X: Double): BivariateFunction { + return BivariateFunction { Ei: Double, Ef: Double -> getTotalLossDeriv(X, Ei, Ef) } + } + + /** + * Значение полной функции потерь с учетом всех неисчезающих порядков + * + * @param x + * @param Ei + * @param Ef + * @return + */ + fun getTotalLossValue(x: Double, Ei: Double, Ef: Double): Double { + if (x == 0.0) { + return 0.0 + } else { + val probs = getLossProbabilities(x) + return (0 until probs.size).sumByDouble { i -> + probs[i] * getLossValue(i, Ei, Ef) + } + } + } + + + /** + * порог по вероятности, до которого вычисляются компоненты функции потерь + */ + private const val SCATTERING_PROBABILITY_THRESHOLD = 1e-3 + private val integrator = GaussRuleIntegrator(100) + private val lossProbCache = Misc.getLRUCache>(100) + + + private val A1 = 0.204 + private val A2 = 0.0556 + private val b = 14.0 + private val pos1 = 12.6 + private val pos2 = 14.3 + private val w1 = 1.85 + private val w2 = 12.5 + + val singleScatterFunction = UnivariateFunction { eps: Double -> + when { + eps <= 0 -> 0.0 + eps <= b -> { + val z = eps - pos1 + A1 * exp(-2.0 * z * z / w1 / w1) + } + else -> { + val z = 4.0 * (eps - pos2) * (eps - pos2) + A2 / (1 + z / w2 / w2) + } + } + } + + + /** + * A generic loss function for numass experiment in "Lobashev" + * parameterization + * + * @param exPos + * @param ionPos + * @param exW + * @param ionW + * @param exIonRatio + * @return + */ + 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 { + Math.max(ex, ion) + } + } + } + + val cutoff = 25.0 + //caclulating lorentz integral analythically + val tailNorm = (Math.atan((ionPos - cutoff) * 2.0 / ionW) + 0.5 * Math.PI) * ionW / 2.0 + val norm = integrator.integrate(0.0, cutoff, func)!! + tailNorm + return UnivariateFunction { e -> func.value(e) / norm } + } + + fun getSingleScatterFunction(set: Values): UnivariateFunction { + + val exPos = set.getDouble("exPos") + val ionPos = set.getDouble("ionPos") + val exW = set.getDouble("exW") + val ionW = set.getDouble("ionW") + val exIonRatio = set.getDouble("exIonRatio") + + return getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio) + } + + val trapFunction: BivariateFunction + get() = BivariateFunction { Ei: Double, Ef: Double -> + val eps = Ei - Ef + if (eps > 10) { + 1.86e-04 * exp(-eps / 25.0) + 5.5e-05 + } else { + 0.0 + } + } + + fun plotScatter(frame: PlotFrame, set: Values) { + //"X", "shift", "exPos", "ionPos", "exW", "ionW", "exIonRatio" + + // JFreeChartFrame frame = JFreeChartFrame.drawFrame("Differential scattering crosssection", null); + val X = set.getDouble("X") + + val exPos = set.getDouble("exPos") + + val ionPos = set.getDouble("ionPos") + + val exW = set.getDouble("exW") + + val ionW = set.getDouble("ionW") + + val exIonRatio = set.getDouble("exIonRatio") + + val scatterFunction = getSingleScatterFunction(exPos, ionPos, exW, ionW, exIonRatio) + + if (set.names.contains("X")) { + val probs = LossCalculator.getGunLossProbabilities(set.getDouble("X")) + val single = { e: Double -> probs[1] * scatterFunction.value(e) } + frame.add(XYFunctionPlot.plot("Single scattering", 0.0, 100.0, 1000) { x: Double -> single(x) }) + + for (i in 2 until probs.size) { + val scatter = { e: Double -> probs[i] * LossCalculator.getLossValue(i, e, 0.0) } + frame.add(XYFunctionPlot.plot(i.toString() + " scattering", 0.0, 100.0, 1000) { x: Double -> scatter(x) }) + } + + val total = UnivariateFunction { eps -> + if (probs.size == 1) { + return@UnivariateFunction 0.0 + } + var sum = probs[1] * scatterFunction.value(eps) + for (i in 2 until probs.size) { + sum += probs[i] * LossCalculator.getLossValue(i, eps, 0.0) + } + return@UnivariateFunction sum + } + + frame.add(XYFunctionPlot.plot("Total loss", 0.0, 100.0, 1000) { x: Double -> total.value(x) }) + + } else { + + frame.add(XYFunctionPlot.plot("Differential cross-section", 0.0, 100.0, 2000) { x: Double -> scatterFunction.value(x) }) + } + + } + +} diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt b/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt index e77ff2d5..3bbc9dec 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt @@ -66,7 +66,7 @@ fun main(args: Array) { for (hv in arrayOf(14000.0, 14500.0, 15000.0, 15500.0, 16050.0)) { val frame = displayChart("integral[$hv]").apply { - this.plots.descriptor = Descriptors.forObject(DataPlot::class) + this.plots.descriptor = Descriptors.forType(DataPlot::class) this.plots.configureValue("showLine", true) } diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/InversedChainProto.kt b/numass-main/src/main/kotlin/inr/numass/scripts/InversedChainProto.kt index 6ade2fd4..e6fdc90c 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/InversedChainProto.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/InversedChainProto.kt @@ -47,7 +47,7 @@ fun main(args: Array) { val point = ProtoNumassPoint.readFile(Paths.get("D:\\Work\\Numass\\data\\2017_05_frames\\Fill_3_events\\set_33\\p36(30s)(HV1=17000).df")) val frame = displayChart("integral").apply { - this.plots.descriptor = Descriptors.forObject(DataPlot::class) + this.plots.descriptor = Descriptors.forType(DataPlot::class) this.plots.configureValue("showLine", true) } diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt b/numass-main/src/main/kotlin/inr/numass/scripts/models/DifferentialSpectrum.kt similarity index 92% rename from numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt rename to numass-main/src/main/kotlin/inr/numass/scripts/models/DifferentialSpectrum.kt index 4fc97e79..dba51c4d 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/models/DifferentialSpectrum.kt @@ -1,5 +1,5 @@ /* - * Copyright 2017 Alexander Nozik. + * Copyright 2018 Alexander Nozik. * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. @@ -14,7 +14,7 @@ * limitations under the License. */ -package inr.numass.scripts +package inr.numass.scripts.models import hep.dataforge.description.Descriptors import hep.dataforge.kodex.buildContext @@ -60,7 +60,7 @@ fun main(args: Array) { } val frame = displayChart("differential").apply { - this.plots.descriptor = Descriptors.forObject(DataPlot::class) + this.plots.descriptor = Descriptors.forType(DataPlot::class) this.plots.configureValue("showLine", true) } diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt b/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt index 7e7aab33..425268cc 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/models/IntegralSpectrum.kt @@ -17,6 +17,7 @@ package inr.numass.scripts.models import hep.dataforge.fx.output.FXOutputManager +import hep.dataforge.io.output.stream import hep.dataforge.kodex.buildContext import hep.dataforge.kodex.configure import hep.dataforge.kodex.step @@ -25,10 +26,19 @@ import hep.dataforge.plots.Plot import hep.dataforge.plots.data.DataPlot import hep.dataforge.plots.jfreechart.JFreeChartPlugin import hep.dataforge.plots.output.plot +import hep.dataforge.stat.fit.FitManager +import hep.dataforge.stat.fit.FitStage +import hep.dataforge.stat.fit.FitState import hep.dataforge.stat.fit.ParamSet +import hep.dataforge.stat.models.XYModel +import hep.dataforge.tables.Adapters +import hep.dataforge.tables.ListTable import inr.numass.NumassPlugin +import inr.numass.data.SpectrumAdapter import inr.numass.models.NBkgSpectrum import inr.numass.models.sterile.SterileNeutrinoSpectrum +import kotlinx.coroutines.experimental.launch +import java.io.PrintWriter import kotlin.math.sqrt @@ -40,11 +50,9 @@ fun main(args: Array) { dataDir = "D:\\Work\\Numass\\data\\2018_04" } - val sp = SterileNeutrinoSpectrum(context, Meta.empty()) - //beta.setCaching(false); + val spectrum = NBkgSpectrum(SterileNeutrinoSpectrum(context, Meta.empty())) - val spectrum = NBkgSpectrum(sp) - //val model = XYModel(Meta.empty(), SpectrumAdapter(Meta.empty()), spectrum) + val t = 30 * 50 // time in seconds per point val params = ParamSet().apply { setPar("N", 8e5, 6.0, 0.0, Double.POSITIVE_INFINITY) @@ -57,9 +65,9 @@ fun main(args: Array) { setPar("trap", 1.0, 0.01) } - fun ParamSet.update(vararg override: Pair): ParamSet = this.copy().apply { + fun ParamSet.update(vararg override: Pair): ParamSet = this.copy().also { set -> override.forEach { - setParValue(it.first, it.second) + set.setParValue(it.first, it.second) } } @@ -70,19 +78,52 @@ fun main(args: Array) { } fun plotResidual(name: String, vararg override: Pair): Plot { + val paramsMod = params.update(*override) + val x = (14000.0..18600.0).step(100.0).toList() val y = x.map { val base = spectrum.value(it, params) - val mod = spectrum.value(it, params.update(*override)) - val err = sqrt(base/1e6) - return@map (mod - base) / err + val mod = spectrum.value(it, paramsMod) + val err = sqrt(base / t) + (mod - base) / err + } + return DataPlot.plot(name, x.toDoubleArray(), y.toDoubleArray()) + } + + val adapter = SpectrumAdapter(Meta.empty()) + val fm = context.get() + + fun plotFitResidual(name: String, vararg override: Pair): Plot { + val paramsMod = params.update(*override) + + val x = (14000.0..18400.0).step(100.0).toList() + + val table = ListTable.Builder(Adapters.getFormat(adapter)).apply { + x.forEach { u -> + row(adapter.buildSpectrumDataPoint(u, t * spectrum.value(u, params).toLong(), t.toDouble())) + } + }.build() + + val model = XYModel(Meta.empty(), adapter, spectrum) + val state = FitState(table, model, params) + val res = fm.runStage(state, "QOW", FitStage.TASK_RUN, "N", "E0","bkg") + + res.printState(PrintWriter(System.out)) + res.printState(PrintWriter(context.output["fitResult", name].stream)) + //context.output["fitResult",name].stream + + val y = x.map { u -> + val base = spectrum.value(u, params) + val mod = spectrum.value(u, res.parameters) + val err = sqrt(base / t) + (mod - base) / err } return DataPlot.plot(name, x.toDoubleArray(), y.toDoubleArray()) } - - context.plot("default") { +/* + context.plot("trap") { plots.configure { "showLine" to true "showSymbol" to false @@ -100,10 +141,37 @@ fun main(args: Array) { "showErrors" to false } plots.setType() - +plotResidual("sterile_1","U2" to 1e-3) - +plotResidual("sterile_3","msterile2" to (3000*3000).toDouble(),"U2" to 1e-3) - +plotResidual("X","X" to 0.11) + +plotResidual("sterile_1", "U2" to 1e-3) + +plotResidual("sterile_3", "msterile2" to (3000 * 3000).toDouble(), "U2" to 1e-3) + +plotResidual("X", "X" to 0.11) +plotResidual("trap", "trap" to 0.99) + }*/ + + context.plot("fit", stage = "plots") { + plots.configure { + "showLine" to true + "showSymbol" to false + "showErrors" to false + "thickness" to 4.0 + } + plots.setType() + launch { + +plotResidual("trap", "trap" to 0.99) + +plotFitResidual("trap_fit", "trap" to 0.99) + } + launch { + +plotResidual("X", "X" to 0.11) + +plotFitResidual("X_fit", "X" to 0.11) + } + launch { + +plotResidual("sterile_1", "U2" to 1e-3) + +plotFitResidual("sterile_1_fit", "U2" to 1e-3) + } + launch { + +plotResidual("sterile_3", "msterile2" to (3000 * 3000).toDouble(), "U2" to 1e-3) + +plotFitResidual("sterile_3_fit", "msterile2" to (3000 * 3000).toDouble(), "U2" to 1e-3) + } + } } \ No newline at end of file