# 1 Program Summary

## 1 Program Summary

| Title of program: | SOPHIA 2.0 |
| --- | --- |
| Catalog number: |  |
| Program obtainable | from A. Mücke |
|  | e-mail: amuecke@physics.adelaide.edu.au |
| Computer on which the program | DEC-Alpha and Intel-Pentium based |
| has been thoroughly tested: | workstations |
| Operating system: | UNIX, Linux, Open-VMS |
| Programming language used: | FORTRAN 77 |
| Memory required to execute | $`<`$1 megabyte |
| No. of bits in a word | 64 |
| Has the code been vectorized? | no |
| Number of lines | 16500 |
| in distributed program: |  |
| Other programs used in SOPHIA | Rndm |
| in modified form | Processor independent random number |
|  | generator based on Ref.  |
|  | Jetset 7.4 |
|  | Lund Monte Carlo for |
|  | jet fragmentation  |
|  | DECSIB |
|  | Sibyll routine which decays |
|  | unstable particles  |
| Nature of physical problem: | Simulation of minimum bias photo- |
|  | production for astrophysical applications |

| Method of solution: | Monte Carlo simulation of individual events |
| --- | --- |
|  | for given nucleon and photon energies, |
|  | photon energies are sampled |
|  | from various distributions. |
| Restrictions on the complexity | Incident ambient photon fields limited |
| of the problem | to power law and blackbody spectra so far. |
|  | No tests were done for center-of-mass |
|  | energies $`\sqrt{s}1000`$GeV. |
| Typical running time: | 10000 events at a center-of-mass energy |
|  | of 1.5 GeV require a typical |
|  | CPU time of about $`75`$ seconds. |

## 2 Introduction

The cosmic ray spectrum extends to extremely high energies. Giant air showers have been observed with energy exceeding $`10^{11}`$ GeV . Energy losses due to interactions with ambient photons can become important, even dominant for such energetic nucleons, above the threshold for pion production. Photoproduction of hadrons is expected to cause a distortion of the ultra-high energy cosmic ray (CR) spectrum by interactions of the nucleons with the microwave background (the Greisen-Zatsepin-Kuzmin cutoff ; see also  for additional references), but it may also be relevant to the observed high energy gamma ray emission from jets of Active Galactic Nuclei (AGN) (e.g. ) or Gamma-Ray Bursts (GRB) . Moreover, it is the major source process for the predicted fluxes of very high energy cosmic neutrinos (e.g. ).

The photohadronic cross section at low interaction energies is dominated by the $`\mathrm{\Delta }`$(1232) resonance. Since the low energy region of the cross section is emphasized in many astrophysical applications, the cross section and decay properties of the prominent $`\mathrm{\Delta }`$-resonance have often been used as an approximation for photopion production, and the subsequent production of gamma rays and neutrinos . As discussed in , this approximation is only valid for a restricted number of cases, and does not describe sufficiently well the whole energy range of photohadronic interactions. A more sophisticated photoproduction simulation code is needed to cover the center-of-mass energy range of about $`\sqrt{s}110^3`$ GeV, which is important in many astrophysical applications.

In this paper we present the newly developed Monte-Carlo event generator SOPHIA (Simulations Of Photo Hadronic Interactions in Astrophysics), which we wrote as a tool for solving problems connected to photohadronic processes in astrophysical environments. The philosophy of the development of SOPHIA has been to implement well established phenomenological models, symmetries of hadronic interactions in a way that describes correctly the available exclusive and inclusive photohadronic cross section data obtained at fixed target and collider experiments.

The paper is organized as follows. After introducing the kinematics of $`N\gamma `$-interactions (Sect. 3) we give a brief physical description of the relevant photohadronic interaction processes (Sect. 4). The implementation of these processes into the SOPHIA event generator, together with the method of cross section decomposition and parametrization, is described in Sect. 5. The structure of the program is outlined in Sec. 6, and a comparison of the model results with experimental data is provided in Sec. 7. The definitions of special functions and tables of parameters used in the cross section and final state parametrization, as well as a compilation of the important routines and functions used in the code, are given in the appendices.

Unless noted otherwise, natural units ($`\mathrm{}=c=e=1`$) are used throughout this paper, with GeV as the general unit. In this notation, cross sections will be in GeV<sup>-2</sup>. A general exception is Section 4, where numerical parametrizations of cross sections are given in $`\mu `$barn. The relevant conversion constant is $`(\mathrm{}c)^2=389.37966`$ GeV$`{}_{}{}^{2}\mu `$barn.

## 3 The physics of photohadronic interactions

### 3.1 Kinematics of $`N`$-$`\gamma `$ collisions

There are three reference frames involved in the description of an astrophysical photohadronic interaction: (i) the astrophysical lab frame (LF), (ii) the rest frame of the nucleon<sup>5</sup><sup>5</sup>5In the following the subscript $`N`$ is used if no distinction between protons and neutrons is made. Interactions for both protons and neutrons are implemented in SOPHIA. (NRF), and (iii) the center-of-mass frame (CMF) of the interaction. For example, in the LF, the initial state can be characterized by the nucleon energy $`E_N`$, the photon energy $`ϵ`$, and the interaction angle $`\theta `$

$$\mathrm{cos}\theta =(\stackrel{}{p}_N\stackrel{}{p}_\gamma )/\beta _NE_Nϵ.$$
(1)   
where $`\stackrel{}{p}_N`$ and $`\stackrel{}{p}_\gamma `$ denote the nucleon and photon momenta. The Lorentz factor is $`\gamma _N=E_N/m_N=(1\beta _N^2)^{1/2}`$ with $`m_N`$ being the nucleon mass. The corresponding quantities in the NRF and CMF are marked with a prime () and an asterisk (), respectively. Fixed target accelerator experiments where a photon beam interacts with a proton target are performed in the NRF. In astronomical applications we assume that the LF can be chosen such that a the photon distribution function is isotropic. The LF may therefore be different from the astronomical observer’s frame if, for example the emission region is moving relative to the observer such as in AGN jets or GRB. The interaction rate of the nucleon in the LF is given by

$$R(E_N)=\frac{1}{8E_N^2\beta _N}_{ϵ_{\mathrm{th}}}^{\mathrm{}}𝑑ϵ\frac{n(ϵ)}{ϵ^2}_{s_{\mathrm{th}}}^{s_{\mathrm{max}}}𝑑s(sm_N^2)\sigma _{N\gamma }(s),$$
(2)   
where $`\sigma _{N\gamma }`$ is the total photohadronic cross section and

$$s=m_N^2+2E_Nϵ(1\beta _N\mathrm{cos}\theta )=m_N^2+2m_Nϵ^{},$$
(3)   
is the square of the center-of-mass energy. The lowest threshold energy for photomeson production is $`\sqrt{s_{\mathrm{th}}}=m_N+m_{\pi ^0}`$. The remaining quantities in Eq. (2) are $`ϵ_{\mathrm{th}}=(s_{\mathrm{th}}m_N^2)/2(E_N+p_N)`$ and $`s_{\mathrm{max}}=m_N^2+2E_Nϵ(1+\beta _N)`$.

The final state of the interaction is described by a number $`N_c`$ of possible channels, each of which has $`N_{f,c}`$ particles in the final state. At threshold, the phase space volume vanishes, which requires kinematically for the partial cross section $`\sigma _c0`$ for $`ss_{\mathrm{th},c}=[_im_i]^2`$. Above threshold, each final state channel has $`3N_{f,c}4`$ degrees of freedom given by the 3-momentum components $`(p_i,\chi _i,\varphi _i)`$ of the particles, constrained by energy and momentum conservation. Here $`p_i`$, $`\chi _i`$, and $`\varphi _i`$ are the particle momentum, and it’s polar and azimuthal angles with respect to the initial nucleon momentum, respectively. One of the $`\varphi _i`$ angles can be chosen to be distributed isotropically since we consider only the scattering of unpolarized photons and nucleons, all other variables are determined by the interaction physics through the differential cross sections. A distinguished role in the final state is played by the “leading-baryon”, which is considered to carry the baryonic quantum numbers of the incoming nucleon. For this particle, the Lorentz invariant 4-momentum transfer $`t=(P_NP_{\mathrm{final}})^2`$ is often used as a final state variable. At small $`s`$, many interaction channels can be reduced to 2-particle final states, for which $`d\sigma /dt`$ gives a complete description.

### 3.2 Interaction processes

Photon-proton interactions are dominated by resonance production at low energies. The incoming baryon is excited to a baryonic resonance due to the absorption of the photon. Such resonances have very short life times and decay immediately into other hadrons. Consequently, the $`N\gamma `$ cross section exhibits a strong energy dependence with clearly visible resonance peaks. Another process being important at low energy is the incoherent interaction of photons with the virtual structure of the nucleon. This process is called direct meson production. Eventually, at high interaction energies ($`\sqrt{s}>2\mathrm{GeV}`$) the total interaction cross section becomes approximately energy-independent, while the contributions from resonances and the direct interaction channels decrease. In this energy range, photon-hadron interactions are dominated by inelastic multiparticle production (also called multipion production).

#### 3.2.1 Baryon resonance excitation and decay

The energy range from the photopion threshold energy $`\sqrt{s}_{\mathrm{th}}1.08`$ GeV for $`\gamma N`$-interactions up to $`\sqrt{s}2`$ GeV is dominated by the process of resonant absorption of a photon by the nucleon with the subsequent emission of particles, i.e. the excitation and decay of baryon resonances. The cross section for the production of a resonance with angular momentum $`J`$ is given by the Breit-Wigner formula

$$\sigma _{\mathrm{bw}}(s;M,\mathrm{\Gamma },J)=\frac{s}{(sm_N^2)^2}\frac{4\pi b_\gamma (2J+1)s\mathrm{\Gamma }^2}{(sM^2)^2+s\mathrm{\Gamma }^2},$$
(4)   
where $`M`$ and $`\mathrm{\Gamma }`$ are the nominal mass and the width of the resonance. $`b_\gamma `$ is the branching ratio for photo-decay of the resonance, which is identical to the probability of photoexcitation. The decay of baryon resonances is generally dominated by hadronic channels. The exclusive cross sections for the resonant contribution to a hadronic channel with branching ratio $`b_c`$ can be written as

$$\sigma _c(s;M,\mathrm{\Gamma },J)=b_c\sigma _{\mathrm{bw}}(s;M,\mathrm{\Gamma },J),$$
(5)   
with $`_cb_c=1b_\gamma 1`$. Most decay channels produce two-particle intermediate or final states, some of them again involving resonances. For the pion-nucleon decay channel, $`N\pi `$, the angular distribution of the final state is given by

$$\frac{d\sigma _{N\pi }}{d\mathrm{cos}\chi ^{}}\underset{\lambda =J}{\overset{J}{}}\left|f_{\frac{1}{2},\lambda }^Jd_{\lambda ,\frac{1}{2}}^J(\chi ^{})\right|^2,$$
(6)   
where $`\chi ^{}`$ denotes the scattering angle in the CMF and $`f_{\frac{1}{2},\lambda }^J`$ are the $`N\pi `$-helicity amplitudes. The functions $`d_{\lambda ,\frac{1}{2}}^J(\chi ^{})`$ are commonly used angular distribution functions which are defined on the basis of spherical harmonics. The $`N\pi `$ helicity amplitudes can be determined from the helicity amplitudes $`A_{\frac{1}{2}}`$ and $`A_{\frac{3}{2}}`$ for photoexcitation (see Ref.  for details), which are measured for many baryon resonances . The same expression applies to other final states involving a nucleon and an isospin-0 meson (e.g., $`N\eta `$). For decay channels with other spin parameters, however, the situation is more complex, and we assume for simplicity an isotropic decay of the resonance.

Baryon resonances are distinguished by their isospin into $`N`$-resonances ($`I=\frac{1}{2}`$, as for the unexcited nucleon) and $`\mathrm{\Delta }`$-resonances ($`I=\frac{3}{2}`$). The charge branching ratios $`b_{\mathrm{iso}}`$ of the resonance decay follow from isospin symmetry. For example, the branching ratios for the decay into a two-particle final state involving a $`N`$\- or $`\mathrm{\Delta }`$-baryon and an $`I=1`$ meson ($`\pi `$ or $`\rho `$) are given in Table 1. Here $`\mathrm{\Delta }I_3`$ is the difference in the isospin 3-component of the baryon between initial and final state (the baryon charge is $`Q_B=I_3+\frac{1}{2}`$). In contrast to the strong decay channels, the electromagnetic excitation of the resonance does not conserve isospin. Hence, the resonance excitation strengths for $`p\gamma `$ and $`n\gamma `$ interactions are not related to each other by isospin symmetry and have to be determined experimentally.

#### 3.2.2 Direct pion production

Direct pion production can be considered as electromagnetic scattering by virtual charged mesons, which are the quantum–mechanical representation of the (color-neutral) strong force field around the baryon. The interacting virtual meson gains enough momentum to materialize. Experimentally the direct production of charged pions is observed as a relatively structureless background in the $`N\pi ^\pm `$ and $`\mathrm{\Delta }\pi ^\pm `$ final states in photon-nucleon interactions.

In terms of Feynman graphs, this process is represented by the $`t`$-channel exchange of a meson. Here, $`t`$ is the squared 4-momentum transfer from the initial to the final state baryon. The graph has a strong vertex at the baryon branch and an electromagnetic vertex for the photon interaction. At the strong vertex, the baryon may be excited and change its isospin. Isospin combination rules determine the iso-branching ratios in the same way as for resonance decay (Table 1 for $`I_{\mathrm{res}}=\frac{1}{2}`$). The presence of the electromagnetic vertex requires that the particle the photon couples to is charged. Thus direct processes with $`\mathrm{\Delta }I_3=0`$ branches (e.g. $`\gamma p\pi ^0p`$) are strongly suppressed.

The low energy structure of the direct cross section is not well constrained. At high energies, Regge theory of the pion exchange implies that $`\sigma _{\mathrm{dir}}(s)s^2`$ . The angular distribution of the process is strongly forward peaked and can be parametrized for small $`|t|`$ by

$$\frac{d\sigma _{\mathrm{dir}}}{dt}\mathrm{exp}(b_{\mathrm{dir}}t).$$
(7)   
with an experimentally determined slope of $`b_{\mathrm{dir}}12\mathrm{GeV}^2`$ .

The total cross section for a direct scattering process is roughly $`m_t^2`$, where $`m_t`$ is the (nominal) mass of the exchanged virtual particle. Therefore, the direct production of pions is dominant, while the contributions from the exchange of heavier mesons are suppressed. The same applies to direct reactions which involve the exchange of a virtual baryon ($`u`$-channel exchange). However, with increasing energy, more and more channels add to the direct cross section, and this makes an explicit treatment difficult.

#### 3.2.3 High energy processes

Phenomenologically, high energy interactions can be interpreted as reggeon and pomeron exchange processes. Both the reggeon and the pomeron are quasi-particles which correspond to sums of certain Feynman diagrams in the Regge limit ($`|t|s`$) . The cross sections for reggeon and pomeron exchange have different but universal energy dependences and account for all of the total cross section  at high energy. There are many different Regge theory-based cross section parametrizations possible. Here we use a recent cross section fit  based on the Donnachie-Landshoff model 

$$\sigma _{\mathrm{reg}}\left(\frac{sm_p^2}{s_0}\right)^{0.34}\sigma _{\mathrm{pom}}\left(\frac{sm_p^2}{s_0}\right)^{0.095},$$
(8)   
with the reference scale $`s_0=1`$ GeV<sup>2</sup>.

Concerning high energy processes, it is convenient to distinguish between diffractive and non-diffractive interactions. Diffractive interactions are characterized by the production of very few secondaries along the direction of the incoming particles. They correspond to the quasi-elastic exchange of a reggeon or pomeron between virtual hadronic states the photon (mainly the vector mesons $`\rho ^0`$, $`\omega `$, and $`\varphi `$) and the nucleon. Because of the spacelike nature of the interaction, the angular distribution is strongly forward peaked, and can be parametrized by Eq. (7) with an energy-dependent slope $`b_{\mathrm{diff}}=6\mathrm{G}\mathrm{e}\mathrm{V}^2+0.5\mathrm{GeV}^2\mathrm{ln}(s/s_0)`$ . At high energies, the cross section of diffractive interactions is approximately a constant fraction of the total cross section. The relative contribution of the different vector mesons is predicted by theory  to $`\rho ^0:\omega =9:1`$. The diffractive production of $`\varphi `$ or heavier mesons is neglected in SOPHIA.

Our treatment of non-diffractive multiparticle production is based on the Dual Parton Model . This model can be considered as a phenomenological realization of the expansion of QCD for large numbers of colors and flavours  in connection with general ideas of Duality and Gribov’s Regge theory . It provides a well developed basic scheme for the simulation of high energy hadronic interactions. The model can be visualized as follows: (i) the incoming nucleon and photon are split into colored quark and diquark constituents, (ii) in the course of the interaction these constituents exchange color quantum numbers, and (iii) confinement and the color field forces result in color strings which fragment to hadrons.

To relate the contributions of reggeon and pomeron exchange to parton configurations, we use the correspondence of their respective amplitudes to certain color flow topologies  which are shown in Figs. 2 and 2. The pomeron exchange topology involves the formation of two color neutral strings, while in case of a reggeon topology only one string is stretched from the diquark to the quark of the photon. The quark and diquark flavors at the string ends are determined by the the spin and valence flavor statistics for the nucleon. For photons the charge difference between $`u`$ and $`d`$ quarks increases the probability that the photon couples to a $`u`$-$`\overline{u}`$ pair instead of a $`d`$-$`\overline{d}`$ pair. In the model we use the theoretically predicted ratio of 4:1 between these two combinations.

The longitudinal momentum fractions $`x`$, $`1x`$ of the partons connected to the string ends are given by Regge asymptotics . One gets for the valence quark ($`x`$) and diquark ($`1x`$) distribution inside the nucleon

$$\rho (x)\frac{1}{\sqrt{x}}(1x)^{1.5}$$
(9)   
and for the quark antiquark distribution inside the photon

$$\rho (x)\frac{1}{\sqrt{x(1x)}}.$$
(10)   
The relatively small transverse momentum of the partons at the string ends are neglected.

In the string fragmentation process, the kinetic energy of the initial partons is reduced by creating new quark-antiquark pairs which are color field-connected to the parent partons. This process continues until the available kinetic energy drops below the particle production threshold causing the newly produced quarks to combine with the valence quarks to form hadrons (see, for example, ).

## 4 Implementation

In this section, all numerical expressions for cross sections are measured in units of $`\mu \mathrm{barn}`$, unless noted otherwise.

### 4.1 Method of cross section parametrization

The basic models of photohadronic interaction processes described in Sect. 3 are used to obtain robust theoretical predictions used for the parametrization of cross sections and final state distribution functions. For theoretically unpredictable parameters we use, if possible, the estimates given in the Review of Particle Properties (RPP. Remaining parameters are determined in fits to available exclusive and inclusive data on $`\gamma p`$ and $`\gamma n`$ interactions, as compiled in standard reference series (, and references therein). Since the parameters published in the RPP generally allow some variations within a given error range, these parameters are then optimized in comparison to data in an iterative process, until a reasonable agreement with the data for a large set of interaction channels is obtained.

This method has been previously described by Rachen , but was considerably improved in the development of SOPHIA. It provides a minimum bias description of photoproduction, which reproduces a large set of available data while reducing a possible bias due to data selection, since data are used only to fix model parameters. Considering the intended applications of SOPHIA for (i) astrophysical applications and (ii) the determination of background spectra in high energy experiments, we put particular emphasis on a good representation of inclusive cross sections and average final state properties in a wide range of interaction energies, while a good representation of complex exclusive channels is generally not expected.

### 4.2 Resonance production

Using Eq. (4), the contribution to the cross section from a resonance with mass $`M`$ and width $`\mathrm{\Gamma }`$ can be written as a function of the NRF photon energy $`ϵ^{}`$ as

$$\sigma (ϵ^{})=\frac{s}{ϵ^2}\frac{\sigma _0\mathrm{\Gamma }^2s}{(sM^2)^2+\mathrm{\Gamma }^2s}.$$
(11)   
The reduced cross section $`\sigma _0`$ is entirely determined by the resonance angular momentum and the electromagnetic excitation strength $`b_\gamma `$. We selected all baryon resonances listed in the RPP with certain existence (overall status: \****) and a well determined minimal photo–excitation strength of $`b_\gamma >10^4`$ for either the $`p\gamma `$ or the $`n\gamma `$ excitation. The resonances fulfilling these criteria and their parameters, as implemented in SOPHIA after iterative optimization, are given in Table 2. The phase-space reduction close to the $`N\pi `$ threshold is heuristically taken into account by multiplying Eq. (11) with the linear quenching function $`\mathrm{Qf}(ϵ^{};0.152,0.17)`$ for the $`\mathrm{\Delta }(1232)`$-resonance, and with $`\mathrm{Qf}(ϵ^{};0.152,0.38)`$ for all other resonances. The function $`\mathrm{Qf}(ϵ^{};ϵ_{\mathrm{th}}^{},w)`$ is defined in Appendix A. The quenching width $`w`$ has been determined from comparison with the data of the total $`p\gamma `$ cross section, and of the exclusive channels $`p\pi ^0`$, $`n\pi ^+`$ and $`\mathrm{\Delta }^{++}\pi ^{}`$ where most of the resonances contribute.

The major hadronic decay channels of these baryon resonances are $`N\pi `$, $`\mathrm{\Delta }\pi `$, and $`N\rho `$; for the $`N(1535)`$, there is also a strong decay into $`N\eta `$, and the $`N(1650)`$ contributes to the $`\mathrm{\Lambda }K`$ channel. The hadronic decay branching ratios $`b_c`$ are all well determined for these resonances and given in the RPP. However, a difficulty arises from the fact that branching ratios can be expected to be energy dependent because of the different masses of the decay products in different branches. In SOPHIA, we consider all secondary particles, including hadronic resonances, as particles of a fixed mass. This implies that, for example, the decay channel $`\mathrm{\Delta }\pi `$ is energetically forbidden for $`\sqrt{s}<m_\mathrm{\Delta }+m_\pi 1.37`$ GeV. To accommodate this problem, we have developed a scheme of energy dependent branching ratios, which change at the thresholds for additional decay channels and are constant in between. The requirements are that (i) the branching ratio $`b_c=0`$ for $`ϵ^{}<ϵ_{\mathrm{th},c}^{}`$, and (ii) the average of the branching ratio over energy, weighted with the Breit-Wigner function, correspond to the average branching ratio given in the RPP for this channel. For all resonances, we considered not more than three decay channels leading to a unique solution to this scheme. No fits to data are required. In practice, however, the experimental error on many branching ratios allows for some freedom, which we have used to generate a scheme that optimizes the agreement with the data on different exclusive channels.

The hadronic branching ratios are given in Tab. 4 in Appendix B. To obtain the contribution to a channel with given particle charges, e.g. $`\mathrm{\Delta }^{++}\pi ^{}`$, the hadronic branching ratio $`b_{\mathrm{\Delta }\pi }`$ has to be multiplied with the iso-branching ratios as given in Tab. 1. We note that with the parameters $`b_\gamma `$, $`b_c`$ and $`b_{\mathrm{iso}}`$, the resonant contribution to all exclusive decay channels is completely determined.

The angular decay distributions for the resonances follow from Eq. (6). In SOPHIA, the kinematics of the decay channels into $`N\pi `$ is implemented in full detail (see Tab. 3). For other decay channels, we assume isotropic decay according to the phase space. Furthermore, there might be some mixing of the different scattering angular distributions since the sampled resonance mass, in general, does not coincide with its nominal mass. This effect is neglected in our work. Instead, we use the angular distributions applying to resonance decay at its nominal mass $`M`$.

The two decay products of a resonance may also decay subsequently. This decay is simulated to occur isotropically according to the available phase space.

### 4.3 Direct pion production

The cross section for direct meson production, unlike those of resonances, is not completely determined by well known parameters. The low and high energy constraints suggest the phenomenological parametrization

$$\sigma _{\mathrm{dir}}(ϵ^{})=\sigma _{\mathrm{max}}\mathrm{Pl}(ϵ^{};ϵ_{\mathrm{th}}^{},ϵ_{\mathrm{max}}^{},2),$$
(12)   
where the function $`\mathrm{Pl}(ϵ^{};ϵ_{\mathrm{th}}^{},ϵ_{\mathrm{max}}^{},\alpha )`$ approaches zero for $`ϵ^{}=ϵ_{\mathrm{th}}^{}`$, goes through a maximum at $`ϵ^{}=ϵ_{\mathrm{max}}^{}`$ and follows an asymptotic behaviour $`(ϵ^{})^\alpha `$. The definition of this function is given in Appendix A.

In SOPHIA, we consider explicitly direct channels with charged pion exchange which are dominating at low energies. The selection is further constrained by the fact that sufficient data are only available for the channels $`p\gamma n\pi ^+`$, $`n\gamma p\pi ^{}`$, and $`p\gamma \mathrm{\Delta }^{++}\pi ^{}`$. We note that proton and neutron induced direct reactions are strictly isospin-symmetric. Both proton and neutron data sets (when available) can be used in the fitting procedure. The high energy data fits on the $`\mathrm{\Delta }\pi `$ and $`N\pi `$ channel, i.e. $`\sigma _\pi 18(ϵ^{})^2`$  and $`\sigma _\mathrm{\Delta }26.4(ϵ^{})^2`$  for $`ϵ^{}>1`$, were primarily used to fix $`\sigma _{\mathrm{max}}`$, while a best fit of $`ϵ_{\mathrm{max}}^{}`$ was obtained by comparing with the residuals of the low energy data after subtracting the resonance contribution. The adopted cross sections are

$`\sigma _{N\pi }(ϵ^{})`$ $`=`$ $`92.7\mathrm{Pl}(ϵ^{};0.152,0.25,2)+40\mathrm{exp}\left({\displaystyle \frac{(ϵ^{}0.29)^2}{0.002}}\right)`$   
$`15\mathrm{exp}\left({\displaystyle \frac{(ϵ^{}0.37)^2}{0.002}}\right),`$   
$`\sigma _{\mathrm{\Delta }\pi }(ϵ^{})`$ $`=`$ $`37.7\mathrm{Pl}(ϵ^{};0.4,0.6,2).`$ (14)   
The two Gaussian-shaped functions included in the direct $`N\pi `$ cross section have been added to improve the representation of the total cross section in the energy region $`0.152`$ GeV $`<ϵ^{}<0.4`$ GeV, where otherwise only the well constrained $`\mathrm{\Delta }(1232)`$ resonance contributes significantly. For $`p\gamma `$\- ($`n\gamma `$-) interactions $`\sigma _{\mathrm{\Delta }\pi }`$ contributes to the $`\mathrm{\Delta }^{++}\pi ^{}`$ ($`\mathrm{\Delta }^+\pi ^{}`$) and $`\mathrm{\Delta }^0\pi ^+`$ ($`\mathrm{\Delta }^{}\pi ^+`$) final states with a ratio 3:1 according to isospin combination rules (see Tab. 1).

By comparison with the total cross section data we find that the resonant and direct interaction channels account for all of the total interaction cross section below the $`3\pi `$ threshold at $`ϵ^{}0.51`$ GeV. Above this threshold, and below the threshold for diffractive interactions at $`ϵ^{}1`$ GeV, where high energy processes set in, we find a residual cross section which can be parametrized as

$$\sigma _{\mathrm{lf}}=80.3(60.2)\mathrm{Qf}(x;0.51;0.1)(ϵ^{})^{0.34},$$
(15)   
where the number 60.2 given in brackets belongs to $`n\gamma `$–interactions while the number 80.3 refers to $`p\gamma `$–collisions. The normalization cross section and the quenching width has been determined by a $`\chi ^2`$ minimization method to the total cross section data for $`p\gamma `$ ($`n\gamma `$) interactions after subtraction of the respective resonant and direct contributions. By analogy, the power law index for this contribution is taken from the high energy parametrization for reggeon exchange (note that $`ϵ^{}sm_N^2`$). Physically, this cross section represents the joint contribution of all $`t`$-channel scattering processes at low energies not considered so far. This is in principle similar to interactions at high energies. Consequently, we use an adapted string fragmentation model to simulate this contribution, and refer to it as low energy fragmentation hereafter.

### 4.4 High energy multipion production

In SOPHIA, we assume that the cross sections for diffractive and non-diffractive high energy interactions are proportional to each other at all energies. This assumption fixes the threshold for high energy interactions to the threshold of the $`N\rho `$ final state, which is nominally at $`ϵ^{}1.1`$ GeV. Because of the large width of the $`\rho `$ there should be some contribution also at lower energies. From comparison with exclusive data of the $`N\rho `$ final state, and the residuals of the total cross section data, with the sum of the contributions of all low energy channels, we find a common threshold for high energy interactions of $`ϵ_{\mathrm{th},\mathrm{high}}^{}=0.85`$ GeV.

We restrict the diffractive channel to the non-resonant production of $`N\rho `$ and $`N\omega `$ final states, for which we assume the theoretically predicted relation $`\sigma _\rho =9\sigma _\omega `$. The ratio between diffractive and non-diffractive interactions is derived from the comparison with exclusive $`N\rho `$ data and total cross section data at high energy,

$$\sigma _{\mathrm{diff}}=0.15\sigma _{\mathrm{frag}}.$$
(16)   
For the parametrization of $`\sigma _{\mathrm{frag}}`$, we use the power law representations of the reggeon and pomeron exchange cross section at high energies, and multiply them by an exponential quenching function $`1\mathrm{exp}([ϵ_{\mathrm{th},\mathrm{high}}^{}ϵ^{}]/a)`$. The relative contributions of the reggeon and pomeron cross sections, and the quenching parameter $`a`$ have been determined by a iterative $`\chi ^2`$ minimization method with respect to the total $`p\gamma `$ ($`n\gamma `$) cross section data after subtraction of all low energy contributions. We find

$$\sigma _{\mathrm{frag}}(ϵ^{})=\left[1\mathrm{exp}\left(\frac{ϵ^{}0.85}{0.69}\right)\right][28.8(26.0)(ϵ^{})^{0.34}+58.3(ϵ^{})^{0.095}],$$
(17)   
where we have used the high-energy behaviour given by Eq. (8).

The string fragmentation is done by the Lund Monte Carlo Jetset 7.4 . This program is well suited for string fragmentation at high energies. Since for our purposes also strings with rather small invariant masses have to be hadronized, several parameters of this fragmentation code had to be tuned to obtain a reasonable description also at low energies. Furthermore, in order to avoid double counting, all final states identical to the processes already considered by resonance production and direct interactions are rejected (note that this is not the case for low-energy fragmentation).

### 4.5 Initial state kinematics and photon radiation fields

The probability for interaction of a proton of energy $`E_N`$ with a photon of energy $`ϵ`$ from a radiation field with the photon density $`n(ϵ)`$ reads

$$𝒫(ϵ)=\frac{1}{R(E_N)}\frac{n(ϵ)}{8E_N^2\beta ϵ^2}_{s_{\mathrm{th}}}^{s_{\mathrm{max}}}𝑑s(sm_N^2)\sigma _{N\gamma }(s),$$
(18)   
where $`R(E_N)`$ is the interaction rate as given in Eq. (2), where also $`s_{\mathrm{th}}`$ and $`s_{\mathrm{max}}`$ are defined.

For a fixed nucleon energy, the CMF energy is sampled from the distribution

$$𝒫(s)=\mathrm{\Phi }^1(sm_N^2)\sigma _{N\gamma }(s)$$
(19)   
with $`\mathrm{\Phi }=_{s_{\mathrm{th}}}^{s_{\mathrm{max}}}𝑑s(sm_N^2)\sigma _{N\gamma }(s)`$. The interaction angle follows from

$$\mathrm{cos}\theta =\frac{1}{\beta }\left(\frac{m_N^2s}{2E_Nϵ}+1\right).$$
(20)   
Currently black body, power law, and broken power law radiation spectra are implemented in SOPHIA. The photon density $`n(ϵ)`$ for a blackbody radiation field of temperature $`T`$ is given in natural units by

$$n(ϵ)=\frac{1}{\pi ^2}\frac{ϵ^2}{\mathrm{exp}(\frac{ϵ}{kT})1},$$
(21)   
where $`k`$ is the Boltzmann constant. For a power law photon spectrum the photon density is given by $`n(ϵ)=ϵ^\alpha `$. The broken power law photon spectrum is given by

$`n(ϵ)=`$ $`ϵ^{\alpha _1}`$ $`\text{for}ϵ<ϵ_b`$ (22)   
$`n(ϵ)=`$ $`ϵ_b^{\alpha _2\alpha _1}ϵ^{\alpha _2}`$ $`\text{for}ϵ>ϵ_b`$ (23)   
where $`ϵ_b`$ is the break energy, and $`\alpha _1`$, $`\alpha _2`$ are the photon indices below and above the break energy, respectively. Note that no absolute normalization of $`n(ϵ)`$ is necessary since it cancels in the definition of $`𝒫(ϵ)`$.

## 5 Structure of the program

### 5.1 Source code

The SOPHIA source code consists of several files which contain a number of routines.

main program containing the routines which organize the various tasks to be performed. Furthermore, the input is handled here and some kinematic transformations needed as input to several routines are performed.

initialization routine for parameter settings.

collection of routines/functions needed for sampling the CMF energy squared $`s`$ and the photon energy $`ϵ`$ in the observer frame.

event generator for photomeson production in $`p\gamma `$ and $`n\gamma `$ collisions.

contains output routines.

### 5.2 The event generator

The simulation of the final state is performed by the photopion production event generator EVENTGEN. Together with the initialization routine INITIAL, EVENTGEN can be used separately for Monte Carlo event simulation. The user has to give the nucleon code number L0, energy E0 (in GeV), the photon of energy $`ϵ`$ (in GeV), and the angle $`\theta `$ in degrees.

EVENTGEN is structured as follows. First the momenta of the incident particles are Lorentz-transformed into the CMF of the interactions. The cross section (in $`\mu `$barn) is calculated in the function CROSSECTION. The cross sections for the various channels considered in this code determine the distribution for the probability of a certain process. The sampling of a process (resonance decay, direct channel, diffractive scattering, fragmentation) at a given NRF energy of the photon is carried out in the routine DEC\_INTER3.

For the resonance decay we sample the resonance at this energy using the Breit-Wigner formula as a probability distribution for a specific resonance (performed in DEC\_RES2). Its branching ratios define the decay mode (in DEC\_PROC2). The subsequent two-particle decay in the CM frame is carried out in RES\_DECAY3, and then decays of all unstable particles are carried out in the Sybill routine DECSIB.

Secondary particle production is simulated in GAMMA\_H for direct and multiparticle production.

Finally, all final state particles are Lorentz-transformed back to the LF.

The output of the final states is organized in the common block /S\_PLIST/ P(2000,5), LLIST(2000), NP, Ideb. Here the array P(i,j) contains the 4-momenta and rest mass of the final state particle i in cartesian coordinates (P(i,1) = $`P_x`$, P(i,2) = $`P_y`$, P(i,3) = $`P_z`$, P(i,4) = energy, P(i,5) = rest mass). LLIST() gives the code numbers of all final state particles and NP is the number of stable final state particles.

### 5.3 Input/Output routines

Using the standard main program, the user provides the following input parameters.

* E0 = energy of incident proton (in GeV), or
* Emin, Emax = low/high energy cutoff of an energy grid of incident protons (in GeV)
* L0 = code number of the incident nucleon (L0 = 13: proton, L0 = 14: neutron)
* ambient photon field:

+ blackbody spectrum: TBB = temperature (in K)
    + straight/broken power law spectrum: ALPHA1, ALPHA2 = power law indices, EPSMIN = low energy cut off (in eV), EPSMAX = high energy cut off (in eV), EPSB = break energy (in eV)
* NTRIAL = number of inelastic interactions
* NBINS = number of bins for output particle spectra ($`200`$ bins)
* DELX = stepsize of output particle spectra

For the calculation of the incident particle momenta we assume that the relativistic nucleon is moving along the positive $`z`$-axis.

The output is organized as follows. All the energy distributions $`(1/N_{\mathrm{eve}})dN_{\mathrm{part}}/d\mathrm{log}x`$ of produced particles are given with logarithmically equal bin sizes in the scaling variable $`x=E_{\mathrm{part}}`$/E0. Here $`N_{\mathrm{eve}}`$ denotes the number of simulated inelastic events and $`N_{\mathrm{part}}`$ is the number of secondary particles of a certain kind. The spectra of photons, protons, neutrons, $`e`$-neutrinos, $`e`$-antineutrinos, $`\mu `$-neutrinos and $`\mu `$-antineutrinos are considered separately. They are stored in a file with name xxxxxx.particle with xxxxxx = input name (chosen by the user):

| particle = | ’gamma’ | $``$ | $`\gamma `$ | spectrum |
| --- | --- | --- | --- | --- |
|  | ’muneu’ | $``$ | $`\nu _\mu `$ | spectrum |
|  | ’muane’ | $``$ | $`\overline{\nu }_\mu `$ | spectrum |
|  | ’e\_neu’ | $``$ | $`\nu _e`$ | spectrum |
|  | ’eaneu’ | $``$ | $`\overline{\nu }_e`$ | spectrum |
|  | ’elect’ | $``$ | $`e^{}`$ | spectrum |
|  | ’posit’ | $``$ | $`e^+`$ | spectrum |
|  | ’proto’ | $``$ | $`p`$ | spectrum |
|  | ’neutr’ | $``$ | $`n`$ | spectrum |

The structure of a typical output file is:   

| 1. line: | high/low energy cutoff of incident nucleon energy grid, |
| --- | --- |
|  | number of energy bins of incident nucleon spectrum (=NINC), |
|  | TBB or ALPHA1, ALPHA2, EPSMIN, EPSB, EPSMAX, incident particles |

| 2. line: | 1. number: | energy of incident nucleon |
| --- | --- | --- |
|  | 2. number: | first number (=$`a`$) of non-zero energy bin of |
|  |  | particle spectrum |
|  | 3. number: | last number (=$`b`$) of non-zero energy bin of |
|  |  | particle spectrum |
| 3. line: | particle spectrum |  |
|  | between $`a`$$`b`$ |  |

## 6 Comparison to data

### 6.1 Total cross section

Fig. 3 (upper panel) shows the total cross section for $`\gamma p`$-interactions with the various contributions considered in SOPHIA. For simplicity, we show both fragmentation contributions together (low energy fragmentation and non-diffractive multipion production). The resonant process $`\gamma p\mathrm{\Delta }^+\pi ^0p`$ is the only one kinematically possible directly at the particle production threshold. Above the $`\pi ^+n`$ threshold at very low energies ($`ϵ^{}<0.25\mathrm{GeV}`$), the direct channel $`\gamma p\pi ^+n`$ constitutes the largest contribution.

To assess the quality of the cross section parametrization, the differences between the experimental data on the total $`\gamma p`$ cross section and the cross section fit are shown in Fig. 3 (lower panel). The total $`\gamma n`$ cross section is overall similar, except at energies of about $`\sqrt{s}1.680`$ GeV ($`ϵ^{}1.035`$ GeV), where it is considerably smaller due to the different excitation strengths of the resonances at this energy.

### 6.2 Exclusive cross sections

Figs. 4 to 7 compare the output of SOPHIA with the data on specific final states as a function of the interaction energy. Such comparisons are important for models that aim to represent correctly photohadronic interactions over a wide energy range.

Fig. 4 compares the total cross sections for $`\pi ^0p`$ and $`\pi ^+n`$ production with experimental data. The major contributions come from the $`\mathrm{\Delta }`$(1232) resonance and the direct channel together with minor contributions from other resonances. The agreement with data in the threshold region is of great importance for many astrophysical applications where this is the dominating energy range in the case of steep proton and ambient photon spectra.

Fig. 5 compares the calculated and measured cross sections for final states involving charged and neutral pions. The general agreement of the model results with data is quite good, although there are some energy ranges that show minor deviations. It is important to note that it is difficult to fit exactly the experimental data without any detailed knowledge of the experimental setups and acceptance constraints, especially in cases like $`\gamma p\pi ^+\pi ^{}p`$ \+ neutrals, where the final state is not well defined.

The number of inelastic $`\gamma p`$ events with 1, 3, and 5 charged particles in the final state are shown in Fig. 6. This comparison is very sensitive to the description of multipion production processes as well as to the smooth transition between different particle production processes. It shows that the different channels are reasonably well modeled in SOPHIA.

Fig. 7 shows reaction cross sections for final states with equal numbers of pions, but having a different isospin of the produced nucleon, together with fixed-target data. This comparison is important for the correct simulation of the ratio of the proton to neutron numbers produced in $`\gamma p`$ collisions. Again the model results agree well with data.

### 6.3 Inclusive distributions

The rapidity distribution tests the kinematic of our simulations. The rapidity of a final state particle with energy $`E`$ is defined as

$$y=\frac{1}{2}\mathrm{ln}\left(\frac{E+p_{}}{Ep_{}}\right)=\mathrm{ln}\left(\frac{E+p_{}}{m_{}}\right),$$
(24)   
where $`p_{}`$ is its momentum component along the direction of the incoming particle. The transverse mass $`m_{}`$ follows from $`m_{}^2=E^2p_{}^2`$. Rapidity is additive under Lorentz transformations, which keeps the rapidity distribution invariant under such transformations.

Moffeit et al.  have measured the rapidity distribution for the interaction $`\gamma p\pi ^{}+`$ anything at three different beam energies $`ϵ^{}`$ = 2.8, 4.7 and 9.3 GeV. In Fig. 8 we compare the calculated rapidity distribution to data. The agreement in the width and the height of the distributions is good.

### 6.4 Multiplicities

In astrophysical environments the production of neutrinos results mainly from the decay of charged secondary pions ($`\pi ^+e^+\nu _\mu \overline{\nu }_\mu \nu _e`$, $`\pi ^{}e^{}\nu _\mu \overline{\nu }_\mu \overline{\nu }_e`$). Most of the gamma rays are produced via neutral pion decay ($`\pi ^0\gamma \gamma `$) and synchrotron/Compton emission from emanating relativistic $`e^{}/e^+`$. Pion multiplicities (see Figs. 9, 10) are therefore instructive to understand the energy dissipation of the initial energy among the $`\nu `$\- and $`\gamma `$-component. For non-astrophysical applications the charged pion multiplicity is a basic interaction parameter that presents a cumulative measure of many interaction channels.

In the resonance region the maximum of the neutral pion multiplicity is reached at the $`\mathrm{\Delta }(1232)`$-resonance (see Fig. 9). At threshold neutral pion production is strongly suppressed in favour of $`\pi ^+`$-production (see Fig. 10) due to the dominance of the direct channel. Multiplicities are approximately growing as $`s^{1/4}`$ in the multipion production region at low energies. The obtained multiplicity distributions from our simulations are in agreement with the data from Moffeit et al. .

By counting the number of protons and neutrons produced in $`\gamma p`$ interactions, one can define a proton-to-neutron ratio. The SOPHIA prediction on the energy dependence of this ratio is shown in Fig. 11. The $`p/n`$ ratio reaches $`2.2`$ at high energy which can be contrasted to an experimentally estimated value of about $`3.8`$ found by Meyer . The ratio derived in  is essentially based on the same data as those our model is compared to. Meyer estimated the experimentally unknown cross sections by isospin symmetry arguments. In our case these cross sections are predicted by the Monte Carlo simulation. We conclude from this that the unmeasured cross sections are important for the $`p/n`$ ratio, and that the difference between the values 2.2 and 3.8 reflect the uncertainty due to the limited experimental data available.

## 7 Conclusions

A newly developed Monte Carlo event generator SOPHIA has been presented. It simulates the interactions of nucleons with photons over a wide range in energy. The simulation of the final state includes all interaction processes which are relevant to astrophysical applications. SOPHIA contains tools, such as for sampling the photon energy from different ambient soft photon spectra, and the nucleon–photon interaction angle, that are needed for such applications. As an event generator, SOPHIA uses all available information on the interaction cross section, the final state particle composition, and kinematics of the interaction processes as provided by particle physics. Comparison with the available accelerator data shows that SOPHIA provides a good description of our current knowledge of photon–nucleon interactions.

Acknowledgements

The work of AM and RJP is supported by the Australian Research Council. RE and TS acknowledge the support by the U.S. Department of Energy under Grant Number DE FG02 01 ER 40626. TS is also supported in part by the NASA grant NAG5–7009. The contribution of JPR was supported by NASA NAG5–2857 and by the EU-TMR network “Astro-Plasma-Physics” under the contract number ERBFMRX-CT98-0168.

## Appendix A Definition of functions

The functions $`\mathrm{Pl}(ϵ^{},ϵ_{\mathrm{th}}^{},ϵ_{\mathrm{max}}^{},\alpha )`$ and $`\mathrm{Qf}(ϵ^{},ϵ_{\mathrm{th}}^{},w)`$ are defined by

$$\mathrm{Pl}(ϵ^{},ϵ_{\mathrm{th}}^{},ϵ_{\mathrm{max}}^{},\alpha )=(\frac{ϵ^{}ϵ_{\mathrm{th}}^{}}{ϵ_{\mathrm{max}}^{}ϵ_{\mathrm{th}}^{}})^{A\alpha }(\frac{ϵ^{}}{ϵ_{\mathrm{max}}^{}})^A$$
(25)   
for $`ϵ^{}>ϵ_{\mathrm{th}}^{}`$, and $`\mathrm{Pl}(ϵ^{},ϵ_{\mathrm{th}}^{},ϵ_{\mathrm{max}}^{},\alpha )=0`$ otherwise, $`A=\alpha ϵ_{\mathrm{max}}^{}/ϵ_{\mathrm{th}}^{}`$, and

$$\mathrm{Qf}(ϵ^{},ϵ_{\mathrm{th}}^{},w)=(ϵ^{}ϵ_{\mathrm{th}}^{})/w$$
(26)   
for $`ϵ_{\mathrm{th}}^{}<ϵ^{}<w+ϵ_{\mathrm{th}}^{}`$, $`\mathrm{Qf}(ϵ^{},ϵ_{\mathrm{th}}^{},w)=0`$ for $`ϵ^{}ϵ_{\mathrm{th}}^{}`$, and $`\mathrm{Qf}(ϵ^{},ϵ_{\mathrm{th}}^{},w)=1`$ for $`ϵ^{}w+ϵ_{\mathrm{th}}^{}`$.

## Appendix B Resonance branching ratios

## Appendix C Compilation of routines/functions

calculates cross section (in $`\mu `$barn) of a resonance with width GAMMA (in GeV), mass DMM (in GeV), maximum cross section SIGMA\_0 (in $`\mu `$barn) and NRF energy of the photon (in GeV) according to the Breit-Wigner formula

collection of functions (SINGLEBACK, TWOBACK) which calculates the cross section (at the NRF energy EPS\_PRIME of the incident photon) of the direct channel (not isospin-corrected)

computes cross section (in $`\mu `$barn) for $`N\gamma `$-interaction at a given energy EPS\_PRIME (=photon energy in GeV in proton’s rest frame); depending on the control variable NDIR it returns the total cross section (NDIR=3) or only a certain part of the cross section (NDIR=1: total resonance cross section, NDIR=4: direct channel cross section NDIR=5: fragmentation cross section, NDIR=11-19: individual resonance cross sections)

returns reaction mode: decay of resonance (IMODE=6), direct pion production (IMODE=2 for $`N\pi `$ final states, IMODE=3 for $`\mathrm{\Delta }\pi `$ final states), fragmentation in resonance region (IMODE=5), diffractive scattering (IMODE=1 for $`N\rho `$ final states, IMODE=4 for $`N\omega `$ final states) and multipion production/fragmentation (IMODE=0) at a given energy EPS\_PRIME (in GeV)

returns in IPROC the decay mode for a given resonance IRES at energy EPS\_PRIME (in GeV) for incident nucleon with code number L0; IRANGE is the number of energy intervals corresponding to the (energy dependent) branching ratios of a specific resonance

returns sampled resonance number IRES (IRES = 1 …IRESMAX) at energy EPS\_PRIME (in GeV) for incident nucleon with code number L0

decay of unstable particles

calculates distribution of the squared CMF energy S (in GeV<sup>2</sup>)

interface routine for hadron-$`\gamma `$ collisions in CM frame;   
E0 is CMF energy of the $`N\gamma `$-system; first particle ($`\gamma `$ or hadron $`N`$) goes in $`+z`$-direction; final state consists of protons, neutrons, gamma, leptons, neutrinos

initialization routine for parameter settings;   
to be called at before calling EVENTGEN

interface to Jetset fragmentation of system with CMF energy SQS (in GeV);   
NP = number of secondaries produced

returns photon density (in photons/cm<sup>3</sup>/eV) at energy EPS (in eV) of blackbody radiation of temperature TBB (in K)

probability distribution for scattering angle of given resonance IRES and incident nucleon with code number L0 ;   
ANGLESCAT is the cosine of scattering angle in CMF frame

calculates probability distribution for thermal photon field with temperature TBB (in K) at energy EPS (in eV)

calculates probability distribution at energy EPS (in eV) for power law radiation n $``$ EPS<sup>-ALPHA</sup> between the limits EPSM1 …EPSM2 (in eV)

carries out 2-particle decay in CM frame of system with mass AMD (in GeV) into particles with code numbers LA,LB (stored in array LRES()) wheras COSTHETA is the cosine of the $`N\gamma `$-CM frame scattering angle ; returns 5-momenta PRES() of the two particles;   
NBAD=1 for kinematically not possible decays, otherwise NBAD=0

determines decay products for the decay of a given resonance IRES by a given decay process IPROC for a incident nucleon with code number L0, and carries out two-particle decay of the resonance with squared CMF energy S; code numbers and (CMF) 5-momenta of produced particles are stored in arrays LLIST and P, respectively, in a common block;   
NBAD=1 for kinematically not possible decays; otherwise NBAD=0

samples incident photon energy; if the blackbody temperature is set to TBB=0 the photon energy is sampled from a power law distribution with index ALPHA between given photon energies EPSMIN, EPSMAX (in eV); otherwise it is sampled from a distribution for a blackbody radiation with temperature TBB (in K)

samples the total CM frame energy S (in GeV<sup>2</sup>) for a given photon with energy EPS (in GeV) and proton with energy E0 (in GeV)

samples the cosine of the scattering angle ANGLESCAT for a given resonance IRES and incident nucleon INC