diff options
Diffstat (limited to 'muon.c')
| -rw-r--r-- | muon.c | 260 |
1 files changed, 0 insertions, 260 deletions
@@ -1,260 +0,0 @@ -#include <stdio.h> -#include <errno.h> -#include <string.h> -#include <stdlib.h> -#include <gsl/gsl_errno.h> -#include <gsl/gsl_spline.h> -#include <math.h> -#include "optics.h" -#include "quantum_efficiency.h" -#include "solid_angle.h" -#include "pdg.h" -#include "vector.h" -#include "muon.h" -#include "sno.h" -#include "scattering.h" - -static int initialized = 0; - -static double *x, *dEdx, *csda_range; -static size_t size; - -static gsl_interp_accel *acc_dEdx; -static gsl_spline *spline_dEdx; - -static gsl_interp_accel *acc_range; -static gsl_spline *spline_range; - -static const double MUON_CRITICAL_ENERGY = 1.029e6; - -static int init() -{ - int i, j; - char line[256]; - char *str; - double value; - int n; - - FILE *f = fopen("muE_water_liquid.txt", "r"); - - if (!f) { - fprintf(stderr, "failed to open muE_water_liquid.txt: %s", strerror(errno)); - return -1; - } - - i = 0; - n = 0; - /* For the first pass, we just count how many values there are. */ - while (fgets(line, sizeof(line), f)) { - size_t len = strlen(line); - if (len && (line[len-1] != '\n')) { - fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line); - goto err; - } - - i += 1; - - /* Skip the first 10 lines since it's just a header. */ - if (i <= 10) continue; - - if (!len) continue; - else if (line[0] == '#') continue; - else if (strstr(line, "Minimum ionization")) continue; - else if (strstr(line, "Muon critical energy")) continue; - - str = strtok(line," \n"); - - while (str) { - value = strtod(str, NULL); - str = strtok(NULL," \n"); - } - - n += 1; - } - - x = malloc(sizeof(double)*n); - dEdx = malloc(sizeof(double)*n); - csda_range = malloc(sizeof(double)*n); - size = n; - - i = 0; - n = 0; - /* Now, we actually store the values. */ - rewind(f); - while (fgets(line, sizeof(line), f)) { - size_t len = strlen(line); - if (len && (line[len-1] != '\n')) { - fprintf(stderr, "got incomplete line on line %i: '%s'\n", i, line); - goto err; - } - - i += 1; - - /* Skip the first 10 lines since it's just a header. */ - if (i <= 10) continue; - - if (!len) continue; - else if (line[0] == '#') continue; - else if (strstr(line, "Minimum ionization")) continue; - else if (strstr(line, "Muon critical energy")) continue; - - str = strtok(line," \n"); - - j = 0; - while (str) { - value = strtod(str, NULL); - switch (j) { - case 0: - x[n] = value; - break; - case 7: - dEdx[n] = value; - break; - case 8: - csda_range[n] = value; - break; - } - j += 1; - str = strtok(NULL," \n"); - } - - n += 1; - } - - fclose(f); - - acc_dEdx = gsl_interp_accel_alloc(); - spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size); - gsl_spline_init(spline_dEdx, x, dEdx, size); - - acc_range = gsl_interp_accel_alloc(); - spline_range = gsl_spline_alloc(gsl_interp_linear, size); - gsl_spline_init(spline_range, x, csda_range, size); - - initialized = 1; - - return 0; - -err: - fclose(f); - - return -1; -} - -double 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 - * based on data from the PDG which uses the continuous slowing down - * approximation. - * - * `T` should be in MeV, and `rho` should be in g/cm^3. - * - * Return value is in cm. - * - * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */ - if (!initialized) { - if (init()) { - exit(1); - } - } - - return gsl_spline_eval(spline_range, T, acc_range)/rho; -} - -double get_T(double T0, double x, double rho) -{ - /* Returns the approximate kinetic energy of a muon in water after - * travelling `x` cm with an initial kinetic energy `T`. - * - * `T` should be in MeV, `x` in cm, and `rho` in g/cm^3. - * - * Return value is in MeV. - * - * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */ - double a, b, range, T; - - if (!initialized) { - if (init()) { - exit(1); - } - } - - range = get_range(T0, rho); - - /* This comes from Equation 33.42 in the PDG Passage of Particles Through - * Matter article. */ - b = log(1 + T0/MUON_CRITICAL_ENERGY)/range; - /* Now we compute the ionization energy loss from the known range and b. */ - a = b*T0/(exp(b*range)-1.0); - - /* Compute the kinetic energy after travelling a distance `x` in the - * continuous slowing down approximation. */ - T = -a/b + (T0+a/b)*exp(-b*x); - - if (T < 0) return 0; - - return T; -} - -double get_dEdx(double T, double rho) -{ - /* Returns the approximate dE/dx for a muon in water with kinetic energy - * `T`. - * - * `T` should be in MeV and `rho` in g/cm^3. - * - * Return value is in MeV/cm. - * - * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */ - if (!initialized) { - if (init()) { - exit(1); - } - } - - return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho; -} - -double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r) -{ - double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R; - - z = 1.0; - - SUB(pmt_dir,pmt_pos,pos); - normalize(pmt_dir); - - if (DOT(pmt_dir,pmt_normal) > 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); - - if (R <= AV_RADIUS) { - rho = HEAVY_WATER_DENSITY; - } else { - rho = WATER_DENSITY; - } - - /* FIXME: I just calculate delta assuming 400 nm light. */ - wavelength0 = 400.0; - n = get_index(rho, wavelength0, 10.0); - - if (beta < 1/n) return 0; - - /* FIXME: is this formula valid for muons? */ - theta0 = get_scattering_rms(x,p,beta,z,rho); - - /* FIXME: add angular response and scattering/absorption. */ - return 2*omega*2*M_PI*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0)/(sqrt(2*M_PI)*theta0); -} |