From 491c6fc921dbb42b83986db7c8577030827a856b Mon Sep 17 00:00:00 2001 From: tlatorre Date: Fri, 19 Oct 2018 10:53:42 -0500 Subject: update path integral to use a fixed number of points I noticed when fitting electrons that the cquad integration routine was not very stable, i.e. it would return different results for *very* small changes in the fit parameters which would cause the fit to stall. Since it's very important for the minimizer that the likelihood function not jump around, I am switching to integrating over the path by just using a fixed number of points and using the trapezoidal rule. This seems to be a lot more stable, and as a bonus I was able to combine the three integrals (direct charge, indirect charge, and time) so that we only have to do a single loop. This should hopefully make the speed comparable since the cquad routine was fairly effective at only using as many function evaluations as needed. Another benefit to this approach is that if needed, it will be easier to port to a GPU. --- src/likelihood.c | 108 +++++++++++++++++++++---------------------------------- 1 file changed, 40 insertions(+), 68 deletions(-) (limited to 'src/likelihood.c') diff --git a/src/likelihood.c b/src/likelihood.c index 1aa1653..ec58d11 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -22,11 +22,6 @@ #include "proton.h" #include "id_particles.h" -typedef struct intParams { - path *p; - int i; -} intParams; - particle *particle_init(int id, double T0, size_t n) { /* Returns a struct with arrays of the particle position and kinetic @@ -160,9 +155,9 @@ void particle_free(particle *p) free(p); } -double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected) +static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o) { - double pmt_dir[3], cos_theta, n, wavelength0, omega, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic; + double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; z = 1.0; @@ -182,31 +177,29 @@ double get_expected_charge(double x, double beta, double theta0, double *pos, do R = NORM(pos); - /* FIXME: I just calculate delta assuming 400 nm light. */ - wavelength0 = 400.0; - if (R <= AV_RADIUS) { - n = get_index_snoman_d2o(wavelength0); + n = n_d2o; } else { - n = get_index_snoman_h2o(wavelength0); + n = n_h2o; } if (beta < 1/n) return 0; - if (reflected) - f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt)); - else - f = get_weighted_pmt_response(acos(-cos_theta_pmt)); + theta_pmt = acos(-cos_theta_pmt); + f_reflec = get_weighted_pmt_reflectivity(theta_pmt); + f = get_weighted_pmt_response(theta_pmt); 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(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o); - l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; - return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); + charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); + + *reflected = f_reflec*charge; + + return f*charge; } double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma) @@ -271,43 +264,41 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)); } -static double gsl_time(double x, void *params) +static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time) { - intParams *pars = (intParams *) params; - double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o; - - path_eval(pars->p, x, pos, dir, &beta, &t, &theta0); + size_t i; + double *qs, *q_indirects, *ts; + double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x; - get_path_length(pos,pmts[pars->i].pos,AV_RADIUS,&l_d2o,&l_h2o); + qs = malloc(sizeof(double)*n); + q_indirects = malloc(sizeof(double)*n); + ts = malloc(sizeof(double)*n); /* FIXME: I just calculate delta assuming 400 nm light. */ wavelength0 = 400.0; n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); - t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; + for (i = 0; i < n; i++) { + x = a + i*(b-a)/(n-1); - return t*get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0); -} + path_eval(p, x, pos, dir, &beta, &t, &theta0); -static double gsl_reflected_charge(double x, void *params) -{ - intParams *pars = (intParams *) params; - double dir[3], pos[3], t, theta0, beta; + get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); - path_eval(pars->p, x, pos, dir, &beta, &t, &theta0); + t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; - return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 1); -} - -static double gsl_charge(double x, void *params) -{ - intParams *pars = (intParams *) params; - double dir[3], pos[3], t, theta0, beta; + qs[i] = get_expected_charge(x, beta, theta0, pos, dir, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o); + ts[i] = t*qs[i]; + } - path_eval(pars->p, x, pos, dir, &beta, &t, &theta0); + *mu_direct = trapz(qs,(b-a)/(n-1),n); + *mu_indirect = trapz(q_indirects,(b-a)/(n-1),n); + *time = trapz(ts,(b-a)/(n-1),n); - return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0); + free(qs); + free(q_indirects); + free(ts); } double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected) @@ -532,7 +523,6 @@ static double getKineticEnergy(double x, void *p) double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast) { size_t i, j, nhit; - intParams params; double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp; particle *p; double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o; @@ -540,19 +530,14 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t double a, b; double tmp[3]; double mu_reflected; + path *path; double mu_direct[MAX_PMTS]; double ts[MAX_PMTS]; double mu[MAX_PMTS]; double mu_noise, mu_indirect; - gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(100); - double result, error; - - size_t nevals; - - gsl_function F; - F.params = ¶ms; + double result; p = particle_init(id, T0, 10000); range = p->range; @@ -565,7 +550,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t /* FIXME: is this formula valid for muons? */ theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2); - params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass); + path = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass); if (beta0 > BETA_MIN) get_smax(p, BETA_MIN, range, &smax); @@ -580,8 +565,6 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t * calculation. */ if (fast && !ev->pmt_hits[i].hit) continue; - params.i = i; - if (fast) { mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected); ts[i] += t0; @@ -640,20 +623,11 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t if (a < 0.0) a = 0.0; if (b > range) b = range; - F.function = &gsl_charge; - - gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); - mu_direct[i] = result; - - F.function = &gsl_reflected_charge; - - gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); - mu_indirect += result; - - F.function = &gsl_time; + double q_indirect; + integrate_path(path, i, a, b, 100, mu_direct+i, &q_indirect, &result); + mu_indirect += q_indirect; if (mu_direct[i] > 1e-9) { - gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); ts[i] = t0 + result/mu_direct[i]; } else { /* If the expected number of PE is very small, our estimate of @@ -678,12 +652,10 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t } } - path_free(params.p); + path_free(path); particle_free(p); - gsl_integration_cquad_workspace_free(w); - mu_noise = DARK_RATE*GTVALID*1e-9; mu_indirect *= CHARGE_FRACTION/10000.0; -- cgit