diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/electron.c | 39 | ||||
| -rw-r--r-- | src/electron.h | 8 | ||||
| -rw-r--r-- | src/likelihood.c | 149 | ||||
| -rw-r--r-- | src/likelihood.h | 12 | ||||
| -rw-r--r-- | src/muon.c | 110 | ||||
| -rw-r--r-- | src/muon.h | 14 | ||||
| -rw-r--r-- | src/proton.h | 9 |
7 files changed, 281 insertions, 60 deletions
diff --git a/src/electron.c b/src/electron.c index 261ba39..c6b5ed8 100644 --- a/src/electron.c +++ b/src/electron.c @@ -161,6 +161,45 @@ static double electron_get_angular_pdf_norm(double alpha, double beta, double mu return result; } +double electron_get_angular_pdf_delta_ray(double cos_theta, double alpha, double beta, double mu) +{ + /* Returns the probability density that a photon in the wavelength range + * 200 nm - 800 nm is emitted at an angle cos_theta. + * + * The angular distribution is modelled by the pdf: + * + * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) + * + * where alpha and beta are constants which depend on the initial energy of + * the particle. + * + * There is no nice closed form solution for the integral of this function, + * so we just compute it on the fly. To make things faster, we keep track + * of the last alpha, beta, and mu parameters that were passed in so we + * avoid recomputing the normalization factor. */ + size_t i; + static double last_alpha, last_beta, last_mu, norm; + static int first = 1; + static double x[10000], y[10000]; + + if (first || alpha != last_alpha || beta != last_beta || mu != last_mu) { + norm = electron_get_angular_pdf_norm(alpha, beta, mu); + + last_alpha = alpha; + last_beta = beta; + last_mu = mu; + + for (i = 0; i < LEN(x); i++) { + x[i] = -1.0 + 2.0*i/(LEN(x)-1); + y[i] = electron_get_angular_pdf_no_norm(x[i], alpha, beta, mu)/norm; + } + + first = 0; + } + + return interp1d(cos_theta,x,y,LEN(x)); +} + double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu) { /* Returns the probability density that a photon in the wavelength range diff --git a/src/electron.h b/src/electron.h index cb3f24e..20e5ea2 100644 --- a/src/electron.h +++ b/src/electron.h @@ -3,9 +3,17 @@ #include <stddef.h> /* for size_t */ +/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy + * lost to radiation for electrons. + * + * FIXME: This is just a rough estimate, should use an energy dependent + * quantity from simulation. */ +#define ELECTRON_PHOTONS_PER_MEV 400.0 + void electron_get_position_distribution_parameters(double T0, double *a, double *b); double electron_get_angular_distribution_alpha(double T0); double electron_get_angular_distribution_beta(double T0); +double electron_get_angular_pdf_delta_ray(double cos_theta, double alpha, double beta, double mu); double electron_get_angular_pdf(double cos_theta, double alpha, double beta, double mu); double electron_get_range(double T, double rho); double electron_get_dEdx_rad(double T, double rho); diff --git a/src/likelihood.c b/src/likelihood.c index 786b221..8e12f57 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -42,16 +42,19 @@ particle *particle_init(int id, double T0, size_t n) * which we solve by dividing the track up into `n` segments and then * numerically computing the energy at each step. */ size_t i; - double dEdx; + double dEdx, rad; particle *p = malloc(sizeof(particle)); p->id = id; p->x = malloc(sizeof(double)*n); p->T = malloc(sizeof(double)*n); p->n = n; - p->rad = 0.0; p->a = 0.0; p->b = 0.0; + p->shower_photons = 0.0; + p->delta_ray_a = 0.0; + p->delta_ray_b = 0.0; + p->delta_ray_photons = 0.0; p->x[0] = 0; p->T[0] = T0; @@ -67,14 +70,17 @@ particle *particle_init(int id, double T0, size_t n) p->a = electron_get_angular_distribution_alpha(T0); p->b = electron_get_angular_distribution_beta(T0); + p->delta_ray_photons = 0.0; + /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ + rad = 0.0; for (i = 1; i < n; i++) { p->x[i] = p->range*i/(n-1); dEdx = electron_get_dEdx(p->T[i-1], WATER_DENSITY); p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]); dEdx = electron_get_dEdx_rad(p->T[i-1], WATER_DENSITY); - p->rad += dEdx*(p->x[i]-p->x[i-1]); + rad += dEdx*(p->x[i]-p->x[i-1]); } /* Make sure that the energy is zero at the last step. This is so that * when we try to bisect the point along the track where the speed of @@ -85,6 +91,8 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + p->shower_photons = rad*ELECTRON_PHOTONS_PER_MEV; + break; case IDP_MU_MINUS: case IDP_MU_PLUS: @@ -96,14 +104,21 @@ particle *particle_init(int id, double T0, size_t n) * load the correct table based on the media. */ p->range = muon_get_range(T0, HEAVY_WATER_DENSITY); + p->a = muon_get_angular_distribution_alpha(T0); + p->b = muon_get_angular_distribution_beta(T0); + + muon_get_delta_ray_distribution_parameters(T0, &p->delta_ray_a, &p->delta_ray_b); + p->delta_ray_photons = muon_get_delta_ray_photons(T0); + /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ + rad = 0.0; for (i = 1; i < n; i++) { p->x[i] = p->range*i/(n-1); dEdx = muon_get_dEdx(p->T[i-1], HEAVY_WATER_DENSITY); p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]); dEdx = muon_get_dEdx_rad(p->T[i-1], WATER_DENSITY); - p->rad += dEdx*(p->x[i]-p->x[i-1]); + rad += dEdx*(p->x[i]-p->x[i-1]); } /* Make sure that the energy is zero at the last step. This is so that * when we try to bisect the point along the track where the speed of @@ -114,6 +129,8 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + p->shower_photons = rad*MUON_PHOTONS_PER_MEV; + break; case IDP_PROTON: p->mass = PROTON_MASS; @@ -121,14 +138,17 @@ particle *particle_init(int id, double T0, size_t n) * tables for heavy water for protons. */ p->range = proton_get_range(T0, WATER_DENSITY); + p->delta_ray_photons = 0.0; + /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ + rad = 0.0; for (i = 1; i < n; i++) { p->x[i] = p->range*i/(n-1); dEdx = proton_get_dEdx(p->T[i-1], WATER_DENSITY); p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]); dEdx = proton_get_dEdx_rad(p->T[i-1], WATER_DENSITY); - p->rad += dEdx*(p->x[i]-p->x[i-1]); + rad += dEdx*(p->x[i]-p->x[i-1]); } /* Make sure that the energy is zero at the last step. This is so that * when we try to bisect the point along the track where the speed of @@ -139,6 +159,8 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + p->shower_photons = rad*PROTON_PHOTONS_PER_MEV; + break; default: fprintf(stderr, "unknown particle id %i\n", id); @@ -170,9 +192,9 @@ void particle_free(particle *p) free(p); } -static double get_expected_charge_shower(double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double alpha, double beta, double rad) +static double get_expected_charge_shower(particle *p, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double *q_delta_ray, double *q_indirect_delta_ray) { - double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter; + double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter, constant; SUB(pmt_dir,pmt_pos,pos); @@ -180,6 +202,9 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p cos_theta_pmt = DOT(pmt_dir,pmt_normal); + charge = 0.0; + *q_delta_ray = 0.0; + *q_indirect_delta_ray = 0.0; *reflected = 0.0; if (cos_theta_pmt > 0) return 0.0; @@ -201,11 +226,26 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; - charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI); + constant = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*omega/(2*M_PI); + + /* Note: We assume here that the peak of the angular distribution is at the + * Cerenkov angle for a particle with beta = 1. This seems to be an OK + * assumption for high energy showers, but is not exactly correct for + * showers from lower energy electrons or for delta rays from lower energy + * muons. It seems good enough, but in the future it would be nice to + * parameterize this. */ + if (p->shower_photons > 0) + charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/n_d2o); + + if (p->delta_ray_photons > 0) + *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o); scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); *reflected = (f_reflec + 1.0 - scatter)*charge; + *q_indirect_delta_ray = (f_reflec + 1.0 - scatter)*(*q_delta_ray); + + *q_delta_ray *= f*scatter; return f*charge*scatter; } @@ -345,14 +385,14 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); } -static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad) +static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma) { size_t i; static double qs[MAX_NPOINTS]; static double q_indirects[MAX_NPOINTS]; static double ts[MAX_NPOINTS]; static double ts2[MAX_NPOINTS]; - double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o; + double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray; if (n > MAX_NPOINTS) { fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n"); @@ -364,6 +404,8 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); + double range = x[n-1] - x[0]; + for (i = 0; i < n; i++) { pos[0] = pos0[0] + x[i]*dir0[0]; pos[1] = pos0[1] + x[i]*dir0[1]; @@ -373,8 +415,9 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 t = x[i]/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; - qs[i] = get_expected_charge_shower(pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, p->a, p->b, rad)*pdf[i]; - q_indirects[i] *= pdf[i]; + qs[i] = get_expected_charge_shower(p, pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, &q_delta_ray, &q_indirect_delta_ray); + qs[i] = qs[i]*pdf[i] + q_delta_ray/range; + q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/range; ts[i] = t*qs[i]; ts2[i] = t*t*qs[i]; } @@ -593,10 +636,11 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa charge = abs_prob*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); - /* Add expected number of photons from electromagnetic shower for electrons - * and positrons. */ - if (p->id == IDP_E_MINUS || p->id == IDP_E_PLUS) - charge += abs_prob*get_weighted_quantum_efficiency()*p->rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); + /* Add expected number of photons from electromagnetic shower. */ + if (p->shower_photons > 0) + charge += abs_prob*get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); + if (p->delta_ray_photons > 0) + charge += abs_prob*get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/n_d2o)*omega/(2*M_PI); *mu_reflected = f_reflected*charge; @@ -891,6 +935,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast npoints = (size_t) (range/dx + 0.5); if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS; + if (npoints > MAX_NPOINTS) npoints = MAX_NPOINTS; /* Calculate total energy */ E0 = v[j].T0 + p->mass; @@ -903,32 +948,43 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast if (!fast) { path = path_init(v[j].pos, v[j].dir, v[j].T0, range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass); - if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) { + switch (v[j].id) { + case IDP_E_MINUS: + case IDP_E_PLUS: electron_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b); + break; + case IDP_MU_MINUS: + case IDP_MU_PLUS: + muon_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b); + break; + default: + fprintf(stderr, "unknown particle id %i\n", v[j].id); + exit(1); + } - /* Lower and upper limits to integrate for the electromagnetic - * shower. */ - xlo = 0.0; - xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); + /* Lower and upper limits to integrate for the electromagnetic + * shower. */ + xlo = 0.0; + xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); - /* Number of points to sample along the longitudinal direction for - * the electromagnetic shower. */ - npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5); + /* Number of points to sample along the longitudinal direction for + * the electromagnetic shower. */ + npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5); - if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS; + if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS; + if (npoints_shower > MAX_NPOINTS) npoints_shower = MAX_NPOINTS; - x = malloc(sizeof(double)*npoints_shower); - pdf = malloc(sizeof(double)*npoints_shower); - for (i = 0; i < npoints_shower; i++) { - x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1); - pdf[i] = gamma_pdf(x[i],pos_a,pos_b); - } + x = malloc(sizeof(double)*npoints_shower); + pdf = malloc(sizeof(double)*npoints_shower); + for (i = 0; i < npoints_shower; i++) { + x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1); + pdf[i] = gamma_pdf(x[i],pos_a,pos_b); + } - /* Make sure the PDF is normalized. */ - double norm = trapz(pdf,x[1]-x[0],npoints_shower); - for (i = 0; i < npoints_shower; i++) { - pdf[i] /= norm; - } + /* Make sure the PDF is normalized. */ + double norm = trapz(pdf,x[1]-x[0],npoints_shower); + for (i = 0; i < npoints_shower; i++) { + pdf[i] /= norm; } } @@ -982,27 +1038,16 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast ts[i][j] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov); } - if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) { - /* If we are fitting for an electron or positron, we need to - * estimate the expected number of photons produced by the - * electromagnetic shower. - * - * Technically we should also have a term for the electromagnetic - * shower produced by muons, but for the energy range I'm - * considering right now, it is very small. */ - integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j],p->rad); - mu_indirect[j] += q_indirect; - ts_shower[i][j] = v[j].t0 + result; - } + integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]); + mu_indirect[j] += q_indirect; + ts_shower[i][j] = v[j].t0 + result; } if (!fast) { path_free(path); - if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) { - free(x); - free(pdf); - } + free(x); + free(pdf); } particle_free(p); diff --git a/src/likelihood.h b/src/likelihood.h index 84fa7ff..4b72742 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -45,13 +45,6 @@ #define GTVALID 400.0 #define BETA_MIN 0.8 -/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy - * lost to radiation for electrons. - * - * FIXME: This is just a rough estimate, should use an energy dependent - * quantity from simulation. */ -#define PHOTONS_PER_MEV 400.0 - /* Maximum number of vertices to fit. */ #define MAX_VERTICES 10 @@ -70,12 +63,15 @@ typedef struct particle { int id; double mass; double range; - double rad; double *x; double *T; size_t n; double a; double b; + double shower_photons; + double delta_ray_a; + double delta_ray_b; + double delta_ray_photons; } particle; particle *particle_init(int id, double T0, size_t n); @@ -40,6 +40,116 @@ static gsl_spline *spline_range; static const double MUON_CRITICAL_ENERGY_H2O = 1029.0e6; static const double MUON_CRITICAL_ENERGY_D2O = 967.0e3; +void muon_get_position_distribution_parameters(double T0, double *a, double *b) +{ + /* Computes the gamma distribution parameters describing the longitudinal + * profile of the photons generated from an electromagnetic shower. + * + * From the PDG "Passage of Particles" section 32.5: + * + * "The mean longitudinal profile of the energy deposition in an + * electromagnetic cascade is reasonably well described by a gamma + * distribution." + * + * Here we use a slightly different form of the gamma distribution: + * + * f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a) + * + * I determined the b parameter by simulating high energy muons using + * RAT-PAC and determined that it's roughly equal to the radiation length. + * To calculate the a parameter we use the formula from the PDG, i.e. + * + * tmax = (a-1)/b = ln(E/E_C) - 0.5 + * + * Therefore, we calculate a as: + * + * a = tmax*b+1. + * + * `T` should be in units of MeV. + * + * Example: + * + * double a, b; + * muon_get_position_distribution_parameters(1000.0, &a, &b); + * double pdf = gamma_pdf(x, a, b); + * + * See http://pdg.lbl.gov/2014/reviews/rpp2014-rev-passage-particles-matter.pdf. + * + * FIXME: Double check that this is correct for muons. */ + double tmax; + + *b = RADIATION_LENGTH; + tmax = log(T0/MUON_CRITICAL_ENERGY_D2O) - 0.5; + *a = fmax(1.1,tmax*(*b)/RADIATION_LENGTH + 1); +} + +double muon_get_angular_distribution_alpha(double T0) +{ + /* To describe the angular distribution of photons in an electromagnetic + * shower I came up with the heuristic form: + * + * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) + * + * I simulated high energy muons using RAT-PAC in heavy water to fit + * for the alpha and beta parameters as a function of energy and determined + * that a reasonably good fit to the data was: + * + * alpha = c0 + c1/log(c2*T0 + c3) + * + * where T0 is the initial energy of the muon in MeV and c0, c1, c2, + * and c3 are constants which I fit out. */ + return 8.238633e-01 + 3.896665e-03/log(1.581060e-05*T0 + 9.991576e-01); +} + +double muon_get_angular_distribution_beta(double T0) +{ + /* To describe the angular distribution of photons in an electromagnetic + * shower I came up with the heuristic form: + * + * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) + * + * I simulated high energy muons using RAT-PAC in heavy water to fit + * for the alpha and beta parameters as a function of energy and determined + * that a reasonably good fit to the data was: + * + * beta = c0 + c1/log(c2*T0 + c3) + * + * where T0 is the initial energy of the muon in MeV and c0, c1, c2, + * and c3 are constants which I fit out. */ + return 2.236058e-01 + 5.376336e-03/log(1.200215e-05*T0 + 1.002540e+00); +} + +double muon_get_delta_ray_photons(double T0) +{ + /* Returns the number of photons in the range 200-800 nm produced by delta + * rays. + * + * This comes from a simple polynomial fit to the number of photons + * generated by simulating muons using RAT-PAC in heavy water. */ + return fmax(0.0,-7532.39 + 39.4548*T0); +} + +void muon_get_delta_ray_distribution_parameters(double T0, double *a, double *b) +{ + /* To describe the angular distribution of photons from delta rays I use + * the same heuristic form as used for electromagnetic showers: + * + * f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) + * + * I simulated high energy muons using RAT-PAC in heavy water to fit + * for the alpha and beta parameters as a function of energy and determined + * that a reasonably good fit to the data was: + * + * beta = c0 + c1/log(c2*T0 + c3) + * + * where T0 is the initial energy of the muon in MeV and c0, c1, c2, + * and c3 are constants which I fit out. */ + *a = 3.463093e-01 + 1.110835e-02/log(5.662948e-06*T0 + 1.009215e+00); + *b = 2.297358e-01 + 4.085721e-03/log(8.218774e-06*T0 + 1.007135e+00); + + return; +} + static int init() { int i, j; @@ -5,6 +5,20 @@ #define EULER_CONSTANT 0.57721 +/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy + * lost to radiation for electrons. + * + * FIXME: This is just a rough estimate, should use an energy dependent + * quantity from simulation. + * + * FIXME: Actually determine what this is. */ +#define MUON_PHOTONS_PER_MEV 7368.0 + +void muon_get_position_distribution_parameters(double T0, double *a, double *b); +double muon_get_angular_distribution_alpha(double T0); +double muon_get_angular_distribution_beta(double T0); +void muon_get_delta_ray_distribution_parameters(double T0, double *a, double *b); +double muon_get_delta_ray_photons(double T0); double muon_get_range(double T, double rho); double muon_get_dEdx_rad(double T, double rho); double muon_get_dEdx(double T, double rho); diff --git a/src/proton.h b/src/proton.h index 008f835..cffeb9f 100644 --- a/src/proton.h +++ b/src/proton.h @@ -3,6 +3,15 @@ #include <stddef.h> /* for size_t */ +/* Number of photons in the range 200 nm - 800 nm generated per MeV of energy + * lost to radiation for electrons. + * + * FIXME: This is just a rough estimate, should use an energy dependent + * quantity from simulation. + * + * FIXME: Actually determine what this is. */ +#define PROTON_PHOTONS_PER_MEV 400.0 + double proton_get_range(double T, double rho); double proton_get_dEdx_rad(double T, double rho); double proton_get_dEdx(double T, double rho); |