diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 43 | ||||
| -rw-r--r-- | src/likelihood.h | 9 | ||||
| -rw-r--r-- | src/muon.c | 51 | ||||
| -rw-r--r-- | src/muon.h | 1 |
4 files changed, 82 insertions, 22 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 000389b..daad9cf 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -105,6 +105,16 @@ static double gsl_muon_time(double x, void *params) return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS); } +static double gsl_muon_reflected_charge(double x, void *params) +{ + intParams *pars = (intParams *) params; + double dir[3], pos[3], T, t, theta0; + + path_eval(pars->p, x, pos, dir, &T, &t, &theta0); + + return get_expected_reflected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS); +} + static double gsl_muon_charge(double x, void *params) { intParams *pars = (intParams *) params; @@ -115,7 +125,7 @@ static double gsl_muon_charge(double x, void *params) return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS); } -double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t) +double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t, double *mu_reflected) { /* Returns the approximate expected number of photons seen by PMT `i` using * an analytic formula. @@ -146,7 +156,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy * * `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, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0; + double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected; /* Calculate beta at the start of the track. */ E0 = T0 + MUON_MASS; @@ -244,6 +254,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy double frac = sqrt(2)*n_d2o*x*beta0*theta0; f = get_weighted_pmt_response(acos(-cos_theta_pmt)); + f_reflected = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt)); absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0); absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0); @@ -256,7 +267,11 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy /* 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; - return f*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); + double charge = f*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; + + return f*charge; } typedef struct betaRootParams { @@ -333,14 +348,13 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl { size_t i, j, nhit; intParams params; - double total_charge; double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp; - double tmean = 0.0; muon_energy *m; double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o; double logp_path; double a, b; double tmp[3]; + double mu_reflected; double mu_direct[MAX_PMTS]; double ts[MAX_PMTS]; @@ -374,15 +388,16 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl else smax = 0.0; - total_charge = 0.0; + mu_indirect = 0.0; for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; params.i = i; if (fast) { - mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i]); + mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i], &mu_reflected); ts[i] += t0; + mu_indirect += mu_reflected; } else { /* 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 @@ -442,6 +457,11 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); mu_direct[i] = result; + F.function = &gsl_muon_reflected_charge; + + gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals); + mu_indirect += result; + F.function = &gsl_muon_time; if (mu_direct[i] > 1e-9) { @@ -468,23 +488,16 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; } } - - tmean += ts[i]*mu_direct[i]; - - total_charge += mu_direct[i]; } path_free(params.p); muon_free_energy(m); - if (total_charge > 0) - tmean /= total_charge; - gsl_integration_cquad_workspace_free(w); mu_noise = DARK_RATE*GTVALID*1e-9; - mu_indirect = total_charge*CHARGE_FRACTION/10000.0; + mu_indirect *= CHARGE_FRACTION/10000.0; for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; diff --git a/src/likelihood.h b/src/likelihood.h index d53b777..b4b3be9 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -18,13 +18,8 @@ * value for n. */ #define MIN_RATIO -10 -/* Roughly equal to the fraction of light which is reflected at each PMT. - * - * Note: The proper thing to do here would be to separately integrate over the - * track to calculate the amount of reflected light using the data from the - * PMTR bank. However, this would make the likelihood function almost twice as - * long, so we just assume a constant fraction of the light is reflected now. */ -#define CHARGE_FRACTION 0.3 +/* The fraction of reflected light which is detected. */ +#define CHARGE_FRACTION 1.0 /* Dark rate of the PMTs (Hz). * @@ -256,6 +256,57 @@ double get_dEdx(double T, double rho) return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho; } +double get_expected_reflected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r) +{ + double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic; + + z = 1.0; + + SUB(pmt_dir,pmt_pos,pos); + + normalize(pmt_dir); + + cos_theta_pmt = DOT(pmt_dir,pmt_normal); + + if (cos_theta_pmt > 0) return 0; + + /* Calculate the cosine of the angle between the track direction and the + * vector to the PMT. */ + cos_theta = DOT(dir,pmt_dir); + + /* Calculate total energy */ + E = T + MUON_MASS; + p = sqrt(E*E - MUON_MASS*MUON_MASS); + beta = p/E; + + omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r); + + R = NORM(pos); + + /* FIXME: I just calculate delta assuming 400 nm light. */ + wavelength0 = 400.0; + + if (R <= AV_RADIUS) { + n = get_index_snoman_d2o(wavelength0); + } else { + n = get_index_snoman_h2o(wavelength0); + } + + if (beta < 1/n) return 0; + + f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt)); + + 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(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o); + + l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; + + return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); +} + double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r) { double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic; @@ -16,6 +16,7 @@ double muon_get_energy(double x, muon_energy *m); void muon_free_energy(muon_energy *m); double get_range(double T, double rho); double get_dEdx(double T, double rho); +double get_expected_reflected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r); double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r); #endif |