move everything to src directory

1 files changed, 0 insertions, 260 deletions

diff --git a/muon.c b/muon.c
deleted file mode 100644
index 4b46769..0000000
--- a/muon.c
+++ /dev/null
@@ -1,260 +0,0 @@
-#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 "muon.h"
-#include "sno.h"
-#include "scattering.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 const double MUON_CRITICAL_ENERGY = 1.029e6;
-
-static int init()
-{
-    int i, j;
-    char line[256];
-    char *str;
-    double value;
-    int n;
-
-    FILE *f = fopen("muE_water_liquid.txt", "r");
-
-    if (!f) {
-        fprintf(stderr, "failed to open muE_water_liquid.txt: %s", 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 10 lines since it's just a header. */
-        if (i <= 10) continue;
-
-        if (!len) continue;
-        else if (line[0] == '#') continue;
-        else if (strstr(line, "Minimum ionization")) continue;
-        else if (strstr(line, "Muon critical energy")) 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 10 lines since it's just a header. */
-        if (i <= 10) continue;
-
-        if (!len) continue;
-        else if (line[0] == '#') continue;
-        else if (strstr(line, "Minimum ionization")) continue;
-        else if (strstr(line, "Muon critical energy")) continue;
-
-        str = strtok(line," \n");
-
-        j = 0;
-        while (str) {
-            value = strtod(str, NULL);
-            switch (j) {
-            case 0:
-                x[n] = value;
-                break;
-            case 7:
-                dEdx[n] = value;
-                break;
-            case 8:
-                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 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
-     * 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 get_T(double T0, double x, double rho)
-{
-    /* Returns the approximate kinetic energy of a muon in water after
-     * travelling `x` cm with an initial kinetic energy `T`.
-     *
-     * `T` should be in MeV, `x` in cm, and `rho` in g/cm^3.
-     *
-     * Return value is in MeV.
-     *
-     * See http://pdg.lbl.gov/2018/AtomicNuclearProperties/adndt.pdf. */
-    double a, b, range, T;
-
-    if (!initialized) {
-        if (init()) {
-            exit(1);
-        }
-    }
-
-    range = get_range(T0, rho);
-
-    /* This comes from Equation 33.42 in the PDG Passage of Particles Through
-     * Matter article. */
-    b = log(1 + T0/MUON_CRITICAL_ENERGY)/range;
-    /* Now we compute the ionization energy loss from the known range and b. */
-    a = b*T0/(exp(b*range)-1.0);
-
-    /* Compute the kinetic energy after travelling a distance `x` in the
-     * continuous slowing down approximation. */
-    T = -a/b + (T0+a/b)*exp(-b*x);
-
-    if (T < 0) return 0;
-
-    return T;
-}
-
-double get_dEdx(double T, double rho)
-{
-    /* Returns the approximate dE/dx for a muon 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);
-        }
-    }
-
-    return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho;
-}
-
-double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
-{
-    double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R;
-
-    z = 1.0;
-
-    SUB(pmt_dir,pmt_pos,pos);
-    normalize(pmt_dir);
-
-    if (DOT(pmt_dir,pmt_normal) > 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);
-
-    if (R <= AV_RADIUS) {
-        rho = HEAVY_WATER_DENSITY;
-    } else {
-        rho = WATER_DENSITY;
-    }
-
-    /* FIXME: I just calculate delta assuming 400 nm light. */
-    wavelength0 = 400.0;
-    n = get_index(rho, wavelength0, 10.0);
-
-    if (beta < 1/n) return 0;
-
-    /* FIXME: is this formula valid for muons? */
-    theta0 = get_scattering_rms(x,p,beta,z,rho);
-
-    /* FIXME: add angular response and scattering/absorption. */
-    return 2*omega*2*M_PI*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0)/(sqrt(2*M_PI)*theta0);
-}