move everything to src directory

1 files changed, 193 insertions, 0 deletions

diff --git a/src/test-likelihood.c b/src/test-likelihood.c
new file mode 100644
index 0000000..f0266b9
--- /dev/null
+++ b/src/test-likelihood.c
@@ -0,0 +1,193 @@
+#include "muon.h"
+#include "random.h"
+#include "optics.h"
+#include "quantum_efficiency.h"
+#include <math.h>
+#include <gsl/gsl_histogram.h>
+#include "sno.h"
+#include "pdg.h"
+#include "vector.h"
+#include "solid_angle.h"
+#include <stdlib.h> /* for atoi() and strtod() */
+#include <unistd.h> /* for exit() */
+#include "scattering.h"
+#include <errno.h> /* for errno */
+#include <string.h> /* for strerror() */
+
+void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double theta0)
+{
+    /* Simulate the cos(theta) distribution around the original track direction
+     * for a muon with kinetic energy T. The angle from the original track
+     * distribution is simulated as a gaussian distribution with standard
+     * deviation `theta0`. */
+    int i;
+    double theta, phi, wavelength, u, qe, index, cerenkov_angle, dir[3], n[3], dest[3], E, p, beta, cos_theta;
+
+    i = 0;
+    while (i < N) {
+        /* Generate a random wavelength in the range 300-600 nm from the
+         * distribution of Cerenkov light. */
+        u = genrand_real2();
+        wavelength = 300.0*600.0/(u*(300.0-600.0) + 600.0);
+
+        qe = get_quantum_efficiency(wavelength);
+
+        /* Check to see if the photon was detected. */
+        if (genrand_real2() > qe) continue;
+
+        index = get_index(HEAVY_WATER_DENSITY, wavelength, 10.0);
+
+        /* Calculate total energy */
+        E = T + MUON_MASS;
+        p = sqrt(E*E - MUON_MASS*MUON_MASS);
+        beta = p/E;
+
+        cerenkov_angle = acos(1/(index*beta));
+
+        /* Assuming the muon track is dominated by small angle scattering, the
+         * angular distribution will be a Gaussian centered around 0 with a
+         * standard deviation of `theta0`. Here, we draw a random angle from
+         * this distribution. */
+        theta = randn()*theta0;
+
+        n[0] = sin(theta);
+        n[1] = 0;
+        n[2] = cos(theta);
+
+        /* To compute the direction of the photon, we start with a vector in
+         * the x-z plane which is offset from the track direction by the
+         * Cerenkov angle and then rotate it around the track direction by a
+         * random angle `phi`. */
+        dir[0] = sin(cerenkov_angle + theta);
+        dir[1] = 0;
+        dir[2] = cos(cerenkov_angle + theta);
+
+        phi = genrand_real2()*2*M_PI;
+
+        rotate(dest,dir,n,phi);
+
+        cos_theta = dest[2];
+
+        gsl_histogram_increment(h, cos_theta);
+
+        i += 1;
+    }
+}
+
+void usage(void)
+{
+    fprintf(stderr,"Usage: ./test-likelihood [options]\n");
+    fprintf(stderr,"  -n     number of events\n");
+    fprintf(stderr,"  -T     kinetic energy of muon (MeV)\n");
+    fprintf(stderr,"  -t     standard deviation of angular distribution\n");
+    fprintf(stderr,"  -b     number of bins\n");
+    fprintf(stderr,"  --xmin lowest value of cos(theta)\n");
+    fprintf(stderr,"  --xmax highest value of cos(theta)\n");
+    fprintf(stderr,"  -h     display this help message\n");
+    exit(1);
+}
+
+int main(int argc, char **argv)
+{
+    size_t i, N, bins;
+    double T, theta0;
+    double E, p, beta;
+    double xmin, xmax;
+
+    N = 100000;
+    bins = 1000;
+    T = 1000.0;
+    theta0 = 0.1;
+    xmin = -1.0;
+    xmax = 1.0;
+
+    for (i = 1; i < argc; i++) {
+        if (!strncmp(argv[i], "--", 2)) {
+            if (!strcmp(argv[i]+2,"xmin")) {
+                xmin = strtod(argv[++i],NULL);
+                continue;
+            } else if (!strcmp(argv[i]+2,"xmax")) {
+                xmax = strtod(argv[++i],NULL);
+                continue;
+            }
+        } else if (argv[i][0] == '-') {
+            switch (argv[i][1]) {
+            case 'n':
+                N = atoi(argv[++i]);
+                break;
+            case 'b':
+                bins = atoi(argv[++i]);
+                break;
+            case 'T':
+                T = strtod(argv[++i],NULL);
+                break;
+            case 't':
+                theta0 = strtod(argv[++i],NULL);
+                break;
+            case 'h':
+                usage();
+            default:
+                fprintf(stderr, "unrecognized option '%s'\n", argv[i]);
+                exit(1);
+            }
+        }
+    }
+
+    gsl_histogram *h = gsl_histogram_alloc(bins);
+    gsl_histogram_set_ranges_uniform(h,xmin,xmax);
+
+    simulate_cos_theta_distribution(N, h, T, theta0);
+
+    gsl_histogram_scale(h, 1.0/gsl_histogram_sum(h));
+
+    FILE *pipe = popen("graph -T X --bitmap-size 2000x2000 -X 'Cos(theta)' -Y Probability", "w");
+
+    if (!pipe) {
+        fprintf(stderr, "error running graph command: %s\n", strerror(errno));
+        exit(1);
+    }
+
+    for (i = 0; i < h->n; i++) {
+        fprintf(pipe, "%g %g\n", h->range[i], h->bin[i]);
+        fprintf(pipe, "%g %g\n", h->range[i+1], h->bin[i]);
+    }
+    fprintf(pipe, "\n\n");
+
+    gsl_histogram_reset(h);
+
+    init_interpolation();
+
+    /* Calculate total energy */
+    E = T + MUON_MASS;
+    p = sqrt(E*E - MUON_MASS*MUON_MASS);
+    beta = p/E;
+
+    for (i = 0; i < bins; i++) {
+        double lo, hi;
+        gsl_histogram_get_range(h, i, &lo, &hi);
+        double cos_theta = (lo+hi)/2.0;
+        h->bin[i] = get_probability(beta, cos_theta, theta0);
+    }
+
+    free_interpolation();
+
+    printf("\n\n");
+
+    gsl_histogram_scale(h, 1.0/gsl_histogram_sum(h));
+
+    for (i = 0; i < h->n; i++) {
+        fprintf(pipe, "%g %g\n", h->range[i], h->bin[i]);
+        fprintf(pipe, "%g %g\n", h->range[i+1], h->bin[i]);
+    }
+    fprintf(pipe, "\n\n");
+
+    if (pclose(pipe)) {
+        fprintf(stderr, "error closing graph command: %s\n", strerror(errno));
+        exit(1);
+    }
+
+    gsl_histogram_free(h);
+
+    return 0;
+}
+