diff --git a/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy b/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy index 1bf4f9f3..87dc8026 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy +++ b/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 secondIteration = new ArrayList<>(); List 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()); diff --git a/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java b/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java index 1e55cfc6..b4c2f98f 100644 --- a/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java +++ b/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 { - double cr; - RealDistribution distribution; - private final RandomGenerator rnd; - private NMEvent prevEvent; + protected double cr; + protected RealDistribution distribution; + protected final RandomGenerator rnd; + 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 shorts = new ArrayList<>(); point.getEvents().stream() @@ -112,7 +118,7 @@ public class NMEventGenerator implements Supplier { 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 { 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); -// } - } diff --git a/numass-main/src/main/java/inr/numass/utils/NMEventGeneratorWithPulser.java b/numass-main/src/main/java/inr/numass/utils/NMEventGeneratorWithPulser.java new file mode 100644 index 00000000..67a55e03 --- /dev/null +++ b/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); + } +} 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 89826a32..a8dd2d67 100644 --- a/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java +++ b/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; }