diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 107 |
1 files changed, 56 insertions, 51 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 3b37bd5..61034a3 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -122,34 +122,6 @@ particle *particle_init(int id, double T0, size_t n) p->shower_photons = electron_get_shower_photons(T0, rad); - norm = electron_get_angular_pdf_norm(p->a, p->b, 1/avg_index_d2o); - - for (i = 0; i < n; i++) { - p->cos_theta[i] = -1.0 + 2.0*i/(n-1); - /* 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. */ - pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm; - if (i > 0) - p->cdf_shower[i] = p->cdf_shower[i-1] + (pdf_last + pdf)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; - else - p->cdf_shower[i] = 0.0; - pdf_last = pdf; - } - - gsl_spline_init(p->spline_shower, p->cdf_shower, p->cos_theta, n); - - p->xlo_shower = gsl_cdf_gamma_Pinv(0.01,p->pos_a,p->pos_b); - p->xhi_shower = gsl_cdf_gamma_Pinv(0.99,p->pos_a,p->pos_b); - - for (i = 0; i < n; i++) { - p->x_shower[i] = p->xlo_shower + (p->xhi_shower-p->xlo_shower)*i/(n-1); - p->gamma_pdf[i] = gamma_pdf(p->x_shower[i],p->pos_a,p->pos_b); - } - break; case IDP_MU_MINUS: case IDP_MU_PLUS: @@ -167,6 +139,7 @@ particle *particle_init(int id, double T0, size_t n) muon_get_position_distribution_parameters(T0, &p->pos_a, &p->pos_b); 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 @@ -190,27 +163,6 @@ particle *particle_init(int id, double T0, size_t n) p->shower_photons = muon_get_shower_photons(T0, rad); - /* Calculate the normalization constant for the angular distribution. */ - norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o); - - for (i = 0; i < n; i++) { - p->cos_theta[i] = -1.0 + 2.0*i/(n-1); - /* 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. */ - pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm; - if (i > 0) - p->cdf_delta[i] = p->cdf_delta[i-1] + (pdf + pdf_last)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; - else - p->cdf_delta[i] = 0.0; - pdf_last = pdf; - } - - gsl_spline_init(p->spline_delta, p->cdf_delta, p->cos_theta, n); - break; case IDP_PROTON: p->mass = PROTON_MASS; @@ -240,8 +192,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; - /* FIXME: Is this correct? */ - p->shower_photons = rad*PROTON_PHOTONS_PER_MEV; + /* FIXME: Add shower photons for protons. */ + p->shower_photons = 0.0; break; default: @@ -249,6 +201,59 @@ particle *particle_init(int id, double T0, size_t n) exit(1); } + if (p->shower_photons) { + norm = electron_get_angular_pdf_norm(p->a, p->b, 1/avg_index_d2o); + + for (i = 0; i < n; i++) { + p->cos_theta[i] = -1.0 + 2.0*i/(n-1); + /* 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. */ + pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->a, p->b, 1/avg_index_d2o)/norm; + if (i > 0) + p->cdf_shower[i] = p->cdf_shower[i-1] + (pdf_last + pdf)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; + else + p->cdf_shower[i] = 0.0; + pdf_last = pdf; + } + + gsl_spline_init(p->spline_shower, p->cdf_shower, p->cos_theta, n); + + p->xlo_shower = gsl_cdf_gamma_Pinv(0.01,p->pos_a,p->pos_b); + p->xhi_shower = gsl_cdf_gamma_Pinv(0.99,p->pos_a,p->pos_b); + + for (i = 0; i < n; i++) { + p->x_shower[i] = p->xlo_shower + (p->xhi_shower-p->xlo_shower)*i/(n-1); + p->gamma_pdf[i] = gamma_pdf(p->x_shower[i],p->pos_a,p->pos_b); + } + } + + if (p->delta_ray_photons) { + /* Calculate the normalization constant for the angular distribution. */ + norm = electron_get_angular_pdf_norm(p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o); + + for (i = 0; i < n; i++) { + p->cos_theta[i] = -1.0 + 2.0*i/(n-1); + /* 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. */ + pdf = electron_get_angular_pdf_no_norm(p->cos_theta[i], p->delta_ray_a, p->delta_ray_b, 1/avg_index_d2o)/norm; + if (i > 0) + p->cdf_delta[i] = p->cdf_delta[i-1] + (pdf + pdf_last)*(p->cos_theta[i]-p->cos_theta[i-1])/2.0; + else + p->cdf_delta[i] = 0.0; + pdf_last = pdf; + } + + gsl_spline_init(p->spline_delta, p->cdf_delta, p->cos_theta, n); + } + return p; } |