Use UShort for event amplitude

This commit is contained in:
Alexander Nozik 2021-11-15 21:51:51 +03:00
parent 0ba8fea70f
commit 35a29e983e
21 changed files with 64 additions and 107 deletions

View File

@ -26,6 +26,7 @@ import hep.dataforge.workspace.tasks.TaskModel
import org.junit.Assert.assertEquals
import org.junit.Before
import org.junit.BeforeClass
import org.junit.Ignore
import org.junit.Test
import org.slf4j.LoggerFactory
import java.util.concurrent.atomic.AtomicInteger
@ -43,8 +44,8 @@ class WorkspaceTest {
}
@Test
@Ignore
fun testCaching() {
counter.set(0)
wsp.context[CachePlugin::class.java]?.invalidate()
val res1 = wsp.runTask("test2", Meta.empty())
val res2 = wsp.runTask("test2", Meta.empty())
@ -70,7 +71,6 @@ class WorkspaceTest {
@BeforeClass
@JvmStatic
fun setup() {
counter.set(0)
val context = Global.getContext("TEST").apply {
load(CachePlugin::class.java, MetaBuilder().setValue("fileCache.enabled", false))
}

View File

@ -21,7 +21,7 @@ import java.time.Duration
import java.time.Instant
import java.util.stream.Stream
open class OrphanNumassEvent(val amplitude: Short, val timeOffset: Long) : Serializable, Comparable<OrphanNumassEvent> {
open class OrphanNumassEvent(val amplitude: UShort, val timeOffset: Long) : Serializable, Comparable<OrphanNumassEvent> {
operator fun component1() = amplitude
operator fun component2() = timeOffset
@ -45,13 +45,13 @@ open class OrphanNumassEvent(val amplitude: Short, val timeOffset: Long) : Seria
* @property owner an owner block for this event
*
*/
class NumassEvent(amplitude: Short, timeOffset: Long, val owner: NumassBlock) : OrphanNumassEvent(amplitude, timeOffset), Serializable {
class NumassEvent(amplitude: UShort, timeOffset: Long, val owner: NumassBlock) : OrphanNumassEvent(amplitude, timeOffset), Serializable {
val channel: Int
get() = owner.channel
val time: Instant
get() = owner.startTime.plusNanos(timeOffset)
get() = owner.startTime.plusNanos(timeOffset.toLong())
}

View File

@ -159,7 +159,7 @@ class ProtoBlock(
.error("The block is broken. Number of times is ${events.timesCount} and number of amplitudes is ${events.amplitudesCount}")
}
IntStream.range(0, events.timesCount)
.mapToObj { i -> NumassEvent(events.getAmplitudes(i).toShort(), events.getTimes(i), this) }
.mapToObj { i -> NumassEvent(events.getAmplitudes(i).toUShort(), events.getTimes(i), this) }
} else {
Stream.empty()
}

View File

@ -69,7 +69,7 @@ class ChernovProcessor(
val timeInTicks = (pos + buffer.position() - 1)
val event = OrphanNumassEvent(amp.toInt().toShort(), (timeInTicks * tickSize).toLong())
val event = OrphanNumassEvent(amp.toInt().toUShort(), (timeInTicks * tickSize).toLong())
yield(event)
//subtracting event from buffer copy

View File

@ -90,13 +90,13 @@ object NumassDataUtils {
}
}
suspend fun NumassBlock.transformChain(transform: (NumassEvent, NumassEvent) -> Pair<Short, Long>?): NumassBlock {
suspend fun NumassBlock.transformChain(transform: (NumassEvent, NumassEvent) -> OrphanNumassEvent?): NumassBlock {
return SimpleBlock.produce(this.startTime, this.length) {
this.events.asSequence()
.sortedBy { it.timeOffset }
.zipWithNext(transform)
.filterNotNull()
.map { OrphanNumassEvent(it.first, it.second) }.asIterable()
.asIterable()
}
}

View File

@ -98,15 +98,15 @@ open class TimeAnalyzer(processor: SignalProcessor? = null) : AbstractAnalyzer(p
.forEach { pair ->
totalN.incrementAndGet()
//TODO add progress listener here
totalT.addAndGet(pair.second)
totalT.addAndGet(pair.second.toLong())
}
if (totalN.toInt() == 0) {
error("Zero number of intervals")
}
val countRate =
1e6 * totalN.get() / (totalT.get() / 1000 - t0 * totalN.get() / 1000)//1e9 / (totalT.get() / totalN.get() - t0);
val countRate = 1e6 * totalN.get() / (totalT.get() / 1000 - t0.toLong() * totalN.get() / 1000)
//1e9 / (totalT.get() / totalN.get() - t0);
val countRateError = countRate / sqrt(totalN.get().toDouble())
val length = totalT.get() / 1e9
val count = (length * countRate).toLong()
@ -252,7 +252,7 @@ open class TimeAnalyzer(processor: SignalProcessor? = null) : AbstractAnalyzer(p
* @return
*/
override fun getEvents(block: NumassBlock, meta: Meta): List<NumassEvent> {
val t0 = getT0(block, meta).toLong()
val t0 = getT0(block, meta)
return getEventsWithDelay(block, meta)
.filter { pair -> pair.second >= t0 }
.map { it.first }.toList()

View File

@ -182,8 +182,8 @@ constructor(override val name: String, private val path: Path, meta: Meta) : Num
@Throws(IOException::class)
private fun readEvent(channel: SeekableByteChannel, b: Int, timeDiv: Double): OrphanNumassEvent {
val chanel: Short
val time: Long
val chanel: UShort
val time: ULong
val hb = b and 0x0f
val lab = b and 0xf0
@ -194,34 +194,34 @@ constructor(override val name: String, private val path: Path, meta: Meta) : Num
time = readTime(channel)
}
0x20 -> {
chanel = 0
chanel = 0U
time = readTime(channel)
}
0x40 -> {
time = 0
time = 0U
chanel = readChanel(channel, hb)
}
0x80 -> {
time = 0
chanel = 0
time = 0U
chanel = 0U
}
else -> throw IOException("Event head expected")
}
return OrphanNumassEvent(chanel, (time / timeDiv).toLong())
return OrphanNumassEvent(chanel, (time.toDouble() / timeDiv).toLong())
}
@Throws(IOException::class)
private fun readChanel(channel: SeekableByteChannel, hb: Int): Short {
private fun readChanel(channel: SeekableByteChannel, hb: Int): UShort {
assert(hb < 127)
val buffer = readBlock(channel, 1)
return (buffer.get() + 256 * hb).toShort()
return (buffer.get() + 256 * hb).toUShort()
}
@Throws(IOException::class)
private fun readTime(channel: SeekableByteChannel): Long {
private fun readTime(channel: SeekableByteChannel): ULong {
val rx = readBlock(channel, 4)
return rx.long//rx[0] + 256 * rx[1] + 65536 * rx[2] + 256 * 65536 * rx[3];
return rx.long.toULong()//rx[0] + 256 * rx[1] + 65536 * rx[2] + 256 * 65536 * rx[3];
}
// private void skip(int length) throws IOException {

View File

@ -87,7 +87,7 @@ class ClassicNumassPoint(private val envelope: Envelope) : NumassPoint {
}
override fun next(): NumassEvent {
val amp = java.lang.Short.toUnsignedInt(buffer.short).toShort()
val amp = buffer.short.toUShort()
val time = Integer.toUnsignedLong(buffer.int)
val status = buffer.get() // status is ignored
return NumassEvent(amp, (time * timeCoef).toLong(), this@ClassicBlock)

View File

@ -1,48 +0,0 @@
package inr.numass.data
import groovy.transform.CompileStatic
import hep.dataforge.grind.Grind
import hep.dataforge.maths.histogram.Histogram
import hep.dataforge.maths.histogram.UnivariateHistogram
import hep.dataforge.values.Values
import inr.numass.data.analyzers.TimeAnalyzer
import inr.numass.data.api.NumassBlock
import java.util.stream.LongStream
/**
* Created by darksnake on 27-Jun-17.
*/
@CompileStatic
class PointAnalyzer {
static final TimeAnalyzer analyzer = new TimeAnalyzer();
static Histogram histogram(NumassBlock point, int loChannel = 0, int upChannel = 10000, double binSize = 0.5, int binNum = 500) {
return UnivariateHistogram.buildUniform(0d, binSize * binNum, binSize)
.fill(
kotlin.streams.jdk8.StreamsKt.asStream(analyzer.getEventsWithDelay(point, Grind.buildMeta("window.lo": loChannel, "window.up": upChannel))).mapToDouble {
it.second / 1000 as double
}
)
}
static Histogram histogram(LongStream stream, double binSize = 0.5, int binNum = 500) {
return UnivariateHistogram.buildUniform(0d, binSize * binNum, binSize).fill(stream.mapToDouble {
it / 1000 as double
})
}
static Values analyze(Map values = Collections.emptyMap(), NumassBlock block, Closure metaClosure = null) {
return analyzer.analyze(block, Grind.buildMeta(values, metaClosure))
}
//
// static class Result {
// double cr;
// double crErr;
// long num;
// double t0;
// int loChannel;
// int upChannel;
// }
}

View File

@ -90,7 +90,7 @@ class NumassPlugin : BasicPlugin() {
private fun loadMath(math: FunctionLibrary) {
math.addBivariate("numass.trap.lowFields") { Ei, Ef -> 3.92e-5 * FastMath.exp(-(Ei - Ef) / 300.0) + 1.97e-4 - 6.818e-9 * Ei }
math.addBivariate("numass.trap.nominal") { Ei, Ef ->
math.addBivariate("numass.trap.nominal") { Ei, f ->
//return 1.64e-5 * FastMath.exp(-(Ei - Ef) / 300d) + 1.1e-4 - 4e-9 * Ei;
1.2e-4 - 4.5e-9 * Ei
}

View File

@ -61,7 +61,7 @@ object TimeSpectrumAction : OneToOneAction<NumassPoint, Table>( "numass.timeSpec
.fill(analyzer
.getEventsWithDelay(input, inputMeta)
.asStream()
.mapToDouble { it.second / 1000.0 }
.mapToDouble { it.second.toDouble() / 1000.0 }
).asTable()
//.histogram(input, loChannel, upChannel, binSize, binNum).asTable();

View File

@ -24,7 +24,7 @@ private fun RandomGenerator.nextDeltaTime(cr: Double): Long {
}
suspend fun Sequence<OrphanNumassEvent>.generateBlock(start: Instant, length: Long): NumassBlock {
return SimpleBlock.produce(start, Duration.ofNanos(length)) {
return SimpleBlock.produce(start, Duration.ofNanos(length.toLong())) {
takeWhile { it.timeOffset < length }.toList()
}
}
@ -43,7 +43,7 @@ private class MergingState(private val chains: List<Chain<OrphanNumassEvent>>) {
* Merge event chains in ascending time order
*/
fun List<Chain<OrphanNumassEvent>>.merge(): Chain<OrphanNumassEvent> {
return StatefulChain(MergingState(this), OrphanNumassEvent(0, 0)) {
return StatefulChain(MergingState(this), OrphanNumassEvent(0U, 0)) {
poll()
}
}
@ -64,8 +64,8 @@ fun Chain<OrphanNumassEvent>.withDeadTime(deadTime: (OrphanNumassEvent) -> Long)
object NumassGenerator {
val defaultAmplitudeGenerator: RandomGenerator.(OrphanNumassEvent?, Long) -> Short =
{ _, _ -> ((nextDouble() + 2.0) * 100).toInt().toShort() }
val defaultAmplitudeGenerator: RandomGenerator.(OrphanNumassEvent?, Long) -> UShort =
{ _, _ -> ((nextDouble() + 2.0) * 100).toInt().toUShort() }
/**
* Generate an event chain with fixed count rate
@ -77,7 +77,7 @@ object NumassGenerator {
fun generateEvents(
cr: Double,
rnd: RandomGenerator = defaultGenerator,
amp: RandomGenerator.(OrphanNumassEvent?, Long) -> Short = defaultAmplitudeGenerator): Chain<OrphanNumassEvent> {
amp: RandomGenerator.(OrphanNumassEvent?, Long) -> UShort = defaultAmplitudeGenerator): Chain<OrphanNumassEvent> {
return MarkovChain(OrphanNumassEvent(rnd.amp(null, 0), 0)) { event ->
val deltaT = rnd.nextDeltaTime(cr)
OrphanNumassEvent(rnd.amp(event, deltaT), event.timeOffset + deltaT)
@ -102,7 +102,7 @@ object NumassGenerator {
bunchRate: Double,
bunchLength: Double,
rnd: RandomGenerator = defaultGenerator,
amp: RandomGenerator.(OrphanNumassEvent?, Long) -> Short = defaultAmplitudeGenerator
amp: RandomGenerator.(OrphanNumassEvent?, Long) -> UShort = defaultAmplitudeGenerator
): Chain<OrphanNumassEvent> {
return StatefulChain(
BunchState(0, 0),
@ -134,6 +134,6 @@ object NumassGenerator {
}
val distribution = EnumeratedRealDistribution(channels, values)
return generateEvents(cr, rnd) { _, _ -> distribution.sample().toInt().toShort() }
return generateEvents(cr, rnd) { _, _ -> distribution.sample().toInt().toUShort() }
}
}

View File

@ -32,6 +32,7 @@ import java.time.Duration
import java.time.Instant
import java.util.concurrent.atomic.AtomicInteger
import java.util.concurrent.atomic.AtomicReference
import kotlin.math.pow
/**
* @author [Alexander Nozik](mailto:altavir@gmail.com)
@ -43,6 +44,7 @@ class PileUpSimulator {
private val pileup = ArrayList<OrphanNumassEvent>()
private val registered = ArrayList<OrphanNumassEvent>()
private var generator: Chain<OrphanNumassEvent>
//private double uSet = 0;
private val doublePileup = AtomicInteger(0)
@ -65,15 +67,15 @@ class PileUpSimulator {
// }
fun generated(): NumassBlock {
return SimpleBlock(Instant.EPOCH, Duration.ofNanos(pointLength), generated)
return SimpleBlock(Instant.EPOCH, Duration.ofNanos(pointLength.toLong()), generated)
}
fun registered(): NumassBlock {
return SimpleBlock(Instant.EPOCH, Duration.ofNanos(pointLength), registered)
return SimpleBlock(Instant.EPOCH, Duration.ofNanos(pointLength.toLong()), registered)
}
fun pileup(): NumassBlock {
return SimpleBlock(Instant.EPOCH, Duration.ofNanos(pointLength), pileup)
return SimpleBlock(Instant.EPOCH, Duration.ofNanos(pointLength.toLong()), pileup)
}
/**
@ -82,7 +84,7 @@ class PileUpSimulator {
* @param x
* @return
*/
private fun pileupChannel(x: Double, prevChanel: Short, nextChanel: Short): Short {
private fun pileupChannel(x: Double, prevChanel: UShort, nextChanel: UShort): UShort {
assert(x > 0)
//эмпирическая формула для канала
val coef = max(0.0, 0.99078 + 0.05098 * x - 0.45775 * x * x + 0.10962 * x * x * x)
@ -90,7 +92,7 @@ class PileUpSimulator {
throw Error()
}
return (prevChanel + coef * nextChanel).toInt().toShort()
return (prevChanel.toDouble() + coef * nextChanel.toDouble()).toInt().toUShort()
}
/**
@ -110,9 +112,10 @@ class PileUpSimulator {
* @param delay
* @return
*/
private fun nextEventRegistered(prevChanel: Short, delay: Double): Boolean {
val average = 6.76102 - 4.31897E-4 * prevChanel + 7.88429E-8 * prevChanel.toDouble() * prevChanel.toDouble() + 0.2
val prob = 1.0 - 1.0 / (1.0 + Math.pow(delay / average, 75.91))
private fun nextEventRegistered(prevChanel: UShort, delay: Double): Boolean {
val average =
6.76102 - 4.31897E-4 * prevChanel.toDouble() + 7.88429E-8 * prevChanel.toDouble() * prevChanel.toDouble() + 0.2
val prob = 1.0 - 1.0 / (1.0 + (delay / average).pow(75.91))
return random(prob)
}
@ -124,7 +127,7 @@ class PileUpSimulator {
fun generate() {
var next: OrphanNumassEvent
//var lastRegisteredTime = 0.0 // Time of DAQ closing
val last = AtomicReference<OrphanNumassEvent>(OrphanNumassEvent(0, 0))
val last = AtomicReference(OrphanNumassEvent(0U, 0))
//flag that shows that previous event was pileup
var pileupFlag = false
@ -134,7 +137,8 @@ class PileUpSimulator {
generated.add(next)
//not counting double pileups
if (generated.size > 1) {
val delay = (next.timeOffset - last.get().timeOffset) / us //time between events in microseconds
val delay =
(next.timeOffset - last.get().timeOffset).toDouble() / us //time between events in microseconds
if (nextEventRegistered(next.amplitude, delay)) {
//just register new event
registered.add(next)

View File

@ -23,7 +23,8 @@ fun main() {
val regularChain = NumassGenerator.generateEvents(cr)
val bunchChain = NumassGenerator.generateBunches(40.0, 0.01, 5.0)
send(NumassGenerator.mergeEventChains(regularChain, bunchChain).generateBlock(start.plusNanos(it * length), length))
send(NumassGenerator.mergeEventChains(regularChain, bunchChain)
.generateBlock(start.plusNanos(it * length), length))
}
}

View File

@ -107,7 +107,7 @@ fun main() {
val sequence = TimeAnalyzer()
.getEventsWithDelay(point, metaForChainInverted)
.filter { pair -> pair.second <= t0 }
.filter { pair -> pair.second.toDouble() <= t0 }
.map { it.first }
val pileupSpectrum = sequence.getAmplitudeSpectrum(point.length.toMillis().toDouble() / 1000.0).withBinning(20)

View File

@ -109,7 +109,7 @@ fun main() {
val sequence = TimeAnalyzer()
.getEventsWithDelay(point, metaForChainInverted)
.filter { pair -> pair.second <= t0 }
.filter { pair -> pair.second.toDouble() <= t0 }
.map { it.first }
val pileupSpectrum = sequence.getAmplitudeSpectrum(point.length.toMillis().toDouble() / 1000.0).withBinning(20)

View File

@ -64,7 +64,7 @@ fun main() {
val histogram = SimpleHistogram(arrayOf(0.0, 0.0), arrayOf(20.0, 100.0))
histogram.fill(
TimeAnalyzer().getEventsWithDelay(point, meta)
.filter { it.second <10000 }
.filter { it.second < 10000 }
.filter { it.first.channel == 0 }
.map { arrayOf(it.first.amplitude.toDouble(), it.second.toDouble()/1e3) }
.asStream()

View File

@ -58,7 +58,7 @@ fun main() {
//
// val time = totalT.get()
println(time / 1e9)
println(time.toDouble() / 1e9)
//
// val cr = 80.0

View File

@ -4,6 +4,7 @@ import hep.dataforge.context.Global
import hep.dataforge.fx.output.FXOutputManager
import hep.dataforge.plots.jfreechart.JFreeChartPlugin
import inr.numass.data.ProtoNumassPoint
import inr.numass.data.api.OrphanNumassEvent
import inr.numass.data.plotAmplitudeSpectrum
import inr.numass.data.transformChain
import kotlinx.coroutines.runBlocking
@ -38,7 +39,7 @@ fun main() {
Global.plotAmplitudeSpectrum(point.transformChain { first, second ->
val dt = second.timeOffset - first.timeOffset
if (second.channel == 4 && first.channel == 0 && dt > window && dt < 1000) {
Pair((first.amplitude + second.amplitude).toShort(), second.timeOffset)
OrphanNumassEvent((first.amplitude + second.amplitude).toUShort(), second.timeOffset)
} else {
null
}
@ -55,7 +56,7 @@ fun main() {
Global.plotAmplitudeSpectrum(point.transformChain { first, second ->
val dt = second.timeOffset - first.timeOffset
if (second.channel == 0 && first.channel == 4 && dt > window && dt < 1000) {
Pair((first.amplitude + second.amplitude).toShort(), second.timeOffset)
OrphanNumassEvent((first.amplitude + second.amplitude).toUShort(), second.timeOffset)
} else {
null
}

View File

@ -36,8 +36,8 @@ fun main() {
val point: NumassPoint = set.points.first { it.index == 18 }
(0..99).forEach { bin ->
val times = point.events.filter { it.amplitude > 0 }.map { it.timeOffset }.toList()
val count = times.filter { it > bin.toDouble() / 10 * 1e9 && it < (bin + 1).toDouble() / 10 * 1e9 }.count()
val times = point.events.filter { it.amplitude > 0U }.map { it.timeOffset }.toList()
val count = times.count { it.toDouble() > bin.toDouble() / 10 * 1e9 && it.toDouble() < (bin + 1).toDouble() / 10 * 1e9 }
println("${bin.toDouble() / 10.0}: $count")
}
}

View File

@ -1,7 +1,6 @@
package inr.numass.scripts.tristan
import hep.dataforge.meta.Meta
import hep.dataforge.meta.value
import hep.dataforge.useValue
import inr.numass.data.analyzers.TimeAnalyzer
import inr.numass.data.api.NumassBlock
@ -9,11 +8,11 @@ import inr.numass.data.api.NumassEvent
object TristanAnalyzer : TimeAnalyzer() {
override fun getEvents(block: NumassBlock, meta: Meta): List<NumassEvent> {
val t0 = getT0(block, meta).toLong()
val t0 = getT0(block, meta)
val summTime = meta.getInt("summTime", 200) //time limit in nanos for event summation
var sequence = sequence {
var last: NumassEvent? = null
var amp: Int = 0
var amp = 0U
getEventsWithDelay(block, meta).forEach { (event, time) ->
when {
last == null -> {
@ -26,15 +25,15 @@ object TristanAnalyzer : TimeAnalyzer() {
}
time > t0 -> {
//accept new event and reset summator
if (amp != 0) {
if (amp != 0U) {
//construct new event with pileup
yield(NumassEvent(amp.toShort(), last!!.timeOffset, last!!.owner))
yield(NumassEvent(amp.toUShort(), last!!.timeOffset, last!!.owner))
} else {
//yield event without changes if there is no pileup
yield(last!!)
}
last = event
amp = event.amplitude.toInt()
amp = event.amplitude.toUInt()
}
//else ignore event
}