SOPHIA Documentation

General Information

SOPHIA is designed to simulate minimum bias photo-production by interactions of unpolarized nucleons with unpolarized photons from the particle production threshold through the resonance region up to s1/2~103 in the multiparticle production region. The implemented interaction processes include: incoherent interaction of photons with the virtual structure of the nucleon at low energies near threshold; resonance excitation/decay (9 resonances with masses 1232-1950 MeV); diffractive scattering (vector meson production) and multipion production on the basis of QCD string fragmentation at high energies s1/2>2GeV. A modified JETSET version 7.4 is used for the latter high energy event generation. SOPHIA simulates individual events for given nucleon and photon energies where photon energies are sampled from various distributions. Corresponding incident ambient photon fields are limited to power law and blackbody spectra in the public code version. A detailed description of all processes and the SOPHIA code can be found in A. Mücke, Ralph Engel, J.P. Rachen, R.J. Protheroe, and Todor Stanev, 2000, Comp. Phys. Commun., 124, 290.


Program structure


SOPHIA is coded in FORTRAN 77 (operating systems : UNIX, Linux, Open-VMS) and has been tested thoroughly on DEC-Alpha and Intel-Pentium based workstations. No tests were done for center-of-mass energies s1/2>1000GeV. Memory required to execute is <1 megabyte. 10000 events at a center-of-mass energy of 1.5 GeV require a typical CPU time of about 75 seconds. Other programs used in SOPHIA in a modified form are: Rndm (a processor independent random number generator based on Marsaglia & Zaman 1987, preprint FSU-SCRI-87-70), Jetset 7.4 (a Lund Monte Carlo for jet fragmentation; Sjostrand 1994, Comp.Phys.Commun., 82, 74), DECSIB (Sibyll-routine which decays unstable particles; Fletcher et al., 1994, Phys.Rev.D 50, 5710).

The SOPHIA source code consists of several files which contain a number of routines. SOPHIA20.f 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. initial.finitialization routine for parameter settings. sampling.f collection of routines/functions needed for sampling the CMF energy squared s and the photon energy eps in the observer frame. eventgen.f event generator for photomeson production in proton-photon and neutron-photon collisions. This is the heart of SOPHIA.

The output of the final states is organized in the common block S_PLIST /P(2000,5), LLIST(2000), NP, Ideb/the array P(i,j) contains the 4-momenta and rest mass of the final state particle 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 NPis the number of stable final state particles. Note that before the initial call of eventgen the initialization routine has to be called. output.fcontains output routines.

Details on each routine can be found in A. Mücke et al., 2000, Comp. Phys. Commun., 124, 290.

The particle identity is given by LLIST with the following numbering scheme (identical to the DECSIB particle identity list):

code particle mass
1 gam 0.0000
2 e+ 0.0005
3 e- 0.0005
4 mu+ 0.1057
5 mu- 0.1057
6 pi0 0.1350
7 pi+ 0.1396
8 pi- 0.1396
9 k+ 0.4937
10 k- 0.4937
11 k0l 0.4977
12 k0s 0.4977
13 p 0.9383
14 n 0.9396
15 nue 0.0000
16 nueb 0.0000
17 num 0.0000
18 numb 0.0000
21 k0 0.4977
22 k0b 0.4977
23 eta 0.5488
24 etap 0.9576
25 rho+ 0.7714
26 rho- 0.7714
27 rho0 0.7717
28 k*+ 0.8921
29 k*- 0.8921
30 k*0 0.8965
31 k*0b 0.8965
32 omeg 0.7826
33 phi 1.1020
34 SIG+ 1.1894
35 SIG0 1.1925
36 SIG- 1.1973
37 XI0 1.3149
38 XI- 1.3213
39 LAM 1.1156
40 DELT++ 1.2300
41 DELT+ 1.2310
42 DELT0 1.2320
43 DELT- 1.2330
44 SIG*+ 1.3828
45 SIG*0 1.3837
46 SIG*- 1.3872
47 XI*0 1.5318
48 XI*- 1.5350
49 OME*- 1.6724

Antibaryons have negative ID numbers correspondingly. Decayed particles are marked by adding 10000 to their ID code.The decay of a particle with code ID can be turned off viaIDB(ID) = -ABS(IDB(ID))and turned on viaIDB(ID) = ABS(IDB(ID))This has to be done only once at the beginning of event generation or might be changed on an event-by-event basis.

Example Programs

The following are example files of input parameter sets and resulting SOPHIA output files: 

History of Updates

Version 2.01 (April 12, 2007)

Last modification on March 27, 2007:

An error in the routine eventgen has been corrected (pointed out by P. Lipari): The problem comes from inconsistent sign conventions for the cos(theta) specified in the parameters. The definition of the explicit particle momenta is only used in the Lorentz transformation. The relevant quantity affected is the boost parameter along the z-axis.

Since the photon momentum is much smaller than the nucleon momentum in "standard" applications in astrophysics, the bug has no effect on simulation results. For a typical UHECR propagation configuration of E_p = 1011GeV, E_ph = 10-10GeV, the relative error would be of the order of 10-20. This applies to isotropic or anisotropic radiation fields.

The situation is different for configurations in which the photon and nucleon energies are of the same order. In this case the total energy would be still fine but the total z-momentum would be wrong.