Author: [log in to unmask]
Date: Thu Dec 4 17:24:25 2014
New Revision: 1649
Log:
Add implementation of landau distribution (in my user area for now but will be moved).
Added:
java/trunk/users/src/main/java/org/hps/users/jeremym/LandauFunction.java
java/trunk/users/src/main/java/org/hps/users/jeremym/LandauPdf.java
java/trunk/users/src/test/java/org/hps/users/jeremym/LandauTest.java
Added: java/trunk/users/src/main/java/org/hps/users/jeremym/LandauFunction.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/jeremym/LandauFunction.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/jeremym/LandauFunction.java Thu Dec 4 17:24:25 2014
@@ -0,0 +1,50 @@
+package org.hps.users.jeremym;
+
+import hep.aida.ref.function.AbstractIFunction;
+
+public class LandauFunction extends AbstractIFunction {
+
+ LandauPdf landauPdf = new LandauPdf();
+
+ public LandauFunction() {
+ this("");
+ }
+
+ public LandauFunction(String title) {
+
+ super();
+
+ variableNames = new String[1];
+ variableNames[0] = "x0";
+
+ parameterNames = new String[2];
+ parameterNames[0] = "mean";
+ parameterNames[1] = "sigma";
+
+ init(title);
+ }
+
+ @Override
+ public double value(double[] v) {
+ return landauPdf.getValue(v[0]);
+ }
+
+ @Override
+ public void setParameter(String key, double value) throws IllegalArgumentException {
+ super.setParameter(key, value);
+ if (key.equals("mean")) {
+ System.out.println("set mean = " + value);
+ landauPdf.setMean(value);
+ } else if (key.equals("sigma")) {
+ System.out.println("set sigma = " + value);
+ landauPdf.setSigma(value);
+ }
+ }
+
+ @Override
+ public void setParameters(double[] parameters) throws IllegalArgumentException {
+ super.setParameters(parameters);
+ landauPdf.setMean(parameters[0]);
+ landauPdf.setSigma(parameters[1]);
+ }
+}
Added: java/trunk/users/src/main/java/org/hps/users/jeremym/LandauPdf.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/jeremym/LandauPdf.java (added)
+++ java/trunk/users/src/main/java/org/hps/users/jeremym/LandauPdf.java Thu Dec 4 17:24:25 2014
@@ -0,0 +1,98 @@
+package org.hps.users.jeremym;
+
+
+/**
+ * <p>
+ * Landau Probability Distribution
+ * <p>
+ * Copied from math/mathcore/src/ProbFuncMathCore.cxx::landau_pdf in ROOT (version 5.34.18).
+ */
+public class LandauPdf {
+
+ static final double[] p1 = { 0.4259894875, -0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253 };
+ static final double[] q1 = { 1.0, -0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063 };
+
+ static final double[] p2 = { 0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211 };
+ static final double[] q2 = { 1.0, 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714 };
+
+ static final double[] p3 = { 0.1788544503, 0.09359161662, 0.006325387654, 0.00006611667319, -0.000002031049101 };
+ static final double[] q3 = { 1.0, 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675 };
+
+ static final double[] p4 = { 0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186 };
+ static final double[] q4 = { 1.0, 106.8615961, 337.6496214, 2016.712389, 1597.063511 };
+
+ static final double[] p5 = { 1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910 };
+ static final double[] q5 = { 1.0, 156.9424537, 3745.310488, 9834.698876, 66924.28357 };
+
+ static final double[] p6 = { 1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109 };
+ static final double[] q6 = { 1.0, 651.4101098, 56974.73333, 165917.4725, -2815759.939 };
+
+ static final double[] a1 = { 0.04166666667, -0.01996527778, 0.02709538966 };
+ static final double[] a2 = { -1.845568670, -4.284640743 };
+
+ protected double sigma = 0;
+ protected double mean = 0;
+
+ public LandauPdf(double mean, double sigma) {
+ this.mean = mean;
+ this.sigma = sigma;
+ }
+
+ public LandauPdf() {
+ }
+
+ public double getValue(double x) {
+ return getValue(x, mean, sigma);
+ }
+
+ public void setSigma(double sigma) {
+ this.sigma = sigma;
+ }
+
+ public void setMean(double mean) {
+ this.mean = mean;
+ }
+
+ /**
+ *
+ * @param x The free variable.
+ * @param sigma The sigma (width) of the distribution.
+ * @param mean The mean of the distribution.
+ * @return The value of the Landau distribution at x.
+ */
+ static double getValue(double x, double mean, double sigma) {
+ // LANDAU pdf : algorithm from CERNLIB G110 denlan same algorithm is used in GSL
+ if (sigma <= 0)
+ return 0;
+ double v = (x - mean) / sigma;
+ double u, ue, us, denlan;
+ if (v < -5.5) {
+ u = Math.exp(v + 1.0);
+ if (u < 1e-10)
+ return 0.0;
+ ue = Math.exp(-1 / u);
+ us = Math.sqrt(u);
+ denlan = 0.3989422803 * (ue / us) * (1 + (a1[0] + (a1[1] + a1[2] * u) * u) * u);
+ } else if (v < -1) {
+ u = Math.exp(-v - 1);
+ denlan = Math.exp(-u) * Math.sqrt(u) * (p1[0] + (p1[1] + (p1[2] + (p1[3] + p1[4] * v) * v) * v) * v) / (q1[0] + (q1[1] + (q1[2] + (q1[3] + q1[4] * v) * v) * v) * v);
+ } else if (v < 1) {
+ denlan = (p2[0] + (p2[1] + (p2[2] + (p2[3] + p2[4] * v) * v) * v) * v) / (q2[0] + (q2[1] + (q2[2] + (q2[3] + q2[4] * v) * v) * v) * v);
+ } else if (v < 5) {
+ denlan = (p3[0] + (p3[1] + (p3[2] + (p3[3] + p3[4] * v) * v) * v) * v) / (q3[0] + (q3[1] + (q3[2] + (q3[3] + q3[4] * v) * v) * v) * v);
+ } else if (v < 12) {
+ u = 1 / v;
+ denlan = u * u * (p4[0] + (p4[1] + (p4[2] + (p4[3] + p4[4] * u) * u) * u) * u) / (q4[0] + (q4[1] + (q4[2] + (q4[3] + q4[4] * u) * u) * u) * u);
+ } else if (v < 50) {
+ u = 1 / v;
+ denlan = u * u * (p5[0] + (p5[1] + (p5[2] + (p5[3] + p5[4] * u) * u) * u) * u) / (q5[0] + (q5[1] + (q5[2] + (q5[3] + q5[4] * u) * u) * u) * u);
+ } else if (v < 300) {
+ u = 1 / v;
+ denlan = u * u * (p6[0] + (p6[1] + (p6[2] + (p6[3] + p6[4] * u) * u) * u) * u) / (q6[0] + (q6[1] + (q6[2] + (q6[3] + q6[4] * u) * u) * u) * u);
+ } else {
+ u = 1 / (v - v * Math.log(v) / (v + 1));
+ denlan = u * u * (1 + (a2[0] + a2[1] * u) * u);
+ }
+ return denlan / sigma;
+ }
+}
Added: java/trunk/users/src/test/java/org/hps/users/jeremym/LandauTest.java
=============================================================================
--- java/trunk/users/src/test/java/org/hps/users/jeremym/LandauTest.java (added)
+++ java/trunk/users/src/test/java/org/hps/users/jeremym/LandauTest.java Thu Dec 4 17:24:25 2014
@@ -0,0 +1,117 @@
+package org.hps.users.jeremym;
+
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IPlotter;
+import hep.aida.IProfile1D;
+import hep.aida.ref.fitter.FitResult;
+
+import java.io.IOException;
+import java.util.Random;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import junit.framework.TestCase;
+
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * <p>
+ * Test of Landau distribution functions ported from C++ and originally defined in the file
+ * <p>
+ * mathcore/src/ProbFuncMathCore.cxx
+ * <p>
+ * within ROOT (version 5.34.18).
+ */
+public class LandauTest extends TestCase {
+
+ static double[] p1 = { 0.4259894875, -0.1249762550, 0.03984243700, -0.006298287635, 0.001511162253 };
+ static double[] q1 = { 1.0, -0.3388260629, 0.09594393323, -0.01608042283, 0.003778942063 };
+
+ static double[] p2 = { 0.1788541609, 0.1173957403, 0.01488850518, -0.001394989411, 0.0001283617211 };
+ static double[] q2 = { 1.0, 0.7428795082, 0.3153932961, 0.06694219548, 0.008790609714 };
+
+ static double[] p3 = { 0.1788544503, 0.09359161662, 0.006325387654, 0.00006611667319, -0.000002031049101 };
+ static double[] q3 = { 1.0, 0.6097809921, 0.2560616665, 0.04746722384, 0.006957301675 };
+
+ static double[] p4 = { 0.9874054407, 118.6723273, 849.2794360, -743.7792444, 427.0262186 };
+ static double[] q4 = { 1.0, 106.8615961, 337.6496214, 2016.712389, 1597.063511 };
+
+ static double[] p5 = { 1.003675074, 167.5702434, 4789.711289, 21217.86767, -22324.94910 };
+ static double[] q5 = { 1.0, 156.9424537, 3745.310488, 9834.698876, 66924.28357 };
+
+ static double[] p6 = { 1.000827619, 664.9143136, 62972.92665, 475554.6998, -5743609.109 };
+ static double[] q6 = { 1.0, 651.4101098, 56974.73333, 165917.4725, -2815759.939 };
+
+ static double[] a1 = { 0.04166666667, -0.01996527778, 0.02709538966 };
+ static double[] a2 = { -1.845568670, -4.284640743 };
+
+ public void testLandauDistribution() {
+ AIDA aida = AIDA.defaultInstance();
+ double mean = 10.0;
+ double sigma = 2.0;
+
+ LandauPdf landauPdf = new LandauPdf();
+ landauPdf.setMean(mean);
+ landauPdf.setSigma(sigma);
+
+ Random random = new Random();
+
+ IProfile1D profile = aida.profile1D("Landau P1D", 50, -0.5, 49.5);
+ for (int i = 0; i < 100; i++) {
+ for (int x = 0; x < 2000; x++) {
+ double xVal = x / 10;
+ double landauVal = landauPdf.getValue(xVal);
+ double smearedLandau = landauVal * (1.0 + 0.1 * random.nextGaussian());
+ profile.fill(xVal, smearedLandau);
+ }
+ }
+
+ IPlotter plotter = aida.analysisFactory().createPlotterFactory().create();
+ plotter.createRegion();
+ plotter.region(0).plot(profile);
+
+ LandauFunction landau = new LandauFunction();
+ //IFunctionFactory functionFactory = aida.analysisFactory().createFunctionFactory(null);
+ IFitFactory fitFactory = aida.analysisFactory().createFitFactory();
+ //functionFactory.catalog().add("landau", landau);
+
+ landau.setParameter("sigma", sigma + 0.1);
+ landau.setParameter("mean", mean + 0.1);
+
+ IFitter fitter = fitFactory.createFitter();
+
+ Logger minuitLogger = Logger.getLogger("org.freehep.math.minuit");
+ minuitLogger.setLevel(Level.ALL);
+ minuitLogger.info("minuit logger test");
+
+ IFitResult fitResult = fitter.fit(profile, landau);
+
+ IFunction fittedFunction = fitResult.fittedFunction();
+
+ plotter.region(0).plot(fittedFunction);
+
+ System.out.println();
+ System.out.println("fitted parameters");
+ for (String fittedParameterName : fitResult.fittedParameterNames()) {
+ System.out.println(fittedParameterName + " = " + fitResult.fittedParameter(fittedParameterName));
+ }
+ System.out.println("fit status = " + fitResult.fitStatus());
+
+ System.out.println();
+ System.out.println("fitted function");
+ for (String parameterName : fittedFunction.parameterNames()) {
+ System.out.println(parameterName + " = " + fittedFunction.parameter(parameterName));
+ }
+
+ ((FitResult)fitResult).printResult();
+
+ try {
+ plotter.writeToFile("Landau.png");
+ } catch (IOException e) {
+ throw new RuntimeException(e);
+ }
+ }
+}
