diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/fit.c | 8 | ||||
| -rw-r--r-- | src/likelihood.c | 132 | ||||
| -rw-r--r-- | src/likelihood.h | 2 | ||||
| -rw-r--r-- | src/path.c | 51 | ||||
| -rw-r--r-- | src/path.h | 3 | ||||
| -rw-r--r-- | src/test-path.c | 2 | ||||
| -rw-r--r-- | src/test.c | 2 |
7 files changed, 125 insertions, 75 deletions
@@ -38,7 +38,7 @@ static size_t iter; typedef struct fitParams { event *ev; - int n; + double dx; int fast; int print; int id; @@ -5003,7 +5003,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params) z2[0] = x[8]; gettimeofday(&tv_start, NULL); - fval = nll_muon(fpars->ev, fpars->id, T, pos, dir, t0, z1, z2, 1, fpars->n, fpars->fast); + fval = nll_muon(fpars->ev, fpars->id, T, pos, dir, t0, z1, z2, 1, fpars->dx, fpars->fast); 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; @@ -5211,7 +5211,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id) * faster, we set the absolute tolerance on the likelihood to 1.0, the * maximum number of function evaluations to 100, and the relative * tolerance on the numerical integration to 10%. */ - fpars.n = 100; + fpars.dx = 1.0; fpars.fast = 1; fpars.print = 0; nlopt_set_ftol_abs(opt, 1.0); @@ -5304,7 +5304,7 @@ int fit_event(event *ev, double *xopt, double *fmin, int id) nlopt_set_upper_bounds(opt, ub); /* Now, we do the "real" minimization. */ - fpars.n = 500; + fpars.dx = 1.0; fpars.fast = 0; fpars.print = 1; nlopt_set_ftol_abs(opt, 1e-5); diff --git a/src/likelihood.c b/src/likelihood.c index 707ff38..8116261 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -320,7 +320,7 @@ 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,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(double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad, double pos_a, double pos_b) +static void integrate_path_shower(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) { size_t i; static double qs[MAX_NPOINTS]; @@ -357,10 +357,11 @@ static void integrate_path_shower(double *x, double *pdf, double T0, double *pos ts2[i] = t*t*qs[i]; } - *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)/(*mu_direct); - double t2 = trapz(ts2,(b-a)/(n-1),n)/(*mu_direct); + double dx = x[1]-x[0]; + *mu_direct = trapz(qs,dx,n); + *mu_indirect = trapz(q_indirects,dx,n); + *time = trapz(ts,dx,n)/(*mu_direct); + double t2 = trapz(ts2,dx,n)/(*mu_direct); *sigma = sqrt(t2-(*time)*(*time)); @@ -372,7 +373,7 @@ static void integrate_path_shower(double *x, double *pdf, double T0, double *pos } } -static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time) +static void integrate_path(path *p, int pmt, double *mu_direct, double *mu_indirect, double *time) { size_t i; static double qs[MAX_NPOINTS]; @@ -380,7 +381,7 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl static double ts[MAX_NPOINTS]; double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x; - if (n > MAX_NPOINTS) { + if (p->len > MAX_NPOINTS) { fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n"); exit(1); } @@ -390,10 +391,27 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); - for (i = 0; i < n; i++) { - x = a + i*(b-a)/(n-1); + for (i = 0; i < p->len; i++) { + x = p->s[i]; + + pos[0] = p->x[i]; + pos[1] = p->y[i]; + pos[2] = p->z[i]; - path_eval(p, x, pos, dir, &beta, &t, &theta0); + dir[0] = p->dx[i]; + dir[1] = p->dy[i]; + dir[2] = p->dz[i]; + + beta = p->beta[i]; + + t = p->t[i]; + + if (p->n > 0) { + theta0 = p->theta0; + } else { + theta0 = fmax(MIN_THETA0,p->theta0*sqrt(x)); + theta0 = fmin(MAX_THETA0,theta0); + } get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); @@ -403,9 +421,10 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl ts[i] = t*qs[i]; } - *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); + double dx = p->s[1]-p->s[0]; + *mu_direct = trapz(qs,dx,p->len); + *mu_indirect = trapz(q_indirects,dx,p->len); + *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) @@ -628,8 +647,19 @@ static double getKineticEnergy(double x, void *p) return particle_get_energy(x, (particle *) p); } -double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, size_t npoints, int fast) +double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double dx, int fast) { + /* Returns the negative log likelihood for event `ev` given a particle with + * id `id`, initial kinetic energy `T0`, position `pos`, direction `dir` and + * initial time `t0`. + * + * `dx` is the distance the track is sampled at to compute the likelihood. + * + * If `fast` is nonzero, returns an approximate faster version of the likelihood. + * + * `z1` and `z2` should be arrays of length `n` and represent coefficients + * in the Karhunen Loeve expansion of the track path. These are only used + * when `fast` is zero. */ size_t i, j, nhit; static double logp[MAX_PE], nll[MAX_PMTS]; double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; @@ -649,6 +679,18 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t static double mu[MAX_PMTS]; double mu_noise, mu_indirect; + double *x, *pdf; + + double result; + + double min_ratio; + + double pos_a, pos_b; + + double xlo, xhi; + + size_t npoints, npoints_shower; + /* Initialize static arrays. */ for (i = 0; i < MAX_PMTS; i++) mu_direct[i] = 0.0; for (i = 0; i < MAX_PMTS; i++) mu_shower[i] = 0.0; @@ -657,28 +699,12 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t for (i = 0; i < MAX_PMTS; i++) ts_sigma[i] = PMT_TTS; for (i = 0; i < MAX_PMTS; i++) mu[i] = 0.0; - double pos_a, pos_b; - - electron_get_position_distribution_parameters(T0, &pos_a, &pos_b); - - /* Upper limit to integrate for the electromagnetic shower. */ - double xlo = 0.0; - double xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); - - double *x = malloc(sizeof(double)*npoints); - double *pdf = malloc(sizeof(double)*npoints); - for (i = 0; i < npoints; i++) { - x[i] = xlo + i*(xhi-xlo)/(npoints-1); - pdf[i] = gamma_pdf(x[i],pos_a,pos_b); - } - - double result; - - double min_ratio; - p = particle_init(id, T0, 10000); range = p->range; + /* Number of points to sample along the particle track. */ + npoints = (size_t) (range/dx + 0.5); + /* Calculate total energy */ E0 = T0 + p->mass; p0 = sqrt(E0*E0 - p->mass*p->mass); @@ -687,8 +713,29 @@ 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); - if (!fast) - path = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass); + if (!fast) { + path = path_init(pos, dir, T0, range, npoints, theta0, getKineticEnergy, p, z1, z2, n, p->mass); + + if (id == IDP_E_MINUS || id == IDP_E_PLUS) { + electron_get_position_distribution_parameters(T0, &pos_a, &pos_b); + + /* Lower and upper limits to integrate for the electromagnetic + * shower. */ + xlo = 0.0; + xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); + + /* Number of points to sample along the longitudinal direction for + * the electromagnetic shower. */ + npoints_shower = (size_t) ((xhi-xlo)/dx + 0.5); + + x = malloc(sizeof(double)*npoints_shower); + pdf = malloc(sizeof(double)*npoints_shower); + for (i = 0; i < npoints_shower; i++) { + x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1); + pdf[i] = gamma_pdf(x[i],pos_a,pos_b); + } + } + } if (beta0 > BETA_MIN) get_smax(p, BETA_MIN, range, &smax); @@ -765,7 +812,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t if (b > range) b = range; double q_indirect; - integrate_path(path, i, a, b, npoints, mu_direct+i, &q_indirect, &result); + integrate_path(path, i, mu_direct+i, &q_indirect, &result); mu_indirect += q_indirect; if (mu_direct[i] > 1e-9) { @@ -792,15 +839,21 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t } if (id == IDP_E_MINUS || id == IDP_E_PLUS) { - integrate_path_shower(x,pdf,T0,pos,dir,i,0,xhi,npoints,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad,pos_a,pos_b); + integrate_path_shower(x,pdf,T0,pos,dir,i,npoints_shower,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad); mu_indirect += q_indirect; ts_shower[i] = t0 + result; } } - if (!fast) + if (!fast) { path_free(path); + if (id == IDP_E_MINUS || id == IDP_E_PLUS) { + free(x); + free(pdf); + } + } + particle_free(p); mu_noise = DARK_RATE*GTVALID*1e-9; @@ -863,8 +916,5 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t for (i = 0; i < n; i++) logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1); - free(x); - free(pdf); - return kahan_sum(nll,nhit) - logp_path; } diff --git a/src/likelihood.h b/src/likelihood.h index 7622ffe..f016b53 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -56,6 +56,6 @@ typedef struct particle { particle *particle_init(int id, double T0, size_t n); double particle_get_energy(double x, particle *p); void particle_free(particle *p); -double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, size_t npoints, int fast); +double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double dx, int fast); #endif @@ -10,8 +10,6 @@ #include <gsl/gsl_sf_psi.h> /* for gsl_sf_psi_n() */ #include "scattering.h" -#define N 1000 - static double foo(double s, double range, double theta0, int k) { return sqrt(2*pow(theta0,2)*range)*sin(M_PI*s*(k-0.5)/range)/(M_PI*(k-0.5)); @@ -32,7 +30,7 @@ double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, return sum*pow(k-0.5,2)*pow(M_PI,2)/pow(theta0*range,2); } -path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m) +path *path_init(double *pos, double *dir, double T0, double range, size_t len, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m) { size_t i, j; double T, E, mom, theta, phi; @@ -49,22 +47,23 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 p->dir[2] = dir[2]; p->n = n; - - double *s = malloc(sizeof(double)*N); - double *theta1 = calloc(N,sizeof(double)); - double *theta2 = calloc(N,sizeof(double)); - double *x = calloc(N,sizeof(double)); - double *y = calloc(N,sizeof(double)); - double *z = calloc(N,sizeof(double)); - double *beta = calloc(N,sizeof(double)); - double *t = calloc(N,sizeof(double)); - double *dx = calloc(N,sizeof(double)); - double *dy = calloc(N,sizeof(double)); - double *dz = calloc(N,sizeof(double)); + p->len = len; + + double *s = calloc(len,sizeof(double)); + double *theta1 = calloc(len,sizeof(double)); + double *theta2 = calloc(len,sizeof(double)); + double *x = calloc(len,sizeof(double)); + double *y = calloc(len,sizeof(double)); + double *z = calloc(len,sizeof(double)); + double *beta = calloc(len,sizeof(double)); + double *t = calloc(len,sizeof(double)); + double *dx = calloc(len,sizeof(double)); + double *dy = calloc(len,sizeof(double)); + double *dz = calloc(len,sizeof(double)); dz[0] = 1.0; - for (i = 0; i < N; i++) { - s[i] = range*i/(N-1); + for (i = 0; i < len; i++) { + s[i] = range*i/(len-1); for (j = 0; j < n; j++) { theta1[i] += z1[j]*foo(s[i],range,theta0,j+1); theta2[i] += z2[j]*foo(s[i],range,theta0,j+1); @@ -125,7 +124,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 /* Now we rotate and translate all the positions and rotate all the * directions. */ - for (i = 0; i < N; i++) { + for (i = 0; i < len; i++) { tmp[0] = x[i]; tmp[1] = y[i]; tmp[2] = z[i]; @@ -186,21 +185,21 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 void path_eval(path *p, double s, double *pos, double *dir, double *beta, double *t, double *theta0) { - pos[0] = interp1d(s,p->s,p->x,N); - pos[1] = interp1d(s,p->s,p->y,N); - pos[2] = interp1d(s,p->s,p->z,N); + pos[0] = interp1d(s,p->s,p->x,p->len); + pos[1] = interp1d(s,p->s,p->y,p->len); + pos[2] = interp1d(s,p->s,p->z,p->len); - *beta = interp1d(s,p->s,p->beta,N); - *t = interp1d(s,p->s,p->t,N); + *beta = interp1d(s,p->s,p->beta,p->len); + *t = interp1d(s,p->s,p->t,p->len); /* Technically we should be interpolating between the direction vectors * using an algorithm like Slerp (see https://en.wikipedia.org/wiki/Slerp), * but since we expect the direction vectors to be pretty close to each * other, we just linearly interpolate and then normalize it. */ - dir[0] = interp1d(s,p->s,p->dx,N); - dir[1] = interp1d(s,p->s,p->dy,N); - dir[2] = interp1d(s,p->s,p->dz,N); + dir[0] = interp1d(s,p->s,p->dx,p->len); + dir[1] = interp1d(s,p->s,p->dy,p->len); + dir[2] = interp1d(s,p->s,p->dz,p->len); normalize(dir); @@ -22,6 +22,7 @@ typedef struct path { double theta0; size_t n; + size_t len; double *s; double *x; @@ -35,7 +36,7 @@ typedef struct path { } path; double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, size_t n); -path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m); +path *path_init(double *pos, double *dir, double T0, double range, size_t len, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m); void path_eval(path *p, double s, double *pos, double *dir, double *beta, double *t, double *theta0); void path_free(path *p); diff --git a/src/test-path.c b/src/test-path.c index 62994d6..f403c37 100644 --- a/src/test-path.c +++ b/src/test-path.c @@ -101,7 +101,7 @@ sim *simulate_path(double *pos0, double *dir0, double range, double theta0, size simvalue->z2[i] = path_get_coefficient(i+1,s,theta2,theta0,N); } - simvalue->p = path_init(pos0,dir0,1.0,range,theta0,getKineticEnergy,NULL,simvalue->z1,simvalue->z2,n,1.0); + simvalue->p = path_init(pos0,dir0,1.0,range,1000,theta0,getKineticEnergy,NULL,simvalue->z1,simvalue->z2,n,1.0); free(s); free(theta1); @@ -498,7 +498,7 @@ int test_path(char *err) rand_sphere(dir0); - path *p = path_init(pos0,dir0,T0,range,0.1,getKineticEnergy,NULL,NULL,NULL,0,m); + path *p = path_init(pos0,dir0,T0,range,1000,0.1,getKineticEnergy,NULL,NULL,NULL,0,m); for (j = 0; j < 100; j++) { s = range*j/99; |