Giu/analistica
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

ex-4: fix multi-line comment warnings

Browse Source
This commit is contained in:
Michele Guerini Rocco 2020-05-30 09:22:29 +00:00
parent 28664215e8
commit 681b1b62b8
1 changed files with 94 additions and 114 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

208
ex-4/main.c
View File

@ -7,18 +7,15 @@
#include <gsl/gsl_deriv.h> #include <gsl/gsl_deriv.h>
// Process CLI arguments. /* Process CLI arguments and print usage */
// int parser(size_t *N, size_t *n, double *p_max,
int parser(size_t *N, size_t *n, double *p_max, char argc, char **argv, size_t *go) size_t argc, char **argv, size_t *go) {
{ for (size_t i = 1; i < argc; i++) {
for (size_t i = 1; i < argc; i++)
{
if (!strcmp(argv[i], "-n")) *N = atol(argv[++i]); if (!strcmp(argv[i], "-n")) *N = atol(argv[++i]);
else if (!strcmp(argv[i], "-b")) *n = atol(argv[++i]); else if (!strcmp(argv[i], "-b")) *n = atol(argv[++i]);
else if (!strcmp(argv[i], "-p")) *p_max = atof(argv[++i]); else if (!strcmp(argv[i], "-p")) *p_max = atof(argv[++i]);
else if (!strcmp(argv[i], "-o")) *go = 1; else if (!strcmp(argv[i], "-o")) *go = 1;
else else {
{
fprintf(stderr, "Usage: %s -[hnbpo]\n", argv[0]); fprintf(stderr, "Usage: %s -[hnbpo]\n", argv[0]);
fprintf(stderr, "\t-h\tShow this message.\n"); fprintf(stderr, "\t-h\tShow this message.\n");
fprintf(stderr, "\t-n N\tThe number of particles to generate. (default: 50000)\n"); fprintf(stderr, "\t-n N\tThe number of particles to generate. (default: 50000)\n");
@ -26,71 +23,63 @@ int parser(size_t *N, size_t *n, double *p_max, char argc, char **argv, size_t *
fprintf(stderr, "\t-p PMAX\tThe maximum value of momentum. (default: 10)\n"); fprintf(stderr, "\t-p PMAX\tThe maximum value of momentum. (default: 10)\n");
fprintf(stderr, "\t-o\tPrint histogram to stdout.\n"); fprintf(stderr, "\t-o\tPrint histogram to stdout.\n");
return 0; return 0;
}
} }
return 1;
} }
return 1;
}
int main(int argc, char **argv) int main(int argc, char **argv) {
{
// Set default options. // Set default options.
//
size_t N = 50000; size_t N = 50000;
size_t n = 50; size_t n = 50;
double p_max = 10; double p_max = 10;
size_t go = 0; size_t go = 0;
// Get eventual CLI arguments. // Get eventual CLI arguments.
//
int res = parser(&N, &n, &p_max, argc, argv, &go); int res = parser(&N, &n, &p_max, argc, argv, &go);
if (go == 0 && res == 1) if (go == 0 && res == 1) {
{ printf("Generating histogram with:\n"
printf("\nGenerating histogram with:\n" "%ld points\n"
"%ld points\n" "%ld bins\n"
"%ld bins\n" "p_max = %.3f\n\n", N, n, p_max);
"p_max = %.3f\n\n", N, n, p_max); }
} else if (res == 0)
else if (res == 0) return EXIT_FAILURE; return EXIT_FAILURE;
// printf("step: \t%.5f\n", step);
// Initialize an RNG. // Initialize an RNG.
//
gsl_rng_env_setup(); gsl_rng_env_setup();
gsl_rng *r = gsl_rng_alloc(gsl_rng_default); gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
// Generate the angle θ uniformly distributed on a sphere using the /* Generate the angle θ uniformly distributed on a sphere using the
// inverse transform: * inverse transform:
// *
// θ = acos(1 - 2X) * θ = acos(1 - 2X)
// *
// where X is a random uniform variable in [0,1), and the module p of * where X is a random uniform variable in [0,1), and the module p of
// the vector: * the vector:
// *
// p² = p_v² + p_h² * p² = p_v² + p_h²
// *
// uniformly distributed between 0 and p_max. The two components are * uniformly distributed between 0 and p_max. The two components are
// then computed as: * then computed as:
// *
// p_v = p⋅cos(θ) * p_v = p⋅cos(θ)
// p_h = p⋅sin(θ) * p_h = p⋅sin(θ)
// *
// The histogram is updated this way. * The histogram is updated this way.
// The j-th bin where p_h goes in is given by: * The j-th bin where p_h goes in is given by:
// *
// step = p_max / n * step = p_max / n
// j = floor(p_h / step) * j = floor(p_h / step)
// *
// Thus an histogram was created and a structure containing the number of * 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 * entries in each bin and the sum of |p_v| in each of them is created and
// filled while generating the events (struct bin). * filled while generating the events (struct bin).
// */
struct bin *histo = calloc(n, sizeof(struct bin)); struct bin *histo = calloc(n, sizeof(struct bin));
// Some useful variables. // Some useful variables.
//
double step = p_max / n; double step = p_max / n;
struct bin *b; struct bin *b;
double theta; double theta;
@ -99,46 +88,40 @@ int main(int argc, char **argv)
double p_h; double p_h;
size_t j; size_t j;
for (size_t i = 0; i < N; i++) for (size_t i = 0; i < N; i++) {
{ // Generate the event
// Generate the event.
//
theta = acos(1 - 2*gsl_rng_uniform(r)); theta = acos(1 - 2*gsl_rng_uniform(r));
p = p_max * gsl_rng_uniform(r); p = p_max * gsl_rng_uniform(r);
// Compute the components. // Compute the components
//
p_v = p * cos(theta); p_v = p * cos(theta);
p_h = p * sin(theta); p_h = p * sin(theta);
// Update the histogram. // Update the histogram
//
j = floor(p_h / step); j = floor(p_h / step);
b = &histo[j]; b = &histo[j];
b -> amo++; b -> amo++;
b -> sum += fabs(p_v); b -> sum += fabs(p_v);
} }
// Compute the mean value of each bin and print it to stodut /* Compute the mean value of each bin and print it to stodut
// together with other useful things to make the histogram. * together with other useful things to make the histogram.
// */
if (go == 1) if (go == 1) {
{
printf("bins: \t%ld\n", n); printf("bins: \t%ld\n", n);
printf("step: \t%.5f\n", step); printf("step: \t%.5f\n", step);
} }
for (size_t i = 0; i < n; i++) for (size_t i = 0; i < n; i++) {
{
histo[i].sum = histo[i].sum / histo[i].amo; // Average P_v histo[i].sum = histo[i].sum / histo[i].amo; // Average P_v
if (go == 1) printf("\n%.5f", histo[i].sum); if (go == 1) printf("\n%.5f", histo[i].sum);
}; };
// Compare the histigram with the expected function: /* Compare the histigram with the expected function:
// *
// x * log(p_max/x)/arctan(sqrt(p_max^2/x^2 - 1)) * x * log(p_max/x)/arctan(sqrt(p_max^2/x^2 - 1))
// *
// using the χ² test. * using the χ² test.
// */
struct parameters params; struct parameters params;
params.histo = histo; params.histo = histo;
params.n = n; params.n = n;
@ -151,8 +134,7 @@ int main(int argc, char **argv)
double min_p = p_max - 5; double min_p = p_max - 5;
double max_p = p_max + 5; double max_p = p_max + 5;
// Initialize minimization. // Initialize minimization
//
double x = p_max; double x = p_max;
int max_iter = 100; int max_iter = 100;
double prec = 1e-7; double prec = 1e-7;
@ -161,71 +143,69 @@ int main(int argc, char **argv)
gsl_min_fminimizer *s = gsl_min_fminimizer_alloc(T); gsl_min_fminimizer *s = gsl_min_fminimizer_alloc(T);
gsl_min_fminimizer_set(s, &func, x, min_p, max_p); gsl_min_fminimizer_set(s, &func, x, min_p, max_p);
// Minimization. // Minimization
// for (int iter = 0; status == GSL_CONTINUE && iter < max_iter; iter++) {
for (int iter = 0; status == GSL_CONTINUE && iter < max_iter; iter++)
{
status = gsl_min_fminimizer_iterate(s); status = gsl_min_fminimizer_iterate(s);
x = gsl_min_fminimizer_x_minimum(s); x = gsl_min_fminimizer_x_minimum(s);
min_p = gsl_min_fminimizer_x_lower(s); min_p = gsl_min_fminimizer_x_lower(s);
max_p = gsl_min_fminimizer_x_upper(s); max_p = gsl_min_fminimizer_x_upper(s);
status = gsl_min_test_interval(min_p, max_p, 0, prec); status = gsl_min_test_interval(min_p, max_p, 0, prec);
} }
double result = x; double result = x;
double res_chi = chi2(result, &params); double res_chi = chi2(result, &params);
if (go == 0) if (go == 0) {
{
printf("Results:\n"); printf("Results:\n");
printf("χ² = %.3f\n", res_chi); printf("χ² = %.3f\n", res_chi);
printf("p_max = %.3f\n", result); printf("p_max = %.3f\n", result);
} }
// Compute the second derivative of χ² in its minimum for the result error. /* Compute the second derivative of χ² in its minimum for the result error.
// *
// p_max = α * p_max = α
// * 2
// (Ei - Oi)² * (Eᵢ - Oᵢ)
// χ² = Σi ---------- * χ² = Σi ──────────
// Ei * Eᵢ
// *
// / Oi² \ * ⎛ 2⎞
// ∂αχ² = Σi | 1 - --- | ∂αE * ⎜ Oᵢ ⎟
// \ Ei² / * ∂αχ² = Σi Eᵢ⋅⎜1 - ───⎟⋅∂αE
// * ⎜ 2⎟
// / Oi² / Oi² \ \ * ⎝ Eᵢ ⎠
// ∂²αχ² = Σi | (∂αE)² 2 --- + ∂²αE | 1 - --- | | *
// \ Ei³ \ Ei² / / * ⎛ 2 ⎛ 2⎞ ⎞
// * ⎜ Oᵢ ⎜ Oᵢ ⎟ ⎟
* ∂²αχ² = Σi ⎜(∂αE)²⋅2⋅─── + ∂²αE ⋅⎜1 - ───⎟ ⎟
* ⎜ 2 ⎜ 2⎟ ⎟
* ⎝ Eᵢ ⎝ Eᵢ ⎠ ⎠
*/
double expecto, A, B; double expecto, A, B;
double error = 0; double error = 0;
for (size_t i = 0; i < n; i++) for (size_t i = 0; i < n; i++) {
{
x = (i + 0.5) * step; x = (i + 0.5) * step;
expecto = expected(x, result); expecto = expected(x, result);
A = 2 * pow(exp1d(x, result) * histo[i].sum / expecto, 2); A = 2 * pow(exp1d(x, result) * histo[i].sum / expecto, 2);
B = exp2d(x, result) * (1 - pow((histo[i].sum / expecto), 2)); B = exp2d(x, result) * (1 - pow((histo[i].sum / expecto), 2));
error = error + A + B; error = error + A + B;
}; };
error = 1/error; error = 1/error;
if (go == 0) printf("ΔP_max = %.3f\n", error); if (go == 0)
printf("ΔP_max = %.3f\n", error);
// Check compatibility // Check compatibility
//
double t = fabs(result - p_max)/error; double t = fabs(result - p_max)/error;
double alpha = 1 - erf(t/sqrt(2)); double alpha = 1 - erf(t/sqrt(2));
if (go == 0) if (go == 0) {
{
printf("\nCompatibility:\n"); printf("\nCompatibility:\n");
printf("t = %.3f\n", t); printf("t = %.3f\n", t);
printf("α = %.3f\n\n", alpha); printf("α = %.3f\n\n", alpha);
} }
// Free memory. // Free memory
//
gsl_min_fminimizer_free(s); gsl_min_fminimizer_free(s);
gsl_rng_free(r); gsl_rng_free(r);
free(histo); free(histo);
return EXIT_SUCCESS; return EXIT_SUCCESS;
} }