Add simulation comparison

build.gradle.kts

@@ -14,7 +14,7 @@ allprojects {
 
 val dataforgeVersion by extra("0.5.2")
 val tablesVersion: String by extra("0.1.2")
-val kmathVersion by extra("0.3.0-dev-17")
+val kmathVersion by extra("0.3.0-dev-18")
 val plotlyVersion: String by extra("0.5.0")
 
 ksciencePublish{

gradle/wrapper/gradle-wrapper.properties

@@ -1,5 +1,5 @@
 distributionBase=GRADLE_USER_HOME
 distributionPath=wrapper/dists
-distributionUrl=https\://services.gradle.org/distributions/gradle-7.3-bin.zip
+distributionUrl=https\://services.gradle.org/distributions/gradle-7.4.1-bin.zip
 zipStoreBase=GRADLE_USER_HOME
 zipStorePath=wrapper/dists

numass-analysis/src/commonMain/kotlin/ru/inr/mass/data/analysis/NumassAmplitudeSpectrum.kt

@@ -55,6 +55,14 @@ public suspend fun NumassBlock.energySpectrum(
     return map.mapValues { it.value.value }
 }
 
+public suspend fun NumassBlock.eventsCount(extractor: NumassEventExtractor = NumassEventExtractor.EVENTS_ONLY): Long {
+    var counter: Long = 0L
+    extractor.extract(this).collect {
+        counter++
+    }
+    return counter
+}
+
 /**
  * Collect events from block in parallel
  */

numass-analysis/src/commonMain/kotlin/ru/inr/mass/data/analysis/NumassEventExtractor.kt

@@ -2,6 +2,7 @@ package ru.inr.mass.data.analysis
 
 import kotlinx.coroutines.flow.Flow
 import kotlinx.coroutines.flow.map
+import kotlinx.coroutines.flow.mapNotNull
 import ru.inr.mass.data.api.NumassBlock
 import ru.inr.mass.data.api.NumassEvent
 
@@ -37,6 +38,41 @@ public fun interface NumassEventExtractor {
                 )
             }
         }
+
+
+        public val TQDC_V2: NumassEventExtractor = NumassEventExtractor { block ->
+            block.frames.mapNotNull { frame ->
+                var max = Short.MIN_VALUE
+                var min = Short.MAX_VALUE
+                var indexOfMax = 0
+
+                // Taking first 8 points as a baseline
+                val baseline = frame.signal.take(8).average()
+
+                frame.signal.forEachIndexed { index, sh: Short ->
+                    if (sh >= max) {
+                        max = sh
+                        indexOfMax = index
+                    }
+                    if (sh <= min) {
+                        min = sh
+                    }
+                }
+
+                /*
+                 * Filtering large negative splashes
+                 */
+                if (baseline - min < 300) {
+                    NumassEvent(
+                        (max - baseline).toInt().toShort(),
+                        frame.timeOffset + frame.tickSize.inWholeNanoseconds * indexOfMax,
+                        block
+                    )
+                } else {
+                    null
+                }
+            }
+        }
     }
 }
 

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/gun_sim.kt

@@ -0,0 +1,63 @@
+package ru.inr.mass.scripts
+
+import ru.inr.mass.workspace.buffer
+import space.kscience.kmath.functions.asFunction
+import space.kscience.kmath.integration.integrate
+import space.kscience.kmath.integration.splineIntegrator
+import space.kscience.kmath.integration.value
+import space.kscience.kmath.interpolation.SplineInterpolator
+import space.kscience.kmath.interpolation.interpolatePolynomials
+import space.kscience.kmath.operations.DoubleField
+import space.kscience.kmath.real.step
+import space.kscience.plotly.Plotly
+import space.kscience.plotly.layout
+import space.kscience.plotly.makeFile
+import space.kscience.plotly.models.AxisType
+import space.kscience.plotly.scatter
+import kotlin.math.PI
+import kotlin.math.exp
+import kotlin.math.pow
+import kotlin.math.sqrt
+
+fun main() {
+    val backScatteringSpectrum: List<Pair<Double, Double>> = {}.javaClass
+        .getResource("/simulation/Gun19_E_back_scatt.dat")!!.readText()
+        .lineSequence().drop(2).mapNotNull {
+            if (it.isBlank()) return@mapNotNull null
+            val (e, p) = it.split('\t')
+            Pair(e.toDouble(), p.toDouble())
+        }.toList()
+
+    val interpolated = SplineInterpolator.double
+        .interpolatePolynomials(backScatteringSpectrum)
+        .asFunction(DoubleField, 0.0)
+
+    val sigma = 0.3
+    val detectorResolution: (Double) -> Double = { x ->
+        1.0 / sqrt(2 * PI) / sigma * exp(-(x / sigma).pow(2) / 2.0)
+    }
+
+    val convoluted: (Double) -> Double = { x ->
+        DoubleField.splineIntegrator.integrate(-2.0..2.0) { y ->
+            detectorResolution(y) * interpolated(x - y)
+        }.value
+    }
+
+    Plotly.plot {
+//        scatter {
+//            name = "simulation"
+//            x.numbers = backScatteringSpectrum.map { 19.0 - it.first }
+//            y.numbers = backScatteringSpectrum.map { it.second }
+//        }
+        scatter {
+            name = "smeared"
+            x.buffer = 0.0..20.0 step 0.1
+            y.numbers = x.doubles.map { convoluted(19.0 - it) * 0.14/0.01 + 0.86 * detectorResolution(it - 19.0) }
+            println(y.doubles.sum()*0.1)//Norm check
+        }
+        layout {
+            yaxis.type = AxisType.log
+        }
+    }.makeFile()
+
+}

numass-workspace/src/main/kotlin/ru/inr/mass/scripts/run_2021_11.kt

@@ -1,50 +1,54 @@
 package ru.inr.mass.scripts
 
-import kotlinx.coroutines.flow.Flow
-import kotlinx.coroutines.flow.map
-import kotlinx.coroutines.flow.toList
 import kotlinx.coroutines.runBlocking
 import kotlinx.html.h2
 import kotlinx.html.p
 import kotlinx.serialization.json.Json
+import ru.inr.mass.data.analysis.NumassEventExtractor
+import ru.inr.mass.data.analysis.amplitudeSpectrum
 import ru.inr.mass.data.api.NumassFrame
-import ru.inr.mass.workspace.Numass.readDirectory
+import ru.inr.mass.data.api.channels
+import ru.inr.mass.workspace.Numass.readPoint
 import ru.inr.mass.workspace.listFrames
 import space.kscience.dataforge.meta.MetaSerializer
 import space.kscience.plotly.*
 
-fun NumassFrame.tqdcAmplitude(): Short {
+//fun NumassFrame.tqdcAmplitude(): Short {
-    var max = Short.MIN_VALUE
+//    var max = Short.MIN_VALUE
-    var min = Short.MAX_VALUE
+//    var min = Short.MAX_VALUE
+//
+//    signal.forEach { sh: Short ->
+//        if (sh >= max) {
+//            max = sh
+//        }
+//        if (sh <= min) {
+//            min = sh
+//        }
+//    }
+//
+//    return (max - min).toShort()
+//}
 
-    signal.forEach { sh: Short ->
+//fun Flow<NumassFrame>.tqdcAmplitudes(): List<Short> = runBlocking {
-        if (sh >= max) {
+//    map { it.tqdcAmplitude() }.toList()
-            max = sh
+//}
-        }
-        if (sh <= min) {
-            min = sh
-        }
-    }
 
-    return (max - min).toShort()
+val IntRange.center: Double get() = (endInclusive + start).toDouble() / 2.0
-}
-
-fun Flow<NumassFrame>.tqdcAmplitudes(): List<Short> = runBlocking {
-    map { it.tqdcAmplitude() }.toList()
-}
 
 suspend fun main() {
     //val repo: DataTree<NumassDirectorySet> = readNumassRepository("D:\\Work\\numass-data\\")
-    val directory = readDirectory("D:\\Work\\Numass\\data\\test\\set_7")
+    //val directory = readDirectory("D:\\Work\\Numass\\data\\2021_11\\Tritium_2\\set_11\\")
-    val point = directory.points.first()
+    val point = readPoint("D:\\Work\\Numass\\data\\2021_11\\Tritium_2\\set_11\\p0(30s)(HV1=14000)")
+    val channel = point.channels[4]!!
 
-    val frames: List<NumassFrame> = point.listFrames()
+    val binning = 16U
+
+    val frames: List<NumassFrame> = channel.listFrames()
     Plotly.page {
         p { +"${frames.size} frames" }
         h2 { +"Random frames" }
         plot {
             val random = kotlin.random.Random(1234)
-
             repeat(10) {
                 val frame = frames.random(random)
                 scatter {
@@ -54,17 +58,31 @@ suspend fun main() {
         }
         h2 { +"Analysis" }
         plot {
-            histogram {
+            scatter {
                 name = "max"
-                x.numbers = frames.map { frame -> frame.signal.maxOrNull() ?: 0 }
+                val spectrum = runBlocking {
+                    channel.amplitudeSpectrum(NumassEventExtractor.EVENTS_ONLY)
+                }.binned(binning)
+                x.numbers = spectrum.keys.map { it.center }
+                y.numbers = spectrum.values.map { it }
             }
 
-            histogram {
+            scatter {
                 name = "max-min"
-                xbins {
+                val spectrum = runBlocking {
-                    size = 2.0
+                    channel.amplitudeSpectrum(NumassEventExtractor.TQDC)
+                }.binned(binning)
+                x.numbers = spectrum.keys.map { it.center }
+                y.numbers = spectrum.values.map { it }
             }
-                x.numbers = frames.map { it.tqdcAmplitude() }
+
+            scatter {
+                name = "max-baseline + filter"
+                val spectrum = runBlocking {
+                    channel.amplitudeSpectrum(NumassEventExtractor.TQDC_V2)
+                }.binned(binning)
+                x.numbers = spectrum.keys.map { it.center }
+                y.numbers = spectrum.values.map { it }
             }
         }
         h2 { +"Meta" }