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);
+ }
+ }
+
+}