update code to use get_index_snoman* functions to calculate the index of refraction

1 files changed, 10 insertions, 4 deletions

diff --git a/src/likelihood.c b/src/likelihood.c
index 4e24a79..a4864a8 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -78,7 +78,7 @@ static double gsl_muon_time(double x, void *params)
 {
     double *params2 = (double *) params;
     double T0 = params2[0];
-    double pos0[3], dir[3], pos[3], pmt_dir[3];
+    double pos0[3], dir[3], pos[3], pmt_dir[3], R, n, wavelength0;
     int i;
     double t;
     i = (int) params2[1];
@@ -93,13 +93,19 @@ static double gsl_muon_time(double x, void *params)
     pos[1] = pos0[1] + dir[1]*x;
     pos[2] = pos0[2] + dir[2]*x;
 
+    R = NORM(pos);
+
     SUB(pmt_dir,pmts[i].pos,pos);
 
     double distance = NORM(pmt_dir);
 
     /* FIXME: I just calculate delta assuming 400 nm light. */
-    double wavelength0 = 400.0;
-    double n = get_index(HEAVY_WATER_DENSITY, wavelength0, 10.0);
+    wavelength0 = 400.0;
+    if (R <= AV_RADIUS) {
+        n = get_index_snoman_d2o(wavelength0);
+    } else {
+        n = get_index_snoman_h2o(wavelength0);
+    }
 
     t = x/SPEED_OF_LIGHT + distance*n/SPEED_OF_LIGHT;
 
@@ -197,7 +203,7 @@ double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
          * just to split up the integral, we only need to find a point along
          * the track close enough such that the integral isn't completely zero.
          */
-        theta_cerenkov = acos(1/get_index(WATER_DENSITY,400.0,10.0));
+        theta_cerenkov = acos(1/get_index_snoman_d2o(400.0));
 
         /* Now, we compute the distance along the track where the PMT is at the
          * Cerenkov angle. */