diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 39 |
1 files changed, 27 insertions, 12 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index fcc8a5f..eb2eb94 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -171,7 +171,7 @@ void particle_free(particle *p) static double get_expected_charge_shower(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 alpha, double beta, double rad) { - double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; + double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter; SUB(pmt_dir,pmt_pos,pos); @@ -195,19 +195,23 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p 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(); + scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); + scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI); - *reflected = f_reflec*charge; + scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); - return f*charge; + *reflected = (f_reflec + 1.0 - scatter)*charge; + + return f*charge*scatter; } 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, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; + 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, scattering_length_h2o, scattering_length_d2o, scatter; z = 1.0; @@ -242,23 +246,34 @@ static double get_expected_charge(double x, double beta, double theta0, double * 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(); + scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); + scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; 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; + scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); - return f*charge; + *reflected = (f_reflec + 1.0 - scatter)*charge; + + return f*charge*scatter; } -double F(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_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the CDF for the time distribution of photons at time `t`. */ size_t i; double p, mu_total; - p = mu_noise*t/GTVALID + mu_indirect*(pow(sigma,2)*norm(tmean,t,sigma) + (t-tmean)*norm_cdf(t,tmean,sigma))/(GTVALID-tmean); + p = mu_noise*t/GTVALID; + + if (t < tmean) + ; + else if (t < tmean + PSUP_REFLECTION_TIME) + p += mu_indirect*(t-tmean)/PSUP_REFLECTION_TIME; + else + p += mu_indirect; mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { @@ -275,7 +290,7 @@ double F(double t, double mu_noise, double mu_indirect, double *mu_direct, doubl return p/mu_total; } -double f(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_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the probability that a photon is detected at time `t`. * @@ -294,7 +309,7 @@ double f(double t, double mu_noise, double mu_indirect, double *mu_direct, doubl size_t i; double p, mu_total; - p = mu_noise/GTVALID + mu_indirect*norm_cdf(t,tmean,sigma)/(GTVALID-tmean); + 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); mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { @@ -324,9 +339,9 @@ 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(f(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_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); else - return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(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_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)); } 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, double rad) |