add epsrel argument to likelihood function

3 files changed, 16 insertions, 6 deletions

diff --git a/src/fit.c b/src/fit.c
index 36de76c..87bf232 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -22,8 +22,14 @@
 
 static size_t iter;
 
+typedef struct fitParams {
+    event *ev;
+    double epsrel;
+} fitParams;
+
 double nll(unsigned int n, const double *x, double *grad, void *params)
 {
+    fitParams *fpars = (fitParams *) params;
     double T, theta, phi, t0;
     double pos[3], dir[3];
     double fval;
@@ -48,7 +54,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
     z2[0] = x[8];
 
     gettimeofday(&tv_start, NULL);
-    fval = nll_muon((event *) params, T, pos, dir, t0, z1, z2, 1);
+    fval = nll_muon(fpars->ev, T, pos, dir, t0, z1, z2, 1, fpars->epsrel);
     gettimeofday(&tv_stop, NULL);
 
     long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
@@ -72,11 +78,12 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
 
 int fit_event(event *ev, double *xopt, double *fmin)
 {
+    fitParams fpars;
     double x[9], ss[9], lb[9], ub[9], fval, n;
     int rv;
 
     nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
-    nlopt_set_min_objective(opt,nll,ev);
+    nlopt_set_min_objective(opt,nll,&fpars);
     nlopt_set_ftol_abs(opt, 1e-5);
 
     x[0] = 1200.0;
@@ -126,6 +133,9 @@ int fit_event(event *ev, double *xopt, double *fmin)
 
     nlopt_set_initial_step(opt, ss);
 
+    fpars.ev = ev;
+    fpars.epsrel = 1e-1;
+
     iter = 0;
     rv = nlopt_optimize(opt,x,&fval);
 
diff --git a/src/likelihood.c b/src/likelihood.c
index c90914b..66e72e4 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -122,7 +122,7 @@ double getKineticEnergy(double x, double T0)
     return get_T(T0, x, HEAVY_WATER_DENSITY);
 }
 
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n)
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel)
 {
     size_t i, j;
     intParams params;
@@ -202,14 +202,14 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
         x = R*sin(theta_cerenkov-theta)/sin(theta_cerenkov);
 
         F.function = &gsl_muon_charge;
-        gsl_integration_cquad(&F, 0, range, 0, 1e-2, w, &result, &error, &nevals);
+        gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
         mu_direct[i] = result;
 
         total_charge += mu_direct[i];
 
         if (mu_direct[i] > 0.001) {
             F.function = &gsl_muon_time;
-            gsl_integration_cquad(&F, 0, range, 0, 1e-2, w, &result, &error, &nevals);
+            gsl_integration_cquad(&F, 0, range, 0, epsrel, w, &result, &error, &nevals);
             ts[i] = result;
 
             ts[i] /= mu_direct[i];
diff --git a/src/likelihood.h b/src/likelihood.h
index 9267ef8..4871f42 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -12,6 +12,6 @@
 /* Event window (ns) */
 #define GTVALID 400.0
 
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n);
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel);
 
 #endif