aboutsummaryrefslogtreecommitdiff
path: root/src/likelihood.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/likelihood.c')
-rw-r--r--src/likelihood.c108
1 files changed, 40 insertions, 68 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 1aa1653..ec58d11 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -22,11 +22,6 @@
#include "proton.h"
#include "id_particles.h"
-typedef struct intParams {
- path *p;
- int i;
-} intParams;
-
particle *particle_init(int id, double T0, size_t n)
{
/* Returns a struct with arrays of the particle position and kinetic
@@ -160,9 +155,9 @@ void particle_free(particle *p)
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)
+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, wavelength0, omega, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_h2o, l_d2o, l_acrylic;
+ 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;
z = 1.0;
@@ -182,31 +177,29 @@ double get_expected_charge(double x, double beta, double theta0, double *pos, do
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);
+ n = n_d2o;
} else {
- n = get_index_snoman_h2o(wavelength0);
+ n = n_h2o;
}
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));
+ 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();
- 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);
+ charge = 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);
+
+ *reflected = f_reflec*charge;
+
+ 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)
@@ -271,43 +264,41 @@ 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_time(double x, void *params)
+static void integrate_path(path *p, int pmt, double a, double b, size_t n, double *mu_direct, double *mu_indirect, double *time)
{
- intParams *pars = (intParams *) params;
- double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o;
-
- path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
+ size_t i;
+ double *qs, *q_indirects, *ts;
+ double dir[3], pos[3], n_d2o, n_h2o, wavelength0, t, theta0, beta, l_d2o, l_h2o, x;
- get_path_length(pos,pmts[pars->i].pos,AV_RADIUS,&l_d2o,&l_h2o);
+ qs = malloc(sizeof(double)*n);
+ q_indirects = malloc(sizeof(double)*n);
+ ts = malloc(sizeof(double)*n);
/* 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);
- t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
+ for (i = 0; i < n; i++) {
+ x = a + i*(b-a)/(n-1);
- return t*get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
-}
+ path_eval(p, x, pos, dir, &beta, &t, &theta0);
-static double gsl_reflected_charge(double x, void *params)
-{
- intParams *pars = (intParams *) params;
- double dir[3], pos[3], t, theta0, beta;
+ get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
- path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
+ t += l_d2o*n_d2o/SPEED_OF_LIGHT + l_h2o*n_h2o/SPEED_OF_LIGHT;
- return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 1);
-}
-
-static double gsl_charge(double x, void *params)
-{
- intParams *pars = (intParams *) params;
- double dir[3], pos[3], t, theta0, beta;
+ qs[i] = get_expected_charge(x, beta, theta0, pos, dir, pmts[pmt].pos, pmts[pmt].normal, PMT_RADIUS, q_indirects+i, n_d2o, n_h2o, l_d2o, l_h2o);
+ ts[i] = t*qs[i];
+ }
- path_eval(pars->p, x, pos, dir, &beta, &t, &theta0);
+ *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);
- return get_expected_charge(x, beta, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS, 0);
+ free(qs);
+ free(q_indirects);
+ free(ts);
}
double get_total_charge_approx(double T0, double *pos, double *dir, particle *p, int i, double smax, double theta0, double *t, double *mu_reflected)
@@ -532,7 +523,6 @@ 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, 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;
particle *p;
double pmt_dir[3], R, cos_theta, theta, wavelength0, n_d2o, n_h2o, theta_cerenkov, s, l_h2o, l_d2o;
@@ -540,19 +530,14 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
double a, b;
double tmp[3];
double mu_reflected;
+ path *path;
double mu_direct[MAX_PMTS];
double ts[MAX_PMTS];
double mu[MAX_PMTS];
double mu_noise, mu_indirect;
- gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(100);
- double result, error;
-
- size_t nevals;
-
- gsl_function F;
- F.params = &params;
+ double result;
p = particle_init(id, T0, 10000);
range = p->range;
@@ -565,7 +550,7 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
/* FIXME: is this formula valid for muons? */
theta0 = get_scattering_rms(range/2,p0,beta0,1.0)/sqrt(range/2);
- params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
+ path = path_init(pos, dir, T0, range, theta0, getKineticEnergy, p, z1, z2, n, p->mass);
if (beta0 > BETA_MIN)
get_smax(p, BETA_MIN, range, &smax);
@@ -580,8 +565,6 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
* calculation. */
if (fast && !ev->pmt_hits[i].hit) continue;
- params.i = i;
-
if (fast) {
mu_direct[i] = get_total_charge_approx(T0, pos, dir, p, i, smax, theta0, &ts[i], &mu_reflected);
ts[i] += t0;
@@ -640,20 +623,11 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
if (a < 0.0) a = 0.0;
if (b > range) b = range;
- F.function = &gsl_charge;
-
- gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
- mu_direct[i] = result;
-
- F.function = &gsl_reflected_charge;
-
- gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
- mu_indirect += result;
-
- F.function = &gsl_time;
+ double q_indirect;
+ integrate_path(path, i, a, b, 100, mu_direct+i, &q_indirect, &result);
+ mu_indirect += q_indirect;
if (mu_direct[i] > 1e-9) {
- gsl_integration_cquad(&F, a, b, 0, epsrel, w, &result, &error, &nevals);
ts[i] = t0 + result/mu_direct[i];
} else {
/* If the expected number of PE is very small, our estimate of
@@ -678,12 +652,10 @@ double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t
}
}
- path_free(params.p);
+ path_free(path);
particle_free(p);
- gsl_integration_cquad_workspace_free(w);
-
mu_noise = DARK_RATE*GTVALID*1e-9;
mu_indirect *= CHARGE_FRACTION/10000.0;
generated by cgit (git 2.34.1) at 2025-05-07 00:09:05 +0000