diff options
| -rw-r--r-- | src/likelihood.c | 2 | ||||
| -rw-r--r-- | src/solid_angle.c | 19 | ||||
| -rw-r--r-- | src/solid_angle.h | 1 | ||||
| -rw-r--r-- | src/test.c | 26 |
4 files changed, 47 insertions, 1 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 2340fb1..ed07826 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -424,7 +424,7 @@ static double get_total_charge_approx(double T0, double *pos, double *dir, parti sin_theta = sqrt(1-pow(cos_theta,2)); /* Get the solid angle of the PMT at the position `s` along the track. */ - omega = get_solid_angle_approx(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS); + omega = get_solid_angle_fast(tmp,pmts[i].pos,pmts[i].normal,PMT_RADIUS); theta0 = fmax(theta0*sqrt(s),MIN_THETA0); diff --git a/src/solid_angle.c b/src/solid_angle.c index a929938..dbac045 100644 --- a/src/solid_angle.c +++ b/src/solid_angle.c @@ -3,6 +3,7 @@ #include <math.h> #include <gsl/gsl_interp2d.h> #include <gsl/gsl_spline2d.h> +#include "vector.h" static double hd[11] = {0.1,0.125,0.150,0.175,0.2,0.3,0.4,0.5,0.6,0.8,1.0}; static double Rd[13] = {0.1,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2.0}; @@ -23,6 +24,24 @@ static double lookupTable[13][11] = { {1.0041,1.0051,1.0061,1.0070,1.0078,1.0106,1.0120,1.0122,1.0116,1.0097,1.0084} }; +double get_solid_angle_fast(double *pos, double *pmt, double *n, double r) +{ + /* Returns a *very* fast approximation of the solid angle subtended by a + * circular disk of radius r at a position `pmt` with a normal vector `n` from + * a position `pos`. */ + double dir[3]; + double R, cos_theta_pmt; + + SUB(dir,pos,pmt); + + /* Distance to the PMT. */ + R = NORM(dir); + + cos_theta_pmt = fabs(DOT(dir,n)/R); + + return pow(r,2)*cos_theta_pmt/pow(R,2); +} + static double A(double u, double a, double h) { return atan(((a*a-h*h)*u - 2*a*a*h*h)/(2*a*h*sqrt((u-h*h)*(u+a*a)))); diff --git a/src/solid_angle.h b/src/solid_angle.h index aacb692..654f605 100644 --- a/src/solid_angle.h +++ b/src/solid_angle.h @@ -3,6 +3,7 @@ #include <math.h> +double get_solid_angle_fast(double *pos, double *pmt, double *n, double r); double get_solid_angle_approx(double *pos, double *pmt, double *n, double r); double get_solid_angle(double *pos, double *pmt, double *n, double r); @@ -970,6 +970,31 @@ err: return 1; } +int test_solid_angle_fast(char *err) +{ + /* Tests the get_solid_angle_fast() function. */ + int i; + double pmt[3] = {0,0,0}; + double pos[3] = {0,0,1}; + double n[3] = {0,0,-1}; + double r = 1.0; + double solid_angle; + + for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) { + pos[0] = solid_angle_results[i].r0*r; + pos[2] = solid_angle_results[i].L*r; + + solid_angle = get_solid_angle_approx(pos,pmt,n,r); + + if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-2, 0)) { + sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega); + return 1; + } + } + + return 0; +} + struct tests { testFunction *test; char *name; @@ -999,6 +1024,7 @@ struct tests { {test_log_norm, "test_log_norm"}, {test_trapz, "test_trapz"}, {test_interp2d, "test_interp2d"}, + {test_solid_angle_fast, "test_solid_angle_fast"}, }; int main(int argc, char **argv) |