analistica/notes/sections/4.md

Exercise 4

Kinematic dip PDF derivation

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 \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 |p_v| follow as a function of p_h?

\begin{figure} \hypertarget{fig:components}{% \centering \begin{tikzpicture}[font=\scriptsize] % Axes \draw [thick, ->] (5,2) -- (5,8); \draw [thick, ->] (5,2) -- (2,1); \draw [thick, ->] (5,2) -- (8,1); \node at (1.5,0.9) {$x$}; \node at (8.5,0.9) {$y$}; \node at (5,8.4) {$z$}; % Momentum \definecolor{cyclamen}{RGB}{146, 24, 43} \draw [ultra thick, ->, cyclamen] (5,2) -- (3.8,6); \draw [thick, dashed, cyclamen] (3.8,0.8) -- (3.8,6); \draw [thick, dashed, cyclamen] (5,7.2) -- (3.8,6); \draw [ultra thick, ->, pink] (5,2) -- (5,7.2); \draw [ultra thick, ->, pink] (5,2) -- (3.8,0.8); \node at (4.8,1.1) {$\vec{p_h}$}; \node at (5.5,6.6) {$\vec{p_v}$}; \node at (3.3,5.5) {$\vec{p}$}; % Angle \draw [thick, cyclamen] (4.4,4) to [out=80,in=210] (5,5); \node at (4.7,4.2) {$\theta$}; \end{tikzpicture} \caption{Momentum components.}\label{fig:components} } \end{figure}

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:

  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_v = P \cos(\theta) \\
    P_h = P \sin(\theta)
  \end{cases}

with \theta \in [0, \pi], the previous ones are written as:

  \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}

which can be shown having Jacobian:

  J = \frac{1}{\sqrt{P_v^2 + P_h^2}}

Hence:

  
  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}}

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}}:

  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)

after a bit of maths, using the identity:

  \sin[ \text{atan}^{\star} ( P_h/P_v )] = \frac{P_h}{\sqrt{P_h^2 + P_v^2}}

and the fact that both the characteristic functions are automatically satisfied within the integral limits, the following result arises:

  I = 2 \, \text{atan} \left( \sqrt{\frac{P_{\text{max}}^2}{x^2} - 1} \right)

from which:

  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)}

Finally, putting all the pieces together, the average value of |P_v| can now be computed:

\begin{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^2_{\text{max}}}{x^2} - 1} \right)} \end{align*}

Namely:

{#fig:plot}

Monte Carlo simulation

The same distribution should be found by generating and binning points in a proper way.
A number of N = 50'000 points were generated as a couple of values (p, \theta), with p evenly distrinuted between 0 and p_{\text{max}} and \theta given by the same procedure described in @sec:3, namely:

  \theta = \text{acos}(1 - 2x)

with x uniformely distributed between 0 and 1.
The data binning turned out to be a bit challenging. Once p was sorted and p_h was computed, the bin in which the latter goes in must be found. If n is the number of bins in which the range [0, $p_{\text{max}}$] is willing to be divided into, then the width w of each bin is given by:

  w = \frac{p_{\text{max}}}{n}

and the j^{th} bin in which p_h goes in is:

  j = \text{floor} \left( \frac{p_h}{w} \right)

where 'floor' is the function which gives the upper integer lesser than its argument and the bins are counted starting from zero.
Then, a vector in which the j^{\text{th}} entry contains the sum S_j of all the $|p_v|$s relative to each p_h fallen into the j^{\text{th}} bin and the number num$_j$ of entries in the j^{\text{th}} bin was reiteratively updated. At the end, the average value of |p_v|_j was computed as $S_j / \text{num}_j$.
For the sake of clarity, for each sorted couple, it works like this:

• the couple [p; \theta] is generated;
• p_h and p_v are computed;
• the j^{\text{th}} bin in which p_h goes in is computed;
• num$_j$ is increased by 1;
• S_j (which is zero at the beginning of everything) is increased by a factor |p_v|.

At the end, \langle |p_h| \rangle_j = \langle |p_h| \rangle (p_h) where:

  p_h = j \cdot w + \frac{w}{2} = w \left( 1 + \frac{1}{2} \right)

The following result was obtained:

In order to check wheter the expected distribution properly metches the produced histogram, a chi-squared minimization was applied. Being a simple one-parameter fit, the \chi^2 was computed without a suitable GSL function and the error of the so obtained estimation of p_{\text{max}} was given as the inverese of the \chi^3 second derivative in its minimum, according to the Cramér-Rao bound.
The following results were obtained:

  p^{\text{oss}}_{\text{max}} = 10.005 \pm 0.018 \with \chi^2 = 0.071

In order to compare p^{\text{oss}}_{\text{max}} with the expected value p_{\text{max}} = 10, the following compatibility t-test was applied:

  p = 1 - \text{erf}\left(\frac{t}{\sqrt{2}}\right)\ \with
  t = \frac{|p^{\text{oss}}_{\text{max}} - p_{\text{max}}|}
           {\Delta p_{\text{max}}}

where \Delta p_{\text{max}} is the absolute error of p_{\text{max}}. At 95% confidence level, the values are compatible if p > 0.05. In this case:

• t = 0.295
• p = 0.768

which allows to assert that the sampled points actually follow the predicted distribution. In @fig:fit, the fit function superimposed on the histogram is shown.

{#fig:fit}