{'fraction_non_alphanumeric': 0.03728619480426537, 'fraction_numerical': 0.011259810587084616, 'mean_word_length': 4.624580705756071, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
STARlight: A Monte Carlo simulation program for ultra-peripheral collisions of relativistic ions

Ultra-peripheral collisions (UPCs) have been a significant source of study at RHIC and the LHC. In these collisions, the two colliding nuclei interact electromagnetically, via two-photon or photonuclear interactions, but not hadronically; they effectively miss each other. Photonuclear interactions produce vector meson states or more general photonuclear final states, while two-photon interactions can produce lepton or meson pairs, or single mesons. In these interactions, the collision geometry plays a major role. We present a program, STARlight, that calculates the cross-sections for a variety of UPC final states and also creates, via Monte Carlo simulation, events for use in determining detector efficiency.

ton fluxes in transverse impact parameter space, subject to the requirement (which is also impact parameter dependent) that the colliding nuclei do not interact hadronically. The program is a two step process. First, it calculates the cross-sections for the reaction of interest, as a function of W (photon-Pomeron or two-photon center of mass energy), Y (final state rapidity) and p T (final state transverse momentum). Second, STARlight generates Monte Carlo events which can be used to determine cross-sections within specific kinematic constraints or for studies of detector efficiencies. The second step includes the decay of any unstable particles produced in the reaction, with appropriate consideration of particle spins and parity. It outputs these events in ASCII format. Solution method: The program generates a two dimensional look-up table of the production cross-section as a function of final state rapidity and mass. The dimen-sions of the table are selectable, allowing the user to choose the desired accuracy. For certain final states, a second two-dimensional look-up table, giving the transverse momentum distribution, as a function of rapidity, is also used. With these look-up tables, the program generates final states. Particle decays and the final angular distributions are calculated for each event.
Restrictions: The program is focused on ultrarelativistic collisions at Brookhaven's RHIC (Relativistic Heavy Ion Collider) and CERN's LHC (Large Hadron Collider), with final states that are visible in a central detector. At lower energies (i.e. at the CERN SPS), caution should be exercised because STARlight does not account for the longitudinal momentum transfer to the nucleus; this is larger at low beam energies. Running time: The running time depends on the binning used for the look-up tables and the number of events requested. Table generation is typically less than 15 minutes on a single core, and event generation is quick -generating 100,000 events takes only a few minutes. Running the provided sample input file to produce 1,000 rho mesons at RHIC energies takes ∼7 minutes; running the provided sample input file to produce 1,000 electron-positron pairs at RHIC energies takes ∼2 minutes. References: http://starlight.hepforge.org and references in this article.

Introduction
Ultra-peripheral collisions (UPCs) [1] are an active field of study at the Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Laboratory [2,3,4,5,6,7], at the Fermilab Tevatron [8] and at the Large Hadron Collider (LHC) [9,10,11,12,13,14,15,16]. In UPCs, the two nuclei interact electromagnetically rather than hadronically by physically missing each other. Two-photon and photonuclear interactions can occur. STARlight is capable of simulating a range of two-photon and photonuclear interactions: two-photon production of single mesons and photoproduction of a variety of vector mesons. It can also interface with the DPMJET III event generator to generate general photonuclear interactions. Short-lived mesons (including vector mesons) are usually decayed within STARlight and take into account the spin state. For more complex cases, such as mesons with multiple final states, PYTHIA is used for the decays, but the spin information is lost.
In ultra-peripheral collisions, two types of interactions can be distinguished: purely electromagnetic, usu-ally referred to as two-photon interactions, or photonuclear.
In a two-photon interaction, the two electromagnetic fields, carried by the two nuclei, collide. These fields couple to any particles that carry charge, so a wide variety of final states are possible. For relativistic nuclei, the electric and magnetic fields are perpendicular, and the configuration may be represented as a flux of almost-real photons, following the Weizsäcker-Williams method. At an impact parameter b from a relativistic nucleus with charge Z and Lorentz boost γ, the photon number density is [1] where k is the photon energy, x = kb/γ c, α ≈ 1/137 is the electromagnetic fine structure constant and K 1 and K 0 are modified Bessel functions. The output N is the number of photons per unit area, per unit energy. The first term gives the flux of transversely polarized photons; it dominates for relativistic nuclei, and it is the only one considered further here. As long as the ions are moderately relativistic, the second term is unimportant. Through these two-photon interactions, pairs of leptons, quarks or electroweak bosons may be produced. If a quark pair is produced, they hadronize into a single meson (with the appropriate spin/parity), a meson pair or a more complex final state. The cross-section depends on the overlap of the two photon fields and on the two-photon coupling to the final state. STARlight can generate a number of single-meson and lepton pair final states.
In photonuclear interactions, the most common exclusive final state is a vector meson. In the vector meson dominance picture, photons fluctuate to qq pairs, which scatter elastically from the target nucleus. For heavy nuclear targets, the amplitudes for scattering from individual nucleons add coherently when the 4−momentum transfer q is small (|q| < /R A ), leading to a large crosssection. The treatment of coherence is handled via a Glauber formalism which accounts for multiple interactions within the nucleus [17].
STARlight also simulates photon scattering at higher p T , where the amplitudes from the multiple target nucleons add incoherently.
In photonuclear UPCs, either nucleus can emit the photon, with the other nucleus serving as a target. These possibilities are indistinguishable. Exchanging the photon emitter and target introduces a propagator exp (−i k · b) to account for the position difference. Here k and b are the photon and impact parameter vectors [18]. For symmetric collisions (e.g. proton-proton or gold-gold), swapping the nucleus is equivalent to a parity transformation. Vector mesons are negative parity, so the amplitudes subtract and, for p T → 0, vector meson production disappears. For proton-antiproton collisions, swapping the nuclei is a charge-parity transformation; since vector mesons are CP positive, the two amplitudes add [19].
One important aspect of ultra-peripheral collisions is that, at least for heavy-ions, multiple interactions are possible between a single ion pair. This has important consequences; it allows for studies involving correlations among the multiple final state particles. For example, one can study the azimuthal angular correlations between the neutrons emitted in mutual Coulomb dissociation [20] or rapidity correlations for states like ρ 0 ρ 0 [21].
Experimentally, UPCs accompanied by single or multiple Coulomb excitation have proven to be very important. The nuclear decays produce neutrons which are easily visible in zero degree calorimeters, greatly simplifying triggering. Many experimental results to-date have relied on these neutrons for triggering. Because of the kinematics of the photon emission and absorption [20], to a large degree, each photon acts independently, and the cross-section for multiple reactions is given by where the P i (b) are the probabilities for the different subreactions, as a function of impact parameter. These probabilities may be substantial; a unitarization procedure corrects for the possibilities of multiple subreactions, such as the production of two e + e − pairs. STARlight can calculate cross-sections and generate events for two-photon and photonuclear interactions accompanied by mutual Coulomb excitation. These are three or four photon interactions: one or two for the photonuclear or two-photon interaction, plus one to excite each nucleus. It can also calculate the cross-sections for reactions accompanied by mutual excitation to a Giant Dipole Resonance (GDR) interactions; the GDR excitations usually decay by single neutron emission. STARlight can simulate the collision of two dissimilar nuclei; both Z and A are selectable by the user. Heavy nuclei are modeled either as hard spheres, or following the Woods-Saxon mass distribution. For selected nuclei (gold, lead, copper), STARlight uses measured parameters for the nuclear radius [22]. Otherwise, the Woods-Saxon radius is set to R A = 1.2 fm · A 1/3 . A skin depth of 0.53 fm is used for all nuclei. For light nuclei, with Z ≤ 6, a Gaussian mass distribution is used.
Each projectile beam energy are defined by the Lorentz factor γ set by the user. The input beam energies can be asymmetric for the two projectiles, but internally all calculations are performed in the center of mass frame. Before the output is written, the final state particles are boosted to the laboratory frame specified by the user. STARlight is optimized for RHIC and LHC energies, although other energy ranges should be usable. STARlight 2.2 should closely reproduce the calculations in the original publications [17,18,23], with the following exceptions. First, for some nuclei, the nuclear radii are set to the measured values, while previously a parameterization assuming R ∝ A 1/3 was always used; this new radius is used exclusively. This reduces the cross-section for gold beams by about 7%, but has almost no effect for lead beams. Second, the parameterizations of the γ + p → J/ψ + p and γ + p → Υ + p cross-sections have been improved, as described below.

Photon spectra and form factors
The photon spectrum is calculated in impact parameter space by integrating Eq. 1 over all impact parameters under the assumption that the beam projectiles do not interact. For photonuclear interactions, the flux is and for two-photon interactions, Here, N(k 1 , b 1 ) is the photon density for photons with energy k 1 at position b 1 from nucleus 1 and N(k 2 , b 2 ) is the corresponding density from nucleus 2. The vectors b 1 and b 2 have their origins at the center of each nucleus. So, | b 1 − b 2 | is the impact parameter between the two nuclei. The requirement that there be no hadronic interactions reduces the usable luminosity. The actual cross-section is sensitive to how this non-interaction is defined [23,24,25]. The probability of not having a hadronic interaction at impact parameter b, P NOHAD ( b), is implemented in different ways for nucleus-nucleus, proton-nucleus and proton-proton collisions.
STARlight uses a trick to lower the dimensionality of the integral in Eq. 4. The integral [24,26]. In nucleus-nucleus collisions, the probability to have no hadronic interaction at impact parameter b is where T AA ( b) is the nuclear overlap function, which is calculated from the nuclear density profiles. The nuclear densities are assumed to follow a Woods-Saxon distribution, and the nucleon-nucleon interaction crosssections follow the particle data group parameterization for pp collisions for proton-proton center of mass energies √ s above 7 GeV [27]: σ = (33.73+0.2838 ln 2 (r)+13.67r −0.412 −7.77r −0.5626 )mb (6) where r = s/1GeV 2 . The use of P NOHAD (b) reduces the γγ luminosities by up to 20% compared to the simpler requirement that b > 2R A [23], because there is still a substantial probability of two nuclei interacting hadronically at impact parameters larger than 2R A .
In proton-nucleus collisions, the hadronic interaction probability is also calculated from the Glauber model. The probability of having no hadronic interactions is where T A ( b) is the nuclear thickness function. For large nuclei, applying this probability function is roughly equivalent to requiring b > R A . This probability is used whether the photon is emitted by the nucleus or the proton. So, the effective photon flux from the proton in a proton-nucleus collision is considerably smaller than if the flux is calculated from just the proton form factor [28].
In proton-proton collisions the probability of having no hadronic interactions is calculated from where Γ(s, b) is the Fourier transform of the pp elastic scattering amplitude [29]. This is modeled by an exponential with b 0 = 19.8 GeV −2 . This interaction probability corresponds roughly to a cut in minimum impact parameter of b > 1.4 fm.
The nuclear form factor is the convolution of a hardsphere form factor with a Yukawa potential with a range of 0.7 fm. This is very close to the the Woods-Saxon distribution, and the form factor may be determined analytically [17]. For a 4-momentum transfer q, (10) Here, ρ 0 is the nuclear density.
For light nuclei, with Z ≤ 6, a Gaussian form factor is used. Proton-proton and proton-nucleus collisions are treated separately, as discussed below, with the dipole form factor for the proton [19,28]

Two-photon interactions
The two photon energies, k 1 and k 2 , are related to the final state invariant mass W and rapidity Y through and The two-photon luminosity d 2 N/dk 1 dk 2 in Eq. 4 can thus be transformed to d 2 N/dWdY. For two-photon interactions, the cross-section to produce a final state X factorizes into a two-photon cross-section and a twophoton luminosity [1]: The limits on |Y| (a symmetric region is assumed) and W in the integrals of Eq. 14 are set by the user.
In differential form, Eq. 14 becomes The look-up table produced during the first phase of the simulation is based on this equation, which thus provides the basis for the invariant mass and rapidity distributions.
The final state transverse momentum p T is the vector sum of the p T of the two photons. The p T distribution of the photons depends on the nuclear form factor [18,26,30,31] as defined by Because the p T distribution depends on the photon energy, it cannot be pre-calculated for every condition. For simplicity, STARlight determines the photon p T using rejection sampling. The azimuthal angle in the transverse plane is assumed to be uniformly distributed. Three types of two-photon interactions are implemented currently: two-photon production of lepton pairs, single meson production and decay, and the reaction γγ → ρ 0 ρ 0 , followed by ρ decays.

Lepton pairs
STARlight creates lepton pairs using the equivalent photon approximation (EPA) approach. This is a lowest order calculation, where the photons are treated as massless. The kinematic distributions are in generally good agreement with the data, except that the EPA approach finds a lower average p T than in the data. A full quantum electrodynamic calculation, including the photon virtuality, found better agreement with the data [4]. Later calculations showed that the STAR result includes room for an additional contribution to the cross-section, due to higher order corrections [32,33]. These corrections are expected, but their effects on the kinematic distributions are not well known. Still, STARlight lepton pairs are useful for most purposes.
The cross-section to produce a pair of leptons with lepton mass m and pair invariant mass W is given by the Breit-Wheeler formula [34] In the region W ≥ 2m, the cross-section drops very rapidly with increasing W, so it is necessary to use a very large number of bins (≥ 100 bins/GeVc −2 ) in W to accurately evaluate the total cross-section. The angular distribution of these lepton pairs is given by [34] (18) where θ is the angle between the beam direction and one of the leptons, in the lepton-lepton center of mass frame. Here, we neglect the effect of the photon p T on the angular distribution.
All three leptons are treated as stable; τ (and, if desired, µ) decays must be handled externally. One weakness of this approach (to be corrected in later versions) is that it loses the effects of the correlations between the lepton spins.
Examples of the kinematic distributions are shown in Figs. 1 and 2. Figure 1 shows the rapidity of each daughter from electron-positron pairs produced in Xn − Xn Au+Au collisions at √ s NN = 200 GeV. Figure   2 shows the transverse momentum from µ − µ + produced in Pb+Pb collisions at LHC energies.

Mesonic final states
The cross-section to produce a single meson from two photons depends on its two-photon width Γ γγ , mass M R , total width Γ and spin J [35] where Eq. 20 applies only for narrow resonances. STARlight can simulate either narrow or wide resonances; in the latter case, the user can select the maximum and minimum W. STARlight simulates the production of a variety of spin 0 and spin 2 mesons, using standard values of Γ γγ .
This distribution is sampled using rejection sampling. The f 2 ′ (1525) is a special case, in that 50% of the time it decays to K s K s . These K s 's are written to the output file, and must be decayed by an external program, such as GEANT.
For some two-photon final states, STARlight uses PYTHIA (8.2 or higher) to handle particle decays. These are the η, η ′ , a 2 (1320) and the η c . The advantage of using PYTHIA is that it can include multiple branching ratios and more complicated final states; the disadvantage is that the spin directional information is lost.
STARlight can simulate the reaction γγ → ρ 0 ρ 0 at threshold (i.e., both ρ 0 are produced with zero momentum in their center of mass frame) with σ ρ 0 ρ 0 = 100 nb. This mode is included to simulate the excess which has been observed in two-photon production, which has been hypothesized to stem from a 4-quark state [36].

Overview of photonuclear interactions
STARlight simulates the photoproduction and decay of a variety of vector meson final states: ρ, ω, φ, J/ψ, ψ ′ , Υ(1S ), Υ(2S ) and the Υ(3S ). It can also simulate ρ+direct π + π − photoproduction, including interference between the two channels. The cross-sections are mostly based on parameterized HERA data for γp → V p. It can also simulate general photonuclear interactions using DPMJET III.
The ρ 0 and ω cross-sections, σ(γ + p → V + p), are parameterized as a function of the γp center of mass, W γp : The first term is for Pomeron exchange and the second for meson exchange. The other mesons are produced solely via Pomeron exchange, so σ M = 0. For the J/ψ, ψ ′ and Υ states, the power law is supplemented with a factor that accounts for the near-threshold decrease in the cross-section: This is shown in Fig. 3a. The cross-sections for ψ(2S ) and Υ(2S , 3S ) states are assumed to have the same shape, with m V replaced by their respective masses. The values of the parameters were obtained from fitting experimental data and are given in Table 1. The rapidity distributions for Υ(1S ), Υ(2S ) and Υ(3S ) produced in pp collisions at √ s = 13 TeV are shown in Fig. 3b. For most of the particles, the masses and widths are the standard Particle Data Group values. The exception is the more complex ρ ′ , with considerable evidence for multiple overlapping resonances [37]. Because key details about individual states and the interference are poorly known, STARlight uses a single resonance with a mass of 1540 MeV and a width of 570 MeV, chosen to match STAR observations [2]. The ρ ′ -photon coupling and nuclear interaction cross-section are arbitrarily taken to be the same as for the ρ.
The vector mesons are decayed assuming that the photon polarization is parallel with the beam axis.
The same nuclear breakup options as for the twophoton interactions are available.

Coherent vector meson production
The cross-sections for coherent production on nuclear targets are determined using a classical Glauber calcu- Table 1: Vector meson parameters for W γp , in GeV. The ρ, ω and φ parameters are from Ref. [38], while the J/ψ values are obtained from a fit of data from [39,40] to Eq. 23. The ρ ′ characteristics are not well known, so it is arbitrarily given the same parameters as the ρ. The Υ(1S ) cross-section uses the parameterization in Eq. 23 and Fig. 3a, while the Υ(2S ) and Υ(3S ) cross-sections are scaled from this, based on their couplings to e + e − (i.e. σ ∝ 1/ f 2 V ), per Eq. 10 of Ref. [17]. The ψ(2S ) cross-section is taken to be 0.166 that of the J/ψ, following Ref. [41]. Previous versions of STARlight used older, slightly different values for the J/ψ [17].
where F(t) is the nuclear form factor. The latter formula is useful because it shows how the p T distribution of the vector mesons is determined. The photon flux dN γ (k)/dk is given by integrating Eq. (1) over impact parameter space, subject to the condition that the two nuclei do not interact hadronically. For each impact parameter, the photon flux striking the nucleus is determined and spread evenly over the nucleus.
The calculation of the cross-section can be done assuming narrow resonances, in which case t min = (M 2 V /4kγ) 2 , or it can be done by convoluting the photon spectrum with the Breit-Wigner shape [17]. This is decided by the user-selected parameter PROD MODE. Setting PROD MODE=2 gives the narrow resonance approximation and setting PROD MODE=3 uses a wide resonance. The difference is significant only for the ρ 0 , where the cross-section is reduced by about 5% in heavy-ion collisions when a wide resonance is used.
The setting of PROD MODE only affects the calculation of the cross section. Regardless of the value of PROD MODE, the event generation is always done with a Breit-Wigner invariant mass distribution.
STARlight calculates a two-dimensional look-up table, covering W and Y, and then generates W, Y pairs for each event by sampling from the table. The W and Y range and number of bins are settable, or STARlight can pick a W range to cover the vector meson in question.
There are multiple options for generating the vector meson p T spectra. STARlight can generate the p T spectra assuming that photoproduction on the two nuclei is independent as shown in Fig. 4 for J/ψ produced in Pb+Pb collisions at LHC energies, or including the interference between the two production sites [18]; this is controlled by the user via the IN-TERFERENCE PARAMETER. When interference is turned on, it can also generate spectra including interference, but with an altered strength, set by the input parameter IF STRENGTH. To simulate proton antiproton collisions, where the interference has the opposite sign to proton-proton collisions [19], one sets IF STRENGTH= −1. When interference is used, the p T spectrum is calculated in another two-dimensional look-up table, in Y and p T ; the range and number of p T bins in this table are settable via the input parameter INT PT N BINS; the number of Y bins is the same as for the cross-section table.

Incoherent vector meson production
The incoherent photonuclear cross-section is calculated under the assumption that it scales the same way as the total inelastic vector meson nucleus cross-section. This can be written Assuming vector meson dominance and using the classical Glauber expression for the total inelastic crosssection, this can also be expressed as (26) The term inside the bracket corresponds to the probability of having at least one vector meson-nucleon interaction at impact parameter b. It is thus assumed in the model that one or more vector meson-nucleon interactions lead to the emergence of a real vector meson and that it is not destroyed by multiple interactions. Other authors have assumed that incoherent production occurs only if the vector meson reacts with exactly one nucleon [44] or calculated the incoherent cross-section based on event-by-event variations in the positions of the nucleons (i.e the variations in the second moment of the amplitudes) [45]. Other authors have made a direct connection between incoherent interactions and nuclear breakup, leading to neutron emission [46]. These approaches may lead to a reduction in the incoherent cross-section at small p T ; this is not present in STARlight. Figure 3d shows the Υ(1S ), Υ(2S ) and Υ(3S ) production in incoherent Pb+Pb collisions at √ s NN = 5.1 GeV.
The transverse momentum is calculated in the same way as for coherent production, but with the nuclear form factor replaced with a nucleon form factor. For the light vector mesons, ρ 0 and ω, a dipole form factor is used: This corresponds to a Fourier transform of an exponential matter distribution. The value of Q 0 depends on the width of this distribution. The width of the matter distribution is chosen such that the RMS radius is equal to the squared sum of the proton and pion RMS radii, leading to Q 2 0 = 0.45 GeV 2 . For the heavier vector mesons, a narrower transverse momentum distribution has been observed in γ + p → V + p data and for these an exponential is used for the form factor [38]: with b = 7.0 GeV −2 for φ, b = 4.0 GeV −2 for J/ψ and Υ, and b = 4.3 GeV −2 for ψ(2S). The distribution used for the light vector mesons is slightly wider than what one gets from an exponential distribution in Q 2 with b ≈ 10 GeV −2 , as has been found previously [47], but it gives a good description of data on ρ 0 production from ALICE [12].

4.3.
π + π − π + π − final states The structure of the 4π final state is somewhat complex in that the bulk of the cross-section likely comes from the decay of two or more excited ρ states. In the absence of a detailed model, the 4π final state is decayed following a phase space distribution. Because the final state depends on the ρ ′ −nucleon cross-section σ ρ−N , if there are two different vector mesons (with different σ ρ−N ), then the actual 4π mass spectrum and substructure could depend on the ion species being collided. The invariant mass distribution of 4π produced in Au+Au collisions at √ s NN = 200 GeV is shown in Fig. 5.

Interference and p T spectra
The vector meson p T is the vector sum of the photon and Pomeron p T . The photon p T is as discussed in the section on two-photon interactions. The Pomeron p T is determined by the form factor of the target nucleus, following the distribution in Eq. 24 with the form factor in Eq. 10. The production amplitude is evenly distributed throughout the nucleus (i.e. absorption is neglected).
Interference between vector meson production on the two nuclear targets is implemented via the cross-section formula: where σ(y, p T ) gives the p T spectrum without interference. Here, c depends on the sign of the combination: c=1 for ion-ion collisions, and c = −1 for proton-antiproton collisions (e. g. at the Tevatron). In STARlight, c is user settable (input parameter IF STRENGTH); it can be any real number. In this approach, Eq. 29 assumes that the production phase is independent of energy. This should be satisfied for Pomeron exchange, but may fail at low energies where production proceeds via meson exchange.
For incoherent photoproduction, this interference is not important, and is therefore not implemented.

General photonuclear interactions
STARlight is interfaced with the DPMJET III [48] Monte Carlo to simulate general photonuclear interactions γ + A → X. This includes photonuclear jet production and production of soft hadrons. Photon emission from one or both nuclei can be handled, and the photon spectrum is calculated for the two cases as described in [49]. Single photon exchange events are characterized by a rapidity gap on the side of the photon-emitting nucleus. The minimum and maximum photon energies (in the rest frame of the target nucleus) must be set by the user through the input parameters MIN GAMMA ENERGY and MAX GAMMA ENERGY (in GeV). DPMJET requires that the photon energy is greater than 6 GeV so MIN GAMMA ENERGY should be set higher than this value. A minimum photon energy of 6 GeV is sufficient to simulate particle production around midrapidity (|η| ≤ 4.5 in Pb+Pb interaction at the LHC), but particle production at very forward rapidities, characterized by low photon energies, cannot be completely described.

Photonuclear breakup
Experimentally, it is usually easier to study ultraperipheral collisions when they are accompanied by photonuclear breakup. The breakup reaction produces neutrons near the beam rapidity; these neutrons are easily detected in forward calorimeters [2,3,4].
Multiple interactions proceed via multiple independent photons [20] and are incorporated following Eq. 2. Since nuclear breakup depends only on the nucleusnucleus separation, it can be included without adding any new integrals. Instead, the two-photon luminosity and the photonuclear cross-section are both multiplied by a probability factor, P X (b), giving the probability of the desired nuclear breakup at that impact parameter.
STARlight considers two types of nuclear breakup: Xn, which is a breakup leading to any number of neutrons, and Giant Dipole Resonance (GDR) excitation, which usually leads to a single neutron emission, 1n [50]. One may require that one or both nuclei break up and the selection is handled by the input parameter BREAKUP MODE.
The cross-section for a nucleus to be excited is [51,52] The subscript '1' shows that this is a non-unitarized probability, and for high enough energy collisions, P 1 (b) can exceed one. The actual probability is given by a unitarization procedure. The excitation crosssection σ γA→A * (k) was determined using experimental data [53]. Because of this, the excitation code works only for gold or lead beams. For proton-proton or proton-nucleus interactions, only the modes that place no restriction on Coulomb excitation are implemented. An interaction with photon energy less than 30 MeV/c (in the target rest frame) is assumed to lead to a Giant-Dipole resonance and, consequently, to a single emitted neutron. Higher energy photons are assumed to produce higher excitations, with multiple neutrons.

Program overview
STARlight runs in two phases. The first phase calculates cross-sections and other kinematic distributions, and stores them in a series of look-up tables. The second phase uses these look-up tables to rapidly simulate nuclear interactions. Because of this two-stage approach, it can generate large numbers of events very rapidly.

Look-up table structure
The first look-up table to be generated is a two-dimensional table for event mass and rapidity, d 2 N/dYdW. For photonuclear vector meson production, the shape of the mass distribution is a relativistic Breit-Wigner distribution [54] with mass and width taken from the PDG. For two-photon production of lepton pairs, the shape of the mass distribution is calculated from the Breit-Wheeler formula (Eq. 18) convoluted with the photon fluxes from each nucleus. For two-photon production of single mesons, the mass distribution is a delta function at the PDG mass.
Additionally, for photonuclear interactions, a second table stores the p T distribution, in the form of d 2 N/d p T dY. In actuality, the p T depends slightly on W, but neglecting this small correlation does not adversely affect the accuracy of STARlight. For γγ interactions, the p T distribution is generated semi-analytically, and no table is needed.
The range of these tables in Y and W, and the number of bins for each dimension, are all user settable. STARlight stores integral tables that are normalized to 1.0. This simplifies the look-up -for each variable, STARlight selects a random number, finds the appropriate locations in the table that straddle that number, and then linearly interpolates between the two nearest values. For two-dimensional tables, an additional onedimensional array holds a collapsed version of the table, just containing dσ/dY. Using this, a value of Y is determined. A second random number is chosen, and the value of W or p T is found by looking at the two rows in the two-dimensional tables that straddle the chosen value of Y. Then, linear interpolation is performed on the two rows, to determine two values of W corresponding to the Y values for those rows. Finally, another linear interpolation is performed, in Y, to determine the final W. Or course, the same value of Y is used to determine W and p T .

Overview of the program structure
STARlight is coded in C++. The package contains files in several folders: Src contains the C++ code Include contains the header files that are needed by the C++ code Config contains example input files Cmake modules contains cmake macros Utils contains root macros that can be used to convert the output to a root tree, and to produce a simple analysis There is also a collection of top-level files outside the folder structure including • CMakeLists, a build file • Makefile created by CMake

• Readme contains instructions for installing
STARlight, a brief description of the physics incorporated in the code, and a description of input parameters and supported channels • starlightDoxyfile.conf gives the settings used by the doxygen documentation system. doxygen documentation may be created (in the doxygen directory) via "doxygen starlightDoxyfile.conf" The driver program Main.cpp instantiates, initializes and runs starlightStandalone.cpp. starlightStandalone.cpp in turns calls methods from STARlight to initialize and produce events of the selected channel and decay them as needed. The final daughter particles are boosted to the lab frame and written out in ASCII format to the file slight.out.
The init method in the starlight class reads the input parameters from the file slight.in and uses these to determine the beam-beam system and the interaction type (e.g., photon-nucleus vs. two-photon, coherent vs. incoherent). It checks whether appropriate differential luminosity tables (dN/dWdY, with granularity set by input parameters specifying the number of bins in W and Y) exist, and calls the appropriate method to create them if they do not. The differential luminosity tables are stored in the file slight.txt. It instantiates the appropriate event channel and calculates the cross-section. Important classes used in the init method of STARlight: • Nucleus: includes methods to determine nuclear radius, thickness and form factor.
• Beam: derives from nucleus; includes methods to calculate photon flux, based on impact parameter and photon energy.
• Beambeamsystem: includes methods to boost to CMS, calculate probability of breakup.
• photonNucleusCrossSection: base class with methods to calculate cross-section for coherent photon-nucleus interactions.
• narrowResonanceCrossSection: derives from pho-tonNucleusCrossSection, includes a method to find the total cross-section for a narrow resonance.
• wideResonanceCrossSection: derives from pho-tonNucleusCrossSection, includes a method to find the total cross-section for a wide resonance.
• incoherentVMCrossSection: derives from pho-tonNucleusCrossSection, includes a method to find the total cross-section for incoherent photonnucleus interactions.
• twophotonluminosity.cpp: calculates the two photon luminosity table, as a function of W and Y.
The run method in starlightStandalone calls the pro-duceEvent method in STARlight. The produceEvent method contains an iterative loop to generate the requested number of events using the method appropriate to the event channel. All methods derive from the eventChannel class. Each method randomly picks an energy and rapidity within the selected range, produces the parent system, and decays it where appropriate. starlightStandalone then boosts the final state particles into the lab frame and writes them into the ASCII file slight.out Important classes used in the produceEvent method of STARlight: • eventChannel: derives from readLuminosity; base class with methods to produce events and keep track of the number attempted and accepted.
• Gammaavectormeson: derives from eventChannel, includes methods to produce vector mesons from photon-nucleus interactions.
• Gammagammaleptonpair: derives from eventChannel, includes methods to create lepton pairs from photon-photon interactions.
• Gammagammasingle: derives from eventChannel; includes methods to create single-meson events in photon-photon interactions.
The utils directory includes a few examples which may be useful for some users. ConvertStarlightAsci-iToTree.C converts STARlight output to a ROOT TTree, while Analyze.cxx reads STARlight output and produces some useful histograms. AnalyzeTree.cxx reads the output from ConvertStarlightAsciiToTree and produces similar histograms.

Description of input data
STARlight runs are controlled through the slight.in file. The input parameters are listed below, with typical values for LHC Pb+Pb running given in parentheses. • PROD MODE=1: Two-photon interaction.
• PROD MODE=2: Coherent photonuclear vector meson production assuming narrow resonances. In pA collisions, this option means that the proton emits the photon and that the gamma-A interaction is coherent.
• PROD MODE=3: Coherent photonuclear vector meson production assuming wide resonances. This option should be used for exclusive ρ 0 production.
• PROD MODE=4: Incoherent photonuclear vector meson production. In pA collisions, this option means that the nucleus emits the photon.
14. N EVENTS: Number of events to be produced. • BREAKUP MODE = 1: hard sphere nuclei (no hadronic break-up if impact parameter is greater than the sum of nuclear radii, no restriction on Coulomb break-up).
• BREAKUP MODE =2: requires Coulomb break-up of both nuclei, with no restriction on the number of neutrons emitted by either nucleus (XnXn).
• BREAKUP MODE = 3: requires Coulomb break-up of both nuclei, but requires that a single neutron is emitted from each nucleus (1n1n).
• BREAKUP MODE = 4: requires Coulomb break-up of neither nucleus. (0n0n) • BREAKUP MODE = 5: requires that there be no hadronic break up, no restriction on Coulomb break-up (This is similar to option 1, but with the actual hadronic interaction probability).
• BREAKUP MODE = 6: requires Coulomb break up of one or both nuclei, with no restriction on the number of neutrons emitted (XnXn + 0nXn + Xn0n).
• BREAKUP MODE = 7: requires Coulomb break up of only one nucleus, with no restriction on the number of neutrons emitted (0nXn+ Xn0n). The following parameters are used only when interfacing with the PYTHIA and/or DPMJET interfaces: 30. MIN GAMMA ENERGY: Allows the user to set the minimum photon energy (in GeV) in the rest frame of the target nucleus. The default is 6.0 GeV and it should never be set below this value since DPMJET was not designed to handle low energy interactions. 31. MAX GAMMA ENERGY: Allows the user to set the maximum photon energy (in GeV) in the rest frame of the target nucleus. The default is 60000.0 GeV. 32. PYTHIA PARAMS: Used to supply input parameters to the PYTHIA interface. This takes a string to pass on semi-colon separated parameters to PYTHIA 6. eg: "mstj(1)=0;paru(13)=0.1" (the default is a blank string " ") 33. PYTHIA FULL EVENT RECORD: Determines whether the full event record from PYTHIA is written to slight.out. true = yes, false = no (false). The additional information added is as follows: daughter production vertex (x [mm], y [mm], z [mm], t [mm/c]), mother1, mother2, daughter1, daughter2, PYTHIA particle status code. The PYTHIA 8 Particle Properties page describes in more detail the properties of mother, daughter, and status code designations.

Description of output
STARlight outputs an ASCII file named slight.out. For each event, a summary line is printed, with the format EVENT: n ntracks nvertices , where n is the event number (starting with 1), ntracks is the number of tracks in the event, and nvertices is the number of vertices in the event (STARlight does not currently produce events with more than one vertex). This is followed by a line describing the vertex, of the format VERTEX: x y z t nv nproc nparent ndaughters , where x, y, z and t are the 4-vector components of the vertex location, nv is the vertex number, nproc is a number intended to represent the physical process (always set to 0), nparent is the track number of the parent track (0 for primary vertex) and ndaughters is the number of daughter tracks from this vertex. This is followed by a series of lines describing each of the daughter tracks emanating from this vertex. Each track line has the format TRACK: GPID px py py nev ntr stopv PDGPID , where GPID is the Geant particle id code, px, py and pz are the three vector components of the track's momentum, nev is the event number, ntr is the number of this track within the vertex (starting with 0), stopv is the vertex number where track ends (0 if track does not terminate within the event), and PDGPID is the Monte Carlo particle ID code endorsed by the Particle Data Group.

Description of test data
The following parameters were used to produce Additional input files used to generate the figures in this paper are available in the code distribution, in the config folder.

Summary
STARlight is a versatile C++ code which can simulate a variety of important ultra-peripheral collisions between relativistic heavy nuclei. The output in terms of cross-sections, rapidity and p T distributions etc. can and has been tested against experimental results. Events generated by STARlight have also been used by several experiments to calculate corrections for acceptance and efficiency; in these calculations, it is important to have an event generator which provide reasonable phase space distributions of the produced particles. STARlight has been found to describe most important features of the data, from both RHIC and the LHC, very well.