a bunch of small changes to speed things up

3 files changed, 47 insertions, 18 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index dcf681a..fcc8a5f 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -323,7 +323,10 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
      * different transit times across the face of the PMT, it seems better to
      * convolve first which is what we do here. In addition, the problem is not
      * analytically tractable if you do things the other way around. */
-    return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
+    if (n == 1)
+        return ln(n) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
+    else
+        return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
 }
 
 static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad)
@@ -430,7 +433,7 @@ static void integrate_path(path *p, int pmt, double *mu_direct, double *mu_indir
     *time = trapz(ts,dx,p->len);
 }
 
-static double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
+static double get_total_charge_approx(double beta0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov)
 {
     /* Returns the approximate expected number of photons seen by PMT `i` using
      * an analytic formula.
@@ -461,18 +464,14 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
      *
      * `smax` is currently calculated as the point where the particle velocity
      * drops to 0.8 times the speed of light. */
-    double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, sin_theta, E0, p0, beta0, f, cos_theta_pmt, theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected, charge;
-
-    /* Calculate beta at the start of the track. */
-    E0 = T0 + p->mass;
-    p0 = sqrt(E0*E0 - p->mass*p->mass);
-    beta0 = p0/E0;
+    double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, sin_theta, f, cos_theta_pmt, theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, f_reflected, charge, prob, frac;
 
     /* First, we find the point along the track where the PMT is at the
      * Cerenkov angle. */
     SUB(pmt_dir,pmts[i].pos,pos);
     /* Compute the distance to the PMT. */
     R = NORM(pmt_dir);
+
     normalize(pmt_dir);
 
     /* Calculate the cosine of the angle between the track direction and the
@@ -530,7 +529,7 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
 
     /* `prob` is the number of photons emitted per cm by the particle at a
      * distance `s` along the track. */
-    double prob = get_probability2(beta);
+    prob = get_probability2(beta);
 
     /* Compute the position of the particle at a distance `s` along the track. */
     tmp[0] = pos[0] + s*dir[0];
@@ -539,7 +538,9 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
 
     SUB(pmt_dir,pmts[i].pos,tmp);
 
-    cos_theta_pmt = DOT(pmt_dir,pmts[i].normal)/NORM(pmt_dir);
+    normalize(pmt_dir);
+
+    cos_theta_pmt = DOT(pmt_dir,pmts[i].normal);
 
     *t = 0.0;
     *mu_reflected = 0.0;
@@ -547,7 +548,7 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
 
     /* Calculate the sine of the angle between the track direction and the PMT
      * at the position `s` along the track. */
-    cos_theta = DOT(dir,pmt_dir)/NORM(pmt_dir);
+    cos_theta = DOT(dir,pmt_dir);
     sin_theta = sqrt(1-pow(cos_theta,2));
 
     /* Get the solid angle of the PMT at the position `s` along the track. */
@@ -555,15 +556,14 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti
 
     theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
 
-    double frac = sqrt(2)*n_d2o*x*beta0*theta0;
+    frac = sqrt(2)*n_d2o*x*beta0*theta0;
 
     theta_pmt = acos(-cos_theta_pmt);
     f = get_weighted_pmt_response(theta_pmt);
     f_reflected = get_weighted_pmt_reflectivity(theta_pmt);
 
-    wavelength0 = 400.0;
-    absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0);
-    absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0);
+    absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
+    absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
     absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
 
     get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o);
@@ -834,7 +834,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
         if (fast && !ev->pmt_hits[i].hit) continue;
 
         if (fast) {
-            mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
+            mu_direct[i] = get_total_charge_approx(beta0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov);
             ts[i] += t0;
             mu_indirect += mu_reflected;
             continue;
@@ -926,9 +926,9 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
         if (ev->pmt_hits[i].hit) {
             for (j = 1; j < MAX_PE; j++) {
                 if (fast)
-                    logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt_fast;
+                    logp[j] = log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt_fast;
                 else
-                    logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]);
+                    logp[j] = log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]);
 
                 if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
 
diff --git a/src/sno_charge.c b/src/sno_charge.c
index 1f39706..a019134 100644
--- a/src/sno_charge.c
+++ b/src/sno_charge.c
@@ -14,6 +14,8 @@ double qhi = 100.0;
 size_t nq = 10000;
 static gsl_spline *splines[MAX_PE];
 static gsl_interp_accel *acc;
+static gsl_spline *log_splines[MAX_PE];
+static gsl_interp_accel *log_acc;
 
 static double qmean, qstd;
 
@@ -84,6 +86,23 @@ double spe_pol2exp(double q)
     return (gmratio*(NG1*funpol1 + NG2*funpol2) + NEXP*exp(-q/Q0))/norm;
 }
 
+double log_pq(double q, int n)
+{
+    if (!initialized) {
+        fprintf(stderr, "charge interpolation hasn't been initialized!\n");
+        exit(1);
+    }
+
+    if (n > MAX_PE) {
+        /* Assume the distribution is gaussian by the central limit theorem. */
+        return log(norm(q,n*qmean,sqrt(n)*qstd));
+    }
+
+    if (q < qlo || q > qhi) return -INFINITY;
+
+    return interp1d(q,log_splines[n-1]->x,log_splines[n-1]->y,log_splines[n-1]->size);
+}
+
 double pq(double q, int n)
 {
     if (!initialized) {
@@ -240,6 +259,15 @@ void init_charge(void)
     }
 
     gsl_integration_cquad_workspace_free(w);
+
+    for (i = 1; i <= MAX_PE; i++) {
+        for (j = 0; j < nq; j++) {
+            y[j] = log(splines[i-1]->y[j]);
+        }
+        log_splines[i-1] = gsl_spline_alloc(gsl_interp_linear,nq);
+        gsl_spline_init(log_splines[i-1],x,y,nq);
+    }
+
     free(x);
     free(y);
 }
diff --git a/src/sno_charge.h b/src/sno_charge.h
index 1004e30..50b018e 100644
--- a/src/sno_charge.h
+++ b/src/sno_charge.h
@@ -2,6 +2,7 @@
 #define SNO_CHARGE_H
 
 void init_charge(void);
+double log_pq(double q, int n);
 double pq(double q, int n);
 double get_log_pmiss(int n);
 double spe_pol2exp(double q);