diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 113 |
1 files changed, 50 insertions, 63 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 8116261..4e44b00 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -642,6 +642,54 @@ static int get_smax(particle *p, double beta_min, double range, double *smax) return status; } +static double guess_time(double *pos, double *dir, double *pmt_pos, double smax, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov) +{ + /* Returns an approximate time at which a PMT is most likely to get hit + * from Cerenkov light. + * + * To do this, we try to compute the distance along the track where the PMT + * is at the Cerenkov angle. */ + double pmt_dir[3], tmp[3]; + double R, cos_theta, sin_theta, s, l_d2o, l_h2o; + + /* First, we find the point along the track where the PMT is at the + * Cerenkov angle. */ + SUB(pmt_dir,pmt_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 + * vector to the PMT at the start of the track. */ + cos_theta = DOT(dir,pmt_dir); + sin_theta = sqrt(1-pow(cos_theta,2)); + + /* Now, we compute the distance along the track where the PMT is at the + * Cerenkov angle. + * + * 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*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; + + /* Make sure that the point is somewhere along the track between 0 and + * `smax`. */ + if (s < 0) s = 0.0; + else if (s > smax) s = smax; + + /* Compute the position of the particle at a distance `s` along the track. */ + tmp[0] = pos[0] + s*dir[0]; + tmp[1] = pos[1] + s*dir[1]; + tmp[2] = pos[2] + s*dir[2]; + + SUB(pmt_dir,pmt_pos,tmp); + + get_path_length(tmp,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o); + + /* Assume the particle is travelling at the speed of light. */ + return s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; +} + static double getKineticEnergy(double x, void *p) { return particle_get_energy(x, (particle *) p); @@ -664,10 +712,8 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t static double logp[MAX_PE], nll[MAX_PMTS]; double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; particle *p; - 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 wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov; double logp_path; - double a, b; - double tmp[3]; double mu_reflected; path *path; @@ -764,53 +810,6 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t continue; } - /* First, we try to compute the distance along the track where the - * PMT is at the Cerenkov angle. The reason for this is because for - * heavy particles like muons which don't scatter much, the - * probability distribution for getting a photon hit along the - * track looks kind of like a delta function, i.e. the PMT is only - * hit over a very narrow window when the angle between the track - * direction and the PMT is *very* close to the Cerenkov angle - * (it's not a perfect delta function since there is some width due - * to dispersion). In this case, it's possible that the numerical - * integration completely skips over the delta function and so - * predicts an expected charge of 0. - * - * To fix this, we only integrate 1 meter on both sides of the - * point along the track where the PMT is at the Cerenkov angle. */ - - /* 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 - * vector to the PMT at the start of the track. */ - cos_theta = DOT(dir,pmt_dir); - sin_theta = sqrt(1-pow(cos_theta,2)); - - /* Now, we compute the distance along the track where the PMT is at the - * Cerenkov angle. - * - * 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*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; - - /* Make sure that the point is somewhere along the track between 0 and - * `smax`. */ - if (s < 0) s = 0.0; - else if (s > smax) s = smax; - - /* Compute the endpoints of the integration. */ - a = s - 100.0; - b = s + 100.0; - - if (a < 0.0) a = 0.0; - if (b > range) b = range; - double q_indirect; integrate_path(path, i, mu_direct+i, &q_indirect, &result); mu_indirect += q_indirect; @@ -823,19 +822,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t * in the denominator. Therefore, we estimate the time using * the same technique as in get_total_charge_approx(). See that * function for more details. */ - n_h2o = get_index_snoman_h2o(wavelength0); - - /* Compute the position of the particle at a distance `s` along the track. */ - tmp[0] = pos[0] + s*dir[0]; - tmp[1] = pos[1] + s*dir[1]; - tmp[2] = pos[2] + s*dir[2]; - - SUB(pmt_dir,pmts[i].pos,tmp); - - get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); - - /* Assume the particle is travelling at the speed of light. */ - ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; + ts[i] = t0 + guess_time(pos,dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov); } if (id == IDP_E_MINUS || id == IDP_E_PLUS) { |