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/electron.c | 195 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 194 insertions(+), 1 deletion(-) (limited to 'src/electron.c') diff --git a/src/electron.c b/src/electron.c index e8e24d3..261ba39 100644 --- a/src/electron.c +++ b/src/electron.c @@ -15,12 +15,17 @@ #include "scattering.h" #include "pmt_response.h" #include "misc.h" +#include +#include /* for fmax() */ 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; @@ -36,6 +41,165 @@ static gsl_spline *spline_range; static const double ELECTRON_CRITICAL_ENERGY_H2O = 78.33; static const double ELECTRON_CRITICAL_ENERGY_D2O = 78.33; +void electron_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 electrons 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; + * electron_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. */ + double tmax; + + *b = RADIATION_LENGTH; + tmax = log(T0/ELECTRON_CRITICAL_ENERGY_D2O) - 0.5; + *a = fmax(1.1,tmax*(*b)/RADIATION_LENGTH + 1); +} + +double electron_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 electrons 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 electron in MeV and c0, c1, c2, + * and c3 are constants which I fit out. */ + return 3.141318e-1 + 2.08198e-01/log(6.33331e-03*T0 + 1.19213e+00); +} + +double electron_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 electrons 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 electron in MeV and c0, c1, c2, + * and c3 are constants which I fit out. */ + return 1.35372e-01 + 2.22344e-01/log(1.96249e-02*T0 + 1.23912e+00); +} + +double electron_get_angular_pdf_no_norm(double cos_theta, double alpha, double beta, double mu) +{ + return exp(-pow(fabs(cos_theta-mu),alpha)/beta); +} + +static double gsl_electron_get_angular_pdf(double cos_theta, void *params) +{ + double alpha = ((double *) params)[0]; + double beta = ((double *) params)[1]; + double mu = ((double *) params)[2]; + return electron_get_angular_pdf_no_norm(cos_theta,alpha,beta,mu); +} + +static double electron_get_angular_pdf_norm(double alpha, double beta, double mu) +{ + double params[3]; + gsl_integration_cquad_workspace *w; + gsl_function F; + double result, error; + int status; + size_t nevals; + + w = gsl_integration_cquad_workspace_alloc(100); + + F.function = &gsl_electron_get_angular_pdf; + F.params = params; + + params[0] = alpha; + params[1] = beta; + params[2] = mu; + + status = gsl_integration_cquad(&F, -1, 1, 0, 1e-9, w, &result, &error, &nevals); + + if (status) { + fprintf(stderr, "error integrating electron angular distribution: %s\n", gsl_strerror(status)); + exit(1); + } + + gsl_integration_cquad_workspace_free(w); + + return result; +} + +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 + * 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)); +} + static int init() { int i, j; @@ -80,6 +244,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; @@ -112,6 +277,9 @@ static int init() case 0: x[n] = value; break; + case 2: + dEdx_rad[n] = value; + break; case 3: dEdx[n] = value; break; @@ -128,6 +296,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); @@ -167,6 +339,27 @@ double electron_get_range(double T, double rho) return gsl_spline_eval(spline_range, T, acc_range)/rho; } +double electron_get_dEdx_rad(double T, double rho) +{ + /* Returns the approximate radiative dE/dx for a electron 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 electron_get_dEdx(double T, double rho) { /* Returns the approximate dE/dx for a electron in water with kinetic energy -- cgit