diff options
| -rw-r--r-- | src/likelihood.c | 36 | ||||
| -rw-r--r-- | src/sno_charge.c | 28 | ||||
| -rw-r--r-- | src/sno_charge.h | 1 |
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); |