File indexing completed on 2024-04-06 12:14:04
0001
0002
0003
0004
0005
0006
0007
0008
0009 subroutine psahot(kcol,ncolp,iret)
0010
0011
0012
0013 include 'epos.inc'
0014 include 'epos.incsem'
0015 include 'epos.incems'
0016 include 'epos.incpar'
0017 double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,wplc
0018 *,wp0,wm0,s,si,smin,xmin,xp1,wpi,wmi,xp,xm,wp1,wm1,wp2,wm2
0019 *,wpq,wmq,wpi1,wmi1,pxh,pyh,pth,xmm,xp2,xg,zg,smax,xmax
0020 *,xmax1,xmin1,zp0,psutz,xpmax0,xpmin0,gb0,tmax0,tmin0,zpm,gb
0021 *,gbyj,tmax,tmin,t,gb7,x,s2,discr,qt,qt2,x1min,x1max,t1min,t1max
0022 *,t1,xq1,qq,qqmax,qqmin,pt2,xmin2,ptnew,xpmin,xpmax,xm0,psuds
0023
0024 dimension ep3(4),ey(3),bx(6)
0025 *,qmin(2),iqc(2),nqc(2),ncc(2,2),amqt(4)
0026 parameter (mjstr=20000)
0027 common /psar7/ delx,alam3p,gam3p
0028 common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
0029 common /psar30/iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
0030 common /testj/ ajeth(4),ajete(5),ajet0(7)
0031 parameter (ntim=1000)
0032 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0033 &,iorprt(ntim),jorprt(ntim),nprtj
0034 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
0035 integer icp(2),ict(2),nemis(2)
0036 integer icp1(2),icp2(2),icm1(2),icm2(2)
0037 integer jcp(nflav,2),jct(nflav,2),jcpr(nflav,2),jctr(nflav,2)
0038 common/cprtx/nprtjx,pprtx(5,2)/ciptl/iptl
0039
0040 call utpri('psahot',ish,ishini,3)
0041
0042 iret=0
0043 alpq=-(alppar+1.)/2.
0044 qqcut=q2min !????????????pt2cut
0045
0046
0047 nptl1=nptl
0048 iptl=nppr(ncolp,kcol)
0049 ip=iproj(kcol)
0050 it=itarg(kcol)
0051 do i=1,2
0052 icp(i)=icproj(i,ip)
0053 ict(i)=ictarg(i,it)
0054 enddo
0055 idpomr=idhpr(ncolp,kcol)
0056 bpomr=bhpr(ncolp,kcol)
0057
0058 q2finsave=q2fin
0059 zzzz=zparpro(kcol)+zpartar(kcol)
0060 nnn=ltarg3(it)+lproj3(ip)
0061 zz=1.+zoeinc*zzzz**2 !<-----
0062 q2fin=q2fin*zz
0063
0064
0065
0066 ajeth(idpomr+1)=ajeth(idpomr+1)+1.
0067
0068 idfpomr=idfpr(ncolp,kcol)
0069 if(ish.ge.3)write(ifch,*)'Start psahot (icp,ict):',ip,icp,it,ict
0070 *,ncolp,kcol,iptl,idpomr,idfpomr,bpomr
0071
0072
0073 if(idfpomr.eq.0)stop'idfpomr??????'
0074 if(ish.ge.3)then
0075 write(ifch,20)iptl
0076 * ,sqrt(pptl(1,iptl)**2+pptl(2,iptl)**2),pptl(3,iptl)
0077 * ,pptl(4,iptl),pptl(5,iptl)
0078 * ,istptl(iptl),ityptl(iptl)
0079 20 format(1x,i4,3x,4(e11.5,1x),i2,1x,i3)
0080 endif
0081 istptl(iptl)=31
0082
0083
0084 1 do i=1,nj
0085 ncj(1,i)=0
0086 ncj(2,i)=0
0087 enddo
0088
0089
0090 nj=0
0091 nptl=nptl1
0092 if(iremn.ge.2)then
0093 call iddeco(icp,jcp)
0094 call iddeco(ict,jct)
0095 endif
0096
0097 wp0=dsqrt(xpr(ncolp,kcol))*dexp(ypr(ncolp,kcol))*dble(engy) !???? new
0098 wm0=dsqrt(xpr(ncolp,kcol))*dexp(-ypr(ncolp,kcol))*dble(engy) !double
0099
0100
0101 amqt(1)=sqrt(sngl(xxp1pr(ncolp,kcol)**2+xyp1pr(ncolp,kcol)**2))
0102 amqt(2)=sqrt(sngl(xxp2pr(ncolp,kcol)**2+xyp2pr(ncolp,kcol)**2))
0103 amqt(3)=sqrt(sngl(xxm2pr(ncolp,kcol)**2+xym2pr(ncolp,kcol)**2))
0104 amqt(4)=sqrt(sngl(xxm1pr(ncolp,kcol)**2+xym1pr(ncolp,kcol)**2))
0105 amqpt=amqt(1)+amqt(2)+amqt(3)+amqt(4)
0106
0107 s2min=4.*q2min
0108 if(sngl(wp0*wm0).le.(sqrt(s2min)+amqpt)**2)then
0109 if(ish.ge.1)then
0110 call utmsg('psahot: insufficient pomeron mass&')
0111 write (ifch,*)'mass:',dsqrt(wp0*wm0),amqpt+sqrt(s2min)
0112 call utmsgf
0113 endif
0114 iret=1
0115 goto 16
0116 endif
0117
0118 ih=iproj(kcol)
0119 jh=itarg(kcol)
0120
0121
0122 rp=r2had(iclpro)+r2had(icltar)+slopom*log(engy**2)
0123 z=exp(-bpomr**2/(4.*.0389*rp)) !coef for rejection
0124
0125 if(z.eq.0)then
0126 write(ifch,*)'psahot : z,ih,jh ! -> ',z,ih,jh
0127 call gakli2(ih,ih)
0128 call gakli2(jh,jh)
0129 call gakli2(iptl,iptl)
0130 stop
0131 endif
0132
0133 do l=1,4
0134 bx(l)=xorptl(l,iptl)
0135 enddo
0136 bx(5)=tivptl(1,iptl)
0137 bx(6)=tivptl(2,iptl)
0138 ity=ityptl(iptl)
0139
0140 if(idpomr.eq.0)then !gg-pomeron
0141 iqq=0 !type of the hard interaction: 0 - gg, 1 - qg, 2 - gq, 3 - qq
0142 pxh=0.d0 !p_x for sh pomeron
0143 pyh=0.d0 !p_y for sh pomeron
0144 elseif(idpomr.eq.1)then !qg-pomeron
0145 iqq=1
0146 pxh=xxp1pr(ncolp,kcol)
0147 pyh=xyp1pr(ncolp,kcol)
0148 amqpt=amqpt-amqt(1)
0149 elseif(idpomr.eq.2)then !gq-pomeron
0150 iqq=2
0151 pxh=xxm2pr(ncolp,kcol)
0152 pyh=xym2pr(ncolp,kcol)
0153 amqpt=amqpt-amqt(3)
0154 elseif(idpomr.eq.3)then !qq-pomeron
0155 iqq=3
0156 pxh=xxp1pr(ncolp,kcol)+xxm2pr(ncolp,kcol)
0157 pyh=xyp1pr(ncolp,kcol)+xym2pr(ncolp,kcol)
0158 amqpt=amqpt-amqt(1)-amqt(3)
0159 else
0160 stop'unknown pomeron'
0161 endif
0162 pth=pxh**2+pyh**2
0163
0164 nj0=nj
0165 if(ish.ge.6)then
0166 write(ifch,*)'iptl,nptl,wp0,wm0,z,iqq,bx:'
0167 write(ifch,*) iptl,nptl,wp0,wm0,z,iqq,bx
0168 endif
0169
0170 s=wp0*wm0 !lc+*lc- for the semihard interaction
0171 smin=dble(s2min)+pth !mass cutoff for the hard pomeron
0172 xmin=smin/s
0173 smax=(dsqrt(s)-dble(amqpt))**2 !max mass for the hard pomeron
0174 xmax=smax/s
0175 !wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
0176 if(iLHC.eq.1.and.iqq.ne.3)then !not for val-val
0177
0178
0179
0180
0181
0182
0183
0184 xmin=max(xmin,min((sqrt(s)-dble(amqpt+ammsdd))**2/s
0185 & ,dble(max(nnn*fegypp*exp(-bpomr**2/b2xscr),zzzz)*zopinc*z)))
0186
0187 smin=xmin*s
0188 endif
0189 !wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
0190 if(smax.le.smin)then
0191 write (ifmt,*)'smax,smin',smax,smin
0192 iret=1
0193 goto 16
0194 endif
0195 xpmax0=psutz(s,smin,dble(amqpt**2)) !max x+ for the hard pomeron
0196 xpmin0=xmin/xpmax0 !min x+ for the hard pomeron
0197 xp1=wp0/dble(engy) !lc+ share for the semihard interaction
0198 xp2=wm0/dble(engy) !lc- share for the semihard interaction
0199
0200
0201
0202
0203
0204
0205 if(ish.ge.4)write(ifch,*)
0206 & 'determine LC momenta wpi,wmi for hard Pomeron'
0207
0208 if(iremn.ge.2)then
0209 if(iclpro.eq.2)then
0210 if(iabs(idptl(ip)).eq.1120)then !proj=proton
0211 nquu1=jcp(1,1)+jcp(1,2)
0212 nqud1=jcp(2,1)+jcp(2,2)
0213 elseif(iabs(idptl(ip)).eq.1220)then !proj=neutron
0214 nquu1=jcp(2,1)+jcp(2,2)
0215 nqud1=jcp(1,1)+jcp(1,2)
0216 else !to avoid flavor problem with exotic projectile (but should not happen (only gg)
0217 nquu1=0
0218 nqud1=0
0219 endif
0220 elseif(iclpro.eq.1)then
0221 if(iabs(idptl(ip)).eq.120)then
0222 nquu1=jcp(1,1)+jcp(1,2)
0223 nqud1=jcp(2,1)+jcp(2,2)
0224 else !to avoid flavor problem with exotic projectile (but should not happen (only gg)
0225 nquu1=0
0226 nqud1=0
0227 endif
0228 elseif(iclpro.eq.3)then
0229 if(iabs(idptl(ip)).eq.130)then !proj=Kch
0230 nquu1=jcp(1,1)+jcp(1,2)
0231 nqud1=jcp(3,1)+jcp(3,2)
0232 elseif(iabs(idptl(ip)).eq.230)then !proj=K0
0233 nquu1=jcp(2,1)+jcp(2,2)
0234 nqud1=jcp(3,1)+jcp(3,2)
0235 else !to avoid flavor problem with exotic projectile (but should not happen (only gg)
0236 nquu1=0
0237 nqud1=0
0238 endif
0239 else !charm
0240 if(iabs(idptl(ip)).eq.140)then
0241 nquu1=jcp(1,1)+jcp(1,2)
0242 nqud1=jcp(4,1)+jcp(4,2)
0243 elseif(iabs(idptl(ip)).eq.240)then
0244 nquu1=jcp(2,1)+jcp(2,2)
0245 nqud1=jcp(4,1)+jcp(4,2)
0246 elseif(iabs(idptl(ip)).eq.340)then
0247 nquu1=jcp(3,1)+jcp(3,2)
0248 nqud1=jcp(4,1)+jcp(4,2)
0249 else
0250 nquu1=jcp(4,1)
0251 nqud1=jcp(4,2)
0252 endif
0253 endif
0254 if(iabs(idptl(maproj+it)).eq.1220)then !targ=neutron
0255 nquu2=jct(2,1)+jct(2,2)
0256 nqud2=jct(1,1)+jct(1,2)
0257 else
0258 nquu2=jct(1,1)+jct(1,2)
0259 nqud2=jct(2,1)+jct(2,2)
0260 endif
0261 endif
0262
0263 iq1=0
0264 iq2=0
0265 wwgg=0.
0266 wwgq=0.
0267 wwqg=0.
0268 wwqq=0.
0269 wwdd=0.
0270
0271 !------------------------------------------
0272 if(iqq.eq.3)then ! val-val
0273 !------------------------------------------
0274
0275 if(ish.ge.4)write(ifch,*)'val-val'
0276 xmin1=xmin**dble(delh+.4)
0277 xmax1=xmax**dble(delh+.4)
0278 zp0=dsqrt(xmin)
0279 if(zp0.ge.1.d0)call utstop('zp0 in sem&')
0280 !........ kinematical bounds
0281 tmax0=dlog((1.d0+dsqrt(1.d0-zp0))/(1.d0-dsqrt(1.d0-zp0)))
0282 tmin0=dlog((1.d0+dsqrt(1.d0-xpmax0))
0283 * /(1.d0-dsqrt(1.d0-xpmax0)))
0284 if(iclpro.ne.4)then
0285 call psjti0(sngl(smax-pth),sqq,sqqb,1,1)
0286 call psjti0(sngl(smax-pth),sqqp,sqqpb,1,2)
0287 call psjti0(sngl(smax-pth),sqaq,sqaqb,-1,1)
0288 else
0289 call psjti0(sngl(smax-pth),sqqp,sqqpb,4,2)
0290 sqq=0.
0291 sqaq=0.
0292 endif
0293 if(iremn.ge.2)then
0294 if(nquu1.gt.nqud1.or.iclpro.ne.2)then
0295 uv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,1)
0296 dv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,2)
0297 else !if nquu<nqud => no u or no d
0298 uv1=psdfh4(sngl(zp0*xp1),q2min,0.,1,1)
0299 dv1=uv1
0300 endif
0301 if(nquu1.eq.0)uv1=0d0
0302 if(nqud1.eq.0)dv1=0d0
0303 if(nquu2.gt.nqud2)then
0304 uv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,1)
0305 dv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,2)
0306 else !if nquu<nqud => no u or no d
0307 uv2=psdfh4(sngl(zp0*xp2),q2min,0.,1,1)
0308 dv2=uv2
0309 endif
0310 if(nquu2.eq.0)uv2=0d0
0311 if(nqud2.eq.0)dv2=0d0
0312 else
0313 uv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,1)
0314 dv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,2)
0315 uv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,1)
0316 dv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,2)
0317 endif
0318 wwuu=uv1*uv2*sqq
0319 if(iclpro.eq.2)then
0320 wwdd=dv1*dv2*sqq
0321 elseif(iclpro.eq.1)then
0322 wwdd=dv1*dv2*sqaq
0323 elseif(iclpro.eq.3)then
0324 wwdd=dv1*dv2*sqqp
0325 elseif(iclpro.eq.4)then
0326 wwuu=uv1*uv2*sqqp
0327 wwdd=0.
0328 endif
0329 wwud=uv1*dv2*sqqp
0330 wwdu=dv1*uv2*sqqp
0331 wudt=wwuu+wwdd+wwud+wwdu
0332 gb0=dble(wudt)/xmax**dble(delh)/xmin**0.4*
0333 * (1.d0-zp0*xp1)**dble(-1.-alpq-alplea(iclpro))*
0334 * (1.d0-zp0*xp2)**dble(-1.-alpq-alplea(icltar))*(tmax0-tmin0)
0335 * *(1.d0-zp0)**dble(.5+alpq)*(1.d0-zp0)**dble(alpq) *5.d0
0336
0337 3 zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh+.4)) !zpm
0338 if(iclpro.ne.4)then
0339 call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
0340 call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
0341 call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
0342 else
0343 call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,4,2)
0344 sqq=0.
0345 sqaq=0.
0346 endif
0347 xpmax=psutz(s,zpm*s,dble(amqpt**2)) !max x+ for sh pomeron
0348 tmax=dlog((1.d0+dsqrt(1.d0-dsqrt(zpm)))
0349 * /(1.d0-dsqrt(1.d0-dsqrt(zpm))))
0350 tmin=dlog((1.d0+dsqrt(1.d0-xpmax))/(1.d0-dsqrt(1.d0-xpmax)))
0351
0352 t=(tmin+dble(rangen())*(tmax-tmin))
0353 xp=1.d0-((1.d0-dexp(-t))/(1.d0+dexp(-t)))**2 !x+_v
0354 xm=zpm/xp !x-_v
0355 if(xm.gt.xp.and.ish.ge.1)write(ifmt,*)'xm,xp',xm,xp
0356 gb=(1.d0-xm)**alpq*(1.d0-xp)**(.5+alpq)*(tmax-tmin)
0357 if(rangen().lt..5)then
0358 xp=xm
0359 xm=zpm/xp
0360 endif
0361
0362 if(iremn.ge.2)then
0363 if(nquu1.gt.nqud1.or.iclpro.ne.2)then
0364 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0365 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0366 else !if nquu<nqud => no u or no d
0367 uv1=psdfh4(sngl(xp*xp1),q2min,0.,1,1)
0368 dv1=uv1
0369 endif
0370 if(nquu1.eq.0)uv1=0d0
0371 if(nqud1.eq.0)dv1=0d0
0372 if(nquu2.gt.nqud2)then
0373 uv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0374 dv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0375 else !if nquu<nqud => no u or no d
0376 uv2=psdfh4(sngl(xm*xp2),q2min,0.,1,1)
0377 dv2=uv2
0378 endif
0379 if(nquu2.eq.0)uv2=0d0
0380 if(nqud2.eq.0)dv2=0d0
0381 else
0382 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0383 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0384 uv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0385 dv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0386 endif
0387
0388 wwuu=uv1*uv2*sqq
0389 if(iclpro.eq.2)then
0390 wwdd=dv1*dv2*sqq
0391 elseif(iclpro.eq.1)then
0392 wwdd=dv1*dv2*sqaq
0393 elseif(iclpro.eq.3)then
0394 wwdd=dv1*dv2*sqqp
0395 elseif(iclpro.eq.4)then
0396 wwuu=uv1*uv2*sqqp
0397 wwdd=0.
0398 endif
0399 wwud=uv1*dv2*sqqp
0400 wwdu=dv1*uv2*sqqp
0401 wudt=wwuu+wwdd+wwud+wwdu
0402 if(wudt.lt.1d-16)then
0403 if(ish.ge.1)write(ifmt,*)'No more valence quark for psahot !'
0404 write(ifch,*)'No more valence quark for psahot !'
0405 & ,ip,it,nquu1,nqud1,nquu2,nqud2
0406 iret=1
0407 goto 16
0408 endif
0409
0410 gb=gb*dble(wudt)/zpm**dble(delh+0.4)
0411 * *(1.d0-xp*xp1)**dble(-1.-alpq-alplea(iclpro))
0412 * *(1.d0-xm*xp2)**dble(-1.-alpq-alplea(icltar))/gb0
0413
0414 if(gb.gt.1.d0.and.ish.ge.1)write (ifch,*)
0415 * 'gb-qq,iclpro,zpm,xp,tmax,tmin,xpmax',
0416 * gb,iclpro,zpm,xp,tmax,tmin,xpmax
0417
0418 if(dble(rangen()).gt.gb)goto 3
0419
0420 aks=rangen()*wudt
0421 if(aks.le.wwuu)then
0422 if(iclpro.le.2)then
0423 iq1=1
0424 elseif(iclpro.eq.3)then
0425 if(iabs(idptl(ip)).eq.130)then !proj=Kch
0426 iq1=1
0427 else !proj=K0
0428 iq1=2
0429 endif
0430 else !charm
0431 if(iabs(idptl(ip)).eq.140)then
0432 iq1=1
0433 elseif(iabs(idptl(ip)).eq.240)then
0434 iq1=2
0435 elseif(iabs(idptl(ip)).eq.340)then
0436 iq1=3
0437 else
0438 iq1=4
0439 endif
0440 endif
0441 iq2=1
0442 elseif(aks.le.wwuu+wwdd)then
0443 if(iclpro.eq.2)then
0444 iq1=2
0445 elseif(iclpro.eq.1)then
0446 iq1=-2
0447 elseif(iclpro.eq.3)then
0448 iq1=-3
0449 else
0450 iq1=-4
0451 endif
0452 iq2=2
0453 elseif(aks.le.wwuu+wwdd+wwud)then
0454 if(iclpro.le.2)then
0455 iq1=1
0456 elseif(iclpro.eq.3)then
0457 if(iabs(idptl(ip)).eq.130)then !proj=Kch
0458 iq1=1
0459 else !proj=K0
0460 iq1=2
0461 endif
0462 else !charm
0463 if(iabs(idptl(ip)).eq.140)then
0464 iq1=1
0465 elseif(iabs(idptl(ip)).eq.240)then
0466 iq1=2
0467 elseif(iabs(idptl(ip)).eq.340)then
0468 iq1=3
0469 else
0470 iq1=4
0471 endif
0472 endif
0473 iq2=2
0474 else
0475 if(iclpro.eq.2)then
0476 iq1=2
0477 elseif(iclpro.eq.1)then
0478 iq1=-2
0479 elseif(iclpro.eq.3)then
0480 iq1=-3
0481 else
0482 iq1=-4
0483 endif
0484 iq2=1
0485 endif
0486
0487 wpi=xp*wp0 !lc+ for the semihard interaction
0488 wmi=xm*wm0 !lc- for the semihard interaction
0489 wp1=(wp0-wpi)
0490 wm1=(wm0-wmi)
0491 wp1=wp1*psutz(wp1*wm1,dble(amqt(2)**2),dble(amqt(4)**2))
0492 wm1=wm1-amqt(2)**2/wp1
0493
0494 !-------------------------------------
0495 else ! sea-sea val-sea sea-val
0496 !-------------------------------------
0497
0498 if(ish.ge.4)write(ifch,*)'sea-sea val-sea sea-val'
0499 xmin1=xmin**(delh-dels)
0500 xmax1=xmax**(delh-dels)
0501
0502 if(iqq.eq.0)then !rejection function normalization
0503 gb0=dlog(xpmax0/xpmin0)*(1.d0-dsqrt(xmin))**(2.*betpom) !y_soft =
0504 else
0505 tmax0=acos(dsqrt(xpmin0)) !kinematical limits for t=cos(x+-)**2
0506 tmin0=acos(dsqrt(xpmax0))
0507 if(iqq.eq.1)then
0508 uv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,1)
0509 * *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
0510 dv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,2)
0511 * *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
0512 else
0513 uv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,1)
0514 * *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
0515 dv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,2)
0516 * *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
0517 endif
0518 gb0=(1.d0-xmin/xpmax0)**betpom*dble(uv1+dv1)
0519 * *xpmin0**(-0.5+dels)
0520 * *(1.d0-xpmin0)**(0.5+alpq)*(tmax0-tmin0)
0521 if(ish.ge.6)write (ifch,*)
0522 * 'gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2',
0523 * gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2
0524 endif
0525 if(iclpro.ne.4.or.iqq.ne.1)then
0526 call psjti0(sngl(smax-pth),sj,sjb,iqq,0) !inclusive (sj) and born
0527 else
0528 call psjti0(sngl(smax-pth),sj,sjb,4,0)
0529 endif
0530 gb0=gb0*dble(sj)/xmax**delh *1.5d0 !rejection function norm.
0531 if(ish.ge.6)write (ifch,*)'gb0,sj,z',gb0,sj,z
0532
0533 if(gb0.le.0.)then
0534 write (ifmt,*)'gb0<0, smax,pth',smax,pth
0535 iret=1
0536 goto 16
0537 endif
0538
0539 ! sharing of light cone momenta between soft preevolution and
0540 ! hard interaction: ( first energy-momentum is shared according to
0541 ! f_hard(yj)~zpm**(delh-dels-1) and then rejected as
0542 ! w_rej ~sigma_hard_tot(yj) / exp(delh*yj)
0543
0544 4 continue
0545 zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh-dels)) !zpm
0546 if(iclpro.ne.4.or.iqq.ne.1)then
0547 call psjti0(sngl(zpm*s-pth),sgg,sggb,0,0)!inclusive (sgg) and born
0548 call psjti0(sngl(zpm*s-pth),sgq,sgqb,0,1)
0549 call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
0550 call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
0551 call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
0552 sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
0553 else
0554 call psjti0(sngl(zpm*s-pth),sgq,sgqb,4,0)
0555 call psjti0(sngl(zpm*s-pth),sqq,sqqb,4,1)
0556 endif
0557 xpmax=psutz(s,zpm*s,dble(amqpt**2)) !max x+ for sh pomeron
0558 xpmin=zpm/xpmax
0559 if(ish.ge.8)write (ifch,*)'zpm,xpmax,xpmin',zpm,xpmax,xpmin
0560
0561 if(iqq.eq.0)then
0562 xp=xpmin*(xpmax/xpmin)**rangen() !lc+ share for the hard interaction
0563 xm=zpm/xp !lc- share for the hard interaction
0564 else
0565 tmax=acos(dsqrt(xpmin)) !kinematical limits for t=cos(x+-)**2
0566 tmin=acos(dsqrt(xpmax))
0567 t=tmin+dble(rangen())*(tmax-tmin)
0568 xp=cos(t)**2
0569 xm=zpm/xp
0570 endif
0571
0572 if(ish.ge.8)write(ifch,*)'zpm,xp,xm,xpmax,xpmin:',
0573 * zpm,xp,xm,xpmax,xpmin
0574
0575
0576 if(iqq.eq.0)then ! --------- sea-sea -----------
0577
0578 dwwgg1=0.
0579 dwwgq1=0.
0580 dwwqg1=0.
0581 dwwqq1=0.
0582
0583 dwwgg2=0.
0584 dwwgq2=0.
0585 dwwqg2=0.
0586 dwwqq2=0.
0587
0588 wwgg=sngl((1.d0-xp)*(1.d0-xm))**betpom
0589 wwgq=sngl(1.d0-xp)**betpom*EsoftQZero(sngl(xm))
0590 wwqg=sngl(1.d0-xm)**betpom*EsoftQZero(sngl(xp))
0591 wwqq=EsoftQZero(sngl(xp))*EsoftQZero(sngl(xm))
0592
0593 if(idfpomr.eq.1)then
0594 rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(zpm))
0595 wwgg=(wwgg*z**(rp/rh)/rh+dwwgg1+dwwgg2)
0596 * *(r2had(iclpro)+r2had(icltar))/z
0597 wwgq=(wwgq*z**(rp/rh)/rh+dwwgq1+dwwgq2)
0598 * *(r2had(iclpro)+r2had(icltar))/z
0599 wwqg=(wwqg*z**(rp/rh)/rh+dwwqg1+dwwqg2)
0600 * *(r2had(iclpro)+r2had(icltar))/z
0601 wwqq=(wwqq*z**(rp/rh)/rh+dwwqq1+dwwqq2)
0602 * *(r2had(iclpro)+r2had(icltar))/z
0603 elseif(idfpomr.eq.2)then
0604 rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(zpm))
0605 wwgg=wwgg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0606 wwqg=wwqg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0607 wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0608 wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0609 elseif(idfpomr.eq.3)then
0610 rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(zpm))
0611 wwgg=wwgg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0612 wwqg=wwqg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0613 wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0614 wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0615 elseif(idfpomr.eq.4)then
0616 rh=alam3p-slopom*sngl(dlog(zpm))
0617 wwgg=wwgg/rh*alam3p*z**(rp/rh-1.)
0618 wwqg=wwqg/rh*alam3p*z**(rp/rh-1.)
0619 wwgq=wwgq/rh*alam3p*z**(rp/rh-1.)
0620 wwqq=wwqq/rh*alam3p*z**(rp/rh-1.)
0621 else
0622 stop'psahot-idfpomr????'
0623 endif
0624
0625 wwgg=wwgg*sgg*(1.-glusea)**2
0626 wwgq=wwgq*sgq*(1.-glusea)*glusea
0627 wwqg=wwqg*sgq*(1.-glusea)*glusea
0628 wwqq=wwqq*sqq*glusea**2
0629 gbyj=dlog(xpmax/xpmin)*dble(wwgg+wwgq+wwqg+wwqq)
0630 wpi=wp0*xp !lc+ for the hard interaction
0631 wmi=wm0*xm !lc+-for the hard interaction
0632 gbyj=gbyj/zpm**dble(delh)/gb0 !rejection fu
0633 if(gbyj.ge.1.d0.and.ish.ge.2)write(ifmt,*)'gbyj',gbyj
0634 if(dble(rangen()).gt.gbyj)goto 4 !rejection
0635 wp1=wp0-wpi
0636 wm1=wm0-wmi
0637 call pslcsh(wp1,wm1,wp2,wm2,amqt,dble(amqpt))
0638
0639 else ! --------- val-sea sea-val -----------
0640
0641 dwwg=0.
0642 dwwq=0.
0643
0644 wwgq=sngl(1.d0-xm)**betpom
0645 wwqq=EsoftQZero(sngl(xm))
0646 if(idfpomr.eq.1)then
0647 rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(xm))
0648 wwgq=(wwgq*z**(rp/rh)/rh+dwwg)
0649 * *(r2had(iclpro)+r2had(icltar))/z
0650 wwqq=(wwqq*z**(rp/rh)/rh+dwwq)
0651 * *(r2had(iclpro)+r2had(icltar))/z
0652 else !tp071031 not used anymore
0653 if(iqq.eq.1)then
0654 rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(xm))
0655 wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0656 wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
0657 else
0658 rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(xm))
0659 wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0660 wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
0661 endif
0662 endif
0663
0664 wwgq=wwgq*sgq*(1.-glusea)*sngl(xp)**(-0.5+dels)
0665 * *sngl(1.d0-xp)**(0.5+alpq)
0666 wwqq=wwqq*sqq*glusea*sngl(xp)**(-0.5+dels)
0667 * *sngl(1.d0-xp)**(0.5+alpq)
0668
0669 if(iqq.eq.1)then !valence quark-gluon hard interaction
0670 if(iremn.ge.2)then
0671 if(nquu1.gt.nqud1.or.iclpro.ne.2)then
0672 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0673 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0674 else !if nquu<nqud => no u or no d
0675 uv1=psdfh4(sngl(xp*xp1),q2min,0.,1,1)
0676 dv1=uv1
0677 endif
0678 if(nquu1.eq.0)uv1=0d0
0679 if(nqud1.eq.0)dv1=0d0
0680 else
0681 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
0682 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
0683 endif
0684 if(uv1+dv1.lt.1d-16)then
0685 if(ish.ge.1)write(ifmt,*)'No more valence quark for psahot !'
0686 write(ifch,*)'No more valence quark in projectile for psahot !'
0687 & ,ip,nquu1,nqud1
0688 iret=1
0689 goto 16
0690 endif
0691 wpi=wp0*xp !lc+ for the hard interaction
0692 wmi=wm0*xm !lc+-for the hard interaction
0693 aks=rangen()
0694 if(aks.le.uv1/(uv1+dv1))then
0695 if(iclpro.le.2)then
0696 iq1=1
0697 elseif(iclpro.eq.3)then
0698 if(iabs(idptl(ip)).eq.130)then !proj=Kch
0699 iq1=1
0700 else !proj=K0
0701 iq1=2
0702 endif
0703 else !charm
0704 if(iabs(idptl(ip)).eq.140)then
0705 iq1=1
0706 elseif(iabs(idptl(ip)).eq.240)then
0707 iq1=2
0708 elseif(iabs(idptl(ip)).eq.340)then
0709 iq1=3
0710 else
0711 iq1=4
0712 endif
0713 endif
0714 else
0715 if(iclpro.eq.2)then
0716 iq1=2
0717 elseif(iclpro.eq.1)then
0718 iq1=-2
0719 elseif(iclpro.eq.3)then
0720 iq1=-3
0721 else
0722 iq1=-4
0723 endif
0724 endif
0725 gbyj=dble((wwgq+wwqq)*(uv1+dv1))*(1.d0-xp*xp1)**
0726 * (-1.-alpq-alplea(iclpro))
0727 else !gluon-valence quark hard interaction
0728 xm=xp
0729 xp=zpm/xm
0730 if(iremn.ge.2)then
0731 if(nquu2.gt.nqud2)then
0732 uv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0733 dv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0734 else !if nquu<nqud => no u or no d
0735 uv1=psdfh4(sngl(xm*xp2),q2min,0.,1,1)
0736 dv1=uv1
0737 endif
0738 if(nquu2.eq.0)uv1=0d0
0739 if(nqud2.eq.0)dv1=0d0
0740 else
0741 uv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
0742 dv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
0743 endif
0744 if(uv1+dv1.lt.1d-16)then
0745 if(ish.ge.1)write(ifmt,*)'No more valence quark for psahot !'
0746 write(ifch,*)'No more valence quark in target for psahot !'
0747 & ,it,nquu2,nqud2
0748 iret=1
0749 goto 16
0750 endif
0751 wpi=wp0*xp !lc+ for the hard interaction
0752 wmi=wm0*xm !lc+-for the hard interaction
0753 aks=rangen()
0754 if(aks.le.uv1/(uv1+dv1))then
0755 iq2=1
0756 else
0757 iq2=2
0758 endif
0759 gbyj=dble(wwgq+wwqq)*dble(uv1+dv1)*
0760 * (1.d0-xm*xp2)**(-1.-alpq-alplea(icltar))
0761 endif
0762
0763 gbyj=gbyj*(tmax-tmin)/zpm**delh /gb0 /2.1d0 !rejection
0764 if(ish.ge.6)write (ifch,*)
0765 * 'gbyj,zpm,xp,tmax,tmin,xpmax,xpmin',
0766 * gbyj,zpm,xp,tmax,tmin,xpmax,xpmin
0767 if(dble(rangen()).gt.gbyj)goto 4
0768
0769 wp1=wp0-wpi
0770 wm1=wm0-wmi
0771 if(ish.ge.8)write (ifch,*)'q_sea mass check',wp1*wm1,amqpt
0772 if(iqq.eq.1)then !valence quark-gluon hard interaction
0773 amq1=amqt(3)**2
0774 s24=(amqt(2)+amqt(4))**2
0775 else
0776 amq1=amqt(1)**2
0777 s24=(amqt(2)+amqt(4))**2
0778 xp=xm
0779 xm=zpm/xp
0780 endif
0781 x1max=psutz(wp1*wm1,dble(amq1),dble(s24))
0782 x1min=dble(amq1)/x1max/wp1/wm1
0783 t1min=(1.d0/x1max-1.d0)
0784 t1max=(1.d0/x1min-1.d0)
0785 5 t1=t1min*(t1max/t1min)**dble(rangen())
0786 if(ish.ge.8)write (ifch,*)'t1,t1min,t1max',t1,t1min,t1max
0787 xq1=1.d0/(1.d0+t1)
0788 if(dble(rangen()).gt.(xq1*(1.d0-xq1))**(1.+(-alpqua)))goto 5
0789 if(iqq.eq.1)then !valence quark-gluon hard interacti
0790 wm2=wm1*(1.d0-xq1)
0791 wm1=wm1*xq1
0792 wp1=wp1-dble(amq1)/wm1
0793 if(ish.ge.8)write (ifch,*)'q_sea+ mass check',
0794 * wp1*wm2,s24
0795 wp1=wp1*psutz(wp1*wm2,dble(amqt(2)**2),dble(amqt(4)**2))
0796 wm2=wm2-dble(amqt(2))**2/wp1
0797 else
0798 wp2=wp1*(1.d0-xq1)
0799 wp1=wp1*xq1
0800 wm1=wm1-dble(amq1)/wp1
0801 if(ish.ge.8)write (ifch,*)'q_sea- mass check',
0802 * wp2*wm1,s24
0803 wm1=wm1*psutz(wm1*wp2,dble(amqt(4)**2),dble(amqt(2)**2))
0804 wp2=wp2-amqt(4)**2/wm1
0805 endif
0806
0807 endif ! -------------------
0808
0809 !-------------------------------
0810 endif
0811 !-------------------------------
0812
0813
0814
0815
0816
0817 if(ish.ge.4)write(ifch,*)
0818 & 'flavor and momenta of end partons of the hard Pomeron'
0819
0820 6 continue
0821 wpi1=wpi
0822 wmi1=wmi
0823 nj=nj0 !initialization for the number of jets
0824 if(ish.ge.8)write (ifch,*)'5-ww,smin',wpi*wmi,smin
0825
0826 rrr=rangen()
0827 jqq=0
0828
0829 call iddeco(icp,jcp) !save valence quarks in jcp and jct
0830 call iddeco(ict,jct)
0831 if(iremn.ge.2)then
0832 do j=1,2
0833 do n=1,nrflav
0834 jcpr(n,j)=jcpref(n,j,ip) !remnant sea quarks in jcpr and jctr
0835 jctr(n,j)=jctref(n,j,it)
0836 enddo
0837 do n=nrflav+1,nflav
0838 jcpr(n,j)=0
0839 jctr(n,j)=0
0840 enddo
0841 enddo
0842 endif
0843 do i=1,2
0844 icp1(i)=0
0845 icp2(i)=0
0846 icm1(i)=0
0847 icm2(i)=0
0848 enddo
0849
0850 iret1=0
0851 iret2=0
0852
0853
0854
0855
0856
0857
0858
0859
0860
0861
0862 if(iqq.eq.1.or.iqq.eq.2)then
0863 if(rrr.lt.wwqq/(wwgq+wwqq))jqq=1
0864 elseif(iqq.eq.0)then
0865 if(rrr.lt.wwqg/(wwgg+wwqg+wwgq+wwqq))then
0866 jqq=1
0867 elseif(rrr.lt.(wwqg+wwgq)/(wwgg+wwqg+wwgq+wwqq))then
0868 jqq=2
0869 elseif(rrr.lt.(wwqg+wwgq+wwqq)/(wwgg+wwqg+wwgq+wwqq))then
0870 jqq=3
0871 endif
0872 endif
0873
0874 if(iLHC.eq.1)then
0875
0876 if((iqq-1)*(iqq-3).eq.0)then !valence quark on projectile side
0877
0878 iqc(1)=iq1 !proj=particle
0879 if(iabs(idptl(ip)).eq.1220)iqc(1)=3-iq1 !proj=neutron
0880 if(idptl(ip).lt.0)iqc(1)=-iqc(1) !proj=antiparticle
0881 ifl1=iabs(iqc(1))
0882
0883
0884
0885 if(iabs(idptl(ip)).lt.1000)then
0886 if(iqc(1).gt.0.and.idp1pr(ncolp,kcol).ne.2)then
0887 idp2pr(ncolp,kcol)=idp1pr(ncolp,kcol)
0888 idp1pr(ncolp,kcol)=2
0889 elseif(iqc(1).lt.0.and.idp2pr(ncolp,kcol).ne.2)then
0890 idp1pr(ncolp,kcol)=idp2pr(ncolp,kcol)
0891 idp2pr(ncolp,kcol)=2
0892 endif
0893 endif
0894
0895
0896 idsp=100
0897 if(iqc(1).gt.0)then
0898 icp1(1)=ifl1
0899 if(idp1pr(ncolp,kcol).ne.2)
0900 & call utstop("psahot: Problem with SE (1)!&")
0901 else
0902 icp2(2)=ifl1
0903 if(idp2pr(ncolp,kcol).ne.2)
0904 & call utstop("psahot: Problem with SE (2)!&")
0905 endif
0906 else
0907 idsp=idsppr(ncolp,kcol)
0908 endif
0909
0910 if((iqq-2)*(iqq-3).eq.0)then !valence quark on target side
0911
0912 iqc(2)=iq2 !tar=particle (can not be pion or kaon !)
0913 if(iabs(idptl(maproj+it)).eq.1220)iqc(2)=3-iq2 !targ=neutron
0914 if(idptl(maproj+it).lt.0)iqc(2)=-iqc(2) !targ=antiparticle
0915 ifl2=iabs(iqc(2))
0916
0917
0918 idst=100
0919 if(iqc(2).gt.0)then
0920 icm1(1)=ifl2
0921 if(idm1pr(ncolp,kcol).ne.2)
0922 & call utstop("psahot: Problem with SE (3)!&")
0923 else
0924 icm2(2)=ifl2
0925 if(idm2pr(ncolp,kcol).ne.2)
0926 & call utstop("psahot: Problem with SE (4)!&")
0927 endif
0928 else
0929 idst=idstpr(ncolp,kcol)
0930 endif
0931
0932
0933
0934 iret=0
0935
0936 call fstrfl(jcpr,jctr,jcp,jct,icp1,icp2,icm1,icm2
0937 * ,idp1pr(ncolp,kcol),idp2pr(ncolp,kcol)
0938 * ,idm1pr(ncolp,kcol),idm2pr(ncolp,kcol)
0939 * ,idsp,idst,iret)
0940
0941 if(iret.ne.0)goto 16
0942
0943
0944
0945 if((iqq-1)*(iqq-3).eq.0)then !valence quark on projectile side
0946
0947
0948
0949 if(iqc(1).gt.0)then
0950 ifl=idtra(icp2,1,0,2)
0951 elseif(iqc(1).lt.0)then
0952 ifl=idtra(icp1,1,0,2)
0953 else
0954 call utstop('No quark for hard Pomeron+ in psahot!&')
0955 endif
0956
0957 if(ish.ge.5)write(ifch,*)'flavor vq+,sqb+',iqc(1)
0958 & ,ifl
0959
0960 nj=nj+1
0961
0962 if(abs(ifl).eq.3)then
0963 iqj(nj)=ifl*4/3
0964 elseif(abs(ifl).ge.4)then
0965 iqj(nj)=ifl*10
0966 else
0967 iqj(nj)=ifl
0968 endif
0969 ioj(nj)=7
0970
0971 eqj(1,nj)=.5*sngl(wp1+dble(amqt(2))**2/wp1)
0972 eqj(2,nj)=wp1-eqj(1,nj)
0973 eqj(3,nj)=xxp2pr(ncolp,kcol)
0974 eqj(4,nj)=xyp2pr(ncolp,kcol)
0975 if(ish.ge.8)write (ifch,*)'q_v+ mass check',(eqj(1,nj)-
0976 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
0977 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
0978 ncc(1,1)=nj
0979 ncc(2,1)=0
0980
0981 else !gluon (soft sea) parton on projectile side
0982
0983
0984
0985
0986 iflq=idtra(icp1,1,0,2)
0987 iflqb=idtra(icp2,1,0,2)
0988
0989
0990 if(ish.ge.5)write(ifch,*)'flavor sq+,sqb+',iflq,iflqb
0991
0992 nj=nj+1
0993
0994 if(abs(iflqb).eq.3)then
0995 iqj(nj)=iflqb*4/3
0996 elseif(abs(iflqb).ge.4)then
0997 iqj(nj)=iflqb*10
0998 else
0999 iqj(nj)=iflqb
1000 endif
1001 ioj(nj)=7
1002 if(abs(iflq).eq.3)then
1003 iqj(nj+1)=iflq*4/3
1004 elseif(abs(iflq).ge.4)then
1005 iqj(nj+1)=iflq*10
1006 else
1007 iqj(nj+1)=iflq
1008 endif
1009 ioj(nj+1)=7
1010
1011 eqj(1,nj)=.5*sngl(wp1+dble(amqt(1))**2/wp1)
1012 eqj(2,nj)=wp1-eqj(1,nj)
1013 eqj(3,nj)=xxp1pr(ncolp,kcol)
1014 eqj(4,nj)=xyp1pr(ncolp,kcol)
1015 if(ish.ge.8)write (ifch,*)'q_s1+ mass check',(eqj(1,nj)-
1016 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1017 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1018 eqj(1,nj+1)=.5*sngl(wp2+dble(amqt(2))**2/wp2)
1019 eqj(2,nj+1)=wp2-eqj(1,nj+1)
1020 eqj(3,nj+1)=xxp2pr(ncolp,kcol)
1021 eqj(4,nj+1)=xyp2pr(ncolp,kcol)
1022 nj=nj+1
1023 if(ish.ge.8)write (ifch,*)'q_s2+ mass check',(eqj(1,nj)-
1024 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1025 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1026
1027
1028 if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.2)then
1029 iqc(1)=0
1030 ncc(1,1)=nj-1
1031 ncc(2,1)=nj
1032
1033 else
1034
1035 iqc(1)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1036 if(1.-1./(wp0*(wmi1-smin/wpi1)).le.xp)goto 1
1037 7 zg=1d0-dble(rangen())*(1.d0-xp)
1038 if(ish.ge.8)write (ifch,*)'6-zg,xp',zg,xp
1039 if(dble(rangen()).gt.zg**dels*((1.d0-xp/zg)/(1.d0-xp))
1040 * **betpom)goto 7
1041 xg=xp/zg
1042 wpq=wp0*(xg-xp)
1043 wmq=1.d0/wpq
1044 wmi1=wmi1-wmq
1045 if(wmi1*wpi1.le.smin)goto 1
1046
1047
1048
1049 nj=nj+1
1050 if(iabs(iqc(1)).eq.3)then
1051 iqj(nj)=-iqc(1)*4/3
1052 elseif(iabs(iqc(1)).ge.3)then
1053 iqj(nj)=-iqc(1)*10
1054 else
1055 iqj(nj)=-iqc(1)
1056 endif
1057 ioj(nj)=-7
1058 eqj(1,nj)=.5*wmq
1059 eqj(2,nj)=-eqj(1,nj)
1060 eqj(3,nj)=0.
1061 eqj(4,nj)=0.
1062 if(ish.ge.8)write (ifch,*)'q_s3+ mass check',eqj(1,nj)**2-
1063 * eqj(2,nj)**2-eqj(3,nj)**2-eqj(4,nj)**2
1064 if(iqc(1).gt.0)then
1065 ncj(1,nj)=nj-1
1066 ncj(1,nj-1)=nj
1067 ncj(2,nj)=0
1068 ncj(2,nj-1)=0
1069 ncc(1,1)=nj-2
1070 ncc(2,1)=0
1071 else
1072 ncj(1,nj)=nj-2
1073 ncj(1,nj-2)=nj
1074 ncj(2,nj)=0
1075 ncj(2,nj-2)=0
1076 ncc(1,1)=nj-1
1077 ncc(2,1)=0
1078 endif
1079 endif
1080 endif
1081
1082
1083
1084 if((iqq-2)*(iqq-3).eq.0)then !valence quark on target side
1085
1086
1087 if(iqc(2).gt.0)then
1088 ifl=idtra(icm2,1,0,2)
1089 elseif(iqc(2).lt.0)then
1090 ifl=idtra(icm1,1,0,2)
1091 else
1092 call utstop('No quark for hard Pomeron- in psahot!&')
1093 endif
1094
1095 if(ish.ge.5)write(ifch,*)'flavor vq-,sqb-',iqc(2)
1096 & ,ifl
1097
1098 nj=nj+1
1099
1100 if(abs(ifl).eq.3)then
1101 iqj(nj)=ifl*4/3
1102 elseif(abs(ifl).ge.4)then
1103 iqj(nj)=ifl*10
1104 else
1105 iqj(nj)=ifl
1106 endif
1107 ioj(nj)=8
1108
1109 eqj(1,nj)=.5*sngl(wm1+dble(amqt(4))**2/wm1)
1110 eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1111 eqj(3,nj)=xxm1pr(ncolp,kcol)
1112 eqj(4,nj)=xym1pr(ncolp,kcol)
1113 if(ish.ge.8)write (ifch,*)'q_v- mass check',(eqj(1,nj)
1114 * +eqj(2,nj))*(eqj(1,nj)-eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1115 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1116 ncc(1,2)=nj
1117 ncc(2,2)=0
1118
1119 else !gluon (soft sea) parton on target side
1120
1121
1122
1123 iflq=idtra(icm1,1,0,2)
1124 iflqb=idtra(icm2,1,0,2)
1125
1126
1127 if(ish.ge.5)write(ifch,*)'flavor sq-,sqb-',iflq,iflqb
1128
1129 nj=nj+1
1130
1131 if(abs(iflqb).eq.3)then
1132 iqj(nj)=iflqb*4/3
1133 elseif(abs(iflqb).ge.4)then
1134 iqj(nj)=iflqb*10
1135 else
1136 iqj(nj)=iflqb
1137 endif
1138 ioj(nj)=8
1139 if(abs(iflq).eq.3)then
1140 iqj(nj+1)=iflq*4/3
1141 elseif(abs(iflq).ge.4)then
1142 iqj(nj+1)=iflq*10
1143 else
1144 iqj(nj+1)=iflq
1145 endif
1146 ioj(nj+1)=8
1147
1148 eqj(1,nj)=.5*sngl(wm1+dble(amqt(3))**2/wm1)
1149 eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1150 eqj(3,nj)=xxm2pr(ncolp,kcol)
1151 eqj(4,nj)=xym2pr(ncolp,kcol)
1152 if(ish.ge.8)write (ifch,*)'q_s1- mass check',(eqj(1,nj)-
1153 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1154 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1155 eqj(1,nj+1)=.5*sngl(wm2+dble(amqt(4))**2/wm2)
1156 eqj(2,nj+1)=eqj(1,nj+1)-wm2
1157 eqj(3,nj+1)=xxm1pr(ncolp,kcol)
1158 eqj(4,nj+1)=xym1pr(ncolp,kcol)
1159 nj=nj+1
1160 if(ish.ge.8)write (ifch,*)'q_s2- mass check',(eqj(1,nj)-
1161 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1162 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1163
1164
1165 if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.1)then
1166 iqc(2)=0
1167 ncc(1,2)=nj-1
1168 ncc(2,2)=nj
1169
1170 else
1171
1172 iqc(2)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1173 if(1.-1./(wm0*(wpi1-smin/wmi1)).le.xm)goto 1
1174 8 zg=1.d0-dble(rangen())*(1.d0-xm)
1175 if(ish.ge.8)write (ifch,*)'7-zg,xm',zg,xm
1176 if(rangen().gt.zg**dels*((1.d0-xm/zg)/(1.d0-xm))**betpom)
1177 * goto 8
1178 xg=xm/zg
1179 wmq=wm0*(xg-xm)
1180 wpq=1.d0/wmq
1181 wpi1=wpi1-wpq
1182 if(wmi1*wpi1.le.smin)goto 1
1183
1184
1185 nj=nj+1
1186 if(iabs(iqc(2)).eq.3)then
1187 iqj(nj)=-iqc(2)*4/3
1188 elseif(iabs(iqc(2)).ge.4)then
1189 iqj(nj)=-iqc(2)*10
1190 else
1191 iqj(nj)=-iqc(2)
1192 endif
1193 ioj(nj)=-8
1194
1195 eqj(1,nj)=.5*sngl(wpq)
1196 eqj(2,nj)=eqj(1,nj)
1197 eqj(3,nj)=0.
1198 eqj(4,nj)=0.
1199 if(ish.ge.8)write (ifch,*)'q_s3- mass check',(eqj(1,nj)-
1200 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1201 if(iqc(2).gt.0)then
1202 ncj(1,nj)=nj-1
1203 ncj(1,nj-1)=nj
1204 ncj(2,nj)=0
1205 ncj(2,nj-1)=0
1206 ncc(1,2)=nj-2
1207 ncc(2,2)=0
1208 else
1209 ncj(1,nj)=nj-2
1210 ncj(1,nj-2)=nj
1211 ncj(2,nj)=0
1212 ncj(2,nj-2)=0
1213 ncc(1,2)=nj-1
1214 ncc(2,2)=0
1215 endif
1216 endif
1217 endif
1218
1219 else
1220
1221
1222
1223 if((iqq-1)*(iqq-3).eq.0)then !valence quark on projectile side
1224
1225 iqc(1)=iq1 !proj=particle
1226 if(iabs(idptl(ip)).eq.1220)iqc(1)=3-iq1 !proj=neutron
1227 if(idptl(ip).lt.0)iqc(1)=-iqc(1) !proj=antiparticle
1228 ifl1=iabs(iqc(1))
1229
1230 nj=nj+1
1231
1232 if(iqc(1).gt.0)then
1233 if(iremn.ne.0)then
1234 if(iremn.eq.3)then
1235 call idsufl3(ifl1,1,jcp) !remove valence quark from remnant
1236 ifl=idrafl(iclpro,jcpr,2,'v',1,iret3) !take sea antiquark from jcpr
1237 elseif(iremn.eq.2)then
1238 call idsufl3(ifl1,1,jcp) !remove valence quark from remnant
1239 ifl=idrafl(iclpro,jcpr,2,'s',0,iret3) !take sea antiquark freely
1240 call idsufl3(ifl,2,jcpr) !remove sea antiquark from remnant
1241 else
1242 call idsufl(ifl1,1,jcp,iret1) !remove valence quark from remnant
1243 if(iret1.ne.1)then
1244 ifl=idrafl(iclpro,jcp,2,'s',1,iret2) !take sea antiquark
1245 endif
1246 endif
1247 else
1248 ifl=ifl1
1249 endif
1250 elseif(iqc(1).lt.0)then
1251 if(iremn.ne.0)then
1252 if(iremn.eq.3)then
1253 call idsufl3(ifl1,2,jcp)
1254 ifl=idrafl(iclpro,jcpr,1,'v',1,iret3) !take sea quark
1255 elseif(iremn.eq.2)then
1256 call idsufl3(ifl1,2,jcp)
1257 ifl=idrafl(iclpro,jcpr,1,'s',0,iret3) !take sea quark
1258 call idsufl3(ifl,1,jcpr) !remove sea quark from remnant
1259 else
1260 call idsufl(ifl1,2,jcp,iret1)
1261 if(iret1.ne.1)then
1262 ifl=idrafl(iclpro,jcp,1,'s',1,iret2)
1263 endif
1264 endif
1265 else
1266 ifl=ifl1
1267 endif
1268 else
1269 call utstop('No quark for hard Pomeron+ in psahot!&')
1270 endif
1271
1272 if(iret1.eq.1.or.iret2.eq.1)then
1273 ifl=ifl1
1274 if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 1)'
1275 if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 1)'
1276 call iddeco(icp,jcp)
1277 iret1=0
1278 iret2=0
1279 endif
1280
1281 if(ish.ge.5)write(ifch,*)'flavor vq+,sqb+',iqc(1)
1282 & ,-isign(ifl,iqc(1))
1283
1284 if(ifl.eq.3)then
1285 iqj(nj)=-isign(ifl,iqc(1))*4/3
1286 elseif(ifl.eq.4)then
1287 iqj(nj)=-isign(ifl,iqc(1))/4
1288 else
1289 iqj(nj)=-isign(ifl,iqc(1))
1290 endif
1291 eqj(1,nj)=.5*sngl(wp1+dble(amqt(2))**2/wp1)
1292 eqj(2,nj)=wp1-eqj(1,nj)
1293 eqj(3,nj)=xxp2pr(ncolp,kcol)
1294 eqj(4,nj)=xyp2pr(ncolp,kcol)
1295 if(ish.ge.8)write (ifch,*)'q_v+ mass check',(eqj(1,nj)-
1296 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1297 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1298 ncc(1,1)=nj
1299 ncc(2,1)=0
1300
1301 else
1302
1303 nj=nj+1
1304 if(idfpomr.lt.3.and.iremn.ne.0)then
1305 if(iremn.eq.3)then
1306 iflq=idrafl(iclpro,jcpr,1,'v',1,iret3) !take sea quark
1307 iflqb=idrafl(iclpro,jcpr,2,'v',1,iret3) !take sea antiquark
1308 elseif(iremn.eq.2)then
1309 iflq=idrafl(iclpro,jcpr,1,'s',0,iret3) !take sea quark
1310 call idsufl3(iflq,1,jcpr) !remove sea quark from remnant
1311 iflqb=idrafl(iclpro,jcpr,2,'s',0,iret3) !take sea antiquark
1312 call idsufl3(iflqb,2,jcpr) !remove sea antiquark from remnant
1313 else
1314 iflq=idrafl(iclpro,jcp,1,'s',1,iret1) !Choose flavor of sea quark.
1315 if(iret1.ne.1)iflqb=idrafl(iclpro,jcp,2,'s',1,iret2) !Choose antiquark.
1316 endif
1317 else
1318 iflq=idrafl(iclpro,jcp,0,'s',0,iret1) !Choose flavor of sea quark.
1319 iflqb=iflq !antiquark=quark (vertex end)
1320 endif
1321
1322 if(iret1.eq.1.or.iret2.eq.1)then
1323 iflqb=iflq
1324 if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 2)'
1325 if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 2)'
1326 call iddeco(icp,jcp)
1327 iret1=0
1328 iret2=0
1329 endif
1330
1331 if(ish.ge.5)write(ifch,*)'flavor sq+,sqb+',iflq,-iflqb
1332
1333 if(iflqb.eq.3)then
1334 iqj(nj)=-iflqb*4/3
1335 elseif(iflqb.eq.4)then
1336 iqj(nj)=-iflqb/4
1337 else
1338 iqj(nj)=-iflqb
1339 endif
1340 ioj(nj)=7
1341 if(iflq.eq.3)then
1342 iqj(nj+1)=iflq*4/3
1343 elseif(iflq.eq.4)then
1344 iqj(nj+1)=iflq/4
1345 else
1346 iqj(nj+1)=iflq
1347 endif
1348 ioj(nj+1)=7
1349
1350 eqj(1,nj)=.5*sngl(wp1+dble(amqt(1))**2/wp1)
1351 eqj(2,nj)=wp1-eqj(1,nj)
1352 eqj(3,nj)=xxp1pr(ncolp,kcol)
1353 eqj(4,nj)=xyp1pr(ncolp,kcol)
1354 if(ish.ge.8)write (ifch,*)'q_s1+ mass check',(eqj(1,nj)-
1355 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1356 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1357 eqj(1,nj+1)=.5*sngl(wp2+dble(amqt(2))**2/wp2)
1358 eqj(2,nj+1)=wp2-eqj(1,nj+1)
1359 eqj(3,nj+1)=xxp2pr(ncolp,kcol)
1360 eqj(4,nj+1)=xyp2pr(ncolp,kcol)
1361 nj=nj+1
1362 if(ish.ge.8)write (ifch,*)'q_s2+ mass check',(eqj(1,nj)-
1363 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1364 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1365
1366
1367 if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.2)then
1368 iqc(1)=0
1369 ncc(1,1)=nj-1
1370 ncc(2,1)=nj
1371
1372 else
1373
1374 iqc(1)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1375 if(1.-1./(wp0*(wmi1-smin/wpi1)).le.xp)goto 1
1376 17 zg=1d0-dble(rangen())*(1.d0-xp)
1377 if(ish.ge.8)write (ifch,*)'6-zg,xp',zg,xp
1378 if(dble(rangen()).gt.zg**dels*((1.d0-xp/zg)/(1.d0-xp))
1379 * **betpom)goto 17
1380 xg=xp/zg
1381 wpq=wp0*(xg-xp)
1382 wmq=1.d0/wpq
1383 wmi1=wmi1-wmq
1384 if(wmi1*wpi1.le.smin)goto 1
1385
1386
1387
1388 nj=nj+1
1389 if(iabs(iqc(1)).eq.3)then
1390 iqj(nj)=-iqc(1)*4/3
1391 else
1392 iqj(nj)=-iqc(1)
1393 endif
1394 ioj(nj)=-7
1395 eqj(1,nj)=.5*wmq
1396 eqj(2,nj)=-eqj(1,nj)
1397 eqj(3,nj)=0.
1398 eqj(4,nj)=0.
1399 if(ish.ge.8)write (ifch,*)'q_s3+ mass check',eqj(1,nj)**2-
1400 * eqj(2,nj)**2-eqj(3,nj)**2-eqj(4,nj)**2
1401 if(iqc(1).gt.0)then
1402 ncj(1,nj)=nj-1
1403 ncj(1,nj-1)=nj
1404 ncj(2,nj)=0
1405 ncj(2,nj-1)=0
1406 ncc(1,1)=nj-2
1407 ncc(2,1)=0
1408 else
1409 ncj(1,nj)=nj-2
1410 ncj(1,nj-2)=nj
1411 ncj(2,nj)=0
1412 ncj(2,nj-2)=0
1413 ncc(1,1)=nj-1
1414 ncc(2,1)=0
1415 endif
1416 endif
1417 endif
1418
1419 if((iqq-2)*(iqq-3).eq.0)then
1420 iqc(2)=iq2 !tar=particle (can not be pion or kaon !)
1421 if(iabs(idptl(maproj+it)).eq.1220)iqc(2)=3-iq2 !targ=neutron
1422 if(idptl(maproj+it).lt.0)iqc(2)=-iqc(2) !targ=antiparticle
1423 ifl2=iabs(iqc(2))
1424
1425 nj=nj+1
1426 if(iqc(2).gt.0)then
1427 if(iremn.ne.0)then
1428 if(iremn.eq.3)then
1429 call idsufl3(ifl2,1,jct) !remove valence quark from remnant
1430 ifl=idrafl(icltar,jctr,2,'v',1,iret3) !take sea antiquark from jctr
1431 elseif(iremn.eq.2)then
1432 call idsufl3(ifl2,1,jct) !remove valence quark from remnant
1433 ifl=idrafl(icltar,jctr,2,'s',0,iret3) !take sea antiquark freely
1434 call idsufl3(ifl,2,jctr) !remove sea antiquark from remnant
1435 else
1436 call idsufl(ifl2,1,jct,iret1) !remove valence quark from remnant
1437 if(iret1.ne.1)then
1438 ifl=idrafl(icltar,jct,2,'s',1,iret2) !take sea antiquark
1439 endif
1440 endif
1441 else
1442 ifl=ifl2
1443 endif
1444 elseif(iqc(2).lt.0)then
1445 if(iremn.ne.0)then
1446 if(iremn.eq.3)then
1447 call idsufl3(ifl2,2,jct)
1448 ifl=idrafl(icltar,jctr,1,'v',1,iret3) !take sea quark
1449 elseif(iremn.eq.2)then
1450 call idsufl3(ifl2,2,jct)
1451 ifl=idrafl(icltar,jctr,1,'s',0,iret3) !take sea quark
1452 call idsufl3(ifl,1,jctr) !remove sea quark from remnant
1453 else
1454 call idsufl(ifl2,2,jct,iret1)
1455 if(iret1.ne.1)then
1456 ifl=idrafl(icltar,jct,1,'s',1,iret2)
1457 endif
1458 endif
1459 else
1460 ifl=ifl2
1461 endif
1462 else
1463 call utstop('No quark for hard Pomeron- in psahot!&')
1464 endif
1465
1466 if(iret1.eq.1.or.iret2.eq.1)then
1467 ifl=ifl2
1468 if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 3)'
1469 if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 3)'
1470 call iddeco(ict,jct)
1471 iret1=0
1472 iret2=0
1473 endif
1474
1475 if(ish.ge.5)write(ifch,*)'flavor vq-,sqb-',iqc(2)
1476 & ,-isign(ifl,iqc(2))
1477
1478 if(ifl.eq.3)then
1479 iqj(nj)=-isign(ifl,iqc(2))*4/3
1480 elseif(ifl.eq.4)then
1481 iqj(nj)=-isign(ifl,iqc(2))/4
1482 else
1483 iqj(nj)=-isign(ifl,iqc(2))
1484 endif
1485
1486 eqj(1,nj)=.5*sngl(wm1+dble(amqt(4))**2/wm1)
1487 eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1488 eqj(3,nj)=xxm1pr(ncolp,kcol)
1489 eqj(4,nj)=xym1pr(ncolp,kcol)
1490 if(ish.ge.8)write (ifch,*)'q_v- mass check',(eqj(1,nj)
1491 * +eqj(2,nj))*(eqj(1,nj)-eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1492 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1493 ncc(1,2)=nj
1494 ncc(2,2)=0
1495 else
1496 nj=nj+1
1497 if(mod(idfpomr,2).ne.0.and.iremn.ne.0)then
1498 if(iremn.eq.3)then
1499 iflq=idrafl(icltar,jctr,1,'v',1,iret3) !take sea quark
1500 iflqb=idrafl(icltar,jctr,2,'v',1,iret3) !take sea antiquark
1501 elseif(iremn.eq.2)then
1502 iflq=idrafl(icltar,jctr,1,'s',0,iret3) !take sea quark
1503 call idsufl3(iflq,1,jctr) !remove sea antiquark from remnant
1504 iflqb=idrafl(icltar,jctr,2,'s',0,iret3) !take sea antiquark
1505 call idsufl3(iflqb,2,jctr) !remove sea antiquark from remnant
1506 else
1507 iflq=idrafl(icltar,jct,1,'s',1,iret2) !Choose flavor of sea quark.
1508 if(iret1.ne.1)iflqb=idrafl(icltar,jct,2,'s',1,iret2) !Choose antiquark.
1509 endif
1510 else
1511 iflq=idrafl(iclpro,jcp,0,'s',0,iret1) !Choose flavor of sea quark.
1512 iflqb=iflq !antiquark=quark (vertex end)
1513 endif
1514
1515 if(iret1.eq.1.or.iret2.eq.1)then
1516 iflqb=iflq
1517 if(ish.ge.3)write(ifmt,*)'Not enough space in rem (psahot 4)'
1518 if(ish.ge.5)write(ifch,*)'Not enough space in rem (psahot 4)'
1519 call iddeco(ict,jct)
1520 iret1=0
1521 iret2=0
1522 endif
1523
1524 if(ish.ge.5)write(ifch,*)'flavor sq-,sqb-',iflq,-iflqb
1525
1526 if(iflqb.eq.3)then
1527 iqj(nj)=-iflqb*4/3
1528 elseif(iflqb.eq.4)then
1529 iqj(nj)=-iflqb/4
1530 else
1531 iqj(nj)=-iflqb
1532 endif
1533 if(iflq.eq.3)then
1534 iqj(nj+1)=iflq*4/3
1535 elseif(iflq.eq.4)then
1536 iqj(nj+1)=iflq/4
1537 else
1538 iqj(nj+1)=iflq
1539 endif
1540 ioj(nj+1)=8
1541
1542 eqj(1,nj)=.5*sngl(wm1+dble(amqt(3))**2/wm1)
1543 eqj(2,nj)=eqj(1,nj)-sngl(wm1)
1544 eqj(3,nj)=xxm2pr(ncolp,kcol)
1545 eqj(4,nj)=xym2pr(ncolp,kcol)
1546 if(ish.ge.8)write (ifch,*)'q_s1- mass check',(eqj(1,nj)-
1547 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1548 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1549 eqj(1,nj+1)=.5*sngl(wm2+dble(amqt(4))**2/wm2)
1550 eqj(2,nj+1)=eqj(1,nj+1)-wm2
1551 eqj(3,nj+1)=xxm1pr(ncolp,kcol)
1552 eqj(4,nj+1)=xym1pr(ncolp,kcol)
1553 nj=nj+1
1554 if(ish.ge.8)write (ifch,*)'q_s2- mass check',(eqj(1,nj)-
1555 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1556 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
1557
1558
1559 if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.1)then
1560 iqc(2)=0
1561 ncc(1,2)=nj-1
1562 ncc(2,2)=nj
1563
1564 else
1565
1566 iqc(2)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
1567 if(1.-1./(wm0*(wpi1-smin/wmi1)).le.xm)goto 1
1568 18 zg=1.d0-dble(rangen())*(1.d0-xm)
1569 if(ish.ge.8)write (ifch,*)'7-zg,xm',zg,xm
1570 if(rangen().gt.zg**dels*((1.d0-xm/zg)/(1.d0-xm))**betpom)
1571 * goto 18
1572 xg=xm/zg
1573 wmq=wm0*(xg-xm)
1574 wpq=1.d0/wmq
1575 wpi1=wpi1-wpq
1576 if(wmi1*wpi1.le.smin)goto 1
1577
1578
1579 nj=nj+1
1580 if(iabs(iqc(2)).eq.3)then
1581 iqj(nj)=-iqc(2)*4/3
1582 else
1583 iqj(nj)=-iqc(2)
1584 endif
1585 ioj(nj)=-8
1586
1587 eqj(1,nj)=.5*sngl(wpq)
1588 eqj(2,nj)=eqj(1,nj)
1589 eqj(3,nj)=0.
1590 eqj(4,nj)=0.
1591 if(ish.ge.8)write (ifch,*)'q_s3- mass check',(eqj(1,nj)-
1592 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
1593 if(iqc(2).gt.0)then
1594 ncj(1,nj)=nj-1
1595 ncj(1,nj-1)=nj
1596 ncj(2,nj)=0
1597 ncj(2,nj-1)=0
1598 ncc(1,2)=nj-2
1599 ncc(2,2)=0
1600 else
1601 ncj(1,nj)=nj-2
1602 ncj(1,nj-2)=nj
1603 ncj(2,nj)=0
1604 ncj(2,nj-2)=0
1605 ncc(1,2)=nj-1
1606 ncc(2,2)=0
1607 endif
1608 endif
1609 endif
1610
1611 endif !iLHC
1612
1613 if(jqq.ne.0)then
1614 if(iqq.ne.0.or.iqq.eq.0.and.jqq.eq.3)then
1615 if(iclpro.ne.4.or.iqq.ne.1)then
1616 call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,1,1)
1617 call psjti0(sngl(wpi1*wmi1-pth),sqaq1,sqaqb1,-1,1)
1618 call psjti0(sngl(wpi1*wmi1-pth),sqqp1,sqqpb1,1,2)
1619 sqq1=(sqq1+sqaq1+2.*(naflav-1)*sqqp1)/naflav/2.
1620 else
1621 call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,4,1)
1622 endif
1623 gbs=sqq1/sqq
1624 else
1625 call psjti0(sngl(wpi1*wmi1-pth),sgq1,sgqb1,0,1)
1626 gbs=sgq1/sgq
1627 endif
1628 if(ish.ge.8)write (ifch,*)'gbs',gbs
1629 if(rangen().gt.gbs)goto 6
1630 endif
1631
1632
1633
1634
1635
1636
1637 if(ish.ge.4)write(ifch,*)
1638 & 'inner partons of the hard Pomeron'
1639 nj00=nj
1640 wpi=wpi1
1641 wmi=wmi1
1642 si=wpi*wmi-pxh**2-pyh**2 !total energy squared for the hard
1643 if(ish.ge.7)write (ifch,*)'si,wpi,wmi',si,wpi,wmi
1644
1645 ept(1)=.5d0*(wpi+wmi)
1646 ept(2)=.5d0*(wpi-wmi)
1647 ept(3)=pxh
1648 ept(4)=pyh
1649 qmin(1)=q2min !effective momentum cutoff above current la
1650 qmin(2)=q2min !effective momentum cutoff below current la
1651 qminn=max(qmin(1),qmin(2)) !effective momentum cutoff for the bo
1652
1653 jfirst=1
1654 jj=int(1.5+rangen())
1655 nemis(1)=0
1656 nemis(2)=0
1657
1658 9 continue ! <<<<----------- ladder rung ---------------------------
1659
1660 if(ish.ge.4)write(ifch,*)'ladder rung'
1661 pt2=ept(3)**2+ept(4)**2
1662 pt=sqrt(sngl(pt2))
1663 if(iabs(iqc(1)).ne.4)then
1664 q2mass=0.
1665 else
1666 q2mass=qcmass**2
1667 si=si-dble(q2mass)
1668 endif
1669 s2min=4.*qminn+q2mass !mass cutoff for born scattering
1670 wwmin=5.*qminn+q2mass-2.*pt*sqrt(q2ini)
1671 wwmin=(wwmin+sqrt(wwmin**2+4.*(q2mass+pt2)*(qminn-q2ini)))
1672 */(1.-q2ini/qminn)/2.
1673 if(ish.ge.5)write(ifch,*)'qminn,q2mass,pt,wwmin,iqc(1):'
1674 if(ish.ge.5)write(ifch,*)qminn,q2mass,pt,wwmin
1675
1676 wpt(1)=ept(1)+ept(2) !lc+ for the current jet emissi
1677 wpt(2)=ept(1)-ept(2) !lc- for the current jet emissi
1678
1679 sjord=0.
1680 if(iabs(iqc(1)).ne.4)then
1681 if(jfirst.eq.1)then
1682 sj=psjti(qmin(jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
1683 sj2=psjti1(qmin(3-jj),qmin(jj),qqcut
1684 * ,sngl(si),iqc(3-jj),iqc(jj),0)
1685 if(ish.ge.5)write(ifch,*)'si,sj,sj2:',si,sj,sj2
1686 if(rangen().gt.sj2/sj.and.si.gt.1.1d0*dble(wwmin))goto 112
1687 jfirst=0
1688 jj=3-jj
1689 sj=sj2
1690 goto 111
1691 endif
1692 sj=psjti1(qmin(jj),qmin(3-jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
1693 111 sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si),iqc(1),iqc(2),0) !born parton-parton
1694 else
1695 sjord=psjci(qmin(2),sngl(si),iqc(2))
1696 sj=sjord
1697 sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si)+q2mass,iqc(1),iqc(2),0)
1698 if(qmin(2).eq.q2min)then
1699 wwmins=2.5*q2min*(1.+sqrt(1.+
1700 * 4.*(pt2+q2mass)/q2min))
1701 if(si.gt.dble(wwmins))call psjti0(sngl(si)+q2mass,sj,sjb
1702 * ,iqc(1),iqc(2))
1703 endif
1704 endif
1705 if(ish.ge.5)write(ifch,*)'si,pt2,wwmin,s2min,sj,sjb,iqc:'
1706 if(ish.ge.5)write(ifch,*)si,pt2,wwmin,s2min,sj,sjb,iqc
1707
1708 if(si.lt.1.1d0*dble(wwmin))goto 12 !------>>>>>>>
1709 if(rangen().lt.sjb/sj)goto 12 !------>>>>>>>
1710
1711 if(iabs(iqc(1)).eq.4)jj=min(2,int(sjord/sj+rangen())+1) !?????????
1712
1713 112 continue
1714
1715 if(iabs(iqc(jj)).ne.4)then
1716
1717 discr=dble((sngl(si)+2.*pt*sqrt(q2ini)-q2mass)**2
1718 * -4.*q2ini*(5.*sngl(si)+q2mass+sngl(pt2)))
1719 if(discr.lt.0.d0.and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
1720 * discr,si,pt,wwmin
1721 discr=dsqrt(discr)
1722 qqmax=(si+2.d0*dble(pt*sqrt(q2ini))-dble(q2mass)+discr)/2.d0
1723 * /(5.d0+(dble(q2mass)+pt2)/si)
1724 qqmin=qqmax-discr/(5.d0+(dble(q2mass)+pt2)/si)
1725 if(s2min.gt.4.d0*qqmin+dble(q2mass))then
1726 xmm=.5d0*(si-s2min+2.d0*dble(pt*sqrt(q2ini)))
1727 discr=xmm**2-si*dble(q2ini)*(1.d0+(dble(q2mass)+pt2)/si)
1728 if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
1729 * discr,si,pt,wwmin
1730 qqmin=(xmm-dsqrt(discr))/(1.d0+(dble(q2mass)+pt2)/si)
1731 endif
1732
1733 xmax=min(1.d0-dble(q2ini)/qqmax,.9999d0)
1734 if(qqmin.lt.dble(qmin(jj)))then
1735 qqmin=dble(qmin(jj))
1736 xmin=max(1.d0-((dble(pt)*dsqrt(qqmin)+dsqrt(pt2*qqmin+
1737 * si*(si-s2min-qqmin*(1.d0+(dble(q2mass)+pt2)/si))))/si)**2,
1738 * (s2min+qqmin*(1.d0+(dble(q2mass)+pt2)/si)-2.d0*dble(pt)
1739 * *dsqrt(qqmin))/si)
1740 if(xmin.le.0.d0)xmin=s2min/si
1741 else
1742 xmin=1.d0-dble(q2ini)/qqmin
1743 endif
1744
1745 qm0=qmin(jj)
1746 xm0=1.0-dble(q2ini/qm0)
1747 if(xm0.gt..999*xmax.or.xm0.lt.1.001*xmin)then
1748 xm0=.5d0*(xmax+xmin)
1749 endif
1750 s2max=sngl(xm0*si-qm0*(dble(q2mass)+pt2)/si)
1751 * +2.*pt*sqrt(q2ini)
1752 xx=xm0
1753
1754 if(iabs(iqc(1)).ne.4)then
1755 if(jfirst.eq.1)then
1756 sj0=psjti(qm0,qqcut,s2max,0,iqc(3-jj),0)*
1757 * psfap(xx,iqc(jj),0)+
1758 * psjti(qm0,qqcut,s2max,7,iqc(3-jj),0)
1759 * *psfap(xx,iqc(jj),1)
1760 else
1761 sj0=psjti1(qm0,qmin(3-jj),qqcut,s2max,0,iqc(3-jj),0)*
1762 * psfap(xx,iqc(jj),0)+
1763 * psjti1(qm0,qmin(3-jj),qqcut,s2max,7,iqc(3-jj),0)
1764 * *psfap(xx,iqc(jj),1)
1765 endif
1766 else
1767 sj0=psjci(qm0,s2max,0)*psfap(xx,iqc(jj),0)+
1768 * psjci(qm0,s2max,1)*psfap(xx,iqc(jj),1)
1769 endif
1770
1771 gb0=dble(sj0/log(q2ini/qcdlam)*qm0*5.)*psuds(qm0,iqc(jj))
1772 if(ish.ge.5)write(ifch,*)'gb0,qm0,xm0,s2max:',
1773 * gb0,qm0,xm0,s2max
1774
1775 if(gb0.le.0.d0)then
1776 write(ifmt,*)'gb0.le.0. si,pt2:',si,pt2
1777 iret=1
1778 goto 16
1779 endif
1780
1781 if(xm0.le..5d0)then
1782 gb0=gb0*xm0**(1.-delh)
1783 else
1784 gb0=gb0*(1.d0-xm0)*2.d0**delh
1785 endif
1786
1787 xmin2=max(.5d0,xmin)
1788 xmin1=xmin**delh
1789 xmax1=min(xmax,.5d0)**delh
1790 if(xmin.ge..5d0)then
1791 djl=1.
1792 elseif(xmax.le..5d0)then
1793 djl=0.
1794 else
1795 djl=1./(1.+((2.*sngl(xmin))**delh-1.)/delh/
1796 * log(2.*(1.-sngl(xmax))))
1797 endif
1798
1799 else
1800
1801 xmin=5.d0*dble(q2min)/si
1802 xmax=min(si/(si+5.0*(pt2+dble(q2mass))),.9999d0)
1803 qqmax=xmax*si/5.d0
1804 qqmin=dble(q2min)
1805
1806 qm0=sngl(qqmin)
1807 xm0=2.d0/(1.d0+sqrt(1.d0+4.d0*(pt2+dble(q2mass))/qm0))
1808 if(xm0.gt..999d0*xmax.or.xm0.lt.1.001d0*xmin)then
1809 xm0=.5d0*(xmax+xmin)
1810 endif
1811 s2max=sngl(xm0*si)
1812 xx=xm0
1813
1814 sj0=psjti(qm0,qmin(3-jj),s2max,0,iqc(3-jj),0)*
1815 * psfap(xx,iqc(jj),0)
1816 gb0=dble(sj0/log(qm0/qcdlam)*qm0 *5.)
1817 gb0=gb0*xm0**(1.-delh)
1818 if(gb0.le.0.d0)then
1819 if(ish.ge.2)write(ifch,*)'gb0.le.0. (charm) si,pt2:',si,pt2
1820 iret=1
1821 goto 16
1822 endif
1823 djl=0.
1824 xmin2=0d0
1825 xmin1=xmin**delh
1826 xmax1=xmax**delh
1827
1828 endif
1829
1830 if(ish.ge.5)write(ifch,*)'xmin,xmax,qqmin,qqmax:',
1831 *xmin,xmax,qqmin,qqmax
1832
1833 ntry=0
1834 10 continue
1835 ntry=ntry+1
1836 if(ntry.gt.5000000)then
1837 if(ish.ge.1)write(*,*)'Reject hard string (too many rejection)'
1838 &,kcol,ncolp,nptl,gb7
1839 iret=1
1840 goto 16
1841 endif
1842 if(rangen().gt.djl)then !lc momentum share in the cur
1843 x=(xmin1+dble(rangen())*(xmax1-xmin1))**(1./delh)
1844 else
1845 x=1.d0-(1.d0-xmin2)*((1.d0-xmax)/(1.d0-xmin2))**rangen()
1846 endif
1847 qq=qqmin/(1.d0+dble(rangen())*(qqmin/qqmax-1.d0))
1848 if(ish.ge.7)write(ifch,*)'x,qq:',x,qq,ntry
1849
1850 if(iabs(iqc(jj)).ne.4)then
1851
1852 qt2=qq*(1.d0-x)
1853 if(qt2.lt.dble(q2ini))then
1854 if(ish.ge.7)write(ifch,*)'qt2:',qt2
1855 goto 10
1856 endif
1857
1858 qmin2=max(qminn,sngl(qq))
1859 qt=dsqrt(qt2)
1860 call pscs(bcos,bsin)
1861 ep3(3)=sngl(qt)*bcos !new parton pt-s
1862 ep3(4)=sngl(qt)*bsin
1863 ptnew=(ept(3)-dble(ep3(3)))**2+(ept(4)-dble(ep3(4)))**2
1864 s2min2=4.*qmin2+q2mass
1865
1866 s2=x*(si-qq)-ptnew-qq*(dble(q2mass)+pt2)/si+pt2 !new ladder mass
1867 if(s2.lt.dble(s2min2))then
1868 if(ish.ge.7)write(ifch,*)'s2,s2min2:',s2,s2min2
1869 goto 10 !rejection in case of too low mass
1870 endif
1871
1872 xx=x
1873 if(iabs(iqc(1)).ne.4)then
1874 if(jfirst.eq.1)then
1875 sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1876 if(iqc(jj).ne.0)then !f1 - f2
1877 sj2=psjti(sngl(qq),qqcut,sngl(s2),iqc(jj),iqc(3-jj),0)
1878 elseif(iqc(3-jj).eq.0)then !f1 - g
1879 sj2=psjti(sngl(qq),qqcut,sngl(s2),1,0,0)
1880 else !g - f2
1881 sj2=psjti(sngl(qq),qqcut,sngl(s2),1,1,0)/naflav/2. !q - q
1882 * +psjti(sngl(qq),qqcut,sngl(s2),-1,1,0)/naflav/2. !q~ - q
1883 * +psjti(sngl(qq),qqcut,sngl(s2),1,2,0)*(1.-1./naflav) !q' - q
1884 endif
1885 else
1886 sj1=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),0,iqc(3-jj),0)
1887 if(iqc(jj).ne.0)then !f1 - f2
1888 sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1889 * ,sngl(s2),iqc(jj),iqc(3-jj),0)
1890 elseif(iqc(3-jj).eq.0)then !f1 - g
1891 sj2=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),1,0,0)
1892 else !g - f2
1893 sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1894 * ,sngl(s2),1,1,0)/naflav/2. !q - q
1895 * +psjti1(sngl(qq),qmin(3-jj),qqcut
1896 * ,sngl(s2),-1,1,0)/naflav/2. !q~ - q
1897 * +psjti1(sngl(qq),qmin(3-jj),qqcut
1898 * ,sngl(s2),1,2,0)*(1.-1./naflav)!q' - q
1899 endif
1900 endif
1901 else
1902 sj1=psjci(sngl(qq),sngl(s2),0)
1903 sj2=psjci(sngl(qq),sngl(s2),1)
1904 endif
1905
1906 !...gb7 - rejection function for x and q**2 simulation
1907 gb7=dble((sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))/
1908 * log(sngl(qt2)/qcdlam))*psuds(sngl(qq),iqc(jj))*qq/gb0
1909
1910 if(x.le..5d0)then
1911 gb7=gb7*x**(1.-delh)
1912 else
1913 gb7=gb7*(1.d0-x)*2.d0**delh
1914 endif
1915 if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1916 * gb7,qq,x,qm0,xm0
1917 if(ish.ge.7)write(ifch,*)'gb7:',gb7
1918 if(dble(rangen()).gt.gb7)goto 10
1919 else
1920 qmin2=max(qminn,sngl(qq))
1921 s2min2=4.*qmin2
1922 s2=x*si-qq !new ladder mass
1923 if(s2.lt.dble(s2min2))goto 10 !rejection in case of too low mass
1924
1925 call pscs(bcos,bsin)
1926 xmm=x*(ept(3)*dble(bcos)+ept(4)*dble(bsin))
1927 discr=xmm**2+qq*(1.d0-x)-x**2*(pt2+dble(q2mass))
1928 if(discr.lt.0.d0)goto 10
1929 qt=xmm+dsqrt(discr)
1930 ep3(3)=sngl(ept(3)-qt*dble(bcos)) !new parton pt-s
1931 ep3(4)=sngl(ept(4)-qt*dble(bsin))
1932 qt2=dble(ep3(3)**2+ep3(4)**2)
1933
1934 xx=x
1935 sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1936 sj2=0.
1937 !.... gb7 - rejection function for x and q**2 simulation
1938 gb7=dble(sj1*psfap(xx,iqc(jj),0)/
1939 * log(sngl(qq)/qcdlam))*qq/gb0
1940 gb7=gb7*x**(1.-delh)
1941 if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1942 * gb7,qq,x,qm0,xm0
1943 if(ish.ge.7)write(ifch,*)'gb7:',gb7
1944 if(dble(rangen()).gt.gb7)goto 10
1945 endif
1946
1947 ! parton type selection iqc -> iqnew (space like) + iq1 (time like)
1948 nqc(2)=0
1949 iqnew=iqc(jj)
1950 if(rangen().lt.sj1/(sj1+sj2))then
1951 if(iqc(jj).eq.0)then !g -> g + g (jt=1)
1952 jt=1
1953 jq=int(1.5+rangen())
1954 nqc(1)=ncc(jq,jj)
1955 else !q -> g + q (jt=2)
1956 jt=2
1957 if(iqc(jj).gt.0)then
1958 jq=1
1959 else
1960 jq=2
1961 endif
1962 nqc(1)=0
1963 iqnew=0
1964 endif
1965 iq1=iqc(jj)
1966 jo=ncc(jq,jj) !parton origin
1967 if(jo.ne.0)then
1968 jo=ioj(jo)
1969 else
1970 jo=ioj(ncc(3-jq,jj))
1971 endif
1972 else
1973 if(iqc(jj).ne.0)then !q -> q + g (jt=3)
1974 iq1=0
1975 jt=3
1976 if(iqc(jj).gt.0)then
1977 jq=1
1978 else
1979 jq=2
1980 endif
1981 nqc(1)=ncc(1,jj)
1982 else !g -> q + q (jt=4)
1983 jt=4
1984 jq=int(1.5+rangen())
1985 iq1=int(naflav*rangen()+1.)*(3-2*jq)
1986 iqnew=-iq1
1987 nqc(1)=ncc(jq,jj)
1988 endif
1989 jo=ncc(1,jj) !parton origin
1990 if(jo.ne.0)then
1991 jo=ioj(jo)
1992 else
1993 jo=ioj(ncc(2,jj))
1994 endif
1995 endif
1996 if(jo.ge.70.and.jo.ne.99)jo=jo/10
1997 if(ish.ge.5)write(ifch,'(a,i3,a4,i3,a4,i3)')' Process :'
1998 * ,iqc(jj),' -> ',iqnew,' + ',iq1
1999
2000 if(iabs(iqc(jj)).ne.4)then
2001 q2part=0.
2002
2003 if(dabs(wpt(3-jj)).gt.1.d-20)then
2004 wplc=wpt(jj)-(pt2+dble(q2mass))/wpt(3-jj)
2005 else
2006 if(ish.gt.1)write(ifmt,*)'Problem with wpt in sem',wpt(3-jj)
2007 wplc=wpt(jj)-(pt2+dble(q2mass))
2008 endif
2009 qp2max=max(0.,sngl(qt2))
2010 else
2011 q2part=q2mass
2012 wplc=wpt(jj)
2013 qp2max=max(0.,sngl(qq))
2014 endif
2015 eprt=max(dsqrt(qt2+dble(q2part)),.5d0*((1.d0-x)*wplc+
2016 *(qt2+dble(q2part))/(1.d0-x)/wplc))
2017 pl=((1.d0-x)*wplc-eprt)*dble(3-2*jj)
2018 zeta=sqrt(qp2max/si)/sqrt(x*(1.-x))
2019 if(iq1.eq.0)then
2020 iq2ini=9
2021 jo=iq2ini
2022 if(zeta.gt.zetacut)jo=-jo
2023 else
2024 iq2ini=iq1
2025 jo=iq2ini
2026 endif
2027 if(ish.ge.5)write(ifch,*)'qq,eprt,iq2ini,jo,E2-Q2',qp2max,eprt
2028 *,iq2ini,jo,eprt**2-qt2
2029 ntest=0
2030 11 ntest=ntest+1
2031 call timsh1(qp2max,sngl(eprt),iq2ini,jo)
2032 amprt=pprt(5,1)**2
2033 plprt=max(0.d0,eprt**2-qt2)-dble(amprt)
2034 if(plprt.lt.0.d0)then
2035 if(ntest.lt.10000)then
2036 goto 11
2037 else
2038 iret=1
2039 goto 16
2040 endif
2041 endif
2042 ep3(1)=sngl(eprt)
2043 ep3(2)=sngl(dsqrt(plprt))
2044 if(pl.lt.0.d0)ep3(2)=-ep3(2)
2045 ey(1)=1.
2046 ey(2)=1.
2047 ey(3)=1.
2048 do i=1,4
2049 ept1(i)=ept(i)-dble(ep3(i))
2050 enddo
2051 call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
2052 s2new=psnorm(ept1)
2053
2054 if(iabs(iqc(jj)).ne.4.and.s2new-q2mass.gt.s2min2.or.
2055 *iabs(iqc(jj)).eq.4.and.s2new.gt.s2min2)then
2056 if(iabs(iqc(1)).ne.4.or.jj.eq.1)then
2057 if(jfirst.eq.1)then
2058 gb=dble(psjti(sngl(qq),qqcut,s2new,iqnew,iqc(3-jj),0))
2059 else
2060 gb=dble(psjti1(sngl(qq),qmin(3-jj),qqcut,s2new,iqnew,
2061 * iqc(3-jj),0))
2062 endif
2063 else
2064 gb=dble(psjci(sngl(qq),s2new-q2mass,iqnew))
2065 endif
2066 if(iqnew.eq.0)then
2067 gb=gb/dble(sj1)
2068 else
2069 gb=gb/dble(sj2)
2070 endif
2071 if(dble(rangen()).gt.gb)goto 10
2072 else
2073 goto 10
2074 endif
2075
2076 if(ish.ge.6)write(ifch,*)'jt,jj,jq,nqc:',jt,jj,jq,nqc
2077 nprtjx=0
2078 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
2079
2080 if(jt.eq.1)then
2081 ncc(jq,jj)=nqc(2)
2082 elseif(jt.eq.2)then
2083 ncc(jq,jj)=ncc(1,jj)
2084 ncc(3-jq,jj)=nqc(1)
2085 elseif(jt.eq.3)then
2086 ncc(1,jj)=nqc(2)
2087 elseif(jt.eq.4)then
2088 ncc(1,jj)=ncc(3-jq,jj)
2089 endif
2090 iqc(jj)=iqnew
2091
2092 do i=1,4
2093 ept(i)=ept1(i)
2094 enddo
2095 ! c.m. energy squared, minimal 4-momentum transfer square and gluon 4-v
2096 ! for the next ladder run
2097 qmin(jj)=sngl(qq)
2098 qminn=qmin2
2099 si=dble(s2new)
2100 nemis(jj)=nemis(jj)+1
2101
2102 goto 9 ! ---------------next ladder rung ------------>>>>>>>>>>>
2103
2104
2105 12 continue !------------------- Born process------------------------
2106
2107 if(ish.ge.4)write(ifch,*)'Born process'
2108 if(ish.ge.6)write(ifch,*)'iqc,si:',iqc,si
2109 qq=dble(qminn)
2110
2111 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2112 xpprbor(ncolp,kcol)=(ept(1)+ept(2))/dble(engy)
2113 xmprbor(ncolp,kcol)=(ept(1)-ept(2))/dble(engy)
2114 nemispr(1,ncolp,kcol)=nemis(1)
2115 nemispr(2,ncolp,kcol)=nemis(2)
2116 ! write(*,'(a,2f8.3,i3,3x,2i3)')'------------'
2117 ! * ,xpprbor(ncolp,kcol),xmprbor(ncolp,kcol),nj-nj00,nemis
2118 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2119
2120 call psabor(sngl(si),sngl(qq),iqc,ncc,ept,0,iptl,bx)
2121
2122 !kkkkkkkkkkkkk out Born partons without timelike cascade
2123 if(nprtjx.eq.2)then
2124 do ii=1,2
2125 ptprboo(ii,ncolp,kcol)=sqrt(pprtx(1,ii)**2+pprtx(2,ii)**2)
2126 enddo
2127 else
2128 stop'psahot: should not happen!!!! '
2129 endif
2130 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2131 ! ptxxx=max( ptprboo(1,ncolp,kcol) , ptprboo(2,ncolp,kcol) )
2132 ! print*,sqrt(wm0*wp0),sqrt(wmi*wpi),ptxxx !++++++++++++++++++++++++++
2133 ! if(ptxxx.lt.10)goto1
2134 ! print*,' ++++++++++++++++++++++++++++++++++++ '
2135 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
2136 if(nj.ne.0.)then
2137 do i=1,nj
2138 do l=1,6
2139 bxj(l,i)=bx(l)
2140 enddo
2141 ityj(i)=ity
2142 iorj(i)=iptl
2143 enddo
2144 endif
2145
2146 call psjarr(jfl) !kinky strings formation
2147
2148 if(jfl.eq.0.and.ish.ge.4)write(ifch,*)
2149 *'jfl,nj,nptl',jfl,nj,nptl
2150 if(jfl.eq.0)goto 1
2151
2152
2153
2154 iret1=0
2155 call idenco(jcp,icp,iret1)
2156 if(iret1.eq.1)call utstop('Problem with proj rem in psahot !&')
2157 iret2=0
2158 call idenco(jct,ict,iret2)
2159 if(iret2.eq.1)call utstop('Problem with targ rem in psahot !&')
2160
2161 do i=1,2
2162 icproj(i,ip)=icp(i)
2163 ictarg(i,it)=ict(i)
2164 enddo
2165
2166 if(ish.ge.3)write(ifch,*)'End psahot (icp,ict):',icp,ict
2167
2168 if(iremn.ge.2)then !uses precalculated flavors
2169
2170 do j=1,2
2171 do n=1,nrflav
2172 jcpref(n,j,ip)=jcpr(n,j)
2173 jctref(n,j,it)=jctr(n,j)
2174 enddo
2175 enddo
2176 if(ish.ge.3)then
2177 write(ifch,'(a,6i3,2x,6i3)')' proj: ',jcpr
2178 write(ifch,'(a,6i3,2x,6i3)')' targ: ',jctr
2179 endif
2180
2181 endif
2182
2183
2184
2185
2186 16 continue
2187 call utprix('psahot',ish,ishini,3)
2188
2189 q2fin=q2finsave
2190 return
2191 end
2192
2193
2194 subroutine psjarr(jfl)
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204 parameter (mjstr=20000)
2205 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2206
2207
2208
2209
2210
2211
2212
2213 dimension mark(mjstr)
2214 double precision ept(4),eptot(2)
2215 include 'epos.inc'
2216 include 'epos.incsem'
2217
2218 if(nj.eq.0)then
2219 jfl=1
2220 return
2221 endif
2222
2223
2224 do k=1,nj
2225 if(iqj(k).eq.0)then !gluon must have two neighbours
2226 if(ncj(1,k).eq.0)then !first neigbour missing
2227
2228 do kk=1,nj !look to which parton he is connected
2229 if(ncj(1,kk).eq.k)then
2230 if(ncj(2,k).ne.kk)ncj(1,k)=kk !if not already connected : connection
2231 elseif(ncj(2,kk).eq.k)then
2232 if(ncj(1,k).ne.kk)ncj(1,k)=kk
2233 endif
2234 enddo
2235 endif
2236 if(ncj(2,k).eq.0)then !second neigbour missing
2237
2238 do kk=1,nj
2239 if(ncj(2,kk).eq.k)then
2240 if(ncj(1,k).ne.kk)ncj(2,k)=kk
2241 elseif(ncj(1,kk).eq.k)then
2242 if(ncj(2,k).ne.kk)ncj(2,k)=kk
2243 endif
2244 enddo
2245 endif
2246 endif
2247 enddo
2248
2249
2250
2251
2252
2253
2254 if(ish.ge.3)then
2255 write (ifch,*)'psjarr: nj',nj
2256 do k=1,nj
2257 if(iqj(k).ne.0)ncj(2,k)=0
2258 write(ifch,'(a,i4)')' parton',k
2259 write(ifch,'(i6,2x,4e10.3,2x,2i3)')iqj(k)
2260 * ,(eqj(j,k),j=1,4),(ncj(j,k),j=1,2)
2261 enddo
2262 endif
2263
2264 jfl=0
2265 do i=1,nj
2266 mark(i)=1
2267 enddo
2268 nptl0=nptl
2269
2270 eptot(1)=0d0
2271 eptot(2)=0d0
2272
2273 1 continue
2274 do ij=1,nj
2275 if(mark(ij).ne.0.and.iqj(ij).ne.0)goto 2
2276 enddo
2277 2 continue
2278 jfirst=1
2279
2280
2281 ij0=ncj(1,ij)
2282 kkk=1
2283 if(ij0.gt.0)then
2284 if(ncj(1,ij0).eq.ij)kkk=1
2285 if(ncj(2,ij0).eq.ij)kkk=2
2286 endif
2287 if(iabs(iqj(ij)).le.2)then
2288 am1=amhadr(1)
2289 elseif(iabs(iqj(ij)).eq.4)then
2290 am1=amhadr(3)
2291 elseif(iabs(iqj(ij)).eq.40)then
2292 am1=qcmass
2293 else
2294 am1=amhadr(2)
2295 endif
2296
2297 do i=1,4
2298 ept(i)=0.
2299 enddo
2300
2301 3 continue
2302 call pspawr(ij,kkk)
2303
2304 mark(ij)=0
2305
2306 do i=1,4
2307 ept(i)=ept(i)+eqj(i,ij)
2308 enddo
2309 eptot(kkk)=eptot(kkk)+eqj(1,ij)
2310
2311 if(iqj(ij).ne.0)then
2312 if(jfirst.ne.1)then
2313 if(iabs(iqj(ij)).le.2)then
2314 am2=amhadr(1)
2315 elseif(iabs(iqj(ij)).eq.4)then
2316 am2=amhadr(3)
2317 elseif(iabs(iqj(ij)).eq.40)then
2318 am2=qcmass
2319 else
2320 am2=amhadr(2)
2321 endif
2322 amj=(am1+am2+stmass)**2
2323 sm=psnorm(ept)
2324 if(sm.lt.amj)then
2325 nptl=nptl0
2326 goto 999
2327 endif
2328
2329 if(nptl-nptl0.lt.nj)then
2330 goto 1
2331 else
2332 if(iLHC.ne.1)then
2333
2334 do k=nptl0+1,nptl
2335 qsqptl(k)=sngl(eptot(nint(qsqptl(k)))**2)
2336 enddo
2337 endif
2338
2339 if(ish.ge.3)then
2340 write (ifch,*)'psjarr: nptl',nptl
2341 do k=nptl0+1,nptl
2342 write(ifch,'(a,i4)')' particle',k
2343 write(ifch,'(i5,2x,6e10.3)')idptl(k)
2344 * ,(pptl(j,k),j=1,5),sqrt(qsqptl(k))
2345 enddo
2346 endif
2347 jfl=1
2348 goto 999
2349 endif
2350 else
2351 jfirst=0
2352 njpar=ij
2353 ij=ncj(1,ij)
2354 if(ij.eq.0)write(ifch,*)
2355 &'IN PSJARR ij=0 :parton,',njpar,'with no connection'
2356 goto 3
2357 endif
2358
2359 else
2360 if(ncj(1,ij).eq.njpar)then
2361 njdau=ncj(2,ij)
2362 else
2363 njdau=ncj(1,ij)
2364 endif
2365 njpar=ij
2366 ij=njdau
2367 if(ij.eq.0)write(ifch,*)
2368 &'IN PSJARR ij=0 :parton,',njpar, 'with no connection'
2369 goto 3
2370 endif
2371
2372
2373 999 continue
2374
2375 end
2376
2377
2378 subroutine pspawr(kj,kkk)
2379
2380
2381
2382
2383 parameter (mjstr=20000)
2384 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2385 common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
2386
2387
2388
2389
2390
2391
2392
2393
2394 include 'epos.inc'
2395 include 'epos.incsem'
2396 common/ciptl/iptl
2397
2398 nptl=nptl+1
2399 if(ish.ge.9)write (ifch,*)'nptl,kj (sto)',nptl,kj
2400 if(nptl.ge.mxptl.or.kj.le.0)then
2401 write (ifmt,*)'nptl,kj',nptl,kj
2402 call alist('Error in pspawr: nptl or kj out of bounds &',1,nptl)
2403 call utstop('nptl or kj out of bounds&')
2404 endif
2405
2406 if(ifrptl(1,iptl).eq.0)ifrptl(1,iptl)=nptl
2407 ifrptl(2,iptl)=nptl
2408
2409 pptl(1,nptl)=eqj(3,kj)
2410 pptl(2,nptl)=eqj(4,kj)
2411 pptl(3,nptl)=eqj(2,kj)
2412 pptl(4,nptl)=eqj(1,kj)
2413 pptl(5,nptl)=0.
2414 idptl(nptl)=psidd(iqj(kj))
2415 iorptl(nptl)=iorj(kj)
2416 jorptl(nptl)=ioj(kj)
2417 istptl(nptl)=20
2418 do i=1,4
2419 xorptl(i,nptl)=bxj(i,kj)
2420 enddo
2421 tivptl(1,nptl)=bxj(5,kj)
2422 tivptl(2,nptl)=bxj(6,kj)
2423 ityptl(nptl)=ityj(kj)
2424
2425 if(iLHC.eq.1)then
2426 qsqptl(nptl)=q2j(kj)
2427 else
2428 qsqptl(nptl)=float(kkk)
2429 endif
2430 !kkkkkkkkkkkkkkkkkkkkkkk
2431 ! write(*,'(a,2i4,i6,f8.3)')'.... ',kj,nptl,idptl(nptl)
2432 !* ,sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2)
2433 return
2434 end
2435
2436
2437 subroutine psabor(si,qq,iqc,ncc,ept,jdis,iptl,coordo)
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447 double precision ept(4),psutz,psuds
2448 dimension ep3(4),ey(3),wsub(5),iqc(2),ncc(2,2),nqc(2),coordo(6)
2449 parameter (mjstr=20000)
2450 common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2451 parameter (ntim=1000)
2452 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
2453 &,iorprt(ntim),jorprt(ntim),nprtj
2454 common/cprtx/nprtjx,pprtx(5,2)
2455 include 'epos.incsem'
2456 include 'epos.inc'
2457 call utpri('psabor',ish,ishini,5)
2458
2459 if(iabs(iqc(1)).ne.4)then !gluon or light quark
2460 q2mass=0.
2461 else !c-quark
2462 q2mass=qcmass**2
2463 endif
2464 p1=si/(1.+q2mass/si)
2465 if(p1.gt.4.*qq)then !|t|-cutoff (p^2_t>qq)
2466 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
2467 else
2468 tmin=2.*qq
2469 endif
2470 tmax=.5*p1
2471
2472 fborn=0.
2473 qt2=tmin*(1.d0-tmin/p1)
2474 if(qt2.lt..999d0*qq.and.ish.ge.2)write(ifmt,*)'qt20,qq',qt2,qq
2475 if(iqc(1).ne.0.or.iqc(2).ne.0)then
2476 do l=1,5 !sum over different subprocesses
2477 wsub(l)=psbori(si,tmin,iqc(1),iqc(2),l) !matrix element
2478 if(l.le.3)then
2479 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
2480 elseif(l.le.4)then
2481 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
2482 else
2483 wsub(l)=wsub(l)*(alfe/2/pi)**2
2484 endif
2485 fborn=fborn+wsub(l)
2486 enddo
2487 fborn=tmin**2*fborn
2488 else
2489 do l=1,5
2490 wsub(l)=psbori(si,.5*si,iqc(1),iqc(2),l)
2491 if(l.le.3)then
2492 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
2493 elseif(l.le.4)then
2494 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
2495 else
2496 wsub(l)=wsub(l)*(alfe/2/pi)**2
2497 endif
2498 fborn=fborn+wsub(l)
2499 enddo
2500 fborn=.25*si**2*fborn
2501 endif
2502 if(jdis.eq.0)then
2503 scale1=qt2
2504 else
2505 scale1=4.*qt2
2506 endif
2507 gb0=dble(fborn)*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))
2508 if(ish.ge.7)write(ifch,*)'tmin,gb0:',tmin,gb0
2509
2510
2511
2512
2513 14 q2=tmin/(1.-rangen()*(1.-tmin/tmax)) !q2=min(|t|,|u|)
2514 qt2=q2*(1.-q2/p1) !qt2=p_t^2 for the process
2515 if(qt2.lt.qq.and.ish.ge.2)write(ifmt,*)'qt2,qq',qt2,qq
2516 if(rangen().lt..5)then !|u|=q2, |t|=p1-q2
2517 jm=2 !first parton to be considered
2518 tq=p1-q2
2519 else !|t|=q2, |u|=p1-q2
2520 jm=1 !first parton to be considered
2521 tq=q2
2522 endif
2523
2524 fborn=0.
2525 do l=1,5 !sum over different subprocesses
2526 wsub(l)=psbori(si,tq,iqc(1),iqc(2),l)
2527 if(l.le.3)then
2528 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
2529 elseif(l.le.4)then
2530 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
2531 else
2532 wsub(l)=wsub(l)*(alfe/2/pi)**2
2533 endif
2534 fborn=fborn+wsub(l)
2535 enddo
2536 if(jdis.eq.0)then
2537 scale1=qt2
2538 else
2539 scale1=4.*qt2
2540 endif
2541 gb=dble(q2**2*fborn)
2542 &*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))/gb0 !rejection function
2543 if(ish.ge.7)write(ifch,*)'q2,qt2,gb:',q2,qt2,gb
2544
2545 if(dble(rangen()).gt.gb)goto 14 !rejection
2546
2547
2548 nqc(2)=0
2549 if(iqc(1).eq.0.and.iqc(2).eq.0)then !g+g
2550 jq=int(1.5+rangen()) !jq=1(2) - transfer of color (anticolor)
2551 nqc(1)=ncc(jq,jm)
2552
2553 if(rangen().lt.wsub(1)/fborn)then !gg->gg
2554 if(rangen().lt..5)then
2555 jt=1 !anticolor-color annihilation
2556 nqc(2)=0
2557 njc1=ncc(3-jq,jm)
2558 njc2=ncc(jq,3-jm)
2559 if(iqj(njc1).ne.0)then
2560 ncj(1,njc1)=njc2
2561 else
2562 ncj(jq,njc1)=njc2
2563 endif
2564 if(iqj(njc2).ne.0)then
2565 ncj(1,njc2)=njc1
2566 else
2567 ncj(3-jq,njc2)=njc1
2568 endif
2569 else
2570 jt=2 !produced gluons get color and
2571 nqc(2)=ncc(3-jq,3-jm) !anticolor from the 2 parents
2572 endif
2573 else !gg->qq~
2574 jt=9 !anticolor-color annihilation
2575 iqc(jm)=int(naflav*rangen()+1)*(3-2*jq) !(anti)quark flavor
2576 iqc(3-jm)=-iqc(jm)
2577 njc1=ncc(3-jq,jm)
2578 njc2=ncc(jq,3-jm)
2579 if(iqj(njc1).ne.0)then
2580 ncj(1,njc1)=njc2
2581 else
2582 ncj(jq,njc1)=njc2
2583 endif
2584 if(iqj(njc2).ne.0)then
2585 ncj(1,njc2)=njc1
2586 else
2587 ncj(3-jq,njc2)=njc1
2588 endif
2589 endif
2590
2591 elseif(iqc(1)*iqc(2).eq.0)then !q(q~)+g
2592 if(iqc(1)+iqc(2).gt.0)then
2593 jq=1 !q
2594 else
2595 jq=2 !q~
2596 endif
2597 if(rangen().lt.wsub(1)/fborn)then !q(q~)g->q(q~)g
2598 if(rangen().lt..5)then !anticolor-color annihilation
2599 if(iqc(jm).eq.0)then
2600 jt=3 !first parton=g
2601 nqc(1)=ncc(jq,jm)
2602 njc1=ncc(3-jq,jm)
2603 njc2=ncc(1,3-jm)
2604 if(iqj(njc1).ne.0)then
2605 ncj(1,njc1)=njc2
2606 else
2607 ncj(jq,njc1)=njc2
2608 endif
2609 if(iqj(njc2).ne.0)then
2610 ncj(1,njc2)=njc1
2611 else
2612 ncj(3-jq,njc2)=njc1
2613 endif
2614
2615 else
2616 jt=4 !first parton=q(q~)
2617 nqc(1)=0
2618 njc1=ncc(1,jm)
2619 njc2=ncc(3-jq,3-jm)
2620 if(iqj(njc1).ne.0)then
2621 ncj(1,njc1)=njc2
2622 else
2623 ncj(3-jq,njc1)=njc2
2624 endif
2625 if(iqj(njc2).ne.0)then
2626 ncj(1,njc2)=njc1
2627 else
2628 ncj(jq,njc2)=njc1
2629 endif
2630 endif
2631
2632 else !color transfer
2633 if(iqc(jm).eq.0)then
2634 jt=5 !first parton=g
2635 nqc(2)=ncc(3-jq,jm)
2636 nqc(1)=ncc(1,3-jm)
2637 else !first parton=q(q~)
2638 jt=6
2639 nqc(1)=ncc(jq,3-jm)
2640 endif
2641 endif
2642
2643 else !q(q~)g->q(q~)-gamma (+-color annihilation)
2644 if(iqc(jm).eq.0)then
2645 jt=11 !first parton=g
2646 nqc(1)=ncc(jq,jm)
2647 njc1=ncc(3-jq,jm)
2648 njc2=ncc(1,3-jm)
2649 if(iqj(njc1).ne.0)then
2650 ncj(1,njc1)=njc2
2651 else
2652 ncj(jq,njc1)=njc2
2653 endif
2654 if(iqj(njc2).ne.0)then
2655 ncj(1,njc2)=njc1
2656 else
2657 ncj(3-jq,njc2)=njc1
2658 endif
2659 iqc(jm)=iqc(3-jm)
2660 iqc(3-jm)=10 !make the second output is gamma.
2661
2662 else
2663 jt=12 !first parton=q(q~)
2664 nqc(1)=ncc(jq,3-jm) !here nqc(1) is gluon.
2665 njc1=ncc(1,jm)
2666 njc2=ncc(3-jq,3-jm)
2667 if(iqj(njc1).ne.0)then
2668 ncj(1,njc1)=njc2
2669 else
2670 ncj(3-jq,njc1)=njc2
2671 endif
2672 if(iqj(njc2).ne.0)then
2673 ncj(1,njc2)=njc1
2674 else
2675 ncj(jq,njc2)=njc1
2676 endif
2677 iqc(3-jm)=10
2678 endif
2679 endif
2680
2681 elseif(iqc(1)*iqc(2).gt.0)then
2682 jt=7 !qq->qq (q~q~->q~q~)
2683 if(iqc(1).gt.0)then
2684 jq=1
2685 else
2686 jq=2
2687 endif
2688 nqc(1)=ncc(1,3-jm)
2689
2690 else ! qq~ ->
2691 if(iqc(jm).gt.0)then
2692 jq=1
2693 else
2694 jq=2
2695 endif
2696 aks=rangen()
2697 if(aks.lt.(wsub(1)+wsub(2))/fborn)then
2698 jt=8 ! qq~->qq~ (anticolor-color annihilation)
2699 if(aks.gt.wsub(1)/fborn)then
2700 iqa=iabs(iqc(jm))
2701 iq=int((naflav-1)*rangen())+1
2702 if(iq.eq.iqa)iq=naflav
2703 iqc(jm)=iq*iqc(jm)/iqa
2704 iqc(3-jm)=-iqc(jm)
2705 endif
2706 nqc(1)=0
2707 njc1=ncc(1,jm)
2708 njc2=ncc(1,3-jm)
2709 if(iqj(njc1).ne.0)then
2710 ncj(1,njc1)=njc2
2711 else
2712 ncj(3-jq,njc1)=njc2
2713 endif
2714 if(iqj(njc2).ne.0)then
2715 ncj(1,njc2)=njc1
2716 else
2717 ncj(jq,njc2)=njc1
2718 endif
2719 elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3))/fborn)then
2720 jt=10 !color transfer qq~->gg
2721 iqc(1)=0
2722 iqc(2)=0
2723 nqc(1)=ncc(1,jm)
2724 nqc(2)=0
2725 elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3)+wsub(4))/fborn)then
2726 jt=13 ! qq~->g+gamma
2727 nqc(1)=ncc(1,jm)
2728 nqc(2)=ncc(1,3-jm)
2729 iqc(jm)=0
2730 iqc(3-jm)=10
2731 else
2732 jt=14 ! qq~->gamma+gamma
2733 njc1=ncc(1,jm)
2734 njc2=ncc(1,3-jm)
2735 if(iqj(njc1).ne.0)then
2736 ncj(jq,njc1)=njc2
2737 else
2738 ncj(3-jq,njc1)=njc2
2739 endif
2740 if(iqj(njc2).ne.0)then
2741 ncj(3-jq,njc2)=njc1
2742 else
2743 ncj(jq,njc2)=njc1
2744 endif
2745 iqc(jm)=10
2746 iqc(3-jm)=10
2747 endif
2748 endif
2749
2750 if(jt.ne.8.and.jt.ne.9)then
2751 jq2=jq
2752 else
2753 jq2=3-jq
2754 endif
2755
2756 call psdeftr(si+q2mass,ept,ey) !lorentz boost to c.m. frame
2757
2758 qt=sqrt(qt2) !p_t
2759 call pscs(bcos,bsin) !cos and sin of the polar angle
2760 if(iabs(iqc(1)).ne.4)then
2761
2762 z=sngl(psutz(dble(si),dble(qt2),dble(qt2)))
2763 if((jt.eq.11.and.jm.eq.1).or.(jt.eq.12.and.jm.eq.2)
2764 $ .or.(jt.eq.13.and.jm.eq.2))z=1-z
2765 wp3=z*sqrt(si)
2766 wm3=qt2/wp3
2767 elseif(jm.eq.1)then
2768 z=sngl(psutz(dble(si),dble(qt2+q2mass),dble(qt2)))
2769 wp3=z*sqrt(si)
2770 wm3=(qt2+q2mass)/wp3
2771 else
2772 z=sngl(psutz(dble(si),dble(qt2),dble(qt2+dble(q2mass))))
2773 wp3=z*sqrt(si)
2774 wm3=qt2/wp3
2775 endif
2776 ep3(1)=.5*(wp3+wm3) !parton 4-momentum
2777 ep3(2)=.5*(wp3-wm3)
2778 ep3(3)=qt*bcos
2779 ep3(4)=qt*bsin
2780 call psdefrot(ep3,s0xh,c0xh,s0h,c0h) !spacial rotation to z-axis
2781
2782 zeta=2. !2=back-to-back emission (angle=pi)
2783 if(iqc(jm).eq.0)then
2784 iq2ini1=9
2785 jo1=iq2ini1
2786 if(zeta.gt.zetacut)jo1=-jo1
2787
2788 else
2789 iq2ini1=iqc(jm)
2790 jo1=iq2ini1
2791 endif
2792 if(iqc(3-jm).eq.0)then
2793 iq2ini2=9
2794 jo2=iq2ini2
2795 if(zeta.gt.zetacut)jo2=-jo2
2796
2797 else
2798 iq2ini2=iqc(3-jm)
2799 jo2=iq2ini2
2800 endif
2801 if(jt.le.10)then
2802 qq1=qt2
2803 qq2=qt2
2804 elseif(jt.le.13)then
2805 qq1=qt2
2806 qq2=0
2807 else
2808 qq1=0
2809 qq2=0
2810 endif
2811
2812 call timsh2(qq1,qq2,sqrt(si),iq2ini1,iq2ini2,jo1,jo2) !final state cascade
2813
2814 if(jt.le.10)then !color connection for the 2nd parton
2815 if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
2816 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2817 if(jt.eq.1)then
2818 nqc(1)=nqc(2)
2819 nqc(2)=ncc(3-jq,3-jm)
2820 elseif(jt.eq.2)then
2821 nqc(2)=ncc(3-jq,jm)
2822 nqc(1)=ncc(jq,3-jm)
2823 elseif(jt.eq.3)then
2824 nqc(1)=nqc(2)
2825 elseif(jt.eq.4)then
2826 nqc(2)=nqc(1)
2827 nqc(1)=ncc(jq,3-jm)
2828 elseif(jt.eq.5)then
2829 nqc(1)=ncc(jq,jm)
2830 elseif(jt.eq.6)then
2831 nqc(2)=ncc(3-jq,3-jm)
2832 nqc(1)=ncc(1,jm)
2833 elseif(jt.eq.7)then
2834 nqc(1)=ncc(1,jm)
2835 elseif(jt.eq.9)then
2836 nqc(1)=ncc(3-jq,3-jm)
2837 elseif(jt.eq.10)then
2838 nqc(1)=nqc(2)
2839 nqc(2)=ncc(1,3-jm)
2840 endif
2841 if(ish.ge.6)write(ifch,*)'jt,jq2,nqc:',jt,jq2,nqc
2842 call psreti(nqc,jq2,2,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstr.
2843 elseif(jt.le.13)then
2844 if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
2845 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2846 ep3(1)=pprt(4,2)
2847 ep3(2)=pprt(3,2)
2848 ep3(3)=pprt(1,2)
2849 ep3(4)=pprt(2,2)
2850 call psrotat(ep3,s0xh,c0xh,s0h,c0h) !special rotation for photon.
2851 call pstrans(ep3,ey,1)
2852 nptl=nptl+1
2853 pptl(1,nptl)=ep3(3)
2854 pptl(2,nptl)=ep3(4)
2855 pptl(3,nptl)=ep3(2)
2856 pptl(4,nptl)=ep3(1)
2857 pptl(5,nptl)=0
2858 idptl(nptl)=10
2859 iorptl(nptl)=iptl
2860 istptl(nptl)=0
2861 jorptl(nptl)=0
2862 do i=1,4
2863 xorptl(i,nptl)=coordo(i)
2864 enddo
2865 tivptl(1,nptl)=coordo(5)
2866 tivptl(2,nptl)=coordo(6)
2867 ityptl(nptl)=71
2868 ifrptl(1,nptl)=0
2869 ifrptl(2,nptl)=0
2870 else
2871 if(ish.ge.6)write(ifch,*)'jt,iqc:',jt,iqc
2872 do j=1,2
2873 ep3(1)=pprt(4,j)
2874 ep3(2)=pprt(3,j)
2875 ep3(3)=pprt(1,j)
2876 ep3(4)=pprt(2,j)
2877 call psrotat(ep3,s0xh,c0xh,s0h,c0h) !special rotation for photon.
2878 call pstrans(ep3,ey,1)
2879 nptl=nptl+1
2880 pptl(1,nptl)=ep3(3)
2881 pptl(2,nptl)=ep3(4)
2882 pptl(3,nptl)=ep3(2)
2883 pptl(4,nptl)=ep3(1)
2884 pptl(5,nptl)=0
2885 idptl(nptl)=10
2886 iorptl(nptl)=iptl
2887 istptl(nptl)=0
2888 jorptl(nptl)=0
2889 do i=1,4
2890 xorptl(i,nptl)=coordo(i)
2891 enddo
2892 tivptl(1,nptl)=coordo(5)
2893 tivptl(2,nptl)=coordo(6)
2894 ityptl(nptl)=72
2895 ifrptl(1,nptl)=0
2896 ifrptl(2,nptl)=0
2897 enddo
2898 endif
2899 call utprix('psabor',ish,ishini,5)
2900 return
2901 end
2902
2903
2904 subroutine psreti(nqc,jort,nfprt,ey,s0xh,c0xh,s0h,c0h)
2905
2906
2907
2908
2909
2910 parameter (ntim=1000)
2911 common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
2912 &,iorprt(ntim),jorprt(ntim),nprtj
2913
2914
2915
2916
2917
2918
2919
2920 parameter (mjstr=20000)
2921 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
2922
2923
2924
2925
2926 common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
2927
2928 dimension ep3(4),nqc(2),ncc(2,ntim),ey(3)
2929 include 'epos.inc'
2930 include 'epos.incsem'
2931 common/cprtx/nprtjx,pprtx(5,2)
2932
2933 if(ish.ge.6)then
2934 write (ifch,*)'nprtj',nprtj
2935 do i=1,nprtj
2936 write (ifch,*)'i,ic,np,ndd',i,idprt(i),iorprt(i),
2937 * idaprt(1,i),idaprt(2,i)
2938 enddo
2939 endif
2940
2941 ncc(1,nfprt)=nqc(1)
2942 if(idprt(nfprt).eq.9)ncc(2,nfprt)=nqc(2)
2943 iprt=nfprt
2944
2945 if(nprtjx.eq.2)then !out Born before timelike cascade
2946 ep3(1)=pprtx(4,iprt)
2947 ep3(2)=pprtx(3,iprt)
2948 ep3(3)=pprtx(1,iprt)
2949 ep3(4)=pprtx(2,iprt)
2950 call psrotat(ep3,s0xh,c0xh,s0h,c0h)
2951 call pstrans(ep3,ey,1)
2952 pprtx(4,iprt)=ep3(1)
2953 pprtx(3,iprt)=ep3(2)
2954 pprtx(1,iprt)=ep3(3)
2955 pprtx(2,iprt)=ep3(4)
2956 endif
2957
2958 1 continue
2959
2960 idau1=idaprt(1,iprt)
2961 idau2=idaprt(2,iprt)
2962 icp=idprt(iprt)
2963
2964 if(ish.ge.6)then
2965 write (ifch,*)'1-iprt,icp,idau1,idau2',iprt,icp,idau1,idau2,
2966 * ncc(1,iprt)
2967 if(icp.eq.9)write (ifch,*)'ncc2',ncc(2,iprt)
2968 endif
2969
2970 if(idau1.ne.0.)then !virtual parton
2971 icd1=idprt(idau1)
2972
2973 if(icp.eq.9)then
2974 if(icd1.ne.9)then !g -> qq~
2975 ncc(1,idau1)=ncc(jort,iprt)
2976 ncc(1,idau2)=ncc(3-jort,iprt)
2977 else !g -> gg
2978 ncc(1,idau1)=ncc(1,iprt)
2979 ncc(2,idau1)=0
2980 ncc(2,idau2)=ncc(2,iprt)
2981 ncc(1,idau2)=0
2982 endif
2983 else !q -> qg
2984 ncc(1,idau1)=0
2985 if(icp*(3-2*jort).gt.0)then
2986 ncc(1,idau2)=ncc(1,iprt)
2987 ncc(2,idau2)=0
2988 else
2989 ncc(1,idau2)=0
2990 ncc(2,idau2)=ncc(1,iprt)
2991 endif
2992 endif
2993 iprt=idau1
2994 goto 1
2995 else
2996
2997 nj=nj+1
2998 ep3(1)=pprt(4,iprt)
2999 ep3(2)=pprt(3,iprt)
3000 ep3(3)=pprt(1,iprt)
3001 ep3(4)=pprt(2,iprt)
3002 call psrotat(ep3,s0xh,c0xh,s0h,c0h)
3003 call pstrans(ep3,ey,1)
3004 do i=1,4
3005 eqj(i,nj)=ep3(i)
3006 enddo
3007
3008 if(icp.eq.9)then
3009 iqj(nj)=0
3010 elseif(iabs(icp).lt.3)then
3011 iqj(nj)=icp
3012 elseif(iabs(icp).eq.3)then
3013 iqj(nj)=icp*4/3
3014 else
3015 iqj(nj)=icp*10
3016 endif
3017
3018 ioj(nj)=jorprt(iprt) !flavor of mother parton
3019 q2j(nj)=q2prt(iprt)
3020
3021 if(iqj(nj).ne.0)then
3022 njc=ncc(1,iprt)
3023 if(njc.ne.0)then
3024 ncj(1,nj)=njc
3025 iqc=iqj(njc)
3026 if(iqc.ne.0)then
3027 ncj(1,njc)=nj
3028 else
3029 if(iqj(nj).gt.0)then
3030 ncj(2,njc)=nj
3031 else
3032 ncj(1,njc)=nj
3033 endif
3034 endif
3035 else
3036 ncc(1,iprt)=nj
3037 endif
3038 else
3039
3040 do m=1,2
3041 if(jort.eq.1)then
3042 m1=m
3043 else
3044 m1=3-m
3045 endif
3046 njc=ncc(m1,iprt)
3047 if(njc.ne.0)then
3048 ncj(m,nj)=njc
3049 iqc=iqj(njc)
3050 if(iqc.ne.0)then
3051 ncj(1,njc)=nj
3052 else
3053 ncj(3-m,njc)=nj
3054 endif
3055 else
3056 ncc(m1,iprt)=nj
3057 endif
3058 enddo
3059 endif
3060 if(ish.ge.6)then
3061 write (ifch,*)'jet-nj,iprt,icp,iqj(nj),ioj(nj),ncj',
3062 * nj,iprt,icp,iqj(nj),ioj(nj),ncj(1,nj)
3063 if(icp.eq.9)write (ifch,*)'ncj2',ncj(2,nj)
3064 endif
3065 endif
3066
3067 2 continue
3068 if(iprt.ne.nfprt)then
3069 icp=idprt(iprt)
3070 ipar=iorprt(iprt)
3071 idau1=idaprt(1,ipar)
3072 idau2=idaprt(2,ipar)
3073 if(ish.ge.6)then
3074 write (ifch,*)'2-iprt,icp,idau1,idau2,ipar',
3075 * iprt,icp,idau1,idau2,ipar,ncc(1,iprt)
3076 if(icp.eq.9)write (ifch,*)ncc(2,iprt)
3077 endif
3078
3079 if(idau1.eq.iprt)then
3080 if(icp.eq.9)then !g -> gg
3081 ncc(1,ipar)=ncc(1,iprt)
3082 ncc(1,idau2)=ncc(2,iprt)
3083 else
3084 icpar=idprt(ipar)
3085 if(icpar.eq.9)then !g -> qq~
3086 ncc(jort,ipar)=ncc(1,iprt)
3087 else !q -> qg
3088 if(icp*(3-2*jort).gt.0)then
3089 ncc(2,idau2)=ncc(1,iprt)
3090 else
3091 ncc(1,idau2)=ncc(1,iprt)
3092 endif
3093 endif
3094 endif
3095 iprt=idau2
3096 goto 1
3097
3098 else
3099 if(icp.eq.9)then
3100 icpar=idprt(ipar)
3101 if(icpar.eq.9)then !g -> gg
3102 ncc(2,ipar)=ncc(2,iprt)
3103 ncc(2,idau1)=ncc(1,iprt)
3104 else !q -> qg
3105 if(icpar*(3-2*jort).gt.0)then
3106 ncc(1,ipar)=ncc(1,iprt)
3107 ncc(1,idau1)=ncc(2,iprt)
3108 else
3109 ncc(1,ipar)=ncc(2,iprt)
3110 ncc(1,idau1)=ncc(1,iprt)
3111 endif
3112 endif
3113 else
3114 ncc(3-jort,ipar)=ncc(1,iprt)
3115 endif
3116 iprt=ipar
3117 goto 2
3118 endif
3119 else
3120 if(ish.ge.6)write (ifch,*)'3-iprt,ncc',iprt,ncc(1,iprt)
3121 nqc(1)=ncc(1,nfprt)
3122 if(idprt(nfprt).eq.9)nqc(2)=ncc(2,nfprt)
3123 endif
3124 return
3125 end
3126
