diff --git a/build.gradle b/build.gradle index cbf22cb7..78048f6b 100644 --- a/build.gradle +++ b/build.gradle @@ -9,6 +9,11 @@ buildscript { classpath "org.jetbrains.kotlin:kotlin-gradle-plugin:$kotlin_version" } } + +//plugins{ +// id 'org.openjfx.javafxplugin' version '0.0.7' apply false +//} + allprojects { apply plugin: 'idea' apply plugin: 'java' diff --git a/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassFrame.kt b/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassFrame.kt index e77abd5a..329190ee 100644 --- a/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassFrame.kt +++ b/numass-core/numass-data-api/src/main/kotlin/inr/numass/data/api/NumassFrame.kt @@ -9,18 +9,19 @@ import java.time.Instant * Created by darksnake on 06-Jul-17. */ class NumassFrame( - /** - * The absolute start time of the frame - */ - val time: Instant, - /** - * The time interval per tick - */ - val tickSize: Duration, - /** - * The buffered signal shape in ticks - */ - val signal: ShortBuffer) { + /** + * The absolute start time of the frame + */ + val time: Instant, + /** + * The time interval per tick + */ + val tickSize: Duration, + /** + * The buffered signal shape in ticks + */ + val signal: ShortBuffer +) { val length: Duration get() = tickSize.multipliedBy(signal.capacity().toLong()) 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 22c30e01..1928d73a 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(frame: NumassFrame): Stream + fun analyze(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 new file mode 100644 index 00000000..d2a314b1 --- /dev/null +++ b/numass-core/numass-signal-processing/build.gradle.kts @@ -0,0 +1,19 @@ +plugins { + idea + kotlin("jvm") +} + + +repositories { + mavenLocal() + mavenCentral() +} + +dependencies { + compile(kotlin("stdlib-jdk8")) + compile(project(":numass-core:numass-data-api")) + + // https://mvnrepository.com/artifact/org.apache.commons/commons-collections4 + compile(group = "org.apache.commons", name = "commons-collections4", version = "4.3") + +} \ No newline at end of file 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 new file mode 100644 index 00000000..fedc5cd2 --- /dev/null +++ b/numass-core/numass-signal-processing/src/main/kotlin/inr/numass/data/ChernovProcessor.kt @@ -0,0 +1,66 @@ +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 org.apache.commons.collections4.queue.CircularFifoQueue +import java.nio.ShortBuffer +import java.util.stream.Stream +import kotlin.streams.asStream + + +private fun ShortBuffer.clone(): ShortBuffer { + val clone = ShortBuffer.allocate(capacity()) + rewind()//copy from the beginning + clone.put(this) + rewind() + clone.flip() + return clone +} + + +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() + + private fun CircularFifoQueue.findMax(): Pair { + TODO() + } + + override fun analyze(parent: NumassBlock, frame: NumassFrame): Stream { + return sequence { + val events = HashMap() + val buffer = frame.signal.clone() + + 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 + + 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) + + //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()) + } + } + } + }.asStream() + } +} + diff --git a/numass-core/src/main/kotlin/inr/numass/data/analyzers/SmartAnalyzer.kt b/numass-core/src/main/kotlin/inr/numass/data/analyzers/SmartAnalyzer.kt index 6b9d59b8..f5ed715c 100644 --- a/numass-core/src/main/kotlin/inr/numass/data/analyzers/SmartAnalyzer.kt +++ b/numass-core/src/main/kotlin/inr/numass/data/analyzers/SmartAnalyzer.kt @@ -41,7 +41,7 @@ class SmartAnalyzer(processor: SignalProcessor? = null) : AbstractAnalyzer(proce "simple" -> simpleAnalyzer "time" -> timeAnalyzer "debunch" -> debunchAnalyzer - else -> throw IllegalArgumentException("Analyzer not found") + else -> throw IllegalArgumentException("Analyzer ${config.getString("type")} not found") } } else { if (config.hasValue("t0") || config.hasMeta("t0")) { diff --git a/numass-core/src/main/proto/inr/numas/numass-proto.proto b/numass-core/src/main/proto/inr/numas/numass-proto.proto deleted file mode 100644 index 6ce6f065..00000000 --- a/numass-core/src/main/proto/inr/numas/numass-proto.proto +++ /dev/null @@ -1,33 +0,0 @@ -syntax = "proto3"; - -package inr.numass.data; - -message Point { - // A single channel for multichannel detector readout - message Channel { - //A continuous measurement block - message Block { - // Raw data frame - message Frame { - uint64 time = 1; // Time in nanos from the beginning of the block - bytes data = 2; // Frame data as an array of int16 mesured in arbitrary channels - } - // Event block obtained directly from device of from frame analysis - // In order to save space, times and amplitudes are in separate arrays. - // Amplitude and time with the same index correspond to the same event - message Events { - repeated uint64 times = 1; // Array of time in nanos from the beginning of the block - repeated uint64 amplitudes = 2; // Array of amplitudes of events in channels - } - - uint64 time = 1; // Block start in epoch nanos - repeated Frame frames = 2; // Frames array - Events events = 3; // Events array - uint64 length = 4; // block size in nanos. If missing, take from meta. - uint64 bin_size = 5; // tick size in nanos. Obsolete, to be removed - } - uint64 id = 1; // The number of measuring channel - repeated Block blocks = 2; // Blocks - } - repeated Channel channels = 1; // Array of measuring channels -} \ No newline at end of file diff --git a/numass-main/build.gradle b/numass-main/build.gradle index bdd6a545..eba39a00 100644 --- a/numass-main/build.gradle +++ b/numass-main/build.gradle @@ -5,6 +5,12 @@ plugins { apply plugin: 'kotlin' +//apply plugin: 'org.openjfx.javafxplugin' +// +//javafx{ +// modules = [ 'javafx.controls'] +//} + //if (!hasProperty('mainClass')) { // ext.mainClass = 'inr.numass.LaunchGrindShell' //} diff --git a/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt b/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt index e6097466..28077c06 100644 --- a/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt +++ b/numass-main/src/main/kotlin/inr/numass/NumassPlugin.kt @@ -72,7 +72,8 @@ class NumassPlugin : BasicPlugin() { plotFitTask, histogramTask, fitScanTask, - sliceTask + sliceTask, + subThresholdTask ) @Provides(Task.TASK_TARGET) diff --git a/numass-main/src/main/kotlin/inr/numass/subthreshold/Threshold.kt b/numass-main/src/main/kotlin/inr/numass/subthreshold/Threshold.kt index 077dfe93..5fd5a023 100644 --- a/numass-main/src/main/kotlin/inr/numass/subthreshold/Threshold.kt +++ b/numass-main/src/main/kotlin/inr/numass/subthreshold/Threshold.kt @@ -10,6 +10,7 @@ import hep.dataforge.nullable import hep.dataforge.storage.Storage import hep.dataforge.tables.ListTable import hep.dataforge.tables.Table +import hep.dataforge.toList import hep.dataforge.values.ValueMap import hep.dataforge.values.Values import inr.numass.data.analyzers.* @@ -25,7 +26,9 @@ import org.apache.commons.math3.analysis.ParametricUnivariateFunction import org.apache.commons.math3.exception.DimensionMismatchException import org.apache.commons.math3.fitting.SimpleCurveFitter import org.apache.commons.math3.fitting.WeightedObservedPoint +import org.slf4j.LoggerFactory import java.util.stream.Collectors +import java.util.stream.StreamSupport object Threshold { @@ -35,13 +38,13 @@ object Threshold { //creating storage instance val storage = NumassDirectory.read(context, meta.getString("data.dir")) as Storage - fun Storage.loaders(): Sequence{ + fun Storage.loaders(): Sequence { return sequence { print("Reading ${this@loaders.fullName}") runBlocking { this@loaders.children }.forEach { - if(it is NumassDataLoader){ + if (it is NumassDataLoader) { yield(it) - } else if (it is Storage){ + } else if (it is Storage) { yieldAll(it.loaders()) } } @@ -51,19 +54,19 @@ object Threshold { //Reading points //Free operation. No reading done val sets = storage.loaders() - .filter { it.fullName.toString().matches(meta.getString("data.mask").toRegex()) } + .filter { it.fullName.toString().matches(meta.getString("data.mask").toRegex()) } val analyzer = TimeAnalyzer(); val data = DataSet.edit(NumassPoint::class).also { dataBuilder -> sets.sortedBy { it.startTime } - .flatMap { set -> set.points.asSequence() } - .groupBy { it.voltage } - .forEach { key, value -> - val point = SimpleNumassPoint(value, key) - val name = key.toInt().toString() - dataBuilder.putStatic(name, point, buildMeta("meta", "voltage" to key)); - } + .flatMap { set -> set.points.asSequence() } + .groupBy { it.voltage } + .forEach { key, value -> + val point = SimpleNumassPoint(value, key) + val name = key.toInt().toString() + dataBuilder.putStatic(name, point, buildMeta("meta", "voltage" to key)); + } }.build() return data.pipe(context, meta) { @@ -90,13 +93,14 @@ object Threshold { // ) return binned.rows - .map { - WeightedObservedPoint( - 1.0,//1d / p.getValue() , //weight - it.getDouble(CHANNEL_KEY), // x - it.getDouble(COUNT_RATE_KEY) / binning) //y - } - .collect(Collectors.toList()) + .map { + WeightedObservedPoint( + 1.0,//1d / p.getValue() , //weight + it.getDouble(CHANNEL_KEY), // x + it.getDouble(COUNT_RATE_KEY) / binning + ) //y + } + .collect(Collectors.toList()) } private fun norm(spectrum: Table, xLow: Int, upper: Int): Double { @@ -132,7 +136,6 @@ object Threshold { } - /** * Exponential function $a e^{\frac{x}{\sigma}}$ */ @@ -149,10 +152,10 @@ object Threshold { val norm = norm(spectrum, xLow, upper) return ValueMap.ofPairs( - "U" to voltage, - "a" to a, - "sigma" to sigma, - "correction" to a * sigma * Math.exp(xLow / sigma) / norm + 1.0 + "U" to voltage, + "a" to a, + "sigma" to sigma, + "correction" to a * sigma * Math.exp(xLow / sigma) / norm + 1.0 ) } @@ -173,14 +176,14 @@ object Threshold { val delta = shift ?: parameters[2] return if (parameters.size > 2) { doubleArrayOf( - Math.pow(x - delta, beta), - a * Math.pow(x - delta, beta) * Math.log(x - delta), - -a * beta * Math.pow(x - delta, beta - 1) + Math.pow(x - delta, beta), + a * Math.pow(x - delta, beta) * Math.log(x - delta), + -a * beta * Math.pow(x - delta, beta - 1) ) } else { doubleArrayOf( - Math.pow(x - delta, beta), - a * Math.pow(x - delta, beta) * Math.log(x - delta) + Math.pow(x - delta, beta), + a * Math.pow(x - delta, beta) * Math.log(x - delta) ) } } @@ -206,11 +209,11 @@ object Threshold { val norm = norm(spectrum, xLow, upper) return ValueMap.ofPairs( - "U" to voltage, - "a" to a, - "beta" to beta, - "delta" to delta, - "correction" to a / (beta + 1) * Math.pow(xLow - delta, beta + 1.0) / norm + 1.0 + "U" to voltage, + "a" to a, + "beta" to beta, + "delta" to delta, + "correction" to a / (beta + 1) * Math.pow(xLow - delta, beta + 1.0) / norm + 1.0 ) } @@ -223,23 +226,33 @@ object Threshold { } fun calculateSubThreshold(set: NumassSet, config: Meta, analyzer: NumassAnalyzer = SmartAnalyzer()): Table { - val reference = config.optNumber("reference").nullable?.let{ + val reference = config.optNumber("reference").nullable?.let { set.getPoints(it.toDouble()).firstOrNull() ?: error("Reference point not found") }?.let { println("Using reference point ${it.voltage}") - analyzer.getAmplitudeSpectrum(it,config) + analyzer.getAmplitudeSpectrum(it, config) } return ListTable.Builder().apply { - set.forEach{ point -> - val spectrum = analyzer.getAmplitudeSpectrum(point,config).let { - if(reference == null){ + StreamSupport.stream(set.spliterator(), true).map { point -> + LoggerFactory.getLogger(Threshold.javaClass).info("Starting point ${point.voltage}") + val spectrum = analyzer.getAmplitudeSpectrum(point, config).let { + if (reference == null) { it - } else{ - subtractAmplitudeSpectrum(it,reference) + } else { + subtractAmplitudeSpectrum(it, reference) } } - row(calculateSubThreshold(spectrum,point.voltage,config)) + LoggerFactory.getLogger(Threshold.javaClass).info("Calculating threshold ${point.voltage}") + try { + calculateSubThreshold(spectrum, point.voltage, config) + } catch (ex: Exception) { + LoggerFactory.getLogger(Threshold.javaClass).error("Failed to fit point ${point.voltage}", ex) + null + } + }.toList().filterNotNull().forEach { + println(it.toString()) + row(it) } }.build() } diff --git a/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt b/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt index d2889193..96676ba9 100644 --- a/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt +++ b/numass-main/src/main/kotlin/inr/numass/tasks/NumassTasks.kt @@ -4,8 +4,6 @@ import hep.dataforge.configure import hep.dataforge.data.* import hep.dataforge.io.output.stream import hep.dataforge.io.render -import hep.dataforge.meta.Meta -import hep.dataforge.meta.MetaUtils import hep.dataforge.meta.buildMeta import hep.dataforge.nullable import hep.dataforge.plots.data.DataPlot @@ -86,11 +84,13 @@ val analyzeTask = task("analyze") { info = "Count the number of events for each voltage and produce a table with the results" } model { meta -> - dependsOn(selectTask, meta); - configure(MetaUtils.optEither(meta, "analyzer", "prepare").orElse(Meta.empty())) + dependsOn(selectTask, meta) + configure { + "analyzer" to meta.getMetaOrEmpty("analyzer") + } } pipe { set -> - SmartAnalyzer().analyzeSet(set, meta).also { res -> + SmartAnalyzer().analyzeSet(set, meta.getMeta("analyzer")).also { res -> val outputMeta = meta.builder.putNode("data", set.meta) context.output.render(res, stage = "numass.analyze", name = name, meta = outputMeta) } @@ -299,7 +299,7 @@ val histogramTask = task("histogram") { value( "binning", types = listOf(ValueType.NUMBER), - defaultValue = 20, + defaultValue = 16, info = "The binning of resulting histogram" ) value( @@ -381,7 +381,7 @@ val histogramTask = task("histogram") { data.toSortedMap().forEach { name, set -> putNode("data", buildMeta { "name" to name - set.meta.useMeta("iteration_info"){"iteration" to it} + set.meta.useMeta("iteration_info") { "iteration" to it } }) } } diff --git a/numass-main/src/main/kotlin/inr/numass/tasks/SpecialNumassTasks.kt b/numass-main/src/main/kotlin/inr/numass/tasks/SpecialNumassTasks.kt new file mode 100644 index 00000000..27b0d231 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/tasks/SpecialNumassTasks.kt @@ -0,0 +1,55 @@ +package inr.numass.tasks + +import hep.dataforge.io.render +import hep.dataforge.plots.data.DataPlot +import hep.dataforge.plots.output.plotFrame +import hep.dataforge.plots.plotData +import hep.dataforge.tables.Adapters +import hep.dataforge.tables.Table +import hep.dataforge.tables.filter +import hep.dataforge.useMeta +import hep.dataforge.values.ValueType +import hep.dataforge.workspace.tasks.task +import inr.numass.data.NumassDataUtils +import inr.numass.data.api.NumassSet +import inr.numass.subthreshold.Threshold + +val subThresholdTask = task("threshold") { + descriptor { + value("plot", types = listOf(ValueType.BOOLEAN), defaultValue = false, info = "Show threshold correction plot") + value( + "binning", + types = listOf(ValueType.NUMBER), + defaultValue = 16, + info = "The binning used for fit" + ) + info = "Calculate sub threshold correction" + } + model { meta -> + dependsOn(selectTask, meta) + configure(meta.getMetaOrEmpty("threshold")) + configure { + meta.useMeta("analyzer") { putNode(it) } + setValue("@target", meta.getString("@target", meta.name)) + } + } + join { data -> + val sum = NumassDataUtils.join(name, data.values) + + val correctionTable = Threshold.calculateSubThreshold(sum, meta).filter { + it.getDouble("correction") in (1.0..1.2) + } + + if (meta.getBoolean("plot", false)) { + context.plotFrame("$name.plot", stage = "numass.threshold") { + plots.setType() + plotData("${name}_cor", correctionTable, Adapters.buildXYAdapter("U", "correction")) + plotData("${name}_a", correctionTable, Adapters.buildXYAdapter("U", "a")) + plotData("${name}_beta", correctionTable, Adapters.buildXYAdapter("U", "beta")) + } + } + + context.output.render(correctionTable, "numass.correction", name, meta = meta) + return@join correctionTable + } +} diff --git a/settings.gradle b/settings.gradle index bbd2f6a1..2096fbcd 100644 --- a/settings.gradle +++ b/settings.gradle @@ -15,6 +15,7 @@ include ":numass-core" include 'numass-core:numass-data-api' include 'numass-core:numass-data-proto' +include 'numass-core:numass-signal-processing' //include ":numass-server" //