[no commit message]

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

@@ -1,18 +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.scripts
-
-inr.numass.context.Main.main("-lc")

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

@@ -1,52 +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.scripts;
-
-import hep.dataforge.plots.jfreechart.JFreeChartFrame;
-import java.io.File;
-import java.io.FileInputStream;
-import java.io.IOException;
-import java.io.ObjectInputStream;
-import org.slf4j.LoggerFactory;
-
-/**
- * A basic plot viewer for serialized JFreeChart plots
- *
- * @author Darksnake
- */
-public class PlotViewer {
-
-    /**
-     * @param args the command line arguments
-     */
-    public static void main(String[] args) {
-        if (args.length > 0) {
-            for (String arg : args) {
-                File ser = new File(arg);
-                try {
-                    FileInputStream stream = new FileInputStream(ser);
-                    ObjectInputStream ostr = new ObjectInputStream(stream);
-                    JFreeChartFrame.deserialize(ostr);
-                } catch (IOException ex) {
-                    LoggerFactory.getLogger(PlotViewer.class).error("IO error during deserialization", ex);
-                } catch (ClassNotFoundException ex) {
-                    LoggerFactory.getLogger(PlotViewer.class).error("Wrong serialized content type", ex);
-                }
-            }
-        }
-    }
-
-}

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

@@ -0,0 +1,57 @@
+/*
+ * 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.scripts
+
+import inr.numass.storage.NMFile
+import inr.numass.storage.NMPoint
+import inr.numass.storage.NumassData
+import inr.numass.storage.NumassDataLoader
+import hep.dataforge.meta.Meta
+import inr.numass.actions.PileupSimulationAction
+import hep.dataforge.grind.GrindMetaBuilder
+
+File dataDir = new File("D:\\Work\\Numass\\data\\2016_04\\T2_data\\Fill_2_2\\set_6_e26d123e54010000")
+if(!dataDir.exists()){
+    println "dataDir directory does not exist"
+}
+int lowerChannel = 400
+int upperChannel = 1700
+
+Meta config = new GrindMetaBuilder().config(lowerChannel: lowerChannel, upperChannel: upperChannel)
+println config
+NumassData data = NumassDataLoader.fromLocalDir(null, dataDir)
+Map<String, NumassData> res = new PileupSimulationAction().simpleRun(data, config)
+
+def keys = res.keySet();
+
+//printing count rate in window
+print "U\t"
+print "Length\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"
+    println keys.collect{res[it].getNMPoints().get(i).getCountInWindow(500,1700)}.join("\t")
+}
+
+//print spectra for selected point
+double u = 15000d;
+
+List<Map> points = res.collect{key, value -> value.getByUset(u).getMapWithBinning(20, false)}
+
+println "\n Spectrum example for U = ${u}\n"
+
+print "channel\t"
+println keys.join("\t")
+
+points.first().keySet().each{
+    print "${it}\t"
+    println points.collect{map-> map[it]}.join("\t")
+}
+
+

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

@@ -0,0 +1,111 @@
+/*
+ * 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.context.Context;
+import hep.dataforge.description.TypedActionDef;
+import hep.dataforge.io.reports.Reportable;
+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 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(Context context, Reportable log, String name, Laminate inputMeta, NumassData input) {
+        int lowerChannel = inputMeta.getInt("lowerChannel", 0);
+        int upperChannel = inputMeta.getInt("upperChannel", RawNMPoint.MAX_CHANEL);
+
+        List<NMPoint> generated = new ArrayList<>();
+        List<NMPoint> registered = new ArrayList<>();
+        List<NMPoint> firstIteration = new ArrayList<>();
+        List<NMPoint> pileup = new ArrayList<>();
+
+        double crScale = inputMeta.getDouble("crScale",1d);
+        
+        input.getNMPoints().forEach(point -> {
+            PileUpSimulator simulator = new PileUpSimulator(crScale * point.getCountRate(lowerChannel, upperChannel, 0), point.getLength())
+                    .withGenerator(point, lowerChannel, upperChannel)
+                    .generate();
+
+            //second iteration to exclude pileup overlap
+            NMPoint pileupPoint = simulator.pileup();
+            firstIteration.add(simulator.registered());
+            simulator = new PileUpSimulator(crScale * point.getCountRate(lowerChannel, upperChannel, 0), point.getLength())
+                    .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("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/BunchGenerator.java

@@ -13,7 +13,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package inr.numass.generators;
+package inr.numass.utils;
 
 import inr.numass.storage.NMEvent;
 import inr.numass.storage.RawNMPoint;

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

@@ -13,7 +13,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package inr.numass.generators;
+package inr.numass.utils;
 
 import inr.numass.storage.NMEvent;
 import inr.numass.storage.RawNMPoint;

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

@@ -13,7 +13,7 @@
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
-package inr.numass.generators;
+package inr.numass.utils;
 
 import inr.numass.storage.NMEvent;
 import inr.numass.storage.NMPoint;
@@ -63,8 +63,8 @@ public class NMEventGenerator {
     }
 
     public void loadSpectrum(Map<Double, Double> spectrum, int minChanel, int maxChanel) {
-        assert minChanel > 0;
+        assert minChanel >= 0;
-        assert maxChanel < RawNMPoint.MAX_CHANEL;
+        assert maxChanel <= RawNMPoint.MAX_CHANEL;
 
         double[] chanels = new double[spectrum.size()];
         double[] values = new double[spectrum.size()];
@@ -78,8 +78,8 @@ public class NMEventGenerator {
     }
 
     public void loadSpectrum(NMPoint point, int minChanel, int maxChanel) {
-        assert minChanel > 0;
+        assert minChanel >= 0;
-        assert maxChanel < RawNMPoint.MAX_CHANEL;
+        assert maxChanel <= RawNMPoint.MAX_CHANEL;
 
         double[] chanels = new double[RawNMPoint.MAX_CHANEL];
         double[] values = new double[RawNMPoint.MAX_CHANEL];
@@ -91,6 +91,20 @@ public class NMEventGenerator {
         distribution = new EnumeratedRealDistribution(chanels, values);
     }
 
+    public void loadSpectrum(NMPoint point, NMPoint reference, int minChanel, int maxChanel) {
+        assert minChanel >= 0;
+        assert maxChanel <= RawNMPoint.MAX_CHANEL;
+
+        double[] chanels = new double[RawNMPoint.MAX_CHANEL];
+        double[] values = new double[RawNMPoint.MAX_CHANEL];
+        for (int i = 0; i < RawNMPoint.MAX_CHANEL; i++) {
+            chanels[i] = i;
+            values[i] = Math.max(0, point.getCountInChanel(i) - reference.getCountInChanel(i));
+            i++;
+        }
+        distribution = new EnumeratedRealDistribution(chanels, values);
+    }
+
     public NMEvent nextEvent(NMEvent prev) {
         short chanel;
 
@@ -100,15 +114,15 @@ public class NMEventGenerator {
             chanel = 1600;
         }
 
-        return new NMEvent(chanel, prev.getTime() + nextExp(1 / cr));
+        return new NMEvent(chanel, prev == null ? 0 : prev.getTime() + nextExp(1 / cr));
     }
 
-    double nextExp(double mean) {
+    public double nextExp(double mean) {
         double rand = this.nextUniform();
         return -mean * Math.log(1 - rand);
     }
 
-    double nextPositiveGaussian(double mean, double sigma) {
+    public double nextPositiveGaussian(double mean, double sigma) {
         double res = -1;
         while (res <= 0) {
             res = mean + generator.nextGaussian() * sigma;
@@ -116,11 +130,11 @@ public class NMEventGenerator {
         return res;
     }
 
-    double nextUniform() {
+    public double nextUniform() {
         return generator.nextDouble();
     }
 
-    void setSeed(int seed) {
+    public void setSeed(int seed) {
         generator.setSeed(seed);
     }
 

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

@@ -5,11 +5,12 @@
  */
 package inr.numass.utils;
 
-import inr.numass.generators.NMEventGenerator;
 import inr.numass.storage.NMEvent;
+import inr.numass.storage.NMPoint;
+import inr.numass.storage.RawNMPoint;
+import static java.lang.Math.max;
 import java.util.ArrayList;
 import java.util.List;
-import java.util.function.Function;
 
 /**
  *
@@ -17,12 +18,119 @@ import java.util.function.Function;
  */
 public class PileUpSimulator {
 
-    private NMEventGenerator generator;
+    private final static double us = 1e-6;//microsecond
+
+    private double uSet = 0;
+    private final double pointLength;
+    private final NMEventGenerator generator;
     private final List<NMEvent> generated = new ArrayList<>();
     private final List<NMEvent> pileup = new ArrayList<>();
     private final List<NMEvent> registred = new ArrayList<>();
 
-    private Function<Double, Double> firstEventSurvivalProb;
+    public PileUpSimulator(double countRate, double length) {
-    private Function<Double, Double> secondEventSurvivalProb;
+        generator = new NMEventGenerator(countRate);
+        this.pointLength = length;
+    }
+
+    public PileUpSimulator withGenerator(NMPoint spectrum, NMPoint reference, int lowerChannel, int upperChannel) {
+        this.uSet = spectrum.getUset();
+        generator.loadSpectrum(spectrum, reference, lowerChannel, upperChannel);
+        return this;
+    }
+
+    public PileUpSimulator withGenerator(NMPoint spectrum, int lowerChannel, int upperChannel) {
+        this.uSet = spectrum.getUset();
+        generator.loadSpectrum(spectrum, lowerChannel, upperChannel);
+        return this;
+    }
+
+    public NMPoint generated() {
+        return new NMPoint(new RawNMPoint(uSet, generated, pointLength));
+    }
+
+    public NMPoint registered() {
+        return new NMPoint(new RawNMPoint(uSet, registred, pointLength));
+    }
+
+    public NMPoint pileup() {
+        return new NMPoint(new RawNMPoint(uSet, pileup, pointLength));
+    }
+
+    /**
+     * The amplitude for pileup event
+     *
+     * @param x
+     * @return
+     */
+    private short pileupChannel(double x, short prevChanel, short nextChanel) {
+        assert x > 0;
+        //эмпирическая формула для канала
+        double coef = max(0, 0.99078 + 0.05098 * x - 0.45775 * x * x + 0.10962 * x * x * x);
+        if (coef < 0 || coef > 1) {
+            throw new Error();
+        }
+
+        return (short) (prevChanel + coef * nextChanel);
+    }
+
+    /**
+     * pileup probability
+     *
+     * @param delay
+     * @return
+     */
+    private boolean pileup(double delay) {
+        double prob = 1d / (1d + Math.pow(delay / 2.5, 42.96));
+        return random(prob);
+    }
+
+    /**
+     * Probability for next event to register
+     *
+     * @param delay
+     * @return
+     */
+    private boolean nextEventRegistered(double delay) {
+        double prob = 1d - 1d / (1d + Math.pow(delay / 6.2, 75.91));
+        return random(prob);
+    }
+
+    private boolean random(double prob) {
+        double r = generator.nextUniform();
+        return r < prob;
+    }
+
+    public synchronized PileUpSimulator generate() {
+        NMEvent current = null;
+        while (true) {
+            NMEvent next = generator.nextEvent(current);
+            if (next.getTime() > pointLength) {
+                break;
+            }
+            generated.add(next.clone());
+            if (current != null) {
+                double delay = (next.getTime() - current.getTime()) / us;
+                if (nextEventRegistered(delay)) {
+                    //just register new event
+                    registred.add(next.clone());
+                } else if (pileup(delay)) {
+                    //pileup event
+                    short newChannel = pileupChannel(delay, current.getChanel(), next.getChanel());
+                    NMEvent newEvent = new NMEvent(newChannel, current.getTime());
+                    //replace already registered event by event with new channel
+                    registred.remove(registred.size() - 1);
+                    registred.add(newEvent.clone());
+                    pileup.add(newEvent.clone());
+                } else {
+                    // second event not registered
+                }
+            } else {
+                //register first event
+                registred.add(next.clone());
+            }
+            current = next;
+        }
+        return this;
+    }
 
 }

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

@@ -186,7 +186,6 @@ public class NMPoint {
                 i++;
             }
             data.add(new MapPoint(dataNames, start + binning / 2d, sum));
-            i++;
         }
         return data;
     }
@@ -211,7 +210,6 @@ public class NMPoint {
                 i++;
             }
             res.put(start + binning / 2d, sum);
-            i++;
         }
         return res;