Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyrand.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyrand.f
1 
2 C*********************************************************************
3 
4 C...PYRAND
5 C...Generates quantities characterizing the high-pT scattering at the
6 C...parton level according to the matrix elements. Chooses incoming,
7 C...reacting partons, their momentum fractions and one of the possible
8 C...subprocesses.
9 
10  SUBROUTINE pyrand
11 
12 C...Double precision and integer declarations.
13  IMPLICIT DOUBLE PRECISION(a-h, o-z)
14  IMPLICIT INTEGER(i-n)
15  INTEGER pyk,pychge,pycomp
16 C...Parameter statement to help give large particle numbers.
17  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
18  &kexcit=4000000,kdimen=5000000)
19 
20 C...User process initialization and event commonblocks.
21  INTEGER maxpup
22  parameter(maxpup=100)
23  INTEGER idbmup,pdfgup,pdfsup,idwtup,nprup,lprup
24  DOUBLE PRECISION ebmup,xsecup,xerrup,xmaxup
25  common/heprup/idbmup(2),ebmup(2),pdfgup(2),pdfsup(2),
26  &idwtup,nprup,xsecup(maxpup),xerrup(maxpup),xmaxup(maxpup),
27  &lprup(maxpup)
28  INTEGER maxnup
29  parameter(maxnup=500)
30  INTEGER nup,idprup,idup,istup,mothup,icolup
31  DOUBLE PRECISION xwgtup,scalup,aqedup,aqcdup,pup,vtimup,spinup
32  common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
33  &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
34  &vtimup(maxnup),spinup(maxnup)
35  SAVE /heprup/,/hepeup/
36 
37 C...Commonblocks.
38  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
39  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
40  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
41  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
42  common/pypars/mstp(200),parp(200),msti(200),pari(200)
43  common/pyint1/mint(400),vint(400)
44  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
45  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
46  common/pyint4/mwid(500),wids(500,5)
47  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
48  common/pyint7/sigt(0:6,0:6,0:5)
49  common/pymssm/imss(0:99),rmss(0:99)
50  SAVE /pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,/pyint1/,
51  &/pyint2/,/pyint3/,/pyint4/,/pyint5/,/pyint7/,/pymssm/
52 C...Local arrays.
53  dimension xpq(-25:25),pmm(2),pdif(4),bhad(4),pmmn(2)
54 
55 C...Parameters and data used in elastic/diffractive treatment.
56  DATA eps/0.0808d0/, alp/0.25d0/, cres/2d0/, pmrc/1.062d0/,
57  &smp/0.880d0/, bhad/2.3d0,1.4d0,1.4d0,0.23d0/
58 
59 C...Initial values, specifically for (first) semihard interaction.
60  mint(10)=0
61  mint(17)=0
62  mint(18)=0
63  vint(143)=1d0
64  vint(144)=1d0
65  vint(157)=0d0
66  vint(158)=0d0
67  mfail=0
68  IF(mstp(171).EQ.1.AND.mstp(172).EQ.2) mfail=1
69  isub=0
70  istsb=0
71  loop=0
72  100 loop=loop+1
73  mint(51)=0
74  mint(143)=1
75  vint(97)=1d0
76 
77 C...Start by assuming incoming photon is entering subprocess.
78  IF(mint(11).EQ.22) THEN
79  mint(15)=22
80  vint(307)=vint(3)**2
81  ENDIF
82  IF(mint(12).EQ.22) THEN
83  mint(16)=22
84  vint(308)=vint(4)**2
85  ENDIF
86  mint(103)=mint(11)
87  mint(104)=mint(12)
88 
89 C...Choice of process type - first event of pileup.
90  inmult=0
91  IF(mint(82).EQ.1.AND.isub.GE.91.AND.isub.LE.96) THEN
92  ELSEIF(mint(82).EQ.1) THEN
93 
94 C...For gamma-p or gamma-gamma first pick between alternatives.
95  iga=0
96  IF(mint(121).GT.1) CALL pysave(4,iga)
97  mint(122)=iga
98 
99 C...For real gamma + gamma with different nature, flip at random.
100  IF(mint(11).EQ.22.AND.mint(12).EQ.22.AND.mint(123).GE.4.AND.
101  & mstp(14).LE.10.AND.pyr(0).GT.0.5d0) THEN
102  mintsv=mint(41)
103  mint(41)=mint(42)
104  mint(42)=mintsv
105  mintsv=mint(45)
106  mint(45)=mint(46)
107  mint(46)=mintsv
108  mintsv=mint(107)
109  mint(107)=mint(108)
110  mint(108)=mintsv
111  IF(mint(47).EQ.2.OR.mint(47).EQ.3) mint(47)=5-mint(47)
112  ENDIF
113 
114 C...Pick process type, possibly by user process machinery.
115 C...(If the latter, also event will be picked here.)
116  IF(mint(111).GE.11.AND.iabs(idwtup).EQ.2.AND.loop.GE.2) THEN
117  CALL upevnt
118  CALL pyupre
119  ELSEIF(mint(111).GE.11.AND.iabs(idwtup).GE.3) THEN
120  CALL upevnt
121  CALL pyupre
122  isub=0
123  110 isub=isub+1
124  IF((iset(isub).NE.11.OR.kfpr(isub,2).NE.idprup).AND.
125  & isub.LT.500) goto 110
126  ELSE
127  rsub=xsec(0,1)*pyr(0)
128  DO 120 i=1,500
129  IF(msub(i).NE.1.OR.i.EQ.96) goto 120
130  isub=i
131  rsub=rsub-xsec(i,1)
132  IF(rsub.LE.0d0) goto 130
133  120 CONTINUE
134  130 IF(isub.EQ.95) isub=96
135  IF(isub.EQ.96) inmult=1
136  IF(iset(isub).EQ.11) THEN
137  idprup=kfpr(isub,2)
138  CALL upevnt
139  CALL pyupre
140  ENDIF
141  ENDIF
142 
143 C...Choice of inclusive process type - pileup events.
144  ELSEIF(mint(82).GE.2.AND.isub.EQ.0) THEN
145  rsub=vint(131)*pyr(0)
146  isub=96
147  IF(rsub.GT.sigt(0,0,5)) isub=94
148  IF(rsub.GT.sigt(0,0,5)+sigt(0,0,4)) isub=93
149  IF(rsub.GT.sigt(0,0,5)+sigt(0,0,4)+sigt(0,0,3)) isub=92
150  IF(rsub.GT.sigt(0,0,5)+sigt(0,0,4)+sigt(0,0,3)+sigt(0,0,2))
151  & isub=91
152  IF(isub.EQ.96) inmult=1
153  ENDIF
154 
155 C...Choice of photon energy and flux factor inside lepton.
156  IF(mint(141).NE.0.OR.mint(142).NE.0) THEN
157  IF (mstp(199).EQ.1) THEN
158  CALL pygaga(5,wtgaga)
159  ELSE
160  CALL pygaga(3,wtgaga)
161  ENDIF
162  IF(isub.GE.131.AND.isub.LE.140) THEN
163  ckin(3)=max(vint(285),vint(154))
164  ckin(1)=2d0*ckin(3)
165  ENDIF
166 C...When necessary set direct/resolved photon by hand.
167  ELSEIF(mint(15).EQ.22.OR.mint(16).EQ.22) THEN
168  IF(mint(15).EQ.22.AND.mint(41).EQ.2) mint(15)=0
169  IF(mint(16).EQ.22.AND.mint(42).EQ.2) mint(16)=0
170  ENDIF
171 
172 C...Restrict direct*resolved processes to pTmin >= Q,
173 C...to avoid doublecounting with DIS.
174  IF(mstp(18).EQ.3.AND.isub.GE.131.AND.isub.LE.136) THEN
175  IF(mint(15).EQ.22) THEN
176  ckin(3)=max(vint(285),vint(154),abs(vint(3)))
177  ELSE
178  ckin(3)=max(vint(285),vint(154),abs(vint(4)))
179  ENDIF
180  ckin(1)=2d0*ckin(3)
181  ENDIF
182 
183 C...Set up for multiple interactions (may include impact parameter).
184  IF(inmult.EQ.1) THEN
185  IF(mint(35).LE.1) CALL pymult(2)
186  IF(mint(35).GE.2) CALL pymign(2)
187  ENDIF
188 
189 C...Loopback point for minimum bias in photon physics.
190  loop2=0
191  140 loop2=loop2+1
192  IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)+mint(143)
193  IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)+mint(143)
194  IF(isub.EQ.96.AND.loop2.EQ.1.AND.mint(82).EQ.1)
195  &ngen(97,1)=ngen(97,1)+mint(143)
196  mint(1)=isub
197  istsb=iset(isub)
198 
199 C...Random choice of flavour for some SUSY processes.
200  IF(isub.GE.201.AND.isub.LE.301) THEN
201 C...~e_L ~nu_e or ~mu_L ~nu_mu.
202  IF(isub.EQ.210) THEN
203  kfpr(isub,1)=ksusy1+11+2*int(0.5d0+pyr(0))
204  kfpr(isub,2)=kfpr(isub,1)+1
205 C...~nu_e ~nu_e(bar) or ~nu_mu ~nu_mu(bar).
206  ELSEIF(isub.EQ.213) THEN
207  kfpr(isub,1)=ksusy1+12+2*int(0.5d0+pyr(0))
208  kfpr(isub,2)=kfpr(isub,1)
209 C...~q ~chi/~g; ~q = ~d, ~u, ~s, ~c or ~b.
210  ELSEIF(isub.GE.246.AND.isub.LE.259.AND.isub.NE.255.AND.
211  & isub.NE.257) THEN
212  IF(isub.GE.258) THEN
213  rkf=4d0
214  ELSE
215  rkf=5d0
216  ENDIF
217  IF(mod(isub,2).EQ.0) THEN
218  kfpr(isub,1)=ksusy1+1+int(rkf*pyr(0))
219  ELSE
220  kfpr(isub,1)=ksusy2+1+int(rkf*pyr(0))
221  ENDIF
222 C...~q1 ~q2; ~q = ~d, ~u, ~s, or ~c.
223  ELSEIF(isub.GE.271.AND.isub.LE.276) THEN
224  IF(isub.EQ.271.OR.isub.EQ.274) THEN
225  ksu1=ksusy1
226  ksu2=ksusy1
227  ELSEIF(isub.EQ.272.OR.isub.EQ.275) THEN
228  ksu1=ksusy2
229  ksu2=ksusy2
230  ELSEIF(pyr(0).LT.0.5d0) THEN
231  ksu1=ksusy1
232  ksu2=ksusy2
233  ELSE
234  ksu1=ksusy2
235  ksu2=ksusy1
236  ENDIF
237  kfpr(isub,1)=ksu1+1+int(4d0*pyr(0))
238  kfpr(isub,2)=ksu2+1+int(4d0*pyr(0))
239 C...~q ~q(bar); ~q = ~d, ~u, ~s, or ~c.
240  ELSEIF(isub.EQ.277.OR.isub.EQ.279) THEN
241  kfpr(isub,1)=ksusy1+1+int(4d0*pyr(0))
242  kfpr(isub,2)=kfpr(isub,1)
243  ELSEIF(isub.EQ.278.OR.isub.EQ.280) THEN
244  kfpr(isub,1)=ksusy2+1+int(4d0*pyr(0))
245  kfpr(isub,2)=kfpr(isub,1)
246 C...~q1 ~q2; ~q = ~d, ~u, ~s, or ~c.
247  ELSEIF(isub.GE.281.AND.isub.LE.286) THEN
248  IF(isub.EQ.281.OR.isub.EQ.284) THEN
249  ksu1=ksusy1
250  ksu2=ksusy1
251  ELSEIF(isub.EQ.282.OR.isub.EQ.285) THEN
252  ksu1=ksusy2
253  ksu2=ksusy2
254  ELSEIF(pyr(0).LT.0.5d0) THEN
255  ksu1=ksusy1
256  ksu2=ksusy2
257  ELSE
258  ksu1=ksusy2
259  ksu2=ksusy1
260  ENDIF
261  IF(isub.EQ.281.OR.isub.LE.283) THEN
262  rkf=5d0
263  ELSE
264  rkf=4d0
265  ENDIF
266  kfpr(isub,2)=ksu2+1+int(rkf*pyr(0))
267  ENDIF
268  ENDIF
269 
270 C...Find resonances (explicit or implicit in cross-section).
271  mint(72)=0
272  kfr1=0
273  IF(istsb.EQ.1.OR.istsb.EQ.3.OR.istsb.EQ.5) THEN
274  kfr1=kfpr(isub,1)
275  ELSEIF(isub.EQ.24.OR.isub.EQ.25.OR.isub.EQ.110.OR.isub.EQ.165.OR.
276  & isub.EQ.171.OR.isub.EQ.176) THEN
277  kfr1=23
278  ELSEIF(isub.EQ.23.OR.isub.EQ.26.OR.isub.EQ.166.OR.isub.EQ.172.OR.
279  & isub.EQ.177) THEN
280  kfr1=24
281  ELSEIF(isub.GE.71.AND.isub.LE.77) THEN
282  kfr1=25
283  IF(mstp(46).EQ.5) THEN
284  kfr1=89
285  pmas(89,1)=parp(45)
286  pmas(89,2)=parp(45)**3/(96d0*paru(1)*parp(47)**2)
287  ENDIF
288  ELSEIF(isub.EQ.194) THEN
289  kfr1=ktechn+113
290  ELSEIF(isub.EQ.195) THEN
291  kfr1=ktechn+213
292  ELSEIF(isub.GE.361.AND.isub.LE.368) THEN
293  kfr1=ktechn+113
294  ELSEIF(isub.GE.370.AND.isub.LE.377) THEN
295  kfr1=ktechn+213
296  ENDIF
297  ckmx=ckin(2)
298  IF(ckmx.LE.0d0) ckmx=vint(1)
299  kcr1=pycomp(kfr1)
300  IF(kfr1.NE.0) THEN
301  IF(ckin(1).GT.pmas(kcr1,1)+20d0*pmas(kcr1,2).OR.
302  & ckmx.LT.pmas(kcr1,1)-20d0*pmas(kcr1,2)) kfr1=0
303  ENDIF
304  IF(kfr1.NE.0) THEN
305  taur1=pmas(kcr1,1)**2/vint(2)
306  IF(kfr1.EQ.ktechn+113) THEN
307  CALL pytecm(s1,s2)
308  taur1=s1/vint(2)
309  ENDIF
310  gamr1=pmas(kcr1,1)*pmas(kcr1,2)/vint(2)
311  mint(72)=1
312  mint(73)=kfr1
313  vint(73)=taur1
314  vint(74)=gamr1
315  ENDIF
316  IF(isub.EQ.141.OR.isub.EQ.194.OR.(isub.GE.364.AND.isub.LE.368))
317  $THEN
318  kfr2=23
319  IF(isub.EQ.194) THEN
320  kfr2=ktechn+223
321  ELSEIF(isub.GE.364.AND.isub.LE.368) THEN
322  kfr2=ktechn+223
323  ENDIF
324  kcr2=pycomp(kfr2)
325  taur2=pmas(kcr2,1)**2/vint(2)
326  IF(kfr2.EQ.ktechn+223) THEN
327  CALL pytecm(s1,s2)
328  taur2=s2/vint(2)
329  ENDIF
330  gamr2=pmas(kcr2,1)*pmas(kcr2,2)/vint(2)
331  IF(ckin(1).GT.pmas(kcr2,1)+20d0*pmas(kcr2,2).OR.
332  & ckmx.LT.pmas(kcr2,1)-20d0*pmas(kcr2,2)) kfr2=0
333  IF(kfr2.NE.0.AND.kfr1.NE.0) THEN
334  mint(72)=2
335  mint(74)=kfr2
336  vint(75)=taur2
337  vint(76)=gamr2
338  ELSEIF(kfr2.NE.0) THEN
339  kfr1=kfr2
340  taur1=taur2
341  gamr1=gamr2
342  mint(72)=1
343  mint(73)=kfr1
344  vint(73)=taur1
345  vint(74)=gamr1
346  ENDIF
347  ENDIF
348 
349 C...Find product masses and minimum pT of process,
350 C...optionally with broadening according to a truncated Breit-Wigner.
351  vint(63)=0d0
352  vint(64)=0d0
353  mint(71)=0
354  vint(71)=ckin(3)
355  IF(mint(82).GE.2) vint(71)=0d0
356  vint(80)=1d0
357  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
358  nbw=0
359  DO 160 i=1,2
360  pmmn(i)=0d0
361  IF(kfpr(isub,i).EQ.0) THEN
362  ELSEIF(mstp(42).LE.0.OR.pmas(pycomp(kfpr(isub,i)),2).LT.
363  & parp(41)) THEN
364  vint(62+i)=pmas(pycomp(kfpr(isub,i)),1)**2
365  ELSE
366  nbw=nbw+1
367 C...This prevents SUSY/t particles from becoming too light.
368  kflw=kfpr(isub,i)
369  IF(kflw/ksusy1.EQ.1.OR.kflw/ksusy1.EQ.2) THEN
370  kcw=pycomp(kflw)
371  pmmn(i)=pmas(kcw,1)
372  DO 150 idc=mdcy(kcw,2),mdcy(kcw,2)+mdcy(kcw,3)-1
373  IF(mdme(idc,1).GT.0.AND.brat(idc).GT.1e-4) THEN
374  pmsum=pmas(pycomp(kfdp(idc,1)),1)+
375  & pmas(pycomp(kfdp(idc,2)),1)
376  IF(kfdp(idc,3).NE.0) pmsum=pmsum+
377  & pmas(pycomp(kfdp(idc,3)),1)
378  pmmn(i)=min(pmmn(i),pmsum)
379  ENDIF
380  150 CONTINUE
381  ELSEIF(kflw.EQ.6) THEN
382  pmmn(i)=pmas(24,1)+pmas(5,1)
383  ENDIF
384  ENDIF
385  160 CONTINUE
386  IF(nbw.GE.1) THEN
387  ckin41=ckin(41)
388  ckin43=ckin(43)
389  ckin(41)=max(pmmn(1),ckin(41))
390  ckin(43)=max(pmmn(2),ckin(43))
391  CALL pyofsh(4,0,kfpr(isub,1),kfpr(isub,2),0d0,pqm3,pqm4)
392  ckin(41)=ckin41
393  ckin(43)=ckin43
394  IF(mint(51).EQ.1) THEN
395  IF(mint(121).GT.1) CALL pysave(2,iga)
396  IF(mfail.EQ.1) THEN
397  msti(61)=1
398  RETURN
399  ENDIF
400  goto 100
401  ENDIF
402  vint(63)=pqm3**2
403  vint(64)=pqm4**2
404  ENDIF
405  IF(min(vint(63),vint(64)).LT.ckin(6)**2) mint(71)=1
406  IF(mint(71).EQ.1) vint(71)=max(ckin(3),ckin(5))
407  ENDIF
408 
409 C...Prepare for additional variable choices in 2 -> 3.
410  IF(istsb.EQ.5) THEN
411  vint(201)=0d0
412  IF(kfpr(isub,2).GT.0) vint(201)=pmas(pycomp(kfpr(isub,2)),1)
413  vint(206)=vint(201)
414  IF(isub.EQ.401.OR.isub.EQ.402) vint(206)=pmas(5,1)
415  vint(204)=pmas(23,1)
416  IF(isub.EQ.124.OR.isub.EQ.174.OR.isub.EQ.179.OR.isub.EQ.351)
417  & vint(204)=pmas(24,1)
418  IF(isub.EQ.352) vint(204)=pmas(pycomp(9900024),1)
419  IF(isub.EQ.121.OR.isub.EQ.122.OR.isub.EQ.181.OR.isub.EQ.182.OR.
420  & isub.EQ.186.OR.isub.EQ.187.OR.isub.EQ.401.OR.isub.EQ.402)
421  & vint(204)=vint(201)
422  vint(209)=vint(204)
423  IF(isub.EQ.401.OR.isub.EQ.402) vint(209)=vint(206)
424  ENDIF
425 
426 C...Select incoming VDM particle (rho/omega/phi/J/psi).
427  IF(istsb.NE.0.AND.(mint(101).GE.2.OR.mint(102).GE.2).AND.
428  &(mint(123).EQ.2.OR.mint(123).EQ.3.OR.mint(123).EQ.7)) THEN
429  vrn=pyr(0)*sigt(0,0,5)
430  IF(mint(101).LE.1) THEN
431  i1mn=0
432  i1mx=0
433  ELSE
434  i1mn=1
435  i1mx=mint(101)
436  ENDIF
437  IF(mint(102).LE.1) THEN
438  i2mn=0
439  i2mx=0
440  ELSE
441  i2mn=1
442  i2mx=mint(102)
443  ENDIF
444  DO 180 i1=i1mn,i1mx
445  kfv1=110*i1+3
446  DO 170 i2=i2mn,i2mx
447  kfv2=110*i2+3
448  vrn=vrn-sigt(i1,i2,5)
449  IF(vrn.LE.0d0) goto 190
450  170 CONTINUE
451  180 CONTINUE
452  190 IF(mint(101).GE.2) mint(103)=kfv1
453  IF(mint(102).GE.2) mint(104)=kfv2
454  ENDIF
455 
456  IF(istsb.EQ.0) THEN
457 C...Elastic scattering or single or double diffractive scattering.
458 
459 C...Select incoming particle (rho/omega/phi/J/psi for VDM) and mass.
460  mint(103)=mint(11)
461  mint(104)=mint(12)
462  pmm(1)=vint(3)
463  pmm(2)=vint(4)
464  IF(mint(101).GE.2.OR.mint(102).GE.2) THEN
465  jj=isub-90
466  vrn=pyr(0)*sigt(0,0,jj)
467  IF(mint(101).LE.1) THEN
468  i1mn=0
469  i1mx=0
470  ELSE
471  i1mn=1
472  i1mx=mint(101)
473  ENDIF
474  IF(mint(102).LE.1) THEN
475  i2mn=0
476  i2mx=0
477  ELSE
478  i2mn=1
479  i2mx=mint(102)
480  ENDIF
481  DO 210 i1=i1mn,i1mx
482  kfv1=110*i1+3
483  DO 200 i2=i2mn,i2mx
484  kfv2=110*i2+3
485  vrn=vrn-sigt(i1,i2,jj)
486  IF(vrn.LE.0d0) goto 220
487  200 CONTINUE
488  210 CONTINUE
489  220 IF(mint(101).GE.2) THEN
490  mint(103)=kfv1
491  pmm(1)=pymass(kfv1)
492  ENDIF
493  IF(mint(102).GE.2) THEN
494  mint(104)=kfv2
495  pmm(2)=pymass(kfv2)
496  ENDIF
497  ENDIF
498  vint(67)=pmm(1)
499  vint(68)=pmm(2)
500 
501 C...Select mass for GVMD states (rejecting previous assignment).
502  q0s=4d0*parp(15)**2
503  q1s=4d0*vint(154)**2
504  loop3=0
505  230 loop3=loop3+1
506  DO 240 jt=1,2
507  IF(mint(106+jt).EQ.3) THEN
508  ps=vint(2+jt)**2
509  pmm(jt)=(q0s+ps)*(q1s+ps)/
510  & (q0s+pyr(0)*(q1s-q0s)+ps)-ps
511  IF(mint(102+jt).GE.333) pmm(jt)=pmm(jt)-
512  & pmas(pycomp(113),1)+pmas(pycomp(mint(102+jt)),1)
513  ENDIF
514  240 CONTINUE
515  IF(pmm(1)+pmm(2)+parp(104).GE.vint(1)) THEN
516  IF(loop3.LT.100.AND.(mint(107).EQ.3.OR.mint(108).EQ.3))
517  & goto 230
518  goto 100
519  ENDIF
520 
521 C...Side/sides of diffractive system.
522  mint(17)=0
523  mint(18)=0
524  IF(isub.EQ.92.OR.isub.EQ.94) mint(17)=1
525  IF(isub.EQ.93.OR.isub.EQ.94) mint(18)=1
526 
527 C...Find masses of particles and minimal masses of diffractive states.
528  DO 250 jt=1,2
529  pdif(jt)=pmm(jt)
530  vint(68+jt)=pdif(jt)
531  IF(mint(16+jt).EQ.1) pdif(jt)=pdif(jt)+parp(102)
532  250 CONTINUE
533  sh=vint(2)
534  if (mint(11).eq.22) then
535  sqm1 = -1.*(vint(3)**2)
536  else
537  sqm1=pmm(1)**2
538  endif
539  sqm2=pmm(2)**2
540  sqm3=pdif(1)**2
541  sqm4=pdif(2)**2
542  smres1=(pmm(1)+pmrc)**2
543  smres2=(pmm(2)+pmrc)**2
544 
545 C...Find elastic slope and lower limit diffractive slope.
546  iha=max(2,iabs(mint(103))/110)
547  IF(iha.GE.5) iha=1
548  ihb=max(2,iabs(mint(104))/110)
549  IF(ihb.GE.5) ihb=1
550  IF(isub.EQ.91) THEN
551  IF(iha.eq.2) then
552  pmvirt=pmas(pycomp(113),1)
553  bmn=5.84/(1+(parp(165))*(vint(307)/(pmvirt**2))**parp(166))
554  & +4.5
555  ELSEIF(iha.eq.3) then
556  pmvirt=pmas(pycomp(333),1)
557  bmn=5.84/(1+(parp(165))*(vint(307)/(pmvirt**2))**parp(166))
558  & +4.5
559  elseif(iha.eq.4) then
560  pmvirt=pmas(pycomp(443),1)
561  bmn=4.5
562  else
563  bmn=2d0*bhad(iha)+2d0*bhad(ihb)+4d0*sh**eps-4.2d0
564  endif
565  ELSEIF(isub.EQ.92) THEN
566  bmn=max(2d0,2d0*bhad(ihb))
567  ELSEIF(isub.EQ.93) THEN
568  bmn=max(2d0,2d0*bhad(iha))
569  ELSEIF(isub.EQ.94) THEN
570  bmn=2d0*alp*4d0
571  ENDIF
572 
573 C...Determine maximum possible t range and coefficient of generation.
574  sqla12=(sh-sqm1-sqm2)**2-4d0*sqm1*sqm2
575  sqla34=(sh-sqm3-sqm4)**2-4d0*sqm3*sqm4
576  tha=sh-(sqm1+sqm2+sqm3+sqm4)+(sqm1-sqm2)*(sqm3-sqm4)/sh
577  thb=sqrt(max(0d0,sqla12))*sqrt(max(0d0,sqla34))/sh
578  thc=(sqm3-sqm1)*(sqm4-sqm2)+(sqm1+sqm4-sqm2-sqm3)*
579  & (sqm1*sqm4-sqm2*sqm3)/sh
580  thl=-0.5d0*(tha+thb)
581  thu=thc/thl
582  thrnd=exp(max(-50d0,bmn*(thl-thu)))-1d0
583 
584 C...Select diffractive mass/masses according to dm^2/m^2.
585  loop3=0
586  260 loop3=loop3+1
587  DO 270 jt=1,2
588  IF(mint(16+jt).EQ.0) THEN
589  pdif(2+jt)=pdif(jt)
590  ELSE
591  pmmin=pdif(jt)
592  pmmax=max(vint(2+jt),vint(1)-pdif(3-jt))
593  pdif(2+jt)=pmmin*(pmmax/pmmin)**pyr(0)
594  ENDIF
595  270 CONTINUE
596  sqm3=pdif(3)**2
597  sqm4=pdif(4)**2
598 
599 C..Additional mass factors, including resonance enhancement.
600  IF(pdif(3)+pdif(4).GE.vint(1)) THEN
601  IF(loop3.LT.100) goto 260
602  goto 100
603  ENDIF
604  IF(isub.EQ.92) THEN
605  fsd=(1d0-sqm3/sh)*(1d0+cres*smres1/(smres1+sqm3))
606  IF(fsd.LT.pyr(0)*(1d0+cres)) goto 260
607  ELSEIF(isub.EQ.93) THEN
608  fsd=(1d0-sqm4/sh)*(1d0+cres*smres2/(smres2+sqm4))
609  IF(fsd.LT.pyr(0)*(1d0+cres)) goto 260
610  ELSEIF(isub.EQ.94) THEN
611  fdd=(1d0-(pdif(3)+pdif(4))**2/sh)*(sh*smp/
612  & (sh*smp+sqm3*sqm4))*(1d0+cres*smres1/(smres1+sqm3))*
613  & (1d0+cres*smres2/(smres2+sqm4))
614  IF(fdd.LT.pyr(0)*(1d0+cres)**2) goto 260
615  ENDIF
616 
617 C...Select t according to exp(Bmn*t) and correct to right slope.
618  th=thu+log(1d0+thrnd*pyr(0))/bmn
619 C if (PMVIRT.eq.PMAS(PYCOMP(443),1)) then
620 C write(*,*) BMN, PMVIRT, VINT(307), MINT(101), PMM(1), PMM(2),
621 C & SQM1, SQM2, "THA: ", THA, "THB: ", THB, "THC: ", THC,
622 C & "THL: ", THL, "THU: ", THU, "THRND: ", THRND,
623 C & "TH: ", TH
624 C write(*,*) "======================="
625 C endif
626  IF(isub.GE.92) THEN
627  IF(isub.EQ.92) THEN
628  badd=2d0*alp*log(sh/sqm3)
629  IF(bhad(ihb).LT.1d0) badd=max(0d0,badd+2d0*bhad(ihb)-2d0)
630  ELSEIF(isub.EQ.93) THEN
631  badd=2d0*alp*log(sh/sqm4)
632  IF(bhad(iha).LT.1d0) badd=max(0d0,badd+2d0*bhad(iha)-2d0)
633  ELSEIF(isub.EQ.94) THEN
634  badd=2d0*alp*(log(exp(4d0)+sh/(alp*sqm3*sqm4))-4d0)
635  ENDIF
636  IF(exp(max(-50d0,badd*(th-thu))).LT.pyr(0)) goto 260
637  ENDIF
638 
639 C...Check whether m^2 and t choices are consistent.
640  sqla34=(sh-sqm3-sqm4)**2-4d0*sqm3*sqm4
641  tha=sh-(sqm1+sqm2+sqm3+sqm4)+(sqm1-sqm2)*(sqm3-sqm4)/sh
642  thb=sqrt(max(0d0,sqla12))*sqrt(max(0d0,sqla34))/sh
643  IF(thb.LE.1d-8) goto 260
644  thc=(sqm3-sqm1)*(sqm4-sqm2)+(sqm1+sqm4-sqm2-sqm3)*
645  & (sqm1*sqm4-sqm2*sqm3)/sh
646  thlm=-0.5d0*(tha+thb)
647  thum=thc/thlm
648  IF(th.LT.thlm.OR.th.GT.thum) goto 260
649 
650 C...Information to output.
651  vint(21)=1d0
652  vint(22)=0d0
653  vint(23)=min(1d0,max(-1d0,(tha+2d0*th)/thb))
654  vint(45)=th
655  vint(59)=2d0*sqrt(max(0d0,-(thc+tha*th+th**2)))/thb
656  vint(63)=pdif(3)**2
657  vint(64)=pdif(4)**2
658  vint(283)=pmm(1)**2/4d0
659  vint(284)=pmm(2)**2/4d0
660 
661 C...Note: in the following, by In is meant the integral over the
662 C...quantity multiplying coefficient cn.
663 C...Choose tau according to h1(tau)/tau, where
664 C...h1(tau) = c1 + I1/I2*c2*1/tau + I1/I3*c3*1/(tau+tau_R) +
665 C...I1/I4*c4*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
666 C...I1/I5*c5*1/(tau+tau_R') +
667 C...I1/I6*c6*tau/((s*tau-m'^2)^2+(m'*Gamma')^2) +
668 C...I1/I7*c7*tau/(1.-tau), and
669 C...c1 + c2 + c3 + c4 + c5 + c6 + c7 = 1.
670  ELSEIF(istsb.GE.1.AND.istsb.LE.5) THEN
671  CALL pyklim(1)
672  IF(mint(51).NE.0) THEN
673  IF(mint(121).GT.1) CALL pysave(2,iga)
674  IF(mfail.EQ.1) THEN
675  msti(61)=1
676  RETURN
677  ENDIF
678  goto 100
679  ENDIF
680  rtau=pyr(0)
681  mtau=1
682  IF(rtau.GT.coef(isub,1)) mtau=2
683  IF(rtau.GT.coef(isub,1)+coef(isub,2)) mtau=3
684  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)) mtau=4
685  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4))
686  & mtau=5
687  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4)+
688  & coef(isub,5)) mtau=6
689  IF(rtau.GT.coef(isub,1)+coef(isub,2)+coef(isub,3)+coef(isub,4)+
690  & coef(isub,5)+coef(isub,6)) mtau=7
691  CALL pykmap(1,mtau,pyr(0))
692 
693 C...2 -> 3, 4 processes:
694 C...Choose tau' according to h4(tau,tau')/tau', where
695 C...h4(tau,tau') = c1 + I1/I2*c2*(1 - tau/tau')^3/tau' +
696 C...I1/I3*c3*1/(1 - tau'), and c1 + c2 + c3 = 1.
697  IF(istsb.GE.3.AND.istsb.LE.5) THEN
698  CALL pyklim(4)
699  IF(mint(51).NE.0) THEN
700  IF(mint(121).GT.1) CALL pysave(2,iga)
701  IF(mfail.EQ.1) THEN
702  msti(61)=1
703  RETURN
704  ENDIF
705  goto 100
706  ENDIF
707  rtaup=pyr(0)
708  mtaup=1
709  IF(rtaup.GT.coef(isub,18)) mtaup=2
710  IF(rtaup.GT.coef(isub,18)+coef(isub,19)) mtaup=3
711  CALL pykmap(4,mtaup,pyr(0))
712  ENDIF
713 
714 C...Choose y* according to h2(y*), where
715 C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
716 C...I0/I3*c3*1/cosh(y*) + I0/I4*c4*1/(1-exp(y*-y*max)) +
717 C...I0/I5*c5*1/(1-exp(-y*-y*min)), I0 = y*max-y*min,
718 C...and c1 + c2 + c3 + c4 + c5 = 1.
719  CALL pyklim(2)
720  IF(mint(51).NE.0) THEN
721  IF(mint(121).GT.1) CALL pysave(2,iga)
722  IF(mfail.EQ.1) THEN
723  msti(61)=1
724  RETURN
725  ENDIF
726  goto 100
727  ENDIF
728  ryst=pyr(0)
729  myst=1
730  IF(ryst.GT.coef(isub,8)) myst=2
731  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
732  IF(ryst.GT.coef(isub,8)+coef(isub,9)+coef(isub,10)) myst=4
733  IF(ryst.GT.coef(isub,8)+coef(isub,9)+coef(isub,10)+
734  & coef(isub,11)) myst=5
735  CALL pykmap(2,myst,pyr(0))
736 
737 C...2 -> 2 processes:
738 C...Choose cos(theta-hat) (cth) according to h3(cth), where
739 C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
740 C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
741 C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
742 C...and c0 + c1 + c2 + c3 + c4 = 1.
743  CALL pyklim(3)
744  IF(mint(51).NE.0) THEN
745  IF(mint(121).GT.1) CALL pysave(2,iga)
746  IF(mfail.EQ.1) THEN
747  msti(61)=1
748  RETURN
749  ENDIF
750  goto 100
751  ENDIF
752  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
753  rcth=pyr(0)
754  mcth=1
755  IF(rcth.GT.coef(isub,13)) mcth=2
756  IF(rcth.GT.coef(isub,13)+coef(isub,14)) mcth=3
757  IF(rcth.GT.coef(isub,13)+coef(isub,14)+coef(isub,15)) mcth=4
758  IF(rcth.GT.coef(isub,13)+coef(isub,14)+coef(isub,15)+
759  & coef(isub,16)) mcth=5
760  CALL pykmap(3,mcth,pyr(0))
761  ENDIF
762 
763 C...2 -> 3 : select pT1, phi1, pT2, phi2, y3 for 3 outgoing.
764  IF(istsb.EQ.5) THEN
765  CALL pykmap(5,0,0d0)
766  IF(mint(51).NE.0) THEN
767  IF(mint(121).GT.1) CALL pysave(2,iga)
768  IF(mfail.EQ.1) THEN
769  msti(61)=1
770  RETURN
771  ENDIF
772  goto 100
773  ENDIF
774  ENDIF
775 
776 C...DIS as f + gamma* -> f process: set dummy values.
777  ELSEIF(istsb.EQ.8) THEN
778  vint(21)=0.9d0
779  vint(22)=0d0
780  vint(23)=0d0
781  vint(47)=0d0
782  vint(48)=0d0
783 
784 C...Low-pT or multiple interactions (first semihard interaction).
785  ELSEIF(istsb.EQ.9) THEN
786  IF(mint(35).LE.1) CALL pymult(3)
787  IF(mint(35).GE.2) CALL pymign(3)
788  isub=mint(1)
789 
790 C...Study user-defined process: kinematics plus weight.
791  ELSEIF(istsb.EQ.11) THEN
792  IF(idwtup.GT.0.AND.xwgtup.LT.0d0) CALL
793  & pyerrm(26,'(PYRAND:) Negative XWGTUP for user process')
794  msti(51)=0
795  IF(nup.LE.0) THEN
796  mint(51)=2
797  msti(51)=1
798  IF(mint(82).EQ.1) THEN
799  ngen(0,1)=ngen(0,1)-1
800  ngen(isub,1)=ngen(isub,1)-1
801  ENDIF
802  IF(mint(121).GT.1) CALL pysave(2,iga)
803  RETURN
804  ENDIF
805 
806 C...Extract cross section event weight.
807  IF(iabs(idwtup).EQ.1.OR.iabs(idwtup).EQ.4) THEN
808  sigs=1d-9*xwgtup
809  ELSE
810  sigs=1d-9*xsecup(kfpr(isub,1))
811  ENDIF
812  IF(iabs(idwtup).GE.1.AND.iabs(idwtup).LE.3) THEN
813  vint(97)=sign(1d0,xwgtup)
814  ELSE
815  vint(97)=1d-9*xwgtup
816  ENDIF
817 
818 C...Construct 'trivial' kinematical variables needed.
819  kfl1=idup(1)
820  kfl2=idup(2)
821  vint(41)=pup(4,1)/ebmup(1)
822  vint(42)=pup(4,2)/ebmup(2)
823  vint(21)=vint(41)*vint(42)
824  vint(22)=0.5d0*log(vint(41)/vint(42))
825  vint(44)=vint(21)*vint(2)
826  vint(43)=sqrt(max(0d0,vint(44)))
827  vint(55)=scalup
828  IF(scalup.LE.0d0) vint(55)=vint(43)
829  vint(56)=vint(55)**2
830  vint(57)=aqedup
831  vint(58)=aqcdup
832 
833 C...Construct other kinematical variables needed (approximately).
834  vint(23)=0d0
835  vint(26)=vint(21)
836  vint(45)=-0.5d0*vint(44)
837  vint(46)=-0.5d0*vint(44)
838  vint(49)=vint(43)
839  vint(50)=vint(44)
840  vint(51)=vint(55)
841  vint(52)=vint(56)
842  vint(53)=vint(55)
843  vint(54)=vint(56)
844  vint(25)=0d0
845  vint(48)=0d0
846  IF(istup(1).NE.-1.OR.istup(2).NE.-1) CALL pyerrm(26,
847  & '(PYRAND:) unacceptable ISTUP code for incoming particles')
848  DO 280 iup=3,nup
849  IF(istup(iup).LT.1.OR.istup(iup).GT.3) CALL pyerrm(26,
850  & '(PYRAND:) unacceptable ISTUP code for particles')
851  IF(istup(iup).EQ.1) vint(25)=vint(25)+2d0*(pup(5,iup)**2+
852  & pup(1,iup)**2+pup(2,iup)**2)/vint(2)
853  IF(istup(iup).EQ.1) vint(48)=vint(48)+0.5d0*(pup(1,iup)**2+
854  & pup(2,iup)**2)
855  280 CONTINUE
856  vint(47)=sqrt(vint(48))
857  ENDIF
858 
859 C...Choose azimuthal angle.
860  vint(24)=0d0
861  IF(istsb.NE.11) vint(24)=paru(2)*pyr(0)
862 
863 C...Check against user cuts on kinematics at parton level.
864  mint(51)=0
865  IF((isub.LE.90.OR.isub.GT.100).AND.istsb.LE.10) CALL pyklim(0)
866  IF(mint(51).NE.0) THEN
867  IF(mint(121).GT.1) CALL pysave(2,iga)
868  IF(mfail.EQ.1) THEN
869  msti(61)=1
870  RETURN
871  ENDIF
872  goto 100
873  ENDIF
874  IF(mint(82).EQ.1.AND.mstp(141).GE.1.AND.istsb.LE.10) THEN
875  mcut=0
876  IF(msub(91)+msub(92)+msub(93)+msub(94)+msub(95).EQ.0)
877  & CALL pykcut(mcut)
878  IF(mcut.NE.0) THEN
879  IF(mint(121).GT.1) CALL pysave(2,iga)
880  IF(mfail.EQ.1) THEN
881  msti(61)=1
882  RETURN
883  ENDIF
884  goto 100
885  ENDIF
886  ENDIF
887 
888 C...Calculate differential cross-section for different subprocesses.
889  IF(istsb.LE.10) CALL pysigh(nchn,sigs)
890  sigsor=sigs
891  siglpt=sigt(0,0,5)*vint(315)*vint(316)
892 
893 C...Multiply cross section by lepton -> photon flux factor.
894  IF(mint(141).NE.0.OR.mint(142).NE.0) THEN
895  sigs=wtgaga*sigs
896  DO 290 ichn=1,nchn
897  sigh(ichn)=wtgaga*sigh(ichn)
898  290 CONTINUE
899  siglpt=wtgaga*siglpt
900  ENDIF
901 
902 C...Multiply cross-section by user-defined weights.
903  IF(mstp(173).EQ.1) THEN
904  sigs=parp(173)*sigs
905  DO 300 ichn=1,nchn
906  sigh(ichn)=parp(173)*sigh(ichn)
907  300 CONTINUE
908  siglpt=parp(173)*siglpt
909  ENDIF
910  wtxs=1d0
911  sigswt=sigs
912  vint(99)=1d0
913  vint(100)=1d0
914  IF(mint(82).EQ.1.AND.mstp(142).GE.1) THEN
915  IF(isub.NE.96.AND.msub(91)+msub(92)+msub(93)+msub(94)+
916  & msub(95).EQ.0) CALL pyevwt(wtxs)
917  sigswt=wtxs*sigs
918  vint(99)=wtxs
919  IF(mstp(142).EQ.1) vint(100)=1d0/wtxs
920  ENDIF
921 
922 C...Calculations for Monte Carlo estimate of all cross-sections.
923  IF(mint(82).EQ.1.AND.isub.LE.90.OR.isub.GE.96) THEN
924  IF(mstp(142).LE.1) THEN
925  xsec(isub,2)=xsec(isub,2)+sigs
926  ELSE
927  xsec(isub,2)=xsec(isub,2)+sigswt
928  ENDIF
929  ELSEIF(mint(82).EQ.1) THEN
930  xsec(isub,2)=xsec(isub,2)+sigs
931  ENDIF
932  IF((isub.EQ.95.OR.isub.EQ.96).AND.loop2.EQ.1.AND.
933  &mint(82).EQ.1) xsec(97,2)=xsec(97,2)+siglpt
934 
935 C...Multiple interactions: store results of cross-section calculation.
936  IF(mint(50).EQ.1.AND.mstp(82).GE.3) THEN
937  vint(153)=sigsor
938  IF(mint(35).LE.1) CALL pymult(4)
939  IF(mint(35).GE.2) CALL pymign(4)
940  ENDIF
941 
942 C...Ratio of actual to maximum cross section.
943  IF(istsb.NE.11) THEN
944  viol=sigswt/xsec(isub,1)
945  IF(isub.EQ.96.AND.mstp(173).EQ.1) viol=viol/parp(174)
946  ELSEIF(idwtup.EQ.1.OR.idwtup.EQ.2) THEN
947  viol=xwgtup/xmaxup(kfpr(isub,1))
948  ELSEIF(idwtup.EQ.-1.OR.idwtup.EQ.-2) THEN
949  viol=abs(xwgtup)/abs(xmaxup(kfpr(isub,1)))
950  ELSE
951  viol=1d0
952  ENDIF
953 
954 C...Check that weight not negative.
955  IF(mstp(123).LE.0) THEN
956  IF(viol.LT.-1d-3) THEN
957  WRITE(mstu(11),5000) viol,ngen(0,3)+1
958  IF(mstp(122).GE.1) WRITE(mstu(11),5100) isub,vint(21),
959  & vint(22),vint(23),vint(26)
960  CALL pystop(2)
961  ENDIF
962  ELSE
963  IF(viol.LT.min(-1d-3,vint(109))) THEN
964  vint(109)=viol
965  IF(mstp(123).LE.2) WRITE(mstu(11),5200) viol,ngen(0,3)+1
966  IF(mstp(122).GE.1) WRITE(mstu(11),5100) isub,vint(21),
967  & vint(22),vint(23),vint(26)
968  ENDIF
969  ENDIF
970 
971 C...Weighting using estimate of maximum of differential cross-section.
972  ratnd=1d0
973  IF(mfail.EQ.0.AND.isub.NE.95.AND.isub.NE.96) THEN
974  IF(viol.LT.pyr(0)) THEN
975  IF(mint(121).GT.1) CALL pysave(2,iga)
976  IF(isub.GE.91.AND.isub.LE.94) isub=0
977  goto 100
978  ENDIF
979  ELSEIF(mfail.EQ.0) THEN
980  ratnd=siglpt/xsec(95,1)
981  viol=viol/ratnd
982  IF(loop2.EQ.1.AND.ratnd.LT.pyr(0)) THEN
983  IF(viol.GT.pyr(0).AND.mint(82).EQ.1.AND.msub(95).EQ.1.AND.
984  & (isub.LE.90.OR.isub.GE.95)) ngen(95,1)=ngen(95,1)+mint(143)
985  IF(mint(121).GT.1) CALL pysave(2,iga)
986  isub=0
987  goto 100
988  ENDIF
989  IF(viol.LT.pyr(0)) THEN
990  goto 140
991  ENDIF
992  ELSEIF(isub.NE.95.AND.isub.NE.96) THEN
993  IF(viol.LT.pyr(0)) THEN
994  msti(61)=1
995  IF(mint(121).GT.1) CALL pysave(2,iga)
996  RETURN
997  ENDIF
998  ELSE
999  ratnd=siglpt/xsec(95,1)
1000  IF(loop.EQ.1.AND.ratnd.LT.pyr(0)) THEN
1001  msti(61)=1
1002  IF(mint(121).GT.1) CALL pysave(2,iga)
1003  RETURN
1004  ENDIF
1005  viol=viol/ratnd
1006  IF(viol.LT.pyr(0)) THEN
1007  IF(mint(121).GT.1) CALL pysave(2,iga)
1008  goto 100
1009  ENDIF
1010  ENDIF
1011 
1012 C...Check for possible violation of estimated maximum of differential
1013 C...cross-section used in weighting.
1014  IF(mstp(123).LE.0) THEN
1015  IF(viol.GT.1d0) THEN
1016  WRITE(mstu(11),5300) viol,ngen(0,3)+1
1017  IF(mstp(122).GE.2) WRITE(mstu(11),5100) isub,vint(21),
1018  & vint(22),vint(23),vint(26)
1019  CALL pystop(2)
1020  ENDIF
1021  ELSEIF(mstp(123).EQ.1) THEN
1022  IF(viol.GT.vint(108)) THEN
1023  vint(108)=viol
1024  IF(viol.GT.1.0001d0) THEN
1025  mint(10)=1
1026  WRITE(mstu(11),5400) viol,ngen(0,3)+1
1027  IF(mstp(122).GE.2) WRITE(mstu(11),5100) isub,vint(21),
1028  & vint(22),vint(23),vint(26)
1029  ENDIF
1030  ENDIF
1031  ELSEIF(viol.GT.vint(108)) THEN
1032  vint(108)=viol
1033  IF(viol.GT.1d0) THEN
1034  mint(10)=1
1035  IF(mstp(123).EQ.2) WRITE(mstu(11),5400) viol,ngen(0,3)+1
1036  IF(istsb.EQ.11.AND.(iabs(idwtup).EQ.1.OR.iabs(idwtup).EQ.2))
1037  & THEN
1038  xmaxup(kfpr(isub,1))=viol*xmaxup(kfpr(isub,1))
1039  IF(kfpr(isub,1).LE.9) THEN
1040  IF(mstp(123).EQ.2) WRITE(mstu(11),5800) kfpr(isub,1),
1041  & xmaxup(kfpr(isub,1))
1042  ELSEIF(kfpr(isub,1).LE.99) THEN
1043  IF(mstp(123).EQ.2) WRITE(mstu(11),5900) kfpr(isub,1),
1044  & xmaxup(kfpr(isub,1))
1045  ELSE
1046  IF(mstp(123).EQ.2) WRITE(mstu(11),6000) kfpr(isub,1),
1047  & xmaxup(kfpr(isub,1))
1048  ENDIF
1049  ENDIF
1050  IF(istsb.NE.11.OR.iabs(idwtup).EQ.1) THEN
1051  xdif=xsec(isub,1)*(viol-1d0)
1052  xsec(isub,1)=xsec(isub,1)+xdif
1053  IF(msub(isub).EQ.1.AND.(isub.LE.90.OR.isub.GT.96))
1054  & xsec(0,1)=xsec(0,1)+xdif
1055  IF(mstp(122).GE.2) WRITE(mstu(11),5100) isub,vint(21),
1056  & vint(22),vint(23),vint(26)
1057  IF(isub.LE.9) THEN
1058  IF(mstp(123).EQ.2) WRITE(mstu(11),5500) isub,xsec(isub,1)
1059  ELSEIF(isub.LE.99) THEN
1060  IF(mstp(123).EQ.2) WRITE(mstu(11),5600) isub,xsec(isub,1)
1061  ELSE
1062  IF(mstp(123).EQ.2) WRITE(mstu(11),5700) isub,xsec(isub,1)
1063  ENDIF
1064  ENDIF
1065  vint(108)=1d0
1066  ENDIF
1067  ENDIF
1068 
1069 C...Multiple interactions: choose impact parameter (if not already done).
1070  IF(mint(39).EQ.0) vint(148)=1d0
1071  IF(mint(50).EQ.1.AND.(isub.LE.90.OR.isub.GE.96).AND.
1072  &mstp(82).GE.3) THEN
1073  IF(mint(35).LE.1) CALL pymult(5)
1074  IF(mint(35).GE.2) CALL pymign(5)
1075  IF(vint(150).LT.pyr(0)) THEN
1076  IF(mint(121).GT.1) CALL pysave(2,iga)
1077  IF(mfail.EQ.1) THEN
1078  msti(61)=1
1079  RETURN
1080  ENDIF
1081  goto 100
1082  ENDIF
1083  ENDIF
1084  IF(mint(82).EQ.1) ngen(0,2)=ngen(0,2)+1
1085  IF(mint(82).EQ.1.AND.msub(95).EQ.1) THEN
1086  IF(isub.LE.90.OR.isub.GE.95) ngen(95,1)=ngen(95,1)+mint(143)
1087  IF(isub.LE.90.OR.isub.GE.96) ngen(96,2)=ngen(96,2)+1
1088  ENDIF
1089  IF(isub.LE.90.OR.isub.GE.96) mint(31)=mint(31)+1
1090 
1091 C...Choose flavour of reacting partons (and subprocess).
1092  IF(istsb.GE.11) goto 320
1093  rsigs=sigs*pyr(0)
1094  qt2=vint(48)
1095  rqqbar=parp(87)*(1d0-(qt2/(qt2+(parp(88)*parp(82)*
1096  &(vint(1)/parp(89))**parp(90))**2))**2)
1097  IF(isub.NE.95.AND.(isub.NE.96.OR.mstp(82).LE.1.OR.
1098  &pyr(0).GT.rqqbar)) THEN
1099  DO 310 ichn=1,nchn
1100  kfl1=isig(ichn,1)
1101  kfl2=isig(ichn,2)
1102  mint(2)=isig(ichn,3)
1103  rsigs=rsigs-sigh(ichn)
1104  IF(rsigs.LE.0d0) goto 320
1105  310 CONTINUE
1106 
1107 C...Multiple interactions: choose qqbar preferentially at small pT.
1108  ELSEIF(isub.EQ.96) THEN
1109  mint(105)=mint(103)
1110  mint(109)=mint(107)
1111  CALL pyspli(mint(11),21,kfl1,kfldum)
1112  mint(105)=mint(104)
1113  mint(109)=mint(108)
1114  CALL pyspli(mint(12),21,kfl2,kfldum)
1115  mint(1)=11
1116  mint(2)=1
1117  IF(kfl1.EQ.kfl2.AND.pyr(0).LT.0.5d0) mint(2)=2
1118 
1119 C...Low-pT: choose string drawing configuration.
1120  ELSE
1121  kfl1=21
1122  kfl2=21
1123  rsigs=6d0*pyr(0)
1124  mint(2)=1
1125  IF(rsigs.GT.1d0) mint(2)=2
1126  IF(rsigs.GT.2d0) mint(2)=3
1127  ENDIF
1128 
1129 C...Reassign QCD process. Partons before initial state radiation.
1130  320 IF(mint(2).GT.10) THEN
1131  mint(1)=mint(2)/10
1132  mint(2)=mod(mint(2),10)
1133  ENDIF
1134  IF(mint(82).EQ.1.AND.mstp(111).GE.0) ngen(mint(1),2)=
1135  &ngen(mint(1),2)+1
1136  mint(15)=kfl1
1137  mint(16)=kfl2
1138  mint(13)=mint(15)
1139  mint(14)=mint(16)
1140  vint(141)=vint(41)
1141  vint(142)=vint(42)
1142  vint(151)=0d0
1143  vint(152)=0d0
1144 
1145 C...Calculate x value of photon for parton inside photon inside e.
1146  DO 350 jt=1,2
1147  mint(18+jt)=0
1148  vint(154+jt)=0d0
1149  mspli=0
1150  IF(jt.EQ.1.AND.mint(43).LE.2) mspli=1
1151  IF(jt.EQ.2.AND.mod(mint(43),2).EQ.1) mspli=1
1152  IF(iabs(mint(14+jt)).LE.8.OR.mint(14+jt).EQ.21) mspli=mspli+1
1153  IF(mspli.EQ.2) THEN
1154  kflh=mint(14+jt)
1155  xhrd=vint(140+jt)
1156  q2hrd=vint(54)
1157  mint(105)=mint(102+jt)
1158  mint(109)=mint(106+jt)
1159  vint(120)=vint(2+jt)
1160  IF(mstp(57).LE.1) THEN
1161  CALL pypdfu(22,xhrd,q2hrd,xpq)
1162  ELSE
1163  CALL pypdfl(22,xhrd,q2hrd,xpq)
1164  ENDIF
1165  wtmx=4d0*xpq(kflh)
1166  IF(mstp(13).EQ.2) THEN
1167  q2pms=q2hrd/pmas(11,1)**2
1168  wtmx=wtmx*log(max(2d0,q2pms*(1d0-xhrd)/xhrd**2))
1169  ENDIF
1170  330 xe=xhrd**pyr(0)
1171  xg=min(1d0-1d-10,xhrd/xe)
1172  IF(mstp(57).LE.1) THEN
1173  CALL pypdfu(22,xg,q2hrd,xpq)
1174  ELSE
1175  CALL pypdfl(22,xg,q2hrd,xpq)
1176  ENDIF
1177  wt=(1d0+(1d0-xe)**2)*xpq(kflh)
1178  IF(mstp(13).EQ.2) wt=wt*log(max(2d0,q2pms*(1d0-xe)/xe**2))
1179  IF(wt.LT.pyr(0)*wtmx) goto 330
1180  mint(18+jt)=1
1181  vint(154+jt)=xe
1182  DO 340 kfls=-25,25
1183  xsfx(jt,kfls)=xpq(kfls)
1184  340 CONTINUE
1185  ENDIF
1186  350 CONTINUE
1187 
1188 C...Pick scale where photon is resolved.
1189  q0s=parp(15)**2
1190  q1s=vint(154)**2
1191  vint(283)=0d0
1192  IF(mint(107).EQ.3) THEN
1193  IF(mstp(66).EQ.1) THEN
1194  vint(283)=q0s*(vint(54)/q0s)**pyr(0)
1195  ELSEIF(mstp(66).EQ.2) THEN
1196  ps=vint(3)**2
1197  q2eff=vint(54)*((q0s+ps)/(vint(54)+ps))*
1198  & exp(ps*(vint(54)-q0s)/((vint(54)+ps)*(q0s+ps)))
1199  q2int=sqrt(q0s*q2eff)
1200  vint(283)=q2int*(vint(54)/q2int)**pyr(0)
1201  ELSEIF(mstp(66).EQ.3) THEN
1202  vint(283)=q0s*(q1s/q0s)**pyr(0)
1203  ELSEIF(mstp(66).GE.4) THEN
1204  ps=0.25d0*vint(3)**2
1205  vint(283)=(q0s+ps)*(q1s+ps)/
1206  & (q0s+pyr(0)*(q1s-q0s)+ps)-ps
1207  ENDIF
1208  ENDIF
1209  vint(284)=0d0
1210  IF(mint(108).EQ.3) THEN
1211  IF(mstp(66).EQ.1) THEN
1212  vint(284)=q0s*(vint(54)/q0s)**pyr(0)
1213  ELSEIF(mstp(66).EQ.2) THEN
1214  ps=vint(4)**2
1215  q2eff=vint(54)*((q0s+ps)/(vint(54)+ps))*
1216  & exp(ps*(vint(54)-q0s)/((vint(54)+ps)*(q0s+ps)))
1217  q2int=sqrt(q0s*q2eff)
1218  vint(284)=q2int*(vint(54)/q2int)**pyr(0)
1219  ELSEIF(mstp(66).EQ.3) THEN
1220  vint(284)=q0s*(q1s/q0s)**pyr(0)
1221  ELSEIF(mstp(66).GE.4) THEN
1222  ps=0.25d0*vint(4)**2
1223  vint(284)=(q0s+ps)*(q1s+ps)/
1224  & (q0s+pyr(0)*(q1s-q0s)+ps)-ps
1225  ENDIF
1226  ENDIF
1227  IF(mint(121).GT.1) CALL pysave(2,iga)
1228 
1229 C...Format statements for differential cross-section maximum violations.
1230  5000 FORMAT(/1x,'Error: negative cross-section fraction',1p,d11.3,1x,
1231  &'in event',1x,i7,'D0'/1x,'Execution stopped!')
1232  5100 FORMAT(1x,'ISUB = ',i3,'; Point of violation:'/1x,'tau =',1p,
1233  &d11.3,', y* =',d11.3,', cthe = ',0p,f11.7,', tau'' =',1p,d11.3)
1234  5200 FORMAT(/1x,'Warning: negative cross-section fraction',1p,d11.3,1x,
1235  &'in event',1x,i7)
1236  5300 FORMAT(/1x,'Error: maximum violated by',1p,d11.3,1x,
1237  &'in event',1x,i7,'D0'/1x,'Execution stopped!')
1238  5400 FORMAT(/1x,'Advisory warning: maximum violated by',1p,d11.3,1x,
1239  &'in event',1x,i7)
1240  5500 FORMAT(1x,'XSEC(',i1,',1) increased to',1p,d11.3)
1241  5600 FORMAT(1x,'XSEC(',i2,',1) increased to',1p,d11.3)
1242  5700 FORMAT(1x,'XSEC(',i3,',1) increased to',1p,d11.3)
1243  5800 FORMAT(1x,'XMAXUP(',i1,') increased to',1p,d11.3)
1244  5900 FORMAT(1x,'XMAXUP(',i2,') increased to',1p,d11.3)
1245  6000 FORMAT(1x,'XMAXUP(',i3,') increased to',1p,d11.3)
1246 
1247  RETURN
1248  END