diff options
| -rw-r--r-- | src/.gitignore | 1 | ||||
| -rw-r--r-- | src/Makefile | 4 | ||||
| -rw-r--r-- | src/fit.c | 11 | ||||
| -rw-r--r-- | src/likelihood.c | 39 | ||||
| -rw-r--r-- | src/likelihood.h | 4 | ||||
| -rw-r--r-- | src/optics.c | 112 | ||||
| -rw-r--r-- | src/optics.h | 10 | ||||
| -rw-r--r-- | src/rsp_rayleigh.dat | 88 | ||||
| -rw-r--r-- | src/test-time-pdf.c | 116 | ||||
| -rw-r--r-- | src/test.c | 135 |
10 files changed, 503 insertions, 17 deletions
diff --git a/src/.gitignore b/src/.gitignore index e87f8da..f613ee3 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -8,3 +8,4 @@ test-charge test-path release.h calculate-csda-range +test-time-pdf diff --git a/src/Makefile b/src/Makefile index 80eef49..c44c172 100644 --- a/src/Makefile +++ b/src/Makefile @@ -3,7 +3,7 @@ release_hdr := $(shell sh -c './mkreleasehdr.sh') CFLAGS=-O2 -Wall -g -DSWAP_BYTES LDLIBS=-lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++ -all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range +all: test test-vector test-likelihood fit test-charge test-path calculate-csda-range test-time-pdf Makefile.dep: -$(CC) -MM *.c > Makefile.dep @@ -24,6 +24,8 @@ test-path: test-path.o mt19937ar.o vector.o path.o random.o misc.o test-charge: test-charge.o sno_charge.o misc.o vector.o +test-time-pdf: test-time-pdf.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o + test-zebra: test-zebra.o zebra.o pack2b.o fit: fit.o zebra.o likelihood.o pmt.o vector.o misc.o muon.o optics.o quantum_efficiency.o solid_angle.o pdg.o scattering.o zdab_utils.o pack2b.o sno_charge.o db.o dqxx.o dict.o siphash.o path.o pmt_response.o release.o electron.o proton.o @@ -5482,7 +5482,6 @@ int main(int argc, char **argv) init_interpolation(); init_charge(); - optics_init(); dict *db = db_init(); if (load_file(db, "DQXX_0000010000.dat")) { @@ -5500,11 +5499,21 @@ int main(int argc, char **argv) exit(1); } + if (load_file(db, "rsp_rayleigh.dat")) { + fprintf(stderr, "failed to load rsp_rayleigh.dat: %s\n", db_err); + exit(1); + } + if (pmt_response_init(db)) { fprintf(stderr, "failed to initialize PMTR bank: %s\n", pmtr_err); exit(1); } + if (optics_init(db)) { + fprintf(stderr, "failed to initialize optics: %s\n", optics_err); + exit(1); + } + int first_ev = 1; int first_mctk = 1; int first_mcvx = 1; diff --git a/src/likelihood.c b/src/likelihood.c index fcc8a5f..eb2eb94 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -171,7 +171,7 @@ void particle_free(particle *p) static double get_expected_charge_shower(double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o, double alpha, double beta, double rad) { - double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; + double pmt_dir[3], cos_theta, omega, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter; SUB(pmt_dir,pmt_pos,pos); @@ -195,19 +195,23 @@ static double get_expected_charge_shower(double *pos, double *dir, double *pmt_p absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic(); + scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); + scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*get_weighted_quantum_efficiency()*rad*PHOTONS_PER_MEV*electron_get_angular_pdf(cos_theta,alpha,beta,1/n_d2o)*omega/(2*M_PI); - *reflected = f_reflec*charge; + scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); - return f*charge; + *reflected = (f_reflec + 1.0 - scatter)*charge; + + return f*charge*scatter; } static double get_expected_charge(double x, double beta, double theta0, double *pos, double *dir, double *pmt_pos, double *pmt_normal, double r, double *reflected, double n_d2o, double n_h2o, double l_d2o, double l_h2o) { - double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge; + double pmt_dir[3], cos_theta, n, omega, z, R, f, f_reflec, cos_theta_pmt, absorption_length_h2o, absorption_length_d2o, absorption_length_acrylic, l_acrylic, theta_pmt, charge, scattering_length_h2o, scattering_length_d2o, scatter; z = 1.0; @@ -242,23 +246,34 @@ static double get_expected_charge(double x, double beta, double theta0, double * absorption_length_d2o = get_weighted_absorption_length_snoman_d2o(); absorption_length_h2o = get_weighted_absorption_length_snoman_h2o(); absorption_length_acrylic = get_weighted_absorption_length_snoman_acrylic(); + scattering_length_d2o = get_weighted_rayleigh_scattering_length_snoman_d2o(); + scattering_length_h2o = get_weighted_rayleigh_scattering_length_snoman_h2o(); l_acrylic = AV_RADIUS_OUTER - AV_RADIUS_INNER; charge = exp(-l_d2o/absorption_length_d2o-l_h2o/absorption_length_h2o-l_acrylic/absorption_length_acrylic)*omega*FINE_STRUCTURE_CONSTANT*z*z*(1-(1/(beta*beta*n*n)))*get_probability(beta, cos_theta, theta0); - *reflected = f_reflec*charge; + scatter = exp(-l_d2o/scattering_length_d2o-l_h2o/scattering_length_h2o); - return f*charge; + *reflected = (f_reflec + 1.0 - scatter)*charge; + + return f*charge*scatter; } -double F(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) +double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the CDF for the time distribution of photons at time `t`. */ size_t i; double p, mu_total; - p = mu_noise*t/GTVALID + mu_indirect*(pow(sigma,2)*norm(tmean,t,sigma) + (t-tmean)*norm_cdf(t,tmean,sigma))/(GTVALID-tmean); + p = mu_noise*t/GTVALID; + + if (t < tmean) + ; + else if (t < tmean + PSUP_REFLECTION_TIME) + p += mu_indirect*(t-tmean)/PSUP_REFLECTION_TIME; + else + p += mu_indirect; mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { @@ -275,7 +290,7 @@ double F(double t, double mu_noise, double mu_indirect, double *mu_direct, doubl return p/mu_total; } -double f(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) +double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma) { /* Returns the probability that a photon is detected at time `t`. * @@ -294,7 +309,7 @@ double f(double t, double mu_noise, double mu_indirect, double *mu_direct, doubl size_t i; double p, mu_total; - p = mu_noise/GTVALID + mu_indirect*norm_cdf(t,tmean,sigma)/(GTVALID-tmean); + p = mu_noise/GTVALID + mu_indirect*(erf((t-tmean)/(sqrt(2)*sigma))-erf((t-tmean-PSUP_REFLECTION_TIME)/(sqrt(2)*sigma)))/(2*PSUP_REFLECTION_TIME); mu_total = mu_noise + mu_indirect; for (i = 0; i < n; i++) { @@ -324,9 +339,9 @@ double log_pt(double t, size_t n, double mu_noise, double mu_indirect, double *m * convolve first which is what we do here. In addition, the problem is not * analytically tractable if you do things the other way around. */ if (n == 1) - return ln(n) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); + return ln(n) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); else - return ln(n) + (n-1)*log1p(-F(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(f(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); + return ln(n) + (n-1)*log1p(-time_cdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)) + log(time_pdf(t,mu_noise,mu_indirect,mu_direct,mu_shower,n2,ts,ts_shower,tmean,sigma,ts_sigma)); } static void integrate_path_shower(particle *p, double *x, double *pdf, double T0, double *pos0, double *dir0, int pmt, size_t n, double *mu_direct, double *mu_indirect, double *time, double *sigma, double rad) diff --git a/src/likelihood.h b/src/likelihood.h index 39b9028..e888d0a 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -4,6 +4,8 @@ #include "event.h" #include <stddef.h> /* for size_t */ +#define PSUP_REFLECTION_TIME 120.0 + /* Minimum number of points in the track integral. */ #define MIN_NPOINTS 100 @@ -64,6 +66,8 @@ typedef struct particle { particle *particle_init(int id, double T0, size_t n); double particle_get_energy(double x, particle *p); void particle_free(particle *p); +double time_pdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma); +double time_cdf(double t, double mu_noise, double mu_indirect, double *mu_direct, double *mu_shower, size_t n, double *ts, double *ts_shower, double tmean, double sigma, double *ts_sigma); double nll_muon(event *ev, int id, double T0, double *pos, double *dir, double t0, double *z1, double *z2, size_t n, double dx, double dx_shower, int fast); #endif diff --git a/src/optics.c b/src/optics.c index 2d28ba4..4a0722a 100644 --- a/src/optics.c +++ b/src/optics.c @@ -7,6 +7,10 @@ #include <gsl/gsl_errno.h> /* for gsl_strerror() */ #include <gsl/gsl_spline.h> #include "misc.h" +#include "db.h" +#include "dict.h" + +char optics_err[256]; /* Absorption coefficients for H2O and D2O as a function of wavelength from * SNOMAN. @@ -15,6 +19,10 @@ static double absorption_coefficient_h2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; static double absorption_coefficient_h2o[6] = {8.410e-05, 2.880e-06, 1.431e-04, 1.034e-04, 5.239e-04, 2.482e-03}; +static double rayleigh_coefficient_wavelengths[51]; +static double rayleigh_coefficient_h2o[51]; +static double rayleigh_coefficient_d2o[51]; + static double absorption_coefficient_d2o_wavelengths[6] = {337.0, 365.0, 386.0, 420.0, 500.0, 620.0}; static double absorption_coefficient_d2o[6] = {6.057e-05, 2.680e-05, 3.333e-05, 3.006e-05, 2.773e-05, 2.736e-05}; @@ -32,7 +40,12 @@ static double absorption_length_h2o_weighted = 0.0; /* Acrylic absorption length weighted by the Cerenkov spectrum and the quantum * efficiency. */ static double absorption_length_acrylic_weighted = 0.0; - +/* D2O Rayleigh scattering length weighted by the Cerenkov spectrum and the + * quantum efficiency. */ +static double rayleigh_scattering_length_d2o_weighted = 0.0; +/* H2O Rayleigh scattering length weighted by the Cerenkov spectrum and the + * quantum efficiency. */ +static double rayleigh_scattering_length_h2o_weighted = 0.0; /* From Table 4 in the paper. */ static double A0 = 0.243905091; @@ -59,6 +72,28 @@ static double RIND_D2O_C3 = 0.32; * From http://pdg.lbl.gov/2018/reviews/rpp2018-rev-phys-constants.pdf. */ static double HC = 1.239841973976e3; +static double gsl_rayleigh_scattering_length_snoman_d2o(double wavelength, void *params) +{ + /* Returns the Rayleigh scattering length in D2O weighted by the quantum efficiency + * and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_rayleigh_scattering_length_snoman_d2o(wavelength)/pow(wavelength,2); +} + +static double gsl_rayleigh_scattering_length_snoman_h2o(double wavelength, void *params) +{ + /* Returns the Rayleigh scattering length in H2O weighted by the quantum efficiency + * and the Cerenkov spectrum. */ + double qe; + + qe = get_quantum_efficiency(wavelength); + + return qe*get_rayleigh_scattering_length_snoman_h2o(wavelength)/pow(wavelength,2); +} + static double gsl_absorption_length_snoman_d2o(double wavelength, void *params) { /* Returns the absorption length in D2O weighted by the quantum efficiency @@ -102,14 +137,16 @@ static double gsl_cerenkov(double wavelength, void *params) return qe/pow(wavelength,2); } -int optics_init(void) +int optics_init(dict *db) { + int i; double norm; double result, error; size_t nevals; int status; gsl_integration_cquad_workspace *w; gsl_function F; + dbval *rspr; w = gsl_integration_cquad_workspace_alloc(100); @@ -156,6 +193,41 @@ int optics_init(void) exit(1); } + rspr = get_bank(db, "RSPR", 1); + + if (!rspr) { + sprintf(optics_err, "failed to load RSPR bank\n"); + return -1; + } + + for (i = 0; i < 51; i++) { + rayleigh_coefficient_wavelengths[i] = rspr[30+i*3].u32; + rayleigh_coefficient_d2o[i] = rspr[30+i*3+1].f; + rayleigh_coefficient_h2o[i] = rspr[30+i*3+2].f; + } + + F.function = &gsl_rayleigh_scattering_length_snoman_d2o; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + rayleigh_scattering_length_d2o_weighted = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + + F.function = &gsl_rayleigh_scattering_length_snoman_h2o; + + status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals); + + rayleigh_scattering_length_h2o_weighted = result/norm; + + if (status) { + fprintf(stderr, "error integrating cerenkov distribution: %s\n", gsl_strerror(status)); + exit(1); + } + gsl_integration_cquad_workspace_free(w); initialized = 1; @@ -175,6 +247,42 @@ double get_weighted_absorption_length_snoman_h2o(void) return absorption_length_h2o_weighted; } +double get_rayleigh_scattering_length_snoman_d2o(double wavelength) +{ + /* Returns the Rayleigh scattering length in D2O in cm. */ + return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_d2o,LEN(rayleigh_coefficient_wavelengths)); +} + +double get_rayleigh_scattering_length_snoman_h2o(double wavelength) +{ + /* Returns the Rayleigh scattering length in H2O in cm. */ + return 1.0/interp1d(wavelength,rayleigh_coefficient_wavelengths, rayleigh_coefficient_h2o,LEN(rayleigh_coefficient_wavelengths)); +} + +double get_weighted_rayleigh_scattering_length_snoman_d2o(void) +{ + /* Returns the weighted Rayleigh scattering length in D2O in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return rayleigh_scattering_length_d2o_weighted; +} + +double get_weighted_rayleigh_scattering_length_snoman_h2o(void) +{ + /* Returns the weighted Rayleigh scattering length in H2O in cm. */ + + if (!initialized) { + fprintf(stderr, "weighted rayleigh scattering length hasn't been initialized!\n"); + exit(1); + } + + return rayleigh_scattering_length_h2o_weighted; +} + double get_absorption_length_snoman_h2o(double wavelength) { /* Returns the absorption length in H2O in cm. */ diff --git a/src/optics.h b/src/optics.h index 06b89ca..64796c6 100644 --- a/src/optics.h +++ b/src/optics.h @@ -1,7 +1,15 @@ #ifndef OPTICS_H #define OPTICS_H -int optics_init(void); +#include "dict.h" + +extern char optics_err[256]; + +int optics_init(dict *db); +double get_rayleigh_scattering_length_snoman_d2o(double wavelength); +double get_rayleigh_scattering_length_snoman_h2o(double wavelength); +double get_weighted_rayleigh_scattering_length_snoman_d2o(void); +double get_weighted_rayleigh_scattering_length_snoman_h2o(void); double get_weighted_absorption_length_snoman_h2o(void); double get_absorption_length_snoman_h2o(double wavelength); double get_weighted_absorption_length_snoman_d2o(void); diff --git a/src/rsp_rayleigh.dat b/src/rsp_rayleigh.dat new file mode 100644 index 0000000..07a3576 --- /dev/null +++ b/src/rsp_rayleigh.dat @@ -0,0 +1,88 @@ +*LOG +*--- Standard titles file rsp_rayleight.dat +*--- This bank contains effective rayleigh scattering +*--- lengths in D2O and H2O for a +/- 10 ns window. +*LOG OFF +*DO RSPR 1 -i(30I / 1I 2F) +#. (This bank contains effective Rayleigh scattering lenghts) +#. +#. Contact: B.A. Moffat Queen's. +#. +#. History:- +#. ======= +#. +#. 4.0183 B.A. Moffat Rayleigh scattering for 10ns window. +#. D2O Rayleigh fraction kd = 0.83 +#. H2O Rayleigh fraction kh = 1.24 +#. D2O Atten. slope md = -2e-06 +#. H2O Atten. slope mh = 0 +#. Atten. slope scaled by (386/lambda)**4 +#. +#. Fri Mar 16 23:21:38 2001 +#. +#. Standard Database Header +#. +19750101 0 20380517 03331900 #. 1..4 Intrinsic validity + 0 0 0 #. 5..7 Data type, Task type, Format no. + 0 0 0 #. 8..10 Creation Date, Time, Source Id. +19750101 0 20380517 03331900 #. 11..14 Effective validity + 0 0 #. 15..16 Entry Date Time +4*0 #. 17..20 Spare +10*0 #. 21..30 Temporary data (not in database) +#. +#. End of Database Header +#. +#. User Data. +#.wl alpha_d2o alpha_h2o [1/cm] for 10 ns window. +220 2.6813e-04 4.4568e-04 +230 2.2291e-04 3.7070e-04 +240 1.8691e-04 3.1100e-04 +250 1.5794e-04 2.6296e-04 +260 1.3441e-04 2.2392e-04 +270 1.1513e-04 1.9192e-04 +280 9.9208e-05 1.6548e-04 +290 8.5955e-05 1.4347e-04 +300 7.4853e-05 1.2502e-04 +310 6.5494e-05 1.0946e-04 +320 5.7557e-05 9.6258e-05 +330 5.0792e-05 8.4999e-05 +340 4.4994e-05 7.5348e-05 +350 4.0003e-05 6.7036e-05 +360 3.5687e-05 5.9844e-05 +370 3.1938e-05 5.3595e-05 +380 2.8671e-05 4.8145e-05 +390 2.5811e-05 4.3373e-05 +400 2.3300e-05 3.9181e-05 +410 2.1088e-05 3.5485e-05 +420 1.9132e-05 3.2217e-05 +430 1.7398e-05 2.9318e-05 +440 1.5857e-05 2.6739e-05 +450 1.4483e-05 2.4439e-05 +460 1.3254e-05 2.2382e-05 +470 1.2154e-05 2.0538e-05 +480 1.1165e-05 1.8880e-05 +490 1.0275e-05 1.7388e-05 +500 9.4719e-06 1.6040e-05 +510 8.7460e-06 1.4821e-05 +520 8.0883e-06 1.3716e-05 +530 7.4913e-06 1.2713e-05 +540 6.9485e-06 1.1800e-05 +550 6.4539e-06 1.0968e-05 +560 6.0027e-06 1.0208e-05 +570 5.5902e-06 9.5136e-06 +580 5.2126e-06 8.8772e-06 +590 4.8663e-06 8.2934e-06 +600 4.5484e-06 7.7571e-06 +610 4.2560e-06 7.2636e-06 +620 3.9867e-06 6.8089e-06 +630 3.7384e-06 6.3894e-06 +640 3.5092e-06 6.0018e-06 +650 3.2973e-06 5.6434e-06 +660 3.1011e-06 5.3114e-06 +670 2.9193e-06 5.0036e-06 +680 2.7507e-06 4.7179e-06 +690 2.5940e-06 4.4524e-06 +700 2.4484e-06 4.2054e-06 +710 2.3128e-06 3.9754e-06 +720 2.1865e-06 3.7609e-06 + diff --git a/src/test-time-pdf.c b/src/test-time-pdf.c new file mode 100644 index 0000000..d0568f8 --- /dev/null +++ b/src/test-time-pdf.c @@ -0,0 +1,116 @@ +#include <stdio.h> +#include <unistd.h> +#include <stdlib.h> +#include <errno.h> /* for errno */ +#include <string.h> /* for strerror() */ +#include "likelihood.h" + +void usage(void) +{ + fprintf(stderr,"Usage: ./test-time-pdf\n"); + fprintf(stderr," -n number of points\n"); + fprintf(stderr," --mu-noise Expected number of noise PE\n"); + fprintf(stderr," --mu-indirect Expected number of PE from scattered/reflected light\n"); + fprintf(stderr," --mu-shower Expected number of PE from the electromagnetic shower\n"); + fprintf(stderr," --ts Expected time of direct light\n"); + fprintf(stderr," --ts-shower Expected time of shower light\n"); + fprintf(stderr," --t-mean Starting time for scattered/reflected light\n"); + fprintf(stderr," --t-sigma Standard deviation of direct light\n"); + fprintf(stderr," --ts-sigma Standard deviation of shower light\n"); + fprintf(stderr," -h Display this help message\n"); + exit(1); +} + +int main(int argc, char **argv) +{ + size_t i, n; + double *x; + double mu_noise, mu_indirect, mu_direct, mu_shower, ts, ts_shower, tmean, tsigma, ts_sigma; + + n = 1000; + mu_noise = 0.1; + mu_indirect = 0.5; + mu_direct = 1.0; + mu_shower = 1.0; + ts = 100.0; + ts_shower = 120.0; + tmean = 100.0; + tsigma = PMT_TTS; + ts_sigma = 10.0; + + for (i = 1; i < argc; i++) { + if (!strncmp(argv[i], "--", 2)) { + if (!strcmp(argv[i]+2,"mu-noise")) { + mu_noise = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"mu-indirect")) { + mu_indirect = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"mu-direct")) { + mu_direct = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"mu-shower")) { + mu_shower = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"ts")) { + ts = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"ts-shower")) { + ts_shower = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"t-mean")) { + tmean = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"t-sigma")) { + tsigma = strtod(argv[++i],NULL); + continue; + } else if (!strcmp(argv[i]+2,"ts-sigma")) { + ts_sigma = strtod(argv[++i],NULL); + continue; + } + } else if (argv[i][0] == '-') { + switch (argv[i][1]) { + case 'n': + n = atoi(argv[++i]); + break; + case 'h': + usage(); + default: + fprintf(stderr, "unrecognized option '%s'\n", argv[i]); + exit(1); + } + } + } + + x = malloc(n*sizeof(double)); + + for (i = 0; i < n; i++) { + x[i] = GTVALID*i/(n-1); + } + + FILE *pipe = popen("graph -T X --bitmap-size 2000x2000 -X 'Time (ns)' -Y Probability", "w"); + + if (!pipe) { + fprintf(stderr, "error running graph command: %s\n", strerror(errno)); + exit(1); + } + + for (i = 0; i < n; i++) { + fprintf(pipe, "%.10f %.10f\n", x[i], time_pdf(x[i],mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma)); + } + fprintf(pipe, "\n\n"); + + for (i = 0; i < n; i++) { + fprintf(pipe, "%.10f %.10f\n", x[i], time_cdf(x[i],mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma)); + } + fprintf(pipe, "\n\n"); + + if (pclose(pipe)) { + fprintf(stderr, "error closing graph command: %s\n", strerror(errno)); + exit(1); + } + + free(x); + + return 0; +} @@ -1227,6 +1227,139 @@ int test_norm_cdf(char *err) return 0; } +static double gsl_time_pdf(double x, void *params) +{ + double mu_noise, mu_indirect, mu_direct, mu_shower, ts, ts_shower, tmean, tsigma, ts_sigma; + + double *pars = (double *) params; + + mu_noise = pars[0]; + mu_indirect = pars[1]; + mu_direct = pars[2]; + mu_shower = pars[3]; + ts = pars[4]; + ts_shower = pars[5]; + tmean = pars[6]; + tsigma = pars[7]; + ts_sigma = pars[8]; + + return time_pdf(x,mu_noise,mu_indirect,&mu_direct,&mu_shower,1,&ts,&ts_shower,tmean,tsigma,&ts_sigma); +} + +int test_time_pdf_norm(char *err) +{ + /* Tests the time_pdf() function. */ + size_t i; + gsl_integration_cquad_workspace *w; + gsl_function F; + double result, error; + int status; + size_t nevals; + double params[9]; + double expected; + + w = gsl_integration_cquad_workspace_alloc(100); + + params[0] = 0.1; + params[1] = 0.5; + params[2] = 1.0; + params[3] = 1.0; + params[4] = 100.0; + params[5] = 120.0; + params[6] = 100.0; + params[7] = PMT_TTS; + params[8] = 10.0; + + F.function = &gsl_time_pdf; + F.params = params; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + params[0] = genrand_real2(); + params[1] = genrand_real2(); + params[2] = genrand_real2(); + params[3] = genrand_real2(); + + status = gsl_integration_cquad(&F, 0, GTVALID, 0, 1e-9, w, &result, &error, &nevals); + + if (status) { + sprintf(err, "error integrating time PDF: %s\n", gsl_strerror(status)); + goto err; + } + + expected = 1.0; + + if (!isclose(result, expected, 1e-2, 0)) { + sprintf(err, "integral of time_pdf = %.5f, but expected %.5f", result, expected); + return 1; + } + } + + return 0; + +err: + return 1; +} + +int test_time_cdf(char *err) +{ + /* Tests the time_cdf() function. */ + size_t i; + gsl_integration_cquad_workspace *w; + gsl_function F; + double result, error; + int status; + size_t nevals; + double params[9]; + double expected; + + w = gsl_integration_cquad_workspace_alloc(100); + + params[0] = 0.1; + params[1] = 0.5; + params[2] = 1.0; + params[3] = 1.0; + params[4] = 100.0; + params[5] = 120.0; + params[6] = 100.0; + params[7] = PMT_TTS; + params[8] = 10.0; + + F.function = &gsl_time_pdf; + F.params = params; + + init_genrand(0); + + for (i = 0; i < 100; i++) { + params[0] = genrand_real2(); + params[1] = genrand_real2(); + params[2] = genrand_real2(); + params[3] = genrand_real2(); + + double t = genrand_real2()*GTVALID; + + status = gsl_integration_cquad(&F, 0, t, 0, 1e-9, w, &result, &error, &nevals); + + if (status) { + sprintf(err, "error integrating time PDF: %s\n", gsl_strerror(status)); + goto err; + } + + expected = time_cdf(t,params[0],params[1],¶ms[2],¶ms[3],1,¶ms[4],¶ms[5],params[6],params[7],¶ms[8]); + + if (!isclose(result, expected, 1e-2, 0)) { + sprintf(err, "integral of time_pdf = %.5f, but expected %.5f", result, expected); + return 1; + } + } + + return 0; + +err: + return 1; +} + struct tests { testFunction *test; char *name; @@ -1263,6 +1396,8 @@ struct tests { {test_solid_angle_fast, "test_solid_angle_fast"}, {test_norm, "test_norm"}, {test_norm_cdf, "test_norm_cdf"}, + {test_time_pdf_norm, "test_time_pdf_norm"}, + {test_time_cdf, "test_time_cdf"}, }; int main(int argc, char **argv) |