Signal processor working.

numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassBlock.kt

@@ -28,6 +28,12 @@ open class OrphanNumassEvent(val amplitude: Short, val timeOffset: Long) : Seria
     override fun compareTo(other: OrphanNumassEvent): Int {
         return this.timeOffset.compareTo(other.timeOffset)
     }
+
+    override fun toString(): String {
+        return "[$amplitude, $timeOffset]"
+    }
+
+
 }
 
 /**

numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/SignalProcessor.kt

@@ -7,5 +7,5 @@ import java.util.stream.Stream
  * Created by darksnake on 07.07.2017.
  */
 interface SignalProcessor {
-    fun analyze(parent: NumassBlock, frame: NumassFrame): Stream<NumassEvent>
+    fun process(parent: NumassBlock, frame: NumassFrame): Stream<NumassEvent>
 }

numass-core/numass-signal-processing/build.gradle.kts

@@ -11,6 +11,7 @@ repositories {
 
 dependencies {
     compile(kotlin("stdlib-jdk8"))
+    compile("hep.dataforge:dataforge-maths")
     compile(project(":numass-core:numass-data-api"))
 
     // https://mvnrepository.com/artifact/org.apache.commons/commons-collections4

numass-core/numass-signal-processing/src/main/kotlin/inr/numass/data/ChernovProcessor.kt

@@ -1,11 +1,10 @@
 package inr.numass.data
 
-import hep.dataforge.meta.Meta
+import inr.numass.data.api.*
-import inr.numass.data.api.NumassBlock
-import inr.numass.data.api.NumassEvent
-import inr.numass.data.api.NumassFrame
-import inr.numass.data.api.SignalProcessor
 import org.apache.commons.collections4.queue.CircularFifoQueue
+import org.apache.commons.math3.fitting.PolynomialCurveFitter
+import org.apache.commons.math3.fitting.WeightedObservedPoint
+import org.slf4j.LoggerFactory
 import java.nio.ShortBuffer
 import java.util.stream.Stream
 import kotlin.streams.asStream
@@ -21,46 +20,81 @@ private fun ShortBuffer.clone(): ShortBuffer {
 }
 
 
-class ChernovProcessor(val meta: Meta) : SignalProcessor {
+class ChernovProcessor(
-    val threshold = meta.getValue("threshold").number.toShort()
+    val threshold: Short,
-    val signalRange: IntRange = TODO()
+    val signalRange: IntRange,
-    val signal: (Double) -> Double = { TODO() }
+    val tickSize: Int = 320,
-    val tickSize: Int = TODO()
+    val signal: (Double) -> Double
+) : SignalProcessor {
 
+    private val fitter = PolynomialCurveFitter.create(2)
+
+    private val signalMax = signal(0.0)
+
+    /**
+     * position an amplitude of peak relative to buffer end (negative)
+     */
     private fun CircularFifoQueue<Short>.findMax(): Pair<Double, Double> {
-        TODO()
+        val data = this.mapIndexed { index, value ->
+            WeightedObservedPoint(
+                1.0,
+                index.toDouble() - size + 1, // final point in zero
+                value.toDouble()
+            )
+        }
+        val (c, b, a) = fitter.fit(data)
+        if (a > 0) error("Minimum!")
+        val x = -b / 2 / a
+        val y = -(b * b - 4 * a * c) / 4 / a
+        return x to y
     }
 
-    override fun analyze(parent: NumassBlock, frame: NumassFrame): Stream<NumassEvent> {
+    fun processBuffer(buffer: ShortBuffer): Sequence<OrphanNumassEvent> {
-        return sequence<NumassEvent> {
-            val events = HashMap<Double, Double>()
-            val buffer = frame.signal.clone()
 
-            val ringBuffer = CircularFifoQueue<Short>(5)
+        val ringBuffer = CircularFifoQueue<Short>(5)
-            while (buffer.remaining() > 0) {
-                ringBuffer.add(buffer.get())
-                val lastValue = ringBuffer[1] ?: -1
-                val currentValue = ringBuffer[0]
-                if (lastValue > threshold && currentValue < lastValue) {
-                    //Found bending, evaluating event
 
-                    ringBuffer.add(buffer.get())//do another step to have 5-points
+        fun roll() {
-                    //TODO check end of frame
+            ringBuffer.add(buffer.get())
-                    val (pos, amp) = ringBuffer.findMax()
+        }
-                    val event = NumassEvent(amp.toShort(), pos.toLong() * tickSize, parent)
-                    yield(event)
 
-                    //subtracting event from buffer copy
+        return sequence<OrphanNumassEvent> {
-                    for (x in signalRange) {
+            while (buffer.remaining() > 1) {
-                        //TODO check all roundings
+                roll()
-                        val position = buffer.position() - x.toShort()
+                if (ringBuffer.isAtFullCapacity) {
-                        val oldValue = buffer.get(position)
+                    if (ringBuffer.all { it > threshold && it <= ringBuffer[2] }) {
-                        val newValue = oldValue - amp * signal(x.toDouble())
+                        //Found bending, evaluating event
-                        buffer.put(position, newValue.toShort())
+                        //TODO check end of frame
+                        try {
+                            val (pos, amp) = ringBuffer.findMax()
+
+                            val timeInTicks = (pos + buffer.position() - 1)
+
+                            val event = OrphanNumassEvent(amp.toShort(), (timeInTicks * tickSize).toLong())
+                            yield(event)
+
+                            //subtracting event from buffer copy
+                            for (x in (signalRange.first + timeInTicks.toInt())..(signalRange.endInclusive + timeInTicks.toInt())) {
+                                //TODO check all roundings
+                                if (x >= 0 && x < buffer.limit()) {
+                                    val oldValue = buffer.get(x)
+                                    val newValue = oldValue - amp * signal(x - timeInTicks) / signalMax
+                                    buffer.put(x, newValue.toShort())
+                                }
+                            }
+                            println(buffer.array().joinToString())
+                        } catch (ex: Exception) {
+                            LoggerFactory.getLogger(javaClass).error("Something went wrong", ex)
+                        }
+                        roll()
                     }
                 }
             }
-        }.asStream()
+        }
+    }
+
+    override fun process(parent: NumassBlock, frame: NumassFrame): Stream<NumassEvent> {
+        val buffer = frame.signal.clone()
+        return processBuffer(buffer).map { it.adopt(parent) }.asStream()
     }
 }
 

numass-core/numass-signal-processing/src/test/kotlin/inr/numass/data/ChernovProcessorTest.kt

@@ -0,0 +1,26 @@
+package inr.numass.data
+
+import org.apache.commons.math3.analysis.function.Gaussian
+import org.junit.Assert.assertTrue
+import org.junit.Test
+import java.nio.ShortBuffer
+
+class ChernovProcessorTest {
+    val gaussian = Gaussian(1000.0, 0.0, 3.0)
+    val processor = ChernovProcessor(10, -12..12, tickSize = 100) { gaussian.value(it) }
+
+    val events = mapOf<Double, Double>(10.0 to 1.0, 16.0 to 0.5)
+
+    val buffer = ShortArray(40) { i ->
+        events.entries.sumByDouble { (pos, amp) -> amp * gaussian.value(pos - i.toDouble()) }.toShort()
+    }
+
+    @Test
+    fun testPeaks() {
+        println(buffer.joinToString())
+        val peaks = processor.processBuffer(ShortBuffer.wrap(buffer)).toList()
+        assertTrue(peaks.isNotEmpty())
+        println(peaks.joinToString())
+    }
+
+}

numass-core/src/main/kotlin/inr/numass/data/analyzers/AbstractAnalyzer.kt

@@ -27,13 +27,13 @@ import inr.numass.data.api.NumassEvent
 import inr.numass.data.api.NumassPoint.Companion.HV_KEY
 import inr.numass.data.api.NumassSet
 import inr.numass.data.api.SignalProcessor
-import java.lang.IllegalArgumentException
 import java.util.stream.Stream
 
 /**
  * Created by darksnake on 11.07.2017.
  */
-abstract class AbstractAnalyzer @JvmOverloads constructor(private val processor: SignalProcessor? = null) : NumassAnalyzer {
+abstract class AbstractAnalyzer @JvmOverloads constructor(private val processor: SignalProcessor? = null) :
+    NumassAnalyzer {
 
     /**
      * Return unsorted stream of events including events from frames.
@@ -58,7 +58,7 @@ abstract class AbstractAnalyzer @JvmOverloads constructor(private val processor:
         return when {
             block.frames.count() == 0L -> block.events
             processor == null -> throw IllegalArgumentException("Signal processor needed to analyze frames")
-            else -> Stream.concat(block.events, block.frames.flatMap { processor.analyze(it) })
+            else -> Stream.concat(block.events, block.frames.flatMap { processor.process(block, it) })
         }
     }
 
@@ -70,14 +70,14 @@ abstract class AbstractAnalyzer @JvmOverloads constructor(private val processor:
      */
     protected open fun getTableFormat(config: Meta): TableFormat {
         return TableFormatBuilder()
-                .addNumber(HV_KEY, X_VALUE_KEY)
+            .addNumber(HV_KEY, X_VALUE_KEY)
-                .addNumber(NumassAnalyzer.LENGTH_KEY)
+            .addNumber(NumassAnalyzer.LENGTH_KEY)
-                .addNumber(NumassAnalyzer.COUNT_KEY)
+            .addNumber(NumassAnalyzer.COUNT_KEY)
-                .addNumber(NumassAnalyzer.COUNT_RATE_KEY, Y_VALUE_KEY)
+            .addNumber(NumassAnalyzer.COUNT_RATE_KEY, Y_VALUE_KEY)
-                .addNumber(NumassAnalyzer.COUNT_RATE_ERROR_KEY, Y_ERROR_KEY)
+            .addNumber(NumassAnalyzer.COUNT_RATE_ERROR_KEY, Y_ERROR_KEY)
-                .addColumn(NumassAnalyzer.WINDOW_KEY)
+            .addColumn(NumassAnalyzer.WINDOW_KEY)
-                .addTime()
+            .addTime()
-                .build()
+            .build()
     }
 
 
@@ -85,18 +85,18 @@ abstract class AbstractAnalyzer @JvmOverloads constructor(private val processor:
         val format = getTableFormat(config)
 
         return ListTable.Builder(format)
-                .rows(set.points.map { point -> analyzeParent(point, config) })
+            .rows(set.points.map { point -> analyzeParent(point, config) })
-                .build()
+            .build()
     }
 
     companion object {
         val NAME_LIST = arrayOf(
-                NumassAnalyzer.LENGTH_KEY,
+            NumassAnalyzer.LENGTH_KEY,
-                NumassAnalyzer.COUNT_KEY,
+            NumassAnalyzer.COUNT_KEY,
-                NumassAnalyzer.COUNT_RATE_KEY,
+            NumassAnalyzer.COUNT_RATE_KEY,
-                NumassAnalyzer.COUNT_RATE_ERROR_KEY,
+            NumassAnalyzer.COUNT_RATE_ERROR_KEY,
-                NumassAnalyzer.WINDOW_KEY,
+            NumassAnalyzer.WINDOW_KEY,
-                NumassAnalyzer.TIME_KEY
+            NumassAnalyzer.TIME_KEY
         )
     }
 }

numass-main/build.gradle

@@ -36,6 +36,7 @@ dependencies {
     compile group: 'commons-cli', name: 'commons-cli', version: '1.+'
     compile group: 'commons-io', name: 'commons-io', version: '2.+'
     compile project(':numass-core')
+    compileOnly "org.jetbrains.kotlin:kotlin-main-kts:1.3.21"
     compile "hep.dataforge:dataforge-minuit" //project(':dataforge-stat:dataforge-minuit')
     compile "hep.dataforge:grind-terminal" //project(':dataforge-grind:grind-terminal')
     compile "hep.dataforge:dataforge-gui"

numass-main/src/main/kotlin/inr/numass/scripts/threshold/threshold-2017_05-CAMAC.kt

@@ -0,0 +1,71 @@
+/*
+ * 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.
+ */
+
+
+import hep.dataforge.buildContext
+import hep.dataforge.fx.output.FXOutputManager
+import hep.dataforge.io.DirectoryOutput
+import hep.dataforge.io.plus
+import hep.dataforge.meta.buildMeta
+import hep.dataforge.nullable
+import hep.dataforge.plots.data.DataPlot
+import hep.dataforge.plots.jfreechart.JFreeChartPlugin
+import hep.dataforge.plots.plotData
+import hep.dataforge.storage.files.FileStorage
+import hep.dataforge.tables.Adapters
+import inr.numass.NumassPlugin
+import inr.numass.data.api.NumassSet
+import inr.numass.data.storage.NumassDirectory
+import inr.numass.displayChart
+import inr.numass.subthreshold.Threshold
+
+fun main(){
+    val context = buildContext("NUMASS", NumassPlugin::class.java, JFreeChartPlugin::class.java) {
+        rootDir = "D:\\Work\\Numass\\sterile\\2017_05"
+        dataDir = "D:\\Work\\Numass\\data\\2017_05"
+        output = FXOutputManager() + DirectoryOutput()
+    }
+
+    val storage = NumassDirectory.read(context, "Fill_3") as? FileStorage ?: error("Storage not found")
+
+    val meta = buildMeta {
+        "delta" to -200
+        "method" to "pow"
+        "t0" to 15e3
+//        "window.lo" to 400
+//        "window.up" to 1600
+        "xLow" to 450
+        "xHigh" to 700
+        "upper" to 3000
+        "binning" to 20
+        //"reference" to 18600
+    }
+
+    val frame = displayChart("correction").apply {
+        plots.setType<DataPlot>()
+    }
+
+    listOf("set_2", "set_3", "set_4", "set_5").forEach { setName ->
+        val set = storage.provide(setName, NumassSet::class.java).nullable ?: error("Set does not exist")
+
+        val correctionTable = Threshold.calculateSubThreshold(set, meta).filter {
+            it.getDouble("correction") in (1.0..1.2)
+        }
+
+        frame.plotData(setName, correctionTable, Adapters.buildXYAdapter("U", "correction"))
+    }
+
+}