analytical function for resolution tail

This commit is contained in:
darksnake 2017-06-15 17:05:17 +03:00
parent 6cc2d5596e
commit 24e48d2d7f
2 changed files with 36 additions and 21 deletions

View File

@ -40,7 +40,6 @@ import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.util.FastMath; import org.apache.commons.math3.util.FastMath;
/** /**
*
* @author Alexander Nozik * @author Alexander Nozik
*/ */
@PluginDef(group = "inr.numass", name = "numass", @PluginDef(group = "inr.numass", name = "numass",
@ -123,6 +122,18 @@ public class NumassPlugin extends BasicPlugin {
double beta = meta.getDouble("tailBeta", 0); double beta = meta.getDouble("tailBeta", 0);
return (double E, double U) -> 1 - (E - U) * (alpha + E / 1000d * beta) / 1000d; 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));
});
} }
/** /**

View File

@ -20,6 +20,7 @@ import hep.dataforge.maths.integration.GaussRuleIntegrator;
import hep.dataforge.maths.integration.UnivariateIntegrator; import hep.dataforge.maths.integration.UnivariateIntegrator;
import hep.dataforge.plots.PlotFrame; import hep.dataforge.plots.PlotFrame;
import hep.dataforge.plots.data.PlottableXYFunction; import hep.dataforge.plots.data.PlottableXYFunction;
import hep.dataforge.utils.Misc;
import hep.dataforge.values.NamedValueSet; import hep.dataforge.values.NamedValueSet;
import org.apache.commons.math3.analysis.BivariateFunction; import org.apache.commons.math3.analysis.BivariateFunction;
import org.apache.commons.math3.analysis.UnivariateFunction; import org.apache.commons.math3.analysis.UnivariateFunction;
@ -48,8 +49,10 @@ public class LossCalculator {
private static final LossCalculator instance = new LossCalculator(); private static final LossCalculator instance = new LossCalculator();
private static final UnivariateIntegrator integrator = new GaussRuleIntegrator(100); 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 final Map<Integer, UnivariateFunction> cache = new HashMap<>();
private LossCalculator() { private LossCalculator() {
cache.put(1, getSingleScatterFunction()); cache.put(1, getSingleScatterFunction());
// cache.put(2, getDoubleScatterFunction()); // cache.put(2, getDoubleScatterFunction());
@ -281,36 +284,37 @@ public class LossCalculator {
* рекурсивно вычисляем все вероятности, котрорые выше порога * рекурсивно вычисляем все вероятности, котрорые выше порога
* <p> * <p>
* дисер, стр.48 * дисер, стр.48
* * </p>
* @param X * @param X
* @return * @return
*/ */
public List<Double> getLossProbabilities(double X) { 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<>(); while (prob > SCATTERING_PROBABILITY_THRESHOLD) {
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) {
/* /*
* prob(n) = prob(n-1)-1/n! * X^n * exp(-X); * prob(n) = prob(n-1)-1/n! * X^n * exp(-X);
*/ */
double delta = Math.exp(-X); double delta = Math.exp(-x);
for (int i = 1; i < res.size() + 1; i++) { for (int i = 1; i < res.size() + 1; i++) {
delta *= X / 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) { public double getLossProbability(int order, double X) {