diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/fit.c | 2 | ||||
| -rw-r--r-- | src/likelihood.c | 4 | ||||
| -rw-r--r-- | src/likelihood.h | 2 | ||||
| -rw-r--r-- | src/muon.c | 11 | ||||
| -rw-r--r-- | src/muon.h | 2 |
5 files changed, 13 insertions, 8 deletions
@@ -134,7 +134,7 @@ int main(int argc, char **argv) /* New event, so we need to fit the old event. */ iter = 0; - gsl_vector_set(x,0,1000.0); + gsl_vector_set(x,0,500.0); gsl_vector_set(x,1,100.0); gsl_vector_set(x,2,100.0); gsl_vector_set(x,3,100.0); diff --git a/src/likelihood.c b/src/likelihood.c index 98a4ad7..4e24a79 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -103,7 +103,7 @@ static double gsl_muon_time(double x, void *params) t = x/SPEED_OF_LIGHT + distance*n/SPEED_OF_LIGHT; - return t*get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS); + return t*get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS); } static double gsl_muon_charge(double x, void *params) @@ -124,7 +124,7 @@ static double gsl_muon_charge(double x, void *params) pos[1] = pos0[1] + dir[1]*x; pos[2] = pos0[2] + dir[2]*x; - return get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS); + return get_expected_charge(x, get_T(T0, x, HEAVY_WATER_DENSITY), T0, pos, dir, pmts[i].pos, pmts[i].normal, PMT_RADIUS); } double nll_muon(event *ev, double T, double *pos, double *dir, double t0) diff --git a/src/likelihood.h b/src/likelihood.h index ebfeb97..43b0823 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -4,7 +4,7 @@ #include "event.h" #define MAX_PE 100 -#define CHARGE_FRACTION 150000.0 +#define CHARGE_FRACTION 60000.0 #define DARK_RATE 500.0 /* Single PE transit time spread (ns). */ #define PMT_TTS 1.5 @@ -216,9 +216,9 @@ double get_dEdx(double T, double rho) return gsl_spline_eval(spline_dEdx, T, acc_dEdx)/rho; } -double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r) +double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r) { - double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R; + double pmt_dir[3], cos_theta, n, wavelength0, omega, theta0, E, p, beta, z, rho, R, E0, p0, beta0; z = 1.0; @@ -253,8 +253,13 @@ double get_expected_charge(double x, double T, double *pos, double *dir, double if (beta < 1/n) return 0; + /* Calculate total energy */ + E0 = T0 + MUON_MASS; + p0 = sqrt(E0*E0 - MUON_MASS*MUON_MASS); + beta0 = p0/E0; + /* FIXME: is this formula valid for muons? */ - theta0 = get_scattering_rms(x,p,beta,z,rho); + 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); @@ -6,6 +6,6 @@ double get_range(double T, double rho); double get_T(double T0, double x, double rho); double get_dEdx(double T, double rho); -double get_expected_charge(double x, double T, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r); +double get_expected_charge(double x, double T, double T0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r); #endif |