aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Makefile4
-rw-r--r--src/e_water_liquid.txt89
-rw-r--r--src/electron.c189
-rw-r--r--src/electron.h9
-rw-r--r--src/fit.c3
-rw-r--r--src/likelihood.c239
-rw-r--r--src/likelihood.h14
-rw-r--r--src/muon.c169
-rw-r--r--src/muon.h15
-rw-r--r--src/path.c25
-rw-r--r--src/path.h4
-rw-r--r--src/pdg.c3
-rw-r--r--src/pdg.h2
-rw-r--r--src/proton.c180
-rw-r--r--src/proton.h9
-rw-r--r--src/proton_water_liquid.txt140
-rw-r--r--src/test-path.c8
-rw-r--r--src/test.c174
18 files changed, 1013 insertions, 263 deletions
diff --git a/src/Makefile b/src/Makefile
index 377388c..62a95af 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -14,7 +14,7 @@ calculate_limits: calculate_limits.c
solid_angle.o: solid_angle.c
-test: test.o solid_angle.o optics.o muon.o vector.o quantum_efficiency.o pdg.o scattering.o misc.o mt19937ar.o sno_charge.o path.o random.o pmt_response.o db.o dict.o siphash.o
+test: test.o solid_angle.o optics.o muon.o vector.o quantum_efficiency.o pdg.o scattering.o misc.o mt19937ar.o sno_charge.o path.o random.o pmt_response.o db.o dict.o siphash.o electron.o proton.o likelihood.o pmt.o
test-vector: test-vector.o vector.o mt19937ar.o
@@ -26,7 +26,7 @@ test-charge: test-charge.o sno_charge.o misc.o vector.o
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
+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
clean:
rm -f *.o calculate_limits test Makefile.dep
diff --git a/src/e_water_liquid.txt b/src/e_water_liquid.txt
new file mode 100644
index 0000000..4e29104
--- /dev/null
+++ b/src/e_water_liquid.txt
@@ -0,0 +1,89 @@
+ESTAR: Stopping Powers and Range Tables for Electrons
+
+WATER, LIQUID
+
+Kinetic Collision Radiative Total CSDA Radiation D. Effect
+Energy Stp. Pow. Stp. Pow. Stp. Pow. Range Yield Parameter
+MeV MeV cm2/g MeV cm2/g MeV cm2/g g/cm2
+
+1.000E-02 2.256E+01 3.898E-03 2.256E+01 2.515E-04 9.408E-05 0.000E+00
+1.250E-02 1.897E+01 3.927E-03 1.898E+01 3.728E-04 1.133E-04 0.000E+00
+1.500E-02 1.647E+01 3.944E-03 1.647E+01 5.147E-04 1.316E-04 0.000E+00
+1.750E-02 1.461E+01 3.955E-03 1.461E+01 6.762E-04 1.493E-04 0.000E+00
+2.000E-02 1.317E+01 3.963E-03 1.318E+01 8.566E-04 1.663E-04 0.000E+00
+2.500E-02 1.109E+01 3.974E-03 1.110E+01 1.272E-03 1.990E-04 0.000E+00
+3.000E-02 9.653E+00 3.984E-03 9.657E+00 1.756E-03 2.301E-04 0.000E+00
+3.500E-02 8.592E+00 3.994E-03 8.596E+00 2.306E-03 2.599E-04 0.000E+00
+4.000E-02 7.777E+00 4.005E-03 7.781E+00 2.919E-03 2.886E-04 0.000E+00
+4.500E-02 7.130E+00 4.018E-03 7.134E+00 3.591E-03 3.165E-04 0.000E+00
+5.000E-02 6.603E+00 4.031E-03 6.607E+00 4.320E-03 3.435E-04 0.000E+00
+5.500E-02 6.166E+00 4.046E-03 6.170E+00 5.103E-03 3.698E-04 0.000E+00
+6.000E-02 5.797E+00 4.062E-03 5.801E+00 5.940E-03 3.955E-04 0.000E+00
+7.000E-02 5.207E+00 4.098E-03 5.211E+00 7.762E-03 4.453E-04 0.000E+00
+8.000E-02 4.757E+00 4.138E-03 4.761E+00 9.773E-03 4.931E-04 0.000E+00
+9.000E-02 4.402E+00 4.181E-03 4.407E+00 1.196E-02 5.393E-04 0.000E+00
+1.000E-01 4.115E+00 4.228E-03 4.119E+00 1.431E-02 5.842E-04 0.000E+00
+1.250E-01 3.591E+00 4.355E-03 3.596E+00 2.083E-02 6.912E-04 0.000E+00
+1.500E-01 3.238E+00 4.494E-03 3.242E+00 2.817E-02 7.926E-04 0.000E+00
+1.750E-01 2.984E+00 4.643E-03 2.988E+00 3.622E-02 8.894E-04 0.000E+00
+2.000E-01 2.793E+00 4.801E-03 2.798E+00 4.488E-02 9.826E-04 0.000E+00
+2.500E-01 2.528E+00 5.141E-03 2.533E+00 6.372E-02 1.161E-03 0.000E+00
+3.000E-01 2.355E+00 5.514E-03 2.360E+00 8.421E-02 1.331E-03 0.000E+00
+3.500E-01 2.235E+00 5.914E-03 2.241E+00 1.060E-01 1.496E-03 0.000E+00
+4.000E-01 2.148E+00 6.339E-03 2.154E+00 1.288E-01 1.658E-03 0.000E+00
+4.500E-01 2.083E+00 6.787E-03 2.090E+00 1.523E-01 1.818E-03 0.000E+00
+5.000E-01 2.034E+00 7.257E-03 2.041E+00 1.766E-01 1.976E-03 0.000E+00
+5.500E-01 1.995E+00 7.747E-03 2.003E+00 2.013E-01 2.134E-03 1.103E-02
+6.000E-01 1.963E+00 8.254E-03 1.972E+00 2.265E-01 2.292E-03 2.938E-02
+7.000E-01 1.917E+00 9.313E-03 1.926E+00 2.778E-01 2.608E-03 7.435E-02
+8.000E-01 1.886E+00 1.042E-02 1.896E+00 3.302E-01 2.928E-03 1.267E-01
+9.000E-01 1.864E+00 1.159E-02 1.876E+00 3.832E-01 3.251E-03 1.835E-01
+1.000E+00 1.849E+00 1.280E-02 1.862E+00 4.367E-01 3.579E-03 2.428E-01
+1.250E+00 1.829E+00 1.600E-02 1.845E+00 5.717E-01 4.416E-03 3.944E-01
+1.500E+00 1.822E+00 1.942E-02 1.841E+00 7.075E-01 5.281E-03 5.437E-01
+1.750E+00 1.821E+00 2.303E-02 1.844E+00 8.432E-01 6.171E-03 6.866E-01
+2.000E+00 1.824E+00 2.678E-02 1.850E+00 9.785E-01 7.085E-03 8.218E-01
+2.500E+00 1.834E+00 3.468E-02 1.868E+00 1.247E+00 8.969E-03 1.069E+00
+3.000E+00 1.846E+00 4.299E-02 1.889E+00 1.514E+00 1.092E-02 1.288E+00
+3.500E+00 1.858E+00 5.164E-02 1.910E+00 1.777E+00 1.291E-02 1.484E+00
+4.000E+00 1.870E+00 6.058E-02 1.931E+00 2.037E+00 1.495E-02 1.660E+00
+4.500E+00 1.882E+00 6.976E-02 1.951E+00 2.295E+00 1.702E-02 1.821E+00
+5.000E+00 1.892E+00 7.917E-02 1.971E+00 2.550E+00 1.911E-02 1.967E+00
+5.500E+00 1.902E+00 8.876E-02 1.991E+00 2.802E+00 2.123E-02 2.102E+00
+6.000E+00 1.911E+00 9.854E-02 2.010E+00 3.052E+00 2.336E-02 2.227E+00
+7.000E+00 1.928E+00 1.185E-01 2.047E+00 3.545E+00 2.766E-02 2.453E+00
+8.000E+00 1.943E+00 1.391E-01 2.082E+00 4.030E+00 3.200E-02 2.652E+00
+9.000E+00 1.956E+00 1.601E-01 2.116E+00 4.506E+00 3.636E-02 2.831E+00
+1.000E+01 1.968E+00 1.814E-01 2.149E+00 4.975E+00 4.072E-02 2.992E+00
+1.250E+01 1.993E+00 2.362E-01 2.230E+00 6.117E+00 5.163E-02 3.341E+00
+1.500E+01 2.014E+00 2.926E-01 2.306E+00 7.219E+00 6.243E-02 3.633E+00
+1.750E+01 2.031E+00 3.501E-01 2.381E+00 8.286E+00 7.309E-02 3.885E+00
+2.000E+01 2.046E+00 4.086E-01 2.454E+00 9.320E+00 8.355E-02 4.107E+00
+2.500E+01 2.070E+00 5.277E-01 2.598E+00 1.130E+01 1.039E-01 4.487E+00
+3.000E+01 2.089E+00 6.489E-01 2.738E+00 1.317E+01 1.233E-01 4.806E+00
+3.500E+01 2.105E+00 7.716E-01 2.876E+00 1.496E+01 1.418E-01 5.082E+00
+4.000E+01 2.118E+00 8.955E-01 3.013E+00 1.665E+01 1.594E-01 5.326E+00
+4.500E+01 2.129E+00 1.021E+00 3.150E+00 1.828E+01 1.762E-01 5.544E+00
+5.000E+01 2.139E+00 1.146E+00 3.286E+00 1.983E+01 1.923E-01 5.741E+00
+5.500E+01 2.148E+00 1.273E+00 3.421E+00 2.132E+01 2.076E-01 5.921E+00
+6.000E+01 2.156E+00 1.400E+00 3.556E+00 2.276E+01 2.222E-01 6.087E+00
+7.000E+01 2.170E+00 1.656E+00 3.827E+00 2.547E+01 2.496E-01 6.383E+00
+8.000E+01 2.182E+00 1.914E+00 4.096E+00 2.799E+01 2.747E-01 6.641E+00
+9.000E+01 2.193E+00 2.173E+00 4.366E+00 3.035E+01 2.978E-01 6.871E+00
+1.000E+02 2.202E+00 2.434E+00 4.636E+00 3.258E+01 3.192E-01 7.077E+00
+1.250E+02 2.222E+00 3.089E+00 5.311E+00 3.761E+01 3.662E-01 7.516E+00
+1.500E+02 2.238E+00 3.749E+00 5.987E+00 4.204E+01 4.060E-01 7.876E+00
+1.750E+02 2.251E+00 4.412E+00 6.663E+00 4.600E+01 4.401E-01 8.182E+00
+2.000E+02 2.263E+00 5.078E+00 7.341E+00 4.957E+01 4.698E-01 8.447E+00
+2.500E+02 2.282E+00 6.416E+00 8.698E+00 5.582E+01 5.190E-01 8.891E+00
+3.000E+02 2.297E+00 7.760E+00 1.006E+01 6.116E+01 5.584E-01 9.254E+00
+3.500E+02 2.311E+00 9.107E+00 1.142E+01 6.583E+01 5.908E-01 9.561E+00
+4.000E+02 2.322E+00 1.046E+01 1.278E+01 6.996E+01 6.180E-01 9.827E+00
+4.500E+02 2.332E+00 1.181E+01 1.414E+01 7.368E+01 6.412E-01 1.006E+01
+5.000E+02 2.341E+00 1.317E+01 1.551E+01 7.706E+01 6.613E-01 1.027E+01
+5.500E+02 2.349E+00 1.453E+01 1.688E+01 8.014E+01 6.789E-01 1.046E+01
+6.000E+02 2.357E+00 1.589E+01 1.824E+01 8.299E+01 6.945E-01 1.064E+01
+7.000E+02 2.370E+00 1.861E+01 2.098E+01 8.810E+01 7.209E-01 1.094E+01
+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
diff --git a/src/electron.c b/src/electron.c
new file mode 100644
index 0000000..e8e24d3
--- /dev/null
+++ b/src/electron.c
@@ -0,0 +1,189 @@
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_spline.h>
+#include <math.h>
+#include "optics.h"
+#include "quantum_efficiency.h"
+#include "solid_angle.h"
+#include "pdg.h"
+#include "vector.h"
+#include "electron.h"
+#include "sno.h"
+#include "scattering.h"
+#include "pmt_response.h"
+#include "misc.h"
+
+static int initialized = 0;
+
+static double *x, *dEdx, *csda_range;
+static size_t size;
+
+static gsl_interp_accel *acc_dEdx;
+static gsl_spline *spline_dEdx;
+
+static gsl_interp_accel *acc_range;
+static gsl_spline *spline_range;
+
+/* Electron critical energy in H2O and D2O.
+ *
+ * These values come from
+ * http://pdgprod.lbl.gov/~deg/AtomicNuclearProperties/HTML/deuterium_oxide_liquid.html for D2O and
+ * http://pdg.lbl.gov/2018/AtomicNuclearProperties/HTML/water_liquid.html for H2O.
+ */
+static const double ELECTRON_CRITICAL_ENERGY_H2O = 78.33;
+static const double ELECTRON_CRITICAL_ENERGY_D2O = 78.33;
+
+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 = 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 3:
+ dEdx[n] = value;
+ break;
+ case 4:
+ csda_range[n] = value;
+ break;
+ }
+ j += 1;
+ str = strtok(NULL," \n");
+ }
+
+ n += 1;
+ }
+
+ fclose(f);
+
+ 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);
+
+ initialized = 1;
+
+ return 0;
+
+err:
+ fclose(f);
+
+ return -1;
+}
+
+double electron_get_range(double T, double rho)
+{
+ /* Returns the approximate range a electron with kinetic energy `T` will travel
+ * in water before losing all of its energy. This range is interpolated
+ * based on data from the PDG which uses the continuous slowing down
+ * approximation.
+ *
+ * `T` should be in MeV, and `rho` should be in g/cm^3.
+ *
+ * Return value is in cm.
+ *
+ * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+ if (!initialized) {
+ if (init()) {
+ exit(1);
+ }
+ }
+
+ return gsl_spline_eval(spline_range, T, acc_range)/rho;
+}
+
+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 (!initialized) {
+ if (init()) {
+ exit(1);
+ }
+ }
+
+ if (T < spline_dEdx->x[0]) return spline_dEdx->y[0];
+
+ return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
+}
diff --git a/src/electron.h b/src/electron.h
new file mode 100644
index 0000000..1bb589f
--- /dev/null
+++ b/src/electron.h
@@ -0,0 +1,9 @@
+#ifndef ELECTRON_H
+#define ELECTRON_H
+
+#include <stddef.h> /* for size_t */
+
+double electron_get_range(double T, double rho);
+double electron_get_dEdx(double T, double rho);
+
+#endif
diff --git a/src/fit.c b/src/fit.c
index 3ab0ffe..dc4f88b 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -21,6 +21,7 @@
#include "pmt_response.h"
#include <signal.h> /* for signal() */
#include "release.h"
+#include "id_particles.h"
char *GitSHA1(void);
char *GitDirty(void);
@@ -4995,7 +4996,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
z2[0] = x[8];
gettimeofday(&tv_start, NULL);
- fval = nll_muon(fpars->ev, T, pos, dir, t0, z1, z2, 1, fpars->epsrel, fpars->fast);
+ fval = nll_muon(fpars->ev, IDP_MU_MINUS, T, pos, dir, t0, z1, z2, 1, fpars->epsrel, fpars->fast);
gettimeofday(&tv_stop, NULL);
long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000;
diff --git a/src/likelihood.c b/src/likelihood.c
index 5f66a72..7da81bb 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -18,12 +18,164 @@
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
#include "pmt_response.h"
+#include "electron.h"
+#include "proton.h"
+#include "id_particles.h"
typedef struct intParams {
path *p;
int i;
} intParams;
+particle *particle_init(int id, double T0, double rho, 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);
+ * double T = particle_get_energy(x, p);
+ *
+ * To compute the kinetic energy as a function of distance we need to solve
+ * the following integral equation:
+ *
+ * T(x) = T0 - \int_0^x dT(T(x))/dx
+ *
+ * which we solve by dividing the track up into `n` segments and then
+ * numerically computing the energy at each step. */
+ size_t i;
+ double dEdx;
+
+ particle *p = malloc(sizeof(particle));
+ p->x = malloc(sizeof(double)*n);
+ p->T = malloc(sizeof(double)*n);
+ p->n = n;
+
+ p->x[0] = 0;
+ p->T[0] = T0;
+
+ switch (id) {
+ case IDP_E_MINUS:
+ case IDP_E_PLUS:
+ p->mass = ELECTRON_MASS;
+ p->range = electron_get_range(T0, rho);
+
+ /* 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);
+ p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ }
+
+ break;
+ case IDP_MU_MINUS:
+ case IDP_MU_PLUS:
+ p->mass = MUON_MASS;
+ p->range = muon_get_range(T0, rho);
+
+ /* 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);
+ 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);
+
+ /* 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);
+ p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ }
+
+ break;
+ default:
+ fprintf(stderr, "unknown particle id %i\n", id);
+ exit(1);
+ }
+
+ return p;
+}
+
+double particle_get_energy(double x, particle *p)
+{
+ /* Returns the approximate kinetic energy of a particle in water after
+ * travelling `x` cm with an initial kinetic energy `T`.
+ *
+ * Return value is in MeV. */
+ double T;
+
+ T = interp1d(x,p->x,p->T,p->n);
+
+ if (T < 0) return 0;
+
+ return T;
+}
+
+void particle_free(particle *p)
+{
+ free(p->x);
+ free(p->T);
+ free(p);
+}
+
+double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected)
+{
+ double pmt_dir[3], cos_theta, n, wavelength0, omega, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
+
+ z = 1.0;
+
+ SUB(pmt_dir,pmt_pos,pos);
+
+ normalize(pmt_dir);
+
+ cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+
+ if (cos_theta_pmt > 0) return 0;
+
+ /* Calculate the cosine of the angle between the track direction and the
+ * vector to the PMT. */
+ cos_theta = DOT(dir,pmt_dir);
+
+ omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
+
+ R = NORM(pos);
+
+ /* FIXME: I just calculate delta assuming 400 nm light. */
+ wavelength0 = 400.0;
+
+ if (R <= AV_RADIUS) {
+ n = get_index_snoman_d2o(wavelength0);
+ } else {
+ n = get_index_snoman_h2o(wavelength0);
+ }
+
+ if (beta < 1/n) return 0;
+
+ if (reflected)
+ f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt));
+ else
+ f = get_weighted_pmt_response(acos(-cos_theta_pmt));
+
+ absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
+ absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
+ absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
+
+ get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+ l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
+
+ return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+}
+
double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
{
/* Returns the CDF for the time distribution of photons at time `t`. */
@@ -86,12 +238,12 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma));
}
-static double gsl_muon_time(double x, void *params)
+static double gsl_time(double x, void *params)
{
intParams *pars = (intParams *) params;
- double dir[3], pos[3], n_d2o, n_h2o, wavelength0, T, t, theta0, l_d2o, l_h2o;
+ double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o;
- path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+ path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
get_path_length(pos,pmts[pars->i].pos,AV_RADIUS,&l_d2o,&l_h2o);
@@ -102,30 +254,30 @@ static double gsl_muon_time(double x, void *params)
t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
- return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+ return t*get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
}
-static double gsl_muon_reflected_charge(double x, void *params)
+static double gsl_reflected_charge(double x, void *params)
{
intParams *pars = (intParams *) params;
- double dir[3], pos[3], T, t, theta0;
+ double dir[3], pos[3], t, theta0, beta;
- path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+ path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
- return get_expected_reflected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+ return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 1);
}
-static double gsl_muon_charge(double x, void *params)
+static double gsl_charge(double x, void *params)
{
intParams *pars = (intParams *) params;
- double dir[3], pos[3], T, t, theta0;
+ double dir[3], pos[3], t, theta0, beta;
- path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+ path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
- return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+ return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
}
-double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t, double *mu_reflected)
+double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected)
{
/* Returns the approximate expected number of photons seen by PMT `i` using
* an analytic formula.
@@ -156,11 +308,11 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
*
* `smax` is currently calculated as the point where the particle velocity
* drops to 0.8 times the speed of light. */
- double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected;
+ double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected;
/* Calculate beta at the start of the track. */
- E0 = T0 + MUON_MASS;
- p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+ E0 = T0 + p->mass;
+ p0 = sqrt(E0*E0 - p->mass*p->mass);
beta0 = p0/E0;
/* First, we find the point along the track where the PMT is at the
@@ -219,12 +371,12 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
b = (s-z)/a;
/* Compute the kinetic energy at the point `s` along the track. */
- T = muon_get_energy(s,m);
+ T = particle_get_energy(s,p);
/* Calculate the particle velocity at the point `s`. */
- E = T + MUON_MASS;
- p = sqrt(E*E - MUON_MASS*MUON_MASS);
- beta = p/E;
+ E = T + p->mass;
+ mom = sqrt(E*E - p->mass*p->mass);
+ beta = mom/E;
if (beta < 1/n_d2o) return 0.0;
@@ -275,7 +427,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
}
typedef struct betaRootParams {
- muon_energy *m;
+ particle *p;
double beta_min;
} betaRootParams;
@@ -283,22 +435,22 @@ static double beta_root(double x, void *params)
{
/* Function used to find at what point along a track a particle crosses
* some threshold in beta. */
- double T, E, p, beta;
+ double T, E, mom, beta;
betaRootParams *pars;
pars = (betaRootParams *) params;
- T = muon_get_energy(x, pars->m);
+ T = particle_get_energy(x, pars->p);
/* Calculate total energy */
- E = T + MUON_MASS;
- p = sqrt(E*E - MUON_MASS*MUON_MASS);
- beta = p/E;
+ E = T + pars->p->mass;
+ mom = sqrt(E*E - pars->p->mass*pars->p->mass);
+ beta = mom/E;
return beta - pars->beta_min;
}
-static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
+static int get_smax(particle *p, double beta_min, double range, double *smax)
{
/* Find the point along the track at which the particle's velocity drops to
* `beta_min`. */
@@ -311,7 +463,7 @@ static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
- pars.m = m;
+ pars.p = p;
pars.beta_min = beta_min;
F.function = &beta_root;
@@ -339,17 +491,17 @@ static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
return status;
}
-double getKineticEnergy(double x, void *params)
+static double getKineticEnergy(double x, void *p)
{
- return muon_get_energy(x,(muon_energy *) params);
+ return particle_get_energy(x, (particle *) p);
}
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
+double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
{
size_t i, j, nhit;
intParams params;
double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
- muon_energy *m;
+ particle *p;
double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o;
double logp_path;
double a, b;
@@ -369,22 +521,21 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
gsl_function F;
F.params = &params;
- range = get_range(T0, HEAVY_WATER_DENSITY);
+ p = particle_init(id, T0, HEAVY_WATER_DENSITY, 10000);
+ range = p->range;
/* Calculate total energy */
- E0 = T0 + MUON_MASS;
- p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+ E0 = T0 + p->mass;
+ p0 = sqrt(E0*E0 - p->mass*p->mass);
beta0 = p0/E0;
/* FIXME: is this formula valid for muons? */
theta0 = get_scattering_rms(range/2,p0,beta0,1.0,HEAVY_WATER_DENSITY)/sqrt(range/2);
- m = muon_init_energy(T0,HEAVY_WATER_DENSITY,10000);
-
- params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, m, z1, z2, n, MUON_MASS);
+ params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
if (beta0 > BETA_MIN)
- get_smax(m, BETA_MIN, range, &smax);
+ get_smax(p, BETA_MIN, range, &smax);
else
smax = 0.0;
@@ -399,7 +550,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
params.i = i;
if (fast) {
- mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i], &mu_reflected);
+ mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected);
ts[i] += t0;
mu_indirect += mu_reflected;
} else {
@@ -456,17 +607,17 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
if (a < 0.0) a = 0.0;
if (b > range) b = range;
- F.function = &gsl_muon_charge;
+ F.function = &gsl_charge;
gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
mu_direct[i] = result;
- F.function = &gsl_muon_reflected_charge;
+ F.function = &gsl_reflected_charge;
gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
mu_indirect += result;
- F.function = &gsl_muon_time;
+ F.function = &gsl_time;
if (mu_direct[i] > 1e-9) {
gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
@@ -496,7 +647,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
path_free(params.p);
- muon_free_energy(m);
+ particle_free(p);
gsl_integration_cquad_workspace_free(w);
diff --git a/src/likelihood.h b/src/likelihood.h
index d708f19..945fc46 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -33,6 +33,18 @@
#define GTVALID 400.0
#define BETA_MIN 0.8
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast);
+typedef struct particle {
+ double mass;
+ double range;
+ double *x;
+ double *T;
+ size_t n;
+} particle;
+
+particle *particle_init(int id, double T0, double rho, 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);
+double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast);
#endif
diff --git a/src/muon.c b/src/muon.c
index f0ede57..02cf55a 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -151,7 +151,7 @@ err:
return -1;
}
-double get_range(double T, double rho)
+double muon_get_range(double T, double rho)
{
/* Returns the approximate range a muon with kinetic energy `T` will travel
* in water before losing all of its energy. This range is interpolated
@@ -172,70 +172,7 @@ double get_range(double T, double rho)
return gsl_spline_eval(spline_range, T, acc_range)/rho;
}
-muon_energy *muon_init_energy(double T0, double rho, size_t n)
-{
- /* Returns a struct with arrays of the muon position and kinetic energy.
- * This struct can then be passed to muon_get_energy() to interpolate the
- * muon's kinetic energy at any point along the track. For example:
- *
- * muon_energy *m = muon_init_energy(1000.0, HEAVY_WATER_DENSITY, 100);
- * double T = muon_get_energy(x, m);
- *
- * To compute the kinetic energy as a function of distance we need to solve
- * the following integral equation:
- *
- * T(x) = T0 - \int_0^x dT(T(x))/dx
- *
- * which we solve by dividing the track up into `n` segments and then
- * numerically computing the energy at each step. */
- size_t i;
- double range, dEdx;
-
- /* Get the range of the muon. */
- range = get_range(T0, rho);
-
- muon_energy *m = malloc(sizeof(muon_energy));
- m->x = malloc(sizeof(double)*n);
- m->T = malloc(sizeof(double)*n);
- m->n = n;
-
- m->x[0] = 0;
- m->T[0] = T0;
-
- /* Now we loop over the points along the track and assume dE/dx is constant
- * between points. */
- for (i = 1; i < n; i++) {
- m->x[i] = range*i/(n-1);
- dEdx = get_dEdx(m->T[i-1], rho);
- m->T[i] = m->T[i-1] - dEdx*(m->x[i]-m->x[i-1]);
- }
-
- return m;
-}
-
-double muon_get_energy(double x, muon_energy *m)
-{
- /* Returns the approximate kinetic energy of a muon in water after
- * travelling `x` cm with an initial kinetic energy `T`.
- *
- * Return value is in MeV. */
- double T;
-
- T = interp1d(x,m->x,m->T,m->n);
-
- if (T < 0) return 0;
-
- return T;
-}
-
-void muon_free_energy(muon_energy *m)
-{
- free(m->x);
- free(m->T);
- free(m);
-}
-
-double get_dEdx(double T, double rho)
+double muon_get_dEdx(double T, double rho)
{
/* Returns the approximate dE/dx for a muon in water with kinetic energy
* `T`.
@@ -255,105 +192,3 @@ double get_dEdx(double T, double rho)
return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
}
-
-double get_expected_reflected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
-{
- double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
-
- z = 1.0;
-
- SUB(pmt_dir,pmt_pos,pos);
-
- normalize(pmt_dir);
-
- cos_theta_pmt = DOT(pmt_dir,pmt_normal);
-
- if (cos_theta_pmt > 0) return 0;
-
- /* Calculate the cosine of the angle between the track direction and the
- * vector to the PMT. */
- cos_theta = DOT(dir,pmt_dir);
-
- /* Calculate total energy */
- E = T + MUON_MASS;
- p = sqrt(E*E - MUON_MASS*MUON_MASS);
- beta = p/E;
-
- omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
-
- R = NORM(pos);
-
- /* FIXME: I just calculate delta assuming 400 nm light. */
- wavelength0 = 400.0;
-
- if (R <= AV_RADIUS) {
- n = get_index_snoman_d2o(wavelength0);
- } else {
- n = get_index_snoman_h2o(wavelength0);
- }
-
- if (beta < 1/n) return 0;
-
- f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt));
-
- absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
- absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
- absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
-
- get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
-
- l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
-
- return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
-}
-
-double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
-{
- double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
-
- z = 1.0;
-
- SUB(pmt_dir,pmt_pos,pos);
-
- normalize(pmt_dir);
-
- cos_theta_pmt = DOT(pmt_dir,pmt_normal);
-
- if (cos_theta_pmt > 0) return 0;
-
- /* Calculate the cosine of the angle between the track direction and the
- * vector to the PMT. */
- cos_theta = DOT(dir,pmt_dir);
-
- /* Calculate total energy */
- E = T + MUON_MASS;
- p = sqrt(E*E - MUON_MASS*MUON_MASS);
- beta = p/E;
-
- omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r);
-
- R = NORM(pos);
-
- /* FIXME: I just calculate delta assuming 400 nm light. */
- wavelength0 = 400.0;
-
- if (R <= AV_RADIUS) {
- n = get_index_snoman_d2o(wavelength0);
- } else {
- n = get_index_snoman_h2o(wavelength0);
- }
-
- if (beta < 1/n) return 0;
-
- f = get_weighted_pmt_response(acos(-cos_theta_pmt));
-
- absorption_length_d2o = get_weighted_absorption_length_snoman_d2o();
- absorption_length_h2o = get_weighted_absorption_length_snoman_h2o();
- absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic();
-
- get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
-
- l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
-
- return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
-}
diff --git a/src/muon.h b/src/muon.h
index ce3db29..7717501 100644
--- a/src/muon.h
+++ b/src/muon.h
@@ -5,18 +5,7 @@
#define EULER_CONSTANT 0.57721
-typedef struct muon_energy {
- double *x;
- double *T;
- size_t n;
-} muon_energy;
-
-muon_energy *muon_init_energy(double T0, double rho, size_t n);
-double muon_get_energy(double x, muon_energy *m);
-void muon_free_energy(muon_energy *m);
-double get_range(double T, double rho);
-double get_dEdx(double T, double rho);
-double get_expected_reflected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r);
-double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r);
+double muon_get_range(double T, double rho);
+double muon_get_dEdx(double T, double rho);
#endif
diff --git a/src/path.c b/src/path.c
index 0aeb20b..340f8c7 100644
--- a/src/path.c
+++ b/src/path.c
@@ -35,7 +35,7 @@ double path_get_coefficient(unsigned int k, double *s, double *x, double theta0,
path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m)
{
size_t i, j;
- double E, mom, beta, theta, phi;
+ double T, E, mom, theta, phi;
double normal[3], k[3], tmp[3], tmp2[3];
path *p = malloc(sizeof(path));
@@ -56,7 +56,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
double *x = calloc(N,sizeof(double));
double *y = calloc(N,sizeof(double));
double *z = calloc(N,sizeof(double));
- double *T = calloc(N,sizeof(double));
+ double *beta = calloc(N,sizeof(double));
double *t = calloc(N,sizeof(double));
double *dx = calloc(N,sizeof(double));
double *dy = calloc(N,sizeof(double));
@@ -69,7 +69,12 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
theta1[i] += z1[j]*foo(s[i],range,theta0,j+1);
theta2[i] += z2[j]*foo(s[i],range,theta0,j+1);
}
- T[i] = fun(s[i],params);
+
+ T = fun(s[i],params);
+ E = T + m;
+ mom = sqrt(E*E - m*m);
+ beta[i] = mom/E;
+
if (i > 0) {
theta = sqrt(theta1[i]*theta1[i] + theta2[i]*theta2[i]);
phi = atan2(theta2[i],theta1[i]);
@@ -79,11 +84,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
dz[i] = (s[i]-s[i-1])*cos(theta);
/* Calculate total energy */
- E = T[i] + m;
- mom = sqrt(E*E - m*m);
- beta = mom/E;
-
- t[i] = t[i-1] + (s[i]-s[i-1])/(beta*SPEED_OF_LIGHT);
+ t[i] = t[i-1] + (s[i]-s[i-1])/(beta[i]*SPEED_OF_LIGHT);
x[i] = x[i-1] + dx[i];
y[i] = y[i-1] + dy[i];
@@ -168,7 +169,7 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
p->x = x;
p->y = y;
p->z = z;
- p->T = T;
+ p->beta = beta;
p->t = t;
p->dx = dx;
p->dy = dy;
@@ -180,13 +181,13 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
return p;
}
-void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0)
+void path_eval(path *p, double s, double *pos, double *dir, double *beta, double *t, double *theta0)
{
pos[0] = interp1d(s,p->s,p->x,N);
pos[1] = interp1d(s,p->s,p->y,N);
pos[2] = interp1d(s,p->s,p->z,N);
- *T = interp1d(s,p->s,p->T,N);
+ *beta = interp1d(s,p->s,p->beta,N);
*t = interp1d(s,p->s,p->t,N);
/* Technically we should be interpolating between the direction vectors
@@ -214,7 +215,7 @@ void path_free(path *p)
free(p->x);
free(p->y);
free(p->z);
- free(p->T);
+ free(p->beta);
free(p->t);
free(p->dx);
free(p->dy);
diff --git a/src/path.h b/src/path.h
index 449088f..aa7fa97 100644
--- a/src/path.h
+++ b/src/path.h
@@ -27,7 +27,7 @@ typedef struct path {
double *x;
double *y;
double *z;
- double *T;
+ double *beta;
double *t;
double *dx;
double *dy;
@@ -36,7 +36,7 @@ typedef struct path {
double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, size_t n);
path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, void *params, double *z1, double *z2, size_t n, double m);
-void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0);
+void path_eval(path *p, double s, double *pos, double *dir, double *beta, double *t, double *theta0);
void path_free(path *p);
#endif
diff --git a/src/pdg.c b/src/pdg.c
index d6a2161..780fbc8 100644
--- a/src/pdg.c
+++ b/src/pdg.c
@@ -16,6 +16,9 @@ double get_scattering_rms(double x, double p, double beta, double z, double rho)
* since the radiation length is only discussed in terms of an
* electromagnetic shower induced by electrons (see Section 33.4.2).
*
+ * Update: MicroBooNE uses this formula for muons in this paper:
+ * https://arxiv.org/abs/1703.06187.
+ *
* See Equation 33.15 in
* http://pdg.lbl.gov/2018/reviews/rpp2018-rev-passage-particles-matter.pdf. */
if (x == 0.0) return 0.0;
diff --git a/src/pdg.h b/src/pdg.h
index ff62f7a..98046f5 100644
--- a/src/pdg.h
+++ b/src/pdg.h
@@ -4,7 +4,9 @@
#define SPEED_OF_LIGHT 29.9792458 /* cm/ns */
/* From http://pdg.lbl.gov/2017/AtomicNuclearProperties/HTML/water_liquid.html */
#define RADIATION_LENGTH 36.08 /* g/cm^2 */
+#define ELECTRON_MASS 0.5109989461 /* MeV */
#define MUON_MASS 105.6583745 /* MeV */
+#define PROTON_MASS 938.272081 /* MeV */
#define FINE_STRUCTURE_CONSTANT 7.297352566417e-3
double get_scattering_rms(double x, double p, double beta, double z, double rho);
diff --git a/src/proton.c b/src/proton.c
new file mode 100644
index 0000000..922c552
--- /dev/null
+++ b/src/proton.c
@@ -0,0 +1,180 @@
+#include <stdio.h>
+#include <errno.h>
+#include <string.h>
+#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_spline.h>
+#include <math.h>
+#include "optics.h"
+#include "quantum_efficiency.h"
+#include "solid_angle.h"
+#include "pdg.h"
+#include "vector.h"
+#include "proton.h"
+#include "sno.h"
+#include "scattering.h"
+#include "pmt_response.h"
+#include "misc.h"
+
+static int initialized = 0;
+
+static double *x, *dEdx, *csda_range;
+static size_t size;
+
+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("proton_water_liquid.txt", "r");
+
+ if (!f) {
+ fprintf(stderr, "failed to open proton_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 = 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 3:
+ dEdx[n] = value;
+ break;
+ case 4:
+ csda_range[n] = value;
+ break;
+ }
+ j += 1;
+ str = strtok(NULL," \n");
+ }
+
+ n += 1;
+ }
+
+ fclose(f);
+
+ 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);
+
+ initialized = 1;
+
+ return 0;
+
+err:
+ fclose(f);
+
+ return -1;
+}
+
+double proton_get_range(double T, double rho)
+{
+ /* Returns the approximate range a proton with kinetic energy `T` will travel
+ * in water before losing all of its energy. This range is interpolated
+ * based on data from the PDG which uses the continuous slowing down
+ * approximation.
+ *
+ * `T` should be in MeV, and `rho` should be in g/cm^3.
+ *
+ * Return value is in cm.
+ *
+ * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
+ if (!initialized) {
+ if (init()) {
+ exit(1);
+ }
+ }
+
+ return gsl_spline_eval(spline_range, T, acc_range)/rho;
+}
+
+double proton_get_dEdx(double T, double rho)
+{
+ /* Returns the approximate dE/dx for a proton 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->x[0]) return spline_dEdx->y[0];
+
+ return gsl_spline_eval(spline_dEdx, T, acc_dEdx)*rho;
+}
diff --git a/src/proton.h b/src/proton.h
new file mode 100644
index 0000000..8e19e9c
--- /dev/null
+++ b/src/proton.h
@@ -0,0 +1,9 @@
+#ifndef PROTON_H
+#define PROTON_H
+
+#include <stddef.h> /* for size_t */
+
+double proton_get_range(double T, double rho);
+double proton_get_dEdx(double T, double rho);
+
+#endif
diff --git a/src/proton_water_liquid.txt b/src/proton_water_liquid.txt
new file mode 100644
index 0000000..ef774db
--- /dev/null
+++ b/src/proton_water_liquid.txt
@@ -0,0 +1,140 @@
+PSTAR: Stopping Powers and Range Tables for Protons
+
+WATER, LIQUID
+
+Kinetic Electron. Nuclear Total CSDA Projected Detour
+Energy Stp. Pow. Stp. Pow. Stp. Pow. Range Range Factor
+MeV MeV cm2/g MeV cm2/g MeV cm2/g g/cm2 g/cm2
+
+1.000E-03 1.337E+02 4.315E+01 1.769E+02 6.319E-06 2.878E-06 0.4555
+1.500E-03 1.638E+02 3.460E+01 1.984E+02 8.969E-06 4.400E-06 0.4906
+2.000E-03 1.891E+02 2.927E+01 2.184E+02 1.137E-05 5.909E-06 0.5197
+2.500E-03 2.114E+02 2.557E+01 2.370E+02 1.357E-05 7.380E-06 0.5440
+3.000E-03 2.316E+02 2.281E+01 2.544E+02 1.560E-05 8.811E-06 0.5647
+4.000E-03 2.675E+02 1.894E+01 2.864E+02 1.930E-05 1.155E-05 0.5986
+5.000E-03 2.990E+02 1.631E+01 3.153E+02 2.262E-05 1.415E-05 0.6254
+6.000E-03 3.276E+02 1.439E+01 3.420E+02 2.567E-05 1.661E-05 0.6473
+7.000E-03 3.538E+02 1.292E+01 3.667E+02 2.849E-05 1.896E-05 0.6656
+8.000E-03 3.782E+02 1.175E+01 3.900E+02 3.113E-05 2.121E-05 0.6813
+9.000E-03 4.012E+02 1.080E+01 4.120E+02 3.363E-05 2.337E-05 0.6950
+1.000E-02 4.229E+02 1.000E+01 4.329E+02 3.599E-05 2.545E-05 0.7070
+1.250E-02 4.660E+02 8.485E+00 4.745E+02 4.150E-05 3.037E-05 0.7318
+1.500E-02 5.036E+02 7.400E+00 5.110E+02 4.657E-05 3.499E-05 0.7514
+1.750E-02 5.372E+02 6.581E+00 5.437E+02 5.131E-05 3.938E-05 0.7674
+2.000E-02 5.673E+02 5.939E+00 5.733E+02 5.578E-05 4.356E-05 0.7808
+2.250E-02 5.946E+02 5.421E+00 6.001E+02 6.005E-05 4.757E-05 0.7923
+2.500E-02 6.195E+02 4.993E+00 6.245E+02 6.413E-05 5.144E-05 0.8022
+2.750E-02 6.421E+02 4.633E+00 6.467E+02 6.806E-05 5.519E-05 0.8109
+3.000E-02 6.628E+02 4.325E+00 6.671E+02 7.187E-05 5.883E-05 0.8187
+3.500E-02 6.989E+02 3.826E+00 7.028E+02 7.916E-05 6.585E-05 0.8319
+4.000E-02 7.290E+02 3.437E+00 7.324E+02 8.613E-05 7.259E-05 0.8429
+4.500E-02 7.538E+02 3.126E+00 7.569E+02 9.284E-05 7.911E-05 0.8522
+5.000E-02 7.740E+02 2.870E+00 7.768E+02 9.935E-05 8.547E-05 0.8602
+5.500E-02 7.901E+02 2.655E+00 7.927E+02 1.057E-04 9.169E-05 0.8673
+6.000E-02 8.026E+02 2.473E+00 8.050E+02 1.120E-04 9.782E-05 0.8735
+6.500E-02 8.119E+02 2.316E+00 8.142E+02 1.182E-04 1.039E-04 0.8791
+7.000E-02 8.183E+02 2.178E+00 8.205E+02 1.243E-04 1.099E-04 0.8842
+7.500E-02 8.223E+02 2.058E+00 8.243E+02 1.303E-04 1.159E-04 0.8889
+8.000E-02 8.241E+02 1.951E+00 8.260E+02 1.364E-04 1.218E-04 0.8931
+8.500E-02 8.239E+02 1.855E+00 8.258E+02 1.425E-04 1.278E-04 0.8971
+9.000E-02 8.222E+02 1.769E+00 8.239E+02 1.485E-04 1.338E-04 0.9007
+9.500E-02 8.190E+02 1.691E+00 8.206E+02 1.546E-04 1.398E-04 0.9041
+1.000E-01 8.145E+02 1.620E+00 8.161E+02 1.607E-04 1.458E-04 0.9073
+1.250E-01 7.801E+02 1.343E+00 7.814E+02 1.920E-04 1.767E-04 0.9207
+1.500E-01 7.360E+02 1.152E+00 7.371E+02 2.249E-04 2.094E-04 0.9310
+1.750E-01 6.959E+02 1.010E+00 6.969E+02 2.598E-04 2.440E-04 0.9393
+2.000E-01 6.604E+02 9.016E-01 6.613E+02 2.966E-04 2.806E-04 0.9460
+2.250E-01 6.286E+02 8.152E-01 6.294E+02 3.354E-04 3.191E-04 0.9515
+2.500E-01 5.999E+02 7.447E-01 6.006E+02 3.761E-04 3.596E-04 0.9562
+2.750E-01 5.737E+02 6.855E-01 5.744E+02 4.186E-04 4.019E-04 0.9601
+3.000E-01 5.497E+02 6.351E-01 5.504E+02 4.631E-04 4.462E-04 0.9635
+3.500E-01 5.075E+02 5.545E-01 5.080E+02 5.577E-04 5.404E-04 0.9689
+4.000E-01 4.714E+02 4.928E-01 4.719E+02 6.599E-04 6.422E-04 0.9731
+4.500E-01 4.401E+02 4.439E-01 4.406E+02 7.697E-04 7.515E-04 0.9764
+5.000E-01 4.128E+02 4.043E-01 4.132E+02 8.869E-04 8.683E-04 0.9790
+5.500E-01 3.888E+02 3.715E-01 3.891E+02 1.012E-03 9.926E-04 0.9811
+6.000E-01 3.676E+02 3.438E-01 3.680E+02 1.144E-03 1.124E-03 0.9829
+6.500E-01 3.489E+02 3.201E-01 3.492E+02 1.283E-03 1.263E-03 0.9844
+7.000E-01 3.322E+02 2.996E-01 3.325E+02 1.430E-03 1.410E-03 0.9857
+7.500E-01 3.172E+02 2.817E-01 3.175E+02 1.584E-03 1.563E-03 0.9868
+8.000E-01 3.037E+02 2.658E-01 3.039E+02 1.745E-03 1.724E-03 0.9877
+8.500E-01 2.914E+02 2.516E-01 2.917E+02 1.913E-03 1.891E-03 0.9886
+9.000E-01 2.803E+02 2.390E-01 2.805E+02 2.088E-03 2.066E-03 0.9893
+9.500E-01 2.700E+02 2.276E-01 2.702E+02 2.270E-03 2.247E-03 0.9899
+1.000E+00 2.606E+02 2.173E-01 2.608E+02 2.458E-03 2.435E-03 0.9905
+1.250E+00 2.228E+02 1.775E-01 2.229E+02 3.499E-03 3.472E-03 0.9925
+1.500E+00 1.955E+02 1.504E-01 1.957E+02 4.698E-03 4.669E-03 0.9938
+1.750E+00 1.748E+02 1.307E-01 1.749E+02 6.052E-03 6.020E-03 0.9946
+2.000E+00 1.585E+02 1.157E-01 1.586E+02 7.555E-03 7.519E-03 0.9952
+2.250E+00 1.453E+02 1.038E-01 1.454E+02 9.203E-03 9.164E-03 0.9957
+2.500E+00 1.343E+02 9.428E-02 1.344E+02 1.099E-02 1.095E-02 0.9960
+2.750E+00 1.250E+02 8.637E-02 1.251E+02 1.292E-02 1.288E-02 0.9963
+3.000E+00 1.171E+02 7.972E-02 1.172E+02 1.499E-02 1.494E-02 0.9965
+3.500E+00 1.041E+02 6.916E-02 1.042E+02 1.952E-02 1.946E-02 0.9968
+4.000E+00 9.398E+01 6.113E-02 9.404E+01 2.458E-02 2.451E-02 0.9971
+4.500E+00 8.580E+01 5.481E-02 8.586E+01 3.015E-02 3.007E-02 0.9973
+5.000E+00 7.906E+01 4.970E-02 7.911E+01 3.623E-02 3.613E-02 0.9974
+5.500E+00 7.339E+01 4.549E-02 7.343E+01 4.279E-02 4.268E-02 0.9975
+6.000E+00 6.854E+01 4.195E-02 6.858E+01 4.984E-02 4.972E-02 0.9976
+6.500E+00 6.434E+01 3.894E-02 6.438E+01 5.737E-02 5.724E-02 0.9977
+7.000E+00 6.068E+01 3.634E-02 6.071E+01 6.537E-02 6.522E-02 0.9977
+7.500E+00 5.744E+01 3.407E-02 5.747E+01 7.384E-02 7.368E-02 0.9978
+8.000E+00 5.456E+01 3.208E-02 5.460E+01 8.277E-02 8.259E-02 0.9978
+8.500E+00 5.199E+01 3.031E-02 5.202E+01 9.215E-02 9.196E-02 0.9979
+9.000E+00 4.966E+01 2.873E-02 4.969E+01 1.020E-01 1.018E-01 0.9979
+9.500E+00 4.756E+01 2.731E-02 4.759E+01 1.123E-01 1.120E-01 0.9979
+1.000E+01 4.564E+01 2.603E-02 4.567E+01 1.230E-01 1.228E-01 0.9980
+1.250E+01 3.813E+01 2.111E-02 3.815E+01 1.832E-01 1.828E-01 0.9981
+1.500E+01 3.290E+01 1.778E-02 3.292E+01 2.539E-01 2.535E-01 0.9982
+1.750E+01 2.904E+01 1.538E-02 2.905E+01 3.350E-01 3.344E-01 0.9982
+2.000E+01 2.605E+01 1.356E-02 2.607E+01 4.260E-01 4.252E-01 0.9983
+2.500E+01 2.174E+01 1.098E-02 2.175E+01 6.370E-01 6.359E-01 0.9983
+2.750E+01 2.012E+01 1.003E-02 2.013E+01 7.566E-01 7.553E-01 0.9984
+3.000E+01 1.875E+01 9.239E-03 1.876E+01 8.853E-01 8.839E-01 0.9984
+3.500E+01 1.656E+01 7.983E-03 1.656E+01 1.170E+00 1.168E+00 0.9984
+4.000E+01 1.487E+01 7.034E-03 1.488E+01 1.489E+00 1.486E+00 0.9985
+4.500E+01 1.353E+01 6.290E-03 1.354E+01 1.841E+00 1.839E+00 0.9985
+5.000E+01 1.244E+01 5.691E-03 1.245E+01 2.227E+00 2.224E+00 0.9985
+5.500E+01 1.154E+01 5.199E-03 1.154E+01 2.644E+00 2.641E+00 0.9985
+6.000E+01 1.078E+01 4.786E-03 1.078E+01 3.093E+00 3.089E+00 0.9986
+6.500E+01 1.012E+01 4.435E-03 1.013E+01 3.572E+00 3.567E+00 0.9986
+7.000E+01 9.555E+00 4.134E-03 9.559E+00 4.080E+00 4.075E+00 0.9986
+7.500E+01 9.059E+00 3.871E-03 9.063E+00 4.618E+00 4.611E+00 0.9986
+8.000E+01 8.622E+00 3.641E-03 8.625E+00 5.184E+00 5.176E+00 0.9986
+8.500E+01 8.233E+00 3.437E-03 8.236E+00 5.777E+00 5.769E+00 0.9986
+9.000E+01 7.884E+00 3.255E-03 7.888E+00 6.398E+00 6.389E+00 0.9986
+9.500E+01 7.570E+00 3.092E-03 7.573E+00 7.045E+00 7.035E+00 0.9986
+1.000E+02 7.286E+00 2.944E-03 7.289E+00 7.718E+00 7.707E+00 0.9987
+1.250E+02 6.190E+00 2.381E-03 6.192E+00 1.146E+01 1.144E+01 0.9987
+1.500E+02 5.443E+00 2.001E-03 5.445E+00 1.577E+01 1.576E+01 0.9987
+1.750E+02 4.901E+00 1.728E-03 4.903E+00 2.062E+01 2.060E+01 0.9988
+2.000E+02 4.491E+00 1.522E-03 4.492E+00 2.596E+01 2.593E+01 0.9988
+2.250E+02 4.169E+00 1.361E-03 4.170E+00 3.174E+01 3.171E+01 0.9988
+2.500E+02 3.910E+00 1.231E-03 3.911E+00 3.794E+01 3.790E+01 0.9989
+2.750E+02 3.697E+00 1.124E-03 3.698E+00 4.452E+01 4.447E+01 0.9989
+3.000E+02 3.519E+00 1.035E-03 3.520E+00 5.145E+01 5.139E+01 0.9989
+3.500E+02 3.240E+00 8.936E-04 3.241E+00 6.628E+01 6.621E+01 0.9989
+4.000E+02 3.031E+00 7.870E-04 3.032E+00 8.225E+01 8.217E+01 0.9990
+4.500E+02 2.870E+00 7.037E-04 2.871E+00 9.921E+01 9.912E+01 0.9990
+5.000E+02 2.743E+00 6.367E-04 2.743E+00 1.170E+02 1.169E+02 0.9990
+5.500E+02 2.640E+00 5.816E-04 2.640E+00 1.356E+02 1.355E+02 0.9991
+6.000E+02 2.555E+00 5.355E-04 2.556E+00 1.549E+02 1.547E+02 0.9991
+6.500E+02 2.485E+00 4.963E-04 2.485E+00 1.747E+02 1.746E+02 0.9991
+7.000E+02 2.426E+00 4.626E-04 2.426E+00 1.951E+02 1.949E+02 0.9991
+7.500E+02 2.375E+00 4.333E-04 2.376E+00 2.159E+02 2.158E+02 0.9991
+8.000E+02 2.333E+00 4.076E-04 2.333E+00 2.372E+02 2.370E+02 0.9992
+8.500E+02 2.296E+00 3.849E-04 2.296E+00 2.588E+02 2.586E+02 0.9992
+9.000E+02 2.264E+00 3.646E-04 2.264E+00 2.807E+02 2.805E+02 0.9992
+9.500E+02 2.236E+00 3.464E-04 2.236E+00 3.029E+02 3.027E+02 0.9992
+1.000E+03 2.211E+00 3.300E-04 2.211E+00 3.254E+02 3.252E+02 0.9992
+1.500E+03 2.070E+00 2.249E-04 2.070E+00 5.605E+02 5.601E+02 0.9993
+2.000E+03 2.021E+00 1.715E-04 2.021E+00 8.054E+02 8.049E+02 0.9994
+2.500E+03 2.004E+00 1.390E-04 2.004E+00 1.054E+03 1.053E+03 0.9995
+3.000E+03 2.001E+00 1.171E-04 2.001E+00 1.304E+03 1.303E+03 0.9995
+4.000E+03 2.012E+00 8.939E-05 2.012E+00 1.802E+03 1.802E+03 0.9996
+5.000E+03 2.031E+00 7.251E-05 2.031E+00 2.297E+03 2.296E+03 0.9996
+6.000E+03 2.052E+00 6.112E-05 2.052E+00 2.787E+03 2.786E+03 0.9997
+7.000E+03 2.072E+00 5.291E-05 2.072E+00 3.272E+03 3.271E+03 0.9997
+8.000E+03 2.091E+00 4.669E-05 2.091E+00 3.752E+03 3.751E+03 0.9997
+9.000E+03 2.109E+00 4.181E-05 2.109E+00 4.228E+03 4.227E+03 0.9997
+1.000E+04 2.126E+00 3.788E-05 2.126E+00 4.700E+03 4.699E+03 0.9998
diff --git a/src/test-path.c b/src/test-path.c
index 523ebd4..62994d6 100644
--- a/src/test-path.c
+++ b/src/test-path.c
@@ -139,7 +139,7 @@ void usage(void)
int main(int argc, char **argv)
{
size_t i, j, n, N, z, bins;
- double pos0[3], dir0[3], theta0, range, pos2[3], dir[3], T, t, tmp;
+ double pos0[3], dir0[3], theta0, range, pos2[3], dir[3], beta, t, tmp;
double *theta1, *theta2;
z = 2;
@@ -211,7 +211,7 @@ int main(int argc, char **argv)
gsl_histogram_increment(h2[j], simvalue->z2[j]);
}
for (j = 0; j < N; j++) {
- path_eval(simvalue->p,range*j/(N-1),pos2,dir,&T,&t,&tmp);
+ path_eval(simvalue->p,range*j/(N-1),pos2,dir,&beta,&t,&tmp);
theta1[j*n+i] = simvalue->dir[j*3+1] - dir[1];
theta2[j*n+i] = simvalue->dir[j*3+2] - dir[2];
}
@@ -282,7 +282,7 @@ int main(int argc, char **argv)
fprintf(pipe, "e\n");
sim *simvalue = simulate_path(pos0,dir0,range,theta0,N,z);
for (i = 0; i < N; i++) {
- path_eval(simvalue->p,range*i/(N-1),pos2,dir,&T,&t,&tmp);
+ path_eval(simvalue->p,range*i/(N-1),pos2,dir,&beta,&t,&tmp);
fprintf(pipe, "%g %g\n", i*range/(N-1), tmp);
}
@@ -307,7 +307,7 @@ int main(int argc, char **argv)
fprintf(pipe,"e\n");
for (i = 0; i < N; i++) {
- path_eval(simvalue->p,i*range/(N-1),pos2,dir,&T,&t,&theta0);
+ path_eval(simvalue->p,i*range/(N-1),pos2,dir,&beta,&t,&theta0);
fprintf(pipe,"%.10g %.10g %.10g\n", pos2[0], pos2[1], pos2[2]);
}
fprintf(pipe,"e\n");
diff --git a/src/test.c b/src/test.c
index 1a30976..2645982 100644
--- a/src/test.c
+++ b/src/test.c
@@ -14,6 +14,10 @@
#include <gsl/gsl_spline.h>
#include "vector.h"
#include <gsl/gsl_statistics.h>
+#include "electron.h"
+#include "proton.h"
+#include "likelihood.h"
+#include "id_particles.h"
typedef int testFunction(char *err);
@@ -111,17 +115,89 @@ struct solid_angle_results {
{2.0,2.0,0.282707}
};
+int test_proton_get_energy(char *err)
+{
+ /* A very simple test to make sure the energy as a function of distance
+ * along the track makes sense. Should include more detailed tests later. */
+ double T0, T;
+ particle *p;
+
+ /* Assume initial kinetic energy is 1 GeV. */
+ T0 = 1000.0;
+ p = particle_init(IDP_PROTON, T0,1.0,10000);
+ T = particle_get_energy(1e-9,p);
+
+ /* At the beginning of the track we should have roughly the same energy. */
+ if (!isclose(T, T0, 1e-5, 0)) {
+ sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
+ goto err;
+ }
+
+ T = particle_get_energy(p->range,p);
+
+ /* At the end of the track the energy should be approximately 0. */
+ if (!isclose(T, 0, 1e-5, 1e-5)) {
+ sprintf(err, "KE = %.5f, but expected %.5f", T, 0.0);
+ goto err;
+ }
+
+ particle_free(p);
+
+ return 0;
+
+err:
+ particle_free(p);
+
+ return 1;
+}
+
+int test_electron_get_energy(char *err)
+{
+ /* A very simple test to make sure the energy as a function of distance
+ * along the track makes sense. Should include more detailed tests later. */
+ double T0, T;
+ particle *p;
+
+ /* Assume initial kinetic energy is 1 GeV. */
+ T0 = 1000.0;
+ p = particle_init(IDP_E_MINUS, T0,1.0,10000);
+ T = particle_get_energy(1e-9,p);
+
+ /* At the beginning of the track we should have roughly the same energy. */
+ if (!isclose(T, T0, 1e-5, 0)) {
+ sprintf(err, "KE = %.5f, but expected %.5f", T, T0);
+ goto err;
+ }
+
+ T = particle_get_energy(p->range,p);
+
+ /* At the end of the track the energy should be approximately 0. */
+ if (!isclose(T, 0, 1e-5, 1e-5)) {
+ sprintf(err, "KE = %.5f, but expected %.5f", T, 0.0);
+ goto err;
+ }
+
+ particle_free(p);
+
+ return 0;
+
+err:
+ particle_free(p);
+
+ return 1;
+}
+
int test_muon_get_energy(char *err)
{
/* A very simple test to make sure the energy as a function of distance
* along the track makes sense. Should include more detailed tests later. */
- double T0, T, range;
- muon_energy *m;
+ double T0, T;
+ particle *p;
/* Assume initial kinetic energy is 1 GeV. */
T0 = 1000.0;
- m = muon_init_energy(T0,1.0,10000);
- T = muon_get_energy(1e-9,m);
+ p = particle_init(IDP_MU_MINUS, T0,1.0,10000);
+ T = particle_get_energy(1e-9,p);
/* At the beginning of the track we should have roughly the same energy. */
if (!isclose(T, T0, 1e-5, 0)) {
@@ -129,8 +205,7 @@ int test_muon_get_energy(char *err)
goto err;
}
- range = get_range(T0,1.0);
- T = muon_get_energy(range,m);
+ T = particle_get_energy(p->range,p);
/* At the end of the track the energy should be approximately 0. */
if (!isclose(T, 0, 1e-5, 1e-5)) {
@@ -138,15 +213,44 @@ int test_muon_get_energy(char *err)
goto err;
}
- muon_free_energy(m);
+ particle_free(p);
return 0;
err:
- muon_free_energy(m);
+ particle_free(p);
return 1;
+}
+int test_proton_get_range(char *err)
+{
+ /* A very simple test to make sure we read in the PDG table correctly. */
+ double value;
+
+ value = proton_get_range(1.0,1.0);
+
+ if (!isclose(value, 2.458E-03, 1e-5, 0)) {
+ sprintf(err, "range = %.5f, but expected %.5f", value, 2.458E-03);
+ return 1;
+ }
+
+ return 0;
+}
+
+int test_electron_get_range(char *err)
+{
+ /* A very simple test to make sure we read in the PDG table correctly. */
+ double value;
+
+ value = electron_get_range(1.0,1.0);
+
+ if (!isclose(value, 4.367e-01, 1e-5, 0)) {
+ sprintf(err, "range = %.5f, but expected %.5f", value, 4.367e-01);
+ return 1;
+ }
+
+ return 0;
}
int test_muon_get_range(char *err)
@@ -154,7 +258,7 @@ int test_muon_get_range(char *err)
/* A very simple test to make sure we read in the PDG table correctly. */
double value;
- value = get_range(1.0,1.0);
+ value = muon_get_range(1.0,1.0);
if (!isclose(value, 2.071e-3, 1e-5, 0)) {
sprintf(err, "range = %.5f, but expected %.5f", value, 1.863e-3);
@@ -164,12 +268,42 @@ int test_muon_get_range(char *err)
return 0;
}
+int test_proton_get_dEdx(char *err)
+{
+ /* A very simple test to make sure we read in the PDG table correctly. */
+ double value;
+
+ value = proton_get_dEdx(1.0,1.0);
+
+ if (!isclose(value, 2.608E+02, 1e-5, 0)) {
+ sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 2.608E+02);
+ return 1;
+ }
+
+ return 0;
+}
+
+int test_electron_get_dEdx(char *err)
+{
+ /* A very simple test to make sure we read in the PDG table correctly. */
+ double value;
+
+ value = electron_get_dEdx(1.0,1.0);
+
+ if (!isclose(value, 1.862, 1e-5, 0)) {
+ sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 1.862);
+ return 1;
+ }
+
+ return 0;
+}
+
int test_muon_get_dEdx(char *err)
{
/* A very simple test to make sure we read in the PDG table correctly. */
double value;
- value = get_dEdx(1.0,1.0);
+ value = muon_get_dEdx(1.0,1.0);
if (!isclose(value, 5.471, 1e-5, 0)) {
sprintf(err, "dE/dx = %.5f, but expected %.5f", value, 6.097);
@@ -341,10 +475,10 @@ int test_path(char *err)
/* Tests the logsumexp() function. */
int i, j, k;
double pos0[3], dir0[3], T0, range, m;
- double T, t;
+ double beta, t;
double pos_expected[3], t_expected;
double pos[3], dir[3];
- double E, mom, beta;
+ double E, mom, beta0;
double s;
double theta0;
@@ -372,11 +506,11 @@ int test_path(char *err)
/* Calculate total energy */
E = T0 + m;
mom = sqrt(E*E - m*m);
- beta = mom/E;
+ beta0 = mom/E;
- t_expected = s/(beta*SPEED_OF_LIGHT);
+ t_expected = s/(beta0*SPEED_OF_LIGHT);
- path_eval(p,s,pos,dir,&T,&t,&theta0);
+ path_eval(p,s,pos,dir,&beta,&t,&theta0);
for (k = 0; k < 3; k++) {
if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {
@@ -392,8 +526,8 @@ int test_path(char *err)
}
}
- if (!isclose(T, 1.0, 1e-9, 1e-9)) {
- sprintf(err, "path_eval(%.2g) returned T = %.5g, but expected %.5g", s, T, 1.0);
+ if (!isclose(beta, beta0, 1e-9, 1e-9)) {
+ sprintf(err, "path_eval(%.2g) returned beta = %.5g, but expected %.5g", s, beta, beta0);
return 1;
}
@@ -737,6 +871,12 @@ struct tests {
{test_get_path_length, "test_get_path_length"},
{test_mean, "test_mean"},
{test_std, "test_std"},
+ {test_electron_get_dEdx, "test_electron_get_dEdx"},
+ {test_electron_get_range, "test_electron_get_range"},
+ {test_electron_get_energy, "test_electron_get_energy"},
+ {test_proton_get_dEdx, "test_proton_get_dEdx"},
+ {test_proton_get_range, "test_proton_get_range"},
+ {test_proton_get_energy, "test_proton_get_energy"},
};
int main(int argc, char **argv)
generated by cgit (git 2.34.1) at 2025-05-06 22:22:38 +0000