pileup generator with pulser

numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy

@@ -6,10 +6,11 @@
 
 package inr.numass.scripts
 
+import hep.dataforge.grind.Grind
 import inr.numass.storage.NMPoint
 import inr.numass.storage.NumassData
 import inr.numass.storage.NumassDataLoader
-import inr.numass.utils.NMEventGenerator
+import inr.numass.utils.NMEventGeneratorWithPulser
 import inr.numass.utils.PileUpSimulator
 import inr.numass.utils.TritiumUtils
 import org.apache.commons.math3.random.JDKRandomGenerator
@@ -33,14 +34,23 @@ List<NMPoint> secondIteration = new ArrayList<>();
 List<NMPoint> pileup = new ArrayList<>();
 
 lowerChannel = 400;
-upperChannel = 3800;
+upperChannel = 1800;
 
 PileUpSimulator buildSimulator(NMPoint point, double cr, NMPoint reference = null, double scale = 1d) {
-    NMEventGenerator generator = new NMEventGenerator(cr, rnd)
+    def cfg = Grind.buildMeta(cr: cr) {
+        pulser(mean: 3450, sigma: 86.45, freq: 66.43)
+    }
+    NMEventGeneratorWithPulser generator = new NMEventGeneratorWithPulser(rnd, cfg)
     generator.loadSpectrum(point, reference, lowerChannel, upperChannel);
     return new PileUpSimulator(point.length * scale, rnd, generator).withUset(point.uset).generate();
 }
 
+double adjustCountRate(PileUpSimulator simulator, NMPoint point) {
+    double generatedInChannel = simulator.generated().getCountInWindow(lowerChannel, upperChannel);
+    double registeredInChannel = simulator.registered().getCountInWindow(lowerChannel, upperChannel);
+    return (generatedInChannel / registeredInChannel) * (point.getCountInWindow(lowerChannel, upperChannel) / point.getLength());
+}
+
 data.NMPoints.forEach { point ->
     double cr = TritiumUtils.countRateWithDeadTime(point, lowerChannel, upperChannel, 6.2e-6);
 
@@ -49,11 +59,15 @@ data.NMPoints.forEach { point ->
     //second iteration to exclude pileup overlap
     NMPoint pileupPoint = simulator.pileup();
     firstIteration.add(simulator.registered());
+
+    //updating count rate
+    cr = adjustCountRate(simulator, point);
     simulator = buildSimulator(point, cr, pileupPoint);
 
     pileupPoint = simulator.pileup();
     secondIteration.add(simulator.registered());
 
+    cr = adjustCountRate(simulator, point);
     simulator = buildSimulator(point, cr, pileupPoint);
 
     generated.add(simulator.generated());

numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java

@@ -15,6 +15,7 @@
  */
 package inr.numass.utils;
 
+import hep.dataforge.meta.Meta;
 import inr.numass.storage.NMEvent;
 import inr.numass.storage.NMPoint;
 import inr.numass.storage.RawNMPoint;
@@ -35,16 +36,21 @@ import java.util.function.Supplier;
  */
 public class NMEventGenerator implements Supplier<NMEvent> {
 
-    double cr;
+    protected final RandomGenerator rnd;
-    RealDistribution distribution;
+    protected double cr;
-    private final RandomGenerator rnd;
+    protected RealDistribution distribution;
-    private NMEvent prevEvent;
+    protected NMEvent prevEvent;
 
-    public NMEventGenerator(double cr, RandomGenerator rnd) {
+    public NMEventGenerator(RandomGenerator rnd, double cr) {
         this.cr = cr;
         this.rnd = rnd;
     }
 
+    public NMEventGenerator(RandomGenerator rnd, Meta meta) {
+        this.cr = meta.getDouble("cr");
+        this.rnd = rnd;
+    }
+
     public void loadSpectrum(RawNMPoint point, int minChanel, int maxChanel) {
         List<Short> shorts = new ArrayList<>();
         point.getEvents().stream()
@@ -112,7 +118,7 @@ public class NMEventGenerator implements Supplier<NMEvent> {
         distribution = new EnumeratedRealDistribution(chanels, values);
     }
 
-    private NMEvent nextEvent(NMEvent prev) {
+    protected NMEvent nextEvent(NMEvent prev) {
         short chanel;
 
         if (distribution != null) {
@@ -134,21 +140,4 @@ public class NMEventGenerator implements Supplier<NMEvent> {
         return -mean * Math.log(1 - rnd.nextDouble());
     }
 
-
-//    public double nextPositiveGaussian(double mean, double sigma) {
-//        double res = -1;
-//        while (res <= 0) {
-//            res = mean + generator.nextGaussian() * sigma;
-//        }
-//        return res;
-//    }
-//
-//    public double nextUniform() {
-//        return generator.nextDouble();
-//    }
-//
-//    public void setSeed(int seed) {
-//        generator.setSeed(seed);
-//    }
-
 }

numass-main/src/main/java/inr/numass/utils/NMEventGeneratorWithPulser.java

@@ -0,0 +1,57 @@
+package inr.numass.utils;
+
+import hep.dataforge.meta.Meta;
+import inr.numass.storage.NMEvent;
+import org.apache.commons.math3.distribution.NormalDistribution;
+import org.apache.commons.math3.distribution.RealDistribution;
+import org.apache.commons.math3.random.RandomGenerator;
+
+/**
+ * Created by darksnake on 25-Nov-16.
+ */
+public class NMEventGeneratorWithPulser extends NMEventGenerator {
+    private RealDistribution pulserChanelDistribution;
+    private double pulserDist;
+    private NMEvent pulserEvent;
+    private NMEvent nextEvent;
+
+    public NMEventGeneratorWithPulser(RandomGenerator rnd, Meta meta) {
+        super(rnd, meta);
+        pulserChanelDistribution = new NormalDistribution(
+                meta.getDouble("pulser.mean", 3450),
+                meta.getDouble("pulser.sigma", 86.45)
+        );
+        pulserDist = 1d / meta.getDouble("pulser.freq", 66.43);
+        pulserEvent = generatePulserEvent();
+    }
+
+    @Override
+    public synchronized NMEvent get() {
+        //expected next event
+        if (nextEvent == null) {
+            nextEvent = nextEvent(prevEvent);
+        }
+        //if pulser event is first, then leave next event as is and return pulser event
+        if (pulserEvent.getTime() < nextEvent.getTime()) {
+            NMEvent res = pulserEvent;
+            pulserEvent = generatePulserEvent();
+            return res;
+        } else {
+            //else return saved next event and generate next one
+            prevEvent = nextEvent;
+            nextEvent = nextEvent(prevEvent);
+            return prevEvent;
+        }
+    }
+
+    private NMEvent generatePulserEvent() {
+        short channel = (short) pulserChanelDistribution.sample();
+        double time;
+        if (pulserEvent == null) {
+            time = rnd.nextDouble() * pulserDist;
+        } else {
+            time = pulserEvent.getTime() + pulserDist;
+        }
+        return new NMEvent(channel, time);
+    }
+}

numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java

@@ -23,11 +23,11 @@ public class PileUpSimulator {
 
     private final static double us = 1e-6;//microsecond
     private final double pointLength;
-    private Supplier<NMEvent> generator;
     private final RandomGenerator rnd;
     private final List<NMEvent> generated = new ArrayList<>();
     private final List<NMEvent> pileup = new ArrayList<>();
     private final List<NMEvent> registred = new ArrayList<>();
+    private Supplier<NMEvent> generator;
     private double uSet = 0;
 
     public PileUpSimulator(double length, RandomGenerator rnd, Supplier<NMEvent> sup) {
@@ -38,7 +38,7 @@ public class PileUpSimulator {
 
     public PileUpSimulator(double length, RandomGenerator rnd, double countRate) {
         this.rnd = rnd;
-        generator = new NMEventGenerator(countRate, rnd);
+        generator = new NMEventGenerator(rnd, countRate);
         this.pointLength = length;
     }
 

numass-storage/src/main/java/inr/numass/storage/NMPoint.java

@@ -30,13 +30,11 @@ import static java.util.Arrays.sort;
 public class NMPoint {
     //TODO transform to annotated and move some parameters to meta
     static final String[] dataNames = {"chanel", "count"};
+    private final int[] spectrum;
     private Instant startTime;
     private long eventsCount;
-
     private int overflow;
-
     private double pointLength;
-    private final int[] spectrum;
     private double uread;
     private double uset;
 
@@ -138,12 +136,12 @@ public class NMPoint {
     }
 
     /**
-     * Events count - overflow
+     * Events count + overflow
      *
      * @return
      */
     public long getEventsCount() {
-        return eventsCount - getOverflow();
+        return eventsCount + getOverflow();
     }
 
     public List<DataPoint> getData(int binning, boolean normalize) {