From 5b170da7d43cf02c99d5a055f347e7f38085724d Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Tue, 19 Feb 2019 16:28:07 +0300 Subject: [PATCH] Signal processor working. --- .../kotlin/inr/numass/data/api/NumassBlock.kt | 6 + .../inr/numass/data/api/SignalProcessor.kt | 2 +- .../numass-signal-processing/build.gradle.kts | 1 + .../inr/numass/data/ChernovProcessor.kt | 104 ++++++++++++------ .../inr/numass/data/ChernovProcessorTest.kt | 26 +++++ .../numass/data/analyzers/AbstractAnalyzer.kt | 38 +++---- numass-main/build.gradle | 1 + .../threshold/threshold-2017_05-CAMAC.kt | 71 ++++++++++++ 8 files changed, 194 insertions(+), 55 deletions(-) create mode 100644 numass-core/numass-signal-processing/src/test/kotlin/inr/numass/data/ChernovProcessorTest.kt create mode 100644 numass-main/src/main/kotlin/inr/numass/scripts/threshold/threshold-2017_05-CAMAC.kt diff --git a/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassBlock.kt b/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassBlock.kt index 751b22af..950b3683 100644 --- a/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassBlock.kt +++ b/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]" + } + + } /** diff --git a/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/SignalProcessor.kt b/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/SignalProcessor.kt index 1928d73a..7689f75a 100644 --- a/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/SignalProcessor.kt +++ b/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 + fun process(parent: NumassBlock, frame: NumassFrame): Stream } diff --git a/numass-core/numass-signal-processing/build.gradle.kts b/numass-core/numass-signal-processing/build.gradle.kts index d2a314b1..8b972e3f 100644 --- a/numass-core/numass-signal-processing/build.gradle.kts +++ b/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 diff --git a/numass-core/numass-signal-processing/src/main/kotlin/inr/numass/data/ChernovProcessor.kt b/numass-core/numass-signal-processing/src/main/kotlin/inr/numass/data/ChernovProcessor.kt index fedc5cd2..a433bb45 100644 --- a/numass-core/numass-signal-processing/src/main/kotlin/inr/numass/data/ChernovProcessor.kt +++ b/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.NumassBlock -import inr.numass.data.api.NumassEvent -import inr.numass.data.api.NumassFrame -import inr.numass.data.api.SignalProcessor +import inr.numass.data.api.* 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 { - val threshold = meta.getValue("threshold").number.toShort() - val signalRange: IntRange = TODO() - val signal: (Double) -> Double = { TODO() } - val tickSize: Int = TODO() +class ChernovProcessor( + val threshold: Short, + val signalRange: IntRange, + val tickSize: Int = 320, + 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.findMax(): Pair { - 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 { - return sequence { - val events = HashMap() - val buffer = frame.signal.clone() + fun processBuffer(buffer: ShortBuffer): Sequence { - val ringBuffer = CircularFifoQueue(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 + val ringBuffer = CircularFifoQueue(5) - ringBuffer.add(buffer.get())//do another step to have 5-points - //TODO check end of frame - val (pos, amp) = ringBuffer.findMax() - val event = NumassEvent(amp.toShort(), pos.toLong() * tickSize, parent) - yield(event) + fun roll() { + ringBuffer.add(buffer.get()) + } - //subtracting event from buffer copy - for (x in signalRange) { - //TODO check all roundings - val position = buffer.position() - x.toShort() - val oldValue = buffer.get(position) - val newValue = oldValue - amp * signal(x.toDouble()) - buffer.put(position, newValue.toShort()) + return sequence { + while (buffer.remaining() > 1) { + roll() + if (ringBuffer.isAtFullCapacity) { + if (ringBuffer.all { it > threshold && it <= ringBuffer[2] }) { + //Found bending, evaluating event + //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 { + val buffer = frame.signal.clone() + return processBuffer(buffer).map { it.adopt(parent) }.asStream() } } diff --git a/numass-core/numass-signal-processing/src/test/kotlin/inr/numass/data/ChernovProcessorTest.kt b/numass-core/numass-signal-processing/src/test/kotlin/inr/numass/data/ChernovProcessorTest.kt new file mode 100644 index 00000000..846808ef --- /dev/null +++ b/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(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()) + } + +} \ No newline at end of file diff --git a/numass-core/src/main/kotlin/inr/numass/data/analyzers/AbstractAnalyzer.kt b/numass-core/src/main/kotlin/inr/numass/data/analyzers/AbstractAnalyzer.kt index c43ccf66..5c53649a 100644 --- a/numass-core/src/main/kotlin/inr/numass/data/analyzers/AbstractAnalyzer.kt +++ b/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(NumassAnalyzer.LENGTH_KEY) - .addNumber(NumassAnalyzer.COUNT_KEY) - .addNumber(NumassAnalyzer.COUNT_RATE_KEY, Y_VALUE_KEY) - .addNumber(NumassAnalyzer.COUNT_RATE_ERROR_KEY, Y_ERROR_KEY) - .addColumn(NumassAnalyzer.WINDOW_KEY) - .addTime() - .build() + .addNumber(HV_KEY, X_VALUE_KEY) + .addNumber(NumassAnalyzer.LENGTH_KEY) + .addNumber(NumassAnalyzer.COUNT_KEY) + .addNumber(NumassAnalyzer.COUNT_RATE_KEY, Y_VALUE_KEY) + .addNumber(NumassAnalyzer.COUNT_RATE_ERROR_KEY, Y_ERROR_KEY) + .addColumn(NumassAnalyzer.WINDOW_KEY) + .addTime() + .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) }) - .build() + .rows(set.points.map { point -> analyzeParent(point, config) }) + .build() } companion object { val NAME_LIST = arrayOf( - NumassAnalyzer.LENGTH_KEY, - NumassAnalyzer.COUNT_KEY, - NumassAnalyzer.COUNT_RATE_KEY, - NumassAnalyzer.COUNT_RATE_ERROR_KEY, - NumassAnalyzer.WINDOW_KEY, - NumassAnalyzer.TIME_KEY + NumassAnalyzer.LENGTH_KEY, + NumassAnalyzer.COUNT_KEY, + NumassAnalyzer.COUNT_RATE_KEY, + NumassAnalyzer.COUNT_RATE_ERROR_KEY, + NumassAnalyzer.WINDOW_KEY, + NumassAnalyzer.TIME_KEY ) } } diff --git a/numass-main/build.gradle b/numass-main/build.gradle index eba39a00..2c075954 100644 --- a/numass-main/build.gradle +++ b/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" diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/threshold/threshold-2017_05-CAMAC.kt b/numass-main/src/main/kotlin/inr/numass/scripts/threshold/threshold-2017_05-CAMAC.kt new file mode 100644 index 00000000..1bd6ef4a --- /dev/null +++ b/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() + } + + 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")) + } + +} \ No newline at end of file