diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/misc.c | 27 | ||||
| -rw-r--r-- | src/misc.h | 2 | ||||
| -rw-r--r-- | src/path.c | 10 | ||||
| -rw-r--r-- | src/path.h | 2 | ||||
| -rw-r--r-- | src/scattering.c | 2 | ||||
| -rw-r--r-- | src/scattering.h | 5 | ||||
| -rw-r--r-- | src/test-path.c | 197 | ||||
| -rw-r--r-- | src/test.c | 62 |
8 files changed, 275 insertions, 32 deletions
@@ -397,3 +397,30 @@ double norm_cdf(double x, double mu, double sigma) * standard deviation `sigma`. */ return erfc(-(x-mu)/(sqrt(2)*sigma))/2.0; } + +double mean(const double *x, size_t n) +{ + /* Returns the mean of the array `x`. */ + size_t i; + double sum = 0.0; + + for (i = 0; i < n; i++) + sum += x[i]; + + return sum/n; +} + +double std(const double *x, size_t n) +{ + /* Returns the standard deviation of the array `x`. */ + size_t i; + double sum, mu; + + mu = mean(x,n); + + sum = 0.0; + for (i = 0; i < n; i++) + sum += pow(x[i]-mu,2); + + return sqrt(sum/n); +} @@ -16,5 +16,7 @@ int allclose(double *a, double *b, size_t n, double rel_tol, double abs_tol); double logsumexp(double *a, size_t n); double norm(double x, double mu, double sigma); double norm_cdf(double x, double mu, double sigma); +double mean(const double *x, size_t n); +double std(const double *x, size_t n); #endif @@ -8,6 +8,7 @@ #include <stdlib.h> /* for malloc(), calloc(), etc. */ #include "misc.h" #include <gsl/gsl_sf_psi.h> /* for gsl_sf_psi_n() */ +#include "scattering.h" #define N 1000 @@ -156,8 +157,9 @@ path *path_init(double *pos, double *dir, double T0, double range, double theta0 * Since there is no nice closed form solution for this, we estimate it by * evaluating the result at the middle of the track. */ if (p->n > 0) { - p->theta0 = 4*range*pow(theta0,2)*gsl_sf_psi_n(1,p->n-0.5)/pow(M_PI,2); + p->theta0 = theta0*sqrt(range)*sqrt(gsl_sf_psi_n(1,p->n+0.5))/M_PI; p->theta0 = fmax(MIN_THETA0, p->theta0); + p->theta0 = fmin(MAX_THETA0, p->theta0); } else { p->theta0 = theta0; } @@ -198,10 +200,12 @@ void path_eval(path *p, double s, double *pos, double *dir, double *T, double *t normalize(dir); - if (p->n > 0) + if (p->n > 0) { *theta0 = p->theta0; - else + } else { *theta0 = fmax(MIN_THETA0,p->theta0*sqrt(s)); + *theta0 = fmin(MAX_THETA0,*theta0); + } } void path_free(path *p) @@ -12,7 +12,7 @@ * distribution such that it effectively averages across the face of a PMT. * * FIXME: Should do some tests to figure out what is the best value. */ -#define MIN_THETA0 0.02 +#define MIN_THETA0 0.01 typedef double getKineticEnergyFunc(double x, void *params); diff --git a/src/scattering.c b/src/scattering.c index c95f31e..0247c62 100644 --- a/src/scattering.c +++ b/src/scattering.c @@ -18,7 +18,7 @@ static double xlo = -1.0; static double xhi = 1.0; static size_t nx = 1000; static double ylo = 0.0; -static double yhi = 1.0; +static double yhi = MAX_THETA0; static size_t ny = 1000; static double *x; diff --git a/src/scattering.h b/src/scattering.h index 308f8ef..65765b0 100644 --- a/src/scattering.h +++ b/src/scattering.h @@ -1,6 +1,11 @@ #ifndef SCATTERING_H #define SCATTERING_H +/* Maximum value for the scattering RMS `theta0`. + * + * This is to prevent the interpolation from going out of bounds. */ +#define MAX_THETA0 1.0 + void init_interpolation(void); double get_probability(double beta, double cos_theta, double theta0); double get_probability2(double beta); diff --git a/src/test-path.c b/src/test-path.c index bca82a4..523ebd4 100644 --- a/src/test-path.c +++ b/src/test-path.c @@ -7,10 +7,15 @@ #include <errno.h> /* for errno */ #include <string.h> /* for strerror() */ #include "random.h" +#include <gsl/gsl_histogram.h> +#include "misc.h" typedef struct sim { double *pos; + double *dir; path *p; + double *z1; + double *z2; } sim; static double getKineticEnergy(double x, void *params) @@ -21,7 +26,7 @@ static double getKineticEnergy(double x, void *params) sim *simulate_path(double *pos0, double *dir0, double range, double theta0, size_t N, size_t n) { size_t i; - double *s, *theta1, *theta2, *x, *y, *z, *z1, *z2, tmp1[3], tmp2[3], k[3], normal[3]; + double *s, *theta1, *theta2, *x, *y, *z, tmp[3], tmp2[3], k[3], normal[3]; double dx, dy, dz, theta, phi; sim *simvalue = malloc(sizeof(sim)); @@ -33,10 +38,9 @@ sim *simulate_path(double *pos0, double *dir0, double range, double theta0, size y = calloc(N,sizeof(double)); z = calloc(N,sizeof(double)); simvalue->pos = calloc(N*3,sizeof(double)); - z1 = calloc(n,sizeof(double)); - z2 = calloc(n,sizeof(double)); - - init_genrand(0); + simvalue->dir = calloc(N*3,sizeof(double)); + simvalue->z1 = calloc(n,sizeof(double)); + simvalue->z2 = calloc(n,sizeof(double)); simvalue->pos[0] = pos0[0]; simvalue->pos[1] = pos0[1]; @@ -59,9 +63,9 @@ sim *simulate_path(double *pos0, double *dir0, double range, double theta0, size y[i] = y[i-1] + dy; z[i] = z[i-1] + dz; - tmp1[0] = x[i]; - tmp1[1] = y[i]; - tmp1[2] = z[i]; + tmp[0] = x[i]; + tmp[1] = y[i]; + tmp[2] = z[i]; k[0] = 0.0; k[1] = 0.0; @@ -73,19 +77,31 @@ sim *simulate_path(double *pos0, double *dir0, double range, double theta0, size phi = acos(DOT(k,dir0)); - rotate(tmp2,tmp1,normal,phi); + rotate(tmp2,tmp,normal,phi); simvalue->pos[i*3] = pos0[0] + tmp2[0]; simvalue->pos[i*3+1] = pos0[1] + tmp2[1]; simvalue->pos[i*3+2] = pos0[2] + tmp2[2]; + + tmp[0] = dx; + tmp[1] = dy; + tmp[2] = dz; + + normalize(tmp); + + rotate(tmp2,tmp,normal,phi); + + simvalue->dir[i*3] = tmp2[0]; + simvalue->dir[i*3+1] = tmp2[1]; + simvalue->dir[i*3+2] = tmp2[2]; } for (i = 0; i < n; i++) { - z1[i] = path_get_coefficient(i+1,s,theta1,theta0,N); - z2[i] = path_get_coefficient(i+1,s,theta2,theta0,N); + simvalue->z1[i] = path_get_coefficient(i+1,s,theta1,theta0,N); + simvalue->z2[i] = path_get_coefficient(i+1,s,theta2,theta0,N); } - simvalue->p = path_init(pos0,dir0,1.0,range,theta0,getKineticEnergy,NULL,z1,z2,n,1.0); + simvalue->p = path_init(pos0,dir0,1.0,range,theta0,getKineticEnergy,NULL,simvalue->z1,simvalue->z2,n,1.0); free(s); free(theta1); @@ -93,32 +109,45 @@ sim *simulate_path(double *pos0, double *dir0, double range, double theta0, size free(x); free(y); free(z); - free(z1); - free(z2); return simvalue; } +void sim_free(sim *simvalue) +{ + free(simvalue->pos); + free(simvalue->dir); + free(simvalue->z1); + free(simvalue->z2); + path_free(simvalue->p); + free(simvalue); +} + void usage(void) { fprintf(stderr,"Usage: ./test-path [options]\n"); fprintf(stderr," -N number of points along track\n"); - fprintf(stderr," -n number of terms in KL expansion\n"); + fprintf(stderr," -z number of terms in KL expansion\n"); + fprintf(stderr," -n number of simulations\n"); fprintf(stderr," -t standard deviation of angular distribution\n"); fprintf(stderr," --range range of track\n"); + fprintf(stderr," -b number of bins\n"); fprintf(stderr," -h display this help message\n"); exit(1); } int main(int argc, char **argv) { - size_t i, n, N; - double pos0[3], dir0[3], theta0, range, pos2[3], dir[3], T, t; + size_t i, j, n, N, z, bins; + double pos0[3], dir0[3], theta0, range, pos2[3], dir[3], T, t, tmp; + double *theta1, *theta2; - n = 2; + z = 2; + n = 1000; N = 1000; - theta0 = 0.1; + theta0 = 0.01; range = 100.0; + bins = 100; for (i = 1; i < argc; i++) { if (!strncmp(argv[i], "--", 2)) { @@ -134,9 +163,15 @@ int main(int argc, char **argv) case 'n': n = atoi(argv[++i]); break; + case 'z': + z = atoi(argv[++i]); + break; case 't': theta0 = strtod(argv[++i],NULL); break; + case 'b': + bins = atoi(argv[++i]); + break; case 'h': usage(); default: @@ -154,28 +189,136 @@ int main(int argc, char **argv) dir0[1] = 0.0; dir0[2] = 0.0; - sim *simvalue = simulate_path(pos0,dir0,range,theta0,N,n); + theta1 = calloc(N*n,sizeof(double)); + theta2 = calloc(N*n,sizeof(double)); + + gsl_histogram **h1 = malloc(sizeof(gsl_histogram *)*z); + gsl_histogram **h2 = malloc(sizeof(gsl_histogram *)*z); - FILE *pipe = popen("graph -T X --bitmap-size 2000x2000 -X X -Y Y", "w"); + for (i = 0; i < z; i++) { + h1[i] = gsl_histogram_alloc(bins); + gsl_histogram_set_ranges_uniform(h1[i],-10.0,10.0); + h2[i] = gsl_histogram_alloc(bins); + gsl_histogram_set_ranges_uniform(h2[i],-10.0,10.0); + } + + init_genrand(0); + + for (i = 0; i < n; i++) { + sim *simvalue = simulate_path(pos0,dir0,range,theta0,N,z); + for (j = 0; j < z; j++) { + gsl_histogram_increment(h1[j], simvalue->z1[j]); + gsl_histogram_increment(h2[j], simvalue->z2[j]); + } + for (j = 0; j < N; j++) { + path_eval(simvalue->p,range*j/(N-1),pos2,dir,&T,&t,&tmp); + theta1[j*n+i] = simvalue->dir[j*3+1] - dir[1]; + theta2[j*n+i] = simvalue->dir[j*3+2] - dir[2]; + } + sim_free(simvalue); + } + + FILE *pipe = popen("gnuplot -p", "w"); if (!pipe) { - fprintf(stderr, "error running graph command: %s\n", strerror(errno)); + fprintf(stderr, "error running gnuplot command: %s\n", strerror(errno)); exit(1); } + fprintf(pipe, "set title 'Path Coefficients'\n"); + fprintf(pipe, "plot "); + for (i = 0; i < z; i++) { + fprintf(pipe, "'-' title 'z1[%zu]' with steps, \\\n",i); + fprintf(pipe, "'-' title 'z2[%zu]' with steps",i); + if (i != z-1) fprintf(pipe, ", \\"); + fprintf(pipe, "\n"); + } + + for (i = 0; i < z; i++) { + for (j = 0; j < h1[i]->n; j++) { + fprintf(pipe, "%g %g\n", h1[i]->range[j], h1[i]->bin[j]); + } + fprintf(pipe, "e\n"); + for (j = 0; j < h2[i]->n; j++) { + fprintf(pipe, "%g %g\n", h2[i]->range[j], h2[i]->bin[j]); + } + if (i != z-1) fprintf(pipe, "e\n"); + } + + for (i = 0; i < z; i++) { + gsl_histogram_free(h1[i]); + gsl_histogram_free(h2[i]); + } + + free(h1); + free(h2); + + if (pclose(pipe)) { + fprintf(stderr, "error closing gnuplot command: %s\n", strerror(errno)); + exit(1); + } + + pipe = popen("gnuplot -p", "w"); + + if (!pipe) { + fprintf(stderr, "error running gnuplot command: %s\n", strerror(errno)); + exit(1); + } + + fprintf(pipe, "set title 'Standard Deviation of the Path Direction along the Track'\n"); + fprintf(pipe, "set xlabel 'Distance (cm)'\n"); + fprintf(pipe, "set ylabel '{/Symbol s}({/Symbol q})'\n"); + fprintf(pipe, "plot '-' title 'x' with lines, \\\n"); + fprintf(pipe, " '-' title 'y' with lines, \\\n"); + fprintf(pipe, " '-' title 'predicted' with lines\n"); + fprintf(pipe, "\n"); + for (i = 0; i < N; i++) { + fprintf(pipe, "%g %g\n", i*range/(N-1), std(theta1+n*i,n)); + } + fprintf(pipe, "e\n"); + for (i = 0; i < N; i++) { + fprintf(pipe, "%g %g\n", i*range/(N-1), std(theta2+n*i,n)); + } + fprintf(pipe, "e\n"); + sim *simvalue = simulate_path(pos0,dir0,range,theta0,N,z); for (i = 0; i < N; i++) { - fprintf(pipe,"%.10g %.10g\n", simvalue->pos[i*3], simvalue->pos[i*3+1]); + path_eval(simvalue->p,range*i/(N-1),pos2,dir,&T,&t,&tmp); + fprintf(pipe, "%g %g\n", i*range/(N-1), tmp); } - fprintf(pipe,"\n\n"); + + if (pclose(pipe)) { + fprintf(stderr, "error closing gnuplot command: %s\n", strerror(errno)); + exit(1); + } + + pipe = popen("gnuplot -p", "w"); + + if (!pipe) { + fprintf(stderr, "error running gnuplot command: %s\n", strerror(errno)); + exit(1); + } + + fprintf(pipe, "set title 'Simulated Path'\n"); + fprintf(pipe, "splot '-' u 1:2:3 title 'simulation' with lines, \\\n"); + fprintf(pipe, " '-' u 1:2:3 title 'path' with lines\n"); + for (i = 0; i < N; i++) { + fprintf(pipe,"%.10g %.10g %.10g\n", simvalue->pos[i*3], simvalue->pos[i*3+1], simvalue->pos[i*3+2]); + } + fprintf(pipe,"e\n"); for (i = 0; i < N; i++) { path_eval(simvalue->p,i*range/(N-1),pos2,dir,&T,&t,&theta0); - fprintf(pipe,"%.10g %.10g\n", pos2[0], pos2[1]); + fprintf(pipe,"%.10g %.10g %.10g\n", pos2[0], pos2[1], pos2[2]); } - fprintf(pipe,"\n\n"); + fprintf(pipe,"e\n"); + + /* Pause so that you can rotate the 3D graph. */ + fprintf(pipe,"pause mouse keypress\n"); + + sim_free(simvalue); if (pclose(pipe)) { - fprintf(stderr, "error closing graph command: %s\n", strerror(errno)); + fprintf(stderr, "error closing gnuplot command: %s\n", strerror(errno)); exit(1); } @@ -13,6 +13,7 @@ #include "pdg.h" #include <gsl/gsl_spline.h> #include "vector.h" +#include <gsl/gsl_statistics.h> typedef int testFunction(char *err); @@ -657,6 +658,65 @@ err: return 1; } +int test_mean(char *err) +{ + /* Tests the mean() function. */ + size_t i, j; + double x[100]; + double mu, expected; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + for (j = 0; j < sizeof(x)/sizeof(x[0]); j++) { + x[j] = genrand_real2(); + } + + mu = mean(x,sizeof(x)/sizeof(x[0])); + expected = gsl_stats_mean(x,1,sizeof(x)/sizeof(x[0])); + + if (!isclose(mu, expected, 0, 1e-9)) { + sprintf(err, "mean() returned %.5g, but expected %.5g", mu, expected); + goto err; + } + } + + return 0; + +err: + return 1; +} + +int test_std(char *err) +{ + /* Tests the mean() function. */ + size_t i, j; + double x[100]; + double sigma, expected, mu; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + for (j = 0; j < sizeof(x)/sizeof(x[0]); j++) { + x[j] = genrand_real2(); + } + + sigma = std(x,sizeof(x)/sizeof(x[0])); + mu = gsl_stats_mean(x,1,sizeof(x)/sizeof(x[0])); + expected = gsl_stats_sd_with_fixed_mean(x,1,sizeof(x)/sizeof(x[0]),mu); + + if (!isclose(sigma, expected, 0, 1e-9)) { + sprintf(err, "std() returned %.5g, but expected %.5g", sigma, expected); + goto err; + } + } + + return 0; + +err: + return 1; +} + struct tests { testFunction *test; char *name; @@ -675,6 +735,8 @@ struct tests { {test_lnfact, "test_lnfact"}, {test_ln, "test_ln"}, {test_get_path_length, "test_get_path_length"}, + {test_mean, "test_mean"}, + {test_std, "test_std"}, }; int main(int argc, char **argv) |