pileup redone

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

@@ -6,28 +6,73 @@
 
 package inr.numass.scripts
 
-import hep.dataforge.grind.GrindMetaBuilder
+import inr.numass.storage.NMPoint
-import hep.dataforge.meta.Meta
-import inr.numass.actions.PileupSimulationAction
 import inr.numass.storage.NumassData
 import inr.numass.storage.NumassDataLoader
+import inr.numass.utils.NMEventGenerator
+import inr.numass.utils.PileUpSimulator
+import inr.numass.utils.TritiumUtils
+import org.apache.commons.math3.random.JDKRandomGenerator
 
-File dataDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_1\\set_10")
+rnd = new JDKRandomGenerator();
+
+//Loading data
+File dataDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_1\\set_28")
 if (!dataDir.exists()) {
     println "dataDir directory does not exist"
 }
-
-Meta config = new GrindMetaBuilder().config(lowerChannel: 500, upperChannel: 1800)
-//println config
 NumassData data = NumassDataLoader.fromLocalDir(null, dataDir)
-Map<String, NumassData> res = new PileupSimulationAction().simpleRun(data,config)
+
+//Simulation process
+Map<String, List<NMPoint>> res = [:]
+
+List<NMPoint> generated = new ArrayList<>();
+List<NMPoint> registered = new ArrayList<>();
+List<NMPoint> firstIteration = new ArrayList<>();
+List<NMPoint> secondIteration = new ArrayList<>();
+List<NMPoint> pileup = new ArrayList<>();
+
+lowerChannel = 400;
+upperChannel = 3800;
+
+PileUpSimulator buildSimulator(NMPoint point, double cr, NMPoint reference = null, double scale = 1d) {
+    NMEventGenerator generator = new NMEventGenerator(cr, rnd)
+    generator.loadSpectrum(point, reference, lowerChannel, upperChannel);
+    return new PileUpSimulator(point.length * scale, rnd, generator).withUset(point.uset).generate();
+}
+
+data.NMPoints.forEach { point ->
+    double cr = TritiumUtils.countRateWithDeadTime(point, lowerChannel, upperChannel, 6.2e-6);
+
+    PileUpSimulator simulator = buildSimulator(point, cr);
+
+    //second iteration to exclude pileup overlap
+    NMPoint pileupPoint = simulator.pileup();
+    firstIteration.add(simulator.registered());
+    simulator = buildSimulator(point, cr, pileupPoint);
+
+    pileupPoint = simulator.pileup();
+    secondIteration.add(simulator.registered());
+
+    simulator = buildSimulator(point, cr, pileupPoint);
+
+    generated.add(simulator.generated());
+    registered.add(simulator.registered());
+    pileup.add(simulator.pileup());
+}
+res.put("original", data.NMPoints);
+res.put("generated", generated);
+res.put("registered", registered);
+//    res.put("firstIteration", new SimulatedPoint("firstIteration", firstIteration));
+//    res.put("secondIteration", new SimulatedPoint("secondIteration", secondIteration));
+res.put("pileup", pileup);
 
 def keys = res.keySet();
 
 //print spectra for selected point
 double u = 16500d;
 
-List<Map> points = res.collect{key, value -> value.getByUset(u).getMapWithBinning(20, false)}
+List<Map> points = res.values().collect { it.find { it.uset == u }.getMapWithBinning(20, false) }
 
 println "\n Spectrum example for U = ${u}\n"
 
@@ -39,15 +84,16 @@ points.first().keySet().each{
     println points.collect { map -> map[it] }.join("\t")
 }
 
-
 //printing count rate in window
 print "U\tLength\t"
 print keys.collect { it + "[total]" }.join("\t") + "\t"
+print keys.collect { it + "[pulse]" }.join("\t") + "\t"
 println keys.join("\t")
 
 for (int i = 0; i < data.getNMPoints().size(); i++) {
     print "${data.getNMPoints().get(i).getUset()}\t"
     print "${data.getNMPoints().get(i).getLength()}\t"
-    print keys.collect{res[it].getNMPoints().get(i).getEventsCount()}.join("\t") + "\t"
+    print keys.collect { res[it].get(i).getEventsCount() }.join("\t") + "\t"
-    println keys.collect { res[it].getNMPoints().get(i).getCountInWindow(500, 1800) }.join("\t")
+    print keys.collect { res[it].get(i).getCountInWindow(3100, 3800) }.join("\t") + "\t"
+    println keys.collect { res[it].get(i).getCountInWindow(400, 3100) }.join("\t")
 }

numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java

@@ -1,122 +0,0 @@
-/*
- * To change this license header, choose License Headers in Project Properties.
- * To change this template file, choose Tools | Templates
- * and open the template in the editor.
- */
-package inr.numass.actions;
-
-import hep.dataforge.actions.OneToOneAction;
-import hep.dataforge.description.TypedActionDef;
-import hep.dataforge.meta.Laminate;
-import hep.dataforge.meta.Meta;
-import inr.numass.storage.NMPoint;
-import inr.numass.storage.NumassData;
-import inr.numass.storage.RawNMPoint;
-import inr.numass.utils.PileUpSimulator;
-import inr.numass.utils.TritiumUtils;
-
-import java.time.Instant;
-import java.util.ArrayList;
-import java.util.LinkedHashMap;
-import java.util.List;
-import java.util.Map;
-
-/**
- * Simulate pileup
- *
- * @author Alexander Nozik
- */
-@TypedActionDef(name = "simulatePileup", inputType = NumassData.class, outputType = Map.class)
-public class PileupSimulationAction extends OneToOneAction<NumassData, Map<String, NumassData>> {
-
-    @Override
-    protected Map<String, NumassData> execute(String name, Laminate inputMeta, NumassData input) {
-        int lowerChannel = inputMeta.getInt("lowerChannel", 1);
-        int upperChannel = inputMeta.getInt("upperChannel", RawNMPoint.MAX_CHANEL - 1);
-
-        List<NMPoint> generated = new ArrayList<>();
-        List<NMPoint> registered = new ArrayList<>();
-        List<NMPoint> firstIteration = new ArrayList<>();
-        List<NMPoint> secondIteration = new ArrayList<>();
-        List<NMPoint> pileup = new ArrayList<>();
-
-        double scale = inputMeta.getDouble("scale", 1d);
-
-        input.getNMPoints().forEach(point -> {
-            double length = point.getLength() * scale;
-            double cr = TritiumUtils.countRateWithDeadTime(point, lowerChannel, upperChannel, 6.4e-6);
-
-            PileUpSimulator simulator = new PileUpSimulator(cr, length)
-                    .withGenerator(point, null, lowerChannel, upperChannel)
-                    .generate();
-
-            //second iteration to exclude pileup overlap
-            NMPoint pileupPoint = simulator.pileup();
-            firstIteration.add(simulator.registered());
-            simulator = new PileUpSimulator(cr, length)
-                    .withGenerator(point, pileupPoint, lowerChannel, upperChannel)
-                    .generate();
-
-            pileupPoint = simulator.pileup();
-            secondIteration.add(simulator.registered());
-            simulator = new PileUpSimulator(cr, length)
-                    .withGenerator(point, pileupPoint, lowerChannel, upperChannel)
-                    .generate();
-
-            generated.add(simulator.generated());
-            registered.add(simulator.registered());
-            pileup.add(simulator.pileup());
-        });
-        Map<String, NumassData> res = new LinkedHashMap<>();
-        res.put("original", input);
-        res.put("generated", new SimulatedPoint("generated", generated));
-        res.put("registered", new SimulatedPoint("registered", registered));
-        res.put("firstIteration", new SimulatedPoint("firstIteration", firstIteration));
-        res.put("secondIteration", new SimulatedPoint("secondIteration", secondIteration));
-        res.put("pileup", new SimulatedPoint("pileup", pileup));
-        return res;
-    }
-
-    private static class SimulatedPoint implements NumassData {
-
-        private final String name;
-        private final List<NMPoint> points;
-
-        public SimulatedPoint(String name, List<NMPoint> points) {
-            this.name = name;
-            this.points = points;
-        }
-
-        @Override
-        public String getDescription() {
-            return name;
-        }
-
-        @Override
-        public Meta meta() {
-            return Meta.empty();
-        }
-
-        @Override
-        public List<NMPoint> getNMPoints() {
-            return points;
-        }
-
-        @Override
-        public boolean isEmpty() {
-            return points.isEmpty();
-        }
-
-        @Override
-        public Instant startTime() {
-            return Instant.EPOCH;
-        }
-
-        @Override
-        public String getName() {
-            return name;
-        }
-
-    }
-
-}

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

@@ -1,184 +0,0 @@
-/* 
- * Copyright 2015 Alexander Nozik.
- *
- * Licensed under the Apache License, Version 2.0 (the "License");
- * you may not use this file except in compliance with the License.
- * You may obtain a copy of the License at
- *
- *      http://www.apache.org/licenses/LICENSE-2.0
- *
- * Unless required by applicable law or agreed to in writing, software
- * distributed under the License is distributed on an "AS IS" BASIS,
- * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
- * See the License for the specific language governing permissions and
- * limitations under the License.
- */
-package inr.numass.utils;
-
-import inr.numass.storage.NMEvent;
-import inr.numass.storage.RawNMPoint;
-import static java.lang.Math.max;
-import java.util.ArrayList;
-import java.util.List;
-import java.util.Map;
-import static java.lang.Math.max;
-
-/**
- * Utility class for pile-up simulation
- * @author Darksnake
- */
-public final class EventChainGenerator {
-
-    private final static double us = 1e-6;//microsecond
-
-    private double blockStartTime = 0;
-
-    private final List<NMEvent> generatedChain = new ArrayList<>();
-    private final NMEventGenerator generator;
-    private final double length;
-    private final List<NMEvent> pileupChain = new ArrayList<>();
-
-    private final List<NMEvent> registredChain = new ArrayList<>();
-
-    public EventChainGenerator(double cr, double length) {
-        generator = new NMEventGenerator(cr);
-        this.length = length;
-        
-        run();
-    }
-
-    public EventChainGenerator(double cr, double length, RawNMPoint source, int minChanel, int maxChanel) {
-        generator = new NMEventGenerator(cr);
-        this.generator.loadSpectrum(source, minChanel, maxChanel);
-        this.length = length;
-        
-        run();
-    }
-    
-    public EventChainGenerator(double cr, double length, Map<Double,Double> spectrum, int minChanel, int maxChanel) {
-        generator = new NMEventGenerator(cr);
-        this.generator.loadSpectrum(spectrum, minChanel, maxChanel);
-        this.length = length;
-        
-        run();
-    }    
-
-    /**
-     * Амлпитуда второго сигнала в зависимости от амплитуд наложенных сигналов и
-     * задержки
-     *
-     * @param delay
-     * @return
-     */
-    private short getNewChanel(double delay, short prevChanel, short newChanel) {
-        assert delay > 0;
-        //эмпирическая формула для канала
-        double x = delay / us;
-        double coef = max(0, 0.99078 + 0.05098 * x - 0.45775 * x * x + 0.10962 * x * x * x);
-        
-        return (short) (prevChanel + coef * newChanel);
-    }
-
-    public RawNMPoint getPileUp() {
-        return new RawNMPoint(2, pileupChain, length);
-    }
-
-    public RawNMPoint getPointAsGenerated() {
-        return new RawNMPoint(0, generatedChain, length);
-    }
-
-    public RawNMPoint getPointAsRegistred() {
-        return new RawNMPoint(1, registredChain, length);
-    }
-
-    /**
-     * Имеется второй сигнал
-     *
-     * @param delay
-     * @return
-     */
-    private boolean hasNew(double delay) {
-        if (delay > 2.65 * us) {
-            return false;
-        } else if (delay < 2.35 * us) {
-            return true;
-        } else {
-            return heads((2.65 * us - delay) / 0.3 / us);
-        }
-    }
-
-    private boolean heads(double prob) {
-        double r = generator.nextUniform();
-        return r < prob;
-    }
-
-    NMEvent nextEvent(NMEvent prev) {
-        if (prev == null) {
-            return generator.nextEvent(new NMEvent((short)0, 0));
-        }
-
-        NMEvent event = generator.nextEvent(prev);
-        generatedChain.add(event);
-        double delay = event.getTime() - blockStartTime;
-        if (notDT(delay)) {
-            //Если система сбора данных успела переварить предыдущие события
-            registredChain.add(event);
-            blockStartTime = event.getTime();
-            return event;
-        } else {
-            if ((!prevSurvived(delay)) && (!registredChain.isEmpty())) {
-                //если первое событие не выжило, а ушло в наложения
-                registredChain.remove(registredChain.size() - 1);
-            }
-            if (hasNew(delay)) {
-                // Если есть событие с увеличенной амлитудой
-                NMEvent pileup = new NMEvent(getNewChanel(delay, prev.getChanel(), event.getChanel()), event.getTime());
-                registredChain.add(pileup);
-                pileupChain.add(pileup);
-            }
-            //возвращаем предыдущее событие, чтобы отсчитывать мертвое время от него
-            return prev;
-        }
-    }
-
-    /**
-     * Не попал в мертвое время и наложения
-     *
-     * @param delay
-     * @return
-     */
-    private boolean notDT(double delay) {
-        if (delay > 7.0 * us) {
-            return true;
-        } else if (delay < 6.5 * us) {
-            return false;
-        } else {
-            return heads((delay - 6.5 * us) / 0.5 / us);
-        }
-    }
-
-    /**
-     * Выжило предыдушее событие
-     *
-     * @param delay
-     * @return
-     */
-    private boolean prevSurvived(double delay) {
-        if (delay > 2.65 * us) {
-            return true;
-        } else if (delay < 2.35 * us) {
-            return false;
-        } else {
-            return heads((delay - 2.35 * us) / 0.3 / us);
-        }
-    }
-
-    private void run() {
-        NMEvent next = null;
-        
-        do {
-            next = nextEvent(next);
-        } while (next.getTime() < length);
-    }
-
-}

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

@@ -18,31 +18,31 @@ package inr.numass.utils;
 import inr.numass.storage.NMEvent;
 import inr.numass.storage.NMPoint;
 import inr.numass.storage.RawNMPoint;
-import java.util.ArrayList;
-import java.util.List;
-import java.util.Map;
 import org.apache.commons.math3.distribution.EnumeratedRealDistribution;
 import org.apache.commons.math3.distribution.RealDistribution;
 import org.apache.commons.math3.random.EmpiricalDistribution;
-import org.apache.commons.math3.random.JDKRandomGenerator;
 import org.apache.commons.math3.random.RandomGenerator;
 
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Map;
+import java.util.function.Supplier;
+
 /**
  * A generator for Numass events with given energy spectrum
  *
  * @author Darksnake
  */
-public class NMEventGenerator {
+public class NMEventGenerator implements Supplier<NMEvent> {
 
     double cr;
-//    UnivariateFunction signalShape;
     RealDistribution distribution;
+    private final RandomGenerator rnd;
+    private NMEvent prevEvent;
 
-    private final RandomGenerator generator;
+    public NMEventGenerator(double cr, RandomGenerator rnd) {
-
-    public NMEventGenerator(double cr) {
         this.cr = cr;
-        generator = new JDKRandomGenerator();
+        this.rnd = rnd;
     }
 
     public void loadSpectrum(RawNMPoint point, int minChanel, int maxChanel) {
@@ -97,7 +97,6 @@ public class NMEventGenerator {
     }
 
     /**
-     *
      * @param point
      * @param reference
      * @param lower     lower channel for spectrum generation
@@ -113,7 +112,7 @@ public class NMEventGenerator {
         distribution = new EnumeratedRealDistribution(chanels, values);
     }
 
-    public NMEvent nextEvent(NMEvent prev) {
+    private NMEvent nextEvent(NMEvent prev) {
         short chanel;
 
         if (distribution != null) {
@@ -122,28 +121,34 @@ public class NMEventGenerator {
             chanel = 1600;
         }
 
-        return new NMEvent(chanel, prev == null ? 0 : prev.getTime() + nextExpDecay(1d / cr));
+        return new NMEvent(chanel, (prev == null ? 0 : prev.getTime()) + nextExpDecay(1d / cr));
     }
 
-    public double nextExpDecay(double mean) {
+
-        double rand = this.nextUniform();
+    @Override
-        return -mean * Math.log(1 - rand);
+    public synchronized NMEvent get() {
+        return prevEvent = nextEvent(prevEvent);
     }
 
-    public double nextPositiveGaussian(double mean, double sigma) {
+    private double nextExpDecay(double mean) {
-        double res = -1;
+        return -mean * Math.log(1 - rnd.nextDouble());
-        while (res <= 0) {
-            res = mean + generator.nextGaussian() * sigma;
-        }
-        return res;
     }
 
-    public double nextUniform() {
-        return generator.nextDouble();
-    }
 
-    public void setSeed(int seed) {
+//    public double nextPositiveGaussian(double mean, double sigma) {
-        generator.setSeed(seed);
+//        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/PileUpSimulator.java

@@ -8,9 +8,11 @@ package inr.numass.utils;
 import inr.numass.storage.NMEvent;
 import inr.numass.storage.NMPoint;
 import inr.numass.storage.RawNMPoint;
+import org.apache.commons.math3.random.RandomGenerator;
 
 import java.util.ArrayList;
 import java.util.List;
+import java.util.function.Supplier;
 
 import static java.lang.Math.max;
 
@@ -21,32 +23,32 @@ public class PileUpSimulator {
 
     private final static double us = 1e-6;//microsecond
     private final double pointLength;
-    private final NMEventGenerator generator;
+    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 double uSet = 0;
 
-    public PileUpSimulator(double countRate, double length) {
+    public PileUpSimulator(double length, RandomGenerator rnd, Supplier<NMEvent> sup) {
-        generator = new NMEventGenerator(countRate);
+        this.rnd = rnd;
+        generator = sup;//new NMEventGenerator(countRate, rnd);
         this.pointLength = length;
     }
 
-    public PileUpSimulator withGenerator(NMPoint spectrum, NMPoint reference) {
+    public PileUpSimulator(double length, RandomGenerator rnd, double countRate) {
-        this.uSet = spectrum.getUset();
+        this.rnd = rnd;
-        generator.loadSpectrum(spectrum, reference);
+        generator = new NMEventGenerator(countRate, rnd);
+        this.pointLength = length;
+    }
+
+    public PileUpSimulator withGenerator(Supplier<NMEvent> sup){
+        this.generator = sup;
         return this;
     }
 
-    public PileUpSimulator withGenerator(NMPoint spectrum, NMPoint reference, int from, int to) {
+    public PileUpSimulator withUset(double uset){
-        this.uSet = spectrum.getUset();
+        this.uSet = uset;
-        generator.loadSpectrum(spectrum, reference, from, to);
-        return this;
-    }
-
-    public PileUpSimulator withGenerator(NMPoint spectrum) {
-        this.uSet = spectrum.getUset();
-        generator.loadSpectrum(spectrum);
         return this;
     }
 
@@ -103,33 +105,32 @@ public class PileUpSimulator {
     }
 
     private boolean random(double prob) {
-        double r = generator.nextUniform();
+        return rnd.nextDouble() <= prob;
-        return r <= prob;
     }
 
     public synchronized PileUpSimulator generate() {
-        NMEvent last = null;// last event
+        NMEvent next;
         double lastRegisteredTime = 0; // Time of DAQ closing
         //flag that shows that previous event was pileup
         boolean pileupFlag = false;
         while (true) {
-            NMEvent next = generator.nextEvent(last);
+            next = generator.get();
             if (next.getTime() > pointLength) {
                 break;
             }
             generated.add(next);
             //not counting double pileups
-            if (last != null) {
+            if (generated.size() > 1) {
                 double delay = (next.getTime() - lastRegisteredTime) / us; //time between events in microseconds
-                if (nextEventRegistered(last.getChanel(), delay)) {
+                if (nextEventRegistered(next.getChanel(), delay)) {
                     //just register new event
                     registred.add(next);
                     lastRegisteredTime = next.getTime();
                     pileupFlag = false;
                 } else if (pileup(delay) && !pileupFlag) {
                     //pileup event
-                    short newChannel = pileupChannel(delay, last.getChanel(), next.getChanel());
+                    short newChannel = pileupChannel(delay, next.getChanel(), next.getChanel());
-                    NMEvent newEvent = new NMEvent(newChannel, last.getTime());
+                    NMEvent newEvent = new NMEvent(newChannel, next.getTime());
                     //replace already registered event by event with new channel
                     registred.remove(registred.size() - 1);
                     registred.add(newEvent);
@@ -145,7 +146,6 @@ public class PileUpSimulator {
                 registred.add(next);
                 lastRegisteredTime = next.getTime();
             }
-            last = next;
         }
         return this;
     }