From ea2aea0f542819b59d40b9d3729b2c0c3ca849d7 Mon Sep 17 00:00:00 2001 From: darksnake Date: Wed, 20 Jul 2016 16:04:42 +0300 Subject: [PATCH] [no commit message] --- .../inr/numass/scripts/ReadTextFile.groovy | 15 ++++ .../main/java/inr/numass/NumassPlugin.java | 2 + .../inr/numass/actions/MergeDataAction.java | 4 +- .../numass/actions/MonitorCorrectAction.java | 78 ++++++++++++++----- .../actions/SubstractSpectrumAction.java | 61 +++++++++++++++ 5 files changed, 138 insertions(+), 22 deletions(-) create mode 100644 numass-main/src/main/groovy/inr/numass/scripts/ReadTextFile.groovy create mode 100644 numass-main/src/main/java/inr/numass/actions/SubstractSpectrumAction.java diff --git a/numass-main/src/main/groovy/inr/numass/scripts/ReadTextFile.groovy b/numass-main/src/main/groovy/inr/numass/scripts/ReadTextFile.groovy new file mode 100644 index 00000000..229b6bfa --- /dev/null +++ b/numass-main/src/main/groovy/inr/numass/scripts/ReadTextFile.groovy @@ -0,0 +1,15 @@ +/* + * 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 hep.dataforge.io.ColumnedDataReader +import hep.dataforge.io.ColumnedDataWriter +import hep.dataforge.tables.Table + +File file = new File("D:\\Work\\Numass\\sterile2016\\empty.dat" ) +Table referenceTable = new ColumnedDataReader(file).toDataSet(); +ColumnedDataWriter.writeDataSet(System.out, referenceTable,"") \ No newline at end of file diff --git a/numass-main/src/main/java/inr/numass/NumassPlugin.java b/numass-main/src/main/java/inr/numass/NumassPlugin.java index 6dff00b6..b282d425 100644 --- a/numass-main/src/main/java/inr/numass/NumassPlugin.java +++ b/numass-main/src/main/java/inr/numass/NumassPlugin.java @@ -39,6 +39,7 @@ import inr.numass.actions.ReadNumassStorageAction; import inr.numass.actions.ShowEnergySpectrumAction; import inr.numass.actions.ShowLossSpectrumAction; import inr.numass.actions.SlicingAction; +import inr.numass.actions.SubstractSpectrumAction; import inr.numass.actions.SummaryAction; import inr.numass.models.BetaSpectrum; import inr.numass.models.CustomNBkgSpectrum; @@ -87,6 +88,7 @@ public class NumassPlugin extends BasicPlugin { actions.registerAction(AdjustErrorsAction.class); actions.registerAction(ReadNumassStorageAction.class); actions.registerAction(ShowEnergySpectrumAction.class); + actions.registerAction(SubstractSpectrumAction.class); } @Override diff --git a/numass-main/src/main/java/inr/numass/actions/MergeDataAction.java b/numass-main/src/main/java/inr/numass/actions/MergeDataAction.java index 38a0b03a..e2788c53 100644 --- a/numass-main/src/main/java/inr/numass/actions/MergeDataAction.java +++ b/numass-main/src/main/java/inr/numass/actions/MergeDataAction.java @@ -17,7 +17,6 @@ package inr.numass.actions; import hep.dataforge.actions.GroupBuilder; import hep.dataforge.actions.ManyToOneAction; -import hep.dataforge.context.Context; import hep.dataforge.data.DataNode; import hep.dataforge.description.NodeDef; import hep.dataforge.description.TypedActionDef; @@ -30,6 +29,7 @@ import hep.dataforge.tables.ListTable; import hep.dataforge.tables.MapPoint; import hep.dataforge.tables.PointSource; import hep.dataforge.tables.Table; +import hep.dataforge.tables.TableFormat; import java.io.OutputStream; import java.util.ArrayList; import java.util.Collection; @@ -169,7 +169,7 @@ public class MergeDataAction extends ManyToOneAction { res.add(curPoint); }); - return new ListTable(res); + return new ListTable(TableFormat.forNames(parnames),res); } diff --git a/numass-main/src/main/java/inr/numass/actions/MonitorCorrectAction.java b/numass-main/src/main/java/inr/numass/actions/MonitorCorrectAction.java index acf72378..6054d5b5 100644 --- a/numass-main/src/main/java/inr/numass/actions/MonitorCorrectAction.java +++ b/numass-main/src/main/java/inr/numass/actions/MonitorCorrectAction.java @@ -16,7 +16,6 @@ package inr.numass.actions; import hep.dataforge.actions.OneToOneAction; -import hep.dataforge.context.Context; import hep.dataforge.description.TypedActionDef; import hep.dataforge.description.ValueDef; import hep.dataforge.exceptions.ContentException; @@ -36,6 +35,9 @@ import java.util.List; import java.util.Map.Entry; import java.util.TreeMap; import java.util.concurrent.CopyOnWriteArrayList; +import javafx.util.Pair; +import org.apache.commons.math3.analysis.interpolation.SplineInterpolator; +import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; /** * @@ -81,28 +83,15 @@ public class MonitorCorrectAction extends OneToOneAction { MapPoint.Builder pb = new MapPoint.Builder(dp); pb.putValue("Monitor", 1.0); if (!isMonitorPoint(monitor, dp) || index.isEmpty()) { - Instant time = getTime(dp); - Entry previousMonitor = index.floorEntry(time); - Entry nextMonitor = index.ceilingEntry(time); - - if (previousMonitor == null) { - previousMonitor = nextMonitor; - } - - if (nextMonitor == null) { - nextMonitor = previousMonitor; - } - - double p; - if (nextMonitor.getKey().isAfter(time) && time.isAfter(previousMonitor.getKey())) { - p = 1.0 * (time.toEpochMilli() - previousMonitor.getKey().toEpochMilli()) - / (nextMonitor.getKey().toEpochMilli() - previousMonitor.getKey().toEpochMilli()); + Pair corr; + if (meta.getBoolean("spline", false)) { + corr = getSplineCorrection(index, dp, norm); } else { - p = 0.5; + corr = getLinearCorrection(index, dp, norm); } + double corrFactor = corr.getKey(); + double corrErr = corr.getValue(); - double corrFactor = (getCR(previousMonitor.getValue()) * (1 - p) + getCR(nextMonitor.getValue()) * p) / norm; - double corrErr = previousMonitor.getValue().getValue("CRerr").doubleValue() / getCR(previousMonitor.getValue()); double pointErr = dp.getValue("CRerr").doubleValue() / getCR(dp); double err = Math.sqrt(corrErr * corrErr + pointErr * pointErr) * getCR(dp); @@ -151,6 +140,55 @@ public class MonitorCorrectAction extends OneToOneAction { return data; } + private Pair getSplineCorrection(TreeMap index, DataPoint dp, double norm) { + double time = getTime(dp).toEpochMilli(); + + double[] xs = new double[index.size()]; + double[] ys = new double[index.size()]; + + int i = 0; + + for (Entry entry : index.entrySet()) { + xs[i] = (double) entry.getKey().toEpochMilli(); + ys[i] = getCR(entry.getValue()) / norm; + i++; + } + + PolynomialSplineFunction corrFunc = new SplineInterpolator().interpolate(xs, ys); + if (corrFunc.isValidPoint(time)) { + double averageErr = index.values().stream().mapToDouble(p -> p.getDouble("CRerr")).average().getAsDouble(); + return new Pair<>(corrFunc.value(time), averageErr / norm / 2d); + } else { + return new Pair<>(1d, 0d); + } + } + + private Pair getLinearCorrection(TreeMap index, DataPoint dp, double norm) { + Instant time = getTime(dp); + Entry previousMonitor = index.floorEntry(time); + Entry nextMonitor = index.ceilingEntry(time); + + if (previousMonitor == null) { + previousMonitor = nextMonitor; + } + + if (nextMonitor == null) { + nextMonitor = previousMonitor; + } + + double p; + if (nextMonitor.getKey().isAfter(time) && time.isAfter(previousMonitor.getKey())) { + p = 1.0 * (time.toEpochMilli() - previousMonitor.getKey().toEpochMilli()) + / (nextMonitor.getKey().toEpochMilli() - previousMonitor.getKey().toEpochMilli()); + } else { + p = 0.5; + } + + double corrFactor = (getCR(previousMonitor.getValue()) * (1 - p) + getCR(nextMonitor.getValue()) * p) / norm; + double corrErr = previousMonitor.getValue().getValue("CRerr").doubleValue() / getCR(previousMonitor.getValue()) / Math.sqrt(2); + return new Pair<>(corrFactor, corrErr); + } + @Override protected void afterAction(String name, Table res, Laminate meta) { printMonitorData(meta); diff --git a/numass-main/src/main/java/inr/numass/actions/SubstractSpectrumAction.java b/numass-main/src/main/java/inr/numass/actions/SubstractSpectrumAction.java new file mode 100644 index 00000000..e761d965 --- /dev/null +++ b/numass-main/src/main/java/inr/numass/actions/SubstractSpectrumAction.java @@ -0,0 +1,61 @@ +/* + * 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.description.TypedActionDef; +import hep.dataforge.io.ColumnedDataReader; +import hep.dataforge.io.ColumnedDataWriter; +import hep.dataforge.io.reports.Reportable; +import hep.dataforge.meta.Laminate; +import hep.dataforge.tables.DataPoint; +import hep.dataforge.tables.ListTable; +import hep.dataforge.tables.MapPoint; +import hep.dataforge.tables.Table; +import java.io.File; +import java.io.IOException; +import java.io.OutputStream; +import java.util.Optional; + +/** + * + * @author Alexander Nozik + */ +@TypedActionDef(name = "substractSpectrum", inputType = Table.class, outputType = Table.class, info = "Substract reference spectrum (background)") +public class SubstractSpectrumAction extends OneToOneAction { + + @Override + protected Table execute(Reportable log, String name, Laminate inputMeta, Table input) { + try { + String referencePath = inputMeta.getString("file", "empty.dat"); + File referenceFile = getContext().io().getFile(referencePath); + Table referenceTable = new ColumnedDataReader(referenceFile).toDataSet(); + ListTable.Builder builder = new ListTable.Builder(input.getFormat()); + input.stream().forEach(point -> { + MapPoint.Builder pointBuilder = new MapPoint.Builder(point); + Optional referencePoint = referenceTable.stream() + .filter(p -> { + return Math.abs(p.getDouble("Uset") - point.getDouble("Uset")) < 0.1; + }).findFirst(); + if (referencePoint.isPresent()) { + pointBuilder.putValue("CR", Math.max(0, point.getDouble("CR") - referencePoint.get().getDouble("CR"))); + pointBuilder.putValue("CRerr", Math.sqrt(Math.pow(point.getDouble("CRerr"), 2d) + Math.pow(referencePoint.get().getDouble("CRerr"), 2d))); + } else { + log.report("No reference point found for Uset = {}", point.getDouble("Uset")); + } + builder.row(pointBuilder.build()); + }); + + Table res = builder.build(); + OutputStream stream = buildActionOutput(name); + ColumnedDataWriter.writeDataSet(stream, res, inputMeta.toString()); + return res; + } catch (IOException ex) { + throw new RuntimeException("Could not read reference file", ex); + } + } + +}