diff --git a/numass-core/src/main/java/inr/numass/data/NumassDataUtils.java b/numass-core/src/main/java/inr/numass/data/NumassDataUtils.java index 167b8273..c58dfd0c 100644 --- a/numass-core/src/main/java/inr/numass/data/NumassDataUtils.java +++ b/numass-core/src/main/java/inr/numass/data/NumassDataUtils.java @@ -73,8 +73,8 @@ public class NumassDataUtils { // if (first.getVoltage() != second.getVoltage()) { // throw new RuntimeException("Voltage mismatch"); // } -// int[] newArray = new int[first.getSpectrum().length]; -// Arrays.setAll(newArray, i -> first.getSpectrum()[i] + second.getSpectrum()[i]); +// int[] newArray = new int[first.getAmplitudeSpectrum().length]; +// Arrays.setAll(newArray, i -> first.getAmplitudeSpectrum()[i] + second.getAmplitudeSpectrum()[i]); // return new NumassPointImpl( // first.getVoltage(), // Instant.EPOCH, @@ -84,8 +84,8 @@ public class NumassDataUtils { // } // // public static NumassPoint substractPoint(NumassPoint point, NumassPoint reference) { -// int[] array = new int[point.getSpectrum().length]; -// Arrays.setAll(array, i -> Math.max(0, point.getSpectrum()[i] - reference.getSpectrum()[i])); +// int[] array = new int[point.getAmplitudeSpectrum().length]; +// Arrays.setAll(array, i -> Math.max(0, point.getAmplitudeSpectrum()[i] - reference.getAmplitudeSpectrum()[i])); // return new NumassPointImpl( // point.getVoltage(), // point.getTime(), diff --git a/numass-core/src/main/kotlin/inr/numass/data/analyzers/NumassAnalyzer.kt b/numass-core/src/main/kotlin/inr/numass/data/analyzers/NumassAnalyzer.kt index ad2e996c..f4c477b4 100644 --- a/numass-core/src/main/kotlin/inr/numass/data/analyzers/NumassAnalyzer.kt +++ b/numass-core/src/main/kotlin/inr/numass/data/analyzers/NumassAnalyzer.kt @@ -99,9 +99,9 @@ interface NumassAnalyzer { return analyze(block, config).getValue(LENGTH_KEY).numberValue().toLong() } - fun getSpectrum(block: NumassBlock, config: Meta): Table { + fun getAmplitudeSpectrum(block: NumassBlock, config: Meta): Table { val seconds = block.length.toMillis().toDouble() / 1000.0 - return getSpectrum(seconds, getEvents(block, config).asSequence(), config) + return getAmplitudeSpectrum(getEvents(block, config).asSequence(), seconds, config) } companion object { @@ -177,11 +177,12 @@ fun countInWindow(spectrum: Table, loChannel: Short, upChannel: Short): Long { /** * Calculate the amplitude spectrum for a given block. The s * - * @param block + * @param events + * @param length length in seconds, used for count rate calculation * @param config * @return */ -fun getSpectrum(length: Double, events: Sequence, config: Meta = Meta.empty()): Table { +fun getAmplitudeSpectrum(events: Sequence, length: Double, config: Meta = Meta.empty()): Table { val format = TableFormatBuilder() .addNumber(NumassAnalyzer.CHANNEL_KEY, X_VALUE_KEY) .addNumber(NumassAnalyzer.COUNT_KEY) diff --git a/numass-core/src/main/kotlin/inr/numass/data/analyzers/TimeAnalyzer.kt b/numass-core/src/main/kotlin/inr/numass/data/analyzers/TimeAnalyzer.kt index 6c3e3c03..2952459a 100644 --- a/numass-core/src/main/kotlin/inr/numass/data/analyzers/TimeAnalyzer.kt +++ b/numass-core/src/main/kotlin/inr/numass/data/analyzers/TimeAnalyzer.kt @@ -157,12 +157,12 @@ class TimeAnalyzer @JvmOverloads constructor(private val processor: SignalProces return block.events.count().toDouble() / block.length.toMillis() * 1000 } - fun getEventsPairs(block: NumassBlock, config: Meta): Sequence> { - return Sequence { getEvents(block, config).iterator() }.zipWithNext() + fun zipEvents(block: NumassBlock, config: Meta): Sequence> { + return getAllEvents(block).asSequence().zipWithNext() } /** - * The chain of event times in nanos + * The chain of event with delays in nanos * * @param block * @param config @@ -173,21 +173,6 @@ class TimeAnalyzer @JvmOverloads constructor(private val processor: SignalProces val delay = Math.max(next.timeOffset - prev.timeOffset, 0) Pair(prev, delay) }.asStream() -// val lastEvent = AtomicReference(null) -// -// val eventStream = super.getEvents(block, config)//using super implementation -// -// return eventStream.map { event -> -// var res = if (lastEvent.get() == null) 0L else event.timeOffset - lastEvent.get().timeOffset -// -// if (res < 0) { -// res = 0L -// } -// -// lastEvent.set(event) -// //TODO remove autoboxing somehow -// Pair(event, res) -// } } /** diff --git a/numass-main/src/main/groovy/inr/numass/scripts/temp/JoinSpectra.groovy b/numass-main/src/main/groovy/inr/numass/scripts/temp/JoinSpectra.groovy index c04a9946..26a0fc64 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/temp/JoinSpectra.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/temp/JoinSpectra.groovy @@ -67,8 +67,8 @@ new GrindShell(ctx).eval { frame.plots.configure(showErrors: false, showSymbol: false, showLine: true, connection: "step") joined.points.filter { it.voltage in [14000d, 15000d, 16000d, 17000d, 18000d] }.forEach { - //Table spectrum = analyzer.getSpectrum(it, Meta.empty()).withBinning(20).withDeadTime() - Table spectrum = analyzer.getSpectrum(it, Grind.buildMeta(t0: t0*1000)).withBinning(20).withDeadTime(t0) + //Table spectrum = analyzer.getAmplitudeSpectrum(it, Meta.empty()).withBinning(20).withDeadTime() + Table spectrum = analyzer.getAmplitudeSpectrum(it, Grind.buildMeta(t0: t0*1000)).withBinning(20).withDeadTime(t0) frame.add(DataPlot.plot(it.voltage.toString(), adapter, spectrum)) } @@ -79,7 +79,7 @@ new GrindShell(ctx).eval { // pointFrame.plots.configure(showErrors: false, showSymbol: false, showLine: true, connection: "step") // // [0, 5, 10,15,20].forEach{ -// Table spectrum = analyzer.getSpectrum(point, Grind.buildMeta(t0: it*1000)).withBinning(20).withDeadTime(it) +// Table spectrum = analyzer.getAmplitudeSpectrum(point, Grind.buildMeta(t0: it*1000)).withBinning(20).withDeadTime(it) // pointFrame.add(DataPlot.plot(it.toString(), adapter, spectrum)) // } diff --git a/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowUtils.groovy b/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowUtils.groovy index 1fa17eac..7c5aa5a6 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowUtils.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/underflow/UnderflowUtils.groovy @@ -58,7 +58,7 @@ class UnderflowUtils { def generate = GrindPipe.build("generate") { result { input -> - return analyzer.getSpectrum(input as NumassPoint, delegate.meta) + return analyzer.getAmplitudeSpectrum(input as NumassPoint, delegate.meta) } } diff --git a/numass-main/src/main/java/inr/numass/utils/UnderflowCorrection.java b/numass-main/src/main/java/inr/numass/utils/UnderflowCorrection.java index b145b0fb..dabcd384 100644 --- a/numass-main/src/main/java/inr/numass/utils/UnderflowCorrection.java +++ b/numass-main/src/main/java/inr/numass/utils/UnderflowCorrection.java @@ -74,7 +74,7 @@ public class UnderflowCorrection { // } public Values fitPoint(NumassPoint point, int xLow, int xHigh, int upper, int binning) { - Table spectrum = analyzer.getSpectrum(point, Meta.empty()); + Table spectrum = analyzer.getAmplitudeSpectrum(point, Meta.empty()); double norm = spectrum.getRows().filter(row -> { int channel = row.getInt(CHANNEL_KEY); diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/Correlation.kt b/numass-main/src/main/kotlin/inr/numass/scripts/Correlation.kt new file mode 100644 index 00000000..e4a81277 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/scripts/Correlation.kt @@ -0,0 +1,73 @@ +/* + * Copyright 2017 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.scripts + +import hep.dataforge.kodex.buildContext +import hep.dataforge.kodex.buildMeta +import inr.numass.NumassPlugin +import inr.numass.data.NumassDataUtils +import inr.numass.data.analyzers.SmartAnalyzer +import inr.numass.data.api.NumassEvent +import inr.numass.data.api.NumassSet +import inr.numass.data.storage.NumassStorageFactory +import org.apache.commons.math3.stat.correlation.PearsonsCorrelation +import java.util.stream.Stream + + +private fun correlation(sequence: Stream): Double { + val amplitudes: MutableList = ArrayList() + val times: MutableList = ArrayList() + sequence.forEach { + amplitudes.add(it.chanel.toDouble()) + times.add(it.timeOffset.toDouble()) + } + + return PearsonsCorrelation().correlation(amplitudes.toDoubleArray(), times.toDoubleArray()) +} + +fun main(args: Array) { + + val context = buildContext("NUMASS", NumassPlugin::class.java) { + rootDir = "D:\\Work\\Numass\\sterile\\2017_05" + dataDir = "D:\\Work\\Numass\\data\\2017_05" + } + //val rootDir = File("D:\\Work\\Numass\\data\\2017_05\\Fill_2") + + val storage = NumassStorageFactory.buildLocal(context, "Fill_2", true, false); + + val sets = (2..14).map { "set_$it" } + + val loaders = sets.mapNotNull { set -> + storage.provide("loader::$set", NumassSet::class.java).orElse(null) + } + + val set = NumassDataUtils.join("sum", loaders) + + val analyzer = SmartAnalyzer(); + + val meta = buildMeta { + "window.lo" to 400 + "window.up" to 2500 + } + + println("Correlation between amplitudes and delays:") + set.points.filter { it.voltage < 16000.0 }.forEach { + val cor = correlation(analyzer.getEvents(it, meta)) + println("${it.voltage}: $cor") + } + +} \ No newline at end of file diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt b/numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt index e3027f72..483e8b8c 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/DifferentialSpectrum.kt @@ -21,17 +21,21 @@ import hep.dataforge.kodex.buildContext import hep.dataforge.kodex.buildMeta import hep.dataforge.kodex.replaceColumn import hep.dataforge.plots.data.DataPlot +import hep.dataforge.tables.Table import inr.numass.NumassPlugin import inr.numass.data.NumassDataUtils import inr.numass.data.analyzers.NumassAnalyzer -import inr.numass.data.analyzers.SmartAnalyzer +import inr.numass.data.analyzers.TimeAnalyzer +import inr.numass.data.analyzers.getAmplitudeSpectrum +import inr.numass.data.api.NumassPoint import inr.numass.data.api.NumassSet import inr.numass.data.storage.NumassStorageFactory + fun main(args: Array) { val context = buildContext("NUMASS", NumassPlugin::class.java, PlotManager::class.java) { - rootDir = "D:\\Work\\Numass\\sterile2017_05" + rootDir = "D:\\Work\\Numass\\sterile\\2017_05" dataDir = "D:\\Work\\Numass\\data\\2017_05" } //val rootDir = File("D:\\Work\\Numass\\data\\2017_05\\Fill_2") @@ -44,14 +48,21 @@ fun main(args: Array) { storage.provide("loader::$set", NumassSet::class.java).orElse(null) } - val analyzer = SmartAnalyzer() - val all = NumassDataUtils.join("sum", loaders) + val t0 = 20e3 val meta = buildMeta { "window.lo" to 400 "window.up" to 1800 + "t0" to t0 + } + + val filter: (NumassPoint) -> Table = { point -> + val analyzer = TimeAnalyzer() + val sequence = analyzer.zipEvents(point, meta).filter { it.second.timeOffset - it.first.timeOffset > t0 }.map { it.second } + //val sequence = analyzer.getEvents(point,meta).asSequence() + getAmplitudeSpectrum(sequence, point.length.toMillis().toDouble() / 1000.0, meta) } val plots = context.getFeature(PlotManager::class.java) @@ -60,22 +71,22 @@ fun main(args: Array) { val integralFrame = plots.getPlotFrame("integral") - for (hv in arrayOf(14000.0, 14200.0, 14400.0, 14600.0, 14800.0, 15000.0)) { - val point1 = all.optPoint(hv).get() + for (hv in arrayOf(14000.0, 14500.0, 15000.0, 15500.0)) { + val basePoint = all.optPoint(hv).get() - val point0 = all.optPoint(hv + 200.0).get() + val subPoint = all.optPoint(hv + 200.0).get() with(NumassAnalyzer) { - val spectrum1 = analyzer.getSpectrum(point1, meta).withBinning(20) + val baseSpectrum = filter(basePoint).withBinning(50) - val spectrum0 = analyzer.getSpectrum(point0, meta).withBinning(20) + val subSpectrum = filter(subPoint).withBinning(50) - val res = subtractAmplitudeSpectrum(spectrum1, spectrum0) + val res = subtractAmplitudeSpectrum(baseSpectrum, subSpectrum) val norm = res.getColumn(COUNT_RATE_KEY).stream().mapToDouble { it.doubleValue() }.sum() - integralFrame.add(DataPlot.plot("point_$hv", AMPLITUDE_ADAPTER, spectrum0)) + integralFrame.add(DataPlot.plot("point_$hv", AMPLITUDE_ADAPTER, baseSpectrum)) frame.add(DataPlot.plot("point_$hv", AMPLITUDE_ADAPTER, res.replaceColumn(COUNT_RATE_KEY) { getDouble(COUNT_RATE_KEY) / norm })) } diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt b/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt index 9cf2bd0f..9a032392 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/InversedChain.kt @@ -27,7 +27,7 @@ import inr.numass.NumassPlugin import inr.numass.data.NumassDataUtils import inr.numass.data.analyzers.NumassAnalyzer import inr.numass.data.analyzers.TimeAnalyzer -import inr.numass.data.analyzers.getSpectrum +import inr.numass.data.analyzers.getAmplitudeSpectrum import inr.numass.data.api.NumassSet import inr.numass.data.storage.NumassStorageFactory import kotlin.streams.asSequence @@ -71,16 +71,16 @@ fun main(args: Array) { } with(NumassAnalyzer) { - val events = getSpectrum(seconds, analyzer.getEvents(point).asSequence(),meta) + val events = getAmplitudeSpectrum(analyzer.getEvents(point).asSequence(), seconds, meta) .withBinning(binning) val eventsNorming = events.getColumn(COUNT_RATE_KEY).stream().mapToDouble{it.doubleValue()}.sum() println("The norming factor for unfiltered count rate is $eventsNorming") - val filtered = getSpectrum( + val filtered = getAmplitudeSpectrum( + analyzer.zipEvents(point, Meta.empty()).filter { it.second.timeOffset - it.first.timeOffset > t0 }.map { it.second }, seconds, - analyzer.getEventsPairs(point, Meta.empty()).filter { it.second.timeOffset - it.first.timeOffset > t0 }.map { it.second }, meta ).withBinning(binning) @@ -88,9 +88,9 @@ fun main(args: Array) { println("The norming factor for filtered count rate is $filteredNorming") - val defaultFiltered = getSpectrum( - seconds, + val defaultFiltered = getAmplitudeSpectrum( analyzer.getEvents(point, buildMeta {"t0" to t0}).asSequence(), + seconds, meta ).withBinning(binning) diff --git a/numass-viewer/src/main/kotlin/inr/numass/viewer/AmplitudeView.kt b/numass-viewer/src/main/kotlin/inr/numass/viewer/AmplitudeView.kt index dc77a3bb..1215ac72 100644 --- a/numass-viewer/src/main/kotlin/inr/numass/viewer/AmplitudeView.kt +++ b/numass-viewer/src/main/kotlin/inr/numass/viewer/AmplitudeView.kt @@ -105,7 +105,7 @@ class AmplitudeView( * Calculate or get spectrum from the cache */ private suspend fun getSpectrum(point: NumassPoint): Table { - return cache.computeIfAbsent(point) { analyzer.getSpectrum(point, Meta.empty()) } + return cache.computeIfAbsent(point) { analyzer.getAmplitudeSpectrum(point, Meta.empty()) } } /** diff --git a/numass-viewer/src/main/kotlin/inr/numass/viewer/SpectrumView.kt b/numass-viewer/src/main/kotlin/inr/numass/viewer/SpectrumView.kt index 0f22ce55..1d9e0bd9 100644 --- a/numass-viewer/src/main/kotlin/inr/numass/viewer/SpectrumView.kt +++ b/numass-viewer/src/main/kotlin/inr/numass/viewer/SpectrumView.kt @@ -154,7 +154,7 @@ class SpectrumView( } private fun getSpectrum(point: NumassPoint): Table { - return cache.computeIfAbsent(point) { analyzer.getSpectrum(point, Meta.empty()) } + return cache.computeIfAbsent(point) { analyzer.getAmplitudeSpectrum(point, Meta.empty()) } }