From 24e48d2d7f581e83058d23f60391ae1b49a88de4 Mon Sep 17 00:00:00 2001
From: darksnake <darksnake@LUNA.localdomain>
Date: Thu, 15 Jun 2017 17:05:17 +0300
Subject: [PATCH] analytical function for resolution tail

---
 .../main/java/inr/numass/NumassPlugin.java    | 13 +++++-
 .../inr/numass/models/LossCalculator.java     | 44 ++++++++++---------
 2 files changed, 36 insertions(+), 21 deletions(-)

diff --git a/numass-main/src/main/java/inr/numass/NumassPlugin.java b/numass-main/src/main/java/inr/numass/NumassPlugin.java
index 561ff11a..8dff6cc1 100644
--- a/numass-main/src/main/java/inr/numass/NumassPlugin.java
+++ b/numass-main/src/main/java/inr/numass/NumassPlugin.java
@@ -40,7 +40,6 @@ import org.apache.commons.math3.analysis.UnivariateFunction;
 import org.apache.commons.math3.util.FastMath;
 
 /**
- *
  * @author Alexander Nozik
  */
 @PluginDef(group = "inr.numass", name = "numass",
@@ -123,6 +122,18 @@ public class NumassPlugin extends BasicPlugin {
             double beta = meta.getDouble("tailBeta", 0);
             return (double E, double U) -> 1 - (E - U) * (alpha + E / 1000d * beta) / 1000d;
         });
+
+        math.registerBivariate("numass.resolutionTail.2017", meta ->
+                (double E, double U) -> {
+                    double D = E - U;
+                    return 0.99797 - 3.05346E-7 * D - 5.45738E-10 * Math.pow(D, 2) - 6.36105E-14 * Math.pow(D, 3);
+                });
+
+        math.registerBivariate("numass.resolutionTail.2017.mod", meta ->
+                (double E, double U) -> {
+                    double D = E - U;
+                    return (0.99797 - 3.05346E-7 * D - 5.45738E-10 * Math.pow(D, 2) - 6.36105E-14 * Math.pow(D, 3));
+                });
     }
 
     /**
diff --git a/numass-main/src/main/java/inr/numass/models/LossCalculator.java b/numass-main/src/main/java/inr/numass/models/LossCalculator.java
index fe25de38..45ba0ac6 100644
--- a/numass-main/src/main/java/inr/numass/models/LossCalculator.java
+++ b/numass-main/src/main/java/inr/numass/models/LossCalculator.java
@@ -20,6 +20,7 @@ import hep.dataforge.maths.integration.GaussRuleIntegrator;
 import hep.dataforge.maths.integration.UnivariateIntegrator;
 import hep.dataforge.plots.PlotFrame;
 import hep.dataforge.plots.data.PlottableXYFunction;
+import hep.dataforge.utils.Misc;
 import hep.dataforge.values.NamedValueSet;
 import org.apache.commons.math3.analysis.BivariateFunction;
 import org.apache.commons.math3.analysis.UnivariateFunction;
@@ -48,8 +49,10 @@ public class LossCalculator {
 
     private static final LossCalculator instance = new LossCalculator();
     private static final UnivariateIntegrator integrator = new GaussRuleIntegrator(100);
+    private static final Map<Double, List<Double>> lossProbCache = Misc.getLRUCache(100);
     private final Map<Integer, UnivariateFunction> cache = new HashMap<>();
 
+
     private LossCalculator() {
         cache.put(1, getSingleScatterFunction());
 //        cache.put(2, getDoubleScatterFunction());
@@ -281,36 +284,37 @@ public class LossCalculator {
      * рекурсивно вычисляем все вероятности, котрорые выше порога
      * <p>
      * дисер, стр.48
-     *
+     * </p>
      * @param X
      * @return
      */
     public List<Double> getLossProbabilities(double X) {
+        return lossProbCache.computeIfAbsent(X, (x) -> {
+            List<Double> res = new ArrayList<>();
+            double prob;
+            if (x > 0) {
+                prob = 1 / x * (1 - Math.exp(-x));
+            } else {
+                // если x ==0, то выживает только нулевой член, первый равен нулю
+                res.add(1d);
+                return res;
+            }
+            res.add(prob);
 
-        List<Double> res = new ArrayList<>();
-        double prob;
-        if (X > 0) {
-            prob = 1 / X * (1 - Math.exp(-X));
-        } else {
-            // если x ==0, то выживает только нулевой член, первый равен нулю
-            res.add(1d);
-            return res;
-        }
-        res.add(prob);
-
-        while (prob > SCATTERING_PROBABILITY_THRESHOLD) {
+            while (prob > SCATTERING_PROBABILITY_THRESHOLD) {
             /*
             * prob(n) = prob(n-1)-1/n! * X^n * exp(-X);
              */
-            double delta = Math.exp(-X);
-            for (int i = 1; i < res.size() + 1; i++) {
-                delta *= X / i;
+                double delta = Math.exp(-x);
+                for (int i = 1; i < res.size() + 1; i++) {
+                    delta *= x / i;
+                }
+                prob -= delta / x;
+                res.add(prob);
             }
-            prob -= delta / X;
-            res.add(prob);
-        }
 
-        return res;
+            return res;
+        });
     }
 
     public double getLossProbability(int order, double X) {