aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/likelihood.c108
-rw-r--r--src/likelihood.h1
-rw-r--r--src/misc.c20
-rw-r--r--src/misc.h6
-rw-r--r--src/optics.c6
-rw-r--r--src/test.c31
6 files changed, 98 insertions, 74 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;
diff --git a/src/likelihood.h b/src/likelihood.h
index 71f77f0..8328322 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -44,7 +44,6 @@ typedef struct particle {
particle *particle_init(int id, double T0, size_t n);
double particle_get_energy(double x, particle *p);
void particle_free(particle *p);
-double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, int reflected);
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);
#endif
diff --git a/src/misc.c b/src/misc.c
index ba0a090..d20bc26 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -218,6 +218,26 @@ static struct {
{100,363.73937555556347},
};
+double trapz(const double *y, double dx, size_t n)
+{
+ /* Returns the integral of `y` using the trapezoidal rule assuming a
+ * constant grid spacing `dx`.
+ *
+ * See https://en.wikipedia.org/wiki/Trapezoidal_rule. */
+ size_t i;
+ double sum = 0.0;
+
+ if (n < 2) return 0.0;
+
+ sum = y[0];
+ for (i = 1; i < n-1; i++) {
+ sum += 2*y[i];
+ }
+ sum += y[n-1];
+
+ return sum*dx/2.0;
+}
+
void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2)
{
/* Returns the path length inside and outside a circle of radius `R` for a
diff --git a/src/misc.h b/src/misc.h
index dadbbce..80a580f 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -3,9 +3,15 @@
#include <stdlib.h> /* for size_t */
+/* Macro to compute the size of a static C array.
+ *
+ * See https://stackoverflow.com/questions/1598773. */
+#define LEN(x) ((sizeof(x)/sizeof(0[x]))/((size_t)(!(sizeof(x) % sizeof(0[x])))))
+
#define LN_MAX 100
#define LNFACT_MAX 100
+double trapz(const double *y, double dx, size_t n);
void get_path_length(double *pos1, double *pos2, double R, double *l1, double *l2);
double ln(unsigned int n);
double lnfact(unsigned int n);
diff --git a/src/optics.c b/src/optics.c
index 1286490..2d28ba4 100644
--- a/src/optics.c
+++ b/src/optics.c
@@ -6,11 +6,7 @@
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h> /* for gsl_strerror() */
#include <gsl/gsl_spline.h>
-
-/* Macro to compute the size of a static C array.
- *
- * See https://stackoverflow.com/questions/1598773. */
-#define LEN(x) ((sizeof(x)/sizeof(0[x]))/((size_t)(!(sizeof(x) % sizeof(0[x])))))
+#include "misc.h"
/* Absorption coefficients for H2O and D2O as a function of wavelength from
* SNOMAN.
diff --git a/src/test.c b/src/test.c
index c1da272..242e728 100644
--- a/src/test.c
+++ b/src/test.c
@@ -879,6 +879,36 @@ err:
return 1;
}
+int test_trapz(char *err)
+{
+ /* Tests the trapz() function. */
+ size_t i;
+ double y[100], integral, expected;
+
+ init_genrand(0);
+
+ for (i = 0; i < LEN(y); i++) {
+ y[i] = genrand_real2();
+ }
+
+ expected = 0.0;
+ for (i = 1; i < LEN(y); i++) {
+ expected += (y[i-1] + y[i])/2;
+ }
+
+ integral = trapz(y, 1.0, LEN(y));
+
+ if (!isclose(integral, expected, 0, 1e-9)) {
+ sprintf(err, "trapz() returned %.5g, but expected %.5g", integral, expected);
+ goto err;
+ }
+
+ return 0;
+
+err:
+ return 1;
+}
+
struct tests {
testFunction *test;
char *name;
@@ -906,6 +936,7 @@ struct tests {
{test_proton_get_range, "test_proton_get_range"},
{test_proton_get_energy, "test_proton_get_energy"},
{test_log_norm, "test_log_norm"},
+ {test_trapz, "test_trapz"},
};
int main(int argc, char **argv)
generated by cgit (git 2.34.1) at 2025-05-06 21:58:34 +0000