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 s^1/2~10³ 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 s^1/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

Coding:

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 s^1/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 = 10¹¹GeV, 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.