From b47e5d91a96e9ad6b61e35bcb5605f3e18847b92 Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Fri, 18 Nov 2016 20:12:46 +0300 Subject: [PATCH] Fixed mathematics for dead time --- .../inr/numass/scripts/Underflow.groovy | 2 +- .../inr/numass/actions/PrepareDataAction.java | 168 +++++++++++++----- .../inr/numass/utils/ExpressionUtils.java | 4 +- 3 files changed, 130 insertions(+), 44 deletions(-) diff --git a/numass-main/src/main/groovy/inr/numass/scripts/Underflow.groovy b/numass-main/src/main/groovy/inr/numass/scripts/Underflow.groovy index cb605dbd..a3b41f74 100644 --- a/numass-main/src/main/groovy/inr/numass/scripts/Underflow.groovy +++ b/numass-main/src/main/groovy/inr/numass/scripts/Underflow.groovy @@ -14,7 +14,7 @@ import inr.numass.utils.UnderflowCorrection //File dataDir = new File("D:\\Work\\Numass\\data\\2016_04\\T2_data\\Fill_2_2\\set_7_b2a3433e54010000") //File dataDir = new File("D:\\Work\\Numass\\data\\2016_04\\T2_data\\Fill_2_2\\set_6_e26d123e54010000") -File dataDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_1\\set_26") +File dataDir = new File("D:\\Work\\Numass\\data\\2016_10\\Fill_1\\set_28") if(!dataDir.exists()){ println "dataDir directory does not exist" } diff --git a/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java b/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java index fe2c283e..b0d3b0a1 100644 --- a/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java +++ b/numass-main/src/main/java/inr/numass/actions/PrepareDataAction.java @@ -16,31 +16,29 @@ package inr.numass.actions; import hep.dataforge.actions.OneToOneAction; +import hep.dataforge.description.NodeDef; import hep.dataforge.description.TypedActionDef; import hep.dataforge.description.ValueDef; import hep.dataforge.exceptions.ContentException; import hep.dataforge.io.ColumnedDataWriter; import hep.dataforge.io.XMLMetaWriter; -import hep.dataforge.io.reports.Logable; import hep.dataforge.meta.Laminate; import hep.dataforge.meta.Meta; import hep.dataforge.tables.*; import inr.numass.storage.NMPoint; import inr.numass.storage.NumassData; import inr.numass.storage.RawNMPoint; -import inr.numass.utils.TritiumUtils; -import inr.numass.utils.UnderflowCorrection; import java.io.OutputStream; import java.time.Instant; import java.util.ArrayList; import java.util.List; import java.util.function.Function; +import java.util.stream.Collectors; import static inr.numass.utils.TritiumUtils.evaluateExpression; /** - * * @author Darksnake */ @TypedActionDef(name = "prepareData", inputType = NumassData.class, outputType = Table.class) @@ -50,6 +48,7 @@ import static inr.numass.utils.TritiumUtils.evaluateExpression; @ValueDef(name = "deadTime", type = "[NUMBER, STRING]", def = "0", info = "Dead time in s. Could be an expression.") @ValueDef(name = "correction", info = "An expression to correct count number depending on potential `U`, point length `T` and point itself as `point`") +@NodeDef(name = "correction", multiple = true, target = "method::inr.numass.actions.PrepareDataAction.makeCorrection") public class PrepareDataAction extends OneToOneAction { public static String[] parnames = {"Uset", "Uread", "Length", "Total", "Window", "Corr", "CR", "CRerr", "Timestamp"}; @@ -67,40 +66,62 @@ public class PrepareDataAction extends OneToOneAction { int upper = meta.getInt("upperWindow", RawNMPoint.MAX_CHANEL - 1); - Function deadTimeFunction; + List corrections = new ArrayList<>(); if (meta.hasValue("deadTime")) { - deadTimeFunction = point -> evaluateExpression(point, meta.getString("deadTime")); - } else { - deadTimeFunction = point -> 0.0; + corrections.add(new DeadTimeCorrection(meta.getString("deadTime"))); + } + + if (meta.hasNode("correction")) { + corrections.addAll(meta.getNodes("correction").stream() + .map((Function) this::makeCorrection) + .collect(Collectors.toList())); + } + + if (meta.hasValue("correction")) { + final String correction = meta.getString("correction"); + corrections.add((point) -> evaluateExpression(point, correction)); } -// double bkg = source.meta().getDouble("background", this.meta().getDouble("background", 0)); List dataList = new ArrayList<>(); for (NMPoint point : dataFile.getNMPoints()) { long total = point.getEventsCount(); - double Uset = point.getUset(); - double Uread = point.getUread(); + double uset = point.getUset(); + double uread = point.getUread(); double time = point.getLength(); - int a = getLowerBorder(meta, Uset); + int a = getLowerBorder(meta, uset); int b = Math.min(upper, RawNMPoint.MAX_CHANEL); // count in window long wind = point.getCountInWindow(a, b); - // count rate after all corrections - double cr = TritiumUtils.countRateWithDeadTime(point, a, b, deadTimeFunction.apply(point)); - // count rate error after all corrections - double crErr = TritiumUtils.countRateWithDeadTimeErr(point, a, b, deadTimeFunction.apply(point)); + double correctionFactor = corrections.stream() + .mapToDouble(cor -> cor.corr(point)) + .reduce((d1, d2) -> d1 * d2).getAsDouble(); + double relativeCorrectionError = Math.sqrt( + corrections.stream() + .mapToDouble(cor -> cor.relativeErr(point)) + .reduce((d1, d2) -> d1 * d1 + d2 * d2).getAsDouble() + ); - double correctionFactor = correction(getReport(name), point, meta); +// // count rate after all corrections +// double cr = TritiumUtils.countRateWithDeadTime(point, a, b, deadTimeFunction.apply(point)); +// // count rate error after all corrections +// double crErr = TritiumUtils.countRateWithDeadTimeErr(point, a, b, deadTimeFunction.apply(point)); +// +// double correctionFactor = correction(getReport(name), point, meta); - cr = cr * correctionFactor; - crErr = crErr * correctionFactor; + double cr = wind / point.getLength() * correctionFactor; + double crErr; + if (relativeCorrectionError == 0) { + crErr = Math.sqrt(wind) / point.getLength() * correctionFactor; + } else { + crErr = Math.sqrt(1d / wind + Math.pow(relativeCorrectionError, 2)) * cr; + } Instant timestamp = point.getStartTime(); - dataList.add(new MapPoint(parnames, new Object[]{Uset, Uread, time, total, wind, correctionFactor, cr, crErr, timestamp})); + dataList.add(new MapPoint(parnames, new Object[]{uset, uread, time, total, wind, correctionFactor, cr, crErr, timestamp})); } TableFormat format; @@ -129,32 +150,97 @@ public class PrepareDataAction extends OneToOneAction { return data; } - @Override - protected String getResultName(String dataName, Meta actionMeta) { - return super.getResultName(dataName, actionMeta); + +// /** +// * The factor to correct for count below detector threshold +// * +// * @param log +// * @param point +// * @param meta +// * @return +// */ +// private double correction(Logable log, NMPoint point, Laminate meta) { +// if (meta.hasValue("correction")) { +//// log.report("Using correction from formula: {}", meta.getString("correction")); +// return evaluateExpression(point, meta.getString("correction")); +// } else if (meta.hasMeta("underflow")) { +// return new UnderflowCorrection().get(log, meta.getMeta("underflow"), point); +// } else { +// return 1; +// } +// } + + @ValueDef(name = "value", type = "[NUMBER, STRING]", info = "Value or function to multiply count rate") + @ValueDef(name = "err", type = "[NUMBER, STRING]", info = "error of the value") + private Correction makeCorrection(Meta corrMeta) { + final String expr = corrMeta.getString("value"); + final String errExpr = corrMeta.getString("err", ""); + return new Correction() { + @Override + public double corr(NMPoint point) { + return evaluateExpression(point, expr); + } + + @Override + public double corrErr(NMPoint point) { + if (errExpr.isEmpty()) { + return 0; + } else { + return evaluateExpression(point, errExpr); + } + } + }; } - /** - * The factor to correct for count below detector threshold - * - * @param log - * @param point - * @param meta - * @return - */ - private double correction(Logable log, NMPoint point, Laminate meta) { - if (meta.hasValue("correction")) { -// log.report("Using correction from formula: {}", meta.getString("correction")); - return evaluateExpression(point, meta.getString("correction")); - } else if (meta.hasMeta("underflow")) { - return new UnderflowCorrection().get(log, meta.getMeta("underflow"), point); - } else { - return 1; + + private interface Correction { + /** + * correction coefficient + * + * @param point + * @return + */ + double corr(NMPoint point); + + /** + * correction coefficient uncertainty + * + * @param point + * @return + */ + default double corrErr(NMPoint point) { + return 0; } - } - + default double relativeErr(NMPoint point) { + double corrErr = corrErr(point); + if (corrErr == 0) { + return 0; + } else { + return corrErr / corr(point); + } + } + } + private class DeadTimeCorrection implements Correction { + + private final Function deadTimeFunction; + + public DeadTimeCorrection(String expr) { + deadTimeFunction = point -> evaluateExpression(point, expr); + } + + @Override + public double corr(NMPoint point) { + double deadTime = deadTimeFunction.apply(point); + double factor = deadTime / point.getLength() * point.getEventsCount(); +// double total = point.getEventsCount(); +// double time = point.getLength(); +// return 1d/(1d - factor); + + return (1d - Math.sqrt(1d - 4d * factor)) / 2d / factor; + } + } } 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 424a7878..091cdcdb 100644 --- a/numass-main/src/main/java/inr/numass/utils/ExpressionUtils.java +++ b/numass-main/src/main/java/inr/numass/utils/ExpressionUtils.java @@ -18,8 +18,8 @@ import java.util.Map; * @author Alexander Nozik */ public class ExpressionUtils { - private static Map cache = Utils.getLRUCache(100); - private static GroovyShell shell; + private static final Map cache = Utils.getLRUCache(100); + private static final GroovyShell shell; static { // Add imports for script.