File indexing completed on 2023-10-25 09:48:24
0001
0002
0003
0004 double precision function totfun(zup,papawt)
0005 implicit double precision (a-h,o-z)
0006 implicit integer (i-n)
0007
0008 #include "inclcon.h"
0009
0010
0011 common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
0012 common/pypars/mstp(200),parp(200),msti(200),pari(200)
0013 parameter (maxnup=500)
0014 common/hepeup/nup,idprup,xwgtup,scalup,aqedup,aqcdup,idup(maxnup),
0015 &istup(maxnup),mothup(2,maxnup),icolup(2,maxnup),pup(5,maxnup),
0016 &vtimup(maxnup),spinup(maxnup)
0017 save /hepeup/
0018
0019
0020 double complex colmat,bundamp
0021 common/upcom/ecm,pmbc,pmb,pmc,fbcc,pmomup(5,8),
0022 & colmat(10,64),bundamp(4),pmomzero(5,8)
0023
0024 common/counter/ibcstate,nev
0025 common/rconst/pi
0026 common/ptpass/ptmin,ptmax,crossmax,etamin,etamax,
0027 & smin,smax,ymin,ymax,psetamin,psetamax
0028 common/confine/ptcut,etacut
0029 common/colflow/amp2cf(10),smatval
0030
0031
0032 common/funtrans/nq2,npdfu
0033
0034
0035 common/subopen/subfactor,subenergy,isubonly
0036
0037
0038 logical generate
0039 common/genefull/generate
0040
0041
0042 common/extraz/zfactor,zmin,zmax
0043 common/outpdf/ioutpdf,ipdfnum
0044 common/intinip/iinip
0045 common/intinif/iinif
0046
0047
0048
0049 common/qqbar/iqqbar,iqcode
0050
0051 dimension xpp(-25:25),xppbar(-25:25),zup(7),pboo(4),pc(4),pl(4)
0052 data conv/3.8938573d+8/ !pb
0053 common /ppp/ pp(4,40),guv(4)
0054
0055
0056
0057 totfun=0.0d0
0058 phase = 0.
0059
0060
0061 if(isubonly.eq.1) then
0062 x1=subfactor
0063 x2=subfactor
0064 else
0065 taumin =((pmbc+pmb+pmc)/ecm)**2
0066 taumax =1.0d0
0067 tau=(taumax-taumin)*zup(6)+taumin
0068 yymin= dlog(dsqrt(tau))
0069 yymax=-dlog(dsqrt(tau))
0070 yy =(yymax-yymin)*zup(7)+yymin
0071 x1 =dsqrt(tau)*exp(yy)
0072 x2 =dsqrt(tau)*exp(-yy)
0073 end if
0074
0075
0076
0077
0078 pup(1,1)= 0.0d0
0079 pup(2,1)= 0.0d0
0080 pup(3,1)= ecm*x1/2.0d0
0081 pup(4,1)= ecm*x1/2.0d0
0082 pup(5,1)= 0.0d0
0083
0084
0085 pup(1,2)= 0.0d0
0086 pup(2,2)= 0.0d0
0087 pup(3,2)=-ecm*x2/2.0d0
0088 pup(4,2)= ecm*x2/2.0d0
0089 pup(5,2)= 0.0d0
0090
0091
0092
0093 do i=1,4
0094 pboo(i)=pup(i,1)+pup(i,2)
0095 end do
0096
0097 do 101, i=1,3
0098 do j=1,4
0099 pc(j)=pmomup(j,i+2)
0100 end do
0101 call lorentz(pboo,pc,pl)
0102 do j=1,4
0103 pmomup(j,i+2)=pl(j)
0104 end do
0105 101 continue
0106
0107
0108 do i=3,5
0109 do j=1,5
0110 pup(j,i)=pmomup(j,i)
0111 end do
0112 end do
0113
0114
0115 do i=1,5
0116 pmomup(i,6)=pmb/(pmb+pmc)*pmomup(i,3)
0117 pmomup(i,7)=pmc/(pmb+pmc)*pmomup(i,3)
0118 end do
0119
0120
0121 do i=1,2
0122 do j=1,5
0123 pmomup(j,i)=pup(j,i)
0124 end do
0125 end do
0126
0127
0128 ptbc =dsqrt(pup(1,3)**2+pup(2,3)**2)
0129 pr =max(1.0d-16,pup(5,3)**2+pup(1,3)**2+pup(2,3)**2)
0130 prs =max(1.0d-16,pup(1,3)**2+pup(2,3)**2)
0131 eta =sign(dlog(min((dsqrt(pr+pup(3,3)**2)+dabs(pup(3,3)))
0132 & /dsqrt(pr),1.0d+20)),pup(3,3))
0133 pseta =sign(dlog(min((dsqrt(prs+pup(3,3)**2)+dabs(pup(3,3)))
0134 & /dsqrt(prs),1.0d+20)),pup(3,3))
0135
0136
0137 if(ptbc.lt.ptcut .or. abs(eta).gt.etacut) then
0138 if (generate) then
0139 do ii=1,10
0140 amp2cf(ii)=0.0d0
0141 end do
0142 end if
0143 smatval=0.0d0
0144 return
0145 end if
0146
0147
0148 if(nq2.eq.1) q2 =x1*x2*ecm**2/4.0d0
0149 if(nq2.eq.2) q2 =x1*x2*ecm**2
0150 if(nq2.eq.3.or.nq2.eq.8) q2 =ptbc**2.0d0+pup(5,3)**2
0151 if(nq2.eq.4) then
0152 q=0.0d0
0153 do i=3,5
0154 q=q+dsqrt(pup(1,i)**2+pup(2,i)**2+pup(5,i)**2)
0155 end do
0156 q2=q**2
0157 end if
0158 if(nq2.eq.5) then
0159 q=0.0d0
0160 do i=3,5
0161 q=q+dsqrt(pup(1,i)**2.0d0+pup(2,i)**2.0d0+pup(5,i)**2.0d0)
0162 end do
0163 q2=(q/3.0d0)**2.0d0
0164 end if
0165 if(nq2.eq.6.or.nq2.eq.7) then
0166 q2=pmb**2+pup(1,4)**2+pup(2,4)**2
0167 end if
0168
0169 if(nq2.eq.9) then
0170 q2=4.0d0*pmb**2
0171 end if
0172 alps = 0.00
0173 alps2 = 0.00
0174
0175 if(ioutpdf.eq.1) then
0176 if(ipdfnum.eq.100) alps =alpgrv(q2,1)*4*pi
0177 if(ipdfnum.eq.200) alps =alpmsrt(dsqrt(q2),0.220d0,0)
0178 if(ipdfnum.eq.300) alps =alpcteq(q2,1)*4*pi
0179 else
0180 alps =pyalps(q2)
0181 end if
0182
0183
0184
0185 if(nq2.eq.6.or.nq2.eq.8) then
0186 alps1=alps
0187 if(ioutpdf.eq.1) then
0188 q2=4.0d0*pmc**2.0d0
0189 if(ipdfnum.eq.100) alps2 =alpgrv(q2,1)*4*pi
0190 if(ipdfnum.eq.200) alps2 =alpmsrt(dsqrt(q2),0.22d0,0)
0191 if(ipdfnum.eq.300) alps2 =alpcteq(q2,1)*4*pi
0192 else
0193 alps2 =pyalps(4.0d0*pmc**2.0d0)
0194 end if
0195 alps =dsqrt(alps1*alps2)
0196 end if
0197
0198
0199 scalup=dsqrt(q2)
0200 aqcdup=alps
0201
0202 if(isubonly.eq.0) then
0203 if(ioutpdf.eq.0) then
0204
0205
0206 if(npdfu.eq.1) then
0207 call pypdfu(2212,x1,q2,xpp)
0208 call pypdfu(-2212,x2,q2,xppbar)
0209 end if
0210
0211
0212 if(npdfu.eq.2) then
0213 call pypdfu(2212,x1,q2,xpp)
0214 call pypdfu(2212,x2,q2,xppbar)
0215 end if
0216 else
0217 if(ipdfnum.eq.100) then
0218 call grv98pa(1, x1, q2, uv, dv, us, ds, ss, gl1)
0219 call grv98pa(1, x2, q2, uv, dv, us, ds, ss, gl2)
0220 if(iqqbar.eq.0) then
0221 xpp(21)=gl1
0222 xppbar(21)=gl2
0223 else
0224 if(iqcode.eq.1) then
0225 if(npdfu.eq.1) then !tevatron
0226 xpp(iqcode)=uv+us
0227 xpp(-iqcode)=us
0228 xppbar(iqcode)=uv+us
0229 xppbar(-iqcode)=us
0230 end if
0231 if(npdfu.eq.2) then !lhc
0232 xpp(iqcode)=uv+us
0233 xpp(-iqcode)=us
0234 xppbar(iqcode)=us
0235 xppbar(-iqcode)=uv+us
0236 end if
0237 end if
0238 if(iqcode.eq.2) then
0239 if(npdfu.eq.1) then !tevatron (u-p,~u-~p and u-~p,~u-p)
0240 xpp(iqcode)=dv+ds
0241 xpp(-iqcode)=ds
0242 xppbar(iqcode)=dv+ds
0243 xppbar(-iqcode)=ds
0244 end if
0245 if(npdfu.eq.2) then !lhc (u-p,~u-p and u-p,~u-p)
0246 xpp(iqcode)=dv+ds
0247 xpp(-iqcode)=ds
0248 xppbar(iqcode)=ds
0249 xppbar(-iqcode)=dv+ds
0250 end if
0251 end if
0252 if(iqcode.eq.3) then ! the same for tevatron or lhc
0253 xpp(3)=2.0d0*ss
0254 xppbar(3)=2.0d0*ss
0255 end if
0256 end if
0257 end if
0258 if(ipdfnum.eq.200) then
0259 qq=dsqrt(q2)
0260 call mrstlo(x1,qq,1,upv,dnv,usea,dsea,str,chm,bot,glu1)
0261 call mrstlo(x2,qq,1,upv,dnv,usea,dsea,str,chm,bot,glu2)
0262 if(iqqbar.eq.0) then
0263 xpp(21)=glu1
0264 xppbar(21)=glu2
0265 else
0266 if(iqcode.eq.1) then
0267 if(npdfu.eq.1) then !tevatron (u-p,~u-~p and u-~p,~u-p)
0268 xpp(iqcode)=upv+usea
0269 xpp(-iqcode)=usea
0270 xppbar(iqcode)=upv+usea
0271 xppbar(-iqcode)=usea
0272 end if
0273 if(npdfu.eq.2) then !lhc (u-p,~u-p and u-p,~u-p)
0274 xpp(iqcode)=upv+usea
0275 xpp(-iqcode)=usea
0276 xppbar(iqcode)=usea
0277 xppbar(-iqcode)=upv+usea
0278 end if
0279 end if
0280 if(iqcode.eq.2) then
0281 if(npdfu.eq.1) then !tevatron (u-p,~u-~p and u-~p,~u-p)
0282 xpp(iqcode)=dnv+dsea
0283 xpp(-iqcode)=dsea
0284 xppbar(iqcode)=dnv+dsea
0285 xppbar(-iqcode)=dsea
0286 end if
0287 if(npdfu.eq.2) then !lhc (u-p,~u-p and u-p,~u-p)
0288 xpp(iqcode)=dnv+dsea
0289 xpp(-iqcode)=dsea
0290 xppbar(iqcode)=dsea
0291 xppbar(-iqcode)=dnv+dsea
0292 end if
0293 end if
0294 if(iqcode.eq.3) then
0295 xpp(3)=str*2.0d0
0296 xppbar(3)=str*2.0d0
0297 end if
0298 end if
0299 end if
0300 if(ipdfnum.eq.300) then
0301 qq=dsqrt(q2)
0302
0303 if(iqqbar.eq.0) then
0304 xpp(21) =ctq6pdf(0,x1,qq)
0305 xppbar(21)=ctq6pdf(0,x2,qq)
0306 else
0307 if(npdfu.eq.1) then !tevatron
0308 xpp(iqcode)=ctq6pdf(iqcode,x1,qq)
0309 xpp(-iqcode)=ctq6pdf(-iqcode,x1,qq)
0310 xppbar(iqcode)=ctq6pdf(iqcode,x2,qq)
0311 xppbar(-iqcode)=ctq6pdf(-iqcode,x2,qq)
0312 end if
0313 if(npdfu.eq.2) then !lhc
0314 xpp(iqcode)=ctq6pdf(iqcode,x1,qq)
0315 xpp(-iqcode)=ctq6pdf(-iqcode,x1,qq)
0316 xppbar(iqcode)=ctq6pdf(-iqcode,x2,qq)
0317 xppbar(-iqcode)=ctq6pdf(iqcode,x2,qq)
0318 end if
0319 end if
0320 end if
0321 end if
0322 end if
0323
0324
0325
0326 if(xpp(21).lt.1.0d-16) xpp(21)=0.0d0
0327 if(xppbar(21).lt.1.0d-16) xppbar(21)=0.0d0
0328
0329 if(iqqbar.eq.1) then
0330 if(xpp(iqcode).lt.1.0d-16) xpp(iqcode)=0.0d0
0331 if(xppbar(iqcode).lt.1.0d-16) xppbar(iqcode)=0.0d0
0332 if(xpp(-iqcode).lt.1.0d-16) xpp(-iqcode)=0.0d0
0333 if(xppbar(-iqcode).lt.1.0d-16) xppbar(-iqcode)=0.0d0
0334 end if
0335
0336
0337 if(isubonly.eq.1) alps=0.20d0
0338
0339
0340 if(ibcstate.eq.3) goto 1005
0341 if(ibcstate.eq.4) goto 1005
0342 if(ibcstate.eq.5) goto 1005
0343 if(ibcstate.eq.6) goto 1005
0344
0345
0346 phase =papawt*alps**4.0d0/(2.0d0**11*pi*3.0d0*dotup(1,2))
0347
0348
0349
0350
0351
0352
0353
0354
0355
0356
0357
0358
0359
0360 call xsection(sigscl,sigvct)
0361
0362
0363
0364
0365 if(iqqbar.eq.1) then
0366 phase=phase*(2.0d0**6/3.0d0**2)
0367 end if
0368
0369
0370
0371 sigscl=conv*phase*sigscl
0372 sigvct=conv*phase*sigvct
0373 sigcross = 0.0
0374 if(ibcstate.eq.1) sigcross=sigscl
0375 if(ibcstate.eq.2) sigcross=sigvct
0376
0377
0378 if(isubonly.eq.0) then
0379 if(iqqbar.eq.0) then
0380 totfun =sigcross*xpp(21)*xppbar(21)/x1/x2
0381 if(ioutpdf.eq.1 .and. ipdfnum.eq.300) then
0382 totfun=sigcross*xpp(21)*xppbar(21)
0383 end if
0384 else
0385 totfun =sigcross*(xpp(iqcode)*xppbar(iqcode)+
0386 & xpp(-iqcode)*xppbar(-iqcode))/x1/x2
0387 if(ioutpdf.eq.1 .and. ipdfnum.eq.300) then
0388 totfun =sigcross*xpp(iqcode)*xppbar(iqcode)
0389 end if
0390 end if
0391 else
0392 totfun =sigcross
0393 end if
0394
0395
0396 zfactor=2.0d0*(dotup(1,3)+dotup(2,3))/(x1*x2*ecm**2.0d0)
0397
0398
0399
0400
0401 if(totfun.lt.1.0d-16) totfun=1.0d-16
0402
0403 return
0404
0405
0406
0407
0408
0409 1005 continue
0410
0411 if(ibcstate.ne.1.and.ibcstate.ne.2) then
0412
0413
0414
0415
0416
0417
0418 do i=1,3
0419 pp(2,i)=pup(1,i)
0420 pp(3,i)=pup(2,i)
0421 pp(4,i)=pup(3,i)
0422 pp(1,i)=pup(4,i)
0423 end do
0424
0425 pp(2,4)=pup(1,5)
0426 pp(3,4)=pup(2,5)
0427 pp(4,4)=pup(3,5)
0428 pp(1,4)=pup(4,5)
0429
0430 pp(2,5)=pup(1,4)
0431 pp(3,5)=pup(2,4)
0432 pp(4,5)=pup(3,4)
0433 pp(1,5)=pup(4,4)
0434
0435
0436
0437
0438
0439 cfact=0.2578125d0
0440
0441
0442 g3=dsqrt(4.0d0*pi*alps)
0443 ampofpw = 0.0
0444 if(ibcstate.eq.3) then
0445 ampofpw=amps2_1p1()
0446 phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0447 end if
0448 if(ibcstate.eq.4) then
0449 ampofpw=amps2_3p0()
0450 phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0451 end if
0452 if(ibcstate.eq.5) then
0453 ampofpw=amps2_3p1()
0454 phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0455 end if
0456 if(ibcstate.eq.6) then
0457 ampofpw=amps2_3p2()
0458 phase=papawt*cfact*g3**8/(2.0d0**9*pi**5*dotup(1,2))
0459 end if
0460
0461 if(isubonly.eq.1) then
0462 totfun=conv*phase*ampofpw
0463 else
0464 totfun =conv*phase*ampofpw*xpp(21)*xppbar(21)/x1/x2
0465 if(ioutpdf.eq.1 .and. ipdfnum.eq.300) then
0466 totfun=conv*phase*ampofpw*xpp(21)*xppbar(21)
0467 end if
0468 end if
0469 end if
0470
0471
0472 zfactor=2.0d0*(dotup(1,3)+dotup(2,3))/(x1*x2*ecm**2.0d0)
0473
0474
0475
0476
0477 if(totfun.lt.1.0d-16) totfun=1.0d-16
0478
0479 return
0480 end