aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/Makefile2
-rw-r--r--src/fit.c30
-rw-r--r--src/likelihood.c92
-rw-r--r--src/likelihood.h3
-rw-r--r--src/muon.c14
-rw-r--r--src/path.c10
-rw-r--r--src/path.h2
-rw-r--r--src/test.c3
8 files changed, 82 insertions, 74 deletions
diff --git a/src/Makefile b/src/Makefile
index 1ea1685..75a75cb 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -24,7 +24,7 @@ test-charge: test-charge.o sno_charge.o misc.o
test-zebra: test-zebra.o zebra.o pack2b.o
-fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o
+fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o
clean:
rm -f *.o calculate_limits test Makefile.dep
diff --git a/src/fit.c b/src/fit.c
index 2a86456..7deec9d 100644
--- a/src/fit.c
+++ b/src/fit.c
@@ -20,6 +20,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
double pos[3], dir[3];
static size_t iter;
double fval;
+ double z1[1], z2[1];
T = x[0];
@@ -35,9 +36,12 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
t0 = x[6];
- fval = nll_muon((event *) params, T, pos, dir, t0);
+ z1[0] = x[7];
+ z2[0] = x[8];
- printf("%5zu %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f f() = %7.3e\n",
+ fval = nll_muon((event *) params, T, pos, dir, t0, z1, z2, 1);
+
+ printf("%5zu %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f %4.2f %4.2f f() = %7.3e\n",
iter++,
x[0],
x[1],
@@ -46,6 +50,8 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
x[4],
x[5],
x[6],
+ x[7],
+ x[8],
fval);
return fval;
@@ -53,20 +59,22 @@ double nll(unsigned int n, const double *x, double *grad, void *params)
int fit_event(event *ev, double *fmin)
{
- double x[7], ss[7], lb[7], ub[7];
+ double x[9], ss[9], lb[9], ub[9];
int rv;
- nlopt_opt opt = nlopt_create(NLOPT_LN_SBPLX, 7);
+ nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9);
nlopt_set_min_objective(opt,nll,ev);
nlopt_set_ftol_abs(opt, 1e-5);
- x[0] = 1000.0;
- x[1] = 100.0;
- x[2] = 100.0;
- x[3] = 100.0;
+ x[0] = 1200.0;
+ x[1] = 0.0;
+ x[2] = 0.0;
+ x[3] = 0.0;
x[4] = 1.57;
x[5] = 0.0;
x[6] = 130.0;
+ x[7] = 0.0;
+ x[8] = 0.0;
ss[0] = 10.0;
ss[1] = 10.0;
@@ -75,6 +83,8 @@ int fit_event(event *ev, double *fmin)
ss[4] = 0.1;
ss[5] = 0.1;
ss[6] = 1.0;
+ ss[7] = 0.1;
+ ss[8] = 0.1;
lb[0] = 1.0;
lb[1] = -1000.0;
@@ -83,6 +93,8 @@ int fit_event(event *ev, double *fmin)
lb[4] = 0.0;
lb[5] = 0.0;
lb[6] = 0.0;
+ lb[7] = -10.0;
+ lb[8] = -10.0;
ub[0] = 100000.0;
ub[1] = 1000.0;
@@ -91,6 +103,8 @@ int fit_event(event *ev, double *fmin)
ub[4] = M_PI;
ub[5] = 2*M_PI;
ub[6] = 400.0;
+ ub[7] = 10.0;
+ ub[8] = 10.0;
nlopt_set_lower_bounds(opt, lb);
nlopt_set_upper_bounds(opt, ub);
diff --git a/src/likelihood.c b/src/likelihood.c
index a4864a8..f06fade 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -11,6 +11,13 @@
#include "optics.h"
#include "sno_charge.h"
#include "pdg.h"
+#include "path.h"
+#include <stddef.h> /* for size_t */
+
+typedef struct intParams {
+ path *p;
+ int i;
+} intParams;
double F(double t, double mu_noise, double mu_indirect, double *mu_direct, size_t n, double *ts, double tmean, double sigma)
{
@@ -76,28 +83,16 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
static double gsl_muon_time(double x, void *params)
{
- double *params2 = (double *) params;
- double T0 = params2[0];
- double pos0[3], dir[3], pos[3], pmt_dir[3], R, n, wavelength0;
- int i;
- double t;
- i = (int) params2[1];
- pos0[0] = params2[2];
- pos0[1] = params2[3];
- pos0[2] = params2[4];
- dir[0] = params2[5];
- dir[1] = params2[6];
- dir[2] = params2[7];
-
- pos[0] = pos0[0] + dir[0]*x;
- pos[1] = pos0[1] + dir[1]*x;
- pos[2] = pos0[2] + dir[2]*x;
+ intParams *pars = (intParams *) params;
+ double dir[3], pos[3], pmt_dir[3], R, n, wavelength0, T, t, theta0, distance;
+
+ path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
R = NORM(pos);
- SUB(pmt_dir,pmts[i].pos,pos);
+ SUB(pmt_dir,pmts[pars->i].pos,pos);
- double distance = NORM(pmt_dir);
+ distance = NORM(pmt_dir);
/* FIXME: I just calculate delta assuming 400 nm light. */
wavelength0 = 400.0;
@@ -107,38 +102,32 @@ static double gsl_muon_time(double x, void *params)
n = get_index_snoman_h2o(wavelength0);
}
- t = x/SPEED_OF_LIGHT + distance*n/SPEED_OF_LIGHT;
+ t += distance*n/SPEED_OF_LIGHT;
- return t*get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
+ return t*get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
}
static double gsl_muon_charge(double x, void *params)
{
- double *params2 = (double *) params;
- double T0 = params2[0];
- double pos0[3], dir[3], pos[3];
- int i;
- i = (int) params2[1];
- pos0[0] = params2[2];
- pos0[1] = params2[3];
- pos0[2] = params2[4];
- dir[0] = params2[5];
- dir[1] = params2[6];
- dir[2] = params2[7];
-
- pos[0] = pos0[0] + dir[0]*x;
- pos[1] = pos0[1] + dir[1]*x;
- pos[2] = pos0[2] + dir[2]*x;
-
- return get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS);
+ intParams *pars = (intParams *) params;
+ double dir[3], pos[3], T, t, theta0;
+
+ path_eval(pars->p, x, pos, dir, &T, &t, &theta0);
+
+ return get_expected_charge(x, T, theta0, pos, dir, pmts[pars->i].pos, pmts[pars->i].normal, PMT_RADIUS);
}
-double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
+double getKineticEnergy(double x, double T0)
+{
+ return get_T(T0, x, HEAVY_WATER_DENSITY);
+}
+
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n)
{
size_t i, j;
- double params[8];
+ intParams params;
double total_charge;
- double logp[MAX_PE], nll, range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov;
+ double logp[MAX_PE], nll, range, pmt_dir[3], R, x, cos_theta, theta, theta_cerenkov, theta0, E0, p0, beta0;
double tmean = 0.0;
int npmt = 0;
@@ -155,21 +144,24 @@ double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
gsl_function F;
F.params = &params;
- range = get_range(T, HEAVY_WATER_DENSITY);
+ range = get_range(T0, HEAVY_WATER_DENSITY);
+
+ /* Calculate total energy */
+ E0 = T0 + MUON_MASS;
+ p0 = sqrt(E0*E0 - MUON_MASS*MUON_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);
+
+ params.p = path_init(pos, dir, T0, range, theta0, getKineticEnergy, z1, z2, n, MUON_MASS);
total_charge = 0.0;
npmt = 0;
for (i = 0; i < MAX_PMTS; i++) {
if (ev->pmt_hits[i].flags || (pmts[i].pmt_type != PMT_NORMAL && pmts[i].pmt_type != PMT_OWL)) continue;
- params[0] = T;
- params[1] = i;
- params[2] = pos[0];
- params[3] = pos[1];
- params[4] = pos[2];
- params[5] = dir[0];
- params[6] = dir[1];
- params[7] = dir[2];
+ params.i = i;
/* First, we try to compute the distance along the track where the
* PMT is at the Cerenkov angle. The reason for this is because for
@@ -245,6 +237,8 @@ double nll_muon(event *ev, double T, double *pos, double *dir, double t0)
}
}
+ path_free(params.p);
+
tmean /= npmt;
gsl_integration_cquad_workspace_free(w);
diff --git a/src/likelihood.h b/src/likelihood.h
index 43b0823..9267ef8 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -2,6 +2,7 @@
#define LIKELIHOOD_H
#include "event.h"
+#include <stddef.h> /* for size_t */
#define MAX_PE 100
#define CHARGE_FRACTION 60000.0
@@ -11,6 +12,6 @@
/* Event window (ns) */
#define GTVALID 400.0
-double nll_muon(event *ev, double T, double *pos, double *dir, double t0);
+double nll_muon(event *ev, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n);
#endif
diff --git a/src/muon.c b/src/muon.c
index 7012214..8db6237 100644
--- a/src/muon.c
+++ b/src/muon.c
@@ -216,9 +216,9 @@ double get_dEdx(double T, double rho)
return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho;
}
-double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
+double get_expected_charge(double x, double T, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r)
{
- double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R, E0, p0, beta0;
+ double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R;
z = 1.0;
@@ -244,23 +244,13 @@ double get_expected_charge(double x, double T, double T0, double *pos, double *d
wavelength0 = 400.0;
if (R <= AV_RADIUS) {
- rho = HEAVY_WATER_DENSITY;
n = get_index_snoman_d2o(wavelength0);
} else {
- rho = WATER_DENSITY;
n = get_index_snoman_h2o(wavelength0);
}
if (beta < 1/n) return 0;
- /* Calculate total energy */
- E0 = T0 + MUON_MASS;
- p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS);
- beta0 = p0/E0;
-
- /* FIXME: is this formula valid for muons? */
- theta0 = get_scattering_rms(x,p0,beta0,z,rho);
-
/* FIXME: add angular response and scattering/absorption. */
return omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0);
}
diff --git a/src/path.c b/src/path.c
index f09ec45..dff1dc7 100644
--- a/src/path.c
+++ b/src/path.c
@@ -36,6 +36,8 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
path *p = malloc(sizeof(path));
+ p->theta0 = theta0;
+
p->pos[0] = pos[0];
p->pos[1] = pos[1];
p->pos[2] = pos[2];
@@ -114,10 +116,13 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0
return p;
}
-void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t)
+void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0)
{
double normal[3], k[3], tmp[3], phi;
+ if (s > p->spline[0]->x[p->spline[0]->size-1])
+ s = p->spline[0]->x[p->spline[0]->size-1];
+
tmp[0] = gsl_spline_eval(p->spline[0],s,p->acc[0]);
tmp[1] = gsl_spline_eval(p->spline[1],s,p->acc[1]);
tmp[2] = gsl_spline_eval(p->spline[2],s,p->acc[2]);
@@ -146,6 +151,9 @@ void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t
normalize(tmp);
rotate(dir,tmp,normal,phi);
+
+ /* FIXME: This should be the *residual* scattering RMS. */
+ *theta0 = p->theta0*sqrt(s);
}
void path_free(path *p)
diff --git a/src/path.h b/src/path.h
index 8aaa024..6e73d10 100644
--- a/src/path.h
+++ b/src/path.h
@@ -15,7 +15,7 @@ typedef struct path {
double path_get_coefficient(unsigned int k, double *s, double *x, double theta0, size_t n);
path *path_init(double *pos, double *dir, double T0, double range, double theta0, getKineticEnergyFunc *fun, double *z1, double *z2, size_t n, double m);
-void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t);
+void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t, double *theta0);
void path_free(path *p);
#endif
diff --git a/src/test.c b/src/test.c
index 3ef6531..baf0914 100644
--- a/src/test.c
+++ b/src/test.c
@@ -342,6 +342,7 @@ int test_path(char *err)
double pos[3], dir[3];
double E, mom, beta;
double s;
+ double theta0;
T0 = 1.0;
range = 1.0;
@@ -371,7 +372,7 @@ int test_path(char *err)
t_expected = s/(beta*SPEED_OF_LIGHT);
- path_eval(p,s,pos,dir,&T,&t);
+ path_eval(p,s,pos,dir,&T,&t,&theta0);
for (k = 0; k < 3; k++) {
if (!isclose(pos[k], pos_expected[k], 1e-9, 1e-9)) {
generated by cgit (git 2.34.1) at 2025-05-07 12:24:21 +0000