118 lines
2.9 KiB
C
118 lines
2.9 KiB
C
|
#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;
|
||
|
|
}
|