From 0d21135c5922e1b453acb8cfaa1c9e685cc0c28a Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 24 Jun 2016 14:29:34 +0300 Subject: [PATCH] [no commit message] --- .../inr/numass/scripts/ListActions.groovy | 18 --- .../inr/numass/scripts/PlotViewer.groovy | 52 -------- .../inr/numass/scripts/SimulatePileup.groovy | 57 +++++++++ .../actions/PileupSimulationAction.java | 111 ++++++++++++++++ .../{generators => utils}/BunchGenerator.java | 2 +- .../EventChainGenerator.java | 2 +- .../NMEventGenerator.java | 34 +++-- .../inr/numass/utils/PileUpSimulator.java | 120 +++++++++++++++++- .../main/java/inr/numass/storage/NMPoint.java | 2 - 9 files changed, 308 insertions(+), 90 deletions(-) delete mode 100644 numass-main/src/main/groovy/inr/numass/scripts/ListActions.groovy delete mode 100644 numass-main/src/main/groovy/inr/numass/scripts/PlotViewer.groovy create mode 100644 numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy create mode 100644 numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java rename numass-main/src/main/java/inr/numass/{generators => utils}/BunchGenerator.java (96%) rename numass-main/src/main/java/inr/numass/{generators => utils}/EventChainGenerator.java (96%) rename numass-main/src/main/java/inr/numass/{generators => utils}/NMEventGenerator.java (74%) diff --git a/numass-main/src/main/groovy/inr/numass/scripts/ListActions.groovy b/numass-main/src/main/groovy/inr/numass/scripts/ListActions.groovy deleted file mode 100644 index a92d45d2..00000000 --- a/numass-main/src/main/groovy/inr/numass/scripts/ListActions.groovy +++ /dev/null @@ -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") diff --git a/numass-main/src/main/groovy/inr/numass/scripts/PlotViewer.groovy b/numass-main/src/main/groovy/inr/numass/scripts/PlotViewer.groovy deleted file mode 100644 index d09049c5..00000000 --- a/numass-main/src/main/groovy/inr/numass/scripts/PlotViewer.groovy +++ /dev/null @@ -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); - } - } - } - } - -} diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy new file mode 100644 index 00000000..bc9bc309 --- /dev/null +++ b/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 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 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") +} + + diff --git a/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java b/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java new file mode 100644 index 00000000..efcde0a6 --- /dev/null +++ b/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> { + + @Override + protected Map 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 generated = new ArrayList<>(); + List registered = new ArrayList<>(); + List firstIteration = new ArrayList<>(); + List 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 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 points; + + public SimulatedPoint(String name, List points) { + this.name = name; + this.points = points; + } + + @Override + public String getDescription() { + return name; + } + + @Override + public Meta meta() { + return Meta.empty(); + } + + @Override + public List getNMPoints() { + return points; + } + + @Override + public boolean isEmpty() { + return points.isEmpty(); + } + + @Override + public Instant startTime() { + return Instant.EPOCH; + } + + @Override + public String getName() { + return name; + } + + } + +} diff --git a/numass-main/src/main/java/inr/numass/generators/BunchGenerator.java b/numass-main/src/main/java/inr/numass/utils/BunchGenerator.java similarity index 96% rename from numass-main/src/main/java/inr/numass/generators/BunchGenerator.java rename to numass-main/src/main/java/inr/numass/utils/BunchGenerator.java index f9ebe4da..463ceb02 100644 --- a/numass-main/src/main/java/inr/numass/generators/BunchGenerator.java +++ b/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; diff --git a/numass-main/src/main/java/inr/numass/generators/EventChainGenerator.java b/numass-main/src/main/java/inr/numass/utils/EventChainGenerator.java similarity index 96% rename from numass-main/src/main/java/inr/numass/generators/EventChainGenerator.java rename to numass-main/src/main/java/inr/numass/utils/EventChainGenerator.java index 78a5201b..ab42e046 100644 --- a/numass-main/src/main/java/inr/numass/generators/EventChainGenerator.java +++ b/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; diff --git a/numass-main/src/main/java/inr/numass/generators/NMEventGenerator.java b/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java similarity index 74% rename from numass-main/src/main/java/inr/numass/generators/NMEventGenerator.java rename to numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java index e8d453b7..27e825dc 100644 --- a/numass-main/src/main/java/inr/numass/generators/NMEventGenerator.java +++ b/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 spectrum, int minChanel, int maxChanel) { - assert minChanel > 0; - assert maxChanel < RawNMPoint.MAX_CHANEL; + assert minChanel >= 0; + 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 maxChanel < RawNMPoint.MAX_CHANEL; + assert minChanel >= 0; + 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); } diff --git a/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java b/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java index 5c3a83c2..b28ecfe6 100644 --- a/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java +++ b/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 generated = new ArrayList<>(); private final List pileup = new ArrayList<>(); private final List registred = new ArrayList<>(); - private Function firstEventSurvivalProb; - private Function secondEventSurvivalProb; - + public PileUpSimulator(double countRate, double length) { + 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; + } + } diff --git a/numass-storage/src/main/java/inr/numass/storage/NMPoint.java b/numass-storage/src/main/java/inr/numass/storage/NMPoint.java index dacd142f..448b3e3a 100644 --- a/numass-storage/src/main/java/inr/numass/storage/NMPoint.java +++ b/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;