[no commit message]

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,"")

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

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<Table, Table> {
             res.add(curPoint);
         });
         
-        return new ListTable(res);
+        return new ListTable(TableFormat.forNames(parnames),res);
 
     }
 

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<Table, Table> {
             MapPoint.Builder pb = new MapPoint.Builder(dp);
             pb.putValue("Monitor", 1.0);
             if (!isMonitorPoint(monitor, dp) || index.isEmpty()) {
-                Instant time = getTime(dp);
+                Pair<Double, Double> corr;
-                Entry<Instant, DataPoint> previousMonitor = index.floorEntry(time);
+                if (meta.getBoolean("spline", false)) {
-                Entry<Instant, DataPoint> nextMonitor = index.ceilingEntry(time);
+                    corr = getSplineCorrection(index, dp, norm);
-
-                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;
+                    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<Table, Table> {
         return data;
     }
 
+    private Pair<Double, Double> getSplineCorrection(TreeMap<Instant, DataPoint> 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<Instant, DataPoint> 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<Double, Double> getLinearCorrection(TreeMap<Instant, DataPoint> index, DataPoint dp, double norm) {
+        Instant time = getTime(dp);
+        Entry<Instant, DataPoint> previousMonitor = index.floorEntry(time);
+        Entry<Instant, DataPoint> 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);

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 <altavir@gmail.com>
+ */
+@TypedActionDef(name = "substractSpectrum", inputType = Table.class, outputType = Table.class, info = "Substract reference spectrum (background)")
+public class SubstractSpectrumAction extends OneToOneAction<Table, Table> {
+
+    @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<DataPoint> 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);
+        }
+    }
+
+}