1 C*********************************************************************
2 C*********************************************************************
6 C* The Lund Monte Carlo for Hadronic Processes **
8 C* PYTHIA version 5.7 **
10 C* Torbjorn Sjostrand **
11 C* Department of theoretical physics 2 **
12 C* University of Lund **
13 C* Solvegatan 14A, S-223 62 Lund, Sweden **
14 C* E-mail torbjorn@thep.lu.se **
15 C* phone +46 - 46 - 222 48 16 **
17 C* Several parts are written by Hans-Uno Bengtsson **
18 C* CTEQ 2 parton distributions are by the CTEQ collaboration **
19 C* SaS photon parton distributions together with Gerhard Schuler **
20 C* g + g -> Z + b + bbar matrix element code by Ronald Kleiss **
21 C* g + g and q + qbar -> t + tbar + H code by Zoltan Kunszt **
23 C* The latest program version and documentation is found on WWW **
24 C* http://thep.lu.se/tf2/staff/torbjorn/Welcome.html **
26 C* Copyright Torbjorn Sjostrand and CERN, Geneva 1993 **
28 C*********************************************************************
29 C*********************************************************************
31 C List of subprograms in order of appearance, with main purpose *
32 C (S = subroutine, F = function, B = block data) *
34 C S PYINIT to administer the initialization procedure *
35 C S PYEVNT to administer the generation of an event *
36 C S PYSTAT to print cross-section and other information *
37 C S PYINRE to initialize treatment of resonances *
38 C S PYINBM to read in beam, target and frame choices *
39 C S PYINKI to initialize kinematics of incoming particles *
40 C S PYINPR to set up the selection of included processes *
41 C S PYXTOT to give total, elastic and diffractive cross-sect. *
42 C S PYMAXI to find differential cross-section maxima *
43 C S PYPILE to select multiplicity of pileup events *
44 C S PYSAVE to save alternatives for gamma-p and gamma-gamma *
45 C S PYRAND to select subprocess and kinematics for event *
46 C S PYSCAT to set up kinematics and colour flow of event *
47 C S PYSSPA to simulate initial state spacelike showers *
48 C S PYRESD to perform resonance decays *
49 C S PYMULT to generate multiple interactions *
50 C S PYREMN to add on target remnants *
51 C S PYDIFF to set up kinematics for diffractive events *
52 C S PYDOCU to compute cross-sections and handle documentation *
53 C S PYFRAM to perform boosts between different frames *
54 C S PYWIDT to calculate full and partial widths of resonances *
55 C S PYOFSH to calculate partial width into off-shell channels *
56 C S PYKLIM to calculate borders of allowed kinematical region *
57 C S PYKMAP to construct value of kinematical variable *
58 C S PYSIGH to calculate differential cross-sections *
59 C S PYSTFU to evaluate structure functions *
60 C S PYSTFL to evaluate structure functions at low x and Q^2 *
61 C S PYSTEL to evaluate electron structure function *
62 C S PYSTGA to evaluate photon structure function (generic) *
63 C S PYGGAM to evaluate photon structure function (SaS sets) *
64 C S PYGVMD to evaluate VMD part of photon structure functions *
65 C S PYGANO to evaluate anomalous part of photon str. func. *
66 C S PYGBEH to evaluate Bethe-Heitler part of photon str. func. *
67 C S PYGDIR to evaluate direct contribution to photon str. func. *
68 C S PYSTPI to evaluate pion structure function *
69 C S PYSTPR to evaluate proton structure function *
70 C F PYCTQ2 to evaluate the CTEQ 2 proton structure function *
71 C F PYHFTH to evaluate threshold factor for heavy flavour *
72 C S PYSPLI to find flavours left in hadron when one removed *
73 C F PYGAMM to evaluate ordinary Gamma function Gamma(x) *
74 C S PYWAUX to evaluate auxiliary functions W1(s) and W2(s) *
75 C S PYI3AU to evaluate auxiliary function I3(s,t,u,v) *
76 C F PYSPEN to evaluate Spence (dilogarithm) function Sp(x) *
77 C S PYQQBH to evaluate matrix element for g + g -> Q + Q~ + H *
78 C S PYTEST to test the proper functioning of the package *
79 C B PYDATA to contain all default values *
80 C S PYKCUT to provide dummy routine for user kinematical cuts *
81 C S PYEVWT to provide dummy routine for weighting events *
82 C S PYUPIN to initialize a user process *
83 C S PYUPEV to generate a user process event (dummy routine) *
84 C S PDFSET dummy routine to be removed when using PDFLIB *
85 C S STRUCTM dummy routine to be removed when using PDFLIB *
87 C*********************************************************************
89 SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN)
91 C...Initializes the generation procedure; finds maxima of the
92 C...differential cross-sections to be used for weighting.
93 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
94 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
95 COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
96 COMMON/LUDAT4/CHAF(500)
98 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
99 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
100 COMMON/PYINT1/MINT(400),VINT(400)
101 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
102 COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
103 COMMON/PYINT9/DXSEC(0:200)
104 DOUBLE PRECISION DXSEC
105 SAVE /LUDAT1/,/LUDAT2/,/LUDAT3/,/LUDAT4/
106 SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT5/,/PYINT9/
107 DIMENSION ALAMIN(20),NFIN(20)
108 CHARACTER*(*) FRAME,BEAM,TARGET
109 CHARACTER CHFRAM*8,CHBEAM*8,CHTARG*8,CHLH(2)*6
111 C...Interface to PDFLIB.
112 COMMON/W50512/QCDL4,QCDL5
114 DOUBLE PRECISION VALUE(20),QCDL4,QCDL5
115 CHARACTER*20 PARM(20)
116 DATA VALUE/20*0D0/,PARM/20*' '/
118 C...Data:Lambda and n_f values for structure functions; months.
119 DATA ALAMIN/0.20,0.29,0.20,0.40,0.213,0.208,0.208,0.322,
120 &0.190,0.235,10*0.2/,NFIN/20*4/
121 DATA CHLH/'lepton','hadron'/
123 C...Reset MINT and VINT arrays. Write headers.
128 IF(MSTU(12).GE.1) CALL LULIST(0)
129 IF(MSTP(122).GE.1) WRITE(MSTU(11),5100)
131 C...Maximum 4 generations; set maximum number of allowed flavours.
132 MSTP(1)=MIN(4,MSTP(1))
133 MSTU(114)=MIN(MSTU(114),2*MSTP(1))
134 MSTP(58)=MIN(MSTP(58),2*MSTP(1))
136 C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton.
140 IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN
145 IF(MDME(IDC,1).EQ.1.OR.MDME(IDC,1).EQ.IPM) VINT(180+I)=
146 & VINT(180+I)+VCKM((IA+1)/2,(IB+1)/2)
148 ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN
153 C...Initialize structure functions: PDFLIB.
154 IF(MSTP(52).EQ.2) THEN
158 VALUE(2)=MSTP(51)/1000
160 VALUE(3)=MOD(MSTP(51),1000)
163 CALL PDFSET(PARM,VALUE)
164 MINT(93)=1000000+MSTP(51)
167 C...Choose Lambda value to use in alpha-strong.
169 IF(MSTP(3).GE.2) THEN
172 IF(MSTP(52).EQ.1.AND.MSTP(51).GE.1.AND.MSTP(51).LE.10) THEN
173 ALAM=ALAMIN(MSTP(51))
175 ELSEIF(MSTP(52).EQ.2) THEN
184 IF(MSTP(3).EQ.3) PARJ(81)=ALAM
187 C...Initialize widths and partial widths for resonances.
190 C...Identify beam and target particles and frame of process.
194 CALL PYINBM(CHFRAM,CHBEAM,CHTARG,WIN)
195 IF(MINT(65).EQ.1) GOTO 170
197 C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives.
198 C...For e-gamma allow 2 alternatives.
201 IF(MSTP(14).EQ.10.AND.(MSEL.EQ.1.OR.MSEL.EQ.2)) THEN
202 IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
203 & (IABS(MINT(11)).GE.28.OR.IABS(MINT(12)).GE.28)) MINT(121)=3
204 IF(MINT(11).EQ.22.AND.MINT(12).EQ.22) MINT(121)=6
205 IF((MINT(11).EQ.22.OR.MINT(12).EQ.22).AND.
206 & (IABS(MINT(11)).EQ.11.OR.IABS(MINT(12)).EQ.11)) MINT(121)=2
209 C...Set up kinematics of process.
212 C...Loop over gamma-p or gamma-gamma alternatives.
213 DO 160 IGA=1,MINT(121)
216 C...Select partonic subprocesses to be included in the simulation.
219 C...Count number of subprocesses on.
222 IF(MINT(50).EQ.0.AND.ISUB.GE.91.AND.ISUB.LE.96.AND.
223 &MSUB(ISUB).EQ.1) THEN
224 WRITE(MSTU(11),5200) ISUB,CHLH(MINT(41)),CHLH(MINT(42))
226 ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN
227 WRITE(MSTU(11),5300) ISUB
229 ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN
230 WRITE(MSTU(11),5400) ISUB
232 ELSEIF(MSUB(ISUB).EQ.1) THEN
236 IF(MINT(48).EQ.0) THEN
240 MINT(49)=MINT(48)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94)
242 C...Reset variables for cross-section calculation.
251 C...Find parametrized total cross-sections.
254 C...Maxima of differential cross-sections.
255 IF(MSTP(121).LE.1) CALL PYMAXI
257 C...Initialize possibility of pileup events.
258 IF(MINT(121).GT.1) MSTP(131)=0
259 IF(MSTP(131).NE.0) CALL PYPILE(1)
261 C...Initialize multiple interactions with variable impact parameter.
262 IF(MINT(50).EQ.1.AND.(MINT(49).NE.0.OR.MSTP(131).NE.0).AND.
263 &MSTP(82).GE.2) CALL PYMULT(1)
265 C...Save results for gamma-p and gamma-gamma alternatives.
266 IF(MINT(121).GT.1) CALL PYSAVE(1,IGA)
269 C...Initialization finished.
270 170 IF(MSTP(122).GE.1) WRITE(MSTU(11),5600)
272 C...Formats for initialization information.
273 5100 FORMAT('1',18('*'),1X,'PYINIT: initialization of PYTHIA ',
274 &'routines',1X,17('*'))
275 5200 FORMAT(1X,'Error: process number ',I3,' not meaningful for ',A6,
276 &'-',A6,' interactions.'/1X,'Execution stopped!')
277 5300 FORMAT(1X,'Error: requested subprocess',I4,' not implemented.'/
278 &1X,'Execution stopped!')
279 5400 FORMAT(1X,'Error: requested subprocess',I4,' not existing.'/
280 &1X,'Execution stopped!')
281 5500 FORMAT(1X,'Error: no subprocess switched on.'/
282 &1X,'Execution stopped.')
283 5600 FORMAT(/1X,22('*'),1X,'PYINIT: initialization completed',1X,