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 probability of
getting a fixed value of $P_v$ given $x$ over the total probability of
getting that $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 $I$ runs over all the possible values of $P_v$ given a certain
$P_h$.  
$f_{P_h , P_v}$ can simply 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 can be written as:
$$
  \begin{cases}
    P = \sqrt{P_v^2 + P_h^2} \\
    \theta = \text{atan}_2 ( 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}_2 ( P_h/P_v )]
    \chi_{[0, \pi]} (\text{atan}_2 ( 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 by 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}_2 ( P_h/P_v )] = \frac{P_h}{\sqrt{P_h^2 + P_v^2}}
$$

and the fact that both the characteristic functions play no role 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 be
computed:
$$
  \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)}
$$ {#eq:dip}

Namely:

![Plot of the expected distribution with
  $P_{\text{max}} = 10$.](images/4-expected.pdf){#fig:plot}


## Monte Carlo simulation

The same distribution should be found by generating and binning points in a
proper way. A number of $N = 50000$ points were generated as a couple of values
($P$, $\theta$), with $P$ evenly distributed 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$ uniformly distributed between 0 and 1.  
The data binning turned out to be a bit challenging. Once $P$ was sampled and
$P_h$ was computed, the bin containing the latter's value must be found. If $n$
is the number of bins in which the range $[0, P_{\text{max}}]$ is divided into,
then the width $w$ of each bin is given by:
$$
  w = \frac{P_{\text{max}}}{n}
$$

and the $i^{th}$ bin in which $P_h$ goes in is:
$$
  i = \text{floor} \left( \frac{P_h}{w} \right)
$$

where 'floor' is the function which gives the bigger integer smaller than its
argument and the bins are counted starting from zero.  
Then, a vector in which the $j^{\text{th}}$ entry contains both the sum $S_j$
of all the $|P_v|$s relative to each $P_h$ fallen into the $j^{\text{th}}$ bin
itself and the number num$_j$ of the bin entries was iteratively 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 sampled couple the procedure is the
following. At first $S_j = 0 \wedge \text{num}_j = 0 \, \forall \, j$, then:

  - the couple $(P, \theta)$ is generated,
  - $P_h$ and $P_v$ are computed,
  - the $j^{\text{th}}$ bin containing $P_h$ is found,
  - num$_j$ is increased by 1,
  - $S_j$ is increased by $|P_v|$.

For $P_{\text{max}} = 10$ and $n = 50$, the following result was obtained:

![Sampled points histogram.](images/4-dip.pdf)

In order to check whether the expected distribution (@eq:dip) properly matches
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 inverse of the $\chi^2$ 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_r^2 = 0.071
$$

where $\chi_r^2$ is the $\chi^2$ per degree of freedom, proving a good
convergence.  
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{oss}}_{\text{max}}}
$$

where $\Delta P^{\text{oss}}_{\text{max}}$ is the $P^{\text{oss}}_{\text{max}}$
uncertainty. 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.

![Fitted sampled data. $P^{\text{oss}}_{\text{max}} = 10.005
  \pm 0.018$, $\chi_r^2 = 0.071$.](images/4-fit.pdf){#fig:fit}