From 8447870e721dd738bce12b45e732c9cc01dc2595 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Sun, 11 Nov 2018 13:22:18 -0600 Subject: update likelihood function to fit electrons! To characterize the angular distribution of photons from an electromagnetic shower I came up with the following functional form: f(cos_theta) ~ exp(-abs(cos_theta-mu)^alpha/beta) and fit this to data simulated using RAT-PAC at several different energies. I then fit the alpha and beta coefficients as a function of energy to the functional form: alpha = c0 + c1/log(c2*T0 + c3) beta = c0 + c1/log(c2*T0 + c3). where T0 is the initial energy of the electron in MeV and c0, c1, c2, and c3 are parameters which I fit. The longitudinal distribution of the photons generated from an electromagnetic shower is described by a gamma distribution: f(x) = x**(a-1)*exp(-x/b)/(Gamma(a)*b**a). This parameterization comes from the PDG "Passage of particles through matter" section 32.5. I also fit the data from my RAT-PAC simulation, but currently I am not using it, and instead using a simpler form to calculate the coefficients from the PDG (although I estimated the b parameter from the RAT-PAC data). I also sped up the calculation of the solid angle by making a lookup table since it was taking a significant fraction of the time to compute the likelihood function. --- src/proton.c | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) (limited to 'src/proton.c') diff --git a/src/proton.c b/src/proton.c index 922c552..f5a6b61 100644 --- a/src/proton.c +++ b/src/proton.c @@ -18,9 +18,12 @@ static int initialized = 0; -static double *x, *dEdx, *csda_range; +static double *x, *dEdx_rad, *dEdx, *csda_range; static size_t size; +static gsl_interp_accel *acc_dEdx_rad; +static gsl_spline *spline_dEdx_rad; + static gsl_interp_accel *acc_dEdx; static gsl_spline *spline_dEdx; @@ -71,6 +74,7 @@ static int init() } x = malloc(sizeof(double)*n); + dEdx_rad = malloc(sizeof(double)*n); dEdx = malloc(sizeof(double)*n); csda_range = malloc(sizeof(double)*n); size = n; @@ -103,6 +107,9 @@ static int init() case 0: x[n] = value; break; + case 2: + dEdx_rad[n] = value; + break; case 3: dEdx[n] = value; break; @@ -119,6 +126,10 @@ static int init() fclose(f); + acc_dEdx_rad = gsl_interp_accel_alloc(); + spline_dEdx_rad = gsl_spline_alloc(gsl_interp_linear, size); + gsl_spline_init(spline_dEdx_rad, x, dEdx_rad, size); + acc_dEdx = gsl_interp_accel_alloc(); spline_dEdx = gsl_spline_alloc(gsl_interp_linear, size); gsl_spline_init(spline_dEdx, x, dEdx, size); @@ -158,6 +169,27 @@ double proton_get_range(double T, double rho) return gsl_spline_eval(spline_range, T, acc_range)/rho; } +double proton_get_dEdx_rad(double T, double rho) +{ + /* Returns the approximate radiative dE/dx for a proton 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); + } + } + + if (T < spline_dEdx_rad->x[0]) return spline_dEdx_rad->y[0]; + + return gsl_spline_eval(spline_dEdx_rad, T, acc_dEdx_rad)*rho; +} + double proton_get_dEdx(double T, double rho) { /* Returns the approximate dE/dx for a proton in water with kinetic energy -- cgit