11 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
16 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
17 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
19 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
22 common/
pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
30 IF(isub.EQ.96) goto 100
34 IF(abs(sqm3).LT.1d-4.AND.abs(sqm4).LT.1d-4)
THEN
35 ckin09=
max(ckin(9),ckin(13))
36 ckin10=min(ckin(10),ckin(14))
37 ckin11=
max(ckin(11),ckin(15))
38 ckin12=min(ckin(12),ckin(16))
40 ckin09=
max(ckin(9),min(0d0,ckin(13)))
41 ckin10=min(ckin(10),
max(0d0,ckin(14)))
42 ckin11=
max(ckin(11),min(0d0,ckin(15)))
43 ckin12=min(ckin(12),
max(0d0,ckin(16)))
50 be34=sqrt(
max(1d-20,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
53 IF(min(sqm3,sqm4).LT.ckin(6)**2.AND.istsb.NE.1.AND.istsb.NE.3)
54 &pthmin=
max(ckin(3),ckin(5))
63 IF(istsb.GE.3.AND.istsb.LE.5) taue=taup
64 x1=sqrt(taue)*exp(yst)
65 x2=sqrt(taue)*exp(-yst)
67 IF(
mint(47).NE.1)
THEN
69 IF(ckin(2).GE.0d0.AND.
tau*
vint(2).GT.ckin(2)**2)
mint(51)=1
70 IF(yst.LT.ckin(7).OR.yst.GT.ckin(8))
mint(51)=1
71 IF(xf.LT.ckin(25).OR.xf.GT.ckin(26))
mint(51)=1
73 IF(
mint(45).NE.1)
THEN
74 IF(x1.LT.ckin(21).OR.x1.GT.ckin(22))
mint(51)=1
76 IF(
mint(46).NE.1)
THEN
77 IF(x2.LT.ckin(23).OR.x2.GT.ckin(24))
mint(51)=1
79 IF(
mint(45).EQ.2)
THEN
82 IF(
mint(46).EQ.2)
THEN
85 IF(istsb.EQ.2.OR.istsb.EQ.4)
THEN
86 pth=0.5d0*be34*sqrt(
tau*
vint(2)*
max(0d0,1d0-cth**2))
87 expy3=
max(1d-20,(1d0+rm3-rm4+be34*cth)/
88 &
max(1d-20,(1d0+rm3-rm4-be34*cth)))
89 expy4=
max(1d-20,(1d0-rm3+rm4-be34*cth)/
90 &
max(1d-20,(1d0-rm3+rm4+be34*cth)))
91 y3=yst+0.5d0*
log(expy3)
92 y4=yst+0.5d0*
log(expy4)
97 sth=sqrt(
max(0d0,1d0-cth**2))
98 exsq3=sqrt(
max(1d-20,((1d0+rm3-rm4)*cosh(yst)+be34*sinh(yst)*
100 exsq4=sqrt(
max(1d-20,((1d0-rm3+rm4)*cosh(yst)-be34*sinh(yst)*
102 IF(sth.GE.1d-10)
THEN
103 expet3=((1d0+rm3-rm4)*sinh(yst)+be34*cosh(yst)*cth+exsq3)/
105 expet4=((1d0-rm3+rm4)*sinh(yst)-be34*cosh(yst)*cth+exsq4)/
107 eta3=
log(min(1d10,
max(1d-10,expet3)))
108 eta4=
log(min(1d10,
max(1d-10,expet4)))
109 etalar=
max(eta3,eta4)
110 etasma=min(eta3,eta4)
112 cts3=((1d0+rm3-rm4)*sinh(yst)+be34*cosh(yst)*cth)/exsq3
113 cts4=((1d0-rm3+rm4)*sinh(yst)-be34*cosh(yst)*cth)/exsq4
114 ctslar=min(1d0,
max(-1d0,cts3,cts4))
115 ctssma=
max(-1d0,min(1d0,cts3,cts4))
117 rpts=4d0*
vint(71)**2/sh
118 be34l=sqrt(
max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4-rpts))
119 rm34=
max(1d-20,2d0*rm3*rm4)
122 rthm=(4d0*rm3*rm4+rpts)/(1d0-rm3-rm4+be34l)
123 tha=0.5d0*sh*
max(rthm,1d0-rm3-rm4-be34*cth)
124 uha=0.5d0*sh*
max(rthm,1d0-rm3-rm4+be34*cth)
125 IF(pth.LT.pthmin)
mint(51)=1
126 IF(ckin(4).GE.0d0.AND.pth.GT.ckin(4))
mint(51)=1
127 IF(ylarge.LT.ckin(9).OR.ylarge.GT.ckin(10))
mint(51)=1
128 IF(ysmall.LT.ckin(11).OR.ysmall.GT.ckin(12))
mint(51)=1
129 IF(etalar.LT.ckin(13).OR.etalar.GT.ckin(14))
mint(51)=1
130 IF(etasma.LT.ckin(15).OR.etasma.GT.ckin(16))
mint(51)=1
131 IF(ctslar.LT.ckin(17).OR.ctslar.GT.ckin(18))
mint(51)=1
132 IF(ctssma.LT.ckin(19).OR.ctssma.GT.ckin(20))
mint(51)=1
133 IF(cth.LT.ckin(27).OR.cth.GT.ckin(28))
mint(51)=1
134 IF(tha.LT.ckin(35))
mint(51)=1
135 IF(ckin(36).GE.0d0.AND.tha.GT.ckin(36))
mint(51)=1
136 IF(uha.LT.ckin(37))
mint(51)=1
137 IF(ckin(38).GE.0d0.AND.uha.GT.ckin(38))
mint(51)=1
139 IF(istsb.GE.3.AND.istsb.LE.5)
THEN
140 IF(taup*
vint(2).LT.ckin(31)**2)
mint(51)=1
141 IF(ckin(32).GE.0d0.AND.taup*
vint(2).GT.ckin(32)**2)
mint(51)=1
145 IF(isub.EQ.10.AND.
mint(43).GE.2)
THEN
147 IF(iabs(
mint(12)).LT.20) xbj=x1
149 w2bj=q2bj*(1d0-xbj)/xbj
150 IF(w2bj.LT.ckin(39))
mint(51)=1
151 IF(ckin(40).GT.0d0.AND.w2bj.GT.ckin(40))
mint(51)=1
154 ELSEIF(ilim.EQ.1)
THEN
160 taumn1=ckin(1)**2/
vint(2)
162 IF(ckin(2).GE.0d0) taumx1=ckin(2)**2/
vint(2)
164 tm3=sqrt(sqm3+pthmin**2)
165 tm4=sqrt(sqm4+pthmin**2)
167 IF(ckin09.GT.ckin12) ydcosh=cosh(ckin09-ckin12)
168 taumn2=(tm3**2+2d0*tm3*tm4*ydcosh+tm4**2)/
vint(2)
171 cth2mn=min(ckin(27)**2,ckin(28)**2)
172 cth2mx=
max(ckin(27)**2,ckin(28)**2)
174 IF(ckin(27)*ckin(28).GT.0d0) taumn3=
175 & (sqrt(sqm3+pthmin**2/(1d0-cth2mn))+
176 & sqrt(sqm4+pthmin**2/(1d0-cth2mn)))**2/
vint(2)
178 IF(ckin(4).GE.0d0.AND.cth2mx.LT.1d0) taumx3=
179 & (sqrt(sqm3+ckin(4)**2/(1d0-cth2mx))+
180 & sqrt(sqm4+ckin(4)**2/(1d0-cth2mx)))**2/
vint(2)
182 taumn4=ckin(21)*ckin(23)
183 taumx4=ckin(22)*ckin(24)
186 taumx5=
max(1d0-ckin(25),1d0+ckin(26))
188 taumn6=(sqm3+sqm4+ckin(35)+ckin(37))/
vint(2)
190 IF(ckin(36).GT.0d0.AND.ckin(38).GT.0d0) taumx6=
191 & (sqm3+sqm4+ckin(36)+ckin(38))/
vint(2)
194 vint(11)=
max(taumn0,taumn1,taumn2,taumn3,taumn4,taumn5,taumn6)
195 vint(31)=min(taumx0,taumx1,taumx2,taumx3,taumx4,taumx5,taumx6)
196 IF(
mint(47).EQ.1.AND.(istsb.EQ.1.OR.istsb.EQ.2))
THEN
199 ELSEIF(
mint(47).EQ.5)
THEN
201 ELSEIF(
mint(47).GE.6)
THEN
206 ELSEIF(ilim.EQ.2)
THEN
209 IF(istsb.GE.3.AND.istsb.LE.5) taue=
vint(26)
218 ystmn2=
log(
max(taue,ckin(21))/taurt)
219 ystmx2=
log(
max(taue,ckin(22))/taurt)
221 ystmn3=-
log(
max(taue,ckin(24))/taurt)
222 ystmx3=-
log(
max(taue,ckin(23))/taurt)
224 yepmn4=0.5d0*abs(ckin(25))/taurt
225 ystmn4=
sign(
log(
max(1d-20,sqrt(1d0+yepmn4**2)+yepmn4)),ckin(25))
226 yepmx4=0.5d0*abs(ckin(26))/taurt
227 ystmx4=
sign(
log(
max(1d-20,sqrt(1d0+yepmx4**2)+yepmx4)),ckin(26))
229 yepsmn=(rm3-rm4)*sinh(ckin09-ckin11)
230 yepsmx=(rm3-rm4)*sinh(ckin10-ckin12)
231 ydifmn=abs(
log(
max(1d-20,sqrt(1d0+yepsmn**2)-yepsmn)))
232 ydifmx=abs(
log(
max(1d-20,sqrt(1d0+yepsmx**2)-yepsmx)))
233 ystmn5=0.5d0*(ckin09+ckin11-ydifmn)
234 ystmx5=0.5d0*(ckin10+ckin12+ydifmx)
237 cthlim=sqrt(
max(0d0,1d0-4d0*pthmin**2/(be34**2*taue*
vint(2))))
238 rzmn=be34*
max(ckin(27),-cthlim)
239 rzmx=be34*min(ckin(28),cthlim)
240 yex3mx=(1d0+rm3-rm4+rzmx)/
max(1d-10,1d0+rm3-rm4-rzmx)
241 yex4mx=(1d0+rm4-rm3-rzmn)/
max(1d-10,1d0+rm4-rm3+rzmn)
242 yex3mn=
max(1d-10,1d0+rm3-rm4+rzmn)/(1d0+rm3-rm4-rzmn)
243 yex4mn=
max(1d-10,1d0+rm4-rm3-rzmx)/(1d0+rm4-rm3+rzmx)
244 ystmn6=ckin09-0.5d0*
log(
max(yex3mx,yex4mx))
245 ystmx6=ckin12-0.5d0*
log(min(yex3mn,yex4mn))
248 vint(12)=
max(ystmn0,ystmn1,ystmn2,ystmn3,ystmn4,ystmn5,ystmn6)
249 vint(32)=min(ystmx0,ystmx1,ystmx2,ystmx3,ystmx4,ystmx5,ystmx6)
250 IF(
mint(47).EQ.1)
THEN
253 ELSEIF(
mint(47).EQ.2.OR.
mint(47).EQ.6)
THEN
254 vint(12)=(1d0-1d-9)*ystmx0
255 vint(32)=(1d0+1d-9)*ystmx0
256 ELSEIF(
mint(47).EQ.3.OR.
mint(47).EQ.7)
THEN
257 vint(12)=-(1d0+1d-9)*ystmx0
258 vint(32)=-(1d0-1d-9)*ystmx0
259 ELSEIF(
mint(47).EQ.5)
THEN
260 ystee=
log((1d0-1d-10)/taurt)
266 ELSEIF(ilim.EQ.3)
THEN
275 ctnmn1=min(0d0,ckin(27))
276 ctnmx1=min(0d0,ckin(28))
277 ctpmn1=
max(0d0,ckin(27))
278 ctpmx1=
max(0d0,ckin(28))
280 ctnmn2=-sqrt(
max(0d0,1d0-4d0*pthmin**2/(be34**2*
tau*
vint(2))))
284 IF(ckin(4).GE.0d0)
THEN
285 ctnmx2=-sqrt(
max(0d0,1d0-4d0*ckin(4)**2/
290 ctnmn3=min(0d0,
max((1d0+rm3-rm4)/be34*tanh(ckin11-yst),
291 & -(1d0-rm3+rm4)/be34*tanh(ckin10-yst)))
292 ctnmx3=min(0d0,(1d0+rm3-rm4)/be34*tanh(ckin12-yst),
293 & -(1d0-rm3+rm4)/be34*tanh(ckin09-yst))
294 ctpmn3=
max(0d0,(1d0+rm3-rm4)/be34*tanh(ckin09-yst),
295 & -(1d0-rm3+rm4)/be34*tanh(ckin12-yst))
296 ctpmx3=
max(0d0,min((1d0+rm3-rm4)/be34*tanh(ckin10-yst),
297 & -(1d0-rm3+rm4)/be34*tanh(ckin11-yst)))
304 IF(ckin(35).GT.0d0)
THEN
305 ctlim=(1d0-rm3-rm4-2d0*ckin(35)/sh)/be34
306 IF(ctlim.GT.0d0)
THEN
313 IF(ckin(36).GT.0d0)
THEN
314 ctlim=(1d0-rm3-rm4-2d0*ckin(36)/sh)/be34
315 IF(ctlim.LT.0d0)
THEN
327 IF(ckin(37).GT.0d0)
THEN
328 ctlim=(2d0*ckin(37)/sh-(1d0-rm3-rm4))/be34
329 IF(ctlim.LT.0d0)
THEN
336 IF(ckin(38).GT.0d0)
THEN
337 ctlim=(2d0*ckin(38)/sh-(1d0-rm3-rm4))/be34
338 IF(ctlim.GT.0d0)
THEN
347 vint(13)=
max(ctnmn0,ctnmn1,ctnmn2,ctnmn3,ctnmn4,ctnmn5)
348 vint(33)=min(ctnmx0,ctnmx1,ctnmx2,ctnmx3,ctnmx4,ctnmx5)
349 vint(14)=
max(ctpmn0,ctpmn1,ctpmn2,ctpmn3,ctpmn4,ctpmn5)
350 vint(34)=min(ctpmx0,ctpmx1,ctpmx2,ctpmx3,ctpmx4,ctpmx5)
356 ELSEIF(ilim.EQ.4)
THEN
360 IF(istsb.EQ.5.AND.
vint(201).GT.0d0)
THEN
362 tapmn0=(sqrt(
tau)+pqrat)**2
366 tapmn1=ckin(31)**2/
vint(2)
368 IF(ckin(32).GE.0d0) tapmx1=ckin(32)**2/
vint(2)
372 vint(36)=min(tapmx0,tapmx1)
373 IF(
mint(47).EQ.1)
THEN
376 ELSEIF(
mint(47).EQ.5)
THEN
378 ELSEIF(
mint(47).EQ.6.OR.
mint(47).EQ.7)
THEN
388 100
IF(ilim.EQ.0)
THEN
389 ELSEIF(ilim.EQ.1)
THEN
390 IF(
mstp(82).LE.1)
THEN
397 ELSEIF(ilim.EQ.2)
THEN
400 ELSEIF(ilim.EQ.3)
THEN
401 IF(
mstp(82).LE.1)
THEN
408 vint(13)=-sqrt(
max(0d0,1d0-st2eff))