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 aa3d37c0..4f645621 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 @@ -26,13 +26,19 @@ import hep.dataforge.values.Value import hep.dataforge.values.ValueMap import hep.dataforge.values.ValueType import hep.dataforge.values.Values +import inr.numass.data.analyzers.NumassAnalyzer.Companion.COUNT_KEY +import inr.numass.data.analyzers.NumassAnalyzer.Companion.COUNT_RATE_ERROR_KEY +import inr.numass.data.analyzers.NumassAnalyzer.Companion.COUNT_RATE_KEY +import inr.numass.data.analyzers.NumassAnalyzer.Companion.LENGTH_KEY import inr.numass.data.api.* import inr.numass.data.api.NumassPoint.Companion.HV_KEY import java.util.* import java.util.concurrent.atomic.AtomicLong import java.util.stream.Stream +import kotlin.math.sqrt import kotlin.streams.asSequence import kotlin.streams.asStream +import kotlin.streams.toList /** @@ -56,12 +62,7 @@ class TimeAnalyzer @JvmOverloads constructor(private val processor: SignalProces val res = getEventsWithDelay(block, config).asSequence().chunked(chunkSize) { analyzeSequence(it.asSequence(), t0) - }.reduce(this::combineBlockResults) ?: ValueMap.ofPairs( - NumassAnalyzer.LENGTH_KEY to 0, - NumassAnalyzer.COUNT_KEY to 0, - NumassAnalyzer.COUNT_RATE_KEY to 0, - NumassAnalyzer.COUNT_RATE_ERROR_KEY to 0 - ) + }.toList().average() return ValueMap.Builder(res) .putValue(NumassAnalyzer.WINDOW_KEY, arrayOf(loChannel, upChannel)) @@ -92,15 +93,7 @@ class TimeAnalyzer @JvmOverloads constructor(private val processor: SignalProces NumassAnalyzer.COUNT_RATE_KEY to countRate, NumassAnalyzer.COUNT_RATE_ERROR_KEY to countRateError ) -// ValueMap.of(NAME_LIST, -// length, -// count, -// countRate, -// countRateError, -// arrayOf(loChannel, upChannel), -// block.startTime, -// t0.toDouble() / 1000.0 -// ) + } override fun analyzeParent(point: ParentBlock, config: Meta): Values { @@ -108,7 +101,7 @@ class TimeAnalyzer @JvmOverloads constructor(private val processor: SignalProces val res = point.blocks.stream() .filter { it.events.findAny().isPresent }// filter for empty blocks .map { it -> analyze(it, config) } - .reduce(null) { v1, v2 -> this.combineBlockResults(v1, v2) } + .toList().average() val map = HashMap(res.asMap()) if (point is NumassPoint) { @@ -118,38 +111,21 @@ class TimeAnalyzer @JvmOverloads constructor(private val processor: SignalProces } /** - * Combine two blocks from the same point into one + * Combine multiple blocks from the same point into one * - * @param v1 - * @param v2 * @return */ - private fun combineBlockResults(v1: Values?, v2: Values?): Values? { - if (v1 == null) { - return v2 - } - if (v2 == null) { - return v1 - } + private fun List.average(): Values { + val totalTime = sumByDouble { it.getDouble(LENGTH_KEY) } + val countRate = sumByDouble { it.getDouble(COUNT_RATE_KEY) * it.getDouble(LENGTH_KEY) } / totalTime + val countRateDispersion = sumByDouble { Math.pow(it.getDouble(COUNT_RATE_ERROR_KEY) * it.getDouble(LENGTH_KEY) / totalTime, 2.0) } - val t1 = v1.getDouble(NumassAnalyzer.LENGTH_KEY) - val t2 = v2.getDouble(NumassAnalyzer.LENGTH_KEY) - val cr1 = v1.getDouble(NumassAnalyzer.COUNT_RATE_KEY) - val cr2 = v2.getDouble(NumassAnalyzer.COUNT_RATE_KEY) - val err1 = v1.getDouble(NumassAnalyzer.COUNT_RATE_ERROR_KEY) - val err2 = v2.getDouble(NumassAnalyzer.COUNT_RATE_ERROR_KEY) - - val countRate = (t1 * cr1 + t2 * cr2) / (t1 + t2) - - val countRateErr = Math.sqrt(Math.pow(t1 * err1 / (t1 + t2), 2.0) + Math.pow(t2 * err2 / (t1 + t2), 2.0)) - - - return ValueMap.Builder(v1) - .putValue(NumassAnalyzer.LENGTH_KEY, t1 + t2) - .putValue(NumassAnalyzer.COUNT_KEY, v1.getInt(NumassAnalyzer.COUNT_KEY) + v2.getInt(NumassAnalyzer.COUNT_KEY)) - .putValue(NumassAnalyzer.COUNT_RATE_KEY, countRate) - .putValue(NumassAnalyzer.COUNT_RATE_ERROR_KEY, countRateErr) + return ValueMap.Builder(first()) + .putValue(LENGTH_KEY, totalTime) + .putValue(COUNT_KEY, sumBy { it.getInt(COUNT_KEY) }) + .putValue(COUNT_RATE_KEY, countRate) + .putValue(COUNT_RATE_ERROR_KEY, sqrt(countRateDispersion)) .build() } diff --git a/numass-core/src/main/kotlin/inr/numass/data/api/NumassBlock.kt b/numass-core/src/main/kotlin/inr/numass/data/api/NumassBlock.kt index 0cb5ea79..601ecfd5 100644 --- a/numass-core/src/main/kotlin/inr/numass/data/api/NumassBlock.kt +++ b/numass-core/src/main/kotlin/inr/numass/data/api/NumassBlock.kt @@ -82,7 +82,7 @@ class SimpleBlock( private val eventList = runBlocking { producer(this@SimpleBlock).toList()} - override val frames: Stream = Stream.empty() + override val frames: Stream get() = Stream.empty() override val events: Stream get() = eventList.stream() diff --git a/numass-main/src/main/groovy/inr/numass/scripts/times/AnalyzePoint.groovy b/numass-main/src/main/groovy/inr/numass/scripts/times/AnalyzePoint.groovy deleted file mode 100644 index 1f95f56c..00000000 --- a/numass-main/src/main/groovy/inr/numass/scripts/times/AnalyzePoint.groovy +++ /dev/null @@ -1,76 +0,0 @@ -package inr.numass.scripts.times - -import hep.dataforge.context.Context -import hep.dataforge.context.Global -import hep.dataforge.data.DataSet -import hep.dataforge.fx.plots.FXPlotManager -import hep.dataforge.grind.Grind -import hep.dataforge.grind.GrindShell -import hep.dataforge.meta.Meta -import inr.numass.NumassPlugin -import inr.numass.actions.TimeAnalyzerAction -import inr.numass.data.NumassDataUtils -import inr.numass.data.api.NumassPoint -import inr.numass.data.api.NumassSet -import inr.numass.data.api.SimpleNumassPoint -import inr.numass.data.storage.NumassStorageFactory - -/** - * Created by darksnake on 27-Jun-17. - */ - - -Context ctx = Global.instance() -ctx.getPluginManager().load(FXPlotManager) -ctx.getPluginManager().load(NumassPlugin) - -new GrindShell(ctx).eval { - - def storage = NumassStorageFactory.buildLocal(ctx, "D:\\Work\\Numass\\data\\2017_05\\Fill_2", true, false); - - Meta meta = Grind.buildMeta(binNum: 200, inverted: true) { - window(lo: 500, up: 3000) - plot(showErrors: false) - } - - //def sets = ((2..14) + (22..31)).collect { "set_$it" } - def sets = (2..14).collect { "set_$it" } - //def sets = (16..31).collect { "set_$it" } - //def sets = (20..28).collect { "set_$it" } - - def loaders = sets.collect { set -> - storage.provide("loader::$set", NumassSet.class).orElse(null) - }.findAll { it != null } - - def hvs = [14000d]//[14000d, 14200d, 14600d, 14800d, 15000d, 15200d, 15400d, 15600d, 15800d] - - def all = NumassDataUtils.join("sum", loaders) - - def builder = DataSet.builder(NumassPoint) - - hvs.each { hv -> - builder.putStatic("point_${hv as int}", new SimpleNumassPoint(hv, all.points.filter { - it.voltage == hv - }.collect())); - } - - def data = builder.build() - -// def hv = 14000; -// def dataBuilder = DataSet.builder(NumassPoint) -// -// StorageUtils.loaderStream(storage, false) -// .filter { it.value instanceof NumassSet } -// .forEach { pair -> -// (pair.value as NumassSet).optPoint(hv).ifPresent { -// dataBuilder.putData(pair.key, it, it.meta); -// } -// } -// def data = dataBuilder.build() - - def result = new TimeAnalyzerAction().run(ctx, data, meta); - - result.computeAll(); - -// storage.close() -} \ No newline at end of file diff --git a/numass-main/src/main/groovy/inr/numass/scripts/times/TestAnalyzer.groovy b/numass-main/src/main/groovy/inr/numass/scripts/times/TestAnalyzer.groovy deleted file mode 100644 index ca60ae0e..00000000 --- a/numass-main/src/main/groovy/inr/numass/scripts/times/TestAnalyzer.groovy +++ /dev/null @@ -1,44 +0,0 @@ -package inr.numass.scripts.times - -import hep.dataforge.context.Context -import hep.dataforge.context.Global -import hep.dataforge.fx.plots.FXPlotManager -import hep.dataforge.grind.GrindShell -import hep.dataforge.meta.Meta -import inr.numass.NumassPlugin -import inr.numass.actions.TimeAnalyzerAction -import inr.numass.data.GeneratorKt -import inr.numass.data.api.SimpleNumassPoint -import org.apache.commons.math3.random.JDKRandomGenerator - -import java.time.Instant - -/** - * Created by darksnake on 27-Jun-17. - */ - - -Context ctx = Global.INSTANCE -ctx.getPluginManager().load(FXPlotManager) -ctx.getPluginManager().load(NumassPlugin.class) - -new GrindShell(ctx).eval { - - double cr = 30e3; - long length = 30e9; - def num = 5; - def dt = 6.5 - - - - def blocks = (1..num).collect { - def chain = GeneratorKt.generateEvents(cr, new JDKRandomGenerator(),{10000}) - GeneratorKt.generateBlock(Instant.now().plusNanos(it * length), length, chain) - } - - def point = new SimpleNumassPoint(10000, blocks) - - def meta = Meta.empty()//Grind.buildMeta(plotHist: false) - - new TimeAnalyzerAction().simpleRun(point, meta); -} \ No newline at end of file diff --git a/numass-main/src/main/kotlin/inr/numass/actions/TimeAnalyzerAction.kt b/numass-main/src/main/kotlin/inr/numass/actions/TimeAnalyzerAction.kt index 2b2b37ad..56205432 100644 --- a/numass-main/src/main/kotlin/inr/numass/actions/TimeAnalyzerAction.kt +++ b/numass-main/src/main/kotlin/inr/numass/actions/TimeAnalyzerAction.kt @@ -13,6 +13,7 @@ import hep.dataforge.tables.Table import hep.dataforge.values.ValueType import inr.numass.data.analyzers.NumassAnalyzer import inr.numass.data.analyzers.TimeAnalyzer +import inr.numass.data.analyzers.TimeAnalyzer.Companion.T0_KEY import inr.numass.data.api.NumassPoint /** @@ -39,9 +40,10 @@ class TimeAnalyzerAction : OneToOneAction() { val pm = context[PlotPlugin::class.java]; - val trueCR = analyzer.analyze(input, inputMeta).getDouble("cr") + val initialEstimate = analyzer.analyze(input, inputMeta) + val trueCR = initialEstimate.getDouble("cr") - log.report("The expected count rate for 30 us delay is $trueCR") + log.report("The expected count rate for ${initialEstimate.getDouble(T0_KEY)} us delay is $trueCR") val binNum = inputMeta.getInt("binNum", 1000); val binSize = inputMeta.getDouble("binSize", 1.0 / trueCR * 10 / binNum * 1e6) @@ -92,7 +94,16 @@ class TimeAnalyzerAction : OneToOneAction() { configure(inputMeta.getMetaOrEmpty("plot")) } - pm.getPlotFrame(name, "stat-method").add(statPlot) + pm.getPlotFrame(name, "stat-method") + .configure { + "xAxis" to { + "title" to "delay" + "units" to "us" + } + "yAxis" to { + "title" to "Relative count rate" + } + }.add(statPlot) (1..100).map { inputMeta.getDouble("t0Step", 1000.0) * it }.map { t -> val result = analyzer.analyze(input, inputMeta.builder.setValue("t0", t)) @@ -103,7 +114,7 @@ class TimeAnalyzerAction : OneToOneAction() { } else { 1.0 } - if(Thread.currentThread().isInterrupted){ + if (Thread.currentThread().isInterrupted) { throw InterruptedException() } statPlot.append( diff --git a/numass-main/src/main/kotlin/inr/numass/data/Generator.kt b/numass-main/src/main/kotlin/inr/numass/data/Generator.kt index c6b63c24..f5a53edd 100644 --- a/numass-main/src/main/kotlin/inr/numass/data/Generator.kt +++ b/numass-main/src/main/kotlin/inr/numass/data/Generator.kt @@ -87,7 +87,7 @@ fun buildBunchChain( BunchState(0, 0), OrphanNumassEvent(rnd.amp(null, 0), 0)) { event -> if (event.timeOffset >= bunchEnd) { - bunchStart = bunchEnd + (rnd.nextDeltaTime(bunchRate)).toLong() + bunchStart = bunchEnd + rnd.nextDeltaTime(bunchRate) bunchEnd = bunchStart + (bunchLength * 1e9).toLong() OrphanNumassEvent(rnd.amp(null, 0), bunchStart) } else { @@ -111,106 +111,4 @@ fun mergeEventChains(vararg chains: Chain): Chain { -// -// protected abstract fun next(event: NumassEvent?, state: S = buildState()): NumassEvent -// -// fun buildSequence(): Sequence { -// val state = buildState() -// return generateSequence(seed = null) { event: NumassEvent? -> -// next(event, state) -// } -// } -// -// protected abstract fun buildState(): S -// -// fun generateBlock(start: Instant, length: Long): NumassBlock { -// val events = buildSequence().takeWhile { it.timeOffset < length }.toList() -// return SimpleBlock(start, Duration.ofNanos(length), events) -// } -//} -// -// -//class SimpleChainGenerator( -// val cr: Double, -// private val rnd: RandomGenerator = JDKRandomGenerator(), -// private val amp: RandomGenerator.(NumassEvent?, Long) -> Short = { _, _ -> ((nextDouble() + 2.0) * 100).toShort() } -//) : ChainGenerator() { -// -// override fun buildState() { -// return Unit -// } -// -// override fun next(event: NumassEvent?, state: Unit): NumassEvent { -// return if (event == null) { -// NumassEvent(rnd.amp(null, 0), Instant.EPOCH, 0) -// } else { -// val deltaT = rnd.nextDeltaTime(cr) -// NumassEvent(rnd.amp(event, deltaT), event.blockTime, event.timeOffset + deltaT) -// } -// } -// -// fun next(event: NumassEvent?): NumassEvent { -// return next(event, Unit) -// } -//} -// -//class BunchGenerator( -// private val cr: Double, -// private val bunchRate: Double, -// private val bunchLength: RandomGenerator.() -> Long, -// private val rnd: RandomGenerator = JDKRandomGenerator(), -// private val amp: RandomGenerator.(NumassEvent?, Long) -> Short = { _, _ -> ((nextDouble() + 2.0) * 100).toShort() } -//) : ChainGenerator() { -// -// private val internalGenerator = SimpleChainGenerator(cr, rnd, amp) -// -// class BunchState(var bunchStart: Long = 0, var bunchEnd: Long = 0) -// -// override fun next(event: NumassEvent?, state: BunchState): NumassEvent { -// if (event?.timeOffset ?: 0 >= state.bunchEnd) { -// state.bunchStart = state.bunchEnd + (rnd.nextExp(bunchRate) * 1e9).toLong() -// state.bunchEnd = state.bunchStart + rnd.bunchLength() -// return NumassEvent(rnd.amp(null, 0), Instant.EPOCH, state.bunchStart) -// } else { -// return internalGenerator.next(event) -// } -// } -// -// override fun buildState(): BunchState { -// return BunchState(0, 0) -// } -//} -// -// -//class MergingGenerator(private vararg val generators: ChainGenerator<*>) : ChainGenerator() { -// -// inner class MergingState { -// val queue: PriorityQueue, NumassEvent>> = -// PriorityQueue(Comparator.comparing, NumassEvent>, Long> { it.second.timeOffset }) -// -// init { -// generators.forEach { generator -> -// val sequence = generator.buildSequence() -// queue.add(Pair(sequence, sequence.iterator().next())) -// } -// } -// } -// -// override fun next(event: NumassEvent?, state: MergingState): NumassEvent { -// val pair = state.queue.poll() -// state.queue.add(Pair(pair.first, pair.first.iterator().next())) -// return pair.second -// } -// -// override fun buildState(): MergingState { -// return MergingState() -// } -//} -// +} \ No newline at end of file diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/AnalyzePoint.kt b/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/AnalyzePoint.kt index 14ff5083..dcb1eace 100644 --- a/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/AnalyzePoint.kt +++ b/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/AnalyzePoint.kt @@ -19,12 +19,13 @@ fun main(args: Array) { dataDir = "D:\\Work\\Numass\\data\\2018_04" } - val storage = NumassStorageFactory.buildLocal(context, "Fill_2", true, false); + val storage = NumassStorageFactory.buildLocal(context, "Fill_3", true, false); val meta = buildMeta { + "t0" to 3000 "binNum" to 200 - "t0Step" to 1000 - "chunkSize" to 1000 + "t0Step" to 100 + "chunkSize" to 3000 node("window") { "lo" to 0 "up" to 4000 @@ -33,7 +34,7 @@ fun main(args: Array) { } //def sets = ((2..14) + (22..31)).collect { "set_$it" } - val sets = (5..6).map { "set_$it" } + val sets = (2..12).map { "set_$it" } //def sets = (16..31).collect { "set_$it" } //def sets = (20..28).collect { "set_$it" } @@ -41,7 +42,7 @@ fun main(args: Array) { storage.provide("loader::$set", NumassSet::class.java).orElse(null) }.filter { it != null } - val hvs = listOf(15000.0)//, 15000d, 15200d, 15400d, 15600d, 15800d] + val hvs = listOf(12000.0, 14000.0, 16000.0)//, 15000d, 15200d, 15400d, 15600d, 15800d] val all = NumassDataUtils.join("sum", loaders) diff --git a/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/TestAnalyzer.kt b/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/TestAnalyzer.kt new file mode 100644 index 00000000..d7341459 --- /dev/null +++ b/numass-main/src/main/kotlin/inr/numass/scripts/timeanalysis/TestAnalyzer.kt @@ -0,0 +1,56 @@ +package inr.numass.scripts.timeanalysis + +import hep.dataforge.fx.plots.FXPlotManager +import hep.dataforge.kodex.buildMeta +import hep.dataforge.maths.chain.MarkovChain +import inr.numass.NumassPlugin +import inr.numass.actions.TimeAnalyzerAction +import inr.numass.data.api.OrphanNumassEvent +import inr.numass.data.api.SimpleNumassPoint +import inr.numass.data.api.timeOffset +import inr.numass.data.generateBlock +import org.apache.commons.math3.random.JDKRandomGenerator +import org.apache.commons.math3.random.RandomGenerator +import java.lang.Math.exp +import java.time.Instant + +fun main(args: Array) { + FXPlotManager().startGlobal() + NumassPlugin().startGlobal() + + val cr = 10e3; + val length = 30e9.toLong(); + val num = 20; + val dt = 6.5 + + val rnd = JDKRandomGenerator() + + fun RandomGenerator.nextExp(mean: Double): Double { + return -mean * Math.log(1 - nextDouble()) + } + + fun RandomGenerator.nextDeltaTime(cr: Double): Long { + return (nextExp(1.0 / cr) * 1e9).toLong() + } + + + val chain = MarkovChain(OrphanNumassEvent(1000, 0)) { event -> + val deltaT = rnd.nextDeltaTime(cr * exp(- event.timeOffset / 1e11)) + OrphanNumassEvent(1000, event.timeOffset + deltaT) + } + + val blocks = (1..num).map { + generateBlock(Instant.now().plusNanos(it * length), length, chain) + } + + val point = SimpleNumassPoint(blocks, 12000.0) + + val meta = buildMeta { + "t0" to 3000 + "binNum" to 200 + "t0Step" to 200 + "chunkSize" to 5000 + } + + TimeAnalyzerAction().simpleRun(point, meta); +} \ No newline at end of file