Study of resonant double Higgs production in the NMSSM model with $\gamma\gamma +b\bar{b}$γγ+bbˉ as the final state, using MadGraph5
Abstract
After the discovery of the Standard Model (SM) Higgs Boson 'h' in 2012, there are ongoing rigourous studies to understand the couplings of h with fermions, massive vector bosons (Z, W ^{+}, W ^{}), massless bosons (Photons ( $\gamma$γ), gluons) and with itself. The crosssection of the SM double Higgs (hh) production process is quite low and hence, difficult to detect. We therefore look at Beyond the Standard Model (BSM) double Higgs production to understand the selfcoupling which predicts an increased crosssection. One of the most popular Supersymmetry (SUSY) models – the NexttoMinimal Supersymmetric Standard Model (NMSSM) predicts the existence of seven Higgs bosons. Typically, the Higgs bosons in the NMSSM are denoted in an order with increasing masses  H_{1}, H_{2}, ...., H_{7}. In the special case of a charge parity  conserving Higgs potential, we have three CP even Higgs bosons  H_{1}, H_{2}, H_{3}, two CP odd ones  A_{1}, A_{2}, and a pair of charged Higgs bosons  H^{+}, H^{ }. In this project, it is assumed that H_{1} is the SMlike Higgs boson. The H_{1}H_{2} pair production processes is discussed via gluongluon fusion for a collision energy of 13 TeV, and the cases in which one Higgs boson decays to $\gamma\gamma$γγ and the other one decays to $b\bar{b}$bbˉare considered. Thus, the emphasis is on the resonant process: gg→H_{3} →H_{1}H_{2} → $\gamma\gamma b\bar{b}$γγbbˉwhich is a Next to Leading Order process (NLO). This process has been simulated using MadGraph5 v2.6.5. MadGraph5 is a Monte Carlo (MC) event generator aimed at simulating events for SM and BSM phenomenology. The result of the simulation is stored in a format, called an LHE (Les Houches Event) file. Anaylsing the data includes studying the kinematics ( $p_{T}, \eta, \phi$pT,η,ϕ) of different particles. Various plots for studying the kinematics of the resonance particle H_{3} and all the final state particles are made and the possibility of resonant double Higgs (H_{1}H_{2} in this case) production using the NMSSM model has been studied.
Abbreviations
NLO  Next to Leading Order 
BEH  BroutEnglertHiggs 
SUSY  Supersymmetry 
NMSSM  NexttoMinimalistic Supersymmetric Standard Model 
LHE  Les Houches Event 
BR  Branching Ratio 
INTRODUCTION
Background
The standard model^{[1]}^{ }
The Standard Model (SM) of high energy physics has been the most successful theory so far, capable of explaining three of the four known fundamental forces of nature – the strong force, the weak force and the electromagnetic force and also in classifying all the rudimentary particles of the universe. Developed in about the latter half of the 20^{th} Century, the standard Model of particle physics was quite successful in predicting the existence of certain particles like the vector bosons – Z, W ^{±}, the gluon, the top and the charm quarks and the Higgs boson. The Standard Model is both strikingly simple and robust. Almost every quantity that has been measured in particle physics experiments over the past five decades falls right on the predicted value, within experimental error limits.
Of the many predictions by the standard model, the discovery of the SM Higgs boson with a mass of about 125 GeV, in 2012, from data collected with CMS and ATLAS detector at CERN, Switzerland, has proven the Standard Model to be the most successful theory of the particle physics. The Higgs boson, predicted as a consequence of the spontaneous symmetry breaking in the electroweak sector of the SM, explains how SM particles get masses, where Higgs boson couples to itself to become massive. Nevertheless, the standard Model does have certain limitations.
One major drawback of the Standard Model is that it does not include the gravitational force, one of the four fundamental forces. Also, it does not address the Hierachy problem. The Hierarchy problem refers to the huge difference in the strengths of the fundamental forces. It also refers to the wide fluctutations in mass for the elementary particles. We see the electron is about 200 times lighter than the muon and 3500 times lighter than the tau. Same thing for the quarks: the top quark is 75,000 times heavier than the up quark. Why is there such a wide spectrum of masses among the building blocks of matter? The hierarchy problem is also related to the mass of the Higgs boson. The equations of the Standard Model establish relations between the fundamental particles. For example, in the equations, the Higgs boson has a basic mass to which theorists add a correction for each particle that interact with it. The heavier the particle, the larger the correction. The top quark being the heaviest particle, it adds such a large correction to the theoretical Higgs boson mass that theorists wonder how the measured Higgs boson mass can be as small as it was found.
In the context of this project, we are studying the double Higgs production. The crosssection of the SM double Higgs (hh) production process is quite low and hence, difficult to detect. We therefore look at Beyond the Standard Model (BSM) double Higgs production to understand the selfcoupling which predicts an increased crosssection. One of the most popular Supersymmetry (SUSY) models – the NexttoMinimal Supersymmetric Standard Model (NMSSM).
Double Higgs production^{[2]}^{ }^{ }
The BroutEnglertHiggs (BEH) mechanism is at the heart of the Standard Model. It introduces the Higgs field which interacts with other particles and gives them mass. The BEH mechanism also prognosticate that Higgs field can interact with itself; in other words, a single (virtual) Higgs boson can decay into two Higgs bosons. Observing and measuring this selfinteraction, or 'Higgs selfcoupling', would be the final validation of the theory of mass generation, while any inconsistency with Standard Model predictions would hint at new physics.
Unfortunately, Higgs boson pairs are predicted to be very rare in proton–proton collisions, with a production rate roughly a thousand times smaller than the single Higgs boson. To make matters worse, not all double Higgs boson production occurs through Higgs selfcoupling. The Higgs boson pairs have various decay channels. The most sensitive of these involve one Higgs boson decaying into a pair of 'b' quarks and one decaying into either another pair of 'b' quarks (hh→bbbb), two tau (τ) leptons (h→bbττ) or two photons (hh→bbγγ). These three searches when statistically combined, the production rate of HH pairs could be excluded beyond 6.7 times the Standard Model prediction, at a 95% confidence level. New physics could be indicated by a Higgs selfcoupling which differs from the Standard Model prediction. This would affect the production rate and the kinematic distributions of the Higgs boson pairs and, as such, is an excellent probe for new physics.
Higgs boson pairs are also a key signature of heavy new particle decays in several scenarios beyond the Standard Model. These might include additional Higgs bosons in models that extend the Higgs sector of the Standard Model, or the excitation of a graviton in models with extra spatial dimensions.
Nexttominimal supersymmetric standard model (NMSSM)^{[3]}^{ }
The NexttoMinimal Supersymmetric Standard Model is a supersymmetric extension of the Standard Model, motivated by a solution to the hierarchy problem. The NMSSM predicts the existence of seven Higgs bosons. Typically, the Higgs bosons in the NMSSM are denoted in an order with increasing masses  H_{1}, H_{2}, ...., H_{7}. In the special case of a charge parity (CP)  conserving Higgs potential, we have three CP even Higgs bosons  H_{1}, H_{2}, H_{3}, two CP odd ones  A_{1}, A_{2}, and a pair of charged Higgs bosons  H^{+}, H^{ }. The A_{1 }, A_{2} are pseudoscalars, meaning these are particles with spin 0 and odd parity, that is, a particle with no intrinsic spin with wave function that changes sign under parity inversion.
In this project, it is assumed that H_{1} is the SMlike Higgs boson. We discuss the H_{1}H_{2} pair production processes via gluongluon fusion for a collision energy of 13 TeV, and consider the cases in which one Higgs boson decays to $\gamma\gamma$and the other one decays to $b\bar{b}$. Thus, the emphasis is on the resonant process: gg→H_{3} →H_{1}H_{2} → $\gamma\gamma b\bar{b}$which is a Next to Leading Order process (NLO). In this process H_{1} decays to $\gamma\gamma$and H_{2} decays to $b\bar{b}$. We also discuss another process where H_{1} decays to $b\bar{b}$and H_{2} decays to $\gamma\gamma$, i.e, gg→H_{3} →H_{1}H_{2} → $b\bar{b} \gamma\gamma$.
Objective of the Research
The objective of this project is to study the resonant double Higgs production (hh) process in the NMSSM, which is important for understanding the Higgs self coupling and its measurement and also to determine if the analysis for the SM double Higgs production can be applied to this model. The double Higgs production process hasn't been observed yet experimentally, at any one of the experimental facilties. Thus, this project is an attempt to simulate the double higgs production process using a Monte Carlo event generator and analyse the produced data using a data analysis tool.
If the double Higgs coupling related to beyond the standard Model scearios (NMSSM in this case) is observed, then it could open a window to new physics, which we haven't been able to observe yet.
Scope
In this project, we are trying to generate the processes: gg→H_{3} →H_{1}H_{2} → $\gamma\gamma b\bar{b}$ and gg→H_{3} →H_{1}H_{2} → $b\bar{b} \gamma\gamma$ in the NMSSM Model. Throughout this project, all the processes have been simulated with a center of mass energy, s=13 TeV which was the center of mass energy for the protonproton collisions of the Large Hadron Collider for Run 2. About 10,000 events have been produced for all of the processes using a Monte Carlo event generator  MadGraph5. While generating either of the processes, we keep the mass of SM higgs fixed, H_{1} = 125 GeV, and vary the masses of H_{2 }as [60, 100, 150, 200, 300, 400, 600, 800] GeV and H_{3 } as [800, 1000, 1500] GeV. Also, we set the decay widths of H_{1}, H_{2} and H_{3} to be 4 MeV, 1 MeV and 1 MeV respectively. All the events generated are at the generator level and hence, no cuts have been applied to take into account the detector specifications. LHAPDF6 was used as the parton distribution function (pdf) in this process.
Any deviations of the selfcoupling from standard model predictions, i.e, validation of BSM scenarios (NMSSM in this case) would open a window on new physics which we are not able to observe yet experimentally.
SOME BASIC CONCEPTS AND DEFINITIONS
Feynman Diagrams^{[4]}^{ }
Feynman diagrams are used for representation of various kinds of particle interactions. They were invented by Nobel Laureate Richard Feynman, to make calculation of scattering amplitudes, cross sections and decay rates, faster and easier. It is a twodimensional representation in which one axis represents the time axis and other axis represents the space axis. An example of a feynman diagram is given below.
Some salient features of a Feynman diagram are:
 Straight lines are used to represent fermions.
 Wavy lines are used to represent bosons.
 The particles moving backward in time are Antiparticles.
 Each process can be represented by mutliple Feynman diagrams. Thus, to calculate the probability of a process to occur all feynman diagrams must be taken into consideration.
 At each vertex, electric charge must be conserved and except in weak interactions, the quark or the lepton flavor must be conserved.
 A Feynman diagram can be used to represent Annihilation/formation or scattering of particles. Fig 3 is an example of annihilation/formation Feynman diagram. Particle A and B collide to form particle X which later decays to C and D.
 An example of a scattering/exchange diagram is given below.
Here particle A scatters off particle B by exchanging particle X. Particle A becomes C and particle B becomes D. We don't know if A emitted X and B absorbed it or vice versa.
Particle Decays
The following properties of particle decay can be measured experimentally:
 The decay rate ' $\Gamma$' of a particle in its own rest frame is defined as the probability per unit time the particle will decay:$dN = \Gamma Ndt \rightarrow N(t) = N_{0} e^{\Gamma t}$.
 Most particles decay through more than one different route. Add up all decay rates to obtain the total decay rate: $\Gamma_{tot} = \Sigma^{n}_{i=1} \Gamma_{i}$.
 Decay width : It is a measure of the probability of a specific decay process occuring within a given amount of time in the parent particle’s rest frame. Since the dimension of ' $\Gamma$ ' is the inverse of time, in our system of natural units, it has the same dimension as mass (or energy). When the mass of an elementary particle is measured, the total rate shows up as the irreducible “width” of the shape of the distribution. Hence the name decay 'width'.
 The decay rate ' $\Gamma$' of a particle in its own rest frame is defined as the probability per unit time the particle will decay:$dN = \Gamma Ndt \rightarrow N(t) = N_{0} e^{\Gamma t}$.
Particle Collisions
 For particle collisions we measure the cross section (σ) of the process. The cross section is a measure of how often a process happens per unit of incident flux. It has dimensions of area, and is measured often measured in barn, b. 1 b ≡ 10^{−28} m^{2}. At the Large Hadron Collider, typical cross sections range from 0.1 b for the inclusive cross section to 0.1 pb for Higgs boson production.
 An elastic collision (or elastic scattering) occurs when the initial and final state particles are the same. Kinetic energy and threemomentum is conserved.
 An inelastic collision (or inelastic scattering) occurs when the initial and final state particles are different. Kinetic energy and threemomentum will not be conserved. (However fourmomentum will be conserved!)
 The exclusive cross section is one with a given final state, e.g. at the LHC, we are interested in measuring how often Higgs bosons are created along with W bosons, so we would like to measure σ(pp → W H). Exclusive cross sections are easier to calculate.
 Cross sections are measured experimentally by simply counting the number of times that a process occurs:
$N(pp\rightarrow WH) = \sigma(pp\rightarrow WH) \times \int\mathcal{L} dt$
where L is the luminosity of the colliding beams, and $\int \mathcal{L} dt$ is the integrated luminosity for the period of time the experiment has run for.
 Luminosity, also known as instantaneous luminosity, can be (roughly) defined as:
$\mathcal{L} = \frac{N_{a} \times N_{b} \times collision frequency}{Overlap area}$
Where N_{a} and N_{b} are the numbers of particle in the overlap area in colliding bunches a and b. Luminosity is measured in dimensions of inverse area per time unit. At the LHC, luminosity is often quoted in units 10^{32} cm^{−2} s^{−1}, equivalent to 1 inverse barn per second. Integrated luminosity is measured in units of b^{−1}. At the LHC, is it often quoted in units of inverse femtobarns, fb^{−1}.
Coupling
Particles which interact with each other are said to be coupled. This interaction is caused by one of the fundamental forces, whose strengths are usually given by a dimensionless coupling constant. In quantum electrodynamics, this value is known as the finestructure constant α, approximately equal to 1/137. For quantum chromodynamics, the constant changes with respect to the distance between the particles. This phenomenon is known as asymptotic freedom. Forces which have a coupling constant greater than 1 are said to be "strongly coupled" while those with constants less than 1 are said to be "weakly coupled".
Higgs Decay
Though Higgs boson has a mass of 125 GeV/c^{2} and a mean lifetime of about 1.6×10^{22} s, it decays almost instantly. As it interacts with all the massive elementary particles of the Standard Model, the Higgs boson has many different channels through which it can decay. Each of these possible channels has its own branching ratio. Here are some of the decay modes:
 Higgs decaying to photons ( $H\rightarrow \gamma\gamma$)
Only a small fraction of the Higgs boson produced decay into 2 photons. The branching fraction for this mode is approximately 2.28 x 10^{22} for a Higgs of mass 125 GeV/c^{2}. This decay is rare because the Higgs cannot interact directly with the two photons as the photon is massless. However, the decay can still happen because of an exchange of intermediate particles
 Higgs decaying to tau leptons ( $H\rightarrow \tau^{+} \tau^{}$)
𝜏 leptons decay before being detected. They can decay in 3 ways: through e^{} + neutrinos, through 𝜇^{} + neutrinos and Hadronically. In hadronic decay, there can be one or three charged hadrons which will be reconstructed as tracks in the detector and there might also be some neutral pions which are reconstructed as energy deposits in the calorimeters. However, the probability of H→𝜏^{+}𝜏^{} is much larger than that of the decay to photons or Z bosons. This mode helps us to understand the mechanism which allows the Higgs boson to give mass to the fermions.
 Higgs decaying to Z bosons ( $H \rightarrow Z Z$)
The two Z bosons decay to various particles such as e^{} e^{+}, 𝜇^{} 𝜇^{+}. This decay mode happens once in about 10,000 decay of Higgs. It is a clean channel because there are not many backgrounds. Looking at the angles between the directions of the electrons and muons, it was found that the Higgs boson has no spin.
 Higgs decaying to bottom quarks ( $H\rightarrow b\bar{b}$)
Higgs decays to bottomquarks about 58% of the time (Branching ratio is about 0.577 for Higgs mass =125 GeV/c^{2}). Yet, it is difficult to observe this because of the large background.
 Decay to W bosons ( $H\rightarrow W^{+}W^{}$)
W bosons are unstable particles which decay into hadrons or leptons. For example, they may decay into e^{} or $\mu$and their corresponding neutrinos.
Kinematic Variables ( $p_{T}, \eta, \phi, m_{o}$)^{ }^{[6]}^{ }
 Transverse momentum ( $p_{T}$): It is the component of the momentum perpendicular to the beam line. It's importance arises because momentum along the beamline may just be left over from the beam particles, while the transverse momentum is always associated with whatever physics happened at the vertex. It can be calculated as $p_{T} = \sqrt{(p_{x}^{2} + p_{y}^{2})}$ where $p_{x} and p_{y}$ are momentum along the axis
 Pseudorapidity ( $\eta$): It is a commonly used spatial coordinate for describing the angle of a particle relative to the beam axis. It is defined as $\eta =  \ln[tan(\frac{\theta}{2})]$ , where $\theta$ is the angle between the particle threemomentum p and the positive direction of the beam axis. It is used instead of polar angle (θ) in hadron collider experiments because the center of mass frame of the interacting particles is boosted in the lab at hadron colliders.
 Azimuthal angle ( $\phi$): It is given by $\phi = tan^{1} \frac{p_{y}}{p_{x}}$.
 Invariant mass (m_{o}): In high energy physics, the invariant mass is equal to the mass in the rest frame of the particle, and can be calculated by the particle's energy 'E' and its momentum 'p' as measured in any frame, by the energymomentum relation: $m_{o}^{2} c^{2} = (\frac{E}{c})^{2}  p^{2}$. In natural units, when c=1, $m_{o}^{2} = {E}^{2}  p^{2}$.
METHODOLOGY
Event Generation Using MadGraph5
MadGraph5_aMC@NLO
MadGraph5_aMC@NLO is a Monte Carlo event generator that aims at providing all the elements necessary for SM and BSM phenomenology, such as the calculations of cross sections, event generation and the use of a variety of tools relevant to event manipulation and analysis. For this project, we will be using MadGraph5 v2.6.5. The result of the simulation is stored in a format, called an LHE (Les Houches Event) file.
Event generation
For this project, our aim is to generate the processes $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow\gamma\gamma b\bar{b}$and $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow b\bar{b} \gamma\gamma$ for 10,000 events at a center of mass energy of protonproton collisions, s = 13 TeV. Note that all the events generated are at the generator level and hence, no cuts have been applied to take into account the detector specifications.
 First, we generate the process: $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow\gamma\gamma b\bar{b}$ . In order to generate the processes, we need to add our model NMSSM to the 'model' folder in MadGraph5 directory. After that we can execute the following set of commands:
> set group_subprocesses Auto
> set ignore_six_quark_processes False
> set loop_optimized_output True
> set gauge unitary
> set complex_mass_scheme False
> import model NMSSMHET_UFO modelname
> generate g g > h03 , (h03 > h2 h01, h2 > b b~, h01 > a a)
>output gg_h3_h1h2_aabb~
>launch
>0
 After entering '0', the MadGraph window allows you to edit the parameter and run cards. For editing the parameter card, enter
>1
 This will open the param_card.dat file in the 'vim editor'. To insert the text, press the 'i' key. Change the following parameters in the param_card file:
> set param_card mass 25 125
> set param_card mass 35 scan:[60,100,150,200,300,400,600,800]
> set param_card mass 45 scan:[800,1000,1500]
> set param_card DECAY 25 4
> set param_card DECAY 35 1
> set param_card DECAY 45 1
 After changing the parmeters, to save and quite the param_card file, press 'shift+z' twice. This will save the changes in the dat file. For editing the run card, enter
> 2
 This will open the run_card.dat file in the 'vim editor'. Make the following changes:
lhapdf = pdlabel ! PDF set
292600= lhaid ! if pdlabel=lhapdf
 After changing the parmeters, to save and quite the param_card file, press 'shift+z' twice. This will save changes in the dat file. Now after setting all the parameters, Enter
> 0
 This will produce 23 LHE files corresponding to the variations in the masses of H_{2} and H_{3} .
 Repeat the above steps for generating the process: $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow\gamma\gamma b\bar{b}$.
Feynman Diagrams of the processes generated using MadGraph5:
Even Data Analysis using ROOT
ROOT  a Data analysis Framework
Root is a standard scientific software toolkit. It has all the components needed to deal with big data processing, statistical analysis, visualisation and storage. It is primarily written in C++ but integrated with other languages such as Python and R. For this project, the version ROOT v6.14/06 has been used.
Data plotting and analysis
 After the generation of LHE files, the next step is to plot the histograms of kinematic variables of the resonance particle H_{3 }and all the final state particles. We do this using PyROOT which is a Python extension module. It allows us to use any ROOT class from the python interpreter.
 To plot the histograms, the following steps have been used:
 Read the LHE files and scan the four momentum of the particles using TLorentzvector class of ROOT.
 Store the value of the four momentum in a variable.
 Using '.Pt()', '.Eta()', 'Phi()' and '.M()' function of ROOT, retrieve the values of the Transverse momentum ( $p_{T}$), pseudorapidity ( $\eta$), Azimuthal angle ( $\phi$) and Masses of the particles.
 Plot the Normalized histograms for the kinematic variables of all particles produced during the process.
 Use proper bin size, Yaxis scale for plotting for each of the histograms.
 For reconstruction of mass of H_{1 , }combine the four momentum of both the photons (for the process where H_{1} decays to $\gamma\gamma$) and use the '.M()' function.
 Similarly, reconstruct the masses of H_{2 }and H_{3 .}
 Save the histograms in '.png' and '.root' formats.
RESULTS AND DISCUSSIONS
Plot of CrossSections
 Plot of Crosssections for the process: $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow\gamma\gamma b\bar{b}$
Through this plot, we can observe that the crosssection is highest for the process, where mass of H_{3 }= 1500 GeV. There is no trend observed in the plot.
 Plot of Crosssections for the process: $gg\rightarrow H_3\rightarrow H_1H_2\rightarrow b\bar{b} \gamma\gamma$
Through this plot, we can observe that the crosssection is the highest for the mass of H_{3} = 1800 GeV. With the increase in mass of H_{3 }, there is an increase in the cross sections of the process. _{ }
Plots of Kinematic Variables _{ }
We plot the histogram of kinematics variables: transverse momentum( $p_T$), pseudorapidity ( $\eta$) and azimuthal angle ( $\phi$) for all final state particles ( $b, \bar{b} and \gamma, \gamma$) have been made. $p_T$, $\eta$, $\phi$ and invariant mass plots have been made for H_{1} and H_{2}, and invariant mass plot for resonance H_{3}. The plots have been made in such a way that the mass of H_{3} is kept constant and the mass of H_{2} is varied.
Plots for $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow\gamma\gamma b\bar{b}$
 Mass of H_{3} = 800 GeV
We see a resonance peak at M_{γγbb~ }= 800 GeV for all the distributions of Fig 13 , which corresponds to the mass of H_{3 .}
We see a resonance peak at M_{γγ }= 800 GeV for all the distributions of Fig 14, which corresponds to the mass of H_{1 }.
The distribution of pT of γ_{1} is slightly skewed (see Sec. APPENDICES) right. Also, many γ_{1} particles have momentum about 40 GeV. The maximum value of pT is about 400 GeV. Also, the right skewness of the distributions decreases with the increasing mass of H_{2} . But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }
The distribution of pT of γ_{2} is slightly skewed right. Also, many γ_{2} particles have momentum about 30 GeV. The maximum value of pT is about 400 GeV. Also, the right skewness of the distributions decreases with the increasing mass of H_{2} . But the maximum value of p_{T} is decreasing with the increasing mass of H_{2.}
The distributions of pT for b are slightly skewed right for M_{H2} = 60, 100 GeV. The skewness of the distributions for M_{H2} = 300, 400 GeV, increases towards left. The distribution for M_{H2} = 600 GeV peaks almost at 240 GeV. The maximum value of pT is 400 GeV.
The distributions of pT for b~ are slightly skewed right for M_{H2} = 60, 100 GeV. The skewness of the distributions for M_{H2} = 300, 400 GeV, increases towards left. The distribution for M_{H2} = 600 GeV peaks almost at 240 GeV. The maximum value of pT is 400 GeV.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
The azimuthal angle (ϕ) distribution is almost symmetric about the origin. All the particles are equally distributed among all the values of ϕ from π to +π. Since, ϕ distributions for all the particles are almost same, we will not analyse the ϕ plots henceforth.
 Mass of H_{3} = 1000 GeV
We see a resonance peak at M_{γγbb~ }= 1000 GeV for all the distributions of Fig 24 , which corresponds to the mass of H_{3 .}
We see a resonance peak at M_{γγ }= 125 GeV for all the distributions of Fig 25 , which corresponds to the mass of H_{1 .}
The distribution of pT of γ_{1} is slightly skewed right. Also, many γ_{1} particles have momentum about 40 GeV. The maximum value of pT is about 500 GeV. Also, the right skewness of the distributions decrease with the increasing mass of H_{2} . But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }
The distribution of pT of γ_{2} is slightly skewed right. Also, many γ_{2} particles have momentum in the range of 3040 GeV. The maximum value of pT is about 500 GeV. Also, the right skewness of the distributions decrease with the increasing mass of H_{2} . But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }
The distributions of pT for b are slightly skewed right for M_{H2} = 60, 100 GeV. The skewness of the distributions for M_{H2} = 400, 600 GeV, increases towards left. The distribution for M_{H2} = 800 GeV peaks almost at 360 GeV. The maximum value of pT is 500 GeV.
The distributions of pT for b~ are slightly skewed right for M_{H2} = 60, 100 GeV. The skewness of the distributions for M_{H2} = 400, 600 GeV, increases towards left. The distribution for M_{H2} = 800 GeV peaks almost at 350 GeV. The maximum value of pT is 500 GeV.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
The azimuthal angle (ϕ) distribution is almost symmetric about the origin. All the particles are equally distributed among all the values of ϕ from π to +π. Since, ϕ distributions for all the particles are almost same, we will not analyse the ϕ plots henceforth.
 Mass of H_{3} = 1500 GeV
We see a resonance peak at M_{γγbb~ }= 1500 GeV for all the distributions of Fig 35 , which corresponds to the mass of H_{3 .}
We see a resonance peak at M_{γγ }= 125 GeV for all the distributions of Fig 36 , which corresponds to the mass of H_{1 .}
The distribution of pT of γ_{1} is slightly skewed right. Also, many γ_{1} particles have momentum about 40 GeV. The maximum value of pT is about 760 GeV. Also, the right skewness of the distributions decrease with the increasing mass of H_{2} . But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }
The distribution of pT of γ_{2} is slightly skewed right. Also, many γ_{1} particles have momentum about 40 GeV. The maximum value of pT is about 760 GeV. Also, the right skewness of the distributions decrease with the increasing mass of H_{2} . But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }
The distributions of pT for b are slightly skewed right for M_{H2} = 60, 100 GeV. The skewness of the distributions for M_{H2} = 600, 800 GeV, increases towards left. The distribution for M_{H2} = 800 GeV peaks almost at 240 GeV. The maximum value of pT is 760 GeV.
The distributions of pT for b~ are slightly skewed right for M_{H2} = 60, 100 GeV. The skewness of the distributions for M_{H2} = 600, 800 GeV, increases towards left. The distribution for M_{H2} = 800 GeV peaks almost at 240 GeV. The maximum value of pT is 740 GeV.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
The azimuthal angle (ϕ) distribution is almost symmetric about the origin. All the particles are equally distributed among all the values of ϕ from π to +π. Since, ϕ distributions for all the particles are almost same, we will not analyse the ϕ plots henceforth.
Plots for $gg\rightarrow H_{3}\rightarrow H_{1} H_{2}\rightarrow b\bar{b} \gamma\gamma$
 Mass of H_{3} = 800 GeV
We see a resonance peak at M_{γγbb~ }= 800 GeV for all the distributions of Fig 46 , which corresponds to the mass of H_{3 .}
We see a resonance peak at M_{bb~ }= 125 GeV for all the distributions of Fig 47 , which corresponds to the mass of H_{1 .}
The distribution of pT of γ_{1} is slightly skewed right. The maximum value of pT is about 400 GeV. Also, the right skewness of the distributions decrease with the increasing mass of H_{2}.
The distribution of pT of γ_{2} is slightly skewed right. The maximum value of pT is about 400 GeV. Also, the right skewness of the distributions decrease with the increasing mass of H_{2}. The distribution for which the mass of H_{2} = 600 GeV, peaks at 240 GeV.
The distribution of pT of b is slightly skewed right. The maximum value of pT is about 405 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}. But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }_{ }
The distribution of pT of b~ is slightly skewed right. The maximum value of pT is about 400 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}. But the maximum value of p_{T} is decreasing with the increasing mass of H_{2. }
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
The azimuthal angle (ϕ) distribution is almost symmetric about the origin. All the particles are equally distributed among all the values of ϕ from π to +π. Since, ϕ distributions for all the particles are almost same, we will not analyse the ϕ plots henceforth.
 Mass of H_{3} = 1000 GeV
We see a resonance peak at M_{γγbb~ }= 1000 GeV for all the distributions of Fig 57 , which corresponds to the mass of H_{3 }.
We see a resonance peak at M_{bb~ }= 125 GeV for all the distributions of Fig 58, which corresponds to the mass of H_{1 .}
The distribution of pT of γ_{1} is slightly skewed right. The maximum value of pT is about 500 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}.
The distribution of pT of γ_{2} is slightly skewed right. The maximum value of pT is about 500 GeV. Also, the right skewness of the distributions decreases with the increasing mass of H_{2}.
The distribution of pT of b is slightly skewed right. The maximum value of pT is about 500 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}.
The distribution of pT of b~ is slightly skewed right. The maximum value of pT is about 500 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
The azimuthal angle (ϕ) distribution is almost symmetric about the origin. All the particles are equally distributed among all the values of ϕ from π to +π. Since, ϕ distributions for all the particles are almost same, we will not analyse the ϕ plots henceforth.
 Mass of H_{3} = 1500 GeV
We see a resonance peak at M_{γγbb~ }= 1500 GeV for all the distributions of Fig 68 , which corresponds to the mass of H_{3 }.
We see a resonance peak at M_{bb~ }= 125 GeV for all the distributions of Fig 69 , which corresponds to the mass of H_{3 }.
The distribution of pT of γ_{1} is slightly skewed right. The maximum value of pT is about 770 GeV. Also, the right skewness of the distributions decreases with the increasing mass of H_{2}.
The distribution of pT of γ_{2} is slightly skewed right. The maximum value of pT is about 770 GeV. Also, the right skewness of the distributions decreases with the increasing mass of H_{2}.
The distribution of pT of b is slightly skewed right. The maximum value of pT is about 770 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}.
The distribution of pT of b~ is slightly skewed right. The maximum value of pT is about 790 GeV. Also, the right skewness of the distributions increases with the increasing mass of H_{2}.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
From the above distributions. we can observe that as p_{T }increases $\eta$tends to 0. Also, the distributions are quite symmetric about 0.
The azimuthal angle (ϕ) distribution is almost symmetric about the origin. All the particles are equally distributed among all the values of ϕ from π to +π. Since, ϕ distributions for all the particles are almost same, we will not analyse the ϕ plots henceforth.
CONCLUSION AND RECOMMENDATIONS
 Through this project, I developed a good understanding of some of the important concepts of particle physics at the elementary level  Feyman diagrams, Relativistic Kinematics, the Standard Model and few ideas about the physics beyond the standard Model. It also helped me to elevate to a level of proficiency in event simulation and data analysis.
 This project aimed at giving a glimpse of how event generation takes place at the generator level, and how the produced event data is analyzed using research grade software tools. As the project included production of large amount of data, it was succesful at giving a basic idea of how the real event data is tackled in actual research.
 The project was aimed at studying the Resonant double higgs production in the NMSSM using a Monte Carlo even generator which has not been observed yet experimentally. We studied the processes: gg→H_{3} →H_{1}H_{2} → $\gamma\gamma b\bar{b}$ and gg→H_{3} →H_{1}H_{2} → $b\bar{b} \gamma\gamma$for the masses of H_{2}=[60, 100, 150, 200, 300, 400, 600, 800] GeV and and H_{3} =[800, 1000, 1500] GeV. We found that for the masses of H_{2 }= [400, 600] GeV and H_{3} = 800 GeV, the final state particles were very less boosted and thus, we feel that the analysis for SM double higgs process can be applied to this.
 This project will be heplful for further studies on beyond the standard model double higgs production, which has not been observed experimentally yet.
 If the NMSSM double higgs production is observed experimentally, then it would open a window for new physics, might overcome the limitations of the Standard Model and can give answers to some of the unsolved problems. Thus, the emphasis would be on experimental testing of the theory of NMSSM at high luminosity.
ACKNOWLEDGEMENTS
Firstly, I would like to express my sincere gratitude to my guide Dr. Jyothsna Rani Komaragiri for giving me this opportunity to pursue a summer internship under her guidance. Also, I would like to thank the Indian Academy of Sciences for facilitating my summer research internship at IISc, Bengaluru.
I would like to give many thanks to the Experimental High Energy physics PhD group of Lata Panwar, Deepak Kumar and Praveen Tiwari at Indian Institute of Science for their help and support. I would also like to thank my fellow summer interns Shreyas Shukla and Shilpi Jain for the efforts they put into making this programme an eventful and enjoyable experience. Last but not least, I would like to thank my room mates Akella Sravan, Lakshay Jain, Nilendu Das, Vasu Garg and Hariharan for making my stay at IAS Fellows residency, a pleasant experience.
APPENDICES
 Skewness: Skewness, in statistics, is the degree of distortion from the symmetrical bell curve, or normal distribution, in a set of data. Skewness can be negative, positive, zero or undefined. A normal distribution has a skew of zero, while a lognormal distribution, for example, would exhibit some degree of rightskew.
In the above figure is a lognormal distribution, which is rightly skewed.
 Kurtosis: It is a measure of whether the data are heavytailed or lighttailed relative to a normal distribution.
References

https://www.quantumdiaries.org/2014/03/14/thestandardmodelabeautifulbutflawedtheory/

https://atlas.cern/updates/physicsbriefing/atlassearchesdoublehiggsproduction

Ellwanger, Ulrich and Hugonie, Cyril and Teixeira, Ana M. (2010). The NexttoMinimal Supersymmetric Standard Model. 496,

https://www.britannica.com/science/Feynmandiagram

https://arxiv.org/pdf/1604.02651.pdf
Source

Fig 1: https://en.wikipedia.org/wiki/Standard_Model#/media/File:Standard_Model_of_Elementary_Particles.svg

Fig 2: https://cds.cern.ch/record/1638469/files/higgspotential.png

Fig 4: https://feynman.aivazis.com/

Fig 6: https://qph.fs.quoracdn.net/mainqimgb0248c6638b23d2896bcd37ebd2fe68f.webp

Fig 7: https://www.researchgate.net/profile/Oliver_Passon/publication/329862752/figure/fig3/AS:706516638261248@1545458101056/LeftRepresentativeFeynmandiagramsfortheHiggsproductionviagluongluonfusionand_W640.jpg

Fig 8: https://cerncourier.com/objects/ccr/cern/55/5/3/CCnew4_05_15.jpg

Fig 9: https://atlas.physicsmasterclasses.org/wpath_files/img/Feynman/Higgs.png
Post your comments
Please try again.