# Physics 77 Final Project

---

## 1. Introduction

### 1.1. Background information

In 1964, Peter Higgs, along with five other scientists proposed the Higgs mechanism, a way that some particles can acquire mass. In their theories, a particle known as a scalar boson should also exist. This particle was called the Higgs boson, and its existence was believed to be a strong indication that the Higgs field was the correct explanation.

It's not easy to find the Higgs boson in the nature, since it decays so rapidly that it cannot be picked up directly by any detector so far. However, based on quantum mechanics, researchers are able to predict the decay products of the Higgs boson. In other words, certain combination of particles that matches the prediction may indicate the existence of the Higgs boson. There are multiple possible processes, resulting in different products. And here we are going to focus on one particular process, which is $H \rightarrow \gamma \gamma$ decay.

In the experiment carried out by the CERN, two beams of high energy particles were accelerated and collided into each other, and the diphotons events produced in this collision are recorded. Not all of them are corresponding to the $H \rightarrow \gamma \gamma$ event, as diphotons can be produced by many other random processes. (In fact, only a very few of them are produced by the decay of the Higgs boson.) And the probability distribution of the invariant mass should be uniformed for those random processes. Hence, If there is an anomalous peak in the distribution, we verify the existence of the Higgs boson and those data points are corresponding to the Higgs boson decay.

### 1.2. Description of the objectives

In this project, we are going to generate the mass distribution of diphoton system, using the collisions' data from LHC. Then fit the diphoton mass distribution to determine if a signal exists by performing hypothesis testing and determining the significance of the signal.

Once the existence of the signal is verified, it will be visualized by several plots. And our group will also use the data to create the sonification of the signal. The signal strength of the $H \rightarrow \gamma \gamma$ will also be measured.

Lastly, using these data of the diphoton system, we can reconstruct the information of the original Higgs boson by energy and momentum conservation. The direction of the the flight of the Higgs boson will be determined and visualized in an animation.

### 1.3. Brief Description of methods used

To determine the existence of the signal, we are going to use hypothesis testing method. The data will be fitted under different hypotheses through a NLLR function minimization process. The results of the fit will be either used to generate pseudo-experiments or to calculate Profile Likelihood Ratio so that the statistical significance of the different hypotheses can be determined.

Similar minimization method will be used again on the signal stregth measurement to obtain its uncertaties.

The programming language used for this project is python. Several python modules, including numpy, matplotlib.pyplot, scipy etc., are used to conduct task regarding file management, mathematical operations, plotting, fitting and so on.

## 2. Method

### 2.1. Data sample

#### 2.1.1 Data structure

In this project, we are using the data of collisions carried out by LHC from 2015 to 2018. They are stored in a h5 file. The shape of data array is (1178902, 10), indicating that there are 1,178,902 collision events that contain two photons. The axis 1 has 10 entries. They are

- transverse momentum of photon 1
- pseudo rapidity of photon 1
- azimuthal angle of photon 1
- energy of photon 1
- transverse momentum of photon 2
- pseudo rapidity of photon 2
- azimuthal angle of photon 2
- energy of photon 2
- Event Number, which is an index of the collision event
- Run Number, which is an index of a `run`. At LHC, the detector is often run for an extended period of time, ranging from a few minutes to a few hours. Data events collected in the same data taking period are said to be in the same `run`.

Data is retrieved and stored in two numpy arrays, for particle 1 and particle 2 repectively, in our program for further manipulations. The last two columns (event number and run number) are not necessary since they will not be the focus of our project.

#### 2.1.2 Mass distribution

To calculate the mass of the diphoton, we need to first calculate the momentum in three directions using the formula below:

$$
p_x = p_T cos(\phi)\\
p_y = p_T sin(\phi)\\
p_z = p_T sinh(\eta)
$$

By the law of energy and momentum conservation, the energy and momentum of the diphoton will be

$$
E_{\gamma\gamma} = E_1+E_2\\
\vec{p_{\gamma\gamma}} = \vec{p_1}+\vec{p_2}
$$
From energy and momentum, we can get the mass of the diphoton
$$
m_{\gamma\gamma} = \sqrt{E_{\gamma\gamma}^2 - |\vec{p_{\gamma\gamma}}|^2}
$$

The result can be visualized in a histogram. 55 bins are set in the mass range 105-160 GeV. The horizontal coordinate represents the mass while the vertical coordinate represents the number of events observed.

> figure 1 Mass distribution

As shown in the plot, a bump around 125 GeV is already discernible. The event in this bump is believed to be produced by the decay of the Higgs boson and the bump itself is therefore a signal for the Higgs boson. Follow up, we will perform hypothesis testing on this signal and quantify the plausibility of it being a Higgs boson statistically.

### 2.2. Probability density functions under different hypotheses

To carry out the test, we need to first fit the mass distribution for different hypothesis.

For background-only hypotheses, we assume that all the data collected are from processes that were already known. The probability distribution of background processes is a curve inversely proportional to mass since it's easier to produce photons with low energy and it's harder to produce photons with higher energy. Therefore, a 4th order polynomial function is used to fit the mass distribution for this hypothesis.

For signal-plus-background hypothesis, we use a normal distribution with a mean of 125 and a standard deviation of 1.6 to approximate the bump. The fit for this hypothesis will be in the form of the sum of this normal distribution and a 4th order polynomial used for background. The function looks like
$$
f(x=m_{\gamma\gamma}) = n_b \sum_i c_i x^i + n_s \mathrm{Normal}(x,125,1.6)
$$
where $n_b$ and $n_s$ are two coefficients of normalization.

Here, we are using maximum likelihood estimation(MLE) to obtain parameters in the fit. As the name implies, MLE proceeds to maximize a likelihood function, which in turn maximizes the agreement between the model and the data. The product of MLE is generally very small, so it is normally replaced by a log-likelihood function or, in this case, a negative log-likelihood function. Maximizing the log-likelihood or minimizing the negative function yields the same results. The negative log-likelihood for both hypotheses are defined as follows:

$$
NLL_{b-only}(obs|b-only)=-\sum_i \log(Poission(k_i, \lambda_i))\\
NLL_{s+b}(obs|s+b)=-\sum_i \log(Poission(k_i, \lambda_i))
$$
where $k_i$ and $\lambda_i$ refer to observed value from data and expected value from fitted functions.

Since maximizing the log-likelihood and minimizing the negative function yields the same results, so we can then fit our data for two hypotheses by utilizing the minimization function in python scipy module.

The output is unlikely to be satisfactory if we only do the fit procedure once. The output will become more accurate as the time we run the code increases. Usually, this iterative process is done manually until the output value becomes stable. But in our program, we define a new function that controls the iteration instead of repeating it by hands. This function can determine the end of the optimization, following the criteria of which the difference between last two NLLvalue is less than 0.1.(0.1 is also an adjustable number.) In the code, we use a while loop to achieve this goal. Besides the coefficient array, the function will also provide how many times the minimization is used.

### 2.3 Signal significance

Statistical significance refers to the claim that a set of observed data are not the result of chance but can instead be attributed to a specific cause. In order to determine whether we can reject the null hypothesis, which is the background-only hypothesis in this case, we need to calculate the statistical significance for both of the hypotheses.

One way to calculate the significance is by calculating the p-value. Firstly, based on the fitted PDF we obtain in previous part, we can use random number generator in python to run pseudo experiments following the Poisson distribution. Then determine the number of pseudo-experiments generated under two different hypotheses that have a NLLR value less than the NLLR value of the observed $m_{\gamma\gamma}$ distribution. NLLR stands for negative log-likelihood ratio and it is defined as
$$
NLLR = \frac{NLL_{s+b}}{NLL_{b-only}} = -2 (\log L_{s+b}-\log L_b)
$$
The p-value is then the ratio of the above number to the total number of pseudo-experiments, and it can be converted to the significance using the norm.ppf method of scipy.stats.

Besides, Profile Likelihood Ratio (PLR) can also be used for significance calculation. It is defined as:
$$
PLR = \frac{NLL_{s+b}(obs|s+b)}{NLL_{b-only}(obs|b-only)}
$$
The numerator in this ratio come from the signal+background fit to the observed data while the denominator come from the background-only fit. Hence no more pseudo-experiments are needed for this method. Follow the formula, the significance $Z$ will be
$$
Z = \sqrt{PLR}
$$

### 2.4 Signal strength

For signal+background hypothesis, recall that a function in the form
$$
f_{s+b}(m_{\gamma\gamma}) = n_b \sum_i c_i x^i + n_s \mathrm{Normal}(x,125,1.6)
$$
is used to fit the observed data. In this formula, the strength of the signal is governed by a normalization factor $n_s$. In order to find its uncertainties, we can re-parameterize the PDF by introducing a $\mu$ factor in front of the $n_s$:
$$
f_{s+b}(m_{\gamma\gamma}) = n_b \sum_i c_i x^i + \mu n_s \mathrm{Normal}(x,125,1.6)
$$
and perform a scan of the NLL value as a function of $\mu$ and calculate their corresponding $-2\Delta LL$ value, which is defined as
$$
-2\Delta LL = NLLR(\mu)-\min(NLLR({\mu}))
$$

The $\mu$ value that minimizes $-2\Delta LL$ should be the central value of the measured $\mu$, and we denote it as $\hat{\mu}$. The two $\mu$ values that give a $-2\Delta LL$ value that differs from the minimum by 1 should define the $\pm1\sigma$ uncertainties. The $+1 \sigma$ and $-1 \sigma$ uncertainties of the $\mu$ measurement are then $\mu_{lo} - \hat{\mu}$ and $\mu_{hi} - \hat{\mu}$, respectively.

It can be presented in the format:
$$
\mu = \hat{\mu}^{+ \mu_{hi} - \hat{\mu}}_{- \hat{\mu} - \mu_{lo}}
$$

There are two different ways to find the $\mu_{lo}$ and $\mu_{hi}$. One is to search the maximum/minimum $\mu$ on the left/right side of the centeral value that satisfies $-2\Delta LL>1$. The other way to determine the $\mu$ is to find two cloest points to 1 on both sides by setting the absolute value of $-2\Delta LL - 1$ to be lower than a adjustable constant. In our code, we choose the second method due to the accuracy.

### 2.5 Visualization of the Higgs boson's flight direction

First of all, we build a model 3D sphere on which the values would be reflected as the pixel colors, and it can rotate around the axis of $\theta = 90^{\circ}$.

To begin with, we have specified various parameters that will be used when drawing a 3D sphere. Then, we normalize the data, that is, define the relative value according to the maximum and minimum values in the three-dimensional array, so that it could range from 0 to 1.0. After that, we match the normalized size of the data with the color. We only retain one decimal place, because it is difficult for the human eye to distinguish more accurate colors. At the same time, more accuracy means more computation, which will increase the load of the computer. We also carried out the transformation between spherical coordinates and rectangular coordinates, so that the data can be mapped to the surface of the sphere. We use ax.plot_surface function to draw a 3D sphere model. The animation. ArtistAnimation function is also used to realize the rotation around the axis of $\theta = 90^{\circ}$.

Secondly, we process the data in the topic. Previously we have used the data of the diphotons system to calculate the energy and momentum components in x, y, z directions of the Higgs boson. And flight directions of Higgs bosons can be determined by these momentum component in the following relationship:

$$
\theta = \arccos(\frac{p_z}{\sqrt{p_x^2+p_y^2+p_z^2}})\\
\\
\phi = \begin{cases}
\arccos(\frac{p_x}{\sqrt{p_x^2+p_y^2}}),\quad &x>0 \\
2\pi-\arccos(\frac{p_x}{\sqrt{p_x^2+p_y^2}}),\quad &x\leq0
\end{cases}
$$

But not all of the data in the data set are corresponding to $H \rightarrow \gamma \gamma$ decay. So We first need to subtracted the background in the data sample with the help of a simulated Monte Carlo sample. From the required component we could find that the signal of the Higgs boson exists in the 120GeV to 130 GeV interval. So we choose interval [120, 130] GeV to be the signal window, leaving other parts to be side bands. Since side bands is dominated by background events, we therefore can use them to predict the background events in the signal window. The normalized background data distribution can be obtained by multiplying the data distribution between 120GeV and 130GeV by $SF$, defined as:
$$
SF = \frac{n_b^{SW}}{n_b^{SBs}}
$$
where $n_b^{SW}$ and $n_b^{SBs}$ represent the number of background events in the signal window and side bands.

We then draw the distribution map, obtained by subtracting the background plus signal data distribution from the background data distribution, of the Higgs boson's energy and quantity in the graph with $\phi$ as the abscissa and $\theta$ as the ordinate, where every 5.6 degrees is divided into one pixel. The relationship between the flight directions and average energy, the total energy, and the number of Higgs bosons in one flight direction were calculated respectively.

Finally, we substitute the processed data into the constructed function, so that the value and coordinate relationship can be displayed through the surface color. We could also obtain an animated sphere rotating around the axis of $\theta = 90^{\circ}$.

### 2.6 Sonification of the signal

We came up with to ways to translate the signal into audio. We can either treat data points as individual notes or discrete points on a waveform. And the data points used in both method is the residual of the $m_{\gamma \gamma}$ distribution. It is defined as:
$$
R(m_{\gamma \gamma})= f_{obs}(m_{\gamma \gamma})-f_{bkg, s+b}(m_{\gamma \gamma})
$$
This function reflects the distance between the observed data and our expected background value in signal+background hypothesis. By doing so, we hope to extract the signal from the background noises.

#### 2.6.1 Method 1: Waveform

This is relatively straight forward. After obtaining the residual points, we can use the Audio package from IPython to convert the points into an audio easily.

#### 2.6.2 Method 2: Individual notes

In this method, we want to consider each point on the residual plot as a note, so that the data set can be translated into a melody. The transformation from experimental data to sound is performed by the following procedures:

Firstly, we linearly normalize the residuals to the range [0, 1], by mapping the maximum to 1 and the minimum to 0.

Secondly, we divide the interval [0, 1] equally into 12 parts and map the residuals in each part to 12 different notes.

Thirdly, we use the principle of Additive Synthesis to generate the sound of Percussion instruments: We first separate the timeline into frames and generate waves by calculating the value of the wave in every single frame. Then, we determine the fundamental frequency by the note and calculate the wave using the given frequency spectrum of the instrument and Fourier Transform. Finally, we use the ADSR rule to scale the amplitude of the wave in different frames and multiply it by the volume.

At last, we put the wave into a .wav file using the "wave" module in python.

### 2.7 Expected output

#### 2.7.1 Signal significance and signal strength

As shown in part 2.1.2, a bump around 125 GeV is already discernible. Hence we expect that the signal+background hypothesis will give us a better fit to the data. Consequently, for pseudo-experiments generated under different hypotheses, their NLLR distribution should also be different and the NLLR of the observation should be closer to the s+b NLLR distribution. In other words, the result is statistically significant.

For the signal strength $\mu$, the NLLR plot is expected to have a parabolic shape with its central value to be exactly 1, since $n_s$ have been normalized already. The intersections of the y = 1 line and the NLLR curve will give us the lower and upper 1$\sigma$ deviation of the signal.

#### 2.7.2 Visualization of the Higgs boson's flight direction

The result will be an animation of a 3D sphere, on which the values would be reflected as the pixel colors, and it can rotate around the axis of $\theta = 90^{\circ}$.

#### 2.7.3 Sonification of the signal

For method 1, the result audio clip is expected to be the sound of a "chord", which is basically a certain combination of notes. Each note in the chord corresponds to a frequency in the waveform obtained by the Fourier Transform.

For method 2, the result audio clip will be a melody since each point in the residual plots corresponds to a certain note.

## 3. Result / Interpretation

### 3.1 Fitting result

The result of the fit under background-only hypothesis is shown as below:
>figure 2 Backgroud-only Fit

And the fitted function is:

$$
f(x) = 4.31168393 \times 10^{-4}(1+1.64361191 \times 10^7 x-2.98726710 \times 10^5x^2 \\+1.85270718 \times 10^3 x^3 -3.89281905 x^4)
$$

The PDF function (red dashed line) fits the distribution of the data well within the intervals (105, 120) and (130, 160), but does not accurately describe the distribution within the (120, 130) signal window.

And the result of the fit under signal+background hypothesis is shown as below:
>figure 3 Signal+background fit

And the fitted function is:
$$
f(x) = 4.60326967 \times 10^{-4}(1+1.66859428 \times 10^7 \cdot x-3.08848670 \times 10^5 \cdot x^2\\+1.95039671 \times 10^3 \cdot x^3 -4.17125072 \cdot x^4) + 5.42065526 \times 10^3 \cdot \mathrm{Normal}(x,125,1.6)
$$
Unlike previous fit, here the fitted function matches the data distribution throughout the entire interval. The results of the fit have tentatively shown that the signal + background hypothesis is more reliable. And next we will quantify its credibility by measuring the significance and plotting out their NLLR distributions.

### 3.2 Significance

#### 3.2.1 Result

We run 10,000 pseudo-experiments for both hypotheses and their NLLR distribution are shown as below:

> figure 4 NLLR distribtuion of pseudo-experiments, b-only hypothesis

> figure 5 NLLR distribution of pseudo-experiments, s+b hypothesis

The first plot is back-ground only hypothesis. The green line in the plot shows the observed NLLR, which equals to -152.9418472751786 in our calculation, and it is far away from the distribution generated by pseudo-experiments.

The second plot is signal+background hypothesis. The green line is still the observed NLLR and now it lies in the middle of the distribution. And all the bins that are less than the observed NLLR is marked in orange.

> figure 6 NLLR plot

They can be shown in the same figure:

Since the observed NLLR doesn't cut through the background-only hypothesis, we cannot use its p-value to calculate the significance. Therefore, the second method mentioned in 2.3, which we only need to calculate the PLR of the observation and take its square root to obtain the significance, can be deployed. The calculated result for significance Z is 12.367.

#### 3.2.2 Interpretation

The certainty of the Higgs boson particle's existence was based on the 5*σ* criterion, which corresponds to a *p*-value of about 1 in 3.5 million. If the background-only hypothesis is rejected at five standard deviations or beyond, we claim evidence of the a signal is established.

Here our obtained Z value is greater than the $5\sigma$, so there must be a signal in our data, and the signal+background hypothesis is valid.

### 3.3 Signal strength

#### 3.3.1 Result
By performing a scan, the $-2\Delta LL$ - $\mu$ function can be plotted out as below:

> figure 7 Signel strengh plot

As shown, its central value is exactly 1, which corresponds to the normalized coefficient of the signal we obtained in the previous s+b fit. And the lower and upper bound are 0.9312 and 1.0689 respectively, giving a $\pm$ 0.07 error if we round it to two decimal places.

So the signal strength factor is given by:
$$
\mu =1.0^{+0.07}_{-0.07}
$$

#### 3.3.2 Interpretation

For a normal distribution, the probability of falling into the $\pm1\sigma$ interval is around 68%. In this case, we can say that there is a 68% chance to have the signal falling in the interval (0.9312, 1.0689).

### 3.4 Visualization

#### 3.4.1 Result

As discussed in previous part, a normalization coefficient $SF$ is needed. In our calculation, $SF$ approximately equals to 0.25.

We ended up generating two sets of animations. In one set, different colors is used to reflect the total energy of the Higgs boson in the corresponding direction, and in the other set, colors is used to reflect the number of Higgs bosons in the corresponding direction. In the animation, the sphere is slowing rotating around its central axis. Alternatively, the sphere can be viewed in all directions by dragging with the mouse.

Unfortunately, due to the memory limitations of the jupyter notebook, the animation cannot be played in the notebook. We will however upload the results we generate in vscode as gif files.

#### 3.4.2 Interpretation

Below are two demos of the animation we get. These are just two screenshots and the animation will be attached to the submission.

> figure 8 Animation demo 1: total energy

> figure 9 Animation demo 2: total quantities

In the first picture, the color reflects the energy of the Higgs boson. Reflected by the color, we can tell that in the polar region, the total energy of Higgs bosons in this direction is significantly larger than other directions. In other regions of the sphere, the color is mostly uniformly distributed, which indicates that the total energy of Higgs bosons in other direction is comparable.

Similar conclusion can be derived from the second picture, yet the colors in this picture reflects the total numbers of Higgs bosons. We can tell that there are significantly less Higgs bosons observed in the directions corresponding to polar regions, comparing to other regions on the sphere. And other than polar regions, the color is mostly uniformly distributed, which indicates that comparable amounts of Higgs bosons are observed in the most of directions.

### 3.5 Sonification

#### 3.5.1 Result

For method 1, the outcome can be parameterized by two parameters.

The first parameter which need to be adjusted is the number of bins, that is, the number of data points, N. Specifically,
$$
T=\frac{N}{f_{sample}}
$$
where $T$ is the duration of the sound, $f_{sample}$ is the sampling frequency.

f_sample should be in the range of 20 Hz~20 kHz to be audible for human being. If N=55, f = 20Hz ~ 20kHz, T may be too short. Our goal is to balance the three variables to generate a tune with a proper frequency, duration and appropriate number of data points that a laptop can handle. Here we set N = 10000, f_sample = 3000, and the duration is 3 seconds. Then we recalculate the negative log likelihood and minimize it to obtain the new fit coefficients for signal plus background hypothesis. The coefficients are [7.83552487e+06 -1.45541532e+05 9.21799887e+02 -1.97655175e+00 5.59045891e-06 2.95912299e+01]. The residual is therefore the difference between the 10000 points of observation and the 10000 predicted number of background events given by the above signal-plus-background fit. The scatter figure is shown below.

> figure 10 Residual plot, 10000 bins

The name of this audio clip is *Sound of Higgs 1* and it is in the attachment.

For method 2, since points will be tranlated into notes, so we will stick to our previous 55 bins calculation. As a result, the residual plot is shown as below:

> figure 11 Residual plot, 55 bins

And we get a audio clip (named *Sound of Higgs 2* in the attachment) around 13 seconds (It's more than 13s but the information are in the first 13s.) from these residual points.

#### 3.5.2 Interpretation

Since there is a large amount onumbers in this chord, so the result, to some extend will stund more like a clip of "noise".

For method 1, since there is a large amount of notes in this chord, so the result, to some extend will stund more like a clip of "noise". Using the Audio package from IPython, one can easily convert the points into an audioor f_sample = 3000, it sounds like howling wind. As f_sample increases, the sound grows shrill. The .wav files is uploaared together on bCourse. However the result can't really give us the infomation arbout the bump in the residual plot, since the Fourier Transform doesn't preserve the shape of the original function.

This problem is solved in method 2, however. Since the data are translated into a melody and the pitch of notes in the melody reflects the value of the function. Hence, the ascending and descending variation of the melody can reflect the convexity of the function. In this example, the bump of the signal generates a peak in the melodic direction, and its auditory trend sounds like changes from bass to treble and back to bass.

This melodic peak appears at about the fourth second in the audio. We can hear a section consisted of hiing pitch notes, which corresponds to the highest point of the Higgs boson signal. This peak is preceded and followed by two short upward and downward scales. Outside the signal window, we can hear mainly oscillations that are close in pitch and dominated by low tones, corresponding to the signal of the background event.

## 4. Conclusion

In this project, for the required component, we calculate the NLLR of observation by fitting the b-only and signal plus background hypotheses to data. $NLLR_{obs} =
-152.94$.
In order to determine how much we can reject the b-only hypothesis and accept the signal-plus-background hypothesis, we try to generate pseudo experiments but failed to get a meaningful pvalue and significance value. Then we turn to use PLR which requires the signal strength. The signal strength $\mu =1.0^{+0.07}_{-0.07}$ is obtained via minimizing the negative log likelihood $-2\Delta LL$. Finally, we calculate the significance $Z = \sqrt{PLR} = 12.38$, which indicates that the siganl exisna.

For the selective component, we use the signal verfied in iprevious part to generate animations and audios cliips In nimations, we map the flight direnctions f Higgs bosons onto a sphere and visualize and the average, total energy, number of particles into color blocks on the sphere. In the sonifications, we covert the signal into both a chord and a melody.

## 5. Highlights and Potential Improvement
### 5.1 Highlights
- In the required component, we define a new function to determine when the output is satisfactory in the minimization process. Thus the fitting process can be more accurate and automatic.
- Compare the two ways to find the $\mu$ which makes $-2\Delta LL$ minimized.
- For the second selective components (sonifcation), we propose two approaches from different perspectives, based on Fourier transform and music theory respectively. The former method gives a single note with a specific harmonic frequency. In the contrast, the latter generates a melody.
### 5.2 Potential improvements
- Increase the number of bins
- For the expected value calculated during the fitting process, we use the central value to approximate the number of events in the bin. We can improve the accuracy by doing the integral over that bin.
- We use the forth degree polynomials to approximate our data in the fitting process. To imporve theaoccuracy, higher degree of polynomials might be appled during the fit.