diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 18 | ||||
| -rw-r--r-- | src/scattering.c | 6 |
2 files changed, 12 insertions, 12 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 17cde5d..950b646 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -547,18 +547,18 @@ double get_theta0_min(double R, double r, double sin_theta_pmt) static void get_expected_charge(double beta, double theta0, double *pos, double *dir, int pmt, double *q, double *reflected, double *t) { - double pmt_dir[3], cos_theta, sin_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o; - - z = 1.0; + double pmt_dir[3], cos_theta, sin_theta, n, omega, R, f, f_reflec, cos_theta_pmt, sin_theta_pmt, charge, prob_abs, prob_sct, l_d2o, l_h2o, cos_theta_cerenkov; R = NORM(pos); n = (R <= AV_RADIUS) ? get_avg_index_d2o() : get_avg_index_h2o(); + cos_theta_cerenkov = 1/(beta*n); + *q = 0.0; *t = 0.0; *reflected = 0.0; - if (beta < 1/n) return; + if (cos_theta_cerenkov > 1) return; SUB(pmt_dir,pmts[pmt].pos,pos); @@ -569,17 +569,17 @@ static void get_expected_charge(double beta, double theta0, double *pos, double cos_theta = DOT(dir,pmt_dir); sin_theta = fast_sqrt(1-cos_theta*cos_theta); - cos_theta_pmt = DOT(pmt_dir,pmts[pmt].normal); + cos_theta_pmt = -DOT(pmt_dir,pmts[pmt].normal); sin_theta_pmt = fast_sqrt(1-cos_theta_pmt*cos_theta_pmt); theta0 = fmax(theta0,get_theta0_min(R,PMT_RADIUS,sin_theta_pmt)); - if (fabs(cos_theta-1.0/(n*beta))/(sin_theta*theta0) > 5) return; + if (fabs(cos_theta-cos_theta_cerenkov)/(sin_theta*theta0) > 5) return; 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); + 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); @@ -603,7 +603,7 @@ static void get_expected_charge(double beta, double theta0, double *pos, double * pretty low, this is expected to be a very small effect. */ prob_sct = 1.0 - get_fsct_d2o(l_d2o)*get_fsct_h2o(l_h2o)*get_fsct_acrylic(AV_THICKNESS); - charge = omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, sin_theta, theta0); + charge = omega*(1-cos_theta_cerenkov*cos_theta_cerenkov)*get_probability(beta, cos_theta, sin_theta, theta0); *reflected = (1.0-prob_abs)*(1.0-prob_sct)*f_reflec*charge + prob_sct*charge; diff --git a/src/scattering.c b/src/scattering.c index 2512cda..740374e 100644 --- a/src/scattering.c +++ b/src/scattering.c @@ -66,7 +66,7 @@ static double prob_scatter(double wavelength, void *params) delta = (1.0/index - beta_cos_theta)/beta_sin_theta_theta0; - return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI); + return FINE_STRUCTURE_CONSTANT*qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI)/beta_sin_theta_theta0; } static double prob_scatter2(double wavelength, void *params) @@ -194,10 +194,10 @@ void init_interpolation(void) double get_probability(double beta, double cos_theta, double sin_theta, double theta0) { /* Make sure theta0 is less than MAX_THETA0, otherwise it's possible that - * gsl_spline2d_eval() will quit. */ + * interp2d() will quit. */ if (theta0 > MAX_THETA0) theta0 = MAX_THETA0; - return interp2d(beta*cos_theta, beta*sin_theta*theta0, x, y, z, nx, ny)/(theta0*sin_theta); + return beta*interp2d(beta*cos_theta, beta*sin_theta*theta0, x, y, z, nx, ny); } double get_probability2(double beta) |