Add function to generate individual angles from moliere distribution

[Jean1995]

This PR adds a function from which one can sample individual angles from the moliere distribution, as well as a function to sample a two-dimensional version of the angle.

[Jean1995]

I have created a short script where one can see the difference between the "old" and "new" way:

import proposal as pp
import matplotlib.pyplot as plt
import numpy as np
import time

particle = pp.particle.EMinusDef()
medium = pp.medium.Air()
cuts = pp.EnergyCutSettings(np.inf, 1, False)
cross = pp.crosssection.make_std_crosssection(particle, medium, cuts, False)

displacement = pp.make_displacement(cross, False)

moliere = pp.scattering.MoliereInterpol(particle, medium)

E_i = 1
E_f = 0.9
grammage = displacement.solve_track_integral(E_i, E_f)

STAT = 5e5

angle_list_old = []

start_time = time.time()
for i in range(int(STAT)):
    deflection = moliere.scatter(grammage, E_i, E_f, np.random.rand(4))
    initial_direction = pp.Cartesian3D(0, 0, 1)
    _, final_direction = pp.scattering.scatter_initial_direction(initial_direction, deflection)
    theta = np.arccos(initial_direction * final_direction / (initial_direction.magnitude * final_direction.magnitude))
    angle_list_old.append(theta)
end_time = time.time()

print(end_time - start_time)

angle_list_new = []

start_time = time.time()
for i in range(int(STAT)):
    angle = moliere.scattering_angle_2d(grammage, E_i, E_f, np.random.rand(), np.random.rand())
    initial_direction = pp.Cartesian3D(0, 0, 1)
    initial_direction.deflect(np.cos(angle), 0)
    theta = np.arccos((initial_direction * pp.Cartesian3D(0, 0, 1))/initial_direction.magnitude)
    angle_list_new.append(theta)
end_time = time.time()
print(end_time - start_time)

min_bin = min(min(angle_list_old), min(angle_list_new))
max_bin = max(max(angle_list_old), max(angle_list_new))
bins = np.linspace(0, 1.1*np.pi, 20)

plt.hist(angle_list_old, bins=bins, label='Old implementation', histtype='step', color='tab:orange')
plt.hist(angle_list_new, bins=bins, label='New implementation', histtype='step', color='tab:blue')
plt.axvline(x = np.mean(angle_list_old), label='Mean old', linestyle='dashed', color='tab:orange')
plt.axvline(x = np.mean(angle_list_new), label='Mean new', linestyle='dashed', color='tab:blue')
plt.legend()
plt.xlabel('deflection angle / rad')
plt.yscale('log')

Which creates the plot

Note that this is for low energies (1 MeV, electrons)