aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authortlatorre <tlatorre@uchicago.edu>2018-10-18 10:04:53 -0500
committertlatorre <tlatorre@uchicago.edu>2018-10-18 10:04:53 -0500
commite8efaa702e296e69d376639ebc0a70a251aaaed2 (patch)
treeaec9e00f4df29ae2f275b187e8bcbb164e043bc5 /src
parent11a22cf650448e113d4fd16c51822dc23b9c1a33 (diff)
downloadsddm-e8efaa702e296e69d376639ebc0a70a251aaaed2.tar.gz
sddm-e8efaa702e296e69d376639ebc0a70a251aaaed2.tar.bz2
sddm-e8efaa702e296e69d376639ebc0a70a251aaaed2.zip
hardcode the density when computing dE/dx
Since we only have the range and dE/dx tables for light water for electrons and protons it's not correct to use the heavy water density. Also, even though we have both tables for muons, currently we only load the heavy water table, so we hardcode the density to that of heavy water. In the future, it would be nice to load both tables and use the correct one depending on if we are fitting in the heavy or light water.
Diffstat (limited to 'src')
-rw-r--r--src/likelihood.c27
-rw-r--r--src/likelihood.h2
-rw-r--r--src/test.c6
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 = &params;
- 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);
diff --git a/src/test.c b/src/test.c
index 20a66ec..c1da272 100644
--- a/src/test.c
+++ b/src/test.c
@@ -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. */