diff options
author | tlatorre <tlatorre@uchicago.edu> | 2018-11-17 13:13:06 -0600 |
---|---|---|
committer | tlatorre <tlatorre@uchicago.edu> | 2018-11-17 13:13:06 -0600 |
commit | 8089acc95c1f337feaff2b8d3ae2f8ffd431234b (patch) | |
tree | 358ba5988ea85f64d6731a69edc7cc91598a2bc9 /src | |
parent | 497e627419d4d57ef5279e447509aacfcb1e9694 (diff) | |
download | sddm-8089acc95c1f337feaff2b8d3ae2f8ffd431234b.tar.gz sddm-8089acc95c1f337feaff2b8d3ae2f8ffd431234b.tar.bz2 sddm-8089acc95c1f337feaff2b8d3ae2f8ffd431234b.zip |
add some comments
Diffstat (limited to 'src')
-rw-r--r-- | src/likelihood.c | 28 | ||||
-rw-r--r-- | src/likelihood.h | 2 |
2 files changed, 28 insertions, 2 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 4e44b00..a24fb41 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -788,7 +788,10 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t else smax = 0.0; - /* Compute the Cerenkov angle at the start of the track. */ + /* Compute the Cerenkov angle at the start of the track. + * + * These are computed at the beginning of this function and then passed to + * the different functions to avoid recomputing them on the fly. */ wavelength0 = 400.0; n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); @@ -821,11 +824,23 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t * the time can be crazy since the error in the total charge is * in the denominator. Therefore, we estimate the time using * the same technique as in get_total_charge_approx(). See that - * function for more details. */ + * function for more details. + * + * Note: I don't really know how much this affects the likelihood. + * I should test it to see if we can get away with something even + * simpler (like just computing the distance from the start of the + * track to the PMT). */ ts[i] = t0 + guess_time(pos,dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov); } if (id == IDP_E_MINUS || id == IDP_E_PLUS) { + /* If we are fitting for an electron or positron, we need to + * estimate the expected number of photons produced by the + * electromagnetic shower. + * + * Technically we should also have a term for the electromagnetic + * shower produced by muons, but for the energy range I'm + * considering right now, it is very small. */ integrate_path_shower(x,pdf,T0,pos,dir,i,npoints_shower,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad); mu_indirect += q_indirect; ts_shower[i] = t0 + result; @@ -846,6 +861,9 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t mu_noise = DARK_RATE*GTVALID*1e-9; mu_indirect *= CHARGE_FRACTION/10000.0; + /* Compute the expected number of photons reaching each PMT by adding up + * the contributions from the noise hits and the direct, indirect, and + * shower light. */ for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; @@ -864,6 +882,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t else min_ratio = MIN_RATIO; + /* Now, we actually compute the negative log likelihood. */ nhit = 0; for (i = 0; i < MAX_PMTS; i++) { if (ev->pmt_hits[i].flags || pmts[i].pmt_type != PMT_NORMAL) continue; @@ -903,5 +922,10 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t for (i = 0; i < n; i++) logp_path += log_norm(z1[i],0,1) + log_norm(z2[i],0,1); + /* Add up the negative log likelihood terms from each PMT. I use a kahan + * sum here instead of just summing all the values to help prevent round + * off error. I actually don't think this makes any difference here since + * the numbers are all of the same order of magnitude, so I should actually + * test to see if it makes any difference. */ return kahan_sum(nll,nhit) - logp_path; } diff --git a/src/likelihood.h b/src/likelihood.h index f016b53..f5adf2f 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -4,6 +4,8 @@ #include "event.h" #include <stddef.h> /* for size_t */ +/* To avoid having to allocate new arrays every time we evaluate the likelihood + * function, we just allocate static arrays with `MAX_POINTS` elements. */ #define MAX_NPOINTS 1000 /* Maximum number of PE to consider if the PMT was hit. */ |