diff --git a/ex-4/lib.c b/ex-4/lib.c new file mode 100644 index 0000000..c6bc6bd --- /dev/null +++ b/ex-4/lib.c @@ -0,0 +1,86 @@ +#include "lib.h" + + +// Minimization wrapper. +// +double chi2(double p_max, void* params) + { + // Pars parameters. + struct parameters p = *((struct parameters*) params); + struct bin* histo = p.histo; + size_t n = p.n; + double step = p.step; + + // Compute χ². + double chi = 0; + double expecto = 0; + for (size_t i = 0; i < n; i++) + { + expecto = expected((i + 0.5) * step, p_max); + chi += pow(histo[i].sum - expecto, 2)/expecto; + }; + return chi; + } + + +// Expected function. +// +double expected (double x, double p_max) + { + // When p_max < x, the argument under sqrt is negative. Return great number + // to lead away the minimization. + if (p_max < x) + { + double num = 500; + return num; + } + return x * log(p_max/x)/atan(sqrt(fabs(pow(p_max, 2)/pow(x,2) - 1))); + } + + +// First derivative of the expected function. +// +double exp1d (double x, double p_max) + { + // When p_max < x, the argument under sqrt is negative. Return great number + // to lead away the minimization. + if (p_max < x) + { + double num = 500; + return num; + } + double a = p_max; + double A = sqrt(fabs(pow(a, 2)/pow(x,2) - 1)); + double B = atan(A); + double C = log(a/x); + return x/(a*B) - x*C/(a*A*pow(B, 2)); + } + +// +// Second derivative of the expected function. +// +double exp2d (double x, double p_max) + { + // When p_max < x, the argument under sqrt is negative. Return great number + // to lead away the minimization. + if (p_max < x) + { + double num = 500; + return num; + } + double a = p_max; + double a2 = pow(a, 2); + + double A = sqrt(fabs(a2/pow(x,2) - 1)); + double B = atan(A); + double C = log(a/x); + double D = 1 - a2/pow(x,2); + + double E = -2/(a2*D*B); + double F = 1/(a2*A); + double G = 1/(pow(x, 2)*pow(A, 3)); + double H = 1/a2 + 2/(a2 * A * B); + + return x*((E + F + G) * C/B - H)/B; + } + diff --git a/ex-4/lib.h b/ex-4/lib.h new file mode 100644 index 0000000..26c9ac7 --- /dev/null +++ b/ex-4/lib.h @@ -0,0 +1,39 @@ +#pragma once + +#include +#include + +// Histogram struct. +// +struct bin + { + size_t amo; // Amount of events in the bin. + double sum; // Sum of |p_v|s of all the events in the bin. + }; + +// Minimization struct. +// +struct parameters + { + struct bin* histo; // Histogram + size_t n; // Number of bins + double step; // Bin width + }; + +///////////////////////////////////////////////////////////////// + +// Minimization wrapper. +// +double chi2(double p_max, void* params); + +// Expected function. +// +double expected (double x, double p_max); + +// First derivative of the expected function. +// +double exp1d (double x, double p_max); + +// Second derivative of the expected function. +// +double exp2d (double x, double p_max); diff --git a/ex-4/main.c b/ex-4/main.c index 778ea9e..8dda0ec 100644 --- a/ex-4/main.c +++ b/ex-4/main.c @@ -1,24 +1,21 @@ +#include "lib.h" #include #include #include #include +#include +#include -int main(int argc, char **argv) + +// Process CLI arguments. +// +int parser(size_t *N, size_t *n, double *p_max, char argc, char **argv) { - - // Set default options. - // - size_t N = 50000; // number of events. - size_t n = 50; // number of bins. - double p_max = 10; // maximum value of momentum module. - - // Process CLI arguments. - // for (size_t i = 1; i < argc; i++) { - if (!strcmp(argv[i], "-N")) N = atol(argv[++i]); - else if (!strcmp(argv[i], "-n")) n = atol(argv[++i]); - else if (!strcmp(argv[i], "-p")) p_max = atof(argv[++i]); + if (!strcmp(argv[i], "-N")) *N = atol(argv[++i]); + else if (!strcmp(argv[i], "-n")) *n = atol(argv[++i]); + else if (!strcmp(argv[i], "-p")) *p_max = atof(argv[++i]); else { fprintf(stderr, "Usage: %s -[hiIntp]\n", argv[0]); @@ -29,6 +26,25 @@ int main(int argc, char **argv) return EXIT_FAILURE; } } + return EXIT_SUCCESS; + } + + +int main(int argc, char **argv) + { + + // Set default options. + // + size_t N = 50000; // number of events. + size_t n = 50; // number of bins. + double p_max = 10; // maximum value of momentum module. + int res = parser(&N, &n, &p_max, argc, argv); + if (res == 0) printf("\nGenerating histogram with:\n" + "%ld points\n" + "%ld bins\n" + "p_max = %.3f\n\n", N, n, p_max); + + // printf("step: \t%.5f\n", step); // Initialize an RNG. // @@ -51,7 +67,7 @@ int main(int argc, char **argv) // p_v = p⋅cos(θ) // p_h = p⋅sin(θ) // - // The histogram will be updated this way. + // The histogram is updated this way. // The j-th bin where p_h goes in is given by: // // step = p_max / n @@ -59,13 +75,8 @@ int main(int argc, char **argv) // // Thus an histogram was created and a structure containing the number of // entries in each bin and the sum of |p_v| in each of them is created and - // filled while generating the events. + // filled while generating the events (struct bin). // - struct bin - { - size_t amo; // Amount of events in the bin. - double sum; // Sum of |p_v|s of all the events in the bin. - }; struct bin *histo = calloc(n, sizeof(struct bin)); // Some useful variables. @@ -94,22 +105,96 @@ int main(int argc, char **argv) // j = floor(p_h / step); b = &histo[j]; - b->amo++; - b->sum += fabs(p_v); + b -> amo++; + b -> sum += fabs(p_v); } // Compute the mean value of each bin and print it to stodut // together with other useful things to make the histogram. // - printf("bins: \t%ld\n", n); - printf("step: \t%.5f\n", step); + // printf("bins: \t%ld\n", n); + // printf("step: \t%.5f\n", step); for (size_t i = 0; i < n; i++) { - histo[i].sum = histo[i].sum / histo[i].amo; - printf("\n%.5f", histo[i].sum); + histo[i].sum = histo[i].sum / histo[i].amo; // Average P_v + //printf("\n%.5f", histo[i].sum); }; + + // Compare the histigram with the expected function: + // + // x * log(p_max/x)/arctan(sqrt(p_max^2/x^2 - 1)) + // + // using the χ² test. + // + struct parameters params; + params.histo = histo; + params.n = n; + params.step = step; - // free memory + gsl_function func; + func.function = &chi2; + func.params = ¶ms; + + double min_p = 5; + double max_p = 15; + + // Initialize minimization. + // + double x = 10; + int max_iter = 100; + double prec = 1e-7; + int status; + const gsl_min_fminimizer_type *T = gsl_min_fminimizer_brent; + gsl_min_fminimizer *s = gsl_min_fminimizer_alloc(T); + gsl_min_fminimizer_set(s, &func, x, min_p, max_p); + + // Minimization. + // + for (int iter = 0; status == GSL_CONTINUE && iter < max_iter; iter++) + { + status = gsl_min_fminimizer_iterate(s); + x = gsl_min_fminimizer_x_minimum(s); + min_p = gsl_min_fminimizer_x_lower(s); + max_p = gsl_min_fminimizer_x_upper(s); + status = gsl_min_test_interval(min_p, max_p, 0, prec); + } + + double result = x; + printf("p_max: %.7f\n", result); + + // Compute the second derivative of χ² in its minimum for the result error. + // + // p_max = α + // + // (Ei - Oi)² + // χ² = Σi ---------- + // Ei + // + // / Oi² \ + // ∂αχ² = Σi | 1 - --- | ∂αE + // \ Ei² / + // + // / Oi² / Oi² \ \ + // ∂²αχ² = Σi | (∂αE)² 2 --- + ∂²αE | 1 - --- | | + // \ Ei³ \ Ei² / / + // + double expecto, A, B; + double error = 0; + for (size_t i = 0; i < n; i++) + { + x = (i + 0.5) * step; + expecto = expected(x, result); + A = 2 * pow(exp1d(x, result) * histo[i].sum / expecto, 2); + B = exp2d(x, result) * (1 - pow((histo[i].sum / expecto), 2)); + error = error + A + B; + }; + error = 1/error; + printf("ΔP_max: %.7f\n\n", error); + + + // Free memory. + // + gsl_min_fminimizer_free(s); gsl_rng_free(r); free(histo); diff --git a/makefile b/makefile index 80898d4..9df4aea 100644 --- a/makefile +++ b/makefile @@ -21,7 +21,7 @@ ex-3/bin/main: ex-3/main.c ex-3/common.c ex-3/likelihood.c ex-3/chisquared.c $(CCOMPILE) ex-4: ex-4/bin/main -ex-4/bin/main: ex-4/main.c +ex-4/bin/main: ex-4/main.c ex-4/lib.c $(CCOMPILE) ex-5: ex-5/bin/casino ex-5/bin/manual ex-5/bin/trifecta