diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 36 |
1 files changed, 18 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]; |