Remaking descriptors

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 {

numass-control/cryotemp/src/main/kotlin/inr/numass/control/cryotemp/PKT8Display.kt

@@ -146,7 +146,7 @@ class PKT8Display : DeviceDisplayFX<PKT8Device>(), 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")
             }
         }

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)),
+        return ValueMap.of(
-                x, count, countErr, t)
+                arrayOf(getComponentName(X_VALUE_KEY), getComponentName(Y_VALUE_KEY), getComponentName(Y_ERROR_KEY), getComponentName(POINT_LENGTH_NAME)),
+                x, count, countErr, t
+        )
     }
 
 

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<Int, UnivariateFunction>().also {
+        it[1] = singleScatterFunction
+    }
+
+
+    fun getGunLossProbabilities(X: Double): List<Double> {
+        val res = ArrayList<Double>()
+        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<Double> {
+        val res = ArrayList<Double>()
+        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<Double> {
+        return (lossProbCache).getOrPut(x) {
+            val res = ArrayList<Double>()
+            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<Double>, 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<Double, List<Double>>(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) })
+        }
+
+    }
+
+}

numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt

@@ -66,7 +66,7 @@ fun main(args: Array<String>) {
     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)
         }
 

numass-main/src/main/kotlin/inr/numass/scripts/InversedChainProto.kt

@@ -47,7 +47,7 @@ fun main(args: Array<String>) {
     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)
     }
 

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<String>) {
     }
 
     val frame = displayChart("differential").apply {
-        this.plots.descriptor = Descriptors.forObject(DataPlot::class)
+        this.plots.descriptor = Descriptors.forType(DataPlot::class)
         this.plots.configureValue("showLine", true)
     }
 

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<String>) {
         dataDir = "D:\\Work\\Numass\\data\\2018_04"
     }
 
-    val sp = SterileNeutrinoSpectrum(context, Meta.empty())
+    val spectrum = NBkgSpectrum(SterileNeutrinoSpectrum(context, Meta.empty()))
-    //beta.setCaching(false);
 
-    val spectrum = NBkgSpectrum(sp)
+    val t = 30 * 50 // time in seconds per point
-    //val model = XYModel(Meta.empty(), SpectrumAdapter(Meta.empty()), spectrum)
 
     val params = ParamSet().apply {
         setPar("N", 8e5, 6.0, 0.0, Double.POSITIVE_INFINITY)
@@ -57,9 +65,9 @@ fun main(args: Array<String>) {
         setPar("trap", 1.0, 0.01)
     }
 
-    fun ParamSet.update(vararg override: Pair<String, Double>): ParamSet = this.copy().apply {
+    fun ParamSet.update(vararg override: Pair<String, Double>): 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<String>) {
     }
 
     fun plotResidual(name: String, vararg override: Pair<String, Double>): 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 mod = spectrum.value(it, paramsMod)
-            val err = sqrt(base/1e6)
+            val err = sqrt(base / t)
-            return@map (mod - base) / err
+            (mod - base) / err
+        }
+        return DataPlot.plot(name, x.toDoubleArray(), y.toDoubleArray())
+    }
+
+    val adapter = SpectrumAdapter(Meta.empty())
+    val fm = context.get<FitManager>()
+
+    fun plotFitResidual(name: String, vararg override: Pair<String, Double>): 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<String>) {
             "showErrors" to false
         }
         plots.setType<DataPlot>()
-        +plotResidual("sterile_1","U2" to 1e-3)
+        +plotResidual("sterile_1", "U2" to 1e-3)
-        +plotResidual("sterile_3","msterile2" to (3000*3000).toDouble(),"U2" to 1e-3)
+        +plotResidual("sterile_3", "msterile2" to (3000 * 3000).toDouble(), "U2" to 1e-3)
-        +plotResidual("X","X" to 0.11)
+        +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<DataPlot>()
+        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)
+        }
+
     }
 
 }