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

ex-4: renamed file dip.c into main.c as you wish

Browse Source
This commit is contained in:
Giù Marcer 2020-04-23 23:59:19 +02:00 committed by rnhmjoj
parent 87beaa2248
commit ad224b825b
1 changed files with 0 additions and 117 deletions

117
ex-4/dip.c
View File

@ -1,117 +0,0 @@
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_rng.h>
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.
// 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]);
else
{
fprintf(stderr, "Usage: %s -[hiIntp]\n", argv[0]);
fprintf(stderr, "\t-h\tShow this message.\n");
fprintf(stderr, "\t-N integer\tThe number of events to generate.\n");
fprintf(stderr, "\t-n integer\tThe number of bins of the histogram.\n");
fprintf(stderr, "\t-p float\tThe maximum value of momentum.\n");
return EXIT_FAILURE;
}
}
// Initialize an RNG.
//
gsl_rng_env_setup();
gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
// Generate the angle θ uniformly distributed on a sphere using the
// inverse transform:
//
// θ = acos(1 - 2X)
//
// where X is a random uniform variable in [0,1), and the module p of
// the vector:
//
// p² = p_v² + p_h²
//
// uniformly distributed between 0 and p_max. The two components are
// then computed as:
//
// p_v = p⋅cos(θ)
// p_h = p⋅sin(θ)
//
// The histogram will be updated this way.
// The j-th bin where p_h goes in is given by:
//
// step = p_max / n
// j = floor(p_h / step)
//
// 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.
//
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.
//
double step = p_max / n;
struct bin *b;
double theta;
double p;
double p_v;
double p_h;
size_t j;
for (size_t i = 0; i < N; i++)
{
// Generate the event.
//
theta = acos(1 - 2*gsl_rng_uniform(r));
p = p_max * gsl_rng_uniform(r);
// Compute the components.
//
p_v = p * cos(theta);
p_h = p * sin(theta);
// Update the histogram.
//
j = floor(p_h / step);
b = &histo[j];
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);
for (size_t i = 0; i < n; i++)
{
histo[i].sum = histo[i].sum / histo[i].amo;
printf("\n%.5f", histo[i].sum);
};
// free memory
gsl_rng_free(r);
free(histo);
return EXIT_SUCCESS;
}