diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 318 |
1 files changed, 230 insertions, 88 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 2340fb1..f769216 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -21,6 +21,7 @@ #include "electron.h" #include "proton.h" #include "id_particles.h" +#include <gsl/gsl_cdf.h> particle *particle_init(int id, double T0, size_t n) { @@ -46,6 +47,7 @@ particle *particle_init(int id, double T0, size_t n) p->x = malloc(sizeof(double)*n); p->T = malloc(sizeof(double)*n); p->n = n; + p->rad = 0.0; p->x[0] = 0; p->T[0] = T0; @@ -64,6 +66,8 @@ particle *particle_init(int id, double T0, size_t n) p->x[i] = p->range*i/(n-1); dEdx = electron_get_dEdx(p->T[i-1], WATER_DENSITY); p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]); + dEdx = electron_get_dEdx_rad(p->T[i-1], WATER_DENSITY); + p->rad += dEdx*(p->x[i]-p->x[i-1]); } /* Make sure that the energy is zero at the last step. This is so that * when we try to bisect the point along the track where the speed of @@ -91,6 +95,8 @@ particle *particle_init(int id, double T0, size_t n) p->x[i] = p->range*i/(n-1); dEdx = muon_get_dEdx(p->T[i-1], HEAVY_WATER_DENSITY); p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]); + dEdx = muon_get_dEdx_rad(p->T[i-1], WATER_DENSITY); + p->rad += dEdx*(p->x[i]-p->x[i-1]); } /* Make sure that the energy is zero at the last step. This is so that * when we try to bisect the point along the track where the speed of @@ -114,6 +120,8 @@ particle *particle_init(int id, double T0, size_t n) p->x[i] = p->range*i/(n-1); dEdx = proton_get_dEdx(p->T[i-1], WATER_DENSITY); p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]); + dEdx = proton_get_dEdx_rad(p->T[i-1], WATER_DENSITY); + p->rad += dEdx*(p->x[i]-p->x[i-1]); } /* Make sure that the energy is zero at the last step. This is so that * when we try to bisect the point along the track where the speed of @@ -155,6 +163,50 @@ void particle_free(particle *p) free(p); } +static double get_expected_charge_shower(double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double alpha, double beta, double rad) +{ + double pmt_dir[3], cos_theta, n, omega, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; + + SUB(pmt_dir,pmt_pos,pos); + + normalize(pmt_dir); + + cos_theta_pmt = DOT(pmt_dir,pmt_normal); + + *reflected = 0.0; + if (cos_theta_pmt > 0) return 0.0; + + /* Calculate the cosine of the angle between the track direction and the + * vector to the PMT. */ + cos_theta = DOT(dir,pmt_dir); + + omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r); + + R = NORM(pos); + + if (R <= AV_RADIUS) { + n = n_d2o; + } else { + n = n_h2o; + } + + theta_pmt = acos(-cos_theta_pmt); + f_reflec = get_weighted_pmt_reflectivity(theta_pmt); + f = get_weighted_pmt_response(theta_pmt); + + absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); + absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); + absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic(); + + l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; + + charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI); + + *reflected = f_reflec*charge; + + return f*charge; +} + static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o) { double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; @@ -174,7 +226,7 @@ static double get_expected_charge(double x, double beta, double theta0, double * * vector to the PMT. */ cos_theta = DOT(dir,pmt_dir); - omega = get_solid_angle_approx(pos,pmt_pos,pmt_normal,r); + omega = get_solid_angle_fast(pos,pmt_pos,pmt_normal,r); R = NORM(pos); @@ -204,7 +256,7 @@ static double get_expected_charge(double x, double beta, double theta0, double * return f*charge; } -double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma) +double F(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the CDF for the time distribution of photons at time `t`. */ size_t i; @@ -217,11 +269,15 @@ double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_ p += mu_direct[i]*norm_cdf(t,ts[i],sigma); mu_total += mu_direct[i]; } + for (i = 0; i < n; i++) { + p += mu_shower[i]*norm_cdf(t,ts_shower[i],ts_sigma[i]); + mu_total += mu_shower[i]; + } return p/mu_total; } -double f(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma) +double f(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the probability that a photon is detected at time `t`. * @@ -247,11 +303,15 @@ double f(double t, double mu_noise, double mu_indirect, double *mu_direct, size_ p += mu_direct[i]*norm(t,ts[i],sigma); mu_total += mu_direct[i]; } + for (i = 0; i < n; i++) { + p += mu_shower[i]*norm(t,ts_shower[i],ts_sigma[i]); + mu_total += mu_shower[i]; + } return p/mu_total; } -double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, size_t n2, double *ts, double tmean, double sigma) +double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n2, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the first order statistic for observing a PMT hit at time `t` * given `n` hits. @@ -263,18 +323,73 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m * different transit times across the face of the PMT, it seems better to * convolve first which is what we do here. In addition, the problem is not * analytically tractable if you do things the other way around. */ - return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,n2,ts,tmean,sigma)); + return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); +} + +static void integrate_path_shower(double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad, double pos_a, double pos_b) +{ + size_t i; + static double qs[MAX_NPOINTS]; + static double q_indirects[MAX_NPOINTS]; + static double ts[MAX_NPOINTS]; + static double ts2[MAX_NPOINTS]; + double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, alpha, beta; + + if (n > MAX_NPOINTS) { + fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n"); + exit(1); + } + + alpha = electron_get_angular_distribution_alpha(T0); + beta = electron_get_angular_distribution_beta(T0); + + /* FIXME: I just calculate delta assuming 400 nm light. */ + wavelength0 = 400.0; + n_d2o = get_index_snoman_d2o(wavelength0); + n_h2o = get_index_snoman_h2o(wavelength0); + + for (i = 0; i < n; i++) { + pos[0] = pos0[0] + x[i]*dir0[0]; + pos[1] = pos0[1] + x[i]*dir0[1]; + pos[2] = pos0[2] + x[i]*dir0[2]; + + get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o); + + t = x[i]/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; + + qs[i] = get_expected_charge_shower(pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, alpha, beta, rad)*pdf[i]; + q_indirects[i] *= pdf[i]; + ts[i] = t*qs[i]; + ts2[i] = t*t*qs[i]; + } + + *mu_direct = trapz(qs,(b-a)/(n-1),n); + *mu_indirect = trapz(q_indirects,(b-a)/(n-1),n); + *time = trapz(ts,(b-a)/(n-1),n)/(*mu_direct); + double t2 = trapz(ts2,(b-a)/(n-1),n)/(*mu_direct); + + *sigma = sqrt(t2-(*time)*(*time)); + + if (*mu_direct == 0.0) { + *time = 0.0; + *sigma = PMT_TTS; + } else { + *sigma = sqrt((*sigma)*(*sigma) + PMT_TTS*PMT_TTS); + } } static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time) { size_t i; - double *qs, *q_indirects, *ts; + static double qs[MAX_NPOINTS]; + static double q_indirects[MAX_NPOINTS]; + static double ts[MAX_NPOINTS]; double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x; - qs = malloc(sizeof(double)*n); - q_indirects = malloc(sizeof(double)*n); - ts = malloc(sizeof(double)*n); + if (n > MAX_NPOINTS) { + fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n"); + exit(1); + } /* FIXME: I just calculate delta assuming 400 nm light. */ wavelength0 = 400.0; @@ -297,10 +412,6 @@ static void integrate_path(path *p, int pmt, double a, double b, size_t n, doubl *mu_direct = trapz(qs,(b-a)/(n-1),n); *mu_indirect = trapz(q_indirects,(b-a)/(n-1),n); *time = trapz(ts,(b-a)/(n-1),n); - - free(qs); - free(q_indirects); - free(ts); } static double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected, double n_d2o, double n_h2o, double cos_theta_cerenkov, double sin_theta_cerenkov) @@ -424,7 +535,7 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti sin_theta = sqrt(1-pow(cos_theta,2)); /* Get the solid angle of the PMT at the position `s` along the track. */ - omega = get_solid_angle_approx(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS); + omega = get_solid_angle_fast(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS); theta0 = fmax(theta0*sqrt(s),MIN_THETA0); @@ -526,7 +637,8 @@ static double getKineticEnergy(double x, void *p) double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, size_t npoints, int fast) { size_t i, j, nhit; - double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp; + static double logp[MAX_PE], nll[MAX_PMTS]; + double range, theta0, E0, p0, beta0, smax, log_mu, max_logp; particle *p; double pmt_dir[3], R, cos_theta, sin_theta, wavelength0, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov, s, l_h2o, l_d2o; double logp_path; @@ -535,11 +647,31 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t double mu_reflected; path *path; - double mu_direct[MAX_PMTS]; - double ts[MAX_PMTS]; - double mu[MAX_PMTS]; + static double mu_direct[MAX_PMTS]; + static double mu_shower[MAX_PMTS] = {0}; + static double ts[MAX_PMTS]; + static double ts_shower[MAX_PMTS] = {0}; + static double ts_sigma[MAX_PMTS] = {0}; + static double mu[MAX_PMTS]; double mu_noise, mu_indirect; + for (i = 0; i < MAX_PMTS; i++) ts_sigma[i] = PMT_TTS; + + double pos_a, pos_b; + + electron_get_position_distribution_parameters(T0, &pos_a, &pos_b); + + /* Upper limit to integrate for the electromagnetic shower. */ + double xlo = 0.0; + double xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); + + double *x = malloc(sizeof(double)*npoints); + double *pdf = malloc(sizeof(double)*npoints); + for (i = 0; i < npoints; i++) { + x[i] = xlo + i*(xhi-xlo)/(npoints-1); + pdf[i] = gamma_pdf(x[i],pos_a,pos_b); + } + double result; double min_ratio; @@ -582,80 +714,87 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected, n_d2o, n_h2o, cos_theta_cerenkov, sin_theta_cerenkov); ts[i] += t0; mu_indirect += mu_reflected; + continue; + } + + /* First, we try to compute the distance along the track where the + * PMT is at the Cerenkov angle. The reason for this is because for + * heavy particles like muons which don't scatter much, the + * probability distribution for getting a photon hit along the + * track looks kind of like a delta function, i.e. the PMT is only + * hit over a very narrow window when the angle between the track + * direction and the PMT is *very* close to the Cerenkov angle + * (it's not a perfect delta function since there is some width due + * to dispersion). In this case, it's possible that the numerical + * integration completely skips over the delta function and so + * predicts an expected charge of 0. + * + * To fix this, we only integrate 1 meter on both sides of the + * point along the track where the PMT is at the Cerenkov angle. */ + + /* First, we find the point along the track where the PMT is at the + * Cerenkov angle. */ + SUB(pmt_dir,pmts[i].pos,pos); + /* Compute the distance to the PMT. */ + R = NORM(pmt_dir); + normalize(pmt_dir); + + /* Calculate the cosine of the angle between the track direction and the + * vector to the PMT at the start of the track. */ + cos_theta = DOT(dir,pmt_dir); + sin_theta = sqrt(1-pow(cos_theta,2)); + + /* Now, we compute the distance along the track where the PMT is at the + * Cerenkov angle. + * + * Note: This formula comes from using the "Law of sines" where the three + * vertices of the triangle are the starting position of the track, the + * point along the track that we want to find, and the PMT position. */ + s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; + + /* Make sure that the point is somewhere along the track between 0 and + * `smax`. */ + if (s < 0) s = 0.0; + else if (s > smax) s = smax; + + /* Compute the endpoints of the integration. */ + a = s - 100.0; + b = s + 100.0; + + if (a < 0.0) a = 0.0; + if (b > range) b = range; + + double q_indirect; + integrate_path(path, i, a, b, npoints, mu_direct+i, &q_indirect, &result); + mu_indirect += q_indirect; + + if (mu_direct[i] > 1e-9) { + ts[i] = t0 + result/mu_direct[i]; } else { - /* First, we try to compute the distance along the track where the - * PMT is at the Cerenkov angle. The reason for this is because for - * heavy particles like muons which don't scatter much, the - * probability distribution for getting a photon hit along the - * track looks kind of like a delta function, i.e. the PMT is only - * hit over a very narrow window when the angle between the track - * direction and the PMT is *very* close to the Cerenkov angle - * (it's not a perfect delta function since there is some width due - * to dispersion). In this case, it's possible that the numerical - * integration completely skips over the delta function and so - * predicts an expected charge of 0. - * - * To fix this, we only integrate 1 meter on both sides of the - * point along the track where the PMT is at the Cerenkov angle. */ - - /* First, we find the point along the track where the PMT is at the - * Cerenkov angle. */ - SUB(pmt_dir,pmts[i].pos,pos); - /* Compute the distance to the PMT. */ - R = NORM(pmt_dir); - normalize(pmt_dir); - - /* Calculate the cosine of the angle between the track direction and the - * vector to the PMT at the start of the track. */ - cos_theta = DOT(dir,pmt_dir); - sin_theta = sqrt(1-pow(cos_theta,2)); - - /* Now, we compute the distance along the track where the PMT is at the - * Cerenkov angle. - * - * Note: This formula comes from using the "Law of sines" where the three - * vertices of the triangle are the starting position of the track, the - * point along the track that we want to find, and the PMT position. */ - s = R*(sin_theta_cerenkov*cos_theta - cos_theta_cerenkov*sin_theta)/sin_theta_cerenkov; - - /* Make sure that the point is somewhere along the track between 0 and - * `smax`. */ - if (s < 0) s = 0.0; - else if (s > smax) s = smax; - - /* Compute the endpoints of the integration. */ - a = s - 100.0; - b = s + 100.0; - - if (a < 0.0) a = 0.0; - if (b > range) b = range; - - double q_indirect; - integrate_path(path, i, a, b, npoints, mu_direct+i, &q_indirect, &result); - mu_indirect += q_indirect; + /* If the expected number of PE is very small, our estimate of + * 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. */ + n_h2o = get_index_snoman_h2o(wavelength0); - if (mu_direct[i] > 1e-9) { - ts[i] = t0 + result/mu_direct[i]; - } else { - /* If the expected number of PE is very small, our estimate of - * 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. */ - n_h2o = get_index_snoman_h2o(wavelength0); + /* Compute the position of the particle at a distance `s` along the track. */ + tmp[0] = pos[0] + s*dir[0]; + tmp[1] = pos[1] + s*dir[1]; + tmp[2] = pos[2] + s*dir[2]; - /* Compute the position of the particle at a distance `s` along the track. */ - tmp[0] = pos[0] + s*dir[0]; - tmp[1] = pos[1] + s*dir[1]; - tmp[2] = pos[2] + s*dir[2]; + SUB(pmt_dir,pmts[i].pos,tmp); - SUB(pmt_dir,pmts[i].pos,tmp); + get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); - get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); + /* Assume the particle is travelling at the speed of light. */ + ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; + } - /* Assume the particle is travelling at the speed of light. */ - ts[i] = t0 + s/SPEED_OF_LIGHT + l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT; - } + if (id == IDP_E_MINUS || id == IDP_E_PLUS) { + integrate_path_shower(x,pdf,T0,pos,dir,i,0,xhi,npoints,mu_shower+i,&q_indirect,&result,&ts_sigma[i],p->rad,pos_a,pos_b); + mu_indirect += q_indirect; + ts_shower[i] = t0 + result; } } @@ -677,7 +816,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t if (fast) mu[i] = mu_direct[i] + mu_noise; else - mu[i] = mu_direct[i] + mu_indirect + mu_noise; + mu[i] = mu_direct[i] + mu_shower[i] + mu_indirect + mu_noise; } if (fast) @@ -697,7 +836,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t if (ev->pmt_hits[i].hit) { for (j = 1; j < MAX_PE; j++) { - logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], 1, &ts[i], ts[i], PMT_TTS); + logp[j] = log(pq(ev->pmt_hits[i].qhs,j)) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect, &mu_direct[i], &mu_shower[i], 1, &ts[i], &ts_shower[i], ts[i], PMT_TTS, &ts_sigma[i]); if (j == 1 || logp[j] > max_logp) max_logp = logp[j]; if (logp[j] - max_logp < min_ratio*ln(10)) { @@ -724,5 +863,8 @@ 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); + free(x); + free(pdf); + return kahan_sum(nll,nhit) - logp_path; } |