grind and groovymath update

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 double cr;
-    RealDistribution distribution;
+    protected RealDistribution distribution;
-    private final RandomGenerator rnd;
+    protected final RandomGenerator rnd;
-    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

@@ -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;
     }