From 9777b92d16cb3a3434c43b084d88f59321567e2d Mon Sep 17 00:00:00 2001 From: tlatorre Date: Mon, 26 Nov 2018 09:53:28 -0600 Subject: update electron range tables This commit also adds a script to calculate the CSDA range for electrons from a table of the stopping power as a function of energy. We need this script since the NIST ESTAR website will only output the stopping power above a certain energy threshold (1 GeV for electrons). See https://physics.nist.gov/PhysRefData/Star/Text/ESTAR.html. --- src/.gitignore | 1 + src/Makefile | 4 +- src/calculate-csda-range.c | 176 +++++++++++++++++++++++++++++++++++++++++++++ src/e_water_liquid.txt | 13 ++++ 4 files changed, 193 insertions(+), 1 deletion(-) create mode 100644 src/.gitignore create mode 100644 src/calculate-csda-range.c diff --git a/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..2750fd2 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1 @@ +calculate-csda-range diff --git a/src/Makefile b/src/Makefile index 62a95af..80eef49 100644 --- a/src/Makefile +++ b/src/Makefile @@ -3,7 +3,7 @@ release_hdr := $(shell sh -c './mkreleasehdr.sh') CFLAGS=-O2 -Wall -g -DSWAP_BYTES LDLIBS=-lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++ -all: test test-vector test-likelihood fit test-charge test-path +all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range Makefile.dep: -$(CC) -MM *.c > Makefile.dep @@ -28,6 +28,8 @@ test-zebra: test-zebra.o zebra.o pack2b.o fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o +calculate-csda-range: calculate-csda-range.o + clean: rm -f *.o calculate_limits test Makefile.dep diff --git a/src/calculate-csda-range.c b/src/calculate-csda-range.c new file mode 100644 index 0000000..7d71096 --- /dev/null +++ b/src/calculate-csda-range.c @@ -0,0 +1,176 @@ +#include +#include +#include /* for size_t, strtod() */ +#include /* for errno */ +#include /* for strerror(), strtok() */ +#include "sno.h" + +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; + +static gsl_interp_accel *acc_range; +static gsl_spline *spline_range; + +static int init() +{ + int i, j; + char line[256]; + char *str; + double value; + int n; + + FILE *f = fopen("e_water_liquid.txt", "r"); + + if (!f) { + fprintf(stderr, "failed to open e_water_liquid.txt: %s\n", 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 8 lines since it's just a header. */ + if (i <= 8) continue; + + if (!len) continue; + else if (line[0] == '#') continue; + + str = strtok(line," \n"); + + while (str) { + value = strtod(str, NULL); + str = strtok(NULL," \n"); + } + + n += 1; + } + + x = malloc(sizeof(double)*n); + dEdx_rad = 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 8 lines since it's just a header. */ + if (i <= 8) continue; + + if (!len) continue; + else if (line[0] == '#') continue; + + str = strtok(line," \n"); + + j = 0; + while (str) { + value = strtod(str, NULL); + switch (j) { + case 0: + x[n] = value; + break; + case 2: + dEdx_rad[n] = value; + break; + case 3: + dEdx[n] = value; + break; + case 4: + csda_range[n] = value; + break; + } + j += 1; + str = strtok(NULL," \n"); + } + + n += 1; + } + + 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); + + acc_range = gsl_interp_accel_alloc(); + spline_range = gsl_spline_alloc(gsl_interp_linear, size); + gsl_spline_init(spline_range, x, csda_range, size); + + return 0; + +err: + fclose(f); + + return -1; +} + +double electron_get_dEdx(double T, double rho) +{ + /* Returns the approximate 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 (T < spline_dEdx->x[0]) return spline_dEdx->y[0]; + + return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho; +} + +int main(int argc, char **argv) +{ + int i; + double T0, T, dx, range; + + init(); + + dx = 1e-5; + + for (i = 0; i < size; i++) { + T0 = x[i]; + + T = T0; + + range = 0.0; + while (T > 0) { + double dEdx2 = electron_get_dEdx(T, WATER_DENSITY); + T -= dEdx2*dx; + range += dx; + } + + printf("%.3E %.3E\n", T0, range); + } + + return 0; +} diff --git a/src/e_water_liquid.txt b/src/e_water_liquid.txt index 4e29104..a4e3a8f 100644 --- a/src/e_water_liquid.txt +++ b/src/e_water_liquid.txt @@ -87,3 +87,16 @@ MeV MeV cm2/g MeV cm2/g MeV cm2/g g/cm2 8.000E+02 2.381E+00 2.133E+01 2.371E+01 9.258E+01 7.425E-01 1.121E+01 9.000E+02 2.391E+00 2.406E+01 2.645E+01 9.657E+01 7.605E-01 1.145E+01 1.000E+03 2.400E+00 2.679E+01 2.919E+01 1.002E+02 7.759E-01 1.166E+01 +1.500E+03 2.435E+00 4.047E+01 4.291E+01 1.142E+02 - 1.247E+01 +2.000E+03 2.459E+00 5.418E+01 5.664E+01 1.244E+02 - 1.304E+01 +2.500E+03 2.479E+00 6.790E+01 7.038E+01 1.323E+02 - 1.349E+01 +3.000E+03 2.494E+00 8.163E+01 8.413E+01 1.388E+02 - 1.385E+01 +3.500E+03 2.507E+00 9.537E+01 9.788E+01 1.443E+02 - 1.416E+01 +4.000E+03 2.519E+00 1.091E+02 1.116E+02 1.490E+02 - 1.443E+01 +4.500E+03 2.529E+00 1.229E+02 1.254E+02 1.533E+02 - 1.466E+01 +5.000E+03 2.538E+00 1.366E+02 1.391E+02 1.570E+02 - 1.488E+01 +6.000E+03 2.553E+00 1.641E+02 1.667E+02 1.636E+02 - 1.524E+01 +7.000E+03 2.566E+00 1.916E+02 1.942E+02 1.692E+02 - 1.555E+01 +8.000E+03 2.578E+00 2.192E+02 2.217E+02 1.740E+02 - 1.582E+01 +9.000E+03 2.588E+00 2.467E+02 2.493E+02 1.782E+02 - 1.605E+01 +1.000E+04 2.597E+00 2.742E+02 2.768E+02 1.820E+02 - 1.626E+01 -- cgit