diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 35 |
1 files changed, 15 insertions, 20 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index a57346b..2340fb1 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -303,7 +303,7 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl 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) +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) { /* Returns the approximate expected number of photons seen by PMT `i` using * an analytic formula. @@ -334,7 +334,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, * * `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, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, 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; + 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; @@ -353,13 +353,6 @@ double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, cos_theta = DOT(dir,pmt_dir); sin_theta = sqrt(1-pow(cos_theta,2)); - /* Compute the Cerenkov angle at the start of the track. */ - wavelength0 = 400.0; - n_d2o = get_index_snoman_d2o(wavelength0); - n_h2o = get_index_snoman_h2o(wavelength0); - cos_theta_cerenkov = 1/(n_d2o*beta0); - sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2)); - /* Now, we compute the distance along the track where the PMT is at the * Cerenkov angle. * @@ -441,6 +434,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, 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_acrylic = get_weighted_absorption_length_snoman_acrylic(); @@ -452,7 +446,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, /* Assume the particle is travelling at the speed of light. */ *t = s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; - double charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI); + charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI); *mu_reflected = f_reflected*charge; @@ -534,7 +528,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t size_t i, j, nhit; 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; + double pmt_dir[3], R, cos_theta, sin_theta, wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov, s, l_h2o, l_d2o; double logp_path; double a, b; double tmp[3]; @@ -569,6 +563,13 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t else smax = 0.0; + /* Compute the Cerenkov angle at the start of the track. */ + wavelength0 = 400.0; + n_d2o = get_index_snoman_d2o(wavelength0); + n_h2o = get_index_snoman_h2o(wavelength0); + cos_theta_cerenkov = 1/(n_d2o*beta0); + sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2)); + mu_indirect = 0.0; for (i = 0; i < MAX_PMTS; i++) { if (pmts[i].pmt_type != PMT_NORMAL) continue; @@ -578,7 +579,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); + 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); ts[i] += t0; mu_indirect += mu_reflected; } else { @@ -607,13 +608,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t /* Calculate the cosine of the angle between the track direction and the * vector to the PMT at the start of the track. */ cos_theta = DOT(dir,pmt_dir); - /* Compute the angle between the track direction and the PMT. */ - theta = acos(cos_theta); - - /* Compute the Cerenkov angle at the start of the track. */ - wavelength0 = 400.0; - n_d2o = get_index_snoman_d2o(wavelength0); - theta_cerenkov = acos(1/(n_d2o*beta0)); + sin_theta = sqrt(1-pow(cos_theta,2)); /* Now, we compute the distance along the track where the PMT is at the * Cerenkov angle. @@ -621,7 +616,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t * Note: This formula comes from using the "Law of sines" where the three * vertices of the triangle are the starting position of the track, the * point along the track that we want to find, and the PMT position. */ - s = R*sin(theta_cerenkov-theta)/sin(theta_cerenkov); + s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; /* Make sure that the point is somewhere along the track between 0 and * `smax`. */ |