Multiple scattering in array of pixel detectors

PHYS291 project: Multiple scattering through an array of pixel detectors

Description
A c++ program which simulates electrons passing through an array of silicon pixel detectors, and undergoing multiple scattering inside the silicon. Other interactions are assumed to be negligible.

Motivation
When a pixel detector has been created, it needs to be tested before it can be put into use. We need to verify that when the detector detects a hit, the hit is detected in the correct location. The proper approach to doing this is to have a prediction as to which pixel ought to light up, so that we can compare this prediction to the real data coming from the pixel detector in question. One way to do this is by sandwiching the pixel detector in between other pixel detectors (ones which are known to be working properly). When a particle passes through the array of detectors, we can reconstruct its track using the data from the surrounding pixel detectors. This track gives us our prediction for which pixel ought to light up, and we can compare that to which pixel actually lights up.
There is, however, a problem: Particles passing through matter dont necessarily move in straight lines. Therefore, the reconstructed track might be different from the actual track. In my project, I will create a Monte-Carlo simulation of the multiple-scattering process which happens in the system of detectors, so that I can visualize how large a deviation away from the straight line we might expect.
In the future (outside the scope of the PHYS291 project), I also plan to use this simulation to test a track reconstruction algorithm that I will be making. This ultimately ties into my masters thesis about track reconstruction.

The program
MonteCarlo.zip
After unzipping, the program should be ran using ROOT. It must be named "MonteCarlo.cpp".

Virtual experiment set-up
I want my simulation to resemble a real-life setup, so I will use some specs which I found in the article 'Performance of the EUDET-type beam telescopes'
I will assume that the detectors are made entirely of silicon, and my source for the physical properties of silicon is the particle data group.

The physical properties are as follows: 7 silicon pixel detectors of 54 micrometer thickness. They are placed 20 millimeters apart. Electrons with momentum of 4 GeV/c will be passing through the detector. The silicon has a radiation depth of 9.37 cm.

Theory
The scattering angles in my simulation will be based on theory presented in section 33.3 of this pdg article .
This section presents formulae (33.22c together with 33.16) which allows us to generate deflection angles by using normal-distributed random numbers as input. This is particularly convenient, because we can get normal-distributed random numbers from the standard library in C++. The formula approximates scatterings towards the x-axis and y-axis as independent (for particle moving in the z-direction), so we will need to generate two different angles for each scattering process: one for scattering in the xz-plane, and one for scattering in the yz-plane.
There is one additional parameter I would like to introduce, which is not mentioned in the pdg article. The scattering angle depends on the thickness of the pixel detector, but the article does not take the incident angle of the particle into account. In order to account for the incident angle (theta) of the particle, the thickness of the pixel detector is multiplied by 1/cos(theta).

For my program, I decided to reverse-engineer the "Vector3" class used with the Unity game engine. This is simply because I've become used to working with it, and it was a good way of practicing some C++ syntax.
Originally, I wanted to have the velocity of each particle being elegantly represented by a single 3-dimensional vector, but this ran into a problem when I wanted to apply rotations to the vector: If you apply a rotation of A degrees in the xz-plane onto the vector (0,0,1), and then project that vector onto the xz-plane, you will observe that the angle between the projected vector and the Z-axis is indeed A degrees. However, if you then were to apply a second rotation of B degrees in the yz-plane to the vector, and then again project that vector onto the xz-plane, you will see that the angle between the projected vector and the z-axis is no longer equal to A. In other words, the second rotation interfered with the first one. In order to solve this problem, my program uses two velocity vectors: The velocity direction in the xz-plane, and the velocity direction in the yz-plane. These will be put together to form the proper velocity direction, so that the next hit location and the incident angle can be calculated. This is done by scaling the vectors so that both their z-components are equal to 1, and then simply grabbing the x-component from the xz-plane, the y-component from the yz-plane, and z = 1. When the resulting vector is multiplied by the distance between two silicon detectors, it points all the way from one hit location to the next. I call this vector Delta in the program.
The pixel detectors are of course arranged along the z-axis, with their surface normals parallel to the z-axis. Incident angles can then be found by a simple trick using the trigonometric definition of the dot product:
$\vec{Delta} \cdot \vec{z} = |Delta||z|cos(\theta) \rightarrow \theta = acos( \frac{\vec{Delta} \cdot \vec{z}}{|Delta|})$
Rotations can be defined using the following rotation matrices:
$R_{xz} = \begin{bmatrix} cos\theta & 0 & sin\theta\\ 0 & 1 & 0 \\ -sin\theta & 0 & cos\theta \end{bmatrix} , R_{yz} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & cos\theta & -sin\theta \\ 0 & sin\theta & cos\theta \end{bmatrix}$

Implementation
All particles in the simulation will start with a velocity of (0,0,1), (0,0,1), an incident angle of 0, and position (0,0,0). Then, the following process will repeat for each pixel detector in the array: Two scattering angles are generated; one for the xz-plane, and one for the yz-plane. These scattering angles will depend on the properties of the immediate pixel detector, as well as the incident angle of the particle. The angles are applied to the velocity vectors using rotation matrices. The program will then combine the velocity vectors, calculate their incident angle, and then the next hit location.
The electrons are initially moving from the origin along the z-axis. The distance from the ideal hit location to the actual hit location is therefore equal to the distance from the actual hit location to the z-axis. This value will be called delta_r, and is calculated using the Pythagorean theorem. Delta_r is calculated for each of the seven pixel detectors, and filled into seven different histograms. However, the first histogram should turn out pretty boring, because all the delta_rs should be equal to zero. (No scatterings have happened yet in the first detector, so all the electrons hit the bulls eye.)

Results / plots
I obtain the following histograms for delta_r:

As you can see, the particles tend to land farther and farther from the origin for each subsequent pixel detector and scattering. However, they arent too far away. The vast majority of the particles are within one nanometer of the ideal situation.
A funny thing about the method I chose, with the x-axis and y-axis being mostly independent, is that there is almost no chance of having zero scattering. The histograms all tend towards zero in the origin. This is because of the fact that we generate two different random numbers for each scattering process, and the likelihood of both being zero is very low.
If we make histograms for delta only along the x-axis, so that one of the two random numbers is suppressed, we see that scattering angles of 0 is indeed possible:

Finally, lets look at histograms of only the incident angles:

This looks very similar in shape to the delta_r histograms, which is to be expected. We do, however, also see that these histograms have a jagged texture to them. Given that these are very small numbers, I suspect that this might be related to rounding errors, or to the random number generator.

References
Jansen, H., Spannagel, S., Behr, J. et al. Performance of the EUDET-type beam telescopes. EPJ Techn Instrum 3, 7 (2016). https://doi.org/10.1140/epjti/s40485-016-0033-2
Passage of particles through matter: http://pdg.lbl.gov/2019/reviews/rpp2019-rev-passage-particles-matter.pdf
Particle data group: silicon http://pdg.lbl.gov/2019/AtomicNuclearProperties/HTML/silicon_Si.html