diff options
Diffstat (limited to 'muon.c')
| -rw-r--r-- | muon.c | 83 |
1 files changed, 67 insertions, 16 deletions
@@ -5,6 +5,14 @@ #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; @@ -95,8 +103,6 @@ static int init() j = 0; while (str) { value = strtod(str, NULL); - /* According to the file, the values are stored for wavelengths - * between 230 and 700 in 1 nm increments. */ switch (j) { case 0: x[n] = value; @@ -156,17 +162,17 @@ double get_range(double T, double rho) return gsl_spline_eval(spline_range, T, acc_range)/rho; } -double get_E(double T, double x, double rho) +double get_T(double T0, double x, double rho) { - /* Returns the approximate energy of a muon in water after travelling `x` - * cm with an initial kinetic energy `T`. + /* 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, E; + double a, b, range, T; if (!initialized) { if (init()) { @@ -174,19 +180,21 @@ double get_E(double T, double x, double rho) } } - range = gsl_spline_eval(spline_range, T, acc_range)/rho; - /* FIXME: Need to double check if it's kosher to be using kinetic energies - * here instead of the total energy. Equation 1 in the document uses the - * total energy, but here I'm using the critical energy in kinetic energy, - * so I should check to see if I need to convert both. */ - b = log(1 + T/MUON_CRITICAL_ENERGY)/range; - a = MUON_CRITICAL_ENERGY*b; + range = get_range(T0, rho); - E = T + a*(1-exp(b*x))/b; + /* 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); - if (E < 0) return 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); - return E; + if (T < 0) return 0; + + return T; } double get_dEdx(double T, double rho) @@ -207,3 +215,46 @@ double get_dEdx(double T, double rho) 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); +} |