aboutsummaryrefslogtreecommitdiff
path: root/src/likelihood.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/likelihood.c')
-rw-r--r--src/likelihood.c269
1 files changed, 193 insertions, 76 deletions
diff --git a/src/likelihood.c b/src/likelihood.c
index 0228bd8..1919a53 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -208,7 +208,7 @@ void particle_free(particle *p)
free(p);
}
-static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *q_delta_ray, double *q_indirect_delta_ray, double *t)
+static void get_expected_charge_delta_ray(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
{
double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
@@ -221,8 +221,6 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
*t = 0.0;
*q = 0.0;
*reflected = 0.0;
- *q_delta_ray = 0.0;
- *q_indirect_delta_ray = 0.0;
if (cos_theta_pmt > 0) return;
/* Calculate the cosine of the angle between the track direction and the
@@ -247,17 +245,55 @@ static void get_expected_charge_shower(particle *p, double *pos, double *dir, in
* 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. */
- charge = 0.0;
- if (p->shower_photons > 0)
- charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_d2o);
-
- if (p->delta_ray_photons > 0)
- *q_delta_ray = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/avg_index_d2o);
+ charge = constant*p->delta_ray_photons*electron_get_angular_pdf_delta_ray(cos_theta,p->delta_ray_a,p->delta_ray_b,1/avg_index_d2o);
*reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
- *q_indirect_delta_ray = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*(*q_delta_ray) + prob_sct*(*q_delta_ray);
- *q_delta_ray *= (1.0-prob_abs)*(1.0-prob_sct)*f;
+ *t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
+
+ *q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
+}
+
+static void get_expected_charge_shower(particle *p, double *pos, double *dir, int pmt, double *q, double *reflected, double *t)
+{
+ double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, charge, constant, prob_abs, prob_sct, l_d2o, l_h2o;
+
+ SUB(pmt_dir,pmts[pmt].pos,pos);
+
+ normalize(pmt_dir);
+
+ cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal);
+
+ *t = 0.0;
+ *q = 0.0;
+ *reflected = 0.0;
+ if (cos_theta_pmt > 0) return;
+
+ /* 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_fast(pos,pmts[pmt].pos,pmts[pmt].normal,PMT_RADIUS);
+
+ f_reflec = get_weighted_pmt_reflectivity(-cos_theta_pmt);
+ f = get_weighted_pmt_response(-cos_theta_pmt);
+
+ get_path_length(pos,pmts[pmt].pos,AV_RADIUS,&l_d2o,&l_h2o);
+
+ prob_abs = 1.0 - get_fabs_d2o(l_d2o)*get_fabs_h2o(l_h2o)*get_fabs_acrylic(AV_THICKNESS);
+ prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o);
+
+ constant = get_weighted_quantum_efficiency()*omega/(2*M_PI);
+
+ /* 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. */
+ charge = constant*p->shower_photons*electron_get_angular_pdf(cos_theta,p->a,p->b,1/avg_index_d2o);
+
+ *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge;
*t = (l_d2o*avg_index_d2o + l_h2o*avg_index_h2o)/SPEED_OF_LIGHT;
@@ -327,7 +363,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double
*q = (1.0-prob_abs)*(1.0-prob_sct)*f*charge;
}
-double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
+double time_cdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma)
{
/* Returns the CDF for the time distribution of photons at time `t`. */
size_t i;
@@ -344,20 +380,15 @@ double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct
mu_total = mu_noise + mu_indirect;
for (i = 0; i < n; i++) {
- if (mu_direct[i] == 0.0) continue;
- p += mu_direct[i]*norm_cdf(t,ts[i],sigma);
- mu_total += mu_direct[i];
- }
- for (i = 0; i < n; i++) {
- if (mu_shower[i] == 0.0) continue;
- p += mu_shower[i]*norm_cdf(t,ts_shower[i],ts_sigma[i]);
- mu_total += mu_shower[i];
+ if (mu[i] == 0.0) continue;
+ p += mu[i]*norm_cdf(t,ts[i],ts_sigma[i]);
+ mu_total += mu[i];
}
return p/mu_total;
}
-double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
+double time_pdf(double t, double mu_noise, double mu_indirect, double *mu, size_t n, double *ts, double tmean, double *ts_sigma)
{
/* Returns the probability that a photon is detected at time `t`.
*
@@ -376,24 +407,19 @@ double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct
size_t i;
double p, mu_total;
- p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*sigma))-erf((t-tmean-PSUP_REFLECTION_TIME)/(sqrt(2)*sigma)))/(2*PSUP_REFLECTION_TIME);
+ p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*PMT_TTS))-erf((t-tmean-PSUP_REFLECTION_TIME)/(sqrt(2)*PMT_TTS)))/(2*PSUP_REFLECTION_TIME);
mu_total = mu_noise + mu_indirect;
for (i = 0; i < n; i++) {
- if (mu_direct[i] == 0.0) continue;
- p += mu_direct[i]*norm(t,ts[i],sigma);
- mu_total += mu_direct[i];
- }
- for (i = 0; i < n; i++) {
- if (mu_shower[i] == 0.0) continue;
- p += mu_shower[i]*norm(t,ts_shower[i],ts_sigma[i]);
- mu_total += mu_shower[i];
+ if (mu[i] == 0.0) continue;
+ p += mu[i]*norm(t,ts[i],ts_sigma[i]);
+ mu_total += mu[i];
}
return p/mu_total;
}
-double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n2, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma)
+double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *mu, size_t n2, double *ts, double tmean, double *ts_sigma)
{
/* Returns the first order statistic for observing a PMT hit at time `t`
* given `n` hits.
@@ -406,15 +432,71 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m
* convolve first which is what we do here. In addition, the problem is not
* analytically tractable if you do things the other way around. */
if (n == 1)
- return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
+ return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma));
else
- return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma));
+ return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu,n2,ts,tmean,ts_sigma));
+}
+
+static void integrate_path_delta_ray(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
+{
+ size_t i;
+ double pos[3], t, q, r, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
+
+ q_sum = 0.0;
+ r_sum = 0.0;
+ t_sum = 0.0;
+ t2_sum = 0.0;
+ for (i = 0; i < n; i++) {
+ pos[0] = pos0[0] + x[i]*dir0[0];
+ pos[1] = pos0[1] + x[i]*dir0[1];
+ pos[2] = pos0[2] + x[i]*dir0[2];
+
+ get_expected_charge_delta_ray(p, pos, dir0, pmt, &q, &r, &t);
+
+ q = q*pdf[i];
+ r = r*pdf[i];
+ t += x[i]/SPEED_OF_LIGHT;
+
+ if (i == 0 || i == (n - 1)) {
+ q_sum += q;
+ r_sum += r;
+ t_sum += t*q;
+ t2_sum += t*t*q;
+ } else {
+ q_sum += 2*q;
+ r_sum += 2*r;
+ t_sum += 2*t*q;
+ t2_sum += 2*t*t*q;
+ }
+ }
+
+ dx = x[1] - x[0];
+ *mu_direct = q_sum*dx*0.5;
+ *mu_indirect = r_sum*dx*0.5;
+
+ if (*mu_direct == 0.0) {
+ *time = 0.0;
+ *sigma = PMT_TTS;
+ } else {
+ *time = t_sum*dx*0.5/(*mu_direct);
+
+ /* Variance in the time = E(t^2) - E(t)^2. */
+ tmp = t2_sum*dx*0.5/(*mu_direct) - (*time)*(*time);
+
+ if (tmp >= 0) {
+ *sigma = sqrt(tmp + PMT_TTS*PMT_TTS);
+ } else {
+ /* This should never happen but does occasionally due to floating
+ * point rounding error. */
+ *sigma = PMT_TTS;
+ }
+ }
}
static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma)
{
size_t i;
- double pos[3], t, q, r, qd, rd, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
+ double pos[3], t, q, r, dx, tmp, q_sum, r_sum, t_sum, t2_sum;
q_sum = 0.0;
r_sum = 0.0;
@@ -425,10 +507,10 @@ static void integrate_path_shower(particle *p, double *x, double *pdf, double T0
pos[1] = pos0[1] + x[i]*dir0[1];
pos[2] = pos0[2] + x[i]*dir0[2];
- get_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &qd, &rd, &t);
+ get_expected_charge_shower(p, pos, dir0, pmt, &q, &r, &t);
- q = q*pdf[i] + qd/p->range;
- r = r*pdf[i] + rd/p->range;
+ q = q*pdf[i];
+ r = r*pdf[i];
t += x[i]/SPEED_OF_LIGHT;
if (i == 0 || i == (n - 1)) {
@@ -803,7 +885,7 @@ double nll_best(event *ev)
* hypothesis and this ideal case. See Chapter 9 in Jayne's "Probability
* Theory: The Logic of Science" for more details. */
size_t i, j, nhit, maxj;
- static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_shower, ts_sigma, mu_shower;
+ static double logp[MAX_PE], nll[MAX_PMTS], mu[MAX_PMTS], ts[MAX_PMTS], ts_sigma;
double log_mu, max_logp, mu_noise, mu_indirect_total, min_ratio;
mu_noise = DARK_RATE*GTVALID*1e-9;
@@ -840,8 +922,6 @@ double nll_best(event *ev)
ts[i] = ev->pmt_hits[i].t;
}
- mu_shower = 0;
- ts_shower = 0;
ts_sigma = PMT_TTS;
mu_indirect_total = 0.0;
@@ -861,7 +941,7 @@ double nll_best(event *ev)
if (ev->pmt_hits[i].hit) {
for (j = 1; j < MAX_PE; j++) {
- logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], &mu_shower, 1, &ts[i], &ts_shower, ts[i], PMT_TTS, &ts_sigma);
+ logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j) + log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i], 1, &ts[i], ts[i], &ts_sigma);
if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
@@ -897,7 +977,7 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
* `z1` and `z2` should be arrays of length `n` and represent coefficients
* in the Karhunen Loeve expansion of the track path. These are only used
* when `fast` is zero. */
- size_t i, j, nhit;
+ size_t i, j, k, nhit;
static double logp[MAX_PE], nll[MAX_PMTS];
double range, theta0, E0, p0, beta0, smax, log_mu, max_logp;
particle *p;
@@ -906,15 +986,29 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
double mu_reflected;
path *path;
- static double mu_direct[MAX_PMTS][MAX_VERTICES];
- static double mu_shower[MAX_PMTS][MAX_VERTICES];
- static double ts[MAX_PMTS][MAX_VERTICES];
- static double ts_shower[MAX_PMTS][MAX_VERTICES];
- static double ts_sigma[MAX_PMTS][MAX_VERTICES];
- static double mu[MAX_PMTS];
+ /* Array representing the expected number of photons at each PMT for each
+ * vertex. The last axis represents the contributions from direct, shower,
+ * and delta ray photons respectively, i.e.
+ *
+ * mu[0][2][0] - expected number of direct photons from vertex 2 at PMT 0
+ * mu[0][2][1] - expected number of shower photons from vertex 2 at PMT 0
+ * mu[0][2][2] - expected number of delta ray photons from vertex 2 at PMT 0
+ * ...
+ * etc.
+ *
+ */
+ static double mu[MAX_PMTS][MAX_VERTICES][3];
+ /* Array representing the average time of arrival at each PMT for each
+ * vertex. See above for the how the array is indexed. */
+ static double ts[MAX_PMTS][MAX_VERTICES][3];
+ /* Array representing the standard deviation of the time at each PMT for
+ * each vertex. See above for the how the array is indexed. */
+ static double ts_sigma[MAX_PMTS][MAX_VERTICES][3];
+ static double mu_sum[MAX_PMTS];
double mu_noise, mu_indirect[MAX_VERTICES], mu_indirect_total;
- double *x, *pdf;
+ static double x[MAX_NPOINTS], pdf[MAX_NPOINTS];
+ static double x_delta[MAX_NPOINTS], pdf_delta[MAX_NPOINTS];
double result;
@@ -935,13 +1029,13 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
/* Initialize static arrays. */
for (i = 0; i < MAX_PMTS; i++) {
- mu[i] = 0.0;
+ mu_sum[i] = 0.0;
for (j = 0; j < n; j++) {
- mu_direct[i][j] = 0.0;
- mu_shower[i][j] = 0.0;
- ts[i][j] = 0.0;
- ts_shower[i][j] = 0.0;
- ts_sigma[i][j] = PMT_TTS;
+ for (k = 0; k < 3; k++) {
+ mu[i][j][k] = 0.0;
+ ts[i][j][k] = 0.0;
+ ts_sigma[i][j][k] = PMT_TTS;
+ }
}
}
@@ -996,8 +1090,6 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
if (npoints_shower < MIN_NPOINTS) npoints_shower = MIN_NPOINTS;
if (npoints_shower > MAX_NPOINTS) npoints_shower = MAX_NPOINTS;
- x = malloc(sizeof(double)*npoints_shower);
- pdf = malloc(sizeof(double)*npoints_shower);
for (i = 0; i < npoints_shower; i++) {
x[i] = xlo + i*(xhi-xlo)/(npoints_shower-1);
pdf[i] = gamma_pdf(x[i],pos_a,pos_b);
@@ -1008,6 +1100,23 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
for (i = 0; i < npoints_shower; i++) {
pdf[i] /= norm;
}
+
+ for (i = 0; i < npoints_shower; i++) {
+ x_delta[i] = i*p->range/(npoints_shower-1);
+ /* For now, we assume the delta rays are produced uniformly
+ * along the path. This isn't quite correct, but the
+ * distribution is close to flat.
+ *
+ * FIXME: Figure out a good parameterization for the
+ * distribution of delta ray photons along the path. */
+ pdf_delta[i] = 1.0;
+ }
+
+ /* Make sure the PDF is normalized. */
+ norm = trapz(pdf_delta,x_delta[1]-x_delta[0],npoints_shower);
+ for (i = 0; i < npoints_shower; i++) {
+ pdf_delta[i] /= norm;
+ }
}
if (beta0 > BETA_MIN)
@@ -1031,18 +1140,18 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
if (hit_only && !ev->pmt_hits[i].hit) continue;
if (fast) {
- mu_direct[i][j] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j], &mu_reflected, cos_theta_cerenkov, sin_theta_cerenkov);
- ts[i][j] += v[j].t0;
+ mu[i][j][0] = get_total_charge_approx(beta0, v[j].pos, v[j].dir, p, i, smax, theta0, &ts[i][j][0], &mu_reflected, cos_theta_cerenkov, sin_theta_cerenkov);
+ ts[i][j][0] += v[j].t0;
mu_indirect[j] += mu_reflected;
continue;
}
double q_indirect;
- integrate_path(path, i, &mu_direct[i][j], &q_indirect, &result);
+ integrate_path(path, i, &mu[i][j][0], &q_indirect, &result);
mu_indirect[j] += q_indirect;
- if (mu_direct[i][j] > 1e-9) {
- ts[i][j] = v[j].t0 + result/mu_direct[i][j];
+ if (mu[i][j][0] > 1e-9) {
+ ts[i][j][0] = v[j].t0 + result/mu[i][j][0];
} else {
/* If the expected number of PE is very small, our estimate of
* the time can be crazy since the error in the total charge is
@@ -1054,19 +1163,24 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
* I should test it to see if we can get away with something even
* simpler (like just computing the distance from the start of the
* track to the PMT). */
- ts[i][j] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,cos_theta_cerenkov,sin_theta_cerenkov);
+ ts[i][j][0] = v[j].t0 + guess_time(v[j].pos,v[j].dir,pmts[i].pos,smax,cos_theta_cerenkov,sin_theta_cerenkov);
}
- integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu_shower[i][j],&q_indirect,&result,&ts_sigma[i][j]);
- mu_indirect[j] += q_indirect;
- ts_shower[i][j] = v[j].t0 + result;
+ if (p->shower_photons > 0) {
+ integrate_path_shower(p,x,pdf,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][1],&q_indirect,&result,&ts_sigma[i][j][1]);
+ mu_indirect[j] += q_indirect;
+ ts[i][j][1] = v[j].t0 + result;
+ }
+
+ if (p->delta_ray_photons > 0) {
+ integrate_path_delta_ray(p,x_delta,pdf_delta,v[j].T0,v[j].pos,v[j].dir,i,npoints_shower,&mu[i][j][2],&q_indirect,&result,&ts_sigma[i][j][2]);
+ mu_indirect[j] += q_indirect;
+ ts[i][j][2] = v[j].t0 + result;
+ }
}
if (!fast) {
path_free(path);
-
- free(x);
- free(pdf);
}
particle_free(p);
@@ -1089,8 +1203,11 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
* calculation. */
if (hit_only && !ev->pmt_hits[i].hit) continue;
+ mu_sum[i] = mu_indirect_total + mu_noise;
for (j = 0; j < n; j++) {
- mu[i] += mu_direct[i][j] + mu_shower[i][j] + mu_indirect_total + mu_noise;
+ for (k = 0; k < 3; k++) {
+ mu_sum[i] += mu[i][j][k];
+ }
}
}
@@ -1108,20 +1225,20 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
* calculation. */
if (hit_only && !ev->pmt_hits[i].hit) continue;
- log_mu = log(mu[i]);
+ log_mu = log(mu_sum[i]);
if (fast)
- log_pt_fast = log_pt(ev->pmt_hits[i].t, 1, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
+ log_pt_fast = log_pt(ev->pmt_hits[i].t, 1, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]);
if (ev->pmt_hits[i].hit) {
for (j = 1; j < MAX_PE; j++) {
- logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu[i] + j*log_mu - lnfact(j);
+ logp[j] = get_log_pq(ev->pmt_hits[i].qhs,j) - mu_sum[i] + j*log_mu - lnfact(j);
if (!charge_only) {
if (fast)
logp[j] += log_pt_fast;
else
- logp[j] += log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu_direct[i][0], &mu_shower[i][0], n, &ts[i][0], &ts_shower[i][0], ts[i][0], PMT_TTS, &ts_sigma[i][0]);
+ logp[j] += log_pt(ev->pmt_hits[i].t, j, mu_noise, mu_indirect_total, &mu[i][0][0], n*3, &ts[i][0][0], ts[i][0][0], &ts_sigma[i][0][0]);
}
if (j == 1 || logp[j] > max_logp) max_logp = logp[j];
@@ -1134,13 +1251,13 @@ double nll(event *ev, vertex *v, size_t n, double dx, double dx_shower, int fast
nll[nhit++] = -logsumexp(logp+1, j-1);
} else {
- logp[0] = -mu[i];
+ logp[0] = -mu_sum[i];
if (fast) {
nll[nhit++] = -logp[0];
continue;
}
for (j = 1; j < MAX_PE_NO_HIT; j++) {
- logp[j] = get_log_pmiss(j) - mu[i] + j*log_mu - lnfact(j);
+ logp[j] = get_log_pmiss(j) - mu_sum[i] + j*log_mu - lnfact(j);
}
nll[nhit++] = -logsumexp(logp, MAX_PE_NO_HIT);
}
generated by cgit (git 2.34.1) at 2025-05-06 21:36:44 +0000