From 1d77bacaae25d40d160f2bcd14ba3a355921213e Mon Sep 17 00:00:00 2001 From: tlatorre Date: Sun, 27 Jan 2019 21:08:25 -0600 Subject: add photons from delta rays to likelihood calculation This commit updates the likelihood function to take into account Cerenkov light produced from delta rays produced by muons. The angular distribution of this light is currently assumed to be constant along the track and parameterized in the same way as the Cerenkov light from an electromagnetic shower. Currently I assume the light is produced uniformly along the track which isn't exactly correct, but should be good enough. --- src/likelihood.c | 149 ++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 97 insertions(+), 52 deletions(-) (limited to 'src/likelihood.c') diff --git a/src/likelihood.c b/src/likelihood.c index 786b221..8e12f57 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -42,16 +42,19 @@ particle *particle_init(int id, double T0, size_t n) * which we solve by dividing the track up into `n` segments and then * numerically computing the energy at each step. */ size_t i; - double dEdx; + double dEdx, rad; particle *p = malloc(sizeof(particle)); p->id = id; p->x = malloc(sizeof(double)*n); p->T = malloc(sizeof(double)*n); p->n = n; - p->rad = 0.0; p->a = 0.0; p->b = 0.0; + p->shower_photons = 0.0; + p->delta_ray_a = 0.0; + p->delta_ray_b = 0.0; + p->delta_ray_photons = 0.0; p->x[0] = 0; p->T[0] = T0; @@ -67,14 +70,17 @@ particle *particle_init(int id, double T0, size_t n) p->a = electron_get_angular_distribution_alpha(T0); p->b = electron_get_angular_distribution_beta(T0); + p->delta_ray_photons = 0.0; + /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ + rad = 0.0; for (i = 1; i < n; i++) { 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]); + 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 @@ -85,6 +91,8 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + p->shower_photons = rad*ELECTRON_PHOTONS_PER_MEV; + break; case IDP_MU_MINUS: case IDP_MU_PLUS: @@ -96,14 +104,21 @@ particle *particle_init(int id, double T0, size_t n) * load the correct table based on the media. */ p->range = muon_get_range(T0, HEAVY_WATER_DENSITY); + p->a = muon_get_angular_distribution_alpha(T0); + p->b = muon_get_angular_distribution_beta(T0); + + muon_get_delta_ray_distribution_parameters(T0, &p->delta_ray_a, &p->delta_ray_b); + p->delta_ray_photons = muon_get_delta_ray_photons(T0); + /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ + rad = 0.0; for (i = 1; i < n; i++) { 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]); + 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 +129,8 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + p->shower_photons = rad*MUON_PHOTONS_PER_MEV; + break; case IDP_PROTON: p->mass = PROTON_MASS; @@ -121,14 +138,17 @@ particle *particle_init(int id, double T0, size_t n) * tables for heavy water for protons. */ p->range = proton_get_range(T0, WATER_DENSITY); + p->delta_ray_photons = 0.0; + /* Now we loop over the points along the track and assume dE/dx is * constant between points. */ + rad = 0.0; for (i = 1; i < n; i++) { 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]); + 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 @@ -139,6 +159,8 @@ particle *particle_init(int id, double T0, size_t n) * range here using the dE/dx table instead of reading in the range. */ p->T[n-1] = 0; + p->shower_photons = rad*PROTON_PHOTONS_PER_MEV; + break; default: fprintf(stderr, "unknown particle id %i\n", id); @@ -170,9 +192,9 @@ 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) +static double get_expected_charge_shower(particle *p, 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 *q_delta_ray, double *q_indirect_delta_ray) { - double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter; + double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter, constant; SUB(pmt_dir,pmt_pos,pos); @@ -180,6 +202,9 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p cos_theta_pmt = DOT(pmt_dir,pmt_normal); + charge = 0.0; + *q_delta_ray = 0.0; + *q_indirect_delta_ray = 0.0; *reflected = 0.0; if (cos_theta_pmt > 0) return 0.0; @@ -201,11 +226,26 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p 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); + constant = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*omega/(2*M_PI); + + /* Note: We assume here that the peak of the angular distribution is at the + * Cerenkov angle for a particle with beta = 1. This seems to be an OK + * assumption for high energy showers, but is not exactly correct for + * showers from lower energy electrons or for delta rays from lower energy + * muons. It seems good enough, but in the future it would be nice to + * parameterize this. */ + if (p->shower_photons > 0) + charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/n_d2o); + + if (p->delta_ray_photons > 0) + *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/n_d2o); scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); *reflected = (f_reflec + 1.0 - scatter)*charge; + *q_indirect_delta_ray = (f_reflec + 1.0 - scatter)*(*q_delta_ray); + + *q_delta_ray *= f*scatter; return f*charge*scatter; } @@ -345,14 +385,14 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); } -static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad) +static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma) { 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; + double pos[3], n_d2o, n_h2o, wavelength0, t, l_d2o, l_h2o, q_delta_ray, q_indirect_delta_ray; if (n > MAX_NPOINTS) { fprintf(stderr, "number of points is greater than MAX_NPOINTS!\n"); @@ -364,6 +404,8 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 n_d2o = get_index_snoman_d2o(wavelength0); n_h2o = get_index_snoman_h2o(wavelength0); + double range = x[n-1] - x[0]; + for (i = 0; i < n; i++) { pos[0] = pos0[0] + x[i]*dir0[0]; pos[1] = pos0[1] + x[i]*dir0[1]; @@ -373,8 +415,9 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0 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, p->a, p->b, rad)*pdf[i]; - q_indirects[i] *= pdf[i]; + qs[i] = get_expected_charge_shower(p, pos, dir0, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o, &q_delta_ray, &q_indirect_delta_ray); + qs[i] = qs[i]*pdf[i] + q_delta_ray/range; + q_indirects[i] = q_indirects[i]*pdf[i] + q_indirect_delta_ray/range; ts[i] = t*qs[i]; ts2[i] = t*t*qs[i]; } @@ -593,10 +636,11 @@ static double get_total_charge_approx(double beta0, double *pos, double *dir, pa charge = abs_prob*n_d2o*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n_d2o*(smax-z)*beta0)/frac) + erf((-a+b*s+n_d2o*z*beta0)/frac))/(b+n_d2o*beta0)/(4*M_PI); - /* Add expected number of photons from electromagnetic shower for electrons - * and positrons. */ - if (p->id == IDP_E_MINUS || p->id == IDP_E_PLUS) - charge += abs_prob*get_weighted_quantum_efficiency()*p->rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); + /* Add expected number of photons from electromagnetic shower. */ + if (p->shower_photons > 0) + charge += abs_prob*get_weighted_quantum_efficiency()*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1.0/n_d2o)*omega/(2*M_PI); + if (p->delta_ray_photons > 0) + charge += abs_prob*get_weighted_quantum_efficiency()*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1.0/n_d2o)*omega/(2*M_PI); *mu_reflected = f_reflected*charge; @@ -891,6 +935,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast npoints = (size_t) (range/dx + 0.5); if (npoints < MIN_NPOINTS) npoints = MIN_NPOINTS; + if (npoints > MAX_NPOINTS) npoints = MAX_NPOINTS; /* Calculate total energy */ E0 = v[j].T0 + p->mass; @@ -903,32 +948,43 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast if (!fast) { path = path_init(v[j].pos, v[j].dir, v[j].T0, range, npoints, theta0, getKineticEnergy, p, v[j].z1, v[j].z2, v[j].n, p->mass); - if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) { + switch (v[j].id) { + case IDP_E_MINUS: + case IDP_E_PLUS: electron_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b); + break; + case IDP_MU_MINUS: + case IDP_MU_PLUS: + muon_get_position_distribution_parameters(v[j].T0, &pos_a, &pos_b); + break; + default: + fprintf(stderr, "unknown particle id %i\n", v[j].id); + exit(1); + } - /* Lower and upper limits to integrate for the electromagnetic - * shower. */ - xlo = 0.0; - xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); + /* Lower and upper limits to integrate for the electromagnetic + * shower. */ + xlo = 0.0; + xhi = gsl_cdf_gamma_Pinv(0.99,pos_a,pos_b); - /* Number of points to sample along the longitudinal direction for - * the electromagnetic shower. */ - npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5); + /* Number of points to sample along the longitudinal direction for + * the electromagnetic shower. */ + npoints_shower = (size_t) ((xhi-xlo)/dx_shower + 0.5); - if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS; + if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS; + if (npoints_shower > MAX_NPOINTS) npoints_shower = MAX_NPOINTS; - x = malloc(sizeof(double)*npoints_shower); - pdf = malloc(sizeof(double)*npoints_shower); - for (i = 0; i < npoints_shower; i++) { - x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1); - pdf[i] = gamma_pdf(x[i],pos_a,pos_b); - } + x = malloc(sizeof(double)*npoints_shower); + pdf = malloc(sizeof(double)*npoints_shower); + for (i = 0; i < npoints_shower; i++) { + x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1); + pdf[i] = gamma_pdf(x[i],pos_a,pos_b); + } - /* Make sure the PDF is normalized. */ - double norm = trapz(pdf,x[1]-x[0],npoints_shower); - for (i = 0; i < npoints_shower; i++) { - pdf[i] /= norm; - } + /* Make sure the PDF is normalized. */ + double norm = trapz(pdf,x[1]-x[0],npoints_shower); + for (i = 0; i < npoints_shower; i++) { + pdf[i] /= norm; } } @@ -982,27 +1038,16 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast ts[i][j] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,n_d2o,n_h2o,cos_theta_cerenkov,sin_theta_cerenkov); } - if (v[j].id == IDP_E_MINUS || v[j].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(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j],p->rad); - mu_indirect[j] += q_indirect; - ts_shower[i][j] = v[j].t0 + result; - } + integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]); + mu_indirect[j] += q_indirect; + ts_shower[i][j] = v[j].t0 + result; } if (!fast) { path_free(path); - if (v[j].id == IDP_E_MINUS || v[j].id == IDP_E_PLUS) { - free(x); - free(pdf); - } + free(x); + free(pdf); } particle_free(p); -- cgit