4 files changed, 40 insertions, 3 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 7c3d10f..f426271 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -704,7 +704,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
 
     logp_path = 0.0;
     for (i = 0; i < n; i++)
-        logp_path += log(norm(z1[i],0,1)) + log(norm(z2[i],0,1));
+        logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1);
 
-    return kahan_sum(nll,nhit) + logp_path;
+    return kahan_sum(nll,nhit) - logp_path;
 }
diff --git a/src/misc.c b/src/misc.c
index 6f50033..ba0a090 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -384,6 +384,13 @@ double logsumexp(double *a, size_t n)
     return amax + sum;
 }
 
+double log_norm(double x, double mu, double sigma)
+{
+    /* Returns the log of the PDF for a gaussian random variable with mean `mu`
+     * and standard deviation `sigma`. */
+    return -pow(x-mu,2)/(2*pow(sigma,2)) - log(sqrt(2*M_PI)*sigma);
+}
+
 double norm(double x, double mu, double sigma)
 {
     /* Returns the PDF for a gaussian random variable with mean `mu` and
diff --git a/src/misc.h b/src/misc.h
index 46931cb..dadbbce 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -14,6 +14,7 @@ double interp1d(double x, double *xp, double *yp, size_t n);
 int isclose(double a, double b, double rel_tol, double abs_tol);
 int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol);
 double logsumexp(double *a, size_t n);
+double log_norm(double x, double mu, double sigma);
 double norm(double x, double mu, double sigma);
 double norm_cdf(double x, double mu, double sigma);
 double mean(const double *x, size_t n);
diff --git a/src/test.c b/src/test.c
index 2645982..20a66ec 100644
--- a/src/test.c
+++ b/src/test.c
@@ -823,7 +823,7 @@ err:
 
 int test_std(char *err)
 {
-    /* Tests the mean() function. */
+    /* Tests the std() function. */
     size_t i, j;
     double x[100];
     double sigma, expected, mu;
@@ -851,6 +851,34 @@ err:
     return 1;
 }
 
+int test_log_norm(char *err)
+{
+    /* Tests the log_norm() function. */
+    size_t i;
+    double x, mu, sigma, logp, expected;
+
+    init_genrand(0);
+
+    for (i = 0; i < 100; i++) {
+        x = randn();
+        mu = randn();
+        sigma = genrand_real2() + 1.0;
+
+        logp = log_norm(x,mu,sigma);
+        expected = log(norm(x,mu,sigma));
+
+        if (!isclose(logp, expected, 0, 1e-9)) {
+            sprintf(err, "log_norm(%.2g,%.2g,%.2g) returned %.5g, but expected %.5g", x, mu, sigma, logp, expected);
+            goto err;
+        }
+    }
+
+    return 0;
+
+err:
+    return 1;
+}
+
 struct tests {
     testFunction *test;
     char *name;
@@ -877,6 +905,7 @@ struct tests {
     {test_proton_get_dEdx, "test_proton_get_dEdx"},
     {test_proton_get_range, "test_proton_get_range"},
     {test_proton_get_energy, "test_proton_get_energy"},
+    {test_log_norm, "test_log_norm"},
 };
 
 int main(int argc, char **argv)