diff --git a/ex-4/main.c b/ex-4/main.c new file mode 100644 index 0000000..778ea9e --- /dev/null +++ b/ex-4/main.c @@ -0,0 +1,117 @@ +#include +#include +#include +#include + +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; + } diff --git a/notes/images/dip.pdf b/notes/images/dip.pdf index 34d7677..2a6ba88 100644 Binary files a/notes/images/dip.pdf and b/notes/images/dip.pdf differ diff --git a/notes/sections/4.md b/notes/sections/4.md index 7a00642..60adde4 100644 --- a/notes/sections/4.md +++ b/notes/sections/4.md @@ -6,16 +6,16 @@ Consider a great number of non-interacting particles, each of which with a random momentum $\vec{p}$ with module between 0 and $p_{\text{max}}$ randomly angled with respect to a coordinate system {$\hat{x}$, $\hat{y}$, $\hat{z}$}. Once the polar angle $\theta$ is defined, the momentum vertical and horizontal -components of a particle, which will be referred as $p_v$ and $p_h$, are the -ones shown in @fig:components. +components of a particle, which will be referred as $\vec{p_v}$ and $\vec{p_h}$ +respectively, are the ones shown in @fig:components. If $\theta$ is evenly distributed on the sphere and the same holds for the -module $p$, which distribution will the average value of the absolute value of -$p_v$ follow as a function of $p_h$? +module $p$, which distribution will the average value of $|p_v|$ follow as a +function of $p_h$? \begin{figure} \hypertarget{fig:components}{% \centering -\begin{tikzpicture} +\begin{tikzpicture}[font=\scriptsize] % Axes \draw [thick, ->] (5,2) -- (5,8); \draw [thick, ->] (5,2) -- (2,1); @@ -41,178 +41,128 @@ $p_v$ follow as a function of $p_h$? } \end{figure} -The aim is to compute $\langle |p_v| \rangle (p_h) dp_h$. +Since the aim is to compute $\langle |p_v| \rangle (p_h)$, the conditional +distribution probability of $p_v$ given a fixed value of $p_h = x$ must first +be determined. It can be computed as the ratio between the probablity of +getting a fixed value of $P_v$ given $x$ over the total probability of $x$: -Consider all the points with $p_h \in [p_h ; p_h - dp_h]$: the values of -$p_v$ that these points can assume depend on $\theta$ and the total momentum -length $p$. +$$ + f (P_v | P_h = x) = \frac{f_{P_h , P_v} (x, P_v)} + {\int_{\{ P_v \}} d P_v f_{P_h , P_v} (x, P_v)} + = \frac{f_{P_h , P_v} (x, P_v)}{I} +$$ + +where $f_{P_h , P_v}$ is the joint pdf of the two variables $P_v$ and $P_h$ and +the integral runs over all the possible values of $P_v$ given $P_h$. + +This joint pdf can simly be computed from the joint pdf of $\theta$ and $p$ with +a change of variables. For the pdf of $\theta$ $f_{\theta} (\theta)$, the same +considerations done in @sec:3 lead to: + +$$ + f_{\theta} (\theta) = \frac{1}{2} \sin{\theta} \chi_{[0, \pi]} (\theta) +$$ + +whereas, being $p$ evenly distributed: + +$$ + f_p (p) = \chi_{[0, p_{\text{max}}]} (p) +$$ + +where $\chi_{[a, b]} (y)$ is the normalized characteristic function which value +is $1/N$ between $a$ and $b$, where $N$ is the normalization term, and 0 +elsewhere. Being a couple of independent variables, their joint pdf is simply +given by the product of their pdfs: + +$$ + f_{\theta , p} (\theta, p) = f_{\theta} (\theta) f_p (p) + = \frac{1}{2} \sin{\theta} \chi_{[0, \pi]} (\theta) + \chi_{[0, p_{\text{max}}]} (p) +$$ + +Given the new variables: $$ \begin{cases} - p_h = p \sin{\theta} \\ p_v = p \cos{\theta} + P_v = P \cos(\theta) \\ + P_h = P \sin(\theta) \end{cases} - \thus |p_v| = p |\cos{\theta}| = p_h \frac{|\cos{\theta}|}{\sin{\theta}} - = |p_v| (\theta) $$ -It looks like the dependence on $p$ has disappeared, but obviously it has -not. In fact, it lies beneath the limits that one has to put on the possible -values of $\theta$. For the sake of clarity, take a look at @fig:sphere (the -system is rotation-invariant, hence it can be drown at a fixed azimuthal angle). - -\begin{figure} -\hypertarget{fig:sphere}{% -\centering -\begin{tikzpicture} - % p_h slice - \definecolor{cyclamen}{RGB}{146, 24, 43} - \filldraw [cyclamen!15!white] (1.5,-3.15) -- (1.5,3.15) -- (1.75,3.05) - -- (2,2.85) -- (2,-2.85) -- (1.75,-3.05) - -- (1.5,-3.15); - \draw [cyclamen] (1.5,-3.15) -- (1.5,3.15); - \draw [cyclamen] (2,-2.9) -- (2,2.9); - \node [cyclamen, left] at (1.5,-0.3) {$p_h$}; - \node [cyclamen, right] at (2,-0.3) {$p_h + dp_h$}; - % Axes - \draw [thick, ->] (0,-4) -- (0,4); - \draw [thick, ->] (0,0) -- (4,0); - \node at (0.3,3.8) {$z$}; - \node at (4,0.3) {$hd$}; - % p_max semicircle - \draw [thick, cyclamen] (0,-3.5) to [out=0, in=-90] (3.5,0); - \draw [thick, cyclamen] (0,3.5) to [out=0, in=90] (3.5,0); - \node [cyclamen, left] at (-0.2,3.5) {$p_{\text{max}}$}; - \node [cyclamen, left] at (-0.2,-3.5) {$-p_{\text{max}}$}; - % Angles - \draw [thick, cyclamen, ->] (0,1.5) to [out=0, in=140] (0.55,1.2); - \node [cyclamen] at (0.4,2) {$\theta$}; - \draw [thick, cyclamen] (0,0) -- (1.5,3.15); - \node [cyclamen, above right] at (1.5,3.15) {$\theta_x$}; - \draw [thick, cyclamen] (0,0) -- (1.5,-3.15); - \node [cyclamen, below right] at (1.5,-3.15) {$\theta_y$}; - % Vectors - \draw [ultra thick, cyclamen, ->] (0,0) -- (1.7,2.2); - \draw [ultra thick, cyclamen, ->] (0,0) -- (1.9,0.6); - \draw [ultra thick, cyclamen, ->] (0,0) -- (1.6,-2); -\end{tikzpicture} -\caption{Momentum space at fixed azimuthal angle ("$hd$" stands for -"horizontal direction"). Some vectors with $p_h \in [p_h, p_h +dp_h]$ -are evidenced.}\label{fig:sphere} -} -\end{figure} - -As can be seen, $\theta_x$ and $\theta_y$ are the minimum and maximum tilts -angles of these vectors respectively, because if a point had $p_h \in [p_h; p_h -+ dp_h]$ and $\theta < \theta_x$ or $\theta > \theta_y$, it would have $p > -p_{\text{max}}$. Therefore their values are easily computed as follow: +with $\theta \in [0, \pi]$, the previous ones are written as: $$ - p_h = p_{\text{max}} \sin{\theta_x} = p_{\text{max}} \sin{\theta_y} - \thus \sin{\theta_x} = \sin{\theta_y} = \frac{p_h}{p_{\text{max}}} + \begin{cases} + P = \sqrt{P_v^2 + P_h^2} \\ + \theta = \text{atan}^{\star} ( P_h/P_v ) := + \begin{cases} + \text{atan} ( P_h/P_v ) &\incase P_v > 0 \\ + \text{atan} ( P_h/P_v ) + \pi &\incase P_v < 0 + \end{cases} + \end{cases} $$ -Since the average value of a quantity is computed by integrating it over all the -possible quantities it depends on weighted on their probability, one gets: +which can be shown having Jacobian: $$ - \langle |p_v| \rangle (p_h) dp_h = \int_{\theta_x}^{\theta_y} - d\theta P(\theta) \cdot P(p) \, dp \cdot |p_v| (\theta) + J = \frac{1}{\sqrt{P_v^2 + P_h^2}} $$ -where $d\theta P(\theta)$ is the probability of generating a point with $\theta -\in [\theta; \theta + d\theta]$ and $P(p) \, dp$ is the probability of -generating a point with $\vec{p}$ in the pink region in @fig:sphere, given a -fixed $\theta$. -The easiest to deduce is $P(p) \, dp$: since $p$ is evenly distributed, it -follows that: +Hence: -$$ - P(p) \, dp = \frac{1}{p_{\text{max}}} dp +$$ + f_{P_h , P_v} (P_h, P_v) = + \frac{1}{2} \sin[ \text{atan}^{\star} ( P_h/P_v )] + \chi_{[0, \pi]} (\text{atan}^{\star} ( P_h/P_v )) \cdot \\ + \frac{\chi_{[0, p_{\text{max}}]} \left( \sqrt{P_v^2 + P_h^2} \right)} + {\sqrt{P_v^2 + P_h^2}} $$ -with: +from which, the integral $I$ can now be computed. The edges of the integral +are fixed bt the fact that the total momentum can not exceed $P_{\text{max}}$: $$ - dp = p(p_h + dp_h) - p(p_h) - = \frac{p_h + dp_h}{\sin{\theta}} - \frac{p_h}{\sin{\theta}} - = \frac{dp_h}{\sin{\theta}} + I = \int + \limits_{- \sqrt{P_{\text{max}}^2 - P_h}}^{\sqrt{P_{\text{max}}^2 - P_h}} + dP_v \, f_{P_h , P_v} (x, P_v) $$ -hence +after a bit of maths, using the identity: $$ - P(p) \, dp = \frac{1}{p_{\text{max}}} \cdot \frac{1}{\sin{\theta}} \, dp_h + \sin[ \text{atan}^{\star} ( P_h/P_v )] = \frac{P_h}{\sqrt{P_h^2 + P_v^2}} $$ -For $d\theta P(\theta)$, instead, one has to do the same considerations done -in @sec:3, from which: +and the fact that both the characteristic functions are automatically satisfied +within the integral limits, the following result arises: $$ - P(\theta) d\theta = \frac{1}{2} \sin{\theta} d\theta + I = 2 \, \text{atan} \left( \sqrt{\frac{P_{\text{max}}^2}{x^2} - 1} \right) $$ -Ultimately, having found all the pieces, the integral must be computed: +from which: -\begin{align*} - \langle |p_v| \rangle (p_h) dp_h &= \int_{\theta_x}^{\theta_y} - d\theta \frac{1}{2} \sin{\theta} \cdot - \frac{1}{p_{\text{max}}} \frac{1}{\sin{\theta}} \, dp_h \cdot - p_h \frac{|\cos{\theta}|}{\sin{\theta}} - \\ - &= \frac{1}{2} \frac{p_h dp_h}{p_{\text{max}}} \int_{\theta_x}^{\theta_y} - d\theta \frac{|\cos{\theta}|}{\sin{\theta}} - \\ - &= \frac{1}{2} \frac{p_h dp_h}{p_{\text{max}}} \cdot \mathscr{O} -\end{align*} +$$ + f (P_v | P_h = x) = \frac{x}{P_v^2 + x^2} \cdot + \frac{1}{2 \, \text{atan} + \left( \sqrt{\frac{P_{\text{max}}^2}{x^2} - 1} \right)} +$$ -Then, with a bit of math: +Finally, putting all the pieces together, the average value of $|P_v|$ can now +be computed: -\begin{align*} - \mathscr{O} &= \int_{\theta_x}^{\theta_y} d\theta - \frac{|\cos{\theta}|}{\sin{\theta}} - \\ - &= \int_{\theta_x}^{\frac{\pi}{2}} d\theta \frac{\cos{\theta}}{\sin{\theta}} - - \int_{\frac{\pi}{2}}^{\theta_y} d\theta \frac{\cos{\theta}}{\sin{\theta}} - \\ - &= \left[ \ln{(\sin{\theta})} \right]_{\theta_x}^{\frac{\pi}{2}} - - \left[ \ln{(\sin{\theta})} \right]_{\frac{\pi}{2}}^{\theta_y} - \\ - &= \ln{(1)} -\ln{ \left( \frac{p_h}{p_{\text{max}}} \right) } - - \ln{ \left( \frac{p_h}{p_{\text{max}}} \right) } + \ln{(1)} - \\ - &= 2 \ln{ \left( \frac{p_{\text{max}}}{p_h} \right) } -\end{align*} - -\newpage -Hence, in conclusion: - -\begin{align*} - \langle |p_v| \rangle (p_h) dp_h &= \frac{1}{2} \frac{p_h dp_h}{p_{\text{max}}} - \cdot 2 \ln{ \left( \frac{p_{\text{max}}}{p_h} \right) } - \\ - &= \ln{ \left( \frac{p_{\text{max}}}{p_h} \right) } - \frac{p_h}{p_{\text{max}}} dp_h -\end{align*} +$$ + \langle |P_v| \rangle = \int + \limits_{- \sqrt{P_{\text{max}}^2 - P_h}}^{\sqrt{P_{\text{max}}^2 - P_h}} + f (P_v | P_h = x) = [ \dots ] = + x \, \frac{\ln \left( \frac{P_{\text{max}}}{x} \right)} + {\text{atan} \left( \sqrt{ \frac{P_{\text{max}}}{x^2} - 1} \right)} +$$ Namely: -\begin{figure} -\hypertarget{fig:plot}{% -\centering -\begin{tikzpicture} - \definecolor{cyclamen}{RGB}{146, 24, 43} - % Axis - \draw [thick, <->] (0,5) -- (0,0) -- (11,0); - \node [below right] at (11,0) {$p_h$}; - \node [above left] at (0,5) {$\langle |p_v| \rangle$}; - % Plot - \draw [domain=0.001:10, smooth, variable=\x, - cyclamen, ultra thick] plot ({\x},{12*ln(10/\x)*\x/10}); - \node [cyclamen, below] at (10,0) {$p_{\text{max}}$}; -\end{tikzpicture} -\caption{Plot of the expected distribution.}\label{fig:plot} -} -\end{figure} +![Plot of the expected distribution with + $P_{\text{max}} = 10$.](images/expected.pdf){#fig:plot} ## Monte Carlo simulation