diff options
| -rw-r--r-- | src/likelihood.c | 7 | ||||
| -rw-r--r-- | src/muon.c | 7 | ||||
| -rw-r--r-- | src/optics.c | 63 | ||||
| -rw-r--r-- | src/optics.h | 2 |
4 files changed, 75 insertions, 4 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 99258a9..000389b 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -146,7 +146,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, l_h2o, l_d2o, 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; /* Calculate beta at the start of the track. */ E0 = T0 + MUON_MASS; @@ -247,13 +247,16 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy 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(); get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); + l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; + /* 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)*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); + 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); } typedef struct betaRootParams { @@ -258,7 +258,7 @@ double get_dEdx(double T, double rho) 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, l_h2o, l_d2o; + 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; @@ -298,8 +298,11 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl 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); - return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); + 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); } diff --git a/src/optics.c b/src/optics.c index 57c9d57..5f01089 100644 --- a/src/optics.c +++ b/src/optics.c @@ -24,6 +24,10 @@ static double absorption_coefficient_h2o[7] = {1e-5, 1e-5, 1.3e-5, 6.0e-5, 2.0e- static double absorption_coefficient_d2o_wavelengths[7] = {254.0, 313.0, 366.0, 406.0, 436.0, 548.0, 578.0}; static double absorption_coefficient_d2o[7] = {1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5}; +/* From prod/media.dat for Acrylic (good). */ +static double absorption_coefficient_acrylic_wavelengths[10] = {300.0, 310.0, 320.0, 330.0, 340.0, 350.0, 360.0, 380.0, 400.0, 450.0}; +static double absorption_coefficient_acrylic[10] = {0.3641, 0.0966, 0.0667, 0.049, 0.035, 0.0261, 0.0194, 0.0115, 0.0086, 0.0071}; + static int initialized = 0; /* D2O absorption length weighted by the Cerenkov spectrum and the quantum @@ -32,6 +36,10 @@ static double absorption_length_d2o_weighted = 0.0; /* H2O absorption length weighted by the Cerenkov spectrum and the quantum * efficiency. */ static double absorption_length_h2o_weighted = 0.0; +/* Acrylic absorption length weighted by the Cerenkov spectrum and the quantum + * efficiency. */ +static double absorption_length_acrylic_weighted = 0.0; + /* From Table 4 in the paper. */ static double A0 = 0.243905091; @@ -80,6 +88,17 @@ static double gsl_absorption_length_snoman_h2o(double wavelength, void *params) return qe*get_absorption_length_snoman_h2o(wavelength)/pow(wavelength,2); } +static double gsl_absorption_length_snoman_acrylic(double wavelength, void *params) +{ + /* Returns the absorption length in acrylic weighted by the quantum + * efficiency and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_absorption_length_snoman_acrylic(wavelength)/pow(wavelength,2); +} + static double gsl_cerenkov(double wavelength, void *params) { /* Returns the quantum efficiency multiplied by the Cerenkov spectrum. */ @@ -133,6 +152,17 @@ int optics_init(void) exit(1); } + F.function = &gsl_absorption_length_snoman_acrylic; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + absorption_length_acrylic_weighted = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + gsl_integration_cquad_workspace_free(w); initialized = 1; @@ -206,6 +236,39 @@ double get_absorption_length_snoman_d2o(double wavelength) return 1.0/gsl_spline_eval(spline, wavelength, acc); } +double get_weighted_absorption_length_snoman_acrylic(void) +{ + /* Returns the weighted absorption length in acrylic in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted absorption length hasn't been initialized!\n"); + exit(1); + } + + return absorption_length_acrylic_weighted; +} + +double get_absorption_length_snoman_acrylic(double wavelength) +{ + /* Returns the absorption length in acrylic in cm. */ + static gsl_spline *spline; + static gsl_interp_accel *acc; + + if (!spline) { + spline = gsl_spline_alloc(gsl_interp_linear, LEN(absorption_coefficient_acrylic_wavelengths)); + gsl_spline_init(spline, absorption_coefficient_acrylic_wavelengths, absorption_coefficient_acrylic, LEN(absorption_coefficient_acrylic_wavelengths)); + acc = gsl_interp_accel_alloc(); + } + + if (wavelength < absorption_coefficient_acrylic_wavelengths[0]) + wavelength = absorption_coefficient_acrylic_wavelengths[0]; + + if (wavelength > absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1]) + wavelength = absorption_coefficient_acrylic_wavelengths[LEN(absorption_coefficient_acrylic_wavelengths)-1]; + + return 1.0/gsl_spline_eval(spline, wavelength, acc); +} + double get_index(double p, double wavelength, double T) { /* Returns the index of refraction of pure water for a given density, diff --git a/src/optics.h b/src/optics.h index fdeefdd..06b89ca 100644 --- a/src/optics.h +++ b/src/optics.h @@ -6,6 +6,8 @@ double get_weighted_absorption_length_snoman_h2o(void); double get_absorption_length_snoman_h2o(double wavelength); double get_weighted_absorption_length_snoman_d2o(void); double get_absorption_length_snoman_d2o(double wavelength); +double get_weighted_absorption_length_snoman_acrylic(void); +double get_absorption_length_snoman_acrylic(double wavelength); double get_index(double p, double wavelength, double T); double get_index_snoman_h2o(double wavelength); double get_index_snoman_d2o(double wavelength); |