From a0c1f4560325311558eedf39bd191ddaad182367 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Tue, 27 Nov 2018 16:07:56 -0600 Subject: add rayleigh scattering This commit adds Rayleigh scattering to the likelihood function. The Rayleigh scattering lengths come from rsp_rayleigh.dat from SNOMAN which only includes photons which scattered +/- 10 ns around the prompt peak. The fraction of light which scatters is treated the same in the likelihood as reflected light, i.e. it is uniform across all the PMTs in the detector and the time PDF is assumed to be a constant for a fixed amount of time after the prompt peak. --- src/.gitignore | 1 + src/Makefile | 4 +- src/fit.c | 11 ++++- src/likelihood.c | 39 ++++++++++----- src/likelihood.h | 4 ++ src/optics.c | 112 +++++++++++++++++++++++++++++++++++++++++- src/optics.h | 10 +++- src/rsp_rayleigh.dat | 88 +++++++++++++++++++++++++++++++++ src/test-time-pdf.c | 116 +++++++++++++++++++++++++++++++++++++++++++ src/test.c | 135 +++++++++++++++++++++++++++++++++++++++++++++++++++ 10 files changed, 503 insertions(+), 17 deletions(-) create mode 100644 src/rsp_rayleigh.dat create mode 100644 src/test-time-pdf.c (limited to 'src') 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 diff --git a/src/fit.c b/src/fit.c index 5bbe1e6..23eab74 100644 --- a/src/fit.c +++ b/src/fit.c @@ -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 /* 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 /* for gsl_strerror() */ #include #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 +#include +#include +#include /* for errno */ +#include /* 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; +} diff --git a/src/test.c b/src/test.c index dddf235..a3f98a9 100644 --- a/src/test.c +++ b/src/test.c @@ -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) -- cgit