diff options
| -rw-r--r-- | src/likelihood.c | 27 | ||||
| -rw-r--r-- | src/likelihood.h | 2 | ||||
| -rw-r--r-- | src/test.c | 6 |
3 files changed, 22 insertions, 13 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index f426271..1bf9549 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -27,14 +27,14 @@ typedef struct intParams { int i; } intParams; -particle *particle_init(int id, double T0, double rho, size_t n) +particle *particle_init(int id, double T0, size_t n) { /* Returns a struct with arrays of the particle position and kinetic * energy. This struct can then be passed to particle_get_energy() to * interpolate the particle's kinetic energy at any point along the track. * For example: * - * particle *p = particle_init(IDP_MU_MINUS, 1000.0, HEAVY_WATER_DENSITY, 100); + * particle *p = particle_init(IDP_MU_MINUS, 1000.0, 100); * double T = particle_get_energy(x, p); * * To compute the kinetic energy as a function of distance we need to solve @@ -59,13 +59,15 @@ particle *particle_init(int id, double T0, double rho, size_t n) case IDP_E_MINUS: case IDP_E_PLUS: p->mass = ELECTRON_MASS; - p->range = electron_get_range(T0, rho); + /* We use the density of light water here since we don't have the + * tables for heavy water for electrons. */ + p->range = electron_get_range(T0, WATER_DENSITY); /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ for (i = 1; i < n; i++) { p->x[i] = p->range*i/(n-1); - dEdx = electron_get_dEdx(p->T[i-1], rho); + 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]); } @@ -73,26 +75,33 @@ particle *particle_init(int id, double T0, double rho, size_t n) case IDP_MU_MINUS: case IDP_MU_PLUS: p->mass = MUON_MASS; - p->range = muon_get_range(T0, rho); + /* We use the density of heavy water here since we only load the heavy + * water table for muons. + * + * FIXME: It would be nice in the future to load both tables and then + * load the correct table based on the media. */ + p->range = muon_get_range(T0, HEAVY_WATER_DENSITY); /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ for (i = 1; i < n; i++) { p->x[i] = p->range*i/(n-1); - dEdx = muon_get_dEdx(p->T[i-1], rho); + 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]); } break; case IDP_PROTON: p->mass = PROTON_MASS; - p->range = proton_get_range(T0, rho); + /* We use the density of light water here since we don't have the + * tables for heavy water for protons. */ + p->range = proton_get_range(T0, WATER_DENSITY); /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ for (i = 1; i < n; i++) { p->x[i] = p->range*i/(n-1); - dEdx = proton_get_dEdx(p->T[i-1], rho); + 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]); } @@ -521,7 +530,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t gsl_function F; F.params = ¶ms; - p = particle_init(id, T0, HEAVY_WATER_DENSITY, 10000); + p = particle_init(id, T0, 10000); range = p->range; /* Calculate total energy */ diff --git a/src/likelihood.h b/src/likelihood.h index 945fc46..71f77f0 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -41,7 +41,7 @@ typedef struct particle { size_t n; } particle; -particle *particle_init(int id, double T0, double rho, size_t n); +particle *particle_init(int id, double T0, size_t n); double particle_get_energy(double x, particle *p); void particle_free(particle *p); double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected); @@ -124,7 +124,7 @@ int test_proton_get_energy(char *err) /* Assume initial kinetic energy is 1 GeV. */ T0 = 1000.0; - p = particle_init(IDP_PROTON, T0,1.0,10000); + p = particle_init(IDP_PROTON, T0, 10000); T = particle_get_energy(1e-9,p); /* At the beginning of the track we should have roughly the same energy. */ @@ -160,7 +160,7 @@ int test_electron_get_energy(char *err) /* Assume initial kinetic energy is 1 GeV. */ T0 = 1000.0; - p = particle_init(IDP_E_MINUS, T0,1.0,10000); + p = particle_init(IDP_E_MINUS, T0, 10000); T = particle_get_energy(1e-9,p); /* At the beginning of the track we should have roughly the same energy. */ @@ -196,7 +196,7 @@ int test_muon_get_energy(char *err) /* Assume initial kinetic energy is 1 GeV. */ T0 = 1000.0; - p = particle_init(IDP_MU_MINUS, T0,1.0,10000); + p = particle_init(IDP_MU_MINUS, T0, 10000); T = particle_get_energy(1e-9,p); /* At the beginning of the track we should have roughly the same energy. */ |