diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/likelihood.c | 9 | ||||
| -rw-r--r-- | src/muon.c | 9 | ||||
| -rw-r--r-- | src/optics.c | 6 | ||||
| -rw-r--r-- | src/optics.h | 1 |
4 files changed, 19 insertions, 6 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 6881e8f..fca14cc 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -153,7 +153,7 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy * * `smax` is currently calculated as the point where the particle velocity * drops to 0.8 times the speed of light. */ - double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length, wavelength0; + double pmt_dir[3], tmp[3], R, cos_theta, theta, x, z, s, a, b, beta, E, p, T, omega, theta_cerenkov, n, sin_theta, E0, p0, beta0, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, wavelength0; /* Calculate beta at the start of the track. */ E0 = T0 + MUON_MASS; @@ -253,9 +253,12 @@ double get_total_charge_approx(double T0, double *pos, double *dir, muon_energy f = get_weighted_pmt_response(acos(-cos_theta_pmt)); - absorption_length = get_absorption_length_snoman_d2o(wavelength0); + absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0); + absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0); - return f*exp(-a/absorption_length)*n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI); + get_path_length(tmp,pmts[i].pos,AV_RADIUS,&l_d2o,&l_h2o); + + return f*exp(-l_d2o/absorption_length_d2o)*exp(-l_h2o/absorption_length_h2o)*n*x*beta0*prob*(1/sin_theta)*omega*(erf((a+b*(smax-s)+n*(smax-z)*beta0)/frac) + erf((-a+b*s+n*z*beta0)/frac))/(b+n*beta0)/(4*M_PI); } typedef struct betaRootParams { @@ -258,7 +258,7 @@ double get_dEdx(double T, double rho) 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, E, p, beta, z, R, f, cos_theta_pmt, absorption_length, distance; + double pmt_dir[3], cos_theta, n, wavelength0, omega, E, p, beta, z, R, f, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, l_h2o, l_d2o, distance; z = 1.0; @@ -298,8 +298,11 @@ double get_expected_charge(double x, double T, double theta0, double *pos, doubl f = get_weighted_pmt_response(acos(-cos_theta_pmt)); - absorption_length = get_absorption_length_snoman_d2o(wavelength0); + absorption_length_d2o = get_absorption_length_snoman_d2o(wavelength0); + absorption_length_h2o = get_absorption_length_snoman_h2o(wavelength0); + + get_path_length(pos,pmt_pos,AV_RADIUS,&l_d2o,&l_h2o); /* FIXME: add angular response and scattering/absorption. */ - return f*exp(-distance/absorption_length)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); + return f*exp(-l_d2o/absorption_length_d2o)*exp(-l_h2o/absorption_length_h2o)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); } diff --git a/src/optics.c b/src/optics.c index 99c2a51..9569865 100644 --- a/src/optics.c +++ b/src/optics.c @@ -26,6 +26,12 @@ static double RIND_D2O_C3 = 0.32; * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */ static double HC = 1.239841973976e3; +double get_absorption_length_snoman_h2o(double wavelength) +{ + /* Returns the absorption length in H2O in cm. */ + return 1e4; +} + double get_absorption_length_snoman_d2o(double wavelength) { /* Returns the absorption length in D2O in cm. */ diff --git a/src/optics.h b/src/optics.h index 5a6dec9..13b3e0b 100644 --- a/src/optics.h +++ b/src/optics.h @@ -1,6 +1,7 @@ #ifndef OPTICS_H #define OPTICS_H +double get_absorption_length_snoman_h2o(double wavelength); double get_absorption_length_snoman_d2o(double wavelength); double get_index(double p, double wavelength, double T); double get_index_snoman_h2o(double wavelength); |