diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/electron.c | 15 | ||||
| -rw-r--r-- | src/electron.h | 1 | ||||
| -rw-r--r-- | src/fit.c | 17 | ||||
| -rw-r--r-- | src/muon.c | 15 | ||||
| -rw-r--r-- | src/muon.h | 1 | ||||
| -rw-r--r-- | src/proton.c | 15 | ||||
| -rw-r--r-- | src/proton.h | 1 |
7 files changed, 57 insertions, 8 deletions
diff --git a/src/electron.c b/src/electron.c index 5b4133a..3f5d4d8 100644 --- a/src/electron.c +++ b/src/electron.c @@ -365,6 +365,21 @@ err: return -1; } +/* Returns the maximum kinetic energy for an electron in the range tables. + * + * If you call electron_get_range() or electron_get_dEdx() with a kinetic + * energy higher you will get a GSL interpolation error. */ +double electron_get_max_energy(void) +{ + if (!initialized) { + if (init()) { + exit(1); + } + } + + return x[size-1]; +} + double electron_get_range(double T, double rho) { /* Returns the approximate range a electron with kinetic energy `T` will travel diff --git a/src/electron.h b/src/electron.h index cc550a7..01a8915 100644 --- a/src/electron.h +++ b/src/electron.h @@ -19,6 +19,7 @@ #include <stddef.h> /* for size_t */ +double electron_get_max_energy(void); double electron_get_shower_photons(double T0, double rad); void electron_get_position_distribution_parameters(double T0, double *a, double *b); double electron_get_angular_distribution_alpha(double T0); @@ -44,16 +44,14 @@ #include "find_peaks.h" #include "util.h" #include "dc.h" +#include "electron.h" +#include "muon.h" +#include "proton.h" /* Maximum number of fit parameters. Should be at least 4 + 3*MAX_VERTICES. */ #define MAX_PARS 100 /* Maximum number of peaks to search for in Hough transform. */ #define MAX_NPEAKS 10 -/* Maximum kinetic energy for any particle. Currently set to 10^4 GeV which is - * approximately the energy of the most energetic muons in SNO (see the paper - * "Measurement of the Cosmic Ray and Neutrino-Induced Muon Flux at the Sudbury - * Neutrino Observatory".*/ -#define MAX_ENERGY 1e7 char *GitSHA1(void); char *GitDirty(void); @@ -5263,7 +5261,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double * press ctrl-c. */ size_t i, j; fitParams fpars; - double x[100], ss[100], lb[100], ub[100], fval, n_d2o, x0[100], T0, Tmin, mass, pos[3], t0, dir[3]; + double x[100], ss[100], lb[100], ub[100], fval, n_d2o, x0[100], T0, Tmin, mass, pos[3], t0, dir[3], Tmax; struct timeval tv_start, tv_stop; double time_elapsed; double peak_theta[MAX_NPEAKS]; @@ -5371,13 +5369,16 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double case IDP_E_MINUS: case IDP_E_PLUS: mass = ELECTRON_MASS; + Tmax = electron_get_max_energy(); break; case IDP_MU_MINUS: case IDP_MU_PLUS: mass = MUON_MASS; + Tmax = muon_get_max_energy(); break; case IDP_PROTON: mass = PROTON_MASS; + Tmax = proton_get_max_energy(); break; default: fprintf(stderr, "unknown particle id %i\n", id[j]); @@ -5391,7 +5392,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double * threshold. */ double Tmin2 = sqrt(mass*mass/(1-pow(0.9,2)))-mass; if (T0 < Tmin2) T0 = Tmin2; - if (T0 > MAX_ENERGY) T0 = MAX_ENERGY; + if (T0 > Tmax) T0 = Tmax; x[4+3*j] = T0; ss[4+3*j] = x[4+3*j]*0.1; @@ -5400,7 +5401,7 @@ int fit_event2(event *ev, double *xopt, double *fmin, int *id, size_t n, double lb[4+3*j] = Tmin; lb[5+3*j] = -INFINITY; lb[6+3*j] = -INFINITY; - ub[4+3*j] = MAX_ENERGY; + ub[4+3*j] = Tmax; ub[5+3*j] = +INFINITY; ub[6+3*j] = +INFINITY; } @@ -312,6 +312,21 @@ err: return -1; } +/* Returns the maximum kinetic energy for a muon in the range tables. + * + * If you call muon_get_range() or muon_get_dEdx() with a kinetic energy higher + * you will get a GSL interpolation error. */ +double muon_get_max_energy(void) +{ + if (!initialized) { + if (init()) { + exit(1); + } + } + + return x[size-1]; +} + double muon_get_range(double T, double rho) { /* Returns the approximate range a muon with kinetic energy `T` will travel @@ -30,6 +30,7 @@ * FIXME: Actually determine what this is. */ #define MUON_PHOTONS_PER_MEV 7368.0 +double muon_get_max_energy(void); double muon_get_shower_photons(double T0, double rad); void muon_get_position_distribution_parameters(double T0, double *a, double *b); double muon_get_angular_distribution_alpha(double T0); diff --git a/src/proton.c b/src/proton.c index 664d930..4714d76 100644 --- a/src/proton.c +++ b/src/proton.c @@ -165,6 +165,21 @@ err: return -1; } +/* Returns the maximum kinetic energy for a proton in the range tables. + * + * If you call proton_get_range() or proton_get_dEdx() with a kinetic + * energy higher you will get a GSL interpolation error. */ +double proton_get_max_energy(void) +{ + if (!initialized) { + if (init()) { + exit(1); + } + } + + return x[size-1]; +} + double proton_get_range(double T, double rho) { /* Returns the approximate range a proton with kinetic energy `T` will travel diff --git a/src/proton.h b/src/proton.h index f658c09..c402229 100644 --- a/src/proton.h +++ b/src/proton.h @@ -28,6 +28,7 @@ * FIXME: Actually determine what this is. */ #define PROTON_PHOTONS_PER_MEV 400.0 +double proton_get_max_energy(void); 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); |