Contents

1 Introduction

Human beings may never understand the very nature of the universe, but regardless, we are driven by courageous inquiry. For thousands of years we have sought to understand Aristotle’s “atom”- the smallest, indivisible stone that comprises our brains, our bodies, and all the universe known to man. And though that universe continues to broaden with technology and time, the atom continues to shrink. The hunt to understand the nature of particles and their interactions has escalated into an enormous puzzle - with seemingly, just a few pieces lingering out of reach.

The Standard Model is the strongest link humans have found to grasp the nature of the universe. It describes the strength of interaction between all known particles; however, it is far from complete. Not only does it fail to incorporate the effects of gravity, but it also predicts the existence of a new field - the Higgs - which has motivated a new era of physics research. The Higgs field allegedly gives mass to all that it interacts with. It is the primary target of physics today, but if it is found (by the biggest and best linear colliders physicists can dream up), then there will be even more problems to confront.

The Higgs boson has a mass that is described in the Standard Model by an infinite, diverging term. Clearly, an infinite mass cannot be in agreement with experimentation. Supersymmetric theory is a leading candidate that eliminates this complication, which I will elaborate on later in this thesis.

Two massive particle accelerators currently in the works to search for the Higgs and other Supersymmetric particles are the Large Hadronic Collider (LHC) and the International Linear Collider (ILC). Both colliders have a unique approach to searching for these particles; the LHC will fire protons at one another at an energy of 14 TeV and will be completed in 2008, while the ILC is a proposed high-energy linear collider that will create e⁺/e^- collisions at extremely high energies – 0.5 to 1 trillion electron-volts (TeV).

While there are advantages to both approaches, the University of Colorado’s SUSY group is working to optimize a proposed electromagnetic calorimeter for the ILC. Our primary objective in optimizing this detector is to determine if our design can provide the needed spatial and energy resolution to infer the mass of the Higgs boson. One key factor in this search is identifying and measuring the decay products that emerge in interactions with the Higgs.

The Higgs often decays into hadronic jets and emerges with Z⁰and W bosons, which contain lots of charged hadronic particles and π⁰s. The charged hadronic particles can be removed by a tracking device, but the π⁰s decay further into two photons, and these photons hit the detector very close to one another at high energies. When this happens, it is difficult to tell if the “merged cluster” we see is due to a single photon with a large energy, or two smaller photons. It is crucial that this distinction be made accurately because the mass of the Higgs can only be reconstructed as accurately as the momentum of these particles can be measured.

My effort in this project has been to study how well CU’s proposed calorimeter design can measure photons. This thesis will outline key aspects of both the Standard Model and Supersymmetric theory motivating the construction of the ILC, describe our simulation, and detail the necessary steps for a photon to be isolated and measured. In addition, I will comment on our current resolution capabilities in making this measurement.

2 Motivating Photon/Pizero Separation Studies

The current lower limit for the mass of the Higgs boson is 109.7GeV/c². This was confirmed in 2002 experiments at LEP (Large Electron-Positron Collider) at CERN, where the 200GeV e⁺e^-→ Z⁰H process was studied. This means that the Higgs mass may, in fact, be much greater than 110GeV, but experimentally it has been proven that it cannot be less.

To arrive at this mass estimate for the Higgs, decay products of the Z⁰ were collected and its mass was reconstructed to 91.176+∕-0.0021GeV/c²[4]. This was done by summing the energy and momentum of daughter particles in the Z⁰ decay and plugging their total momentum into Einstein’s relativistic equation (Eq. 2.1).

(2.1)

Equation 2.1 can be used to solve for the Higgs mass for this interaction if it is re-written it in terms of properties of Z⁰ boson and the energy in the center of mass frame (Eq. 2.2). This can be expanded (Eq. 2.3) and simplified (Eq. 2.4). The mass of the Z⁰ is now well-known from experimentation, but each time the Z⁰ emerges, we must deduce its energy by measuring and adding the energies of decay particles to obtain the E_z term in equation 2.4. The detector must be very sensitive for this measurement to be made adequately.

(2.2)

(2.3)

(2.4)

Similarly, many Supersymmetric particles decay into the Z⁰and W bosons in the intermediate stage. To discover these particles or the Higgs boson, it is essential that a detector be sensitive enough to reconstruct the mass of the Z⁰ accurately by collecting the decay products (including hadrons, pizeros, and other particles) in this way. Each particle adds a momentum vector to the equation which makes my research – the separation of pizeros and photons- a very small but pivotal task in the search for the Higgs.

2.1 The Large Hadronic and International Linear Colliders

The LHC and ILC are two major particle accelerators that will be built in the next two decades. Each offers different advantages in searching for the Higgs and other new particles.

The LHC is thought to be a “means to discovery.” It operates at great energies and it may be able to see new particles first, but the ILC is the true precision machine; it will be able to take more detailed measurements when these particles are found.

The LHC in Geneva, Switzerland, can rival the high energies at which the ILC will operate. It operates up to 14 TeV, but unlike the ILC, the LHC fires protons at one another. In a proton-proton collision, six quarks are involved but only two collide. In this type of collision, it is not possible to tell in any one event which quarks collided (since quarks have differing charges of -1/3 or +2/3), and it is not possible to tell which passed through with the beam, so the total energy of the collision cannot be determined. It is possible that the Higgs may be discovered at the LHC, but even if it is, it will be extremely difficult to decipher any of its properties because of the complexities in the resultant decay.

In contrast, it is electrons and positrons that will be fired at one other at the ILC, and in every collision the total energy and momentum is known. This will make it possible for a better estimate of the Higgs mass and momentum. These leptons will be fired at one another through a 35-kilometer tunnel that accelerates particles to an energy of 0.5TeV in the first stage (with a completion date set in the 2020s). An upgrade of this project is planned for the 2030s, which will increase the collision energy to 1TeV. It is hoped that SUSY particles, or the elusive Higgs boson can be observed and measured in this massive collider.

3 The Standard Model

The Standard Model describes the interactions of all elementary particles through the three fundamental forces - electromagnetic, strong, and weak nuclear. While it predicts and describes these interactions with extraordinary precision (in agreement with experimental data) there are problems with the theory. For instance, it fails to incorporate the effects of gravity. In this discussion, we will consider only the aspects of the Standard Model that motivate research in Supersymmetry. In doing so, we will best highlight the goals of the SUSY group in optimizing our proposed detector for the ILC.

The Standard Model describes how different particles react based on the principle of Gauge invariance. Theory groups elementary particles into two classes: fermions and bosons.

Fermions have half-integer spin and make up all known matter. They obey Fermi-Dirac statistics and the Pauli-Exclusion principle stating that no two identical fermions can occupy the same quantum state simultaneously. Fermions can be further divided into leptons and quarks, and an identical partner of opposite charge exists for each such that for every lepton there is an anti-lepton, and for every quark there is an anti-quark. There are six known leptons, three of which are neutrinos; these are classified into doublet pairs, each containing one neutrino (Fig. 1).

So far, experimentation has revealed six quarks with charge of -1/3 or +2/3. Quarks bind together to create composite particles known as hadrons including mesons, which consist of quark/antiquark pairs, and baryons, which consist of a triplet of quarks bound together. In addition to this triplet of “valence” (high-energy) quarks, baryons and mesons can contain “sea quarks”, which are referred to as virtual particles, because they exist only for a limited period of time and space. Because of this, the Heisenberg uncertainty principle mandates that the energy and momentum of sea quarks cannot be defined with unlimited precision.

In contrast to fermions, bosons have integer spin and participate in the propagation of the weak, electro-magnetic, and strong nuclear forces. Bosons are exchanged when any of the three fundamental interactions takes place. Through such experiments as the photo-electric effect and Rutherford scattering, the pioneers of quantum physics discovered that electromagnetic forces between charged particles are mediated by a photon, a mass-less particle. In contrast, strong nuclear interactions are mediated by mass-less bosons known as gluons, discovered at the DESY laboratory in Germany in 1979 [1]. The third force predicted and described by the Standard Model is the weak nuclear force, which dictates interactions such as beta decay, meditated by the exchange of charged W+/W- and neutral Z⁰ bosons. These were discovered at CERN in 1983 [1]. See Figure 1 for a table of bosons, and Figure 2 for a diagram of the interactions between bosons and fermions.

---

Figure 2: Standard Model Interactions[14]

---

3.1 Postulation of the Higgs

Since the Standard Model was first postulated, every performed laboratory experiment has been in extraordinary agreement with it. In 1964, Peter Higgs, François Englert and Robert Brout first recognized that a mass term could be obtained from a massless potential if it was written in the form μ²φ² + λφ⁴.

The mass term in the equation for the Higgs was added to the LaGrangian for spin-zero particles (Eq. 3.1). Typically, a mass term would be of the form “+μ²φ²” ; however, the spin-zero LaGrangian is not gauge invariant if a mass term is included, so such a term cannot exist. To investigate this peculiarity, Peter Higgs made the mass term negative and included it in the expression for the potential (Eq. 3.1).

(3.1)

We can find the minimum of this potential by taking the derivative and setting it equal to zero, arriving at Equation 3.2:

(3.2)

Solving this for φ, a minimum can be found at ν = +∕-. This value is also known as the vacuum expectation value of the Higgs field φ. To study this minimum, we let φ = ν + η(x) and a mass terms appears of the form in m_n² = 2λν² = -2μ² [2].

This term is considered to be the mass of the Higgs boson. It is believed that the interaction with this field is what gives particles mass. This transformation has been interpreted by some to suggest a historic phase transformation in the universe[3]. Many physicists believe that at the beginning of the universe (the time of the Big Bang), all forces were one. It is thought that the Higgs field became manifest when the universe cooled, and the universal potential shifted from zero to. At this time, all particles were given mass[3] (Fig. 3).

Field Theory predicts that as fields propagate, they also dissociate and recombine as they move through space. This dissociation of the boson is into particles that usually reunite into the same boson. For example, the W^- may dissociate into an electron and an electron-neutrino; in turn the electron and the electron-neutrino may recombine to form a W^- boson.

The discovery of the Higgs remains a critical link in validating the Standard model, but if it is found, it will bring about a concern for particle physicists because the mass of the Higgs has a term described by an infinite series of fluctuations which contain a quadratic, diverging term - due in part to the fact that the Higgs is spin zero, and can interact with itself [7]. This term is problematic in that the predicted infinite mass cannot possibly agree with experiment. A grander theory is needed to eliminate the diverging mass of the Higgs; Supersymmetry is such a theory and offers an elegant solution.

4 Supersymmetry

Super symmetrical (SUSY) theory introduces a theoretical framework that compliments the Standard Model by introducing a SUSY “super-partner” for every known particle. The spin of the super-partner differs from the spin of the particle by 1/2 ℏ, so that when bosons are exchanged in field interactions, the super-partner is also exchanged. As a result of this difference in spin, the vacuum fluctuations of the Higgs (that lead to a quadratic mass divergence) are cancelled by a negative quadratic divergence produced by the vacuum fluctuations of the Higgs into SUSY particles[6].

Another interesting fallout of SUSY is that the theory introduces a quantum number which demands the existence of a particle that may account for the 25% of unidentified mass in the universe referred to as “dark matter.” This comes from the fact that SUSY theory allows for proton decay.

Since protons are known not to decay experimentally, the theory introduces a new quantum number called “R-parity,” which voids the most likely proton decay. R-parity will be defined shortly, but for now, I will highlight the important consequences of this new quantum number.

R-Parity has a dependence on the spin number (j) , and on the baryon(B) and lepton numbers (L) (Eq. 4.1). For R-Parity to be conserved, SUSY requires that lepton and baryon number also be conserved in all interactions. Lepton number is the relationship between the number of leptons (N_l) and anti-leptons (N_l) in an interaction (Eq. 4.2), while baryon number (Eq. 4.3) is similarly the relationship between the number of quarks (N_q) and anti-quarks (N_q) involved. An example of lepton conservation is provided in Figure 4.

(4.1)

(4.2)

(4.3)

---

Figure 4: Lepton conservation requires that the number of total leptons before a decay matches that after the decay.

---

R-parity introduces the restriction that Supersymmetric particles must always decay into a Standard Model particle and another SUSY particle. This leads to the postulate that there must exist a very weakly interacting particle which cannot decay further - this particle is known as the neutralino and is the super-partner of the photon. Supersymmetry theorists believe that the neutralino is the dark matter particle[9].

There is one final aesthetically pleasing argument supporting the unification of fields through the implementation of Supersymmetry. If Supersymmetry is neglected, and one extrapolates the inverses of the coupling constants for the electromagnetic, strong, and weak field interactions against energy, then the lines do not meet at one point. However, if these values are extrapolated including the effect of SUSY particles, then the three lines meet at a single point (10¹⁶ GeV), which happens to be the energy of the Big Bang (Fig. 5).

If the Higgs is found, a modifying term (which SUSY provides) is needed to deal with the diverging mass term introduced into the Standard Model. Thus, research in Supersymmetry will advance greatly if the Higgs, or any SUSY particle can be observed and provide the first validation of this theory.

5 Cascade Showers and Calorimetry

Electromagnetic calorimeters are built to measure the energy deposited when incident electrons and/or photons penetrate into their volume and produce an electromagnetic shower. This shower results in a multiplicative process with a temporal extent determined by the magnitude of the energy of the incident particles. High energy electrons lose most of their energy through bremsstrahlung radiation, an electro-magnetic radiation produced when a charged, decelerating particle is deflected by another charged particle¹. This effect results in the production of high-energy photons, which then undergo Compton scattering or Bethe Heitler pair production, a process through which an electron-positron pair is created when a high-energy photon interacts with the nucleus of the atom. Figure 6 depicts the processes of pair creation and bremsstrahlung radiation.

The process repeats when these electrons and positrons radiate more high energy photons which also convert to electron positron pairs, so long as their energy is sufficiently large. This cycle is often referred to as a cascade shower. In the final stage of this process, photons produced interact with the calorimetric medium through the photoelectric process, while the electrons dissipate their energy through collisions with the medium [8].

5.1 Properties of Cascade Showers

The radial distribution of a given cascade shower is described by the Moliere radius, given by Equation 5.1. The Moliere radius refers to the radius inside which 91% of the energy of the shower is contained. It is expressed in terms of the multiple scattering energy (E_s), the mean path length of an electron in a material (radiation length X_o), and the critical energy of the shower (E_c), which is the energy at which an electron loses as much energy through collisions as through radiation.

The multiple scattering energy (E_S), used in the calculation for the Moliere radius, is due to several small coulomb scatterings that may take place, wherein an incident charged particle (e⁺∕e^- ) is scattered slightly when it interacts with a proton in the material of the medium². When this interaction occurs, a photon is passed from medium to the charged particle, which causes the charged particle to scatter slightly off its initial path . This exchange might happen several times, and the end deviation of several such scatterings from the natural “unscattered” path can be described by a Gaussian (Fig. 7).

When the transverse momentum of the photon exchanged in this process is very large, we may see a single scattering that is greater than the “random path” summations over successive smaller scatterings. This is known as a single Coulomb scattering and contributes a tail onto the Gaussian for multiple scatterings. Since significant single single scatterings happen so infrequently, they make no accountable contribution to the E_s variable in the Moliere radius. However, it is because of single Coulomb scattering that we occasionally observe a shower with an energy distribution larger than the Moliere radius.

The radiation length (χ₀) is given by Equation 5.2, where N is the Avogadro constant, A is the atomic weight of the medium, r_e is the electron radius, and ζ is the necessary correction for atomic electrons in the overall bremsstrahlung process [15]. The Moliere radius for CU’s proposed tungsten-scintillator design (elaborated on in the section 6.1) is calculated to be 2.5 cm.

(5.1)

(5.2)

6 Optimizing a Detector

There are currently three competing detector designs being considered for implementation at the ILC. The smallest proposed detector is the Si-D design, which incorporates a Tungsten-Silicon electromagnetic calorimeter (ecal). The second is the LDC (Large Detector Concept) design which offers a bigger tracking chamber, and implements the same ecal as the Si-D detector. The third design is the Global Large Detector (GLD) design which is the largest detector, advocating the use of CU’s Scintillator-Tungsten ecal design.

The advantage of a Silicon-Tungsten detector (as in the Si-D and LDC designs) is that it delivers better spatial resolution because the design consists of 1x1 cm tiles, and thus has several points of signal readout. In contrast, the design being researched at the University of Colorado uses a cheaper material - scintillator instead of silicon, and implements 5x5 cm tiles interleaved so that offset layers form a brick-like pattern giving the resolution of a 2.5x2.5 cm tile. CU’s calorimeter design then requires only 1/4 the number of readouts of a 2.5x2.5 tile design. Since GLD design is the largest³ and most costly of the three, the GLD group is pushing to cut back on costs by implementing CU’s cheaper scintillator-tungsten prototype.

One advantage to the Scintillator detector is that it offers better energy resolution than the silicon detector. Experimentally, Silicon detectors have exhibited energy resolution with an error of 17% (Eq. 6.1), while Scintillator detectors typically improve resolution to an error of 11% (Eq. 6.2). ⁴

(6.1)

(6.2)

6.1 The University of Colorado ECAL Design

---

Figure 8: Detector Schematic: The electromagnetic calorimeter is the central cylinder (blue)

---

The Tungsten-Scintillator ecal design being optimized at CU is that of a large, cylindrical detector, consisting of 40 layers of alternating Tungsten/Scintillator material parallel to the beam-pipe (Fig. 8). Each alternating layer of Scintillator tiles is offset in a brick-like pattern (so that the corner of the tiles in the first layer line up with the center of tiles in the next layer) to improve overall spatial resolution of the calorimeter.

The effect of offsetting alternating layers of tiles is to improve the effective sensor resolution from 5x5 cm to 2.5x2.5 cm. We assess the energy distribution of a cluster by examining how the energy spills into several tiles. Therefore, we observe the worst spatial resolution when the hit is in the very center of two tile edges, because most of the cluster energy is read out by one tile measurement instead of a spread of a few different tile measurements. When this happens, we can infer the least about the general shape of the energy distribution of the cluster. Tile geometry and this point of “worst separation” can be seen in Figure 9.

Each Tungsten-Scintillator tile is seven mm thick (Fig. 10). When a layer of Tungsten is struck by high-energy photons, a shower of e⁺e^-pairs is produced. These charged particles propagate though space and hit the scintillator, exciting electrons in orbit in the scintillator molecule. This produces light that can be observed by electronic means. The electronic signal produced is read out by a small Silicon photodetector (SiPD) on each tile. This signal is proportional to the energy in the shower and is used to determine the energy that hits each tile. After being measured, the photons then propagate through the next layer of the detector and strike the Tungsten, knocking off more e⁺e^-pairs.

The cascade continues until all of the energy in the shower is deposited into the calorimeter. The energy deposition patterns created by hits in the electro-magnetic shower are unique for each particle. They can be used to identify the parent particles in the decay, and to determine the incident angle of each particle involved.

---

Figure 10: Calorimeter Tile Design

---

7 Simulation

Simulations of the showers were generated using SLIC (Simulation Tool for the Linear Collider), a GEANT-4⁵ based simulations package developed at the Stanford Linear Accelerator Center (SLAC). All the events were generated using a geometry file which models the Large Detector Concept (LDC) incorporating our design for a scintillator electromagnetic calorimeter (ecal). The reconstruction and analysis software is written in java, incorporating the hep.lcd framework [5]. The analysis is presently carried out using JAVA Analysis Studio Version 1.6 and ROOT.

An extensive database has been constructed using this software by simulating data for photons, and for π⁰events. Each data file contains the energy deposition in each calorimeter tile for 10,000 identical events. For example, one data file might contain 10,000 identical photon events at an energy of 20GeV, fired at the same θ, φ (CU’s ecal design requires that a group of tiles aligned in any radial direction occupy a constant θ and φ value. See Figure 12).

Data files exist for photons and π⁰s with energies ranging from 1-100GeV at 25 separate points on the tiles that are used for reference. These 25 points are spaced equally throughout the bottom left quadrant of the tile. By symmetry, this tile could be any tile in the detector, and this quadrant can be flipped any way so that any of the corresponding 25 points in the other three quadrants can be simulated (Fig. 11).

The energy deposition patterns for points in-between these are simulated using interpolation code written by SUSY researcher Joseph Proulx in 2004. Previous CU undergraduate students Maria Gulda and Sarah Moll laid the groundwork for this database, creating all data files up to 75GeV. I have finished creating 75GeV and 100GeV data for photons and π⁰s at each of the 25 points used for analysis. By utilizing these simulation data files made for various shower energies and reference points, we can use interpolation, rotational, and statistical techniques to reconstruct the expected energy observed in each tile for any photon or π⁰ shower less than 100GeV.

For each event, our simulation code looks at a group of 3x3 simulation squares (with the central tile containing the energy-weighted center of the cluster). We do this for each layer throughout the calorimeter. By studying the mean distribution patterns in these towers, we can begin to separate the hits due to single photons from those due to π⁰ decays. This is accomplished by comparing known mean tile energies (reconstructed from database data) to the energy distribution of any single event.

---

Figure 12: Each radial stack of tiles occupies a constant angular width in θ and φ. The arrows mark the contact point of incident particles on the detector surface. The red dots show the maximum energy tower, and the green dots mark peripheral energy deposition.

---

7.1 Calorimeter Resolution

Before we can begin to ask how well CU’s ecal can resolve individual particles, we must investigate how effective the 40 layers of Tungsten/Scintillator are in capturing all of the energy from the shower.

Using 75GeV photon data simulated by GEANT-4, we can plot the total energy collected in the calorimeter against the number of events (Fig 13 ). A deviation from the expected Gaussian can be seen in the lower energies, resulting in a lower mean energy of 72.309. This deviation, or “tail” can be attributed to “punch-through” energy that made it all the way through the calorimeter’s 40 layers and escaped.

In the current 40-layer configuration, “punch-through” energy loss will therefore become significant in high-energy events. Adding more layers of tiles to the calorimeter would be a costly design alteration, so I investigated the plausibility of correcting for this energy loss without modification of the sensor design.

If the calorimeter is perfect, we should expect 11% error in our energy resolution (Eq. 7.1). If Equation 7.1 is expanded to include a linear constant representing the punch-through error in the detector, then it can be represented by Equation 7.2. I plotted sigma-to-mean ratios for data at 1, 5, 10, 20, 50, 75, and 100 GeV to examine this relation (Fig. 14), and fit a curve to the data points to solve for constants A and B in Equation 7.2.

If we let A equal 0.114, then this curve is best fit when B equals 0.015, suggesting to us that 1.5% of the energy will escape the detector.

(7.1)

(7.2)

Figure 14: The plot below shows the trending of the ration of energy observed against increasing energy. This line is best fit by Equation 7.2 with A=0.114, and B=0.015.

7.1.1 Energy Resolution and Correction

To account for this 1.5% energy loss, we will add a correction factor to the simulated observed energy. If we plot the energy in the last five layers of the detector against the difference between the expected energy and the observed energy, then we get a roughly linear correlation (Fig. 15).

---

Figure 15: Trendline for Missing Energy

---

When this plot (Fig 15) was regenerated for 50, 75, and 100GeV data, all energies were found to be roughly consistent with the linear fit y=0.7x+1.7. This suggests that the missing energy is roughly proportional to the energy in the last five layers divided by 0.7.

I began with an initial correction given by Equation 7.3, and replotted the Gaussian of Figure 13.

(7.3)

The mean best matched the expected 75GeV when the entire correction is multiplied by .98 (Eq. 7.4). The result (Fig. 16), has a new mean corrected energy of 74.998.

(7.4)

This correction gives the same result (a better fitting Gaussian) for data at 1, 5, 10, 20, 50, and 100 GeV . This correction can be applied to data seen by the detector to more accurately report the energy of the cluster. Finally, if we re-plot Figure 14 with sigmas and means for the corrected energies, the constant B (energy unaccounted for) drops from 1.5% to 0.8%. However, the overall energy resolution (constant A) increases from 11.4% error to 11.9%.

Still, it can be said that with this correction implemented, the total energy of the cluster is better represented, and we can begin to analyze energy deposition patterns signature of different particles.

---

Figure 16: Corrected Cluster Energies for 75GeV Data vs. Number of Events

---

8 Particle Separation

Using SLIC (the GEANT-4 based simulation package), we analyze energy distributions on the calorimeter by utilizing a number of tools. This process consists of three stages that will be detailed in the following sections. First, the charged hadron tracks must be identified and removed. Secondly, a χ² calculation must be performed on remaining clusters, comparing them to the energy distributions in the database to decide if these clusters are single photons, or double photons from a π⁰ decay. Once this distinction is made, these photons must be measured.

8.1 Removal of Charged Hadron Tracks

One primary aim of CU’s research with a Scintillator-Tungsten detector is to show that the Scintillator-Tungsten design detector offers spatial resolution that is adequate in reconstructing the masses of the W and Z⁰ bosons. Both of these particles frequently appear in decays where the Higgs boson and SUSY particles emerge, so it is essential that a new detector have an electromagnetic calorimeter that can reconstruct these masses with sufficient precision. As mentioned in Section 2, an e⁺e^- pair can interact and form a Z⁰ and an H⁰. In this interaction, the Higgs mass can only be inferred as well as the energy of the Z⁰ is reconstructed. It is therefore important to be able to measure the daughter particles from a great number of different interactions that might take place.

The Z⁰ boson decays into leptons only 6% of the time. Leptons are easily recognized since their masses and charges are well-known. Three percent of the time, the Z⁰ decays into a τ⁺τ^- pair, which decays further into leptons that are well-known, and neutrinos which escape the detector completely. However, 70% of the time, the Z⁰ decays into hadrons - both charged and neutral. If the energies of the Z⁰ and the W are to be reconstructed accurately, then hadrons need to be collected and measured with great precision.

Charged hadrons have slightly curved trajectories that can be seen by the tracking chamber and then separated from the shower. Figure 17 depicts a general decay shower as it first appears on the detector surface. Charged hadron trajectories tracked by the tracking chamber can be traced out onto the calorimeter and their hits can be removed. Currently, the SUSY programs are able to do this with 88.6% efficiency (meaning that we identify 86.6% of the charged clusters that are present in the Monte Carlo data), and with 87.8% purity (meaning that of these clusters found, 87.8% of them are correctly identified as “charged” clusters. The others are actually neutral clusters wrongly identified).

The neutral particles remaining after this process include K⁰ and π⁰s and neutrons. K⁰s can decay into a combinations of neutral or charged pions and well-known leptons. These charged pions can also be removed by the tracking chamber, so the primary particle left behind is the π⁰. The neutral π⁰ has a lifetime of 8.4x10^-17seconds, thus, it nearly instantaneously decays into two photons. We can only reconstruct the energy of the incident π⁰ if we can accurately separate and measure these two photons.

8.2 Creating Covariance Matrices

After the charged hadron tracks are removed, we need to measure the energies and positions of the two photons produced by π⁰ decays. However, we can’t begin to separate and measure the two photons in these decays until we test to be sure that we are actually looking at two photons and not just one. To decide if an event is a merged cluster or a single photon event, we will make a χ² calculation which compares the measured tile energy depositions to mean tile depositions from known events (calculated from database data). However, covariance matrices - which demonstrate the energy correlations between the tiles - must be created to be factored into this χ² calculation.

When examining a cluster, our software returns the tile that contains the energy-weighted center of that cluster, as well as the eight surrounding tiles. We refer to this three by three array of tiles stacked throughout the 40 layers as an “energy tower”.

Covariance matrices are “trained” over thousands of events so that they best represent the mean energy correlations between tiles. First, we isolate offset stacked layers in groups of three (i.e., the first group - mini-stack one - will be layers 0,2, and 4, and the second group - mini-stack two - will be layers 1,3, and 5, etc. See Figure 18).

Then, each mini-stack is rotated in a standard way (described in Section 8.3). After this rotation is made, covariance matrices can be trained for each group of 9 tiles throughout these three layers, so there exists one covariance matrix for every three layers in the detector. Each element in this 9x9 matrix then represents the energy correlation between two tiles of a 3x3 array (which is a summation over three layers thick).

The elements in the covariance matrix are calculated by Equation 8.1 where μ is one tile’s mean (for the mini-stack), calculated from thousands of photon events. The indices i and j are tile elements, and x is a single energy measurement in one tile.

(8.1)

A covariance matrix is calculated for each “mini-stack” with matrix values that represent the average tile correlations in that mini-stack.

8.3 Rotation

Since only 25 reference points of photon incidence are kept in our database, a rotation algorithm is used to standardize each individual event to look like one of these points.

The energy-weighted center of a 3x3 “tower” will always be in the central tile stack, but it may be in any of its four quadrants. The rotation code flips tile matrix elements either vertically or horizontally (or both) to move this energy-weighted center to the top-left quadrant of the central tile (Fig.19).

The 25 reference points that we use actually fall in the bottom-left corner of the central tile (refer back to Fig. 11). Before standardized means and covariances are calculated, the reference data is also flipped to concentrate all energy in the top-left corner. This process results in a standardized tower matrix that can then be compared to also-standardized means and covariances. The energy distribution observed over the nine tiles after this rotation is made will then look similar to Figure 20, which graphs the energy for each tile in a 3x3 square.

This rotation algorithm makes it possible to generate a finite amount of data (for only 25 points), and use symmetry to map these points into the other three quadrants. With this limited data set, we can now make a χ² comparison between a measured energy cluster (most likely due to two nearby photons from a π⁰ decay) and the energy cluster of a single photon. This comparison will be a sufficient means of determining if the measured cluster is from a single photon, or two photons.

8.4 The Chi-Squared

A χ² calculation is the best known method to decide whether or not an observed cluster is a single photon, or two photons. We utilize this calculation to compare the energy distribution of a single, measured event to a set of means calculated from known photon events. If the measured event is due to a single photon, then the energy distribution will be similar to the means and the χ² value will be small. Conversely, if the measured event was a π⁰ decay (meaning two photons are present in the cluster) then the energy distribution will not resemble the single photon energy means very well, and the χ² calculated will be large.

There are several different ways to go about calculating a χ² , especially for the case we are considering which includes energy measurements in many tiles and layers. Different χ² methods may lead to better or worse photon/π⁰ separation resolution, but from a statistical standpoint, all methods are correct.

A traditional χ² can be found by taking the difference between the measured value and the mean value, and dividing this difference by the square of the standard deviation (Eq. 8.2).

(8.2)

If a measurement is similar to the mean values to which it is compared, then the χ² will have a value equal to the number of degrees of freedom. If we are considering the case with nine degrees of freedom (in a 3x3 block of simulation tiles - for one layer only), a χ² value of nine can be expected.

Going a step beyond this, we calculate a χ² for not just one layer, but for an entire energy tower. For best resolution, the 40 layers have been broken down into groups of three with tile layers that align vertically (for clarification, refer back to the geometry discussed in Section 8.2). For each of these groups, one χ² is calculated. We have 14 different groups of three layers, meaning that we have 14 different χ² calculations to sum. Each group or “mini-stack” contributes nine measurements (for the nine tiles in the matrix), so we can expect a total χ² for the tower of around 120 (9x14).

We can further improve the resolution of the χ² by considering the correlations between the tiles. The best resolution has been attained by calculating the χ² using a Hessian Matrix (H_ij), which is the inverse of the covariance matrix (V^-1)_ij discussed in Section 8.2. For one mini-stack, with an expected χ² of nine, this is given by Equation 8.3. This calculation is made for each of the 14 different groups of layers, and these values are summed together.

(8.3)

Figure 21 demonstrates what a χ² distribution looks like for 10,000 photon events versus 10,000 π⁰ events fired at the same point on the detector. The photons typically have a much smaller χ² than π⁰ events.

It is harder to separate between photon and π⁰ events at higher energies, when photons in a π⁰ event are Lorentz boosted and cause the two photons to hit closer together. This causes the event to look more like it is from a single photon, and these π⁰ χ² values are characteristically much smaller (Fig. 22).

---

Figure 21: Overlain for 10,000 events. Single photon events are in blue and π0events are in red. Left : a 20GeV (e+e- ) collision for the central point in the tile. Single photon events have a mean of 122.95, and π0 events have a mean of 522.91. Middle :a 50GeV (e+e- ) collision for the central point in the tile.

---

8.5 Separation Capabilities

The efficiency of separation we can achieve based on χ² plots (Fig. 21) varies not only with energy, but also with the location of the hit in the tile. The efficiency plot (Fig. 23 ) shows efficiency curves for both photon and π⁰events. The calculated χ² for each event is plotted on the x-axis against the efficiency at that point. For photon events, the efficiency is calculated by dividing number of events below a χ² value by the total number of photon events. Similarly, the efficiency of π⁰ events is found by dividing the number of π⁰ events above a χ² value by the total number of π⁰ events. The cut efficiency varies slightly by point, but remains around 90% for 20GeV photons, and the separation efficiency drops to roughly 70-80% for photons of 50GeV.

It is much more difficult to make this separation for higher energy events. However, over 90% of the photons that emerge from Standard Model and SUSY e⁺e^- events have energies less than 100GeV (Fig. 24).

---

Figure 23: Efficiency plot of photon events (red) vs pizero events (blue) for 20GeV events. At 20GeV(left), photon and π0 events can be separated with 94% efficiency. At 50GeV (right), the efficiency for this same point drops to 67%.

---

---

Figure 24: Plots of the number of of e+e- events (where photons and π0s emerge )vs particle energy. Left: Supersymmetry events. Right: Standard Model Events. Photon events are dark blue and Standard Model events are shown in light blue. Over 90% of Standard Model events do not produce photons with energies greater than 100GeV.

---

9 Merged Cluster Measurements

The χ² separation method implemented in simulations at CU will determine the particle type with 95% efficiency at 20GeV, but it is an entirely different and more difficult task to go on to measure the energies and positions of the these photons. No research group has yet claimed to be able to reconstruct the π⁰ energy without first making the unrealistic assumption that the two photons can be measured independently.

Therefore, the following method devised by myself and CU SUSY researcher Joseph Proulx is an entirely rudimentary and experimental idea. The method employed is still in need of refinement to reduce computation time and error margins, but the results yielded are indicative of a successful means of measuring two photons in a merged cluster.

9.1 Separation Method Currently Being Explored at CU

I have spent the bulk of my research investigating the analysis that follows the photon-or-π⁰ cluster diagnosis discussed in Section 8. If the event is a single photon, it is passed to trajectory reconstruction software written by undergraduate SUSY researcher Jack Gill. If the event is a merged cluster due to a π⁰ decay, it is to be fed into a JAVA program I wrote which uses χ² minimization tactics to try to identify the most likely energies and positions of the two photons present in the cluster.

To try to make this determination, this fitting program constructs a “guess” array energy tower to compare with the observed energy tower. It does this by looking up the mean energy distributions for two single photon towers and adding them together, layer by layer. The program then computes a χ² for comparison, and varies the “guess” photon parameters to make new comparison arrays until the reconstructed energy tower best matches the observed energy tower (when this happens, the χ² will be minimized).

9.1.1 Constructing a Comparison Array

To construct a “guess” energy tower to compare to the observed energy distribution, it is necessary to look up the towers for two single photon events, and add them together. The mean energy distributed in each tile due to a photon decay can be looked up in the database for any one of the 40 layers by specifying an energy value (in GeV) and a θ and φ position of incidence for that photon. In this sense, an energy tower (for 40 layers) is just a function of these variables. For a single layer, the energy matrix (M_π⁰) for a π⁰ decay then can be thought of as a sum of the energy matrices for two separate photons (Eq. 9.3).

(9.1)

(9.2)

(9.3)

The θ and φ specified offer two unique pieces of information. First, these coordinates correspond to a specific tile in the detector, and secondly, they tell us which quadrant of that tile the hit was in (the energy-weighted center of the cluster). From this information, we can look up an array for that layer, and rotate it appropriately so that the center of the hit corresponds with the specified coordinates.

---

Figure 25: Two photons arrays are looked up based on an energy and a position. The position is specific to the tile and quadrant, so they can be rotated correctly, and added together to simulate the array of the merged cluster.

---

9.1.2 χ² Computation to Measure Merged-Cluster Parameters

Once a “guess” energy array is constructed, a variance array is created in the same way - by adding the variance arrays for each individual photon together. Each element of the variance array is then just the combined variances for that tile due to each of the photons used as “guess” parameters. With this information, a χ² for the cluster can be computed by comparing the observed energy tower with a reconstructed tower (Fig. 25). This is done by subtracting the guessed mean from the observed mean for a single tile, squaring the difference, and then dividing by the new variance array (Eq. 9.4). A χ² is computed for every three alternating layers (just as before), and these results are added together to get a χ² for the tower.

(9.4)

The goal of this program is then to systematically vary the energies and positions of the two trial photons to assemble guess arrays until this χ² measurement is minimized when the best parameters are found. The most effective means I have found in doing this (thus far) requires carefully setting limits on the range of energies, θ, and φ values, and then looping though these possible coordinates at increments just large enough to limit computation time within a feasible level.

9.1.3 Loop Restrictions Due to π⁰ Decay Dynamics

The χ² minimization method of my photon measurement program employs a quintet of nested loops that vary the values in θ and φ for both photons as well as distribution of energy split between the two photons. It is essential that limits be placed on both the step size and the angular separation transversed in the coordinate loops to mitigate a costly degree of computation time.

The maximum opening angle of decay can be computed by solving Einstein’s relativistic formula for the mass of the pion (0.135 MeV) and then expanded in terms of the energies and momentum of the photons that it decays into (Eq. 9.5). Since it is the parameters of the photons that are being varied to search for a best χ² fit, we can calculate the opening angle of the photon for any event by working backward and using the energy and spatial parameters of the photons present after the decay. An expansion of Equation 9.5 givesEquation 9.6, which contains a cosine term which can be expressed as Equation 9.7. This θ value is the opening angle of the decay (Eq. 9.8).

(9.5)

(9.6)

(9.7)

(9.8)

We can also look at where the opening angle will be greatest (this will happen when one photon has almost all the energy and the other has only a tiny bit⁶). The tangent of this opening angle in the lab frame will be equal to the transverse momentum of one photon divided by its momentum in the z-direction. This can be written as in Equation 9.11 (with the momentum of the photon written in the π⁰ rest frame as ).

Though it is difficult to estimate the noise background in our calorimeter, preliminary tests have shown that it is unlikely to see a photon smaller than about 300MeV. If we calculate the opening angle for this specific case, we can then use it to set a limit on the maximum separation of the two photons. To do this, we set the total momentum of the photon equal to E² which is then .09GeV (Eq. 9.12). Equation 9.12 can be solved for θ in the center of mass frame, which can then be plugged into 9.11 to solve for the opening angle in the lab frame. This calculation yields a maximum opening angle of 0.076 radians for a 10GeV π⁰, and .035 radians for a 20GeV π⁰ (corresponding to a spatial separation of 7.6 cm and 3.5 cm respectively).

(9.9)

(9.10)

(9.11)

(9.12)

It would be somewhat of a crude assumption to set the loop limits solely to span half this maximum angle radially in θ and φ outward from the energy-weighted center of the cluster because the energies are very rarely distributed evenly between the two photons. It is better to assume that the photon with greater energy hit nearer to the central tile than the smaller photon. With this assumption, the energy in surrounding tiles can be examined, and we can roughly diagnose where the lesser-energy photon hit.

The program searches tiles stacks around the central tile (in the three-by-three square) and identifies stacks with the second, third, and fourth-highest energies. The program then checks to see if any of these “high-energy” tile-stacks are adjacent to one another, and if so, sets a flag corresponding to limits on that area. For example, the top row of three tiles correspond (by geometry to the detector) to a higher φ value. The bottom row of tiles corresponds to a lower φ. The right-most column corresponds to a higher θ; the left, a lower θ. The program is restricted to only search for the best χ² within a total delta angle equal to the maximum opening angle, but the start and end points of this loop are moved around by these “high/low” flags determined by the energy distribution of surrounding tiles.

Typically, a 10 GeV π⁰event has an opening angle of around 0.013 radians, corresponding to a spatial separation of 1.3 cm on the detector surface. In contrast, a 20GeV π⁰ decay has an opening angle on the order of 0.006 radians, or about 6 millimeters spatially, and a 50GeV event has a separation given by 0.002 radians or 2mm.

The step size chosen to step through these angles is also important, and somewhat limited by the computation time constraints. Because so many energies and coordinates must be tried for each photon to get the best χ² fit, it currently takes about a week to run just 50 events.

For the results evaluated in this thesis, the angular step size is limited to 0.005 radians, and the energy step size is set to 0.05 GeV. The results do indicate that smaller steps would likely improve measurement resolution, however, this would increase the computation time to a few weeks to get even a handful of data points. Testing the program’s capabilities with smaller steps is therefore not feasible. Other means to improve resolution will be discussed in Section 9.3.

9.2 Merged-Cluster Resolution Capabilities

While this method of measuring photons in merged clusters appears promising, the resultant program resolution is not yet ideal for implementation. The program is still much too slow for necessary mass-scale testing, and in a typical case the spatial error margins may range to half a centimeter or more. A total of 50 events were used for this analysis, for a spread of energies⁷ including 10, 20, 25GeV. These events occur at a variety of incident tile points - including points from the 25-point database, and interpolated points between. Typical deviations between the program output and actual events values are summarized in Table 1. Spatial deviations are reported in radians, but can be converted to millimeters by multiplying by 1000 (the electromagnetic calorimeter radius is one meter).

The results summarized were created by setting step sizes as specified in Section 9.1.3. The larger photon contributes more energy to the cluster, so a χ² fit based on positions and energies is going to be best when there is less error in the position of the “big” photon. This is why errors reported for the larger photon’s parameters are smaller in magnitude than errors reported for the smaller photon’s parameters.

If this program were finding the “most correct” point each time (as close as we can get with the specified step size), then the average error should be concentrated within half of this radial step size, or to 0.0025 radians. The spread of frequency distributions in error for theta values (Fig. 26) for the large photon shows that 23 of 50 events had errors less than this, while 40 events had errors less than one full step size. Similarly, 20 of 50 errors in phi for the larger photon are contained within half a step size, and 41 of 50 are within one step size. This suggests that spatial resolution is heavily limited by the step size.

A few outliers removed from the tabulated data include decays where the opening angle observed was very near or at the maximum. In these decays, one of the two photons gets most all of the π⁰ energy (the second photon has an energy less than 1GeV) and the program fails to reconstruct an array within the typical error range (error jumps from a few millimeters to several centimeters). In the future, this problem may be combated by limiting the energy distribution and requiring one photon to take on almost all the π⁰ energy when the χ² calculated is for spatial parameters that yield large opening angles.

The χ² value is much more sensitive to spatial accuracy than energy resolution. Unlike the spatial errors discussed above, the errors in energy distribution can not really be pinned on the limiting step size of .05 GeV (Fig. 27). Deviations in energy for each of the two photons occupy a wide spread of values, with an average deviation of just less than 2GeV from the correct energy value. It is my hope that improved spatial resolution will also improve this energy resolution, but there is currently no evidence to suggest this to be the case.

Also, it should be noted that this program does not yet work for higher energies. A 50GeV photon has a typical opening angle of 0.002 radians, or 2mm. The step size is currently set to 0.005 radians, and a smaller step-size is not ideal in terms of computation time. At this energy, both photons will lie in this distant between one step, and the program will return the exact same θ and φ value for each of the two photons. The photons are not actually overlapping, but this point is the closest point that was “tried” to the real photon coordinates. If the resolution could be improved, the photons in high-energy merged clusters could be distinguished from one another.

9.3 Future Work

The first stages of this project have set a solid foundation for future research in the area of merged-cluster photon measurements. Code has been written to look up, rotate, and stack photon arrays, as well as compute a χ² by comparing a guess energy array to an observed energy array. However, the spatial margin of errors needs to be reduced so that the pion momentum can be more accurately reconstructed, and computation time needs to be diminished accordingly. Additionally, modifications to the nature of the utilized “step size” are needed so that higher-energy photons can be measured. There are quite a few ways one could proceed in this. In the following sections, I will highlight a few of the options that I have considered.

9.3.1 Variable Step Size

We can see from the frequency distributions in section 9.2 that maximum spatial errors linger very closely to the asserted angular step size of .005 radians (5mm on the detector surface). Precarious computation time makes it insensible to reduce this step further, but it is possible that a variable step size could be set.

A variable step could be implemented so that as the χ² measurement worsened/improved, the spatial step would accordingly grow and shrink. This would avoid wasting time calculating χ² values when the input parameters are far from the actual parameters. For this to be done successfully, it would be best to make a broad study of the χ² surface, because the χ² surface is very complex. I have found that these χ² measurements are much more sensitive to the positional measurements than to the energies of the “guess” photons specified. For example, when the correct parameters are input, it would be reasonable to observe a χ² of 131, while the same spatial parameters with a vastly wrong energy (10GeV or more) might deliver a χ² of 135. However, when the correct energy is input and the angular parameters are rounded to the hundredths place, the χ² might shoot up to well over 2000.

It would also be a good idea to limit the energy that is “guessed” based on the opening angle assumed. If the two “guess” points are far away from one another, then it is much more likely that one of the two photons will have almost all of the pion energy (this is just by nature of the Lorentz transformation which provides a relativistic “boost” proportional to the photon energies). I have made attempts, but have so far been unsuccessful in setting a self-varying step size that reduces time and does not skip over the correct point.

9.3.2 An Optimized Constrained-Fitting Algorithm

One rather obvious solution to the time-computation problem would be to implement some sort of fitting algorithm to quickly find the absolute minimum of the χ² function.

Originally, I had intended to implement an unconstrained optimized-fitting JAVA program called “UNCMIN” that would very quickly vary the six parameters (energy1, energy2, θ₁,θ₂,φ₁,φ₂) to search for a minimum χ². I chose UNCMIN because the SUSY team works primarily in JAVA, and it is the most widely-available optimized-fitting JAVA routine. UNCMIN uses use a quasi-Newtonian algorithm that is designed for multi-dimensional fitting. (See Appendix A for a discussion of the quasi-Newtonian type algorithm). However, UNCMIN is an unconstrained fitting routine, meaning that it does not permit constraints on the variables (i.e.- one cannot tell the program that maximum energy distributed between the two photons is equal to the energy of the cluster). UNCMIN would search for a minimum χ² by escalating the two energies to infinity.

One solution around this is to let the input parameters be various angles in the pion’s frame of reference (through a series of rotations and a Lorentz transformation matrix, the six photon parameters in the lab frame can then be calculated). The advantage to this is that these angles can be infinite in either direction and they needn’t be constrained. Logistics of these alternative parameters are discussed in Appendix B. However, this problem would be less complicated if a fitting routine that allows constrained fits on the variables were implemented.

There are several different types of algorithms that may be suitable for this problem. One method that will not work is a grid-search algorithm, which works to minimize the χ² by varying only one parameter at a time. I have investigated how the χ² varies by independently varying one energy, θ, or φ value. If a correct parameter point is input and the other five parameters are incorrect, the χ² is not characteristically smaller than the χ² calculated when all six parameters are incorrect. This is because the χ² surface is very intricate and contains several sharp minimums.

A gradient-search algorithm (such as UNCMIN) would be promising if constraints could be placed on the variables. In a gradient-search algorithm, all parameters are incremented simultaneously, with relative magnitudes adjusted so that the resultant direction of travel in parameter space is along the gradient (direction of greatest change in the χ² surface). This type of algorithm is ideal for approaching the minimum from far away, but it does not converge rapidly (with the tiny steps and intricate precision needed) near the minimum. A Marquardt algorithm uses a gradient-search method for the beginning of the search and then implements an analytical search method close to the minimum. This algorithm is known for being relatively insensitive to the starting parameters, and for finding fits more quickly and efficiently than other methods.

The IDL programming language supports a Marquardt algorithm called MPFIT. I have also come across similar algorithms in PYTHON and Fortran. Since the linear collider simulation code used by the SUSY group is all written in JAVA, it may be best to write a wrapper in one of these other languages to do the actual curve fitting. This wrapper could call the code that I have already written to look up energy distributions and set constraints, and then the fitter algorithm could quickly vary the parameters within such bounds to find the best fit.

10 Final Comments

I am extremely grateful to have had the opportunity to step in on this project at such a pivotal time. My efforts would not have been possible without the simulation framework, database construction, and various coding contributions from past SUSY researchers including Jack Gill, and Maria Gulda. I offer extended gratitude to Joseph Proulx, who fostered a myriad of ingenious ideas for quite some time before tendering them to my dominion. Also, I would like to thank my adviser, Uriel Nauenberg, who - on faith alone - indulged me with the opportunity to manifest such ideas. Finally, I offer sincere thanks to any future individual who may take this project further.

A Appendix: The Quasi-Newtonian Method

The Quasi-Newtonian method is a variation on Newton’s method, (also called the Newton-Raphson method), which is a root-finding algorithm that uses a Taylor series approximation to find accurate roots of a system. A first-order Taylor expansion (Eq. A.1) can be used to estimate the offset needed to land closer to the root of a function, starting from an initial guess x_o.

(A.1)

If we set the left side of eq A.1 to zero and solve for the offset ε = ε_o, we arrive at Equation A.2. If the initial guess is chosen close enough to the root, then the series can be approximated by Equation A.3.

(A.2)

(A.3)

This root-search can be generalized to several dimensions by replacing the derivative f’(x) with the gradient of the function, and the second derivative can be replaced by the inverse of the Hessian matrix Hf(x) (Fig. A). Calculating the Hessian matrix can be time costly, so the quasi-Newtonian method finds a means around calculating this matrix at a single point by gradually building up a matrix approximation using gradient information from the previous iterations [13]. If a starting value of x_o is chosen close enough to the root, then convergence is guaranteed.

B Appendix: Alternative Parameters for Un-Constrained Optimized Fitting

If unconstrained parameters must be used by a fitting algorithm such as UNCMIN, one alternative is to make all parameters be angles so that any output returned will offer a significant solution. (An infinite angle can be modded by 2π). One approach I tried (that one may wish to explore further in the future) was to let the fitting algorithm have 4 angular input parameters: θ and φ of the pion in its rest frame, and a θ and φ decay angle of one of the two photons.

With these parameters in the pion’s rest frame, the photon parameters in the lab frame can be calculated. If we let the energy of each photon be equal to the pion mass (.135 MeV)/2, then the momentum components of one of the two resultant photons in the pion’s frame can be given by Equations B.2 through B.4 below. Since the photons will split at angles of 180 degrees from one another in the pion’s frame, the momentum components for the second photon are just Px_γ2= -Px_γ1, Py_γ2 = -Px_γ1, and Py_γ2 = -Py_γ1.

(B.1)

(B.2)

(B.3)

(B.4)

We must perform a standardized rotation on this vector to rotate it into the Z-direction so that a Lorentz transformation can be performed which is independent of the initial momentum vector. Two separate Lorentz transformations of this kind will yield the momentum vectors for each of the photons in the lab frame. This calculation looks like B.5 where R_φ and R_θ are the rotation matrices about those respective coordinates (θ and φ are coordinates of the incident pion in its rest frame), L is the Lorentz transformation matrix, and R_φ^-1 and R_θ^-1 are the inverses of the (equivalent to the matrix transposes).

(B.5)

Expanded, we get the somewhat lengthy quintet of matrices to multiply by the momentum vector below:

=

Once this calculation is worked out, the energy of the first photon is just E′_γ1, and the momentum vector of the photon is known in the lab frame. Therefore, θ_γ1 = Tan^-1(), and φ_γ1 = Tan^-1(). This calculation can be repeated for the second photon with its momentum components in the pion frame just equal to the negative momentum components of the first photon.

The unconstrained fitting algorithm can then vary the parameters θ_π,φ_π,and θ_γ1,φ_γ1 to look for a minimum χ², and the lab measurements of the photons can be calculated no matter how big or negative, the input angles become.

There were a few reasons I decided not to use these parameters. First of all, the fitter still requires an “initial guess” of the best parameters and this guess must be close enough to the correct value that a false minimum is not returned. The best means of “guessing” I devised included looping through energies and positions with a relatively large step size and making χ² calculations. I wanted to then feed this rough guess for the best parameters into the fitter (UNCMIN) and let the routine do the fine-tuning on the point, but this introduced a further problem of then translating the six parameters backward to angles in the pion’s frame of reference that were now the required input.

Unfortunately, the backward calculation proved more difficult than the forward calculation, because a mass term is needed to constrain the Lorentz transformation (M= with M being the known mass of the π⁰ , 0.135 MeV). The squared cluster energy minus the squared momentum vector (the “best initial guess” for the photon parameters) will not necessarily force M to go to 0.135, meaning that yet another χ² minimization loop is required to optimize the photon momentum components to make this true. The prospect of this χ² mass minimization appeared to be so computation time costly that I chose to look for an alternate solution.

Also, at this point in my work, testing began to suggest that the program works relatively well at low energies without the implementation of this fitting routine at all. The above frame translation may be useful if one finds a better way of arriving at an “initial guess” for the best parameters that can be expressed in terms of angles in the pion’s reference frame.

References