aboutsummaryrefslogtreecommitdiff
path: root/src/likelihood.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/likelihood.c')
-rw-r--r--src/likelihood.c239
1 files changed, 195 insertions, 44 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 5f66a72..7da81bb 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -18,12 +18,164 @@
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
#include "pmt_response.h"
+#include "electron.h"
+#include "proton.h"
+#include "id_particles.h"
typedef struct intParams {
path *p;
int i;
} intParams;
+particle *particle_init(int id, double T0, double rho, size_t n)
+{
+ /* Returns a struct with arrays of the particle position and kinetic
+ * energy. This struct can then be passed to particle_get_energy() to
+ * interpolate the particle's kinetic energy at any point along the track.
+ * For example:
+ *
+ * particle *p = particle_init(IDP_MU_MINUS, 1000.0, HEAVY_WATER_DENSITY, 100);
+ * double T = particle_get_energy(x, p);
+ *
+ * To compute the kinetic energy as a function of distance we need to solve
+ * the following integral equation:
+ *
+ * T(x) = T0 - \int_0^x dT(T(x))/dx
+ *
+ * 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;
+
+ particle *p = malloc(sizeof(particle));
+ p->x = malloc(sizeof(double)*n);
+ p->T = malloc(sizeof(double)*n);
+ p->n = n;
+
+ p->x[0] = 0;
+ p->T[0] = T0;
+
+ switch (id) {
+ case IDP_E_MINUS:
+ case IDP_E_PLUS:
+ p->mass = ELECTRON_MASS;
+ p->range = electron_get_range(T0, rho);
+
+ /* Now we loop over the points along the track and assume dE/dx is
+ * constant between points. */
+ for (i = 1; i < n; i++) {
+ p->x[i] = p->range*i/(n-1);
+ dEdx = electron_get_dEdx(p->T[i-1], rho);
+ p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ }
+
+ break;
+ case IDP_MU_MINUS:
+ case IDP_MU_PLUS:
+ p->mass = MUON_MASS;
+ p->range = muon_get_range(T0, rho);
+
+ /* Now we loop over the points along the track and assume dE/dx is
+ * constant between points. */
+ for (i = 1; i < n; i++) {
+ p->x[i] = p->range*i/(n-1);
+ dEdx = muon_get_dEdx(p->T[i-1], rho);
+ p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ }
+
+ break;
+ case IDP_PROTON:
+ p->mass = PROTON_MASS;
+ p->range = proton_get_range(T0, rho);
+
+ /* Now we loop over the points along the track and assume dE/dx is
+ * constant between points. */
+ for (i = 1; i < n; i++) {
+ p->x[i] = p->range*i/(n-1);
+ dEdx = proton_get_dEdx(p->T[i-1], rho);
+ p->T[i] = p->T[i-1] - dEdx*(p->x[i]-p->x[i-1]);
+ }
+
+ break;
+ default:
+ fprintf(stderr, "unknown particle id %i\n", id);
+ exit(1);
+ }
+
+ return p;
+}
+
+double particle_get_energy(double x, particle *p)
+{
+ /* Returns the approximate kinetic energy of a particle in water after
+ * travelling `x` cm with an initial kinetic energy `T`.
+ *
+ * Return value is in MeV. */
+ double T;
+
+ T = interp1d(x,p->x,p->T,p->n);
+
+ if (T < 0) return 0;
+
+ return T;
+}
+
+void particle_free(particle *p)
+{
+ free(p->x);
+ free(p->T);
+ free(p);
+}
+
+double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected)
+{
+ double pmt_dir[3], cos_theta, n, wavelength0, omega, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
+
+ z = 1.0;
+
+ SUB(pmt_dir,pmt_pos,pos);
+
+ normalize(pmt_dir);
+
+ cos_theta_pmt = DOT(pmt_dir,pmt_normal);
+
+ if (cos_theta_pmt > 0) return 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_approx(pos,pmt_pos,pmt_normal,r);
+
+ R = NORM(pos);
+
+ /* FIXME: I just calculate delta assuming 400 nm light. */
+ wavelength0 = 400.0;
+
+ if (R <= AV_RADIUS) {
+ n = get_index_snoman_d2o(wavelength0);
+ } else {
+ n = get_index_snoman_h2o(wavelength0);
+ }
+
+ if (beta < 1/n) return 0;
+
+ if (reflected)
+ f = get_weighted_pmt_reflectivity(acos(-cos_theta_pmt));
+ else
+ f = get_weighted_pmt_response(acos(-cos_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();
+
+ get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+ l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER;
+
+ return f*exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
+}
+
double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
{
/* Returns the CDF for the time distribution of photons at time `t`. */
@@ -86,12 +238,12 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
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));
}
-static double gsl_muon_time(double x, void *params)
+static double gsl_time(double x, void *params)
{
intParams *pars = (intParams *) params;
- double dir[3], pos[3], n_d2o, n_h2o, wavelength0, T, t, theta0, l_d2o, l_h2o;
+ double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o;
- path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+ path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
get_path_length(pos,pmts[pars->i].pos,AV_RADIUS,&l_d2o,&l_h2o);
@@ -102,30 +254,30 @@ static double gsl_muon_time(double x, void *params)
t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
- return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+ return t*get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
}
-static double gsl_muon_reflected_charge(double x, void *params)
+static double gsl_reflected_charge(double x, void *params)
{
intParams *pars = (intParams *) params;
- double dir[3], pos[3], T, t, theta0;
+ double dir[3], pos[3], t, theta0, beta;
- path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+ path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
- return get_expected_reflected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+ return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 1);
}
-static double gsl_muon_charge(double x, void *params)
+static double gsl_charge(double x, void *params)
{
intParams *pars = (intParams *) params;
- double dir[3], pos[3], T, t, theta0;
+ double dir[3], pos[3], t, theta0, beta;
- path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+ path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
- return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
+ return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
}
-double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy *m, int i, double smax, double theta0, double *t, double *mu_reflected)
+double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected)
{
/* Returns the approximate expected number of photons seen by PMT `i` using
* an analytic formula.
@@ -156,11 +308,11 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
*
* `smax` is currently calculated as the point where the particle velocity
* drops to 0.8 times the speed of light. */
- double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, p, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected;
+ double pmt_dir[3], tmp[3], R, cos_theta, x, z, s, a, b, beta, E, mom, T, omega, cos_theta_cerenkov, sin_theta_cerenkov, n_d2o, n_h2o, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic, wavelength0, f_reflected;
/* Calculate beta at the start of the track. */
- E0 = T0 + MUON_MASS;
- p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+ E0 = T0 + p->mass;
+ p0 = sqrt(E0*E0 - p->mass*p->mass);
beta0 = p0/E0;
/* First, we find the point along the track where the PMT is at the
@@ -219,12 +371,12 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
b = (s-z)/a;
/* Compute the kinetic energy at the point `s` along the track. */
- T = muon_get_energy(s,m);
+ T = particle_get_energy(s,p);
/* Calculate the particle velocity at the point `s`. */
- E = T + MUON_MASS;
- p = sqrt(E*E - MUON_MASS*MUON_MASS);
- beta = p/E;
+ E = T + p->mass;
+ mom = sqrt(E*E - p->mass*p->mass);
+ beta = mom/E;
if (beta < 1/n_d2o) return 0.0;
@@ -275,7 +427,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy
}
typedef struct betaRootParams {
- muon_energy *m;
+ particle *p;
double beta_min;
} betaRootParams;
@@ -283,22 +435,22 @@ static double beta_root(double x, void *params)
{
/* Function used to find at what point along a track a particle crosses
* some threshold in beta. */
- double T, E, p, beta;
+ double T, E, mom, beta;
betaRootParams *pars;
pars = (betaRootParams *) params;
- T = muon_get_energy(x, pars->m);
+ T = particle_get_energy(x, pars->p);
/* Calculate total energy */
- E = T + MUON_MASS;
- p = sqrt(E*E - MUON_MASS*MUON_MASS);
- beta = p/E;
+ E = T + pars->p->mass;
+ mom = sqrt(E*E - pars->p->mass*pars->p->mass);
+ beta = mom/E;
return beta - pars->beta_min;
}
-static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
+static int get_smax(particle *p, double beta_min, double range, double *smax)
{
/* Find the point along the track at which the particle's velocity drops to
* `beta_min`. */
@@ -311,7 +463,7 @@ static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
s = gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
- pars.m = m;
+ pars.p = p;
pars.beta_min = beta_min;
F.function = &beta_root;
@@ -339,17 +491,17 @@ static int get_smax(muon_energy *m, double beta_min, double range, double *smax)
return status;
}
-double getKineticEnergy(double x, void *params)
+static double getKineticEnergy(double x, void *p)
{
- return muon_get_energy(x,(muon_energy *) params);
+ return particle_get_energy(x, (particle *) p);
}
-double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
+double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double epsrel, int fast)
{
size_t i, j, nhit;
intParams params;
double logp[MAX_PE], nll[MAX_PMTS], range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
- muon_energy *m;
+ particle *p;
double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o;
double logp_path;
double a, b;
@@ -369,22 +521,21 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
gsl_function F;
F.params = &params;
- range = get_range(T0, HEAVY_WATER_DENSITY);
+ p = particle_init(id, T0, HEAVY_WATER_DENSITY, 10000);
+ range = p->range;
/* Calculate total energy */
- E0 = T0 + MUON_MASS;
- p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
+ E0 = T0 + p->mass;
+ p0 = sqrt(E0*E0 - p->mass*p->mass);
beta0 = p0/E0;
/* FIXME: is this formula valid for muons? */
theta0 = get_scattering_rms(range/2,p0,beta0,1.0,HEAVY_WATER_DENSITY)/sqrt(range/2);
- m = muon_init_energy(T0,HEAVY_WATER_DENSITY,10000);
-
- params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, m, z1, z2, n, MUON_MASS);
+ params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
if (beta0 > BETA_MIN)
- get_smax(m, BETA_MIN, range, &smax);
+ get_smax(p, BETA_MIN, range, &smax);
else
smax = 0.0;
@@ -399,7 +550,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
params.i = i;
if (fast) {
- mu_direct[i] = get_total_charge_approx(T0, pos, dir, m, i, smax, theta0, &ts[i], &mu_reflected);
+ mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected);
ts[i] += t0;
mu_indirect += mu_reflected;
} else {
@@ -456,17 +607,17 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
if (a < 0.0) a = 0.0;
if (b > range) b = range;
- F.function = &gsl_muon_charge;
+ F.function = &gsl_charge;
gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
mu_direct[i] = result;
- F.function = &gsl_muon_reflected_charge;
+ F.function = &gsl_reflected_charge;
gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
mu_indirect += result;
- F.function = &gsl_muon_time;
+ F.function = &gsl_time;
if (mu_direct[i] > 1e-9) {
gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
@@ -496,7 +647,7 @@ double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, doubl
path_free(params.p);
- muon_free_energy(m);
+ particle_free(p);
gsl_integration_cquad_workspace_free(w);
generated by cgit (git 2.34.1) at 2025-05-07 05:10:19 +0000