add get_avg_index_{d2o,h2o} functions

4 files changed, 50 insertions, 28 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 3f27915..25563af 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -202,7 +202,7 @@ particle *particle_init(int id, double T0, size_t n)
     }
 
     if (p->shower_photons) {
-        norm = electron_get_angular_pdf_norm(p->a, p->b, 1/avg_index_d2o);
+        norm = electron_get_angular_pdf_norm(p->a, p->b, 1/get_avg_index_d2o());
 
         for (i = 0; i < n; i++) {
             p->cos_theta[i] = -1.0 + 2.0*i/(n-1);
@@ -215,7 +215,7 @@ particle *particle_init(int id, double T0, size_t n)
              *
              * Note: We add EPSILON to the pdf since the cdf values are
              * required to be strictly increasing in order to use gsl splines. */
-            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm + EPSILON;
+            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/get_avg_index_d2o())/norm + EPSILON;
             if (i > 0)
                 p->cdf_shower[i] = p->cdf_shower[i-1] + (pdf_last + pdf)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0;
             else
@@ -236,7 +236,7 @@ particle *particle_init(int id, double T0, size_t n)
 
     if (p->delta_ray_photons) {
         /* Calculate the normalization constant for the angular distribution. */
-        norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o);
+        norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/get_avg_index_d2o());
 
         for (i = 0; i < n; i++) {
             p->cos_theta[i] = -1.0 + 2.0*i/(n-1);
@@ -249,7 +249,7 @@ particle *particle_init(int id, double T0, size_t n)
              *
              * Note: We add EPSILON to the pdf since the cdf values are
              * required to be strictly increasing in order to use gsl splines. */
-            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm + EPSILON;
+            pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/get_avg_index_d2o())/norm + EPSILON;
             if (i > 0)
                 p->cdf_delta[i] = p->cdf_delta[i-1] + (pdf + pdf_last)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0;
             else
@@ -521,7 +521,7 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
 
     *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
 
-    *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+    *t = (l_d2o*get_avg_index_d2o() + l_h2o*get_avg_index_h2o())/SPEED_OF_LIGHT;
 
     *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
 }
@@ -553,7 +553,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
 
     R = NORM(pos);
 
-    n = (R <= AV_RADIUS) ? avg_index_d2o : avg_index_h2o;
+    n = (R <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o();
 
     *q = 0.0;
     *t = 0.0;
@@ -583,7 +583,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
 
     get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
 
-    *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+    *t = (l_d2o*get_avg_index_d2o() + l_h2o*get_avg_index_h2o())/SPEED_OF_LIGHT;
 
     /* Probability that a photon is absorbed. We calculate this by computing:
      *
@@ -872,7 +872,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     *t = 0.0;
     *mu_reflected = 0.0;
-    if (beta < 1/avg_index_d2o) return 0.0;
+    if (beta < 1/get_avg_index_d2o()) return 0.0;
 
     /* `prob` is the number of photons emitted per cm by the particle at a
      * distance `s` along the track. */
@@ -903,7 +903,7 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
 
     theta0 = fmax(theta0*sqrt(s),MIN_THETA0);
 
-    frac = sqrt(2)*avg_index_d2o*x*beta0*theta0;
+    frac = sqrt(2)*get_avg_index_d2o()*x*beta0*theta0;
 
     f = get_weighted_pmt_response(-cos_theta_pmt);
     f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
@@ -914,15 +914,15 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa
     prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS);
 
     /* Assume the particle is travelling at the speed of light. */
-    *t = s/SPEED_OF_LIGHT + l_d2o*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/SPEED_OF_LIGHT;
+    *t = s/SPEED_OF_LIGHT + l_d2o*get_avg_index_d2o()/SPEED_OF_LIGHT + l_h2o*get_avg_index_h2o()/SPEED_OF_LIGHT;
 
-    charge = avg_index_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+avg_index_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+avg_index_d2o*z*beta0)/frac))/(b+avg_index_d2o*beta0)/(4*M_PI);
+    charge = get_avg_index_d2o()*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+get_avg_index_d2o()*(smax-z)*beta0)/frac) + erf((-a+b*s+get_avg_index_d2o()*z*beta0)/frac))/(b+get_avg_index_d2o()*beta0)/(4*M_PI);
 
     /* Add expected number of photons from electromagnetic shower. */
     if (p->shower_photons > 0)
-        charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/avg_index_d2o)*omega/(2*M_PI);
+        charge += get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/get_avg_index_d2o())*omega/(2*M_PI);
     if (p->delta_ray_photons > 0)
-        charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/avg_index_d2o)*omega/(2*M_PI);
+        charge += get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/get_avg_index_d2o())*omega/(2*M_PI);
 
     *mu_reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
 
@@ -1049,7 +1049,7 @@ static double guess_time(double *pos, double *dir, double *pmt_pos, double smax,
     get_path_length(tmp,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
 
     /* Assume the particle is travelling at the speed of light. */
-    return s/SPEED_OF_LIGHT + l_d2o*avg_index_d2o/SPEED_OF_LIGHT + l_h2o*avg_index_h2o/SPEED_OF_LIGHT;
+    return s/SPEED_OF_LIGHT + l_d2o*get_avg_index_d2o()/SPEED_OF_LIGHT + l_h2o*get_avg_index_h2o()/SPEED_OF_LIGHT;
 }
 
 static double getKineticEnergy(double x, void *p)
@@ -1273,8 +1273,8 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, in
 
         /* Compute the distance at which the particle falls below the cerenkov
          * threshold. */
-        if (beta0 > 1/avg_index_d2o)
-            get_smax(p, 1/avg_index_d2o, range, &cerenkov_range);
+        if (beta0 > 1/get_avg_index_d2o())
+            get_smax(p, 1/get_avg_index_d2o(), range, &cerenkov_range);
         else
             cerenkov_range = range;
 
@@ -1298,7 +1298,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, int ns, const int fast, in
          *
          * These are computed at the beginning of this function and then passed to
          * the different functions to avoid recomputing them on the fly. */
-        cos_theta_cerenkov = 1/(avg_index_d2o*beta0);
+        cos_theta_cerenkov = 1/(get_avg_index_d2o()*beta0);
         sin_theta_cerenkov = sqrt(1-pow(cos_theta_cerenkov,2));
 
         mu_indirect[j] = 0.0;
diff --git a/src/optics.c b/src/optics.c
index 89bee66..7782678 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -29,11 +29,8 @@
 char optics_err[256];
 
 /* Index of refraction in water and heavy water averaged over the Cerenkov
- * spectrum and the PMT quantum efficiency.
- *
- * Set these to 0.0 so that it's obvious if someone forgets to call
- * optics_init() before using them. */
-double avg_index_h2o = 0.0, avg_index_d2o = 0.0;
+ * spectrum and the PMT quantum efficiency. */
+static double avg_index_h2o = 0.0, avg_index_d2o = 0.0;
 
 /* Absorption coefficients for H2O, D2O, and acrylic as a function of
  * wavelength from SNOMAN.
@@ -113,6 +110,26 @@ static double RIND_ACRYLIC_C3 = 0.32;
  * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */
 static double HC = 1.239841973976e3;
 
+double get_avg_index_d2o(void)
+{
+    if (!initialized) {
+        fprintf(stderr, "average index in d2o hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return avg_index_d2o;
+}
+
+double get_avg_index_h2o(void)
+{
+    if (!initialized) {
+        fprintf(stderr, "average index in h2o hasn't been initialized!\n");
+        exit(1);
+    }
+
+    return avg_index_h2o;
+}
+
 double get_index(double p, double wavelength, double T)
 {
     /* Returns the index of refraction of pure water for a given density,
diff --git a/src/optics.h b/src/optics.h
index 0e7740b..26e268e 100644
--- a/src/optics.h
+++ b/src/optics.h
@@ -23,16 +23,15 @@
 /* Global error string when optics_init() returns -1. */
 extern char optics_err[256];
 
-/* Index of refraction in water and heavy water averaged over the Cerenkov
- * spectrum and the PMT quantum efficiency.
- *
- * Note: You must call optics_init() before using these! */
-extern double avg_index_h2o, avg_index_d2o;
-
 /* Initialize the optics data by precomputing the average absorption and
  * scattering tables. */
 int optics_init(void);
 
+/* Functions to get the index of refraction in water and heavy water averaged
+ * over the Cerenkov spectrum and the PMT quantum efficiency. */
+double get_avg_index_d2o(void);
+double get_avg_index_h2o(void);
+
 /* Functions for computing the index of refraction. */
 double get_index(double p, double wavelength, double T);
 double get_index_snoman_h2o(double wavelength);
diff --git a/src/test.c b/src/test.c
index 5513d33..5ab4fba 100644
--- a/src/test.c
+++ b/src/test.c
@@ -147,6 +147,8 @@ int test_proton_get_energy(char *err)
     double T0, T;
     particle *p;
 
+    optics_init();
+
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
     p = particle_init(IDP_PROTON, T0, 10000);
@@ -183,6 +185,8 @@ int test_electron_get_energy(char *err)
     double T0, T;
     particle *p;
 
+    optics_init();
+
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
     p = particle_init(IDP_E_MINUS, T0, 10000);
@@ -219,6 +223,8 @@ int test_muon_get_energy(char *err)
     double T0, T;
     particle *p;
 
+    optics_init();
+
     /* Assume initial kinetic energy is 1 GeV. */
     T0 = 1000.0;
     p = particle_init(IDP_MU_MINUS, T0, 10000);