From 42e762c42184859cb1d8d035a9ef4750a0ef0ca8 Mon Sep 17 00:00:00 2001 From: darksnake Date: Mon, 3 Jul 2017 16:14:37 +0300 Subject: [PATCH] Histogram complete --- .../inr/numass/data/PointAnalyzer.groovy | 37 +++++++++++++---- .../numass/scripts/TestPointAnalyzer.groovy | 41 ++++++++++++++----- 2 files changed, 61 insertions(+), 17 deletions(-) diff --git a/numass-main/src/main/groovy/inr/numass/data/PointAnalyzer.groovy b/numass-main/src/main/groovy/inr/numass/data/PointAnalyzer.groovy index 8354c7ad..604fcbd9 100644 --- a/numass-main/src/main/groovy/inr/numass/data/PointAnalyzer.groovy +++ b/numass-main/src/main/groovy/inr/numass/data/PointAnalyzer.groovy @@ -4,6 +4,8 @@ import groovy.transform.CompileStatic import hep.dataforge.maths.histogram.Histogram import hep.dataforge.maths.histogram.UnivariateHistogram +import java.util.stream.DoubleStream + /** * Created by darksnake on 27-Jun-17. */ @@ -18,30 +20,51 @@ class PointAnalyzer { for (int i = 1; i < point.events.size(); i++) { NMEvent event = point.events[i]; double t = event.time - lastEvent.time; - if (t >= t0 && event.chanel <= upChannel && event.chanel >= loChannel) { + if (t < 0) { + lastEvent = event + } else if (t >= t0 && event.chanel <= upChannel && event.chanel >= loChannel) { totalN++ totalT += t + lastEvent = event } - lastEvent = event } double cr = 1d / (totalT / totalN - t0); return new Result(cr: cr, crErr: cr / Math.sqrt(totalN), num: totalN, t0: t0, loChannel: loChannel, upChannel: upChannel) } - - static Histogram histogram(RawNMPoint point, int loChannel = 0, int upChannel = 4000) { + private static DoubleStream timeChain(RawNMPoint point, int loChannel = 0, int upChannel = 4000) { List ts = new ArrayList<>(); NMEvent lastEvent = point.events[0]; for (int i = 1; i < point.events.size(); i++) { NMEvent event = point.events[i]; double t = event.time - lastEvent.time; - if (t >= 0 && event.chanel <= upChannel && event.chanel >= loChannel) { + if (t < 0) { + lastEvent = event + } else if (t >= 0 && event.chanel <= upChannel && event.chanel >= loChannel) { ts << t + lastEvent = event } - lastEvent = event } - return UnivariateHistogram.buildUniform(0d, 5e-4, 1e-6).fill(ts.stream().mapToDouble { it }) + return ts.stream().mapToDouble { it } + } + + /** + * Calculate the number of events in chain with delay and channel in given regions + * @param point + * @param t1 + * @param t2 + * @param loChannel + * @param upChannel + * @return + */ + static long count(RawNMPoint point, double t1, double t2, int loChannel = 0, int upChannel = 4000) { + return timeChain(point, loChannel, upChannel).filter { it > t1 && it < t2 }.count(); + } + + + static Histogram histogram(RawNMPoint point, int loChannel = 0, int upChannel = 4000, double binSize = 1e-6d, int binNum = 500) { + return UnivariateHistogram.buildUniform(0d, binSize*binNum, binSize).fill(timeChain(point, loChannel, upChannel)) } static class Result { diff --git a/numass-main/src/main/groovy/inr/numass/scripts/TestPointAnalyzer.groovy b/numass-main/src/main/groovy/inr/numass/scripts/TestPointAnalyzer.groovy index 6985de7d..4175b225 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/TestPointAnalyzer.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/TestPointAnalyzer.groovy @@ -25,35 +25,56 @@ GrindShell shell = new GrindShell(ctx) shell.eval { PlotHelper plot = plots - File rootDir = new File("D:\\Work\\Numass\\data\\2017_05\\Fill_1") + File rootDir = new File("D:\\Work\\Numass\\data\\2017_05\\Fill_1C") NumassStorage storage = NumassStorageFactory.buildLocal(rootDir); + def set = "set_6" def hv = 15000; - def point = storage.provide("loader::set_5/rawPoint::$hv", RawNMPoint.class).get(); - def histogram = PointAnalyzer.histogram(point,1000,1300).asTable(); + def loChannel = 3000; + def upChannel = 3600; - plot.configure("histogram"){ - yAxis(type:"log") + def point = storage.provide("loader::$set/rawPoint::$hv", RawNMPoint.class).get(); + + def histogram = PointAnalyzer.histogram(point, loChannel, upChannel).asTable(); + + println "finished histogram calculation..." + + plot.configure("histogram") { + yAxis(type: "log") } - plot.plot(histogram, ["frame": "histogram","showLine": true, "showSymbol": false, "showErrors": false, "connectionType": "step"]){ + plot.plot(name: hv, frame: "histogram", showLine: true, showSymbol: false, showErrors: false, connectionType: "step", histogram, { adapter("x.value": "x", "y.value": "count") - } + }) + def trueCR = PointAnalyzer.analyzePoint(point, 30e-6, loChannel, upChannel).cr + + println "The expected count rate for 30 us delay is $trueCR" def t0 = (1..150).collect { 5.5e-6 + 2e-7 * it } - def plotPoints = t0.collect { - def result = PointAnalyzer.analyzePoint(point, it) + def statPlotPoints = t0.collect { + def result = PointAnalyzer.analyzePoint(point, it, loChannel, upChannel) ValueMap.fromMap("x.value": it, "y.value": result.cr, "y.err": result.crErr); } //def cr = t0.collect { PointAnalyzer.analyzePoint(point, it).cr } - plot.plot(plotPoints, ["name": hv]) + plot.plot(name: hv, frame: "stat-method", statPlotPoints) + + def delta = 5e-6 + def discrepancyPlotPoints = (1..20).collect { delta * it }.collect { + def t1 = it + def t2 = it + delta + def result = PointAnalyzer.count(point, t1, t2, loChannel, upChannel) - (Math.exp(- trueCR * t1) - Math.exp(- trueCR * t2)) * point.length * trueCR + ValueMap.fromMap("x.value": it + delta / 2, "y.value": result); + } + + plot.plot(name: hv, frame: "discrepancy", discrepancyPlotPoints) // plot.plot(title: "dead time", from: 5.5e-6, to: 2e-5) { point.cr * 1d / (1d - 6.55e-6 * point.cr) } + storage.close() } \ No newline at end of file