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 bc9bc309..a10f67c1 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/SimulatePileup.groovy @@ -18,27 +18,14 @@ File dataDir = new File("D:\\Work\\Numass\\data\\2016_04\\T2_data\\Fill_2_2\\set 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 +Meta config = new GrindMetaBuilder().config(lowerChannel: 400, upperChannel: 1700) +//println config NumassData data = NumassDataLoader.fromLocalDir(null, dataDir) -Map res = new PileupSimulationAction().simpleRun(data, config) +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; @@ -55,3 +42,14 @@ points.first().keySet().each{ } +//printing count rate in window +print "U\tLength\t" +print keys.collect{it+"[total]"}.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" + println keys.collect{res[it].getNMPoints().get(i).getCountInWindow(500,1700)}.join("\t") +} diff --git a/numass-main/src/main/java/inr/numass/NumassPlugin.java b/numass-main/src/main/java/inr/numass/NumassPlugin.java index 95cf884e..39028153 100644 --- a/numass-main/src/main/java/inr/numass/NumassPlugin.java +++ b/numass-main/src/main/java/inr/numass/NumassPlugin.java @@ -46,7 +46,6 @@ import inr.numass.models.EmpiricalLossSpectrum; import inr.numass.models.ExperimentalVariableLossSpectrum; import inr.numass.models.GaussSourceSpectrum; import inr.numass.models.GunSpectrum; -import inr.numass.models.LossCalculator; import inr.numass.models.ModularSpectrum; import inr.numass.models.NBkgSpectrum; import inr.numass.models.RangedNamedSetSpectrum; @@ -221,10 +220,10 @@ public class NumassPlugin extends BasicPlugin { //Intercept = 4.95745, B1 = -0.36879, B2 = 0.00827 //sp.setTrappingFunction((Ei,Ef)->LossCalculator.getTrapFunction().value(Ei, Ef)*(4.95745-0.36879*Ei+0.00827*Ei*Ei)); sp.setTrappingFunction((Ei, Ef) -> { - return 4.1e-5 * FastMath.exp(-(Ei - Ef) / 600d) + 7.2e-3 * (0.05349 - 4.36375E-6 * (Ei) + 1.09875E-10 * Ei * Ei); + return 7.12e-5 * FastMath.exp(-(Ei - Ef) / 350d) + 0.1564 * (0.00123-4.2e-8*Ei); }); context.getReport().report("Using folowing trapping formula: {}", - "4e-5 * FastMath.exp(-(Ei - Ef) / 550d) + 7.2e-3 * (0.05349 - 4.36375E-6 * (Ei) + 1.09875E-10 * Ei**2)"); + "7.12e-5 * FastMath.exp(-(Ei - Ef) / 350d) + 0.1564 * (0.00123-4.2e-8*Ei)"); NBkgSpectrum spectrum = new NBkgSpectrum(sp); return new XYModel(spectrum, getAdapter(an)); diff --git a/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java b/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java index efcde0a6..042e0386 100644 --- a/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java +++ b/numass-main/src/main/java/inr/numass/actions/PileupSimulationAction.java @@ -31,25 +31,35 @@ public class PileupSimulationAction extends OneToOneAction 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); + int lowerChannel = inputMeta.getInt("lowerChannel", 1); + int upperChannel = inputMeta.getInt("upperChannel", RawNMPoint.MAX_CHANEL-1); List generated = new ArrayList<>(); List registered = new ArrayList<>(); List firstIteration = new ArrayList<>(); + List secondIteration = new ArrayList<>(); List pileup = new ArrayList<>(); - double crScale = inputMeta.getDouble("crScale",1d); - + double scale = inputMeta.getDouble("scale", 1d); + input.getNMPoints().forEach(point -> { - PileUpSimulator simulator = new PileUpSimulator(crScale * point.getCountRate(lowerChannel, upperChannel, 0), point.getLength()) - .withGenerator(point, lowerChannel, upperChannel) + double length = point.getLength() * scale; + double cr = point.getCountRate(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(crScale * point.getCountRate(lowerChannel, upperChannel, 0), point.getLength()) + 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(); @@ -62,6 +72,7 @@ public class PileupSimulationAction extends OneToOneAction { public static String[] parnames = {"Uset", "Uread", "Length", "Total", "Window", "Corr", "CR", "CRerr", "Timestamp"}; @@ -79,9 +81,14 @@ public class PrepareDataAction extends OneToOneAction { int upper = meta.getInt("upperWindow", RawNMPoint.MAX_CHANEL - 1); - double deadTime = meta.getDouble("deadTime", 0); -// double bkg = source.meta().getDouble("background", this.meta().getDouble("background", 0)); + Function deadTimeFunction; + if (meta.hasValue("deadTime")) { + deadTimeFunction = point -> evaluate(point, meta.getString("deadTime")); + } else { + deadTimeFunction = point -> 0.0; + } +// double bkg = source.meta().getDouble("background", this.meta().getDouble("background", 0)); List dataList = new ArrayList<>(); for (NMPoint point : dataFile.getNMPoints()) { @@ -96,11 +103,11 @@ public class PrepareDataAction extends OneToOneAction { long wind = point.getCountInWindow(a, b); // count rate after all corrections - double cr = point.getCountRate(a, b, deadTime); + double cr = point.getCountRate(a, b, deadTimeFunction.apply(point)); // count rate error after all corrections - double crErr = point.getCountRateErr(a, b, deadTime); + double crErr = point.getCountRateErr(a, b, deadTimeFunction.apply(point)); - double correctionFactor = getUnderflowCorrection(log, point, meta, cr); + double correctionFactor = correction(log, point, meta); cr = cr * correctionFactor; crErr = crErr * correctionFactor; @@ -136,6 +143,22 @@ public class PrepareDataAction extends OneToOneAction { return data; } + /** + * Evaluate groovy expression to number + * + * @param point + * @param expression + * @param countRate + * @return + */ + private double evaluate(NMPoint point, String expression) { + Map exprParams = new HashMap<>(); + exprParams.put("T", point.getLength()); + exprParams.put("U", point.getUread()); + exprParams.put("point", point); + return ExpressionUtils.evaluate(expression, exprParams); + } + /** * The factor to correct for count below detector threshold * @@ -145,30 +168,30 @@ public class PrepareDataAction extends OneToOneAction { * @param countRate precalculated count rate in main window * @return */ - private double getUnderflowCorrection(Reportable log, NMPoint point, Laminate meta, double countRate) { - if (meta.hasNode("underflow") || meta.getBoolean("underflow", true)) { - if (point.getUset() > meta.getDouble("underflow.thresholdU", 17000)) { - if (meta.hasValue("underflow.correctionFunction")) { - Map exprParams = new HashMap<>(); - exprParams.put("U", point.getUread()); - exprParams.put("CR", countRate); - exprParams.put("T", point.getLength()); - return ExpressionUtils.evaluate(meta.getString("underflow.correctionFunction"), exprParams); + private double correction(Reportable log, NMPoint point, Laminate meta) { + if (meta.hasValue("correction")) { +// log.report("Using correction from formula: {}", meta.getString("correction")); + return evaluate(point, meta.getString("correction")); + } else if (meta.hasNode("underflow")) { + if (point.getUset() >= meta.getDouble("underflow.threshold", 17000)) { +// log.report("Using underflow factor from formula: {}", meta.getString("underflow.function")); + if (meta.hasValue("underflow.function")) { + return evaluate(point, meta.getString("underflow.function")); } else { return 1; } } else { try { - log.report("Calculating underflow correction coefficient for point {}", point.getUset()); +// log.report("Calculating underflow correction coefficient for point {}", point.getUset()); int xLow = meta.getInt("underflow.lowerBorder", meta.getInt("lowerWindow")); int xHigh = meta.getInt("underflow.upperBorder", 800); int binning = meta.getInt("underflow.binning", 20); int upper = meta.getInt("upperWindow", RawNMPoint.MAX_CHANEL - 1); long norm = point.getCountInWindow(xLow, upper); double[] fitRes = getUnderflowExpParameters(point, xLow, xHigh, binning); - log.report("Underflow interpolation function: {}*exp(c/{})", fitRes[0], fitRes[1]); +// log.report("Underflow interpolation function: {}*exp(c/{})", fitRes[0], fitRes[1]); double correction = fitRes[0] * fitRes[1] * (Math.exp(xLow / fitRes[1]) - 1d) / norm + 1d; - log.report("Underflow correction factor: {}", correction); +// log.report("Underflow correction factor: {}", correction); return correction; } catch (Exception ex) { log.reportError("Failed to calculate underflow parameters for point {} with message:", point.getUset(), ex.getMessage()); @@ -196,7 +219,7 @@ public class PrepareDataAction extends OneToOneAction { List points = point.getMapWithBinning(binning, false) .entrySet().stream() .filter(entry -> entry.getKey() >= xLow && entry.getKey() <= xHigh) - .map(p -> new WeightedObservedPoint(1d / p.getValue(), p.getKey(), p.getValue()/ binning)) + .map(p -> new WeightedObservedPoint(1d / p.getValue(), p.getKey(), p.getValue() / binning)) .collect(Collectors.toList()); SimpleCurveFitter fitter = SimpleCurveFitter.create(new ExponentFunction(), new double[]{1d, 200d}); return fitter.fit(points); diff --git a/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java b/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java index 71b1aa73..1ed099ef 100644 --- a/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/BetaSpectrum.java @@ -41,7 +41,9 @@ public class BetaSpectrum extends AbstractParametricFunction implements RangedNa public BetaSpectrum(File FSSFile) { super(list); - this.fss = new FSS(FSSFile); + if (FSSFile != null) { + this.fss = new FSS(FSSFile); + } } double derivRoot(int n, double E0, double mnu2, double E) throws NotDefinedException { diff --git a/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java b/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java index fb581163..4bf123a0 100644 --- a/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java +++ b/numass-main/src/main/java/inr/numass/models/ModularSpectrum.java @@ -212,9 +212,9 @@ public class ModularSpectrum extends AbstractParametricFunction { this.caching = caching; this.trappingCache.setCachingEnabled(caching); - for (NamedSpectrumCaching sp : this.cacheList) { + this.cacheList.stream().forEach((sp) -> { sp.setCachingEnabled(caching); - } + }); } /** diff --git a/numass-main/src/main/java/inr/numass/models/SterileNeutrinoSpectrum.java b/numass-main/src/main/java/inr/numass/models/SterileNeutrinoSpectrum.java new file mode 100644 index 00000000..22018981 --- /dev/null +++ b/numass-main/src/main/java/inr/numass/models/SterileNeutrinoSpectrum.java @@ -0,0 +1,41 @@ +/* + * 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.models; + +import hep.dataforge.functions.AbstractParametricFunction; +import hep.dataforge.functions.ParametricFunction; +import hep.dataforge.values.NamedValueSet; + +/** + * + * @author Alexander Nozik + */ +public class SterileNeutrinoSpectrum extends AbstractParametricFunction { + + private static final String[] list = {"X", "trap", "E0", "mnu2", "msterile2", "U2"}; + private ParametricFunction spectrum; + + public SterileNeutrinoSpectrum() { + super(list); + } + + + @Override + public double derivValue(String parName, double x, NamedValueSet set) { + throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates. + } + + @Override + public double value(double x, NamedValueSet set) { + throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates. + } + + @Override + public boolean providesDeriv(String name) { + throw new UnsupportedOperationException("Not supported yet."); //To change body of generated methods, choose Tools | Templates. + } + +} diff --git a/numass-main/src/main/java/inr/numass/utils/ExpressionUtils.java b/numass-main/src/main/java/inr/numass/utils/ExpressionUtils.java index 95140165..1c2349a4 100644 --- a/numass-main/src/main/java/inr/numass/utils/ExpressionUtils.java +++ b/numass-main/src/main/java/inr/numass/utils/ExpressionUtils.java @@ -17,7 +17,7 @@ import org.codehaus.groovy.control.customizers.ImportCustomizer; */ public class ExpressionUtils { - public static Double evaluate(String expression, Map binding) { + public static Double evaluate(String expression, Map binding) { Binding b = new Binding(binding); // Add imports for script. ImportCustomizer importCustomizer = new ImportCustomizer(); 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 27e825dc..f97bf756 100644 --- a/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java +++ b/numass-main/src/main/java/inr/numass/utils/NMEventGenerator.java @@ -72,35 +72,43 @@ public class NMEventGenerator { for (Map.Entry entry : spectrum.entrySet()) { chanels[i] = entry.getKey(); values[i] = entry.getValue(); - i++; } distribution = new EnumeratedRealDistribution(chanels, values); } - public void loadSpectrum(NMPoint point, int minChanel, int maxChanel) { - assert minChanel >= 0; - assert maxChanel <= RawNMPoint.MAX_CHANEL; - + public void loadSpectrum(NMPoint point) { 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] = point.getCountInChanel(i); - i++; } distribution = new EnumeratedRealDistribution(chanels, values); } - public void loadSpectrum(NMPoint point, NMPoint reference, int minChanel, int maxChanel) { - assert minChanel >= 0; - assert maxChanel <= RawNMPoint.MAX_CHANEL; - + public void loadSpectrum(NMPoint point, NMPoint reference) { 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); + } + + /** + * + * @param point + * @param reference + * @param lower lower channel for spectrum generation + * @param upper upper channel for spectrum generation + */ + public void loadSpectrum(NMPoint point, NMPoint reference, int lower, int upper) { + double[] chanels = new double[RawNMPoint.MAX_CHANEL]; + double[] values = new double[RawNMPoint.MAX_CHANEL]; + for (int i = lower; i < upper; i++) { + chanels[i] = i; + values[i] = Math.max(0, point.getCountInChanel(i) - (reference == null ? 0 : reference.getCountInChanel(i))); } distribution = new EnumeratedRealDistribution(chanels, values); } 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 b28ecfe6..4d42d5c9 100644 --- a/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java +++ b/numass-main/src/main/java/inr/numass/utils/PileUpSimulator.java @@ -32,15 +32,21 @@ public class PileUpSimulator { this.pointLength = length; } - public PileUpSimulator withGenerator(NMPoint spectrum, NMPoint reference, int lowerChannel, int upperChannel) { + public PileUpSimulator withGenerator(NMPoint spectrum, NMPoint reference) { this.uSet = spectrum.getUset(); - generator.loadSpectrum(spectrum, reference, lowerChannel, upperChannel); + generator.loadSpectrum(spectrum, reference); return this; } - public PileUpSimulator withGenerator(NMPoint spectrum, int lowerChannel, int upperChannel) { + public PileUpSimulator withGenerator(NMPoint spectrum, NMPoint reference, int from, int to) { this.uSet = spectrum.getUset(); - generator.loadSpectrum(spectrum, lowerChannel, upperChannel); + generator.loadSpectrum(spectrum, reference, from, to); + return this; + } + + public PileUpSimulator withGenerator(NMPoint spectrum) { + this.uSet = spectrum.getUset(); + generator.loadSpectrum(spectrum); return this; } @@ -102,18 +108,22 @@ public class PileUpSimulator { public synchronized PileUpSimulator generate() { NMEvent current = null; + boolean pileupFlag = false; while (true) { NMEvent next = generator.nextEvent(current); if (next.getTime() > pointLength) { break; } generated.add(next.clone()); + //flag that shows that previous event was pileup + //not counting double pileups 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)) { + pileupFlag = false; + } else if (pileup(delay) && !pileupFlag) { //pileup event short newChannel = pileupChannel(delay, current.getChanel(), next.getChanel()); NMEvent newEvent = new NMEvent(newChannel, current.getTime()); @@ -121,6 +131,7 @@ public class PileUpSimulator { registred.remove(registred.size() - 1); registred.add(newEvent.clone()); pileup.add(newEvent.clone()); + pileupFlag = true; } else { // second event not registered } diff --git a/numass-main/src/test/java/inr/numass/actions/PrepareDataActionTest.java b/numass-main/src/test/java/inr/numass/actions/PrepareDataActionTest.java index bdbf40cd..7b83d3bb 100644 --- a/numass-main/src/test/java/inr/numass/actions/PrepareDataActionTest.java +++ b/numass-main/src/test/java/inr/numass/actions/PrepareDataActionTest.java @@ -23,7 +23,7 @@ public class PrepareDataActionTest { @Test public void testExpression() { - Map exprParams = new HashMap<>(); + Map exprParams = new HashMap<>(); exprParams.put("U", 18000d); double correctionFactor = ExpressionUtils.evaluate("1 + 13.265 * exp(- U / 2343.4)", exprParams); Assert.assertEquals("Testing expression calculation", 1.006125, correctionFactor, 1e-5);