diff options
| -rw-r--r-- | src/muon.c | 2 | ||||
| -rw-r--r-- | src/scattering.c | 9 | ||||
| -rw-r--r-- | src/test-likelihood.c | 37 |
3 files changed, 26 insertions, 22 deletions
@@ -262,5 +262,5 @@ double get_expected_charge(double x, double T, double T0, double *pos, double *d theta0 = get_scattering_rms(x,p0,beta0,z,rho); /* FIXME: add angular response and scattering/absorption. */ - return 2*omega*2*M_PI*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0)/(sqrt(2*M_PI)*theta0); + return omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); } diff --git a/src/scattering.c b/src/scattering.c index a1a63df..66e8398 100644 --- a/src/scattering.c +++ b/src/scattering.c @@ -40,12 +40,9 @@ static double prob_scatter(double wavelength, void *params) index = get_index_snoman_d2o(wavelength); - delta = (1.0/index - beta_cos_theta)/(2*beta_sin_theta_theta0); + delta = (1.0/index - beta_cos_theta)/beta_sin_theta_theta0; - /* FIXME: ignore GSL error for underflow here. */ - if (delta*delta > 500) return 0.0; - else if (delta*delta == 0.0) return INFINITY; - return qe*exp(-delta*delta)*gsl_sf_bessel_K0(delta*delta)/pow(wavelength,2)*1e7/(4*M_PI*M_PI); + return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*M_PI); } void init_interpolation(void) @@ -106,7 +103,7 @@ double get_probability(double beta, double cos_theta, double theta0) * we are going to square it everywhere. */ sin_theta = fabs(sin(acos(cos_theta))); - return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/sin_theta; + return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/(theta0*sin_theta); } void free_interpolation(void) diff --git a/src/test-likelihood.c b/src/test-likelihood.c index f0266b9..45bbfe2 100644 --- a/src/test-likelihood.c +++ b/src/test-likelihood.c @@ -21,7 +21,7 @@ void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double t * distribution is simulated as a gaussian distribution with standard * deviation `theta0`. */ int i; - double theta, phi, wavelength, u, qe, index, cerenkov_angle, dir[3], n[3], dest[3], E, p, beta, cos_theta; + double theta, phi, wavelength, u, qe, index, cerenkov_angle, dir[3], n[3], dest[3], E, p, beta, cos_theta, thetax, thetay; i = 0; while (i < N) { @@ -35,7 +35,7 @@ void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double t /* Check to see if the photon was detected. */ if (genrand_real2() > qe) continue; - index = get_index(HEAVY_WATER_DENSITY, wavelength, 10.0); + index = get_index_snoman_d2o(wavelength); /* Calculate total energy */ E = T + MUON_MASS; @@ -45,21 +45,28 @@ void simulate_cos_theta_distribution(int N, gsl_histogram *h, double T, double t cerenkov_angle = acos(1/(index*beta)); /* Assuming the muon track is dominated by small angle scattering, the - * angular distribution will be a Gaussian centered around 0 with a - * standard deviation of `theta0`. Here, we draw a random angle from - * this distribution. */ - theta = randn()*theta0; - - n[0] = sin(theta); - n[1] = 0; + * angular distribution looks like the product of two uncorrelated + * Gaussian distributions with a standard deviation of `theta0` in the + * plane perpendicular to the track direction. Here, we draw two random + * angles and then compute the polar and azimuthal angle for the track + * direction. */ + thetax = randn()*theta0; + thetay = randn()*theta0; + + theta = sqrt(thetax*thetax + thetay*thetay); + phi = atan2(thetay,thetax); + + n[0] = sin(theta)*cos(phi); + n[1] = sin(theta)*sin(phi); n[2] = cos(theta); - /* To compute the direction of the photon, we start with a vector in - * the x-z plane which is offset from the track direction by the - * Cerenkov angle and then rotate it around the track direction by a - * random angle `phi`. */ - dir[0] = sin(cerenkov_angle + theta); - dir[1] = 0; + /* To compute the direction of the photon, we start with a vector which + * has the same azimuthal angle as the track direction but is offset + * from the track direction in the polar angle by the Cerenkov angle + * and then rotate it around the track direction by a random angle + * `phi`. */ + dir[0] = sin(cerenkov_angle + theta)*cos(phi); + dir[1] = sin(cerenkov_angle + theta)*sin(phi); dir[2] = cos(cerenkov_angle + theta); phi = genrand_real2()*2*M_PI; |