]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | C********************************************************************* |
2 | C********************************************************************* | |
3 | C* ** | |
4 | C* December 1993 ** | |
5 | C* ** | |
6 | C* The Lund Monte Carlo for Hadronic Processes ** | |
7 | C* ** | |
8 | C* PYTHIA version 5.7 ** | |
9 | C* ** | |
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 ** | |
16 | C* ** | |
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 ** | |
22 | C* ** | |
23 | C* The latest program version and documentation is found on WWW ** | |
24 | C* http://thep.lu.se/tf2/staff/torbjorn/Welcome.html ** | |
25 | C* ** | |
26 | C* Copyright Torbjorn Sjostrand and CERN, Geneva 1993 ** | |
27 | C* ** | |
28 | C********************************************************************* | |
29 | C********************************************************************* | |
30 | C * | |
31 | C List of subprograms in order of appearance, with main purpose * | |
32 | C (S = subroutine, F = function, B = block data) * | |
33 | C * | |
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 * | |
86 | C * | |
87 | C********************************************************************* | |
88 | ||
89 | SUBROUTINE PYINIT(FRAME,BEAM,TARGET,WIN) | |
90 | ||
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) | |
97 | CHARACTER CHAF*8 | |
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 | |
110 | ||
111 | C...Interface to PDFLIB. | |
112 | COMMON/W50512/QCDL4,QCDL5 | |
113 | SAVE /W50512/ | |
114 | DOUBLE PRECISION VALUE(20),QCDL4,QCDL5 | |
115 | CHARACTER*20 PARM(20) | |
116 | DATA VALUE/20*0D0/,PARM/20*' '/ | |
117 | ||
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'/ | |
122 | ||
123 | C...Reset MINT and VINT arrays. Write headers. | |
124 | DO 100 J=1,400 | |
125 | MINT(J)=0 | |
126 | VINT(J)=0. | |
127 | 100 CONTINUE | |
128 | IF(MSTU(12).GE.1) CALL LULIST(0) | |
129 | IF(MSTP(122).GE.1) WRITE(MSTU(11),5100) | |
130 | ||
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)) | |
135 | ||
136 | C...Sum up Cabibbo-Kobayashi-Maskawa factors for each quark/lepton. | |
137 | DO 120 I=-20,20 | |
138 | VINT(180+I)=0. | |
139 | IA=IABS(I) | |
140 | IF(IA.GE.1.AND.IA.LE.2*MSTP(1)) THEN | |
141 | DO 110 J=1,MSTP(1) | |
142 | IB=2*J-1+MOD(IA,2) | |
143 | IPM=(5-ISIGN(1,I))/2 | |
144 | IDC=J+MDCY(IA,2)+2 | |
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) | |
147 | 110 CONTINUE | |
148 | ELSEIF(IA.GE.11.AND.IA.LE.10+2*MSTP(1)) THEN | |
149 | VINT(180+I)=1. | |
150 | ENDIF | |
151 | 120 CONTINUE | |
152 | ||
153 | C...Initialize structure functions: PDFLIB. | |
154 | IF(MSTP(52).EQ.2) THEN | |
155 | PARM(1)='NPTYPE' | |
156 | VALUE(1)=1 | |
157 | PARM(2)='NGROUP' | |
158 | VALUE(2)=MSTP(51)/1000 | |
159 | PARM(3)='NSET' | |
160 | VALUE(3)=MOD(MSTP(51),1000) | |
161 | PARM(4)='TMAS' | |
162 | VALUE(4)=PMAS(6,1) | |
163 | CALL PDFSET(PARM,VALUE) | |
164 | MINT(93)=1000000+MSTP(51) | |
165 | ENDIF | |
166 | ||
167 | C...Choose Lambda value to use in alpha-strong. | |
168 | MSTU(111)=MSTP(2) | |
169 | IF(MSTP(3).GE.2) THEN | |
170 | ALAM=0.2 | |
171 | NF=4 | |
172 | IF(MSTP(52).EQ.1.AND.MSTP(51).GE.1.AND.MSTP(51).LE.10) THEN | |
173 | ALAM=ALAMIN(MSTP(51)) | |
174 | NF=NFIN(MSTP(51)) | |
175 | ELSEIF(MSTP(52).EQ.2) THEN | |
176 | ALAM=QCDL4 | |
177 | NF=4 | |
178 | ENDIF | |
179 | PARP(1)=ALAM | |
180 | PARP(61)=ALAM | |
181 | PARP(72)=ALAM | |
182 | PARU(112)=ALAM | |
183 | MSTU(112)=NF | |
184 | IF(MSTP(3).EQ.3) PARJ(81)=ALAM | |
185 | ENDIF | |
186 | ||
187 | C...Initialize widths and partial widths for resonances. | |
188 | CALL PYINRE | |
189 | ||
190 | C...Identify beam and target particles and frame of process. | |
191 | CHFRAM=FRAME//' ' | |
192 | CHBEAM=BEAM//' ' | |
193 | CHTARG=TARGET//' ' | |
194 | CALL PYINBM(CHFRAM,CHBEAM,CHTARG,WIN) | |
195 | IF(MINT(65).EQ.1) GOTO 170 | |
196 | ||
197 | C...For gamma-p or gamma-gamma allow many (3 or 6) alternatives. | |
198 | C...For e-gamma allow 2 alternatives. | |
199 | MINT(121)=1 | |
200 | MINT(123)=MSTP(14) | |
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 | |
207 | ENDIF | |
208 | ||
209 | C...Set up kinematics of process. | |
210 | CALL PYINKI(0) | |
211 | ||
212 | C...Loop over gamma-p or gamma-gamma alternatives. | |
213 | DO 160 IGA=1,MINT(121) | |
214 | MINT(122)=IGA | |
215 | ||
216 | C...Select partonic subprocesses to be included in the simulation. | |
217 | CALL PYINPR | |
218 | ||
219 | C...Count number of subprocesses on. | |
220 | MINT(48)=0 | |
221 | DO 130 ISUB=1,200 | |
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)) | |
225 | STOP | |
226 | ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).EQ.-1) THEN | |
227 | WRITE(MSTU(11),5300) ISUB | |
228 | STOP | |
229 | ELSEIF(MSUB(ISUB).EQ.1.AND.ISET(ISUB).LE.-2) THEN | |
230 | WRITE(MSTU(11),5400) ISUB | |
231 | STOP | |
232 | ELSEIF(MSUB(ISUB).EQ.1) THEN | |
233 | MINT(48)=MINT(48)+1 | |
234 | ENDIF | |
235 | 130 CONTINUE | |
236 | IF(MINT(48).EQ.0) THEN | |
237 | WRITE(MSTU(11),5500) | |
238 | STOP | |
239 | ENDIF | |
240 | MINT(49)=MINT(48)-MSUB(91)-MSUB(92)-MSUB(93)-MSUB(94) | |
241 | ||
242 | C...Reset variables for cross-section calculation. | |
243 | DO 150 I=0,200 | |
244 | DO 140 J=1,3 | |
245 | NGEN(I,J)=0 | |
246 | XSEC(I,J)=0. | |
247 | 140 CONTINUE | |
248 | DXSEC(I)=0D0 | |
249 | 150 CONTINUE | |
250 | ||
251 | C...Find parametrized total cross-sections. | |
252 | CALL PYXTOT | |
253 | ||
254 | C...Maxima of differential cross-sections. | |
255 | IF(MSTP(121).LE.1) CALL PYMAXI | |
256 | ||
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) | |
260 | ||
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) | |
264 | ||
265 | C...Save results for gamma-p and gamma-gamma alternatives. | |
266 | IF(MINT(121).GT.1) CALL PYSAVE(1,IGA) | |
267 | 160 CONTINUE | |
268 | ||
269 | C...Initialization finished. | |
270 | 170 IF(MSTP(122).GE.1) WRITE(MSTU(11),5600) | |
271 | ||
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, | |
284 | &22('*')) | |
285 | ||
286 | RETURN | |
287 | END |