Histogram complete

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
             }
+        }
         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<Double> 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
             }
-        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 {

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()
 }