diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/fit.c | 173 | ||||
| -rw-r--r-- | src/path.c | 2 |
2 files changed, 147 insertions, 28 deletions
@@ -17,6 +17,7 @@ #include <errno.h> /* for errno */ #include "pdg.h" #include "optics.h" +#include "vector.h" #define EV_RECORD 0x45562020 // 'EV ' (as written to ZDAB file) @@ -36,15 +37,34 @@ static struct startingParameters { double x; double y; double z; - double theta; - double phi; } startingParameters[] = { - {0.0, 0.0, 0.0, 0.0, 0.0}, - {0.0, 0.0, 0.0, M_PI/2, 0.0}, - {0.0, 0.0, 0.0, M_PI/2, M_PI/2}, - {0.0, 0.0, 0.0, M_PI/2, M_PI}, - {0.0, 0.0, 0.0, M_PI/2, 3*M_PI/2}, - {0.0, 0.0, 0.0, M_PI, 0.0}, + { 0.0, 0.0, 0.0}, + { 100.0, 0.0, 0.0}, + { 100.0, 100.0, 0.0}, + { 0.0, 100.0, 0.0}, + {-100.0, 100.0, 0.0}, + {-100.0, 0.0, 0.0}, + {-100.0,-100.0, 0.0}, + { 0.0,-100.0, 0.0}, + { 100.0,-100.0, 0.0}, + { 0.0, 0.0, 100.0}, + { 100.0, 0.0, 100.0}, + { 100.0, 100.0, 100.0}, + { 0.0, 100.0, 100.0}, + {-100.0, 100.0, 100.0}, + {-100.0, 0.0, 100.0}, + {-100.0,-100.0, 100.0}, + { 0.0,-100.0, 100.0}, + { 100.0,-100.0, 100.0}, + { 0.0, 0.0, -100.0}, + { 100.0, 0.0, -100.0}, + { 100.0, 100.0, -100.0}, + { 0.0, 100.0, -100.0}, + {-100.0, 100.0, -100.0}, + {-100.0, 0.0, -100.0}, + {-100.0,-100.0, -100.0}, + { 0.0,-100.0, -100.0}, + { 100.0,-100.0, -100.0}, }; double nll(unsigned int n, const double *x, double *grad, void *params) @@ -96,25 +116,127 @@ double nll(unsigned int n, const double *x, double *grad, void *params) return fval; } +double guess_t0(event *ev, double *pos) +{ + /* Returns a guess for the t0 of event `ev` given a position `pos`. The t0 + * is calculated by computing the charge weighted average of the time + * difference between each PMT's hit time and the time it would take a + * photon to travel from `pos` to the PMT. */ + size_t i; + double pmt_dir[3], distance, t0, n, qhs_sum; + + /* Compute the index of refraction for heavy water at 400 nm. */ + n = get_index_snoman_d2o(400.0); + + t0 = 0.0; + qhs_sum = 0.0; + for (i = 0; i < MAX_PMTS; i++) { + /* Only look at normal PMTs that were hit and aren't flagged. */ + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue; + + /* Compute the vector between `pos` and the PMT. */ + SUB(pmt_dir,pmts[i].pos,pos); + + /* Compute the distance to the PMT. */ + distance = NORM(pmt_dir); + + /* Add the charge weighted time difference between the PMT hit time and + * the time it would take a photon to hit the PMT from `pos`. */ + t0 += ev->pmt_hits[i].qhs*(ev->pmt_hits[i].t - distance*n/SPEED_OF_LIGHT); + + qhs_sum += ev->pmt_hits[i].qhs; + } + + /* Divide by the total QHS sum. */ + t0 /= qhs_sum; + + return t0; +} + +void guess_direction(event *ev, double *pos, double *theta, double *phi) +{ + /* Returns the approximate direction of the event from a given position. + * The direction is computed by taking the charge weighted average of the + * vectors from `pos` to each hit PMT. + * + * A possible improvement here is to do something like a Hough transform + * where we map PMT hits to cones with a 42 degree angle and then search + * for the highest point in the transformed space. This method might also + * generalize better when searching for a second ring since we can look for + * secondary peaks not near the highest peak. */ + size_t i; + double pmt_dir[3], dir[3]; + + dir[0] = 0.0; + dir[1] = 0.0; + dir[2] = 0.0; + + for (i = 0; i < MAX_PMTS; i++) { + /* Only look at normal PMTs that were hit and aren't flagged. */ + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue; + + /* Compute the vector between `pos` and the PMT. */ + SUB(pmt_dir,pmts[i].pos,pos); + + /* Normalize it. */ + normalize(pmt_dir); + + /* Multiply the vector by the QHS charge in the PMT. */ + MUL(pmt_dir,ev->pmt_hits[i].qhs); + + /* Add this to the estimated direction. */ + ADD(dir,dir,pmt_dir); + } + + /* Normalize the charge weighted sum. */ + normalize(dir); + + /* Compute the spherical coordinates for `dir`. */ + *theta = acos(dir[2]); + *phi = atan2(dir[1],dir[0]); +} + int fit_event(event *ev, double *xopt, double *fmin) { size_t i; fitParams fpars; - double x[9], ss[9], lb[9], ub[9], fval, n; + double x[9], ss[9], lb[9], ub[9], fval, n, qhs_sum, x0[9], T0, Tmin; int rv; nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9); nlopt_set_min_objective(opt,nll,&fpars); - x[0] = 1200.0; - x[1] = 0.0; - x[2] = 0.0; - x[3] = 0.0; - x[4] = 1.57; - x[5] = 0.0; - x[6] = 130.0; - x[7] = 0.0; - x[8] = 0.0; + /* Make a guess as to the energy. Right now we just use a simple + * approximation that the muon produces approximately 6 photons/MeV. + * + * FIXME: Should update this to something better. */ + + qhs_sum = 0.0; + for (i = 0; i < MAX_PMTS; i++) { + if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL || !ev->pmt_hits[i].hit) continue; + qhs_sum += ev->pmt_hits[i].qhs; + } + + T0 = qhs_sum/6.0; + + n = get_index_snoman_d2o(400.0); + + /* Calculate the Cerenkov threshold for a muon. */ + Tmin = MUON_MASS/sqrt(1.0-1.0/(n*n)) - MUON_MASS; + + /* If our guess is below the Cerenkov threshold, start at the Cerenkov + * threshold. */ + if (T0 < Tmin) T0 = Tmin; + + x0[0] = T0; + x0[1] = 0.0; + x0[2] = 0.0; + x0[3] = 0.0; + x0[4] = 1.57; + x0[5] = 0.0; + x0[6] = 130.0; + x0[7] = 0.0; + x0[8] = 0.0; ss[0] = 10.0; ss[1] = 10.0; @@ -126,9 +248,7 @@ int fit_event(event *ev, double *xopt, double *fmin) ss[7] = 0.1; ss[8] = 0.1; - n = get_index_snoman_d2o(400.0); - - lb[0] = MUON_MASS/sqrt(1.0-1.0/(n*n)); + lb[0] = Tmin; lb[1] = -1000.0; lb[2] = -1000.0; lb[3] = -1000.0; @@ -167,15 +287,14 @@ int fit_event(event *ev, double *xopt, double *fmin) nlopt_set_maxeval(opt, 100); for (i = 0; i < sizeof(startingParameters)/sizeof(startingParameters[0]); i++) { - x[0] = 1200.0; + memcpy(x,x0,sizeof(x)); x[1] = startingParameters[i].x; x[2] = startingParameters[i].y; x[3] = startingParameters[i].z; - x[4] = startingParameters[i].theta; - x[5] = startingParameters[i].phi; - x[6] = 130.0; - x[7] = 0.0; - x[8] = 0.0; + + guess_direction(ev,x+1,&x[4],&x[5]); + + x[6] = guess_t0(ev,x+1); rv = nlopt_optimize(opt,x,&fval); @@ -8,7 +8,7 @@ #include <stdlib.h> /* for malloc(), calloc(), etc. */ #include "misc.h" -#define N 100 +#define N 1000 static double foo(double s, double range, double theta0, int k) { |