Monte Carlo Threshold Energy Estimation for A + B → C + D processes : An Educational Resource in Experimental High Energy Physics Research

Specialized documentation envisioned from a pedagogical bases to train scientifically and technologically teachers and researchers, who initiate themselves in the analysis of high energy physics (HEP) experiments, is scarce. The lack of this material makes that young scientists’ learning process be prolonged in time, raising costs in experimental research. In this paper we present the Monte Carlo technique applied to simulate the threshold energy for producing final-state particles of a specific two-body process (A + B → C + D), as pedagogical environment to face both computationally and conceptually an experimental analysis. The active/interactive learning-teaching formative process presented here is expected to be an educational resource for reducing young scientists’ learning curve and saving time and costs in HEP scientific research.


Introduction
This article is the result of a research study titled "Design of an educational resource aimed at the study of relativistic kinematics of a two-by-two process" carried out at the Universidad Distrital Francisco José de Caldas, whose goal was the generation of an educational resource based on ROOT-C++ for the study of a 2 → 2 scattering processes, in particular, with the understanding of the experiment details at the level of both physics and instrumentation; the second one corresponds to a description of the kinematic variables needed for event reconstruction; the third one is linked with computational methods used to correct experimental observables determined through the analysis of experimental data, namely, the Monte Carlo technique; and the last one, related to the implementation of the analysis codes, makes reference to the validation of the software analysis in terms of already-known results.
In accordance with the above and for the purpose of illustrating the Monte Carlo technique applied to a A + B → C + D process, the pedagogical route for an active/interactive learning process of this educational resource is based on several aspects: discussion of the main physics aspects that justify the choice of the γ + p → K + + Λ 0 reaction, the detector settings required to perform the experiment, the description of the reaction in terms of physical variables (energy, momentum, Mandelstam variables), the analytical calculation of the minimum energy (threshold energy) for production of particles in the final state, the method of Monte Carlo along with its computational implementation to estimate the threshold energy, and lastly, the corroboration of the results.

Experimental setup
In this work, the reaction γ + p → K + + Λ 0 was chosen as the reference process.From the experimental point of view, this reaction has several elements to consider.
The production of a photon beam (γ) is well known in the field of accelerators around the world and the modeling of a hydrogen target as a cluster of protons (p) results very familiar in physics.
The kaon K + is a charged particle that can be detected relatively easily due precisely to its electric charge and significant stability (at | p| K + =1 GeV/c, K + decays, on average, after 7.5 m); lambda, Λ 0 , on the contrary, is a neutral and unstable particle (at | p| Λ 0 =1 GeV/c, Λ 0 decays, on average, after 0.95 m) whose direct detection is unfeasible, making it more convenient to reconstruct through the detection of charged decay products (mainly protons and negative pions, π − ) or through its neutral decay mode (neutrons and neutral pions); in both cases, the reconstruction imposes conditions of energy and momentum conservation.
In addition, the γ + p → K + + Λ 0 has been measured in various facilities around the world for which real data are currently available.In particular, we have access to reaction's real data taken at the Hall B of the Thomas Jefferson National Accelerator Laboratory (TJNAF), located at Newport News, VA, USA [1].
An overall description of the experiment begins in the CEBAF linear accelerator (Continuous Electron Beam Accelerator Facility) where a continuous electron beam with a maximum energy of 6 GeV is produced (a 12 GeV upgrade is in progress).The resulting electron beam is driven to the experimental Hall where it interacts with the electromagnetic field of an atomic nucleus with charge Z (radiator, i.e a gold foil).As a consequence of the interaction, an incoming electron, initially with energy E 0 , is accelerated by the field and emerges with lower energy E e (energy-degraded electron), emitting a photon of energy E γ .This process is kinematically possible only if a small amount of momentum is transferred to the nucleus (recoil momentum).
In the relativistic limit, bremsstrahlung is the dominant mode of energy loss for electrons in a material, and the nucleus recoil energy can usually be neglected [2][3][4].Energy and momentum conservation then gives that E 0 = E e + E γ and P 0 = P e + P γ , with P 0 , P e , and P γ , being the momentum of the incident electron, emerging or outgoing electron, and emitted bremsstrahlung photon, respectively.
The spectrum of bremsstrahlung photons is not monoenergetic; however, a beam of incident electrons with a fixed energy produces a beam of photons spread over a range of energies.The polychromatic nature of these photons coming out from this bremsstrahlung process, along with the fact that the electron and photon beams come out from the radiator (in this case a thin-gold foil) nearly parallel (a mixed beam) with respect to the incident electron beam direction requires the use of a tagger system to find the energy of the photons interacting in the target.
The outgoing beams are therefore separated when passing through the tagger-system magnet; the electrons are bent downwards by the magnet while the photons continue towards the target.The tagger magnetic field is properly matched to the incident electron beam energy (E 0 ) such that electrons which do not radiate (full-energy electrons) are directed into a beam dump.Electrons which do generate a bremsstrahlung photon are deflected significantly by the tagger magnet, with a radius of curvature depending on the fraction of incident energy transferred to the photon (E γ /E 0 ).This fractional energy allows the photon-tagging system to tag photons with energies between 20% and 95% of the incident electron energy.The photon energy E γ is then calculated by energy conservation (neglecting the recoil energy The photon beam, typically with an energy range between 1.0 and 2.5 GeV, then faces a collimator1 whose purpose is to reduce the angular spread of the beam, to finally impinge directly on the hydrogen target placed inside the CLAS detector as shown in Fig. 1.CLAS, CEBAF Large Acceptance Spectrometer, corresponds to a spectrometer with a coverage nearly to a 4π in solid angle, which allowed to detect charge and some neutral particles [5,6].This detector allows the particles in the final state of the photoproduction (particle production from a photon beam) to be detected and identified.Due to the statistical nature of this kind of experiments, various final states are produced; for the purpose of this work as an educational resource, we are interested on reconstructing exclusive events with Λ 0 and K + in the final state.

Physical variables
In a relativistic particle dispersion event A + B → C + D, the description of the state of each particle is made in terms of four-momentum vector P = (P x , P y , P z , E), where P x , P y and P z are the components of the three-dimensional momentum vector P , and E is the total energy of the particle. 2ith the above, the initial state of the reaction is described by the sum of its four-momentum vectors: P i = P B + P A , and the final state: Having two particles in the final state, a two-body reaction can be characterized, in general, by six degrees of freedom (NDF).Imposing energy and momentum conservation, which defines four constrains, NDF is then reduced to two.Therefore, the description of a two-body reaction relies on two independent variables.Typically, cos θ and φ fulfill such a requirement; however, rather than working with two variables, it is frequently convenient to use invariants with one constraint.The most widely used choice are the well-known Mandelstam variables, s, t and u, defined by: s = (P A + P B ) 2 = (P C + P D ) 2  (1) where s is the square of the energy of the center of momentum (which also matches the square of the invariant mass of the reaction, a measure of the amount of energy available for the reaction [7,8]), t is the square of the momentum transfer of A to C and B to D (refers to the amount of energy transferred from one particle to another) and u, analogous to t, coincides with the momentum transfer of A to D and B to C. The physical significance of the Mandelstam variables deals with the fact that they are "Lorentz invariants", which means that they have the same value in any inertial frame.These variables allow to verify the consistency of the proposed simulation with the laws of physics, since they must have the same value, for example, in the "laboratory" and in the "center of momentum" of the reaction reference frames.In a real experiment, the Mandelstam variables are an important tool to reduce and limit the reaction events that meet these invariant quantities.
The experiment illustrated in Fig. 1, fixes the zenithal direction on z axis, along the photon beam.Axis y is set perpendicular to the surface of the experimental area where CLAS detector is located.For this case, the azimuthal angle and the polar angle are used as variables to characterize the directions of P and the corresponding angular distributions of final-state particles.
Importantly, the angular distributions can be described from the angles between a fixed reaction plane and decay planes (generated by final-state particles), which turn out to be unknown event by event; in other words, decay planes are randomly generated.Figure 2 outlines the reaction plane for initial-state particles (photon γ and proton p), and for two examples of possible reaction planes of finalstate particles (positive kaon: K + , and lambda: Λ 0 ).

Threshold energy
One of quantities of interest when studying reaction processes is the minimum required energy (or threshold energy) to produce a certain final state in the laboratory.This quantity allows to project out the architecture of a detector system for an specif variety of experiments.
The threshold energy is calculated by using energymomentum invariant P 2 = E 2 − P 2 = m 2 0 , and from  the properties of the particles involved in the process (see Table 1).From a two-body process A + B → C + D, we define A as the photon, γ, which describes the particles that makes up the photon beam in the experiment; B as the proton, p, which specifies the target of hydrogen at the geometric center of the CLAS detector; and C, D, as the final-state particles of the process: K + , Λ 0 .
Assuming the above, the initial-states four-momentum is expressed as P i = P γ + P p , where the photon and proton four-momentum vectors are set to P γ = (0, 0, E γ , E γ ) and P p = (0, 0, 0, m 0,p ): the photon beam is fixed in z direction, the value of P z,γ = E γ since m 0,γ = 0, and, the value of E p = m 0,p due to hydrogen is at rest in the laboratory reference frame; thus, the initial state four-momentum vector and its square s are expressed respectively as with "LAB" ("laboratory") as a label to identify the reference frame where the calculation is performed.
Equally for the center of momentum (CM) reference frame, the final state four-momentum vector and its square s are expressed respectively as Since s is a Lorentz invariant, that is , the energy of the photon beam can be expressed as Given that, the threshold energy of the photon beam that produces final-state particles with nearto-zero kinetic energy should approximately be where the total energy of the final-state particles are reduced to rest energies. e1504-5 From reported mass values in Table 1, we obtain that E LAB γ, threshold ≈ 0.9739 GeV and s = 2.7076 GeV 2 . 3

Monte Carlo method
The kinematic variables describing final-state particles of the reaction depend on the experimental specific circumstances, for example, the kinetic energy of the beam, the rest energy of the particles or the chance of interaction.Therefore, every single event will decay, in general, in a different reaction plane randomly.It makes then necessary, for simulation purposes, to generate randomly reaction planes for each event.
Depending on the reaction and the experimental conditions, different distributions can be used: isotropic distributions, perhaps for corrections to experimental data, or distributions according to some phenomenological model in order to understand what happens in a particular reaction.
The generation of random decay planes allows, among others, to carry out acceptance and efficiency corrections to experimental data, or compare the signal obtained from an experiment with simulated phenomenological models through the generation of random variables or Monte Carlo method [10] .
The Monte Carlo method is a mathematical tool used to find solutions to mathematical, physical and engineering problems, through the generation of random variables.One way to obtain a Monte Carlo dynamics is through the Method of congruences, which is operationally defined as [11,12]: with x being the variable to generate and m, c and a, the modulus, the multiplier and increased recurrence.These constants must satisfy that 0 ≤ a < m, 0 ≤ c < m and 0 ≤ x 0 < m, where x 0 corresponds to an initial value (or seed) required by the algorithm to generate the random variable x n .The quality of these algorithms is based on their repetition periods, depending on the parameters involving.In the case of the method of congruences, the bigger the m, the longer the repetition period.Due to the interest of generating random numbers with long enough periods, we used professional 3 A more general expression for E LAB γ,threshold falls to . Here, j Mj is defined as the sum of the rest masses of all particles in the final state.
algorithms, as the ones implemented in the classes TRandom2 and TRandom3 in ROOT software infrastructure, whose periods of repetition are of the order of 2 88 and 2 19987 , respectively [13].

Implementation of the Monte Carlo method
The generation of random reaction planes has to do with generating four-momentum vectors of finalstate particles living in the same plane, accordingly to momentum conservation.For computational simplicity, and for the sake of this educational resource, the generation of reaction planes is made isotropically.
The ROOT platform has a class called TGenPhaseSpace, which allows to generate events in the center of momentum of the reaction.TGenPhaseSpace allows to generate particles with a three-dimensional momentum vector, uniformly distributed in cos θ and φ, respectively (or phasespace).
Each event occurs at random following the dynamics of the Monte Carlo technique, whose original code, proposed by James F. algorithm [14], called GENBOD, generates final state N-bodies events in phasespace.
TGenPhaseSpace class requires defining the fourmomentum vectors of the initial-state particles, P γ , P p , the energy of the center of momentum √ s, understood as P γ +P p , and the masses of the final-state particles.
Four-momentum vectors are implemented by TLorentzVector P(Px,Py,Pz,E) class, which allows the methods Generate() and GetDecay() to be invoked.It generates, for each event, four-momentum vectors for final-state particles.Figure 3 shows a code snippet of TGenPhaseSpace implementation for a photon beam energy E LAB γ of 2.0 GeV.

Estimation of the threshold energy
With the purpose of illustrating the Monte Carlo technique to estimate the minimum photon beam energy required to produce K + and Λ 0 in the final state, it is necessary to generate events within a certain range of energies, according to the physics describing the reaction.First of all, it is proposed to set initial state conditions to generate events with a uniform distribution for E LAB γ between 0 and 2 GeV.Then, values as total energy of final-state particles, E K + and E Λ 0 , and their three-dimensional momentum components are stored in histograms as a function of E LAB γ .Finally, the results of the simulation are corroborated with s, t and u distributions, and the theoretical value of E LAB γ,threshold , constituting the flowchart in Fig. 4. A fragment of the computational implementation can be seen in Fig. 5 7

. Results and analysis
Several tests to validate the simulation through the conservation of momentum and energy can be chosen.
Generally, distributions or histograms of the initial and final state Lorentz invariants, namely Mandelstam variables s, t and u, are compared in different reference frames.The most common ones are the laboratory frame (LAB) and the center of momentum frame (CM).

s variable
The histograms in Fig. 6 show the distributions of s referred to the initial-state particles in both LAB and CM frames.
As expected, the LAB and CM distributions are identical, which corresponds to a Lorentz invariant.In particular, the distributions for s are partially different for initial and final states.Recall that the simulation was conducted to estimate the required minimum value of the photon beam to produce finalstate particles.The event generation ranges E γ from 0 to 2.0 GeV.For the initial state, the value of s must start at around M 2 p = 0.88 GeV 2 (squared proton mass at rest), which is consistent with the results obtained in Fig. 6.Regarding to the final state, s value should start around ∼ 2.707 GeV 2 , which is consistent with E LAB γ,threshold = 0.9739 GeV.Although the values where s is computed for the initial and final states are expected to be different, it is important to note that if the value of s for the final state starts at ∼ 2.707 GeV 2 , the distributions of the initial and final state must be identical.Figure 8 shows a superposition of s distributions for initial and final states in CM frame, where regions above s ∼ 2.707 GeV 2 match to each other.

t variable
t variable represents the square of the momentum transfer in the reaction, from γ to K + (t γ,K + ) or from p to Λ 0 (t γ,K + ).For a fixed E γ value, event by event, the momentum transferred is not the same and therefore a particular form of t distribution is expected, as shown in Fig. 9.
The value of t turns out to be less than a certain value t max , different from zero.This characterize the physical limits of the reaction, which allows t to have physical meaning [16].
By doing an event counting from Fig. 9, we found that the events number is significantly smaller than those generated as initial-state events.Such drawback allows to consolidate the analysis in the sense of looking for a consistent explanation with the implementation and physics of the problem.
The issue here is that the computational implementation only initialize values of the initial state but not of the final state.Final-state particles are generated by TGenPhaseSpace when the center of momentum is sufficient to produce them, otherwise, the program associates a nan ("not a number") value to each of particle's four-momentum components.Since Mandelstam variables are squared fourmomentum sums, the class TLorentzVector has the "+" operator overloaded in order to do component-bycomponent sums.Adding a nan to a non-nan value,  the result is a vector with nan components.Upon visualization, the Draw() does not display events with nan values.
It is necessary, then, to implement some lines of code to solve the problem and see all the events generated.The fragment of code depicted in Fig. 10 allows to test for nan values and setting them to zero.
This way, event counting is used as a test to be checked in the analysis, and for validating the analysis itself.Figure 11 shows the results after setting nan values to zero.
The left inset panel matches the same t distribution in Fig. 9.

u variable
Analogue to t, the u variable represents the momentum transfer in the reaction from γ to Λ 0 or from p to + .Similarly to t variable, for a fixed E γ value in the event-to-event basis, the u value is not   expected to be the same.Figure 12 illustrates the distribution of u variable for the reaction.
u variable turns out to have, again, certain values to frame physical region of the reaction [16].For the case of γ + p → K + + Λ 0 , computational estimate for the minimum value of E LAB γ to produce final state particles K + and Λ 0 , by Monte Carlo method, is consistent with the calculated: E LAB γ,threshold ≈ 0.9739 GeV.

Conclusions
The results from an exemplification through an estimate of the threshold energy of a two-body reaction, as this article describes methodologically, allows to outline, in general, the process involved in the analysis of experimental data from a real experiment and how to validate it computationally and from the related physics.
Since physical variables of interest as Mandelstam variables are well defined within the simulation (Figs. 7 and 8), these variables protect the validity of the simulation that allows to conceptualize "in- variant" quantities, typical from special relativity.More complex analysis can be derived from Mandelstam variables, for example, to study scattering amplitudes from s, u and t variables.
Furthermore, the threshold energy from the simulation obtained from Fig. 13 reflects consistency with the simulation model γ + p → K + + Λ 0 .It should be noted that these results can be used from junior scientist to test conceptual HEP basis and to get familiar with scientific HEP terms.
In general, the visualization of the variables involved in the reaction allows the scientist to verify, extrapolate, and in addition, associate physical concepts with the output code information.It is extremely necessary to have those active/interactive learning processes realized, not only when facing a data analysis of a real experiment, but also in any learning-teaching process in physics education.
The production of educational resources for the initiation of young scientists in the production of scientific knowledge is sometimes scarce; documentation from dissertations and articles tend to fulfill the required training in research.Unfortunately, experience shows a challenge in the time required by researchers to learn, for example, on how to conduct a data analysis from a real experiment in high-energy physics.
The document presented here happens to be a synthesis that allows young researchers to have conceptual background for facing a first analysis of a particular reaction from the physics of the experiment, to structure a very preliminar data analysis of an experiment, and to expect derived problems (as the nan values in the section where t variable was studied).All this is expected to be materialized in reducing the young scientist's learning curve, which is reversed in saving time and costs for experimental research, and to improve better learning processes in scientific research.

Figure 2 :
Figure 2: Reaction and two examples of possible reaction planes for γp → K + Λ 0 in the laboratory reference frame.

Figure 5 :
Figure 5: Code snippet showing the computational implementation of TRandom3 to generate a photon beam energy distribution in the simulation.

Figure 6 :
Figure 6: s variable in LAB (left panel) and CM (right panel) reference frames for initial-state particles γ and p.

Figure 7
Figure7shows the s distributions for final-state particles in both LAB and CM frames.In particular, the distributions for s are partially different for initial and final states.Recall that the simulation was conducted to estimate the required minimum value of the photon beam to produce finalstate particles.The event generation ranges E γ from 0 to 2.0 GeV.For the initial state, the value of s must start at around M 2 p = 0.88 GeV 2 (squared proton mass at rest), which is consistent with the results obtained in Fig.6.Regarding to the final state, s value should start around ∼ 2.707 GeV 2 , which is consistent with E LAB γ,threshold = 0.9739 GeV.Although the values where s is computed for the initial and final states are expected to be different, it is important to note that if the value of s for the final state starts at ∼ 2.707 GeV 2 , the distributions of the initial and final state must be identical.Figure8shows a superposition of s distributions for initial and final states in CM frame, where regions above s ∼ 2.707 GeV 2 match to each other.

Figure 7 :
Figure 7: s variable on LAB (left panel) and CM (right panel) for final state particles K + and Λ 0 .

Figure 8 :
Figure 8: s variable overlap for initial (solid line) and final (blue dashed line) states in CM frame.

Figure 9 :
Figure 9: Events distribution for t variable.

Figure 10 :
Figure 10: Fragment of code to set nan values to zero.

Figure 11 :
Figure 11: t variable distribution with the computational implementation avoiding nan values.

7. 4 .
Momentum distributions for K + and Λ 0 as a function of the beam energyThe minimum energy particle production can be obtained from various ways in the simulation.The proposal here is to obtain it by studying the threedimensional momentum distributions of the finalstate particles.Figure 13 (top panel) shows the total three-momentum vector of K + , | P K + |, as a function of E LAB γ .As expected, the photoproduction of final-state particles starts around E LAB γ ≈ 0.97 GeV.Also, this behavior is repeated for the histogram associated with the total three-momentum distribution of Λ 0 , | P Λ 0 |, as a function of E LAB γ , as shown in Fig. 13 (bottom panel).

Figure 13 :
Figure 13: Three-dimensional momentum distributions for K + (top panel) and Λ 0 (bottom panel) as a function of E LAB γ