diff options
Diffstat (limited to 'src/muon.c')
| -rw-r--r-- | src/muon.c | 169 |
1 files changed, 2 insertions, 167 deletions
@@ -151,7 +151,7 @@ err: return -1; } -double get_range(double T, double rho) +double muon_get_range(double T, double rho) { /* Returns the approximate range a muon with kinetic energy `T` will travel * in water before losing all of its energy. This range is interpolated @@ -172,70 +172,7 @@ double get_range(double T, double rho) return gsl_spline_eval(spline_range, T, acc_range)/rho; } -muon_energy *muon_init_energy(double T0, double rho, size_t n) -{ - /* Returns a struct with arrays of the muon position and kinetic energy. - * This struct can then be passed to muon_get_energy() to interpolate the - * muon's kinetic energy at any point along the track. For example: - * - * muon_energy *m = muon_init_energy(1000.0, HEAVY_WATER_DENSITY, 100); - * double T = muon_get_energy(x, m); - * - * To compute the kinetic energy as a function of distance we need to solve - * the following integral equation: - * - * T(x) = T0 - \int_0^x dT(T(x))/dx - * - * which we solve by dividing the track up into `n` segments and then - * numerically computing the energy at each step. */ - size_t i; - double range, dEdx; - - /* Get the range of the muon. */ - range = get_range(T0, rho); - - muon_energy *m = malloc(sizeof(muon_energy)); - m->x = malloc(sizeof(double)*n); - m->T = malloc(sizeof(double)*n); - m->n = n; - - m->x[0] = 0; - m->T[0] = T0; - - /* Now we loop over the points along the track and assume dE/dx is constant - * between points. */ - for (i = 1; i < n; i++) { - m->x[i] = range*i/(n-1); - dEdx = get_dEdx(m->T[i-1], rho); - m->T[i] = m->T[i-1] - dEdx*(m->x[i]-m->x[i-1]); - } - - return m; -} - -double muon_get_energy(double x, muon_energy *m) -{ - /* Returns the approximate kinetic energy of a muon in water after - * travelling `x` cm with an initial kinetic energy `T`. - * - * Return value is in MeV. */ - double T; - - T = interp1d(x,m->x,m->T,m->n); - - if (T < 0) return 0; - - return T; -} - -void muon_free_energy(muon_energy *m) -{ - free(m->x); - free(m->T); - free(m); -} - -double get_dEdx(double T, double rho) +double muon_get_dEdx(double T, double rho) { /* Returns the approximate dE/dx for a muon in water with kinetic energy * `T`. @@ -255,105 +192,3 @@ 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; - - 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_response(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); -} |