1 c------------------------------------------------------------------------
2 function ffsigiut(xx1,xx2,jpp,je1,je2)
3 c------------------------------------------------------------------------
5 c \int(dt) \int(du) ffsig *s/sh**3 *2*pi*alpha**2 *delta(uh+th+sh)
7 c-----------------------------------------------------------------------
8 common /ar3/ x1(7),a1(7)
11 double precision tmin,tmax,t,sh2,sqrtq2s
17 if(sh.le.4.*q2min)return
19 c tmin=sh/2-sqrt(sh*sh/4-q2min*sh)
20 sqrtq2s=sqrt(dble(q2min*sh))
21 tmin=sh2-sqrt((sh2-sqrtq2s)*(sh2+sqrtq2s))
25 t=2d0*tmin/(1d0+tmin/tmax-dble(tgss(ig,i)*(2*m-3))
27 qq=sngl(t*(1d0-t/dble(sh)))
28 ft=ffsigj(sngl(t),qq,xx1,xx2,jpp,je1,je2)/sh**3
29 * * (2*pi*pssalf(qq/qcdlam))**2
30 ffsigiut=ffsigiut+wgss(ig,i)*ft*sngl(t)**2
34 * *0.5*sngl(1d0/tmin-1d0/tmax)
36 * /2 !CS for parton pair
40 c-----------------------------------------------------------------------
41 function ffsigj(t,qt,x1,x2,jpp,je1,je2)
42 c-----------------------------------------------------------------------
44 c \sum x1*f_i(x1,qt) * x2*f_k(x2,qt) * B_ik
46 c B_ik = psbori = contribution to Born xsection:
48 c = s/pi * delta(s+t+u) * 2*pi*alpha**2 /s**2 * B_ik
50 c qt = virtuality scale
51 c x1, x2 = light cone momentum fractions
53 c x*f_j(x,qt) = function fparton(x,qt,j)
55 c-----------------------------------------------------------------------
56 c jpp: type of Pomeron
66 c-----------------------------------------------------------------------
75 sea1=pifpartone(x1,qt,-1,je1,ji1)
76 g1= pifpartone(x1,qt, 0,je1,ji1)
77 uv1= pifpartone(x1,qt, 1,je1,ji1)
78 dv1= pifpartone(x1,qt, 2,je1,ji1)
79 sea2=pifpartone(x2,qt,-1,je2,ji2)
80 g2= pifpartone(x2,qt, 0,je2,ji2)
81 uv2= pifpartone(x2,qt, 1,je2,ji2)
82 dv2= pifpartone(x2,qt, 2,je2,ji2)
84 sea1=pifpartone(x1,qt,-1,je1,1)+pifpartone(x1,qt,-1,je1,2)
85 g1= pifpartone(x1,qt, 0,je1,1)+pifpartone(x1,qt, 0,je1,2)
86 uv1= pifpartone(x1,qt, 1,je1,1)+pifpartone(x1,qt, 1,je1,2)
87 dv1= pifpartone(x1,qt, 2,je1,1)+pifpartone(x1,qt, 2,je1,2)
88 sea2=pifpartone(x2,qt,-1,je2,1)+pifpartone(x2,qt,-1,je2,2)
89 g2= pifpartone(x2,qt, 0,je2,1)+pifpartone(x2,qt, 0,je2,2)
90 uv2= pifpartone(x2,qt, 1,je2,1)+pifpartone(x2,qt, 1,je2,2)
91 dv2= pifpartone(x2,qt, 2,je2,1)+pifpartone(x2,qt, 2,je2,2)
94 ffsigj= ffborn(s,t, g1*g2 !gg
96 * ,(uv1+dv1+2.*naflav*sea1)*g2+g1*(uv2+dv2+2.*naflav*sea2) !gq
98 * ,(uv1+sea1)*(uv2+sea2) !qq
99 * +(dv1+sea1)*(dv2+sea2)+sea1*sea2*(naflav-1)*2.
101 * ,(uv1+sea1)*sea2+(uv2+sea2)*sea1 !qa
102 * +(dv1+sea1)*sea2+(dv2+sea2)*sea1+sea1*sea2*(naflav-2)*2.
104 * ,dv1*uv2+dv2*uv1+(uv2+dv2)*sea1*(naflav-1)*2. !qqp
105 * +(uv1+dv1)*sea2*(naflav-1)*2.
106 * +sea1*sea2*naflav*(naflav-1)*4.
111 c-----------------------------------------------------------------------
112 function ffsig(t,qt,x1,x2) !former psjy
113 c-----------------------------------------------------------------------
114 include 'epos.incsem'
119 g1= pifpartone(x1,qt, 0,2,1)+pifpartone(x1,qt, 0,2,2)
120 uv1= pifpartone(x1,qt, 1,2,1)+pifpartone(x1,qt, 1,2,2)
121 dv1= pifpartone(x1,qt, 2,2,1)+pifpartone(x1,qt, 2,2,2)
122 sea1=pifpartone(x1,qt,-1,2,1)+pifpartone(x1,qt,-1,2,2)
123 g2= pifpartone(x2,qt, 0,2,1)+pifpartone(x2,qt, 0,2,2)
124 uv2= pifpartone(x2,qt, 1,2,1)+pifpartone(x2,qt, 1,2,2)
125 dv2= pifpartone(x2,qt, 2,2,1)+pifpartone(x2,qt, 2,2,2)
126 sea2=pifpartone(x2,qt,-1,2,1)+pifpartone(x2,qt,-1,2,2)
128 ffsig= ffborn(s,t, g1*g2 !gg
130 * ,(uv1+dv1+2.*naflav*sea1)*g2+g1*(uv2+dv2+2.*naflav*sea2) !gq
132 * ,(uv1+sea1)*(uv2+sea2) !qq
133 * +(dv1+sea1)*(dv2+sea2)+sea1*sea2*(naflav-1)*2.
135 * ,(uv1+sea1)*sea2+(uv2+sea2)*sea1 !qa
136 * +(dv1+sea1)*sea2+(dv2+sea2)*sea1+sea1*sea2*(naflav-2)*2.
138 * ,dv1*uv2+dv2*uv1+(uv2+dv2)*sea1*(naflav-1)*2. !qqp
139 * +(uv1+dv1)*sea2*(naflav-1)*2.
140 * +sea1*sea2*naflav*(naflav-1)*4.
145 c------------------------------------------------------------------------
146 function ffborn(s,t,gg,gq,qq,qa,qqp)
147 c------------------------------------------------------------------------
150 *( psbori(s,t,0,0,1)+psbori(s,s-t,0,0,1)
151 * +psbori(s,t,0,0,2)+psbori(s,s-t,0,0,2)) /2. *gg !gg
153 *+(psbori(s,t,0,1,1)+psbori(s,s-t,0,1,1)) *gq !gq
155 *+(psbori(s,t,1,1,1)+psbori(s,s-t,1,1,1))/2. *qq !qq
157 *+(psbori(s,t,1,-1,1)+psbori(s,s-t,1,-1,1)+psbori(s,t,1,-1,2)+
158 * psbori(s,s-t,1,-1,2)+psbori(s,t,1,-1,3)+psbori(s,s-t,1,-1,3)) !qa
161 *+(psbori(s,t,1,2,1)+psbori(s,s-t,1,2,1)) *qqp !qq'
165 c-----------------------------------------------------------------------
166 function pifpartone(xx,qq,j,je,ji) ! pol interpolation of partone
167 c-----------------------------------------------------------------------
168 include 'epos.incsem'
170 common/tabfptn/kxxmax,kqqmax,fptn(20,20,-1:2,0:2,2)
172 common /cpifpartone/npifpartone
174 npifpartone=npifpartone+1
175 if(npifpartone.eq.1)call MakeFpartonTable
181 xxk=1.+log(xx/xxmin)/log(xxmax/xxmin)*(kxxmax-1)
182 qqk=1.+log(qq/q2min)/log(qqmax/q2min)*(kqqmax-1)
187 if(kxx.gt.(kxxmax-2))kxx=kxxmax-2
188 if(kqq.gt.(kqqmax-2))kqq=kqqmax-2
191 wi(3)=wi(2)*(wi(2)-1.)*.5
196 wj(3)=wj(2)*(wj(2)-1.)*.5
202 pifpartone=pifpartone+fptn(kxx+kx-1,kqq+kq-1,j,je,ji)
208 c-----------------------------------------------------------------------
209 subroutine MakeFpartonTable
210 c-----------------------------------------------------------------------
211 include 'epos.incsem'
213 common/tabfptn/kxxmax,kqqmax,fptn(20,20,-1:2,0:2,2)
214 write (*,'(a,$)')'(Fparton table'
225 xx=xxmin*(xxmax/xxmin)**((kxx-1.)/(kxxmax-1.))
227 qq=q2min*(qqmax/q2min)**((kqq-1.)/(kqqmax-1.))
228 fptn(kxx,kqq,j,je,ji)= fpartone(xx,qq,j,je,ji)
234 write (*,'(a,$)')'done)'
237 c------------------------------------------------------------------------
238 function fpartone(xx,qq,j,je,ji) !former pspdf0 (sha)
239 c-----------------------------------------------------------------------
241 c parton distribution function for proton ( actually x*f(x) !!!!!!! )
243 c xx = light cone momentum fraction
244 c qq = virtuality scale
246 c -1 ... sea (distribution function per flavor)
254 c ji = initial parton type
257 c-----------------------------------------------------------------------
258 double precision z,xmin,xm,zx,psuds
259 common/ar3/ x1(7),a1(7)
261 include 'epos.incsem'
266 c ...... f_0 * sudakov.........
268 if(j.eq.0.and.ji.eq.1)then
269 fpartone=fzeroGlu(xx,2,1) !hadron class 2, projectile side
270 elseif((j.eq.1.or.j.eq.2).and.ji.eq.2)then
271 fpartone=psdfh4(xx,q2min,0.,2,j)
272 elseif(j.eq.-1.and.ji.eq.1)then
273 fpartone=fzeroSea(xx,2,1)
275 fpartone=fpartone*sngl(psuds(qq,j)/psuds(q2min,j))
278 c......... integral f_0 E_qcd............
281 xmin=dble(xx)/(1.d0-dble(q2ini/qq))
287 !numerical integration xm -> 1
291 zx=1.d0-(1.d0-xm)*(.5d0+(dble(m)-1.5d0)*dble(x1(i)))**.25d0
294 gl=fzeroGlu(sngl(zx),2,1)
295 uv=psdfh4(sngl(zx),q2min,0.,2,1)
296 dv=psdfh4(sngl(zx),q2min,0.,2,2)
297 sea=fzeroSea(sngl(zx),2,1)
302 * fz=gl *psevi(q2min,qq,z,1,1)
303 * +sea*psevi(q2min,qq,z,2,1) !ccccc
305 * fz=(uv+dv)*psevi(q2min,qq,z,2,1)
306 elseif(j.eq.1.and.ji.eq.2)then
307 fz=psevi(q2min,qq,z,3,2)*uv
308 elseif(j.eq.2.and.ji.eq.2)then
309 fz=psevi(q2min,qq,z,3,2)*dv
311 akns=psevi(q2min,qq,z,3,2) !nonsinglet contribution
312 aks=(psevi(q2min,qq,z,2,2)-akns) !singlet contribution
314 * fz=psevi(q2min,qq,z,1,2)*gl
315 * +sea*aks+sea*akns !ccccc
319 dpd1=dpd1+a1(i)*fz/sngl(zx)**2/sngl(1.d0-zx)**3
322 dpd1=dpd1*sngl(1.d0-xm)**4/8.*xx
324 !numerical integration xmin -> xm
330 & *((xmin-xx)/(xm-xx))**(.5d0-(dble(m)-1.5d0)*dble(x1(i)))
333 gl=fzeroGlu(sngl(zx),2,1)
334 uv=psdfh4(sngl(zx),q2min,0.,2,1)
335 dv=psdfh4(sngl(zx),q2min,0.,2,2)
336 sea=fzeroSea(sngl(zx),2,1)
341 * fz=gl *psevi(q2min,qq,z,1,1)
342 * +sea*psevi(q2min,qq,z,2,1) !ccccc
344 * fz=(uv+dv)*psevi(q2min,qq,z,2,1)
345 elseif(j.eq.1.and.ji.eq.2)then
346 fz=psevi(q2min,qq,z,3,2)*uv
347 elseif(j.eq.2.and.ji.eq.2)then
348 fz=psevi(q2min,qq,z,3,2)*dv
350 akns=psevi(q2min,qq,z,3,2) !nonsinglet contribution
351 aks=(psevi(q2min,qq,z,2,2)-akns) !singlet contribution
353 * fz=psevi(q2min,qq,z,1,2)*gl
354 * +sea*aks+sea*akns !ccccc
358 dpd2=dpd2+a1(i)*fz*sngl((1.d0-xx/zx)/zx)
361 dpd2=dpd2*sngl(log((xm-xx)/(xmin-xx))*.5d0*xx)
363 fpartone=fpartone+dpd2+dpd1
367 if(j.lt.0)fpartone=fpartone/naflav/2.
371 c------------------------------------------------------------------------
372 function fparton(xx,qq,j) !former pspdf0 (sha)
373 c-----------------------------------------------------------------------
375 c parton distribution function for proton ( actually x*f(x) !!!!!!! )
377 c xx = light cone momentum fraction
378 c qq = virtuality scale
380 c -1 ... sea (dsistribution fuction per flavor)
385 c-----------------------------------------------------------------------
386 c (see pages 105 - 107 of our report)
388 c fparton(xx) = xx * f(xx) !!!!!
390 c f_j(xx,qq) = \sum_k \int(xx<x<1) dx/x f0_k(x) Eqcd_k_j(xx/x,qq)
392 c f0_k = fzeroGlu or fzeroSea
394 c Eqcd=E~qcd+delta*sudakov, E~qcd: at least one emission
396 c-----------------------------------------------------------------------
397 double precision z,xmin,xm,zx,psuds
398 common/ar3/ x1(7),a1(7)
400 include 'epos.incsem'
402 c ...... f_0 * sudakov.........
405 fparton=fzeroGlu(xx,2,1)
406 elseif(j.eq.1.or.j.eq.2)then
407 fparton=psdfh4(xx,q2min,0.,2,j)
409 fparton=fzeroSea(xx,2,1)
411 fparton=fparton*sngl(psuds(qq,j)/psuds(q2min,j))
413 c......... integral f_0 E_qcd............
415 xmin=xx/(1.d0-dble(q2ini/qq))
421 !numerical integration xm -> 1
425 zx=1.d0-(1.d0-xm)*(.5d0+(dble(m)-1.5d0)*dble(x1(i)))**.25d0
428 gl=fzeroGlu(sngl(zx),2,1)
429 uv=psdfh4(sngl(zx),q2min,0.,2,1)
430 dv=psdfh4(sngl(zx),q2min,0.,2,2)
431 sea=fzeroSea(sngl(zx),2,1)
434 fz=psevi(q2min,qq,z,1,1)*gl
435 * +(uv+dv+sea)*psevi(q2min,qq,z,2,1)
437 fz=psevi(q2min,qq,z,3,2)*uv
439 fz=psevi(q2min,qq,z,3,2)*dv
441 akns=psevi(q2min,qq,z,3,2) !nonsinglet contribution
442 aks=(psevi(q2min,qq,z,2,2)-akns) !singlet contribution
443 fz=(psevi(q2min,qq,z,1,2)*gl+(uv+dv+sea)*aks+sea*akns)
445 dpd1=dpd1+a1(i)*fz/sngl(zx)**2/sngl(1.d0-zx)**3
448 dpd1=dpd1*sngl((1.d0-xm)**4/8.*xx)
450 !numerical integration xmin -> xm
455 zx=xx+(xm-xx)*((xmin-xx)/(xm-xx))
456 * **(.5d0-(dble(m)-1.5)*dble(x1(i)))
459 gl=fzeroGlu(sngl(zx),2,1)
460 uv=psdfh4(sngl(zx),q2min,0.,2,1)
461 dv=psdfh4(sngl(zx),q2min,0.,2,2)
462 sea=fzeroSea(sngl(zx),2,1)
465 fz=psevi(q2min,qq,z,1,1)*gl+(uv+dv+sea)*
466 * psevi(q2min,qq,z,2,1)
468 fz=psevi(q2min,qq,z,3,2)*uv
470 fz=psevi(q2min,qq,z,3,2)*dv
472 akns=psevi(q2min,qq,z,3,2) !nonsinglet contribution
473 aks=(psevi(q2min,qq,z,2,2)-akns) !singlet contribution
474 fz=(psevi(q2min,qq,z,1,2)*gl+(uv+dv+sea)*aks+sea*akns)
476 dpd2=dpd2+a1(i)*fz*sngl((1.d0-xx/zx)/zx)
479 dpd2=dpd2*sngl(log((xm-xx)/(xmin-xx))*.5d0*xx)
481 fparton=fparton+dpd2+dpd1
483 if(j.lt.0)fparton=fparton/naflav/2.
487 c------------------------------------------------------------------------
488 function fzeroGlu(z,k,ipt)
489 c-----------------------------------------------------------------------
493 c f = F & EsoftGluon &=convolution
495 c F(x) = alpff(k)*x**betff(ipt)*(1-x)**alplea(k)
497 c EsoftGluon(x) = x**(-1-dels) * EsoftGluonTil(x)
501 c ipt - 1=proj 2=targ
502 c-----------------------------------------------------------------------
503 double precision xpmin,xp
505 common /ar3/ x1(7),a1(7)
506 include 'epos.incsem'
510 xpmin=xpmin**(1+betff(ipt)+dels)
513 xp=(.5*(1.+xpmin+(2*m-3)*x1(i)*(1.-xpmin)))**(1./
514 * (1+betff(ipt)+dels))
516 fzeroGlu=fzeroGlu+a1(i)*(1.-xp)**alplea(k)*EsoftGluonTil(zz)
519 fzeroGlu=fzeroGlu*.5*(1.-xpmin)/(1+betff(ipt)+dels)
521 fzeroGlu=fzeroGlu *alpff(k) *z**(-dels)
525 c------------------------------------------------------------------------
526 function fzeroSea(z,k,ipt)
527 c-----------------------------------------------------------------------
531 c f = F & EsoftQuark &=convolution
533 c F(x) = alpff(k)*x**betff(ipt)*(1-x)**alplea(k)
535 c EsoftQuark(x) = x**(-1-dels) * EsoftQuarkTil(x)
537 c z - light cone x of the quark,
539 c-----------------------------------------------------------------------
540 double precision xpmin,xp
541 common /ar3/ x1(7),a1(7)
543 include 'epos.incsem'
547 xpmin=xpmin**(1+betff(ipt)+dels)
550 xp=(.5*(1.+xpmin+(2*m-3)*x1(i)*(1.-xpmin)))**(1./
551 * (1+betff(ipt)+dels))
553 fzeroSea=fzeroSea+a1(i)*(1.-xp)**alplea(k)*EsoftQuarkTil(zz)
556 fzeroSea=fzeroSea*.5*(1.-xpmin)/(1+betff(ipt)+dels)
558 fzeroSea=fzeroSea *alpff(k) *z**(-dels)
562 c------------------------------------------------------------------------
563 function EsoftGluonTil(zz)
564 c-----------------------------------------------------------------------
565 c EsoftGluon = zz^(-1-dels) * EsoftGluonTil
566 c-----------------------------------------------------------------------
568 include 'epos.incsem'
569 EsoftGluonTil=gamsoft*(1-glusea)*(1.-zz)**betpom
572 c------------------------------------------------------------------------
573 function EsoftQuarkTil(zz)
574 c-----------------------------------------------------------------------
575 c EsoftQuark = zz^(-1-dels) * EsoftQuarkTil
576 c-----------------------------------------------------------------------
577 double precision zmin,z
578 common /ar3/ x1(7),a1(7)
580 include 'epos.incsem'
587 z=(.5d0*(1.+zmin+(2*m-3)*x1(i)*(1.d0-zmin)))
588 * **(1.d0/(1.d0+dels))
589 EsoftQuarkTil=EsoftQuarkTil+a1(i)*max(1.d-5,(1.d0-zz/z))**betpom
593 EsoftQuarkTil=EsoftQuarkTil*1.5*(1.d0-zmin)/(1.+dels)
595 EsoftQuarkTil=gamsoft*glusea*EsoftQuarkTil
599 c------------------------------------------------------------------------
600 function EsoftQZero(zz) ! former psftilf
601 c-----------------------------------------------------------------------
603 c EsoftQuark = EsoftQZero * wsplit * z^(-1-dels) * gamsoft
605 c zz - ratio of the quark and pomeron light cone x (zz=x_G/x_P)
606 c integration over quark to gluon light cone momentum ratio (z=x/x_G):
608 c EsoftQZero = int(dz) z^dels * (1-zz/z)^betpom * P_qG(z)
610 c-----------------------------------------------------------------------
611 double precision zmin,z
612 common /ar3/ x1(7),a1(7)
613 include 'epos.incsem'
620 z=(.5d0*(1.+zmin+(2*m-3)*x1(i)*(1.d0-zmin)))
621 * **(1.d0/(1.d0+dels))
622 EsoftQZero=EsoftQZero+a1(i)*max(1.d-5,(1.d0-zz/z))**betpom
626 EsoftQZero=EsoftQZero*1.5*(1.d0-zmin)/(1.+dels) !1.5=naflav/2 at Q0
630 c------------------------------------------------------------------------
631 function ffsigi(qq,y0) !former psjx1 (sto)
632 c------------------------------------------------------------------------
634 c dsigma/dpt_jet = \int dy \int dx1 ffsig(x1,x2(x1))
636 c x1=xplus, x2=xminus
637 c x2=x2(x1) due to u+t+s=0
638 c ( s=x1*x2*spp, t/spp=-x1*xt*exp(-y)/2, u/spp=-x2*xt*exp(y)/2 )
640 c qq = pt**2, xt=2.*sqrt(qq/s)
641 c rapidity range: 0 to y0
643 c ffsig = function ffsig(t,qq,x1,x2)
645 c-----------------------------------------------------------------------
646 include 'epos.incsem'
648 double precision xx1,xx2,xt,ymax,ymin,y,xmin,xmax
654 if(qq.lt.q2min)return
655 xt=2d0*sqrt(dble(qq)/dble(s))
656 ymax=min(dble(y0),log(1d0/xt+sqrt((1d0/xt-1d0)*(1d0/xt+1d0))))
657 ymin=-ymax !final result must be divided by 2
660 y=.5d0*(ymax+ymin+(ymin-ymax)*dble((2*m-3)*tgss(ig,i)))
661 !for xx1-integration, use variable x=xx1-xt*exp(y)/2.,with xmin<x<xmax
662 xmin=xt**2/2.d0/(2.d0-xt*exp(-y)) !condition x2<1
663 xmax=1.d0-xt*exp(y)/2.d0 !condition x1<1
667 xx1=xt*exp(y)/2.d0+xmin*(xmax/xmin)**dble(.5
668 & +tgss(ig1,i1)*(m1-1.5))
669 xx2=xt*exp(-y)*xx1/(2.d0*xx1-xt*exp(y))
674 t=sh/2.*(1.-sqrt(aa)) !formula in parton-parton cms
675 ft=ffsig(t,qq,sngl(xx1),sngl(xx2))
676 fx=fx+wgss(ig1,i1)*ft/sh**2
679 fx=fx*0.5*sngl(log(xmax/xmin)) !dx/x=0.5*log(xmax/xmin)dt (gauss)
680 ffsigi=ffsigi+wgss(ig,i)*fx
683 ffsigi=ffsigi*0.5*sngl(ymax-ymin) !dy=0.5*(ymax-ymin)dt (gauss)
684 * *2*pi*(2*pi*pssalf(qq/qcdlam))**2 !alpha = 2*pi*pssalf
685 * *2*sqrt(qq) !d2pt=2*pi*pt*dpt
686 * /2 ! y interval 2 * Delta_y
687 * /2 ! condition t < sqrt(s)/2,
688 ! since t > sqrt(s)/2 is automatically included,
693 c------------------------------------------------------------------------
694 function psbori(s,t,j,l,n)
695 c-----------------------------------------------------------------------
696 c contribution to the born cross-section:
698 c dsigmaBorn/d2pt/dy = s/pi * delta(s+t+u) * 2*pi*alpha**2 /s**2 *psbori
700 c s - c.m. energy squared for the born scattering,
701 c t - invariant variable for the born scattering |(p1-p3)**2|,
702 c j - parton type at current end of the ladder (0 - g, 1,-1,2,... - q)
703 c l - parton type at opposite end of the ladder (0 - g, 1,-1,2,... - q)
704 c n - subprocess number
705 c-----------------------------------------------------------------------
706 include 'epos.incsem'
712 if(iabs(j).ne.4)then !light quarks and gluons
714 if(j.eq.0.and.l.eq.0)then !gg->gg
715 psbori=(3.-t*u/s**2+s*u/t**2+s*t/u**2)*4.5
716 elseif(j*l.eq.0)then !gq->gq
717 psbori=(s**2+u**2)/t**2+(s/u+u/s)/2.25
718 elseif(j.eq.l)then !qq->qq
719 psbori=((s**2+u**2)/t**2+(s**2+t**2)/u**2)/2.25
721 elseif(j.eq.-l)then !qq~->qq~
722 psbori=((s**2+u**2)/t**2+(u**2+t**2)/s**2)/2.25
725 psbori=(s**2+u**2)/t**2/2.25
728 if(j.eq.0.and.l.eq.0)then !gg->qq~
729 psbori=.5*(t/u+u/t)-1.125*(t*t+u*u)/s**2
730 elseif(j.eq.-l)then !qq~->q'q'~
731 psbori=(t*t+u*u)/s**2/1.125
736 if(j.ne.0.and.j.eq.-l)then !qq~->gg
737 psbori=32./27.*(t/u+u/t)-(t*t+u*u)/s**2/.375
742 c............ n=4 for photon product processes, make e_q**2 =2/9.,
743 c the average value of charge squared for all types of quarks.
745 if(j.ne.0.and.j.eq.-l)then !qq~->g+gamma
746 psbori=16*factgam*(u/t+t/u)/81.
747 elseif (j*l.eq.0.and.j+l.ne.0) then !q(q~)g->q(q~)+gamma
748 psbori=2*factgam*(u/s+s/u)/27.
753 if(j.ne.0.and.j.eq.-l)then !qq~->gamma+gamma
754 psbori=4*factgam*(t/u+u/t)/81.
760 elseif(n.eq.1)then !c-quark
762 if(l.eq.0)then !cg->cg
764 psbori=(s**2+u**2)/t**2+(s/u+u/s)/2.25
765 * -4.*qcmass**2/t+xm*(xm*t**2-t)/.5625+4.*qcmass**2*xm
767 psbori=(s**2+u**2)/t**2/2.25-qcmass**2/t/1.125
778 c-----------------------------------------------------------------------
779 double precision function om51p(sy,xh,yp,b,iqq)
780 c-----------------------------------------------------------------------
782 c xh - fraction of the energy squared s for the pomeron;
783 c yp - rapidity for the pomeron;
784 c b - impact parameter between the pomeron ends;
785 c iqq =-1 - 0+1+2+3+4,
786 c iqq = 0 - soft pomeron,
791 c iqq = 5 - soft(int)|b,
792 c iqq = 6 - gg(int)|b,
793 c iqq = 7 - soft(proj)|b,
794 c iqq = 8 - gg(proj)|b,
795 c iqq = 9 - qg(proj)|b,
796 c iqq = 10 - total fro-uncut integrated,
797 c iqq = 11 - total uncut integrated,
798 c iqq = 12 - soft(int),
799 c iqq = 13 - gg(int),
800 c iqq = 14 - <b^2*soft(int)>,
801 c iqq = 15 - <b^2*gg(int)>,
802 c iqq = 16 - soft(proj-int),
803 c iqq = 17 - gg(proj-int),
804 c iqq = 18 - qg(proj-int),
805 c iqq = 19 - <b^2*soft(proj)>,
806 c iqq = 20 - <b^2*gg(proj)>,
807 c iqq = 21 - <b^2*qg(proj)>
808 c-----------------------------------------------------------------------
809 double precision xh,yp!,coefom1,coefom2
810 common /psar7/ delx,alam3p,gam3p
811 common /psar37/ coefom1,coefom2
813 include 'epos.incsem'
821 rp=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy))
822 zb=exp(-b**2/(4.*.0389*rp))
823 rh=r2had(iclpro)+r2had(icltar)
825 if(iqq.eq.0)then !soft
826 c rp=r2hads(iclpro)+r2hads(icltar)+slopoms*log(max(1.,sy))
827 zb=exp(-b**2/(4.*.0389*rp))
828 om51p=chad(iclpro)*chad(icltar)*gamhads(iclpro)
829 * *gamhads(icltar)*sy**dels*(xp*xm)**(-alppar)*zb/rp
830 elseif(iqq.le.4)then !gg,qg,gq,qq
831 om51p=psvin(sy,xp,xm,zb,iqq)
832 elseif(iqq.eq.5)then !soft(int)|b
833 c rh=alam3p+slopoms*log(max(1.,sy))
834 om51p=sy**dels*zb**(rp/rh)/rh
835 elseif(iqq.eq.6)then !gg(int)|b
836 om51p=psvin(sy,xp,xm,zb,14)
837 elseif(iqq.eq.7)then !soft(proj)b
838 c rh=r2hads(iclpro)+.5*alam3p+slopoms*log(max(1.,sy))
839 om51p=chad(iclpro)*gamhads(iclpro)*sy**dels
840 * *xp**(-alppar)*zb**(rp/rh)/rh
841 elseif(iqq.eq.8)then !gg(proj)b
842 om51p=psvin(sy,xp,xm,zb,16)
843 elseif(iqq.eq.9)then !qg(proj)b
844 om51p=psvin(sy,xp,xm,zb,18)
845 elseif(iqq.eq.10)then !total fro-uncut integrated
848 elseif(iqq.eq.11)then !total uncut integrated
849 om51p=psvin(sy,xp,xm,zb,9)
850 c om51p=om51p+dble(coefom1)/2.d0*om51p**2+dble(coefom2)/6.d0*om51p**3 !!!!!!!!!!
851 c if(om51p.gt.100.d0)om51p=100.d0
852 elseif(iqq.eq.12)then !soft(int)
853 om51p=sy**dels*4.*.0389
854 elseif(iqq.eq.13)then !gg(int)
855 om51p=psvin(sy,xp,xm,zb,5)
856 elseif(iqq.eq.14)then !<b^2*soft(int)>
857 c rh=alam3p+slopoms*log(max(1.,sy))
858 om51p=sy**dels*rh*(4.*.0389)**2
859 elseif(iqq.eq.15)then !<b^2*gg(int)>
860 om51p=psvin(sy,xp,xm,zb,15)
861 elseif(iqq.eq.16)then !soft(proj-int)
862 om51p=chad(iclpro)*gamhads(iclpro)*sy**dels
863 * *xp**(-alppar)*4.*.0389
864 elseif(iqq.eq.17)then !gg(proj-int)
865 om51p=psvin(sy,xp,xm,zb,6)
866 elseif(iqq.eq.18)then !qg(proj-int)
867 om51p=psvin(sy,xp,xm,zb,7)
868 elseif(iqq.eq.19)then !<b^2*soft(proj)>
869 c rh=r2hads(iclpro)+.5*alam3p+slopoms*log(max(1.,sy))
870 om51p=chad(iclpro)*gamhads(iclpro)*sy**dels
871 * *xp**(-alppar)*rh*(4.*.0389)**2
872 elseif(iqq.eq.20)then !<b^2*gg(proj)>
873 om51p=psvin(sy,xp,xm,zb,17)
874 elseif(iqq.eq.21)then !<b^2*qg(proj)>
875 om51p=psvin(sy,xp,xm,zb,19)
881 cc-----------------------------------------------------------------------
882 c double precision function om2p(xh,yp,xprem0,xmrem0,b,iqq)
883 cc-----------------------------------------------------------------------
884 cc om2p - chi~(x,y) for cut pomeron
885 cc xh - fraction of the energy squared s for the pomeron;
886 cc yp - rapidity for the pomeron;
887 cc xprem - x+ for the projectile remnant;
888 cc xmrem - x- for the target remnant;
889 cc b - impact parameter between the pomeron ends;
894 cc iqq = 3 - 1-cut(soft),
898 cc iqq = 7 - 1+(difr)
899 cc iqq = -7 - 1-(difr)
900 cc-----------------------------------------------------------------------
901 c double precision xh,yp,xprem0,xmrem0
903 c include 'epos.incsem'
909 c xp=dsqrt(xh)*dexp(yp)
915 c rp=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy))
916 c zb=exp(-b**2/(4.*.0389*rp))
919 c om2p=psvy(xp,xprem,xm,xmrem,b,2)
920 c * +psvy(xp,xprem,xm,xmrem,b,-2)
921 c * +psvy(xp,xprem,xm,xmrem,b,3)
922 c * +psvy(xp,xprem,xm,xmrem,b,-3)
923 c * +psvy(xp,xprem,xm,xmrem,b,9)
924 c * +psvy(xp,xprem,xm,xmrem,b,-9)
925 c * +psvx(xp,xprem,xm,xmrem,b,1)
926 c * +psvx(xp,xprem,xm,xmrem,b,2)
927 c * +psvx(xp,xprem,xm,xmrem,b,-2)
928 c * +psvx(xp,xprem,xm,xmrem,b,6)
929 c * +psvx(xp,xprem,xm,xmrem,b,-6)
930 c om2p=om2p+(chad(iclpro)*chad(icltar)*gamhad(iclpro)
931 c * *gamhad(icltar)*sy**dels*(xp*xm)**(-alppar)*zb/rp
932 c * +psvin(sy,xp,xm,zb,1)+psvin(sy,xp,xm,zb,2)
933 c * +psvin(sy,xp,xm,zb,3)+psvin(sy,xp,xm,zb,4))
934 c elseif(iqq.eq.1)then
935 c om2p=psvy(xp,xprem,xm,xmrem,b,2)+psvy(xp,xprem,xm,xmrem,b,-2)
936 c * +psvx(xp,xprem,xm,xmrem,b,1)
937 c elseif(iqq.eq.2)then
938 c om2p=psvy(xp,xprem,xm,xmrem,b,3)
939 c * +psvx(xp,xprem,xm,xmrem,b,2)
940 c elseif(iqq.eq.-2)then
941 c om2p=psvy(xp,xprem,xm,xmrem,b,-3)
942 c * +psvx(xp,xprem,xm,xmrem,b,-2)
943 c elseif(iqq.eq.3)then
944 c om2p=psvy(xp,xprem,xm,xmrem,b,4)+psvy(xp,xprem,xm,xmrem,b,-4)
945 c * +psvx(xp,xprem,xm,xmrem,b,3)
946 c elseif(iqq.eq.4)then
947 c om2p=psvy(xp,xprem,xm,xmrem,b,5)+psvy(xp,xprem,xm,xmrem,b,7)
948 c * +psvy(xp,xprem,xm,xmrem,b,-5)+psvy(xp,xprem,xm,xmrem,b,-7)
949 c * +psvx(xp,xprem,xm,xmrem,b,4)+psvx(xp,xprem,xm,xmrem,b,-4)
950 c elseif(iqq.eq.5)then
951 c om2p=psvy(xp,xprem,xm,xmrem,b,6)+psvy(xp,xprem,xm,xmrem,b,-8)
952 c * +psvx(xp,xprem,xm,xmrem,b,5)
953 c elseif(iqq.eq.6)then
954 c om2p=psvy(xp,xprem,xm,xmrem,b,-6)+psvy(xp,xprem,xm,xmrem,b,8)
955 c * +psvx(xp,xprem,xm,xmrem,b,-5)
956 c elseif(iqq.eq.7)then
957 c om2p=psvy(xp,xprem,xm,xmrem,b,9)
958 c * +psvx(xp,xprem,xm,xmrem,b,6)
959 c elseif(iqq.eq.-7)then
960 c om2p=psvy(xp,xprem,xm,xmrem,b,-9)
961 c * +psvx(xp,xprem,xm,xmrem,b,-6)
963 c stop'om2p-wrong iqq!!!'
968 cc-----------------------------------------------------------------------
969 c double precision function om3p(xh,yp,xleg,xprem,xmrem,xlrem
971 cc-----------------------------------------------------------------------
972 cc om3p - chi~(x,y) for cut pomeron (nuclear effects)
973 cc xh - fraction of the energy squared s for the pomeron;
974 cc yp - rapidity for the pomeron;
975 cc xleg - x for the pomeron leg;
976 cc xprem - x+ for the projectile remnant;
977 cc xmrem - x- for the target remnant;
978 cc xlrem - x for the leg remnant;
979 cc b1 - impact parameter between the pomeron ends;
980 cc b2 - impact parameter for the second pomeron end;
989 cc iqq = 9 - uncut-h+,
990 cc iqq = 10 - uncut-h-,
991 cc iqq = 11 - uncut-YY+,
992 cc iqq = 12 - uncut-YY-,
993 cc-----------------------------------------------------------------------
994 c double precision xh,yp,xleg,xprem,xmrem,xlrem
997 c return !!!!!!!!!!!!!!!
998 cc if(iqq.ne.1.and.iqq.ne.5.and.iqq.ne.9.and.iqq.ne.10
999 cc *.and.iqq.ne.11.and.iqq.ne.12)return
1001 cc$$$ xp=dsqrt(xh)*exp(yp)
1002 cc$$$ if(xh.ne.0.d0)then
1011 cc-----------------------------------------------------------------------
1012 c double precision function om4p(xx1,xx2,xx3,xx4
1013 c *,b12,b13,b14,b23,b24,b34,iqq)
1014 cc-----------------------------------------------------------------------
1015 cc om4p - chi for 2-leg contributions
1016 cc xx_i - x+- for pomeron ends;
1017 cc b_ij - impact parameter diff. between pomeron ends;
1018 cc iqq = 1 - uncut-H,
1019 cc iqq = 2 - uncut-YY+,
1020 cc iqq = 3 - uncut-YY-
1021 cc-----------------------------------------------------------------------
1022 c double precision xx1,xx2xx3,xx4
1027 cc------------------------------------------------------------------------
1028 c function omi5pp(sy,xpp,xpm,z,iqq) !former psfsh1
1029 cc-----------------------------------------------------------------------
1030 cc omi5pp - integrated semihard interaction eikonal
1031 cc sy - energy squared for the hard interaction,
1032 cc z - impact parameter factor, z=exp(-b**2/rp),
1033 cc iqq - type of the hard interaction:
1034 cc 0 - soft, 1 - gg, 2 - qg, 3 - gq
1035 cc-----------------------------------------------------------------------
1036 c common /ar3/ x1(7),a1(7)
1037 c common /ar9/ x9(3),a9(3)
1038 c include 'epos.inc'
1039 c include 'epos.incsem'
1040 c fsy(zsy)=zsy**dels !*(1.-1./zsy)**betpom
1043 c if(iclpro.eq.4.and.iqq.eq.2.or.icltar.eq.4.and.iqq.eq.3)then
1044 c spmin=4.*q2min+2.*qcmass**2
1045 c elseif(iqq.ne.0)then
1050 c if(sy.le.spmin)return
1052 c rp=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy))
1053 c alpq=(alppar+1.)/2.
1063 c xpmax=(1.-spmin/sy)**(1.+alplea(iclp))
1066 c xp=1.-(xpmax*(.5+x9(i)*(m-1.5)))**(1./(1.+alplea(iclp)))
1067 c xmmax=(1.-spmin/sy/xp)**(1.+alplea(iclt))
1070 c xm=1.-(xmmax*(.5+x9(i1)*(m1-1.5)))**(1./(1.+alplea(iclt)))
1073 c rh=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy1))
1074 c omi5pp=omi5pp+a9(i)*a9(i1)*fsy(sy1)*xmmax*z**(rp/rh)/rh
1075 c * *(xp*xm)**(-alppar)
1080 c omi5pp=omi5pp*xpmax/(1.+alplea(iclp))/(1.+alplea(iclt))
1081 c * *chad(iclpro)*chad(icltar)*gamhad(iclpro)*gamhad(icltar)
1082 c * *(xpp*xpm)**(1.-alppar)/4.
1086 c xmin=(spmin/sy)**(delh-dels)
1089 c zh=(.5*(1.+xmin-(2*m-3)*x9(i)*(1.-xmin)))**(1./(delh-dels))
1090 c if(iclpro.eq.4.and.iqq.eq.2.or.icltar.eq.4.and.iqq.eq.3)then
1091 c call psjti0(zh*sy,sgq,sgqb,4,0)
1092 c call psjti0(zh*sy,sqq,sqqb,4,1)
1094 c call psjti0(zh*sy,sgg,sggb,0,0)
1095 c call psjti0(zh*sy,sgq,sgqb,0,1)
1096 c call psjti0(zh*sy,sqq,sqqb,1,1)
1097 c call psjti0(zh*sy,sqaq,sqaqb,-1,1)
1098 c call psjti0(zh*sy,sqqp,sqqpb,1,2)
1099 c sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
1106 c xx=.5+x9(i1)*(m1-1.5)
1110 c xp1max=(1.-xp)**(1.+alplea(iclp))
1111 c xm1max=(1.-xm)**(1.+alplea(iclt))
1114 c xp1=1.-(xp1max*(.5+x9(i2)*(m2-1.5)))
1115 c * **(1./(1.+alplea(iclp)))
1118 c xm1=1.-(xm1max*(.5+x9(i3)*(m3-1.5)))
1119 c * **(1./(1.+alplea(iclt)))
1120 c if(xp1.lt.xp.or.xm1.lt.xm)write (*,*)'xp1,xm1,xp,xm'
1123 c rh=r2had(iclpro)+r2had(icltar)+slopom
1124 c * *log(xp1*xm1/xp/xm)
1125 c glu1=(1.-xp/xp1)**betpom*(1.-glusea)
1126 c sea1=EsoftQZero(xp/xp1)*glusea
1127 c glu2=(1.-xm/xm1)**betpom*(1.-glusea)
1128 c sea2=EsoftQZero(xm/xm1)*glusea
1129 c stg=stg+a9(i1)*a9(i2)*a9(i3)*(glu1*glu2*sgg
1130 c * +(glu1*sea2+sea1*glu2)*sgq+sea1*sea2*sqq)
1131 c * *xp1max*xm1max*(xp1*xm1)**(dels-alppar)
1139 c omi5pp=omi5pp-a9(i)*log(zh)*stg/zh**delh
1143 c xpmin=zh**(dels+.5)
1146 c xp=(.5*(1.+xpmin-(2*m1-3)*x9(i1)*(1.-xpmin)))
1147 c * **(1./(dels+.5))
1149 c if(xp*xpp.lt..99999)then
1150 c uv1=psdfh4(xp*xpp,q2min,0.,iclp,1)
1151 c dv1=psdfh4(xp*xpp,q2min,0.,iclp,2)
1152 c xm1max=(1.-xm)**(1.+alplea(iclt))
1155 c xm1=1.-(xm1max*(.5+x9(i2)*(m2-1.5)))
1156 c * **(1./(1.+alplea(iclt)))
1158 c rh=r2had(iclpro)+r2had(icltar)+slopom*log(xm1/xm)
1159 c glu2=(1.-xm/xm1)**betpom*(1.-glusea)
1160 c sea2=EsoftQZero(xm/xm1)*glusea
1161 c stq=stq+a9(i1)*a9(i2)*(glu2*sgq+sea2*sqq)*(uv1+dv1)
1162 c * *z**(rp/rh)/rh*xm1max*xm1**(dels-alppar)/sqrt(xp)
1163 c * *((1.-xp)/(1.-xp*xpp))**(1.-alpq+alplea(iclp))
1169 c stq=stq*(1.-xpmin)
1170 c omi5pp=omi5pp+a9(i)*stq/zh**delh
1176 c omi5pp=omi5pp*(1.-xmin)/(delh-dels)
1178 c omi5pp=omi5pp*chad(iclp)*chad(iclt)*gamhad(iclp)
1179 c * *gamhad(iclt)*rr**2*(xpp*xpm)**(1.-alppar)
1180 c * /(1.+alplea(iclp))/(1.+alplea(iclt))*pi/8.*factk
1182 c omi5pp=omi5pp*chad(iclp)*chad(iclt)*rr*gamhad(iclt)
1183 c * *xpp**(1.-alpq)*xpm**(1.-alppar)/(.5+dels)
1184 c * /(1.+alplea(iclt))/16.*factk
1189 c------------------------------------------------------------------------
1190 function om52pi(sy,xpp,xpm,iqq,je1,je2) !modified om51pp
1191 c-----------------------------------------------------------------------
1192 c sy - energy squared for the hard interaction
1194 c iqq = 0 - sea-sea,
1195 c iqq = 1 - val-sea,
1196 c iqq = 2 - sea-val,
1197 c iqq = 3 - val-val,
1199 c je = emission type
1200 c 0 ... no emissions
1204 c already b-averaged (\int d2b /sigine*10)
1205 c-----------------------------------------------------------------------
1206 common /ar3/ x1(7),a1(7)
1207 common /psar7/ delx,alam3p,gam3p
1209 include 'epos.incsem'
1210 if(iqq.lt.0.or.iqq.gt.3)stop'om52pi: unvalid iqq'
1218 if( je1.ge.1 .and. je2.ge.1) ef1=1
1219 if( je1.ge.1 .and.(je2.eq.0.or.je2.eq.2))ef2=1
1220 if((je1.eq.0.or.je1.eq.2).and. je2.ge.1) ef3=1
1221 if((je1.eq.0.or.je1.eq.2).and.(je2.eq.0.or.je2.eq.2))ef4=1
1224 if(sy.le.spmin)goto999
1228 ctp060829 icls=icltar
1229 elseif(iqq.eq.2)then
1230 ctp060829 icls=iclpro
1235 if(iqq.eq.3)delss=-0.5
1237 xmin=xmin**(delh-delss)
1240 c numerical integration over zh
1243 zh=(.5*(1.+xmin-(2*m-3)*x1(i)*(1.-xmin)))**(1./(delh-delss))
1244 sgg= ef1 *pijet(2,q2min,q2min,zh*sy,0,0)
1245 * + (ef2+ef3)*pijet(1,q2min,q2min,zh*sy,0,0)
1246 * + ef4 *pijet(0,q2min,q2min,zh*sy,0,0)
1247 sgq= ef1 *pijet(2,q2min,q2min,zh*sy,0,1)
1248 * + (ef2+ef3)*pijet(1,q2min,q2min,zh*sy,0,1)
1249 * + ef4 *pijet(0,q2min,q2min,zh*sy,0,1)
1250 sqq= ef1 *pijet(2,q2min,q2min,zh*sy,1,1)
1251 * + (ef2+ef3)*pijet(1,q2min,q2min,zh*sy,1,1)
1252 * + ef4 *pijet(0,q2min,q2min,zh*sy,1,1)
1253 sqaq= ef1 *pijet(2,q2min,q2min,zh*sy,-1,1)
1254 * + (ef2+ef3)*pijet(1,q2min,q2min,zh*sy,-1,1)
1255 * + ef4 *pijet(0,q2min,q2min,zh*sy,-1,1)
1256 sqqp= ef1 *pijet(2,q2min,q2min,zh*sy,1,2)
1257 * + (ef2+ef3)*pijet(1,q2min,q2min,zh*sy,1,2)
1258 * + ef4 *pijet(0,q2min,q2min,zh*sy,1,2)
1260 sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
1265 xx=.5+x1(i1)*(m1-1.5)
1268 glu1=EsoftGluonTil(xp)
1269 sea1=EsoftQuarkTil(xp)
1270 glu2=EsoftGluonTil(xm)
1271 sea2=EsoftQuarkTil(xm)
1273 * +(glu1*sea2+sea1*glu2)*sgq !ccccc
1274 * +sea1*sea2*sqq !ccccc
1278 om52pi=om52pi-a1(i)*log(zh)*stg/zh**delh
1279 elseif(iqq.eq.3)then
1280 stq=0. !int^1_(sqrt(z)) dx_p / x_p / sqrt(1-x_p) =int^(tmax)_(0) dt
1281 tmax=sqrt(1.-sqrt(zh)) !t=ln((1+sqrt(1-x_p))/(1-sqrt(1-x_p)))
1282 tmax=log((1.+tmax)/(1.-tmax))
1283 if(tmax.gt.1.e-20)then
1286 t=tmax*(.5+x1(i1)*(m1-1.5))
1287 z01=((1.d0-exp(-1.d0*t))/(1.d0+exp(-1.d0*t)))**2
1290 if(xp*xpp.le..9999.and.xm*xpm.le..9999
1291 * .or.xm*xpp.le..9999.and.xp*xpm.le..9999)then
1293 * *(psharg(xp*xpp,xm*xpm,sqqi,sqqp,sqaq)
1294 * +psharg(xm*xpp,xp*xpm,sqqi,sqqp,sqaq))
1295 * *max(1e-20,1.-xp)**(.5-alpq)
1296 * *max(1e-20,1.-xm)**(-alpq)
1297 * *xp**delss*xm**delss
1298 * *xpp**alppar/gamhad(iclpro) ! Eval
1299 * *xpm**alppar/gamhad(icltar) ! Eval
1305 om52pi=om52pi+a1(i)*stq/zh**delh
1306 elseif(iqq.eq.1.or.iqq.eq.2)then
1311 t=tmax*(.5+x1(i1)*(m1-1.5))
1314 if(xp*xpp.lt..99999)then
1315 uv1=psdfh4(xp*xpp,q2min,0.,iclv,1) ! Eval
1316 dv1=psdfh4(xp*xpp,q2min,0.,iclv,2) ! Eval
1317 glu2=EsoftGluonTil(xm)
1318 sea2=EsoftQuarkTil(xm)
1321 * dstq=(glu2*sgq+sea2*sqq)*(uv1+dv1)
1322 * *(1.-xp*xpp)**(-1.+alpq-alplea(iclv)) ! Eval
1323 * *xp**(delss-.5)*(1.-xp)**(-alpq+.5) ! Eval *sqrt(1-x)/sqrt(x)
1324 * *xpp**alppar/gamhad(iclv) ! Eval
1330 om52pi=om52pi+a1(i)*stq/zh**delh
1332 stop'om52pi: unvalid iqq (2). '
1337 om52pi=om52pi*(1.-xmin)/(delh-delss)
1341 elseif(iqq.eq.3)then
1343 * * utgam1(2.+alplea(iclpro)-alpq) ! Eval
1344 * /utgam1(1.+alplea(iclpro))/utgam1(1.-alpq) ! Eval
1345 * * utgam1(2.+alplea(icltar)-alpq) ! Eval
1346 * /utgam1(1.+alplea(icltar))/utgam1(1.-alpq) ! Eval
1347 * /xpp**alpq/xpm**alpq ! Eval
1348 elseif(iqq.le.2)then
1350 * *utgam1(2.+alplea(iclv)-alpq)/utgam1(1.+alplea(iclv)) ! Eval
1351 * /utgam1(1.-alpq) ! Eval
1356 om52pi=om52pi*factk * .0390 /sigine*10 /2.
1359 c------------------------------------------------------------------------
1360 function psharg(zh1,zh2,sqq,sqqp,sqaq)
1361 c-----------------------------------------------------------------------
1362 include 'epos.incsem'
1366 if(zh1.le..9999.and.zh2.le..9999)then
1367 uv1=psdfh4(zh1,q2min,0.,iclpro,1)
1368 dv1=psdfh4(zh1,q2min,0.,iclpro,2)
1369 uv2=psdfh4(zh2,q2min,0.,icltar,1)
1370 dv2=psdfh4(zh2,q2min,0.,icltar,2)
1371 if(iclpro.eq.2.and.icltar.eq.2)then !proton
1372 fff=sqq*(uv1*uv2+dv1*dv2)+sqqp*(uv1*dv2+dv1*uv2)
1373 elseif(iclpro.eq.1.or.icltar.eq.1)then !pion
1374 fff=sqq*uv1*uv2+sqaq*dv1*dv2+sqqp*(uv1*dv2+dv1*uv2)
1375 elseif(iclpro.eq.3.or.icltar.eq.3)then !kaon
1376 fff=sqq*uv1*uv2+sqqp*(uv1*dv2+dv1*uv2+dv1*dv2)
1377 elseif(iclpro.eq.4.or.icltar.eq.4)then !J/psi
1378 fff=sqq*uv1*(uv2+dv2)
1381 * *(1.-zh1)**(-1.+alpq-alplea(iclpro))
1382 * *(1.-zh2)**(-1.+alpq-alplea(icltar))
1389 c------------------------------------------------------------------------
1390 function om51pp(sy,xpp,z,iqq) !former psfsh
1391 c-----------------------------------------------------------------------
1392 c om51pp - semihard interaction eikonal
1393 c sy - energy squared for the hard interaction,
1394 c z - impact parameter factor, z=exp(-b**2/rp),
1395 c iqq - type of the hard interaction:
1396 c 0 - gg, 1 - qg, 2 - gq, 3 - gg(int), 4 - gg(proj), 5 - qg(proj),
1397 c 6 - gg(int)|b=0, 7 - <b^2*gg(int)>, 8 - gg(proj)|b=0,
1398 c 9 - <b^2*gg(proj)>, 10 - qg(proj)|b=0, 11 - <b^2*qg(proj)>
1399 c-----------------------------------------------------------------------
1400 common /ar3/ x1(7),a1(7)
1401 common /psar7/ delx,alam3p,gam3p
1403 include 'epos.incsem'
1406 if(iqq.eq.0.or.iqq.eq.3.or.iqq.eq.4
1407 *.or.iqq.eq.6.or.iqq.eq.7.or.iqq.eq.8.or.iqq.eq.9
1408 *.or.iclpro.ne.4.and.(iqq.eq.1.or.iqq.eq.5
1409 *.or.iqq.eq.10.or.iqq.eq.11)
1410 *.or.icltar.ne.4.and.iqq.eq.2)then
1413 spmin=4.*q2min+2.*qcmass**2
1415 if(sy.le.spmin)goto999
1417 if(iqq.eq.1.or.iqq.eq.5.or.iqq.eq.10.or.iqq.eq.11)then
1420 elseif(iqq.eq.2)then
1426 xmin=xmin**(delh-dels)
1427 rp=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy))
1430 c numerical integration over zh
1433 zh=(.5*(1.+xmin-(2*m-3)*x1(i)*(1.-xmin)))**(1./
1435 if(iqq.eq.0.or.iqq.eq.3.or.iqq.eq.4
1436 * .or.iqq.eq.6.or.iqq.eq.7.or.iqq.eq.8.or.iqq.eq.9
1437 * .or.iclpro.ne.4.and.(iqq.eq.1.or.iqq.eq.5
1438 * .or.iqq.eq.10.or.iqq.eq.11)
1439 * .or.icltar.ne.4.and.iqq.eq.2)then
1440 call psjti0(zh*sy,sgg,sggb,0,0) !inclusive (sj) and born (sjb)
1441 call psjti0(zh*sy,sgq,sgqb,0,1)
1442 call psjti0(zh*sy,sqq,sqqb,1,1)
1443 call psjti0(zh*sy,sqaq,sqaqb,-1,1)
1444 call psjti0(zh*sy,sqqp,sqqpb,1,2)
1445 sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
1446 c...........test.......
1447 c tgg= psjet(q2min,q2min,q2min,zh*sy,0,0,0)
1448 c * +2*psjet1(q2min,q2min,q2min,zh*sy,0,0,0)
1449 c * + psborn(q2min,q2min,q2min,zh*sy,0,0,0,1)
1450 c tgq= psjet(q2min,q2min,q2min,zh*sy,0,1,0)
1451 c * +2*psjet1(q2min,q2min,q2min,zh*sy,0,1,0)
1452 c * + psborn(q2min,q2min,q2min,zh*sy,0,1,0,1)
1453 c tqq= psjet(q2min,q2min,q2min,zh*sy,1,1,0)
1454 c * +2*psjet1(q2min,q2min,q2min,zh*sy,1,1,0)
1455 c * + psborn(q2min,q2min,q2min,zh*sy,1,1,0,1)
1456 c tqa= psjet(q2min,q2min,q2min,zh*sy,-1,1,0)
1457 c * +2*psjet1(q2min,q2min,q2min,zh*sy,-1,1,0)
1458 c * + psborn(q2min,q2min,q2min,zh*sy,-1,1,0,1)
1459 c tqqp= psjet(q2min,q2min,q2min,zh*sy,1,2,0)
1460 c * +2*psjet1(q2min,q2min,q2min,zh*sy,1,2,0)
1461 c * + psborn(q2min,q2min,q2min,zh*sy,1,2,0,1)
1462 c write(6,'(f12.2,3x,2f7.3,2(3x,2f7.3))')
1463 c * zh*sy,tgg,sgg, tgq,sgq, tqqp,sqqp
1464 c.......................
1466 call psjti0(zh*sy,sgq,sgqb,4,0)
1467 call psjti0(zh*sy,sqq,sqqb,4,1)
1470 if(iqq.eq.0.or.iqq.eq.3.or.iqq.eq.4
1471 * .or.iqq.eq.6.or.iqq.eq.7.or.iqq.eq.8.or.iqq.eq.9)then
1475 xx=.5+x1(i1)*(m1-1.5)
1478 glu1=(1.-xp)**betpom*(1.-glusea)
1479 sea1=EsoftQZero(xp)*glusea
1480 glu2=(1.-xm)**betpom*(1.-glusea)
1481 sea2=EsoftQZero(xm)*glusea
1483 rh=r2had(iclpro)+r2had(icltar)-slopom*log(zh)
1484 elseif(iqq.eq.3.or.iqq.eq.4)then
1486 elseif(iqq.eq.6.or.iqq.eq.7)then
1487 rh=alam3p-slopom*log(zh)
1488 elseif(iqq.eq.8.or.iqq.eq.9)then
1489 rh=r2had(iclpro)+.5*alam3p-slopom*log(zh)
1491 dstg=(glu1*glu2*sgg+
1492 * (glu1*sea2+sea1*glu2)*sgq+sea1*sea2*sqq)
1494 if(iqq.eq.7.or.iqq.eq.9)dstg=dstg*rh**2
1498 om51pp=om51pp-a1(i)*log(zh)*stg/zh**delh
1504 t=tmax*(.5+x1(i1)*(m1-1.5))
1507 if(xp*xpp.lt..99999)then
1508 uv1=psdfh4(xp*xpp,q2min,0.,iclv,1)
1509 dv1=psdfh4(xp*xpp,q2min,0.,iclv,2)
1510 glu2=(1.-xm)**betpom*(1.-glusea)
1511 sea2=EsoftQZero(xm)*glusea
1513 rh=r2had(iclpro)+r2had(icltar)-slopom*log(xm)
1514 elseif(iqq.eq.5)then
1516 elseif(iqq.le.10.or.iqq.le.11)then
1517 rh=r2had(iclpro)+.5*alam3p-slopom*log(xm)
1521 * dstq=(glu2*sgq+sea2*sqq)*(uv1+dv1)
1523 * *(1.-xp*xpp)**(-1.+alpq-alplea(iclv))
1524 * *xp**(dels-.5)*(1.-xp)**(-alpq+.5)
1525 if(iqq.eq.11)dstq=dstq*rh**2
1531 om51pp=om51pp+a1(i)*stq/zh**delh
1536 om51pp=om51pp*(1.-xmin)/(delh-dels)/sy**delh/2.
1538 om51pp=om51pp*chad(iclpro)*chad(icltar)*gamhad(iclpro)
1539 * *gamhad(icltar)*rr**2*pi
1540 elseif(iqq.eq.3)then
1541 om51pp=om51pp*rr**2*pi*4.*.0389
1542 elseif(iqq.eq.6)then
1543 om51pp=om51pp*rr**2*pi
1544 elseif(iqq.eq.7)then
1545 om51pp=om51pp*rr**2*pi*(4.*.0389)**2
1546 elseif(iqq.eq.4.or.iqq.eq.8.or.iqq.eq.9)then
1547 om51pp=om51pp*rr**2*pi*chad(iclpro)*gamhad(iclpro)
1548 if(iqq.eq.4)om51pp=om51pp*4.*.0389
1549 if(iqq.eq.9)om51pp=om51pp*(4.*.0389)**2
1550 elseif(iqq.le.2)then
1551 om51pp=om51pp*chad(iclpro)*chad(icltar)*rr*gamhad(icls)
1552 * *utgam1(2.+alplea(iclv)-alpq)/utgam1(1.+alplea(iclv))
1553 * /utgam1(1.-alpq)/2./xpp**alpq
1554 elseif(iqq.eq.5.or.iqq.eq.10.or.iqq.eq.11)then
1555 om51pp=om51pp*chad(iclv)*rr
1556 * *utgam1(2.+alplea(iclv)-alpq)/utgam1(1.+alplea(iclv))
1557 * /utgam1(1.-alpq)/2./xpp**alpq
1558 if(iqq.eq.5)om51pp=om51pp*4.*.0389
1559 if(iqq.eq.11)om51pp=om51pp*(4.*.0389)**2
1564 c------------------------------------------------------------------------
1565 subroutine psfz(gz2,b)
1566 c-----------------------------------------------------------------------
1567 c hadron-nucleus cross sections calculation
1568 c b - impact parameter squared
1569 c-----------------------------------------------------------------------
1570 double precision PhiExpo
1572 common /ar3/ x1(7),a1(7)
1573 external pttcs,pprcs
1577 rs=r2had(iclpro)+r2had(icltar)+max(slopom,slopoms)*log(engy**2)
1578 & +gwidth*(r2had(iclpro)+r2had(icltar))
1579 & +bmxdif(iclpro,icltar)/4./0.0389
1588 b1=sqrt(-rpom*log(zv1))
1589 b2=sqrt(-rpom*log(zv2))
1591 vv21=sngl(PhiExpo(1.,1.d0,1.d0,engy**2,b1))
1592 vv22=sngl(PhiExpo(1.,1.d0,1.d0,engy**2,b2))
1594 if(maproj.eq.1.and.matarg.eq.1)then
1597 elseif(matarg.eq.1)then
1598 cg1=ptrot(pprcs,b,b1)
1599 cg2=ptrot(pprcs,b,b2)
1601 cg1=ptrot(pttcs,b,b1)
1602 cg2=ptrot(pttcs,b,b2)
1605 gz2=gz2+a1(i1)*(cg1*(1.-vv21)+cg2*(1.-vv22)/z)
1614 c------------------------------------------------------------------------
1615 function ptgau(func,bm,iqq)
1616 c-----------------------------------------------------------------------
1617 c impact parameter integration for impact parameters <bm -
1618 c for nucleus-nucleus and hadron-nucleus cross-sections calculation
1619 c iqq=1 : projectile, iqq=2 : target
1620 c-----------------------------------------------------------------------
1622 common /ar3/ x1(7),a1(7)
1628 b=bm*sqrt(.5+x1(i)*(m-1.5))
1629 ptgau=ptgau+func(b,iqq)*a1(i)
1632 ptgau=ptgau*bm**2*pi*.5
1636 c------------------------------------------------------------------------
1637 function ptgau1(bm,iqq)
1638 c-----------------------------------------------------------------------
1639 c impact parameter integration for impact parameters >bm -
1640 c for hadron-nucleus cross-sections calculation
1641 c iqq=1 : projectile, iqq=2 : target
1642 c-----------------------------------------------------------------------
1644 common /ar5/ x5(2),a5(2)
1654 ptgau1=ptgau1+ptfau(b,iqq)*a5(i)*exp(x5(i))*b*2.*pi*difn
1658 c------------------------------------------------------------------------
1660 c-----------------------------------------------------------------------
1661 c impact parameter integration for impact parameters >bm -
1662 c for nucleus-nucleus cross-sections calculation
1663 c-----------------------------------------------------------------------
1665 common /ar5/ x5(2),a5(2)
1668 difn=difnuc(maproj)+difnuc(matarg)
1671 ptgau2=ptgau2+ptfauAA(b)*a5(i)*exp(x5(i))*b*2.*pi*difn
1677 c------------------------------------------------------------------------
1678 function ptfau(b,iqq)
1679 c-----------------------------------------------------------------------
1680 c ptfau - integrands for hadron-nucleus cross-sections calculation
1681 c iqq=1 : projectile, iqq=2 : target
1682 c-----------------------------------------------------------------------
1684 common /psar35/ anorm,anormp
1689 ptfau=1.-max(0.,(1.-anormp*gz2))**maproj
1691 ptfau=1.-max(0.,(1.-anorm*gz2))**matarg
1697 c------------------------------------------------------------------------
1699 c-----------------------------------------------------------------------
1700 c ptfau - integrands for hadron-nucleus cross-sections calculation
1701 c-----------------------------------------------------------------------
1703 common /ar3/ x1(7),a1(7)
1704 common /psar35/ anorm,anormp
1709 rs=r2had(iclpro)+r2had(icltar)+max(slopom,slopoms)*log(engy**2)
1710 & +gwidth*(r2had(iclpro)+r2had(icltar))
1711 & +bmxdif(iclpro,icltar)/4./0.0389
1718 b1=sqrt(-rpom*log(zv1))
1719 b2=sqrt(-rpom*log(zv2))
1722 ptfau1=max(0.,(1.-anorm*gz21))**matarg
1723 ptfau2=max(0.,(1.-anorm*gz22))**matarg
1724 cg1=ptrot(pprcs,b,b1)
1725 cg2=ptrot(pprcs,b,b2)
1726 ptfauAA=ptfauAA+a1(i1)*(cg1*(1.-ptfau1)+cg2*(1.-ptfau2)/z)
1729 ptfauAA=ptfauAA*rpom/2.
1730 ptfauAA=1.-max(0.,(1.-anormp*ptfauAA))**maproj
1735 c------------------------------------------------------------------------
1736 function ptrot(func,s,b)
1737 c-----------------------------------------------------------------------
1738 c convolution of nuclear profile functions (axial angle integration)
1739 c-----------------------------------------------------------------------
1740 common /ar8/ x2(4),a2
1745 sb1=b**2+s**2-2.*b*s*(2.*x2(i)-1.)
1746 sb2=b**2+s**2-2.*b*s*(1.-2.*x2(i))
1747 ptrot=ptrot+(func(sb1)+func(sb2))
1753 c------------------------------------------------------------------------
1755 c-----------------------------------------------------------------------
1756 c ptt - nuclear profile function value at imp param squared b*difnuc**2
1757 c-----------------------------------------------------------------------
1759 common /psar34/ rrr,rrrm
1760 common /ar5/ x5(2),a5(2)
1761 common /ar9/ x9(3),a9(3)
1763 b=b0/difnuc(matarg)**2
1773 z1=zm*(1.+x9(i))*0.5
1774 z2=zm*(1.-x9(i))*0.5
1775 quq=sqrt(b+z1**2)-rrr
1776 if (quq.lt.85.)pttcs=pttcs+a9(i)/(1.+exp(quq))
1777 quq=sqrt(b+z2**2)-rrr
1778 if (quq.lt.85.)pttcs=pttcs+a9(i)/(1.+exp(quq))
1785 quq=sqrt(b+z1**2)-rrr-x5(i)
1786 if (quq.lt.85.)dt=dt+a5(i)/(exp(-x5(i))+exp(quq))
1793 c------------------------------------------------------------------------
1795 c-----------------------------------------------------------------------
1796 c ptt - nuclear profile function value at imp param squared b*difnuc**2
1797 c-----------------------------------------------------------------------
1799 common /psar41/ rrrp,rrrmp
1800 common /ar5/ x5(2),a5(2)
1801 common /ar9/ x9(3),a9(3)
1803 b=b0/difnuc(maproj)**2
1813 z1=zm*(1.+x9(i))*0.5
1814 z2=zm*(1.-x9(i))*0.5
1815 quq=sqrt(b+z1**2)-rrrp
1816 if (quq.lt.85.)pprcs=pprcs+a9(i)/(1.+exp(quq))
1817 quq=sqrt(b+z2**2)-rrrp
1818 if (quq.lt.85.)pprcs=pprcs+a9(i)/(1.+exp(quq))
1825 quq=sqrt(b+z1**2)-rrrp-x5(i)
1826 if (quq.lt.85.)dt=dt+a5(i)/(exp(-x5(i))+exp(quq))
1833 c------------------------------------------------------------------------------
1834 function pscrse(ek,mapr,matg)
1835 c------------------------------------------------------------------------------
1836 c hadron-nucleus (hadron-proton) and nucl-nucl particle production cross section
1837 c ek - lab kinetic energy for the interaction
1838 c maproj - projec mass number (1<maproj<64)
1839 c matarg - target mass number (1<matarg<64)
1840 c------------------------------------------------------------------------------
1841 dimension wk(3),wa(3),wb(3)
1843 common /psar33/ asect(7,4,7),asectn(7,7,7)
1844 common /psar34/ rrr,rrrm
1845 common /psar35/ anorm,anormp
1846 common /psar41/ rrrp,rrrmp
1847 external ptfau,ptfauAA
1850 call idmass(1120,amt1)
1851 call idmass(1220,amt2)
1852 amtar=0.5*(amt1+amt2)
1853 if(matg.eq.1)amtar=amt1
1855 call idmass(idproj,ampro)
1860 c p=sqrt(max(0.,egy**2-ampro**2))
1861 egy=sqrt( 2*egy*amtar+amtar**2+ampro**2 )
1870 if(matg.eq.1.and.mapr.eq.1)then
1873 elseif(mapr.eq.1)then
1876 rrr=rad/difnuc(matg)
1878 anorm=1.5/pi/rrr**3/(1.+(pi/rrr)**2)/difnuc(matg)**2
1879 gin=(ptgau(ptfau,bm,2)+ptgau1(bm,2))*10. !sig_in
1880 elseif(matg.eq.1)then
1883 rrrp=rad/difnuc(mapr)
1885 anormp=1.5/pi/rrrp**3/(1.+(pi/rrrp)**2)/difnuc(mapr)**2
1886 gin=(ptgau(ptfau,bm,1)+ptgau1(bm,1))*10. !sig_in
1889 radp=radnuc(mapr)+1.
1891 rrr=rad/difnuc(matg)
1893 rrrp=radp/difnuc(mapr)
1895 anorm=1.5/pi/rrr**3/(1.+(pi/rrr)**2)/difnuc(matg)**2
1896 anormp=1.5/pi/rrrp**3/(1.+(pi/rrrp)**2)/difnuc(mapr)**2
1897 gin=(ptgau(ptfauAA,bm,2)+ptgau2(bm))*10.
1904 ye=log10(max(1.,egy/1.5))+1.
1908 wk(3)=wk(2)*(wk(2)-1.)*.5
1909 wk(1)=1.-wk(2)+wk(3)
1910 wk(2)=wk(2)-2.*wk(3)
1913 ya=log(ya)/.69315+1.
1916 wa(3)=wa(2)*(wa(2)-1.)*.5
1917 wa(1)=1.-wa(2)+wa(3)
1918 wa(2)=wa(2)-2.*wa(3)
1924 pscrse=pscrse+asect(je+i-1,iclpro,ja+m-1)*wk(i)*wa(m)
1931 yb=log(yb)/.69315+1.
1934 wb(3)=wb(2)*(wb(2)-1.)*.5
1935 wb(1)=1.-wb(2)+wb(3)
1936 wb(2)=wb(2)-2.*wb(3)
1941 pscrse=pscrse+asectn(je+i-1,jb+n-1,ja+m-1)*wk(i)*wa(m)*wb(n)
1953 c------------------------------------------------------------------------------
1954 function eposcrse(ek,mapro,matar,id)
1955 c------------------------------------------------------------------------------
1956 c inelastic cross section of epos
1957 c (id=0 corresponds to air)
1958 c ek - kinetic energy for the interaction
1959 c maproj - projec mass number (1<maproj<64)
1960 c matarg - target mass number (1<matarg<64)
1961 c------------------------------------------------------------------------------
1968 eposcrse=eposcrse+airwnxs(k)*pscrse(ek,mapro,mt)
1971 eposcrse=pscrse(ek,mapro,matar)
1978 cc------------------------------------------------------------------------
1979 c function pshard1(sy,xpp,xpm,z)
1980 cc-----------------------------------------------------------------------
1981 cc pshard - qq-pomeron eikonal
1982 cc sy - energy squared for the pomeron,
1983 cc xpp - lc+ for the pomeron,
1984 cc xpm - lc- for the pomeron
1985 cc-----------------------------------------------------------------------
1986 c common /ar3/ x1(7),a1(7)
1987 c common /ar9/ x9(3),a9(3)
1988 c include 'epos.inc'
1989 c include 'epos.incsem'
1992 c if(iclpro.ne.4.and.icltar.ne.4)then
1995 c spmin=4.*q2min+2.*qcmass**2
1997 c if(sy.le.spmin)return
1999 c rp=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy))
2000 c alpq=(alppar+1.)/2.
2001 c xmin=spmin/sy !min hard pomeron mass share
2002 c xminl=xmin**(delh+.5)
2006 c zh=(.5*(1.+xminl-(2*m-3)*x9(i)*(1.-xminl)))**(1./(delh+.5))
2007 c if(iclpro.ne.4.and.icltar.ne.4)then
2008 c call psjti0(zh*sy,sqq,sqqb,1,1)
2009 c call psjti0(zh*sy,sqqp,sqqpb,1,2)
2010 c call psjti0(zh*sy,sqaq,sqaqb,-1,1)
2012 c call psjti0(zh*sy,sqq,sqqb,4,1)
2020 c xx=.5+x9(i1)*(m1-1.5)
2023 c if(xp*xpp.le..9999.and.xm*xpm.le..9999.or.
2024 c * xm*xpp.le..9999.and.xp*xpm.le..9999)then
2025 c stq=stq+a9(i1)*psharf(xp*xpp,xm*xpm,sqq,sqqp,sqaq)
2026 c * *(1.-xp)**(1.+alplea(iclpro)-alpq)
2027 c * *(1.-xm)**(1.+alplea(icltar)-alpq)
2031 c pshard1=pshard1-a9(i)*stq/zh**(delh+0.5)*log(zh)
2034 c pshard1=pshard1*(1.-xminl)/(delh+.5)/4.*factk
2035 c **chad(iclpro)*chad(icltar)*(xpp*xpm)**(1.-alpq)
2036 c **z**(rp/(r2had(iclpro)+r2had(icltar)))
2037 c */(8.*pi*(r2had(iclpro)+r2had(icltar)))
2041 c------------------------------------------------------------------------
2042 function pshard(sy,xpp,xpm)
2043 c-----------------------------------------------------------------------
2044 c pshard - qq-pomeron eikonal
2045 c sy - energy squared for the pomeron,
2046 c xpp - lc+ for the pomeron,
2047 c xpm - lc- for the pomeron
2048 c-----------------------------------------------------------------------
2049 double precision z01
2050 common /ar3/ x1(7),a1(7)
2052 include 'epos.incsem'
2055 if(iclpro.ne.4.and.icltar.ne.4)then
2058 spmin=4.*q2min+2.*qcmass**2
2060 if(sy.le.spmin)return
2063 xmin=spmin/sy !min hard pomeron mass share
2064 xminl=xmin**(delh+.5)
2068 zh=(.5*(1.+xminl-(2*m-3)*x1(i)*(1.-xminl)))**(1./(delh+.5))
2069 if(iclpro.ne.4.and.icltar.ne.4)then
2070 call psjti0(zh*sy,sqq,sqqb,1,1)
2071 call psjti0(zh*sy,sqqp,sqqpb,1,2)
2072 call psjti0(zh*sy,sqaq,sqaqb,-1,1)
2074 call psjti0(zh*sy,sqq,sqqb,4,1)
2079 stq=0. !int^1_(sqrt(z)) dx_p / x_p / sqrt(1-x_p) =int^(tmax)_(0) dt
2080 tmax=sqrt(1.-sqrt(zh)) !t=ln((1+sqrt(1-x_p))/(1-sqrt(1-x_p)))
2081 tmax=log((1.+tmax)/(1.-tmax))
2082 if(tmax.gt.1.e-20)then
2085 t=tmax*(.5+x1(i1)*(m1-1.5))
2086 z01=((1.d0-exp(-1.d0*t))/(1.d0+exp(-1.d0*t)))**2
2089 if(xp*xpp.le..9999.and.xm*xpm.le..9999.or.
2090 * xm*xpp.le..9999.and.xp*xpm.le..9999)then
2091 stq=stq+a1(i1)*(psharf(xp*xpp,xm*xpm,sqq,sqqp,sqaq)+
2092 * psharf(xm*xpp,xp*xpm,sqq,sqqp,sqaq))
2093 * *z01**(.5-alpq)/(1.-xm)**alpq
2099 pshard=pshard+a1(i)*stq/zh**(delh+0.5)
2102 pshard=pshard*(1.-xminl)/(delh+.5)/4.*
2103 *utgam1(2.+alplea(iclpro)-alpq)/utgam1(1.+alplea(iclpro))/
2105 *utgam1(2.+alplea(icltar)-alpq)/utgam1(1.+alplea(icltar))/
2107 *chad(iclpro)*chad(icltar)/(8.*pi*(r2had(iclpro)+r2had(icltar)))*
2108 *(xpp*xpm)**(-alpq)/sy**delh
2112 c------------------------------------------------------------------------
2113 function psharf(zh1,zh2,sqq,sqqp,sqaq)
2114 c-----------------------------------------------------------------------
2115 include 'epos.incsem'
2119 if(zh1.le..9999.and.zh2.le..9999)then
2120 uv1=psdfh4(zh1,q2min,0.,iclpro,1)
2121 dv1=psdfh4(zh1,q2min,0.,iclpro,2)
2122 uv2=psdfh4(zh2,q2min,0.,icltar,1)
2123 dv2=psdfh4(zh2,q2min,0.,icltar,2)
2124 if(iclpro.eq.2.and.icltar.eq.2)then !proton
2125 fff=sqq*(uv1*uv2+dv1*dv2)+sqqp*(uv1*dv2+dv1*uv2)
2126 elseif(iclpro.eq.1.or.icltar.eq.1)then !pion
2127 fff=sqq*uv1*uv2+sqaq*dv1*dv2+sqqp*(uv1*dv2+dv1*uv2)
2128 elseif(iclpro.eq.3.or.icltar.eq.3)then !kaon
2129 fff=sqq*uv1*uv2+sqqp*(uv1*dv2+dv1*uv2+dv1*dv2)
2130 elseif(iclpro.eq.4.or.icltar.eq.4)then !J/psi
2131 fff=sqq*uv1*(uv2+dv2)
2133 psharf=fff*(1.-zh1)**(-1.+alpq-alplea(iclpro))*
2134 * (1.-zh2)**(-1.+alpq-alplea(icltar))
2141 c------------------------------------------------------------------------
2142 function psvin(sy,xpp,xpm,z,iqq)
2143 c-----------------------------------------------------------------------
2144 c psvin - contributions to the interaction eikonal
2145 c sy - energy squared for the hard interaction,
2146 c xpp - lc+ for the sh pomeron,
2147 c xpm - lc- for the sh pomeron,
2148 c z - impact parameter factor, z=exp(-b**2/4*rp),
2153 c iqq = 5 - gg(int),
2154 c iqq = 6 - gg(proj),
2155 c iqq = 7 - qg(proj),
2156 c iqq = 9 - total uncut-integrated,
2157 c iqq = 10 - total cut,
2158 c iqq = 14 - gg(int)|b=0,
2159 c iqq = 15 - <b^2*gg(int)>,
2160 c iqq = 16 - gg(proj)|b=0,
2161 c iqq = 17 - <b^2*gg(proj)>,
2162 c iqq = 18 - qg(proj)|b=0,
2163 c iqq = 19 - <b^2*qg(proj)>
2164 c-----------------------------------------------------------------------
2165 dimension wk(3),wi(3),wj(3),wz(3),fa(3)
2166 common /psar2/ edmax,epmax
2167 common /psar4/ fhgg(11,10,8),fhqg(11,10,80)
2168 *,fhgq(11,10,80),fhqq(11,10,80),fhgg0(11,10),fhgg1(11,10,4)
2169 *,fhqg1(11,10,40),fhgg01(11),fhgg02(11),fhgg11(11,4)
2170 *,fhgg12(11,4),fhqg11(11,10,4),fhqg12(11,10,4)
2171 *,ftoint(11,14,2,2,3)
2172 common /psar7/ delx,alam3p,gam3p
2174 include 'epos.incsem'
2187 rp=r2had(iclpro)+r2had(icltar)+slopom*log(max(1.,sy))
2190 if(iqq.eq.1.or.iqq.eq.5.or.iqq.eq.6.or.iqq.eq.14
2191 *.or.iqq.eq.15.or.iqq.eq.16.or.iqq.eq.17
2192 *.or.iclpro.ne.4.and.(iqq.eq.2.or.iqq.eq.7
2193 *.or.iqq.eq.18.or.iqq.eq.19)
2194 *.or.icltar.ne.4.and.iqq.eq.3
2195 *.or.iclpro.ne.4.and.icltar.ne.4)then
2198 spmin=4.*q2min+2.*qcmass**2
2200 if(sy.le.spmin.and.(iqq.le.7.or.iqq.gt.13))return
2202 if(iqq.le.7.or.iqq.gt.13)then
2203 yl=log(sy/spmin)/log(epmax/2./spmin)*10.+1
2207 wk(3)=wk(2)*(wk(2)-1.)*.5
2208 wk(1)=1.-wk(2)+wk(3)
2209 wk(2)=wk(2)-2.*wk(3)
2211 if(iqq.ne.4)then !---------------- not 4 ------------------
2215 psvin=max(0.,exp(fhgg01(k+1))*wk(2)
2216 * +exp(fhgg01(k+2))*wk(3))
2218 psvin=exp(fhgg01(k)*wk(1)+fhgg01(k+1)*wk(2)
2219 * +fhgg01(k+2)*wk(3))
2221 psvin=psvin*factk*sy**delh
2224 elseif(iqq.eq.15)then
2226 psvin=max(0.,exp(fhgg02(k+1))*wk(2)
2227 * +exp(fhgg02(k+2))*wk(3))
2229 psvin=exp(fhgg02(k)*wk(1)+fhgg02(k+1)*wk(2)
2230 * +fhgg02(k+2)*wk(3))
2232 psvin=psvin*factk*sy**delh
2235 elseif(iqq.eq.6)then
2237 psvin=max(0.,exp(fhgg11(k+1,iclpro))*wk(2)
2238 * +exp(fhgg11(k+2,iclpro))*wk(3))
2240 psvin=exp(fhgg11(k,iclpro)*wk(1)+fhgg11(k+1,iclpro)*wk(2)
2241 * +fhgg11(k+2,iclpro)*wk(3))
2243 psvin=psvin*factk*sy**delh*xp**(-alppar)
2246 elseif(iqq.eq.17)then
2248 psvin=max(0.,exp(fhgg12(k+1,iclpro))*wk(2)
2249 * +exp(fhgg12(k+2,iclpro))*wk(3))
2251 psvin=exp(fhgg12(k,iclpro)*wk(1)+fhgg12(k+1,iclpro)*wk(2)
2252 * +fhgg12(k+2,iclpro)*wk(3))
2254 psvin=psvin*factk*sy**delh*xp**(-alppar)
2257 elseif(iqq.eq.7.or.iqq.eq.19)then
2259 xl=log(10.*xp)/log(2.)+5.
2268 wi(3)=wi(2)*(wi(2)-1.)*.5
2269 wi(1)=1.-wi(2)+wi(3)
2270 wi(2)=wi(2)-2.*wi(3)
2276 fhhh=fhqg11(k2,i+i1-1,iclpro)
2277 elseif(iqq.eq.19)then
2278 fhhh=fhqg12(k2,i+i1-1,iclpro)
2280 fa(k1)=fa(k1)+fhhh*wi(i1)
2284 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2286 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2288 psvin=psvin*factk*sy**delh
2296 wz(3)=wz(2)*(wz(2)-1.)*.5
2297 wz(1)=1.-wz(2)+wz(3)
2298 wz(2)=wz(2)-2.*wz(3)
2303 fa(k1)=fhgg0(k2,jz)*wz(1)+fhgg0(k2,jz+1)
2304 * *wz(2)+fhgg0(k2,jz+2)*wz(3)
2307 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2309 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2311 psvin=psvin*z*factk*sy**delh
2313 elseif(iqq.eq.16)then
2316 fa(k1)=fhgg1(k2,jz,iclpro)*wz(1)+fhgg1(k2,jz+1,iclpro)
2317 * *wz(2)+fhgg1(k2,jz+2,iclpro)*wz(3)
2320 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2322 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2324 psvin=psvin*z*factk*sy**delh*xp**(-alppar)
2326 elseif(iqq.eq.18)then
2328 xl=log(10.*xp)/log(2.)+5.
2337 wi(3)=wi(2)*(wi(2)-1.)*.5
2338 wi(1)=1.-wi(2)+wi(3)
2339 wi(2)=wi(2)-2.*wi(3)
2345 l2=jz+l1-1+10*(iclpro-1)
2346 fhhh=fhqg1(k2,i+i1-1,l2)
2347 fa(k1)=fa(k1)+fhhh*wi(i1)*wz(l1)
2352 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2354 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2356 psvin=psvin*z*factk*sy**delh
2358 elseif(iqq.eq.1)then !1111111111111111111111111111111111
2362 iclpt=iclpro+4*(icltar-1)
2363 fa(k1)=fhgg(k2,jz,iclpt)*wz(1)+fhgg(k2,jz+1,iclpt)
2364 * *wz(2)+fhgg(k2,jz+2,iclpt)*wz(3)
2367 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2369 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2371 psvin=psvin*z*factk*sy**delh*(xp*xm)**(-alppar)
2373 else ! 2222222222222222222222 3333333333333333333333 ....
2376 xl=log(10.*xp)/log(2.)+5.
2385 wi(3)=wi(2)*(wi(2)-1.)*.5
2386 wi(1)=1.-wi(2)+wi(3)
2387 wi(2)=wi(2)-2.*wi(3)
2394 l2=jz+l1-1+10*(iclpro+4*(icltar-1)-1)
2395 fhhh=fhqg(k2,i+i1-1,l2)
2396 elseif(iqq.eq.3)then
2397 l2=jz+l1-1+10*(iclpro+4*(icltar-1)-1)
2398 fhhh=fhgq(k2,i+i1-1,l2)
2400 fa(k1)=fa(k1)+fhhh*wi(i1)*wz(l1)
2405 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2407 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2409 psvin=psvin*xm**(-alppar)*z*factk*sy**delh
2412 else ! ------------- 4444444444444444444 -----------------------
2415 xl1=log(10.*xp)/log(2.)+5.
2423 wi(3)=wi(2)*(wi(2)-1.)*.5
2424 wi(1)=1.-wi(2)+wi(3)
2425 wi(2)=wi(2)-2.*wi(3)
2428 xl2=log(10.*xm)/log(2.)+5.
2436 wj(3)=wj(2)*(wj(2)-1.)*.5
2437 wj(1)=1.-wj(2)+wj(3)
2438 wj(2)=wj(2)-2.*wj(3)
2445 j2=j+j1-1+10*(iclp+4*(iclt-1)-1)
2446 fa(k1)=fa(k1)+fhqq(k2,i+i1-1,j2)*wi(i1)*wj(j1)
2451 psvin=max(0.,exp(fa(2))*wk(2)+exp(fa(3))*wk(3))
2453 psvin=exp(fa(1)*wk(1)+fa(2)*wk(2)+fa(3)*wk(3))
2455 psvin=psvin*z**(rp/(r2had(iclpro)+r2had(icltar)))*
2458 endif !--------------------------------------------
2463 yl=log(sy)/log(1.e8)*10.+1
2465 k=min(k,9) !?????????????9
2467 wk(3)=wk(2)*(wk(2)-1.)*.5
2468 wk(1)=1.-wk(2)+wk(3)
2469 wk(2)=wk(2)-2.*wk(3)
2480 wz(3)=wz(2)*(wz(2)-1.)*.5
2481 wz(1)=1.-wz(2)+wz(3)
2482 wz(2)=wz(2)-2.*wz(3)
2489 psvin=psvin+ftoint(k2,l2,icdp,icdt,iclp)*wk(k1)*wz(l1)
2498 c------------------------------------------------------------------------
2499 function psbint(q1,q2,qqcut,ss,m1,l1,jdis)
2500 c-----------------------------------------------------------------------
2501 c psbint - born cross-section interpolation
2502 c q1 - virtuality cutoff at current end of the ladder;
2503 c q2 - virtuality cutoff at opposite end of the ladder;
2504 c qqcut - p_t cutoff for the born process;
2505 c s - total c.m. energy squared for the scattering,
2506 c m1 - parton type at current end of the ladder (0 - g, 1,-1,2,... - q)
2507 c l1 - parton type at opposite end of the ladder (0 - g, 1,-1,2,... - q)
2508 c-----------------------------------------------------------------------
2509 dimension wi(3),wk(3)
2510 common /psar2/ edmax,epmax
2511 common /psar21/ csbor(20,160,2)
2512 include 'epos.incsem'
2513 double precision psuds
2522 if(iabs(m1).ne.4)then
2524 if(m1.ne.0.and.m1.eq.l1)then
2527 elseif(m1.ne.0.and.m1.eq.-l1)then
2530 elseif(m1.ne.0.and.l1.ne.0.and.m1.ne.l1)then
2543 spmin=4.*q2min+q2mass
2545 if(s.le.s2min)return
2549 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
2556 ml=20*(m-1)+80*(l-1)
2557 qli=log(qq/q2min)/log(qmax/q2min)*19.+1.
2558 sl=log(s/spmin)/log(epmax/2./spmin)*19.+1.
2567 wi(3)=wi(2)*(wi(2)-1.)*.5
2568 wi(1)=1.-wi(2)+wi(3)
2569 wi(2)=wi(2)-2.*wi(3)
2572 wk(3)=wk(2)*(wk(2)-1.)*.5
2573 wk(1)=1.-wk(2)+wk(3)
2574 wk(2)=wk(2)-2.*wk(3)
2578 psbint=psbint+csbor(i+i1-1,k+k1+ml-1,jdis+1)
2582 psbint=exp(psbint)*(1./tmin-1./tmax)
2583 if(jdis.eq.0.and.qq.gt.q1)then
2584 psbint=psbint*sngl(psuds(qq,m1)/psuds(q1,m1))
2585 elseif(jdis.eq.1.and.4.*qq.gt.q1)then
2586 psbint=psbint*sngl(psuds(4.*qq,m1)/psuds(q1,m1))
2588 if(qq.gt.q2)psbint=psbint*sngl(psuds(qq,l1)/psuds(q2,l1))
2592 c-----------------------------------------------------------------------
2593 function psborn(q1,q2,qqcut,s,j,l,jdis,md)
2594 c-----------------------------------------------------------------------
2596 c hard 2->2 parton scattering born cross-section
2597 c including sudakov on both sides
2599 c q1 - virtuality cutoff at current end of the ladder;
2600 c q2 - virtuality cutoff at opposite end of the ladder;
2601 c qqcut - p_t cutoff for the born process;
2602 c s - c.m. energy squared for the scattering;
2603 c j - parton type at current end of the ladder (0 - g, 1,2 etc. - q);
2604 c l - parton type at opposite end of the ladder (0 - g, 1,2 etc. - q).
2605 c-----------------------------------------------------------------------
2606 common /ar3/ x1(7),a1(7)
2607 double precision sud0,psbornd,psuds
2609 include 'epos.incsem'
2619 c if(j.ne.3)then !kkkkkkkkkk charm is 3 ???
2629 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
2632 ! return !tmin=2.*qq !kkkkkkk !????????????? tp why not ?
2635 sud0=psuds(q1,j1)*psuds(q2,l)
2640 t=2.*tmin/(1.+tmin/tmax-x1(i)*(2*m-3)
2643 if(qt.lt..999*qq.and.ish.ge.1)write(ifch,*)'psborn:qt,qq,q1,q2'
2651 if(j1.eq.0.and.l.eq.0)then
2652 fb=ffborn(s,t, 1. , 0. , 0. , 0. , 0. ) !gg
2653 elseif(j1*l.eq.0)then
2654 fb=ffborn(s,t, 0. , 1. , 0. , 0. , 0.) !qg
2656 fb=ffborn(s,t, 0. , 0. , 1. , 0. , 0.) !qq
2657 elseif(j1.eq.-l)then
2658 fb=ffborn(s,t, 0. , 0. , 0. , 1. , 0.) !qq
2660 fb=ffborn(s,t, 0. , 0. , 0. , 0. , 1.) !qq
2662 fb=fb*pssalf(qt/qcdlam)**2
2663 psbornd=psbornd+dble(a1(i)*fb)*dble(t)**2
2664 & *psuds(scale,j1)*psuds(qt,l)
2667 psbornd=psbornd*dble(2.*pi**3)/dble(s)**2/sud0*2
2668 * /2 !CS for parton pair
2669 if(md.eq.1)psbornd=psbornd*(1./tmin-1./tmax)
2670 psborn=sngl(psbornd)
2674 c------------------------------------------------------------------------
2675 function psdgh(s,qq,long)
2676 c-----------------------------------------------------------------------
2678 c s - energy squared for the interaction (hadron-hadron),
2679 c-----------------------------------------------------------------------
2680 common/ar3/ x1(7),a1(7)
2681 common /cnsta/ pi,pii,hquer,prom,piom,ainfin
2682 include 'epos.incsem'
2683 double precision psuds
2687 psdgh=(psdfh4(xd,q2min,0.,2,1)/2.25+psdfh4(xd,q2min,0.,2,2)/9.
2688 * +psdfh4(xd,q2min,0.,2,3)/9.+
2689 * 2.*(psdfh4(xd,q2min,0.,2,-1)+psdfh4(xd,q2min,0.,2,-2)+
2690 * psdfh4(xd,q2min,0.,2,-3))/4.5)
2691 * *sngl(psuds(qq,1)/psuds(q2min,1))*4.*pi**2*alfe/qq
2698 s2min=qq/(1.-q2ini/qq)
2700 s2min=4.*max(q2min,qcmass**2)+qq
2701 s2min=s2min/(1.-4.*q2ini/(s2min-qq))
2706 do i=1,7 !numerical integration over z1
2709 z1=qq/s+(xmin-qq/s)*((1.-qq/s)/(xmin-qq/s))
2710 * **(.5+(m-1.5)*x1(i))
2712 z1=.5*(1.+xmin+(2*m-3)*x1(i)*(1.-xmin))
2714 call psdint(z1*s,qq,sds,sdn,sdb,sdt,sdr,1,long)
2715 call psdint(z1*s,qq,sdsg,sdng,sdbg,sdtg,sdrg,0,long)
2716 tu=psdfh4(z1,q2min,0.,2,1)
2717 td=psdfh4(z1,q2min,0.,2,2)
2718 ts=psdfh4(z1,q2min,0.,2,3)
2719 tg=psdfh4(z1,q2min,0.,2,0)
2720 tsea=2.*(psdfh4(z1,q2min,0.,2,-1)+psdfh4(z1,q2min,0.,2,-2)
2721 * +psdfh4(z1,q2min,0.,2,-3))
2722 gy=sdn*(tu/2.25+td/9.+ts/9.+tsea/4.5)+sdtg*tg/4.5
2723 * +sdt*(tu+td+ts+tsea)/4.5
2724 dgh=dgh+a1(i)*gy*(1.-qq/s/z1)
2727 dgh=dgh*log((1.-qq/s)/(xmin-qq/s))*.5
2733 c------------------------------------------------------------------------
2734 function psdh(s,qq,iclpro0,long)
2735 c-----------------------------------------------------------------------
2736 c pshard - hard quark-quark interaction cross-section
2737 c s - energy squared for the interaction (hadron-hadron),
2738 c iclpro0 - type of the primary hadron (nucleon)
2739 c-----------------------------------------------------------------------
2740 common /ar3/ x1(7),a1(7)
2741 include 'epos.incsem'
2743 double precision psuds
2747 if(long.eq.0.and.(idisco.eq.0.or.idisco.eq.1))then
2748 psdh=(psdfh4(xd,qqs,0.,iclpro0,1)/2.25+
2749 * psdfh4(xd,qqs,0.,iclpro0,2)/9.)
2750 * *sngl(psuds(qq,1)/psuds(qqs,1))
2758 s2min=qq/(1.-q2ini/qq)
2760 s2min=4.*max(q2min,qcmass**2)+qq
2761 s2min=s2min/(1.-4.*q2ini/(s2min-qq))
2765 do i=1,7 !numerical integration over z1
2768 z1=qq/s+(xmin-qq/s)*((1.-qq/s)/(xmin-qq/s))
2769 * **(.5+(m-1.5)*x1(i))
2771 z1=.5*(1.+xmin+(2*m-3)*x1(i)*(1.-xmin))
2773 call psdint(z1*s,qq,sds,sdn,sdb,sdt,sdr,1,long)
2774 tu=psdfh4(z1,qqs,0.,iclpro0,1)
2775 td=psdfh4(z1,qqs,0.,iclpro0,2)
2776 gy=sdt*(tu+td)/4.5+sdn*(tu/2.25+td/9.)
2786 dh=dh*log((1.-qq/s)/(xmin-qq/s))*.5
2795 c------------------------------------------------------------------------
2796 function psdsh(s,qq,iclpro0,dqsh,long)
2797 c-----------------------------------------------------------------------
2798 c psdsh - semihard interaction eikonal
2799 c s - energy squared for the interaction (hadron-hadron),
2800 c iclpro0 - hadron class,
2801 c z - impact parameter factor, z=exp(-b**2/rp),
2802 c iqq - type of the hard interaction (0 - gg, 1 - qg, 2 - gq)
2803 c-----------------------------------------------------------------------
2804 common /ar3/ x1(7),a1(7)
2806 include 'epos.incsem'
2807 double precision psuds
2810 if(long.eq.0.and.(idisco.eq.0.or.idisco.eq.1))then
2811 dqsh=fzeroSeaZZ(xd,iclpro0)/xd**dels
2812 * *rr*4.*pi*gamhad(iclpro0)/
2813 * 4.5*sngl(psuds(qq,1)/psuds(q2min,1))
2820 s2min=qq/(1.-q2ini/qq)
2822 s2min=qq+4.*max(q2min,qcmass**2)
2825 xmin=xmin**(delh-dels)
2828 c numerical integration over z1
2831 z1=(.5*(1.+xmin-(2*m-3)*x1(i)*(1.-xmin)))**(1./
2833 call psdint(z1*s,qq,sdsg,sdng,sdbg,sdtg,sdrg,0,long)
2834 call psdint(z1*s,qq,sdsq,sdnq,sdbq,sdtq,sdrq,1,long)
2835 dsh=dsh+a1(i)/z1**delh*(sdtg*fzeroGluZZ(z1,iclpro0)
2836 * +(sdtq+sdnq)*fzeroSeaZZ(z1,iclpro0))
2839 dsh=dsh*(1.-xmin)/(delh-dels)/2.
2841 psdsh=dqsh+dsh*rr*4.*pi*gamhad(iclpro0)/4.5 !*ccorr(1,1,iclpro0)
2845 c------------------------------------------------------------------------
2846 function psdsh1(s,qq,iclpro0,dqsh,long)
2847 c-----------------------------------------------------------------------
2848 c psdsh - semihard interaction eikonal
2849 c s - energy squared for the interaction (hadron-hadron),
2850 c iclpro0 - hadron class,
2851 c z - impact parameter factor, z=exp(-b**2/rp),
2852 c iqq - type of the hard interaction (0 - gg, 1 - qg, 2 - gq)
2853 c-----------------------------------------------------------------------
2854 common /ar3/ x1(7),a1(7)
2856 include 'epos.incsem'
2857 c double precision psuds
2859 psdsh1=0. !only for plotting in psaevp : not use any more
2862 c$$$ write(ifch,*)'Psdsh1 for xd,qq',xd,qq
2863 c$$$ if(long.eq.0.and.(idisco.eq.0.or.idisco.eq.1))then
2864 c$$$ dqsh=psftist(xd)/4.5*sngl(psuds(qq,1)/psuds(q2min,1))
2865 c$$$ * *4.*pi**2*alfe/qq
2870 c$$$ if(long.eq.0)then
2871 c$$$ s2min=qq/(1.-q2ini/qq)
2873 c$$$ s2min=qq+4.*max(q2min,qcmass**2)
2876 c$$$ xmin=xmin**(delh-dels)
2878 c$$$ if(xmin.lt.1.)then
2879 c$$$c numerical integration over z1
2882 c$$$ z1=(.5*(1.+xmin-(2*m-3)*x1(i)*(1.-xmin)))**(1./
2884 c$$$ call psdint(z1*s,qq,sdsg,sdng,sdbg,sdtg,sdrg,0,long)
2885 c$$$ call psdint(z1*s,qq,sdsq,sdnq,sdbq,sdtq,sdrq,1,long)
2886 c$$$ dsh=dsh+a1(i)/z1**delh*(sdtg*psftigt(z1)
2887 c$$$ * +(sdtq+sdnq)*psftist(z1))*z1**dels
2890 c$$$ dsh=dsh*(1.-xmin)/(delh-dels)/2.
2892 c$$$ psdsh1=dqsh+dsh/4.5
2897 c------------------------------------------------------------------------
2898 function psev0(q1,qq,xx,j)
2899 c-----------------------------------------------------------------------
2900 double precision xx,psuds,psev00
2901 common /ar3/ x1(7),a1(7)
2902 include 'epos.incsem'
2908 if(j.eq.1)then !g->q
2909 qi=2.*q1/(1.+q1/qq+(1.-q1/qq)*(2.*m-3.)*x1(i))
2910 psev00=psev00+a1(i)*qi*psuds(qi,0)/psuds(qi,1)
2911 * /log(qi*(1.d0-xx)/qcdlam)
2913 qi=(.5*(q1+qq+(q1-qq)*(2.*m-3.)*x1(i)))
2914 psev00=psev00+a1(i)/qi/psuds(qi,0)*psuds(qi,1)
2915 * /log(qi*(1.d0-xx)/qcdlam)
2921 psev00=psev00*(1.d0/q1-1.d0/qq)*psuds(qq,1)/psuds(qq,0)/2.d0
2923 psev00=psev00*(qq-q1)*psuds(qq,0)/psuds(qq,1)/2.d0
2925 psev00=psev00/log(log(qq*(1.d0-xx)/qcdlam)
2926 & /log(q1*(1.d0-xx)/qcdlam))
2931 c------------------------------------------------------------------------
2932 function psev(q1,qq,xx,j,l,n)
2933 c------------------------------------------------------------------------
2934 double precision xx,zmax,zmax1,zmin,zmin1,z,psuds,fk,fq
2936 common /ar3/ x1(7),a1(7)
2937 include 'epos.incsem'
2945 if(zmin.lt.zmax)then
2946 if(zmin.lt..1d0)then
2947 zmax1=min(.1d0,zmax)
2951 z=xx+(zmin-xx)*((zmax1-xx)/(zmin-xx))**(.5+(m-1.5)*x1(i))
2953 z=zmin*(zmax1/zmin)**(.5+(m-1.5)*x1(i))
2955 z=(.5d0*(zmax1+zmin+(zmax1-zmin)*(2*m-3)*x1(i)))
2957 qmin=max(q2ini/(1.d0-xx/z),q2ini/(1.d0-z))
2965 qi=qmin*(qmax/qmin)**(.5+x1(i1)*(m1-1.5))
2967 qi=(.5*(qmax+qmin+(qmax-qmin)*(2.*m1-3.)*x1(i1)))
2970 if(j.eq.3.and.k.eq.1)then
2974 fk=dble(psevi0(q1,qi,xx/z,min(2,j),k))
2976 fk=dble(psevi(q1,qi,xx/z,j,k)/qi)
2980 fq=fq+a1(i1)*fk/psuds(qi,l-1)*pssalf(qt/qcdlam)
2984 fq=fq*log(qmax/qmin)*(1.d0-xx/z)
2990 fz1=fz1+a1(i)*fq*psfap(z,k-1,l-1)
2995 fz1=fz1*log((zmax1-xx)/(zmin-xx))/4.
2997 fz1=fz1*log(zmax1/zmin)/4.
2999 fz1=fz1*(zmax1-zmin)/4.
3003 if(zmax.gt..1d0)then
3004 zmin1=max(.1d0,zmin)
3007 z=1.d0-(1.d0-zmax)*((1.d0-zmin1)/(1.d0-zmax))**
3008 * (.5+x1(i)*(m-1.5))
3009 qmin=max(q2ini/(1.d0-z),q2ini/(1.d0-xx/z))
3017 qi=qmin*(qmax/qmin)**(.5+x1(i1)*(m1-1.5))
3019 qi=(.5*(qmax+qmin+(qmax-qmin)*(2.*m1-3.)*x1(i1)))
3022 if(j.eq.3.and.k.eq.1)then
3026 fk=dble(psevi0(q1,qi,xx/z,min(2,j),k))
3028 fk=dble(psevi(q1,qi,xx/z,j,k)/qi)
3032 fq=fq+a1(i1)*fk/psuds(qi,l-1)*pssalf(qt/qcdlam)
3036 fq=fq*log(qmax/qmin)
3040 fz2=fz2+a1(i)*fq*psfap(z,k-1,l-1)*(1.d0/z-1.d0)
3044 fz2=fz2*log((1.d0-zmin1)/(1.d0-zmax))/4.
3047 psev=sngl((fz1+fz2)*psuds(qq,l-1))
3051 c------------------------------------------------------------------------
3052 function psevi0(q1,qq,xx,m,l)
3053 c------------------------------------------------------------------------
3054 double precision xx,xmax,psuds
3055 dimension wi(3),wj(3),wk(3)
3056 common /psar2/ edmax,epmax
3057 common /psar31/ evk0(21,21,54)
3059 include 'epos.incsem'
3061 xmax=1.d0-2.d0*q2ini/epmax
3062 qmin=max(1.d0*q2min,q2ini/(1.d0-xx))
3064 if(qq.gt..5001*epmax.and.ish.ge.1)then
3065 write(ifch,*)'0-extrap.:q1,qq,epmax,xx,m,l:',q1,qq,epmax,xx,m,l
3068 if(xx.ge.xmax.or.qq.le.1.000*qm1)then
3070 c write (*,*)'xx,xmax,qq,qm1,qmin,q1',xx,xmax,qq,qm1,qmin,q1
3078 yx=log(10.d0*xx)+13.
3082 elseif(xx.lt..9d0)then
3087 yx=log(10.d0*(1.d0-xx))/log(10.d0*(1.d0-xmax))*6.+21
3092 wk(3)=wk(2)*(wk(2)-1.)*.5
3093 wk(1)=1.-wk(2)+wk(3)
3094 wk(2)=wk(2)-2.*wk(3)
3096 qli=log(qq/qmin)/log(.5*epmax/qmin)*20.+1.
3097 qlj=log(qm1/qmin)/log(qq/qmin)*20.+1.
3102 wi(3)=wi(2)*(wi(2)-1.)*.5
3103 wi(1)=1.-wi(2)+wi(3)
3104 wi(2)=wi(2)-2.*wi(3)
3110 wj(3)=wj(2)*(wj(2)-1.)*.5
3111 wj(1)=1.-wj(2)+wj(3)
3112 wj(2)=wj(2)-2.*wj(3)
3118 psevi0=psevi0+evk0(i+i1-1,j+j1-1,k+k1-1+27*(m-1))
3119 * *wi(i1)*wj(j1)*wk(k1)
3125 psevi0=psevi0*psfap(xx,m-1,l-1)*log(log(qq*(1.d0-xx)/qcdlam)
3126 */log(qm1*(1.d0-xx)/qcdlam))*sngl(psuds(qq,m-1)/psuds(q1,m-1))/4.5
3130 c------------------------------------------------------------------------
3131 function psevi(q1,qq,xx,m,l)
3132 c------------------------------------------------------------------------
3133 c m l: 1 1 ... gluon -> gluon
3134 c 2 1 ... quark -> gluon
3135 c 1 2 ... gluon -> quark
3136 c 3 2 ... quark -> quark non singlet
3137 c 2 2 ... quark -> quark all
3138 c singlet = all - non singlet
3139 c-----------------------------------------------------------------------
3140 double precision xx,xmax,psuds
3141 dimension wi(3),wj(3),wk(3)
3142 common /psar2/ edmax,epmax
3143 common /psar32/ evk(21,21,135)
3145 include 'epos.incsem'
3148 xmax=1.d0-2.d0*q2ini/epmax
3149 if(qq.gt..5001*epmax.and.ish.ge.1)then
3150 write(ifch,*)'1-extrap.:q1,qq,epmax,xx,m,l:',q1,qq,epmax,xx,m,l
3153 qmin=max(1.d0*q2min,q2ini/(1.d0-xx))
3155 if(xx.ge.xmax.or.qq.le.1.0001*qm1)then
3158 qmin1=max(1.d0*qmin,q2ini/(1.d0-dsqrt(xx)))
3159 if(qq.le.1.0001*qmin1)then
3160 psevi=psevi0(q1,qq,xx,min(m,2),l)
3165 yx=log(10.d0*xx)+13.
3169 elseif(xx.lt..9d0)then
3174 yx=log(10.d0*(1.d0-xx))/log(10.d0*(1.d0-xmax))*6.+21
3179 wk(3)=wk(2)*(wk(2)-1.)*.5
3180 wk(1)=1.-wk(2)+wk(3)
3181 wk(2)=wk(2)-2.*wk(3)
3183 qli=log(qq/qmin)/log(.5*epmax/qmin)*20.+1.
3184 qlj=log(qm1/qmin)/log(qq/qmin)*20.+1.
3189 wi(3)=wi(2)*(wi(2)-1.)*.5
3190 wi(1)=1.-wi(2)+wi(3)
3191 wi(2)=wi(2)-2.*wi(3)
3197 wj(3)=wj(2)*(wj(2)-1.)*.5
3198 wj(1)=1.-wj(2)+wj(3)
3199 wj(2)=wj(2)-2.*wj(3)
3207 k2=k+k1-1+27*(m-1)+54*(l-1)
3209 psevi=psevi+evk(i+i1-1,j+j1-1,k2)
3210 * *wi(i1)*wj(j1)*wk(k1)
3214 psevi=exp(psevi)*psfap(xx,m-1,l-1)*log(log(qq*(1.d0-xx)/qcdlam)
3215 */log(qm1*(1.d0-xx)/qcdlam))/4.5
3216 if(q1.lt.qm1)psevi=psevi*sngl(psuds(qm1,m-1)/psuds(q1,m-1))
3220 c------------------------------------------------------------------------
3221 function psjci(q1,s,l1)
3222 c-----------------------------------------------------------------------
3223 c psjci - inclusive ordered ladder cross-section interpolation for c-quark
3224 c q1 - virtuality cutoff at current end of the ladder
3225 c s - total c.m. energy squared for the ladder,
3226 c l1 - parton type at current end of the ladder (0-g, 1,2,etc.-q)
3227 c-----------------------------------------------------------------------
3228 dimension wi(3),wk(3)
3229 common /psar2/ edmax,epmax
3230 common /psar23/ cschar(20,20,2)
3231 include 'epos.incsem'
3235 spmin=4.*q2min+q2mass
3238 if(s.le.s2min)return
3240 smins=s2min/(1.-q2ini/q1)
3241 c if(s.le.smins)goto 1
3242 if(s.le.smins.or.qq.le.q2min)goto 1 !??????? ctp070618
3246 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
3254 qli=log(qq/q2min)/log(qmax/q2min)*19.+1.
3255 sl=log(s/spmin)/log(epmax/2./spmin)*19.+1.
3263 wi(3)=wi(2)*(wi(2)-1.)*.5
3264 wi(1)=1.-wi(2)+wi(3)
3265 wi(2)=wi(2)-2.*wi(3)
3268 wk(3)=wk(2)*(wk(2)-1.)*.5
3269 wk(1)=1.-wk(2)+wk(3)
3270 wk(2)=wk(2)-2.*wk(3)
3274 psjci=psjci+cschar(i+i1-1,k+k1-1,l)*wi(i1)*wk(k1)
3277 psjci=exp(psjci)*(1./tmin-1./tmax)
3279 1 psjci=psbint(q2min,q1,0.,s,4,l1,0)
3283 c-----------------------------------------------------------------------
3285 c-----------------------------------------------------------------------
3286 c psjct - unordered ladder cross-section for c-quark
3287 c s - c.m. energy squared for the scattering;
3288 c l - parton type at opposite end of the ladder (0 - g, 1,2 etc. - q).
3289 c-----------------------------------------------------------------------
3290 double precision xx,zmax,qmax,qmin,qi,zmin,fsj,z,s2,sj
3291 common /ar3/ x1(7),a1(7)
3293 include 'epos.incsem'
3297 zmax=dble(s)/(dble(s)+dble(5.*q2mass))
3298 qmax=zmax**2*dble(q2mass)/(1.d0-zmax)
3301 if(qmax.lt.qmin.and.ish.ge.1)write(ifch,*)'psjct:qmin,qmax'
3305 qi=2.d0*qmin/(1.d0+qmin/qmax+dble((2*m-3)*x1(i))
3306 * *(1.d0-qmin/qmax))
3307 zmax=(2.d0/(1.d0+dsqrt(1.d0+4.d0*dble(q2mass)/qi)))**delh
3308 zmin=(5.d0*qi/dble(s))**delh
3311 if(zmax.lt.zmin.and.ish.ge.1)write(ifch,*)'psjct:zmin,zmax'
3315 z=(.5d0*(zmax+zmin+dble((2*m1-3)*x1(i1))
3316 * *(zmax-zmin)))**(1./delh)
3319 sj=dble(psjti(sngl(qi),q2min,sngl(s2),0,l,0)*psfap(xx,1,0))*z
3320 fsj=fsj+dble(a1(i1))*sj*dble(pssalf(sngl(qi)/qcdlam))/z**delh
3324 psjct=psjct+a1(i)*sngl(fsj*qi)
3327 psjct=psjct*sngl(1./qmin-1./qmax)/delh/4.
3331 c------------------------------------------------------------------------
3332 function psjet1(q1,q2,qqcut,s,j,l,jdis)
3333 c-----------------------------------------------------------------------
3334 c psjet1 - ordered parton ladder cross-section
3335 c q1 - virtuality cutoff at current end of the ladder;
3336 c q2 - virtuality cutoff at opposite end of the ladder;
3337 c qqcut - p_t cutoff for the born process;
3338 c s - c.m. energy squared for the scattering;
3339 c j - parton type at current end of the ladder (0 - g, 1,2 etc. - q);
3340 c l - parton type at opposite end of the ladder (0 - g, 1,2 etc. - q).
3341 c-----------------------------------------------------------------------
3342 double precision xx,z,qq,xmax,xmin,s2min,smin,p1,q2ms,q2inis,xmin1
3343 *,sh,qtmin,t,xmax1,fx1,fx2,psuds
3344 common /ar3/ x1(7),a1(7)
3345 common /ar9/ x9(3),a9(3)
3347 include 'epos.incsem'
3352 elseif(jdis.eq.1)then
3353 qq=dble(max(q1/4.,q2))
3355 qq=dble(max(q1,q2/4.))
3357 qq=max(qq,dble(qqcut))
3363 s2min=dble(q2mass)+4.d0*qq
3364 if(jdis.eq.0.or.jdis.eq.2)then
3365 smin=s2min/(1.d0-dble(q2ini)/qq)
3367 smin=s2min/(1.d0-dble(q2ini)/qq/4.d0)
3369 if(dble(s).le.smin)return
3371 q2ms=dble(q2mass)/dble(s)
3372 q2inis=dble(q2ini)/dble(s)
3373 p1=dble(s)/(1.d0+q2ms)
3376 if(jdis.eq.0.or.jdis.eq.2)then
3377 xmax=.5d0*(1.d0+q2ms)+dsqrt(.25d0*(1.d0-q2ms)**2-4.d0*q2inis)
3379 xmax=.5d0*(1.+q2ms)+dsqrt(.25d0*(1.-q2ms)**2-q2inis)
3381 xmin=max(1.d0+q2ms-xmax,s2min/dble(s))
3382 if(xmin.ge.xmax.and.ish.ge.1)then
3383 write(ifch,*)'jti1,xmin,xmax',xmin,xmax
3389 if(xmax.gt..8d0)then
3390 xmin1=max(xmin,.8d0)
3393 z=1.d0-(1.d0-xmax)*((1.d0-xmin1)/(1.d0-xmax))**
3394 * (.5d0+dble(x9(i)*(m-1.5)))
3397 p1=sh/(1.d0+dble(q2mass)/sh)
3399 if(jdis.eq.0.or.jdis.eq.2)then
3400 qtmin=max(qq,dble(q2ini)/(1.d0-z))
3402 qtmin=max(qq,dble(q2ini)/(1.d0-z)/4.d0)
3404 tmin=2.d0*dble(qtmin)/(1.d0+dsqrt(1.d0-4.d0*dble(qtmin)/p1))
3408 if(tmin.ge.tmax.and.ish.ge.1)write(ifch,*)'psjet1:tmin,tmax'
3412 t=2.d0*tmin/(1.d0+tmin/tmax-dble(x9(i1)*(2*m1-3))
3413 & *(1.d0-tmin/tmax))
3414 qt=sngl(t*(1.d0-t/p1))
3415 c if(qt.lt.qtmin)write (*,*)'psjet1:qt,qq',qt,qq
3420 elseif(jdis.eq.1)then
3423 elseif(jdis.eq.2)then
3429 fb=fb+psjetj(q1,scale1,sngl(t),xx,sngl(sh),j,l,n)
3431 ft=ft+a9(i1)*fb*pssalf(qt/qcdlam)**2*sngl(t**2
3435 fx1=fx1+dble(a9(i)*ft)*(1.d0/tmin-1.d0/tmax)/sh**2*(1.d0-z)
3438 fx1=fx1*dlog((1.d0-xmin1)/(1.d0-xmax))
3441 if(xmin.lt..8d0)then
3442 xmax1=min(xmax,.8d0)**(-delh)
3446 z=(.5d0*(xmax1+xmin1+(xmin1-xmax1)*dble((2*m-3)*x9(i))))
3450 p1=sh/(1.d0+dble(q2mass)/sh)
3452 if(jdis.eq.0.or.jdis.eq.2)then
3453 qtmin=max(qq,dble(q2ini)/(1.d0-z))
3455 qtmin=max(qq,dble(q2ini)/(1.d0-z)/4.d0)
3457 tmin=2.d0*dble(qtmin)/(1.d0+dsqrt(1.d0-4.d0*dble(qtmin)/p1))
3461 if(tmin.ge.tmax.and.ish.ge.1)write(ifch,*)'psjet1:tmin,tmax'
3465 t=2.d0*tmin/(1.d0+tmin/tmax-dble(x9(i1)*(2*m1-3))
3466 & *(1.d0-tmin/tmax))
3467 qt=sngl(t*(1.d0-t/p1))
3468 if(qt.lt.sngl(qtmin).and.ish.ge.1)write(ifch,*)'psjet1:qt,qq'
3474 elseif(jdis.eq.1)then
3477 elseif(jdis.eq.2)then
3483 fb=fb+psjetj(q1,scale1,sngl(t),xx,sngl(sh),j,l,n)
3485 ft=ft+a9(i1)*fb*pssalf(qt/qcdlam)**2*sngl(t**2
3489 fx2=fx2+dble(a9(i)*ft)*(1.d0/tmin-1.d0/tmax)/sh**2*z**(1.+delh)
3492 fx2=fx2*(xmin1-xmax1)/dble(delh)
3494 psjet1=sngl((fx1+fx2)/psuds(q2,l))*pi**3*2
3495 * /2 !CS for parton pair
3499 c-----------------------------------------------------------------------
3500 function psjet(q1,q2,qqcut,s,j,l,jdis)
3501 c-----------------------------------------------------------------------
3502 c parton ladder cross-section
3503 c with at least one emission on each side
3505 c q1 - virtuality cutoff at current end of the ladder;
3506 c q2 - virtuality cutoff at opposite end of the ladder;
3507 c qqcut - p_t cutoff for the born process;
3508 c s - c.m. energy squared for the scattering;
3509 c j - parton type at current end of the ladder (0 - g, 1,2 etc. - q);
3510 c l - parton type at opposite end of the ladder (0 - g, 1,2 etc. - q).
3511 c-----------------------------------------------------------------------
3512 double precision xx1,xx2,qq,s2min,xmin,xmax,xmin1,xmax1,t,tmin
3513 *,tmax,sh,z,qtmin,ft,fx1,fx2
3514 common /ar3/ x1(7),a1(7)
3515 common /ar9/ x9(3),a9(3)
3517 include 'epos.incsem'
3518 common/ccctest/iiitest
3525 qq=dble(max(q1/4.,q2))
3527 qq=max(qq,dble(qqcut))
3529 if(dble(s).le.s2min/(1.d0-dble(q2ini)/qq)**2)return !kkkkkkk
3531 phi=acos(1.-54.*q2ini/s)/3.
3532 zmax=(1.+2.*cos(phi))**2/9. !kkkkkkk
3533 zmin=(1.-cos(phi)+sqrt(3.d0)*sin(phi))/3. !kkkkkkk
3534 zmin=max(zmin**2,sngl(s2min/dble(s)))
3535 if(zmin.gt.zmax.and.ish.ge.1)write(ifch,*)'psjet:zmin,zmax'
3541 z=dble(.5*(zmax+zmin+(zmin-zmax)*(2*m-3)*x9(i)))**(-1./delh)
3545 qtmin=max(qq,dble(q2ini)/(1.d0-dsqrt(z)))
3546 tmin=max(0.d0,1.d0-4.d0*qtmin/sh)
3547 tmin=2.d0*qtmin/(1.d0+dsqrt(tmin)) !kkkkkkk
3551 c if(tmin.gt.tmax)write (*,*)'psjet:tmin,tmax',tmin,tmax
3554 t=2.d0*tmin/(1.d0+tmin/tmax-dble(x9(i1)*(2*m1-3))
3555 & *(1.d0-tmin/tmax))
3557 c if(qt.lt.qtmin)write (*,*)'psjet:qt,qq',qt,qq
3559 xmin=max(dsqrt(z),z/xmax) !xm>xp !!!
3560 if(xmin.gt.xmax.and.ish.ge.1)write(ifch,*)'psjet:xmin,xmax'
3564 if(xmax.gt..8d0)then
3565 xmin1=max(xmin,.8d0)
3568 xx1=1.d0-(1.d0-xmax)*((1.d0-xmin1)/(1.d0-xmax))**
3569 * dble(.5+x9(i2)*(m2-1.5))
3573 fb=fb+psjeti(q1,q2,qt,sngl(t),xx1,xx2,sngl(sh)
3575 * +psjeti(q1,q2,qt,sngl(t),xx2,xx1,sngl(sh)
3577 fx1=fx1+dble(a9(i2)*fb)*(1.d0/xx1-1.d0)
3578 * *pssalf(qt/qcdlam)**2
3581 fx1=fx1*dlog((1.d0-xmin1)/(1.d0-xmax))
3583 if(xmin.lt..8d0)then
3584 xmax1=min(xmax,.8d0)
3587 xx1=xmin*(xmax1/xmin)**dble(.5+x9(i2)*(m2-1.5))
3591 fb=fb+psjeti(q1,q2,qt,sngl(t),xx1,xx2,sngl(sh)
3593 * +psjeti(q1,q2,qt,sngl(t),xx2,xx1,sngl(sh)
3595 fx2=fx2+dble(a9(i2))*fb*pssalf(qt/qcdlam)**2
3598 fx2=fx2*dlog(xmax1/xmin)
3600 ft=ft+dble(a9(i1))*(fx1+fx2)*t**2
3603 ft=ft*(1.d0/tmin-1.d0/tmax)
3604 psjet=psjet+a9(i)*sngl(ft*z**(1.+delh)/sh**2)
3607 psjet=psjet*(zmin-zmax)/delh*pi**3
3608 * /2. !CS for parton pair
3612 c-----------------------------------------------------------------------
3613 function pijet(ii,qi,qq,sk,m1,l1) !polynomial interpol of jet CS
3614 c-----------------------------------------------------------------------
3615 c ii ..... type of CS (2 = bothside, 1 = oneside, 0 = no emission, Born)
3616 c qi ..... virtuality cutoff at current end of the ladder
3617 c qq ..... virtuality cutoff of Born
3618 c sk ..... energy squared for the scattering
3619 c m1,l1 .. parton types
3620 c-----------------------------------------------------------------------
3621 include 'epos.incsem'
3622 common/psar2/edmax,epmax
3623 common/tabcsjet/ksmax,iqmax,jqmax,csjet(0:2,2,20,20,20,3,2)
3624 real wi(3),wj(3),wk(3)
3625 common/cpijet/npijet
3628 if(npijet.eq.1)call MakeCSTable
3630 if(m1.ne.0.and.m1.eq.l1)then
3633 elseif(m1.ne.0.and.m1.eq.-l1)then
3636 elseif(m1.ne.0.and.l1.ne.0.and.m1.ne.l1)then
3647 spmed=spmin*(epmax/2./spmin)**(1./(ksmax-1.))
3656 qli=1.+log(qi/q2min)/log(qmax/q2min)*(iqmax-1)
3657 qlj=1.+log(qq/qqmin)/log(qmax/qqmin)*(jqmax-1)
3658 sl= 1.+log(sk/spmin)/log(spmax/spmin)*(ksmax-1)
3665 if(k.gt.(ksmax-2))k=ksmax-2
3666 if(i.gt.(iqmax-2))i=iqmax-2
3667 if(j.gt.(jqmax-2))j=jqmax-2
3670 wi(3)=wi(2)*(wi(2)-1.)*.5
3671 wi(1)=1.-wi(2)+wi(3)
3672 wi(2)=wi(2)-2.*wi(3)
3675 wj(3)=wj(2)*(wj(2)-1.)*.5
3676 wj(1)=1.-wj(2)+wj(3)
3677 wj(2)=wj(2)-2.*wj(3)
3680 wk(3)=wk(2)*(wk(2)-1.)*.5
3681 wk(1)=1.-wk(2)+wk(3)
3682 wk(2)=wk(2)-2.*wk(3)
3688 pijet=pijet+csjet(ii,kk,k+k1-1,i+i1-1,j+j1-1,m,l)
3689 * *wi(i1)*wj(j1)*wk(k1)
3693 ! if(ii.eq.2)print*,' '
3694 ! write(*,'(i2,f6.0,i2,3x,3(2f5.2,2x),f5.2)')
3695 !* ii,sk,k,(wk(kk1),csjet(ii,kk,k+kk1-1,1,1,m,l),kk1=1,3) ,pijet
3698 c-----------------------------------------------------------------------
3699 subroutine MakeCSTable !tabulates psjet
3700 c-----------------------------------------------------------------------
3701 c last two indices of table: parton types
3708 c-----------------------------------------------------------------------
3709 include 'epos.incsem'
3710 common/psar2/edmax,epmax
3711 common/tabcsjet/ksmax,iqmax,jqmax,csjet(0:2,2,20,20,20,3,2)
3712 write (*,'(a,$)')'(CS table'
3718 if(kk.eq.1)spmax=epmax/2.
3719 if(kk.eq.2)spmax=spmin*(epmax/2./spmin)**(1./(ksmax-1.))
3720 do m=1,3 !parton type at upper end of the ladder
3721 write (*,'(a,$)')'.'
3722 do l=1,2 !parton type at lower end of the ladder
3725 if(m.eq.3.and.l.eq.1)l1=-m1
3727 sk=spmin*(spmax/spmin)**((k-1.)/(ksmax-1.))
3730 qi=q2min*(qmax/q2min)**((i-1.)/(iqmax-1.))
3732 qq=qi*(qmax/qi)**((j-1.)/(jqmax-1.))
3733 !write(*,'(i3,4f8.3,2i4,$)')j, qi,q2min,qq,sk,m1,l1
3734 csjet(2,kk,k,i,j,m,l)= psjet(qi,q2min,qq,sk,m1,l1,0)
3735 csjet(1,kk,k,i,j,m,l)=psjet1(qi,q2min,qq,sk,m1,l1,0)
3736 csjet(0,kk,k,i,j,m,l)=psborn(qi,q2min,qq,sk,m1,l1,0,1)
3737 ! if(i.eq.1.and.j.eq.1.and.m.eq.1.and.l.eq.1)
3738 ! *write(*,'(2f8.2,f13.2,2i3,3x,i3,3f8.3)')
3739 ! * qi,qq,sk,m1,l1,k,csjet(2,kk,k,i,j,m,l)
3740 ! * ,csjet(1,kk,k,i,j,m,l),csjet(0,kk,k,i,j,m,l)
3747 write (*,'(a,$)')'done)'
3750 c-----------------------------------------------------------------------
3751 function psjeti(q1,q2,qt,t,xx1,xx2,s,j,l,jdis)
3752 c-----------------------------------------------------------------------
3754 c E~qcd_ji * E~qcd_lk * B_ik
3756 c B_ik = psbori = contribution to Born xsection:
3757 c dsigmaBorn/d2pt/dy
3758 c = s/pi * delta(s+t+u) * 2*pi*alpha**2 /s**2 * B_ik
3760 c E~qcd: at least one emission
3762 c q1 - virtuality cutoff at current end of the ladder
3763 c q2 - virtuality cutoff at opposite end of the ladder
3764 c xx1 - feinman x for the first parton for the born process
3765 c xx2 - feinman x for the second parton for the born process
3766 c s - c.m. energy squared for the born scattering
3767 c t - invariant variable for the scattering |(p1-p3)**2|,
3768 c j - parton type at current end of the ladder (0 - g, 1,-1,2,... - q)
3769 c l - parton type at opposite end of the ladder (0 - g, 1,-1,2,... - q)
3770 c-----------------------------------------------------------------------
3772 c psevi: 1 1 ... gluon -> gluon
3773 c 2 1 ... quark -> gluon
3774 c 1 2 ... gluon -> quark
3775 c 3 2 ... quark -> quark non singlet
3776 c 2 2 ... quark -> quark all
3777 c singlet = all - non singlet
3778 c-----------------------------------------------------------------------
3779 double precision xx1,xx2
3780 include 'epos.incsem'
3781 common/ccctest/iiitest
3788 if(j.eq.0.and.l.eq.0)then ! gluon-gluon --->
3789 akg1=psevi(q1,scale,xx1,1,1) !gluon contribution
3790 akg2=psevi(q2,qt,xx2,1,1) !gluon contribution
3791 aks1=psevi(q1,scale,xx1,1,2)/naflav/2. !singlet contribution per quark
3792 aks2=psevi(q2,qt,xx2,1,2)/naflav/2. !singlet contribution per quark
3793 psjeti=ffborn(s,t,akg1*akg2
3794 * ,(akg1*aks2+aks1*akg2)*naflav*2. !ccccc
3795 * ,aks1*aks2*naflav*2.
3796 * ,aks1*aks2*naflav*2.
3797 * ,aks1*aks2*naflav*2.*(naflav-1)*2.
3799 elseif(j.eq.0)then ! gluon-quark --->
3800 akg1=psevi(q1,scale,xx1,1,1) !gluon contribution
3801 akg2=psevi(q2,qt,xx2,2,1) !gluon contribution
3802 aks1=psevi(q1,scale,xx1,1,2)/naflav/2. !singlet contribution
3803 akns2=psevi(q2,qt,xx2,3,2) !nonsinglet contribution
3804 aks2=(psevi(q2,qt,xx2,2,2)-akns2)/naflav/2. !singlet contribution
3805 psjeti=ffborn(s,t,akg1*akg2
3806 * ,(akg1*(akns2+aks2*naflav*2.)+aks1*akg2*naflav*2.)
3807 * ,aks1*(akns2+aks2*naflav*2.)
3808 * ,aks1*(akns2+aks2*naflav*2.)
3809 * ,aks1*(akns2+aks2*naflav*2.)*(naflav-1)*2.)
3810 elseif(l.eq.0)then ! quark-gluon --->
3811 akg1=psevi(q1,scale,xx1,2,1) !gluon contribution
3812 akg2=psevi(q2,qt,xx2,1,1) !gluon contribution
3813 akns1=psevi(q1,scale,xx1,3,2) !nonsinglet contribution
3814 aks1=(psevi(q1,scale,xx1,2,2)-akns1)/naflav/2. !singlet contribution
3815 aks2=psevi(q2,qt,xx2,1,2)/naflav/2. !singlet contribution
3816 psjeti=ffborn(s,t,akg1*akg2
3817 * ,(akg2*(akns1+aks1*naflav*2.)+aks2*akg1*naflav*2.)
3818 * ,aks2*(akns1+aks1*naflav*2.)
3819 * ,aks2*(akns1+aks1*naflav*2.)
3820 * ,aks2*(akns1+aks1*naflav*2.)*(naflav-1)*2.)
3821 else ! quark-quark --->
3822 akg1=psevi(q1,scale,xx1,2,1) !gluon contribution
3823 akg2=psevi(q2,qt,xx2,2,1) !gluon contribution
3824 akns1=psevi(q1,scale,xx1,3,2) !nonsinglet contribution
3825 aks1=(psevi(q1,scale,xx1,2,2)-akns1)/naflav/2.!singlet contribution
3826 akns2=psevi(q2,qt,xx2,3,2) !nonsinglet contribution
3827 aks2=(psevi(q2,qt,xx2,2,2)-akns2)/naflav/2.!singlet contribution
3830 psjeti=ffborn(s,t,akg1*akg2
3831 * ,(akg2*(akns1+aks1*naflav*2.)+akg1*(akns2+aks2*naflav*2.))
3832 * ,((akns1+aks1)*(akns2+aks2)+aks1*aks2*(2.*naflav-1.))
3833 * ,(akns1*aks2+akns2*aks1+aks1*aks2*naflav*2.)
3834 * ,(akns1*aks2+akns2*aks1+aks1*aks2*naflav*2.)*(naflav-1)*2.)
3836 psjeti=ffborn(s,t,akg1*akg2
3837 * ,(akg2*(akns1+aks1*naflav*2.)+akg1*(akns2+aks2*naflav*2.))
3838 * ,(akns1*aks2+akns2*aks1+aks1*aks2*naflav*2.)
3839 * ,((akns1+aks1)*(akns2+aks2)+aks1*aks2*(2.*naflav-1.))
3840 * ,(akns1*aks2+akns2*aks1+aks1*aks2*naflav*2.)*(naflav-1)*2.)
3842 psjeti=ffborn(s,t,akg1*akg2
3843 * ,(akg2*(akns1+aks1*naflav*2.)+akg1*(akns2+aks2*naflav*2.))
3844 * ,(akns1*aks2+akns2*aks1+aks1*aks2*naflav*2.)
3845 * ,(akns1*aks2+akns2*aks1+aks1*aks2*naflav*2.)
3846 * ,(akns1*akns2+akns1*aks2*(naflav-1)*2.
3847 * +akns2*aks1*(naflav-1)*2.+aks1*aks2*naflav*2.*(naflav-1)*2.))
3853 c-----------------------------------------------------------------------
3854 function psjetj(q1,scale,t,xx,s,j,l,n)
3855 c-----------------------------------------------------------------------
3856 c psjetj - integrand for the ordered ladder cross-section
3857 c q1 - virtuality cutoff at current end of the ladder,
3858 c scale - born process scale,
3859 c t - invariant variable for the scattering |(p1-p3)**2|,
3860 c xx - feinman x for the first parton for the born process
3861 c s - c.m. energy squared for the born scattering,
3862 c j - parton type at current end of the ladder (0 - g, 1,-1,2,... - q)
3863 c l - parton type at opposite end of the ladder (0 - g, 1,-1,2,... - q)
3864 c n - subprocess number
3865 c-----------------------------------------------------------------------
3867 include 'epos.incsem'
3872 psjetj=psevi(q1,scale,xx,m,1)*(psbori(s,t,0,0,n)+ !gg
3873 * psbori(s,s-t,0,0,n))/2.
3874 * +psevi(q1,scale,xx,m,2)*(psbori(s,t,1,0,n)+ !qg
3875 * psbori(s,s-t,1,0,n))
3877 aks=psevi(q1,scale,xx,1,2)/naflav/2. !singlet contribution per quark
3878 psjetj=psevi(q1,scale,xx,1,1)*(psbori(s,t,0,1,n)+ !gq
3879 * psbori(s,s-t,0,1,n))
3880 * +aks*(psbori(s,t,1,1,n)+psbori(s,s-t,1,1,n))/2. !qq
3881 * +aks*(psbori(s,t,-1,1,n)+psbori(s,s-t,-1,1,n)) !qq~
3882 * +aks*(psbori(s,t,1,2,n)+psbori(s,s-t,1,2,n))*(naflav-1)*2. !qq'
3884 akg=psevi(q1,scale,xx,2,1) !gluon contribution
3885 akns=psevi(q1,scale,xx,3,2) !nonsinglet contribution
3886 aks=(psevi(q1,scale,xx,2,2)-akns)/naflav/2. !singlet contribution
3888 psjetj=akg*(psbori(s,t,0,1,n)+psbori(s,s-t,0,1,n)) !gq
3889 * +(akns+aks)*(psbori(s,t,1,1,n)+psbori(s,s-t,1,1,n))/2. !qq
3890 * +aks*(psbori(s,t,-1,1,n)+psbori(s,s-t,-1,1,n)) !qq~
3891 * +aks*(psbori(s,t,1,2,n)+psbori(s,s-t,1,2,n))*(naflav-1)*2. !qq'
3893 psjetj=akg*(psbori(s,t,0,1,n)+psbori(s,s-t,0,1,n)) !gq
3894 * +aks*(psbori(s,t,1,1,n)+psbori(s,s-t,1,1,n))/2. !qq
3895 * +(akns+aks)*(psbori(s,t,-1,1,n)+psbori(s,s-t,-1,1,n)) !qq~
3896 * +aks*(psbori(s,t,1,2,n)+psbori(s,s-t,1,2,n))*(naflav-1)*2.!qq'
3898 psjetj=akg*(psbori(s,t,0,1,n)+psbori(s,s-t,0,1,n)) !gq
3899 * +aks*(psbori(s,t,1,1,n)+psbori(s,s-t,1,1,n))/2. !qq
3900 * +aks*(psbori(s,t,-1,1,n)+psbori(s,s-t,-1,1,n)) !qq~
3901 * +(akns+aks*(naflav-1)*2.)*
3902 * (psbori(s,t,1,2,n)+psbori(s,s-t,1,2,n)) !qq'
3906 p1=s/(1.+qcmass**2/s)
3907 psjetj=psevi(q1,scale,xx,m,1)*(psbori(s,t,4,0,n)+ !cg
3908 * psbori(s,p1-t,4,0,n))
3909 * +psevi(q1,scale,xx,m,2)*(psbori(s,t,4,1,n)+ !cq
3910 * psbori(s,p1-t,4,1,n))
3917 c------------------------------------------------------------------------
3918 function psjti(q1,qqcut,s,m1,l1,jdis)
3919 c-----------------------------------------------------------------------
3920 c psjti - inclusive hard cross-section interpolation - for any ordering
3922 c q1 - virtuality cutoff at current end of the ladder
3923 c qqcut - p_t cutoff for the born process;
3924 c s - total c.m. energy squared for the ladder
3925 c m1 - parton type at current end of the ladder (0-g, 1,2,etc.-q)
3926 c l1 - parton type at opposite end of the ladder (0-g, 1,2,etc.-q)
3927 c-----------------------------------------------------------------------
3928 dimension wi(3),wj(3),wk(3)
3929 common /psar2/ edmax,epmax
3930 common /psar19/ cstot(20,20,240)
3931 include 'epos.incsem'
3939 qqmin=max(q2min,q1/4.)
3945 if(s.le.s2min)return
3948 smins=s2min/(1.-q2ini/qq)
3950 smins=s2min/(1.-q2ini/qq/4.)
3952 if(s.le.smins)goto 1
3955 tmin=2.*qq/(1.+sqrt(1.-4.*qq/s))
3961 if(m1.ne.0.and.m1.eq.l1)then
3964 elseif(m1.ne.0.and.m1.eq.-l1)then
3967 elseif(m1.ne.0.and.l1.ne.0.and.m1.ne.l1)then
3975 ml=20*(m-1)+60*(l-1)+120*jdis
3976 qli=log(q1/q2min)/log(qmax/q2min)*19.+1.
3977 qlj=log(qq/qqmin)/log(s/4./qqmin)*19.+1.
3978 sl=log(s/spmin)/log(epmax/2./spmin)*19.+1.
3989 wi(3)=wi(2)*(wi(2)-1.)*.5
3990 wi(1)=1.-wi(2)+wi(3)
3991 wi(2)=wi(2)-2.*wi(3)
3994 wj(3)=wj(2)*(wj(2)-1.)*.5
3995 wj(1)=1.-wj(2)+wj(3)
3996 wj(2)=wj(2)-2.*wj(3)
3999 wk(3)=wk(2)*(wk(2)-1.)*.5
4000 wk(1)=1.-wk(2)+wk(3)
4001 wk(2)=wk(2)-2.*wk(3)
4006 psjti=psjti+cstot(i+i1-1,j+j1-1,k+k1+ml-1)
4007 * *wi(i1)*wj(j1)*wk(k1)
4011 psjti=exp(psjti)*(1./tmin-1./tmax)
4014 psjti=psbint(q1,q2min,qqcut,s,m1,l1,jdis)
4018 c------------------------------------------------------------------------
4019 subroutine psjti0(ss,sj,sjb,m1,l1)
4020 c-----------------------------------------------------------------------
4021 c psjti0 - inclusive hard cross-section interpolation -
4022 c for minimal virtuality cutoff in the ladder
4023 c s - total c.m. energy squared for the ladder,
4024 c sj - inclusive jet cross-section,
4025 c sjb - born cross-section,
4026 c m1 - parton type at current end of the ladder (0-g, 1,2,etc.-q)
4027 c l1 - parton type at opposite end of the ladder (0-g, 1,2,etc.-q)
4028 c-----------------------------------------------------------------------
4030 common /psar2/ edmax,epmax
4031 common /psar22/ cstotzero(20,4,2),csborzer(20,4,2)
4032 include 'epos.incsem'
4036 if(iabs(m1).ne.4)then
4038 if(m1.ne.0.and.m1.eq.l1)then
4041 elseif(m1.ne.0.and.m1.eq.-l1)then
4044 elseif(m1.ne.0.and.l1.ne.0.and.m1.ne.l1)then
4059 if(s.le.spmin)return
4063 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
4069 sl=log(s/spmin)/log(epmax/2./spmin)*19.+1.
4073 wk(3)=wk(2)*(wk(2)-1.)*.5
4074 wk(1)=1.-wk(2)+wk(3)
4075 wk(2)=wk(2)-2.*wk(3)
4078 sj=sj+cstotzero(k+k1-1,m,l)*wk(k1)
4079 sjb=sjb+csborzer(k+k1-1,m,l)*wk(k1)
4082 sjb=exp(sjb)*(1./tmin-1./tmax)
4083 sj=max(sjb,exp(sj)*(1./tmin-1./tmax))
4087 c------------------------------------------------------------------------
4088 function psjti1(q1,q2,qqcut,s,m1,l1,jdis)
4089 c-----------------------------------------------------------------------
4090 c psjti1 - inclusive hard cross-section interpolation - for strict order
4092 c q1 - virtuality cutoff at current end of the ladder
4093 c q2 - virtuality cutoff at opposite end of the ladder
4094 c qqcut - p_t cutoff for the born process;
4095 c s - total c.m. energy squared for the ladder,
4096 c m1 - parton type at current end of the ladder (0-g, 1,2,etc.-q)
4097 c l1 - parton type at opposite end of the ladder (0-g, 1,2,etc.-q)
4098 c-----------------------------------------------------------------------
4099 dimension wi(3),wj(3),wk(3)
4100 common /psar2/ edmax,epmax
4101 common /psar20/ csord(20,20,240)
4102 include 'epos.incsem'
4103 double precision psuds
4114 if(s.le.s2min)return
4116 smins=s2min/(1.-q2ini/qq)
4117 if(s.le.smins)goto 1
4120 tmin=2.*qq/(1.+sqrt(1.-4.*qq/s))
4126 if(m1.ne.0.and.m1.eq.l1)then
4129 elseif(m1.ne.0.and.m1.eq.-l1)then
4132 elseif(m1.ne.0.and.l1.ne.0.and.m1.ne.l1)then
4140 ml=20*(m-1)+60*(l-1)+120*jdis
4141 qli=log(q1/q2min)/log(s/4./q2min)*19.+1.
4142 qlj=log(qq/qqmin)/log(s/4./qqmin)*19.+1.
4143 sl=log(s/spmin)/log(epmax/2./spmin)*19.+1.
4154 wi(3)=wi(2)*(wi(2)-1.)*.5
4155 wi(1)=1.-wi(2)+wi(3)
4156 wi(2)=wi(2)-2.*wi(3)
4159 wj(3)=wj(2)*(wj(2)-1.)*.5
4160 wj(1)=1.-wj(2)+wj(3)
4161 wj(2)=wj(2)-2.*wj(3)
4164 wk(3)=wk(2)*(wk(2)-1.)*.5
4165 wk(1)=1.-wk(2)+wk(3)
4166 wk(2)=wk(2)-2.*wk(3)
4172 psjti1=psjti1+csord(i+i1-1,j+j1-1,k2)
4173 * *wi(i1)*wj(j1)*wk(k1)
4177 psjti1=exp(psjti1)*(1./tmin-1./tmax)
4179 if(jdis.eq.0.and.qq.gt.q2)then
4180 psjti1=psjti1*sngl(psuds(qq,l1)/psuds(q2,l1))
4181 elseif(jdis.eq.1.and.4.*qq.gt.q2)then
4182 psjti1=psjti1*sngl(psuds(4.*qq,l1)/psuds(q2,l1))
4187 psjti1=psbint(q1,q2,qqcut,s,m1,l1,0)
4189 psjti1=psbint(q2,q1,qqcut,s,l1,m1,1)
4194 c------------------------------------------------------------------------
4195 function pspdfg(xx,qqs,qq,iclpro0,j)
4196 c-----------------------------------------------------------------------
4197 c pspdf - parton distribution function
4198 c qq - virtuality scale
4199 c qqs - initial virtuality for the input distributions
4200 c iclpro0 - hadron class
4202 c-----------------------------------------------------------------------
4204 common/ar3/ x1(7),a1(7)
4205 include 'epos.incsem'
4206 double precision psuds
4208 pspdfg=psdfh4(xx,qqs,0.,iclpro0,j)
4209 if(j.gt.0)pspdfg=pspdfg+psdfh4(xx,qqs,0.,iclpro0,-j) !+sea contr.
4210 pspdfg=pspdfg*sngl(psuds(qq,j)/psuds(qqs,j))
4212 xmin=xx/(1.-q2ini/qq)
4213 if(xmin.ge.1.)return
4218 do i=1,7 !numerical integration over zx
4220 zx=1.-(1.-xm)*(.5+(m-1.5)*x1(i))**.25
4224 aks=psevi(qqs,qq,z,2,1) !quark contribution
4225 akg=psevi(qqs,qq,z,1,1) !gluon contribution
4228 akg=psevi(qqs,qq,z,1,2)/naflav/2. !gluon contribution
4229 akns=psevi(qqs,qq,z,3,2) !nonsinglet contribution
4230 aks=(psevi(qqs,qq,z,2,2)-akns)/naflav/2. !quark contribution
4233 fz=akg*psdfh4(zx,qqs,0.,iclpro0,0)
4234 * +akns*psdfh4(zx,qqs,0.,iclpro0,j)
4235 * +aks*(psdfh4(zx,qqs,0.,iclpro0,1)+
4236 * 2.*psdfh4(zx,qqs,0.,iclpro0,-1)
4237 * +psdfh4(zx,qqs,0.,iclpro0,2)+2.*psdfh4(zx,qqs,0.,iclpro0,-2)
4238 * +2.*psdfh4(zx,qqs,0.,iclpro0,-3))
4239 if(j.gt.0)fz=fz+akns*psdfh4(zx,qqs,0.,iclpro0,-j)
4241 dpd1=dpd1+a1(i)*fz/zx**2/(1.-zx)**3
4244 dpd1=dpd1*(1.-xm)**4/8.*xx
4247 do i=1,7 !numerical integration
4249 zx=xx+(xm-xx)*((xmin-xx)/(xm-xx))**(.5-(m-1.5)*x1(i))
4253 aks=psevi(qqs,qq,z,2,1) !quark contribution
4254 akg=psevi(qqs,qq,z,1,1) !gluon contribution
4257 akg=psevi(qqs,qq,z,1,2)/naflav/2. !gluon contribution
4258 akns=psevi(qqs,qq,z,3,2) !nonsinglet contribution
4259 aks=(psevi(qqs,qq,z,2,2)-akns)/naflav/2. !quark contribution
4262 fz=akg*psdfh4(zx,qqs,0.,iclpro0,0)
4263 * +akns*psdfh4(zx,qqs,0.,iclpro0,j)
4264 * +aks*(psdfh4(zx,qqs,0.,iclpro0,1)
4265 * +2.*psdfh4(zx,qqs,0.,iclpro0,-1)
4266 * +psdfh4(zx,qqs,0.,iclpro0,2)+2.*psdfh4(zx,qqs,0.,iclpro0,-2)
4267 * +2.*psdfh4(zx,qqs,0.,iclpro0,-3))
4268 if(j.gt.0)fz=fz+akns*psdfh4(zx,qqs,0.,iclpro0,-j)
4270 dpd2=dpd2+a1(i)*fz*(1.-xx/zx)/zx
4273 dpd2=dpd2*log((xm-xx)/(xmin-xx))*.5*xx
4275 pspdfg=pspdfg+dpd2+dpd1
4279 c-----------------------------------------------------------------------
4281 c-----------------------------------------------------------------------
4283 include 'epos.incsem'
4288 if(jmod.eq.0)then !??????????????ttttttt
4289 write(*,*)"no more triple Pomeron, xpar2=0 in psaevp not accepted"
4290 write(*,*)"use xpar2=1 instead"
4296 xx=xminim+(xmaxim-xminim)*(i-.5)/nrbins
4298 xx=xminim*(xmaxim/xminim)**((i-.5)/nrbins)
4302 if(jmod.eq.0)then !evolution+matrix element +3P (ours)
4304 ar(i,3)=(psdh(ww,qq,2,0)+psdh(ww,qq,2,1)+
4305 * psdsh1(ww,qq,2,dqsh,0)+psdsh1(ww,qq,2,dqsh,1)
4306 * )/(4.*pi**2*alfe)*qq
4307 elseif(jmod.eq.1)then !evolution+matrix element (ours)
4309 ar(i,3)=(psdh(ww,qq,2,0)+psdh(ww,qq,2,1)+
4310 * psdsh(ww,qq,2,dqsh,0)+psdsh(ww,qq,2,dqsh,1)
4311 * )/(4.*pi**2*alfe)*qq
4312 elseif(jmod.eq.2)then !just evolution (grv)
4313 ar(i,3)=(pspdfg(xx,q2min,qq,2,1)/2.25+
4314 * pspdfg(xx,q2min,qq,2,2)/9.+
4315 * pspdfg(xx,q2min,qq,2,-1)*2./3.6+
4316 * pspdfg(xx,q2min,qq,2,-3)*2./9.)
4317 if(naflav.eq.4)ar(i,3)=ar(i,3)+pspdfg(xx,q2min,qq,2,-4)
4319 elseif(jmod.eq.3)then !grv
4320 ar(i,3)=(psdfh4(xx,qq,0.,2,1)+2.*psdfh4(xx,qq,0.,2,-1))/2.25
4321 * +(psdfh4(xx,qq,0.,2,2)+2.*psdfh4(xx,qq,0.,2,-2))/9.
4322 * +2.*psdfh4(xx,qq,0.,2,-3)/9. !
4323 elseif(jmod.eq.4)then !just evolution (ours)
4324 ar(i,3)=(fparton(xx,qq,1)/2.25+fparton(xx,qq,2)/9.+
4325 * fparton(xx,qq,-1)*6./4.5) !uv+dv+6*sea
4326 if(naflav.eq.4)ar(i,3)=ar(i,3)+fparton(xx,qq,-4)*2./2.25
4327 elseif(jmod.eq.5)then !grv+res
4329 ar(i,3)=(psdgh(ww,qq,0)+psdgh(ww,qq,1)
4330 * )/(4.*pi**2*alfe)*qq
4337 c------------------------------------------------------------------------
4338 subroutine pscs(c,s)
4339 c-----------------------------------------------------------------------
4340 c pscs - cos (c) and sin (s) generation for uniformly distributed angle
4341 c-----------------------------------------------------------------------
4352 c------------------------------------------------------------------------
4353 subroutine psdefrot(ep,s0x,c0x,s0,c0)
4354 c-----------------------------------------------------------------------
4355 c psdefrot - determination of the parameters the spacial rotation to the
4356 c system for 4-vector ep
4357 c s0, c0 - sin and cos for the zx-rotation;
4358 c s0x, c0x - sin and cos for the xy-rotation
4359 c-----------------------------------------------------------------------
4362 c transverse momentum square for the current parton (ep)
4363 pt2=ep(3)**2+ep(4)**2
4366 c system rotation to get pt=0 - euler angles are determined (c0x = cos t
4367 c s0x = sin theta, c0 = cos phi, s0 = sin phi)
4370 c total momentum for the gluon
4371 pl=sqrt(pt2+ep(2)**2)
4388 c------------------------------------------------------------------------
4389 subroutine psdeftr(s,ep,ey)
4390 c-----------------------------------------------------------------------
4391 c psdeftr - determination of the parameters for the lorentz transform to
4392 c rest frame system for 4-vector ep of mass squared s
4393 c-----------------------------------------------------------------------
4395 double precision ep(4)
4398 if(ep(i+1).eq.0.d0)then
4403 if(wp.gt.1.e-8.and.wm/wp.lt.1.e-8)then
4406 if(l.ne.i)ww=ww+ep(l+1)**2
4409 elseif(wm.gt.1.e-8.and.wp/wm.lt.1.e-8)then
4412 if(l.ne.i)ww=ww+ep(l+1)**2
4421 ep(1)=dsqrt(dble(s))
4425 c------------------------------------------------------------------------
4426 function psdfh4(x,qqs,qq,icq,iq)
4427 c------------------------------------------------------------------------
4428 c psdfh4 - GRV structure functions
4429 c------------------------------------------------------------------------
4430 common /psar36/ alvc
4437 sq=log(log(qqs/.232**2)/log(.23/.232**2))
4438 if(iq.eq.0)then !gluon
4444 bg=16.69-22.74*sq+5.779*sq*sq
4445 cg=-25.59+29.71*sq-7.296*sq*sq
4446 dg=2.792+2.215*sq+.422*sq*sq-.104*sq*sq*sq
4449 psdfh4=(1.-x)**dg*(x**aag*(ag+bg*x+cg*x**2)*log(1./x)**bbg
4450 * +sq**alg*exp(-eg+sqrt(eeg*sq**betg*log(1./x))))
4451 elseif(iq.eq.1.or.iq.eq.2)then !u_v or d_v
4454 auu=2.284+.802*sq+.055*sq*sq
4455 au=-.449-.138*sq-.076*sq*sq
4456 bu=.213+2.669*sq-.728*sq*sq
4457 cu=8.854-9.135*sq+1.979*sq*sq
4458 du=2.997+.753*sq-.076*sq*sq
4459 uv=auu*x**aau*(1.-x)**du*
4460 * (1.+au*x**bbu+bu*x+cu*x**1.5)
4464 add=.371+.083*sq+.039*sq*sq
4465 ad=-.509+3.31*sq-1.248*sq*sq
4466 bd=12.41-10.52*sq+2.267*sq*sq
4467 ccd=6.373-6.208*sq+1.418*sq*sq
4468 dd=3.691+.799*sq-.071*sq*sq
4469 dv=add*x**aad*(1.-x)**dd*
4470 * (1.+ad*x**bbd+bd*x+ccd*x**1.5)
4472 if(iq.eq.1)then !u_v
4474 elseif(iq.eq.2)then !d_v
4477 elseif(iq.eq.-3)then !s_sea
4481 as=-5.548+3.669*sqrt(sq)-.616*sq
4482 bs=18.92-16.73*sqrt(sq)+5.168*sq
4483 ds=6.379-.35*sq+.142*sq*sq
4486 psdfh4=(1.-x)**ds*sq**als/log(1./x)**aas*(1.+as*sqrt(x)
4487 * +bs*x)*exp(-es+sqrt(ees*sq**bets*log(1./x)))
4488 elseif(iabs(iq).lt.3)then !u_sea or d_sea
4491 addel=.082+.014*sq+.008*sq*sq
4492 adel=-38.07+36.13*sq-.656*sq*sq
4493 bdel=90.31-74.15*sq+7.645*sq*sq
4495 ddel=7.486+1.217*sq-.159*sq*sq
4496 delv=addel*x**aadel*(1.-x)**ddel*
4497 * (1.+adel*x**bbdel+bdel*x+ccdel*x**1.5)
4506 dud=4.752+1.164*sq+.286*sq*sq
4509 udsea=(1.-x)**dud*(x**aaud*(aud+bud*x+cud*x**2)
4510 * *log(1./x)**bbud+sq**alud*exp(-eud+sqrt(eeud*sq**betud*
4513 if(iq.eq.-1)then !u_sea
4514 psdfh4=(udsea-delv)/2.
4515 elseif(iq.eq.-2)then !d_sea
4516 psdfh4=(udsea+delv)/2.
4521 elseif(icq.eq.1.or.icq.eq.3)then
4522 sq=log(log(qqs/.204**2)/log(.26/.204**2))
4523 if(iq.eq.1.or.iq.eq.2)then
4528 anorm=1.212+.498*sq+.009*sq**2
4529 psdfh4=.5*anorm*x**aapi*(1.-x)**dpi*
4530 * (1.+api*sqrt(x)+bpi*x)
4534 aapi=2.251-1.339*sqrt(sq)
4535 api=2.668-1.265*sq+.156*sq**2
4538 cpi=-1.014+.92*sq-.101*sq**2
4542 psdfh4=(1.-x)**dpi*(x**aapi*(api+bpi*sqrt(x)+cpi*x)*
4543 * log(1./x)**bbpi+sq**alfpi*
4544 * exp(-epi+sqrt(eppi*sq**betpi*log(1./x))))
4545 elseif(iq.eq.-3)then
4554 psdfh4=sq**alfpi/log(1./x)**aapi*(1.-x)**dpi*
4555 * (1.+api*sqrt(x)+bpi*x)*
4556 * exp(-epi+sqrt(eppi*sq**betpi*log(1./x)))
4557 elseif(iabs(iq).lt.3)then
4560 aapi=.309-.134*sqrt(sq)
4562 bbpi=.893-.264*sqrt(sq)
4568 psdfh4=(1.-x)**dpi*(x**aapi*(api+bpi*sqrt(x)+cpi*x)*
4569 * log(1./x)**bbpi+sq**alfpi*
4570 * exp(-epi+sqrt(eppi*sq**betpi*log(1./x))))
4574 elseif(icq.eq.0)then
4575 sq=log(log(qqs/.204**2)/log(.26/.204**2))
4579 aapi=2.251-1.339*sqrt(sq)
4580 api=2.668-1.265*sq+.156*sq**2
4583 cpi=-1.014+.92*sq-.101*sq**2
4587 psdfh4=(1.-x)**dpi*(x**aapi*(api+bpi*sqrt(x)+cpi*x)*
4588 * log(1./x)**bbpi+sq**alfpi*
4589 * exp(-epi+sqrt(eppi*sq**betpi*log(1./x))))
4600 str=sq**alfpi/log(1./x)**aapi*(1.-x)**dpi*
4601 * (1.+api*sqrt(x)+bpi*x)*
4602 * exp(-epi+sqrt(eppi*sq**betpi*log(1./x)))
4610 anorm=1.212+.498*sq+.009*sq**2
4611 val=.5*anorm*x**aapi*(1.-x)**dpi*
4612 * (1.+api*sqrt(x)+bpi*x)
4616 aapi=.309-.134*sqrt(sq)
4618 bbpi=.893-.264*sqrt(sq)
4624 sea=(1.-x)**dpi*(x**aapi*(api+bpi*sqrt(x)+cpi*x)*
4625 * log(1./x)**bbpi+sq**alfpi*
4626 * exp(-epi+sqrt(eppi*sq**betpi*log(1./x))))
4628 psdfh4=(.836*(val+2.*sea)-.587*str)
4630 psdfh4=(.25*(val+2.*sea)+.587*str)
4636 psdfh4=psdfh4/(1.+qq/.59)**2
4638 elseif(icq.eq.4.and.iq.eq.1)then
4639 psdfh4=x**3*(1.-x)**alvc*(alvc+3.)*(alvc+2.)*(alvc+1.)
4647 c------------------------------------------------------------------------
4648 function psfap(x,j,l)
4649 c-----------------------------------------------------------------------
4650 c psfap - altarelli-parisi function (multiplied by x)
4651 c x - light cone momentum share value,
4652 c j - type of the parent parton (0-g;1,2,etc.-q)
4653 c l - type of the daughter parton (0-g;1,2,etc.-q)
4654 c-----------------------------------------------------------------------
4656 include 'epos.incsem'
4660 psfap=((1.d0-x)/x+x/(1.d0-x)+x*(1.d0-x))*6.d0
4662 psfap=(x**2+(1.d0-x)**2)*naflav
4666 psfap=(1.d0+(1.d0-x)**2)/x/.75d0
4668 psfap=(x**2+1.d0)/(1.d0-x)/.75d0
4674 cc------------------------------------------------------------------------
4675 c function psgen(a1,a2)
4676 cc-----------------------------------------------------------------------
4677 cc psgen - x-values generation according to distribution
4678 cc x1^(-a1) x2^(-0.5)
4679 cc-----------------------------------------------------------------------
4680 c common/lept1/engy,elepti,elepto,angmue,icinpu
4685 c x1=.5*rangen()**(1./(1.-aa))
4686 c elseif(aa.eq.1.)then
4687 c x1=.5/engy**rangen()
4689 c x1=.5*(1.+rangen()*(engy**(aa-1.)-1.))**(1./(1.-aa))
4691 c if(x1.lt.1.e-7.or.x1.gt..999999)then
4694 c if(rangen().lt..5)then
4695 c gb=x1**(aa-a1)*.5**aa/(1.-x1)**a2
4698 c gb=(1.-x1)**(aa-a2)*.5**aa/x1**a1
4700 c if(rangen().gt.gb)goto 1
4705 c------------------------------------------------------------------------
4707 c-----------------------------------------------------------------------
4708 c psidd - kink type decoder
4709 c-----------------------------------------------------------------------
4712 elseif(iabs(icc).le.2)then !u,u~,d,d~
4714 elseif(iabs(icc).eq.4)then !s,s~
4716 elseif(iabs(icc).gt.10)then !c,c~ etc.
4718 elseif(icc.eq.3)then !ud
4720 elseif(icc.eq.-3)then !u~d~
4722 elseif(icc.eq.6)then !uu
4724 elseif(icc.eq.-6)then !u~u~
4726 elseif(icc.eq.7)then !dd
4728 elseif(icc.eq.-7)then !d~d~
4731 write (*,*)'psidd?????????',icc
4736 cc------------------------------------------------------------------------
4737 c function pslam(s,a,b)
4738 cc-----------------------------------------------------------------------
4739 cc kinematical function for two particle decay - maximal pt-value
4740 cc a - first particle mass squared,
4741 cc b - second particle mass squared,
4742 cc s - two particle invariant mass squared
4743 cc-----------------------------------------------------------------------
4744 c pslam=.25/s*(s+a-b)**2-a
4748 c------------------------------------------------------------------------
4749 function psjvrg1(qt,s,y0)
4750 c-----------------------------------------------------------------------
4751 common /ar3/ x1(7),a1(7)
4752 common /cnsta/ pi,pii,hquer,prom,piom,ainfin
4753 include 'epos.incsem'
4754 double precision xt,ymin,ymax,y,xmin,xmax,xx1,xx2
4757 if(s.le.4.*qt)return
4759 xt=2.d0*sqrt(dble(qt)/dble(s))
4760 ymax=min(dble(y0),log(1d0/xt+sqrt((1d0/xt-1d0)*(1d0/xt+1d0))))
4765 y=.5d0*(ymax+ymin+(ymin-ymax)*dble((2*m-3)*x1(i)))
4766 xmin=xt**2/2.d0/(2.d0-xt*exp(-y))
4767 xmax=1.d0-xt*exp(y)/2.d0
4772 xx1=xt*exp(y)/2d0+xmin*(xmax/xmin)**dble(.5+x1(i1)*(m1-1.5))
4773 xx2=xt*exp(-y)*xx1/(2.d0*xx1-xt*exp(y))
4776 t=sngl(dble(sh)/2d0*(1d0
4777 & -sqrt(max(0d0,1d0-4d0*dble(qt)/dble(sh)))))
4778 ft=psjvrx(t,qt,sngl(xx1),sngl(xx2),sh)
4779 fx=fx+a1(i1)*ft/sh**2
4782 fx=fx*sngl(log(xmax/xmin))
4783 psjvrg1=psjvrg1+a1(i)*fx
4786 psjvrg1=psjvrg1*sngl(ymax-ymin)*pi**3
4787 **pssalf(qt/qcdlam)**2*sqrt(qt)
4791 c-----------------------------------------------------------------------
4792 function psjvrx(t,qt,xx1,xx2,s)
4793 c-----------------------------------------------------------------------
4794 include 'epos.incsem'
4796 g1=psdfh4(xx1,qt,0.,2,0)
4797 ub1=psdfh4(xx1,qt,0.,2,-1)
4798 u1=psdfh4(xx1,qt,0.,2,1)+ub1
4799 db1=psdfh4(xx1,qt,0.,2,-2)
4800 d1=psdfh4(xx1,qt,0.,2,2)+db1
4801 sb1=psdfh4(xx1,qt,0.,2,-3)
4803 g2=psdfh4(xx2,qt,0.,2,0)
4804 ub2=psdfh4(xx2,qt,0.,2,-1)
4805 u2=psdfh4(xx2,qt,0.,2,1)+ub2
4806 db2=psdfh4(xx2,qt,0.,2,-2)
4807 d2=psdfh4(xx2,qt,0.,2,2)+db2
4808 sb2=psdfh4(xx2,qt,0.,2,-3)
4811 psjvrx=g1*g2*(psbori(s,t,0,0,1)+psbori(s,s-t,0,0,1)
4812 *+psbori(s,t,0,0,2)+psbori(s,s-t,0,0,2))/2.
4813 *+(psbori(s,t,0,1,1)+psbori(s,s-t,0,1,1))*
4814 *(g2*(u1+ub1+d1+db1+s1+sb1)+g1*(u2+ub2+d2+db2+s2+sb2))
4815 *+(psbori(s,t,1,1,1)+psbori(s,s-t,1,1,1))/2.*
4816 *(u1*u2+ub1*ub2+d1*d2+db1*db2+s1*s2+sb1*sb2)
4817 *+(psbori(s,t,1,-1,1)+psbori(s,s-t,1,-1,1)+psbori(s,t,1,-1,2)+
4818 *psbori(s,s-t,1,-1,2)+psbori(s,t,1,-1,3)+psbori(s,s-t,1,-1,3))*
4819 *(u1*ub2+ub1*u2+d1*db2+db1*d2+s1*sb2+sb1*s2)
4820 *+(psbori(s,t,1,2,1)+psbori(s,s-t,1,2,1))*
4821 *((u1+ub1)*(d2+db2+s2+sb2)+(u2+ub2)*(d1+db1+s1+sb1)+
4822 *(d1+db1)*(u2+ub2+s2+sb2)+(d2+db2)*(u1+ub1+s1+sb1)+
4823 *(s1+sb1)*(u2+ub2+d2+db2)+(s2+sb2)*(u1+ub1+d1+db1))
4827 c------------------------------------------------------------------------
4828 function psjwo1(qt,s,y0)
4829 c-----------------------------------------------------------------------
4830 common /ar3/ x1(7),a1(7)
4831 common /cnsta/ pi,pii,hquer,prom,piom,ainfin
4832 double precision xt,ymax,ymin,y,xmin,xmax,xx1,xx2
4833 include 'epos.incsem'
4836 if(s.le.4.*qt)return
4838 xt=2.d0*sqrt(dble(qt)/dble(s))
4839 ymax=min(dble(y0),log(1d0/xt+sqrt((1d0/xt-1d0)*(1d0/xt+1d0))))
4844 y=.5d0*(ymax+ymin+(ymin-ymax)*dble(2*m-3)*dble(x1(i)))
4845 xmin=xt**2/2.d0/(2.d0-xt*exp(-y))
4846 xmax=1.d0-xt*exp(y)/2.d0
4851 xx1=xt*exp(y)/2.d0+xmin*(xmax/xmin)**dble(.5+x1(i1)*(m1-1.5))
4852 xx2=xt*exp(-y)/(2.d0-xt*exp(y)/xx1)
4855 t=sngl(dble(sh)/2d0*(1d0-sqrt(1d0-4d0*dble(qt)/dble(sh))))
4856 ft=psjwox(t,qt,sngl(xx1),sngl(xx2),sh)
4857 fx=fx+a1(i1)*ft/sh**2
4860 fx=fx*log(xmax/xmin)
4861 psjwo1=psjwo1+a1(i)*fx
4864 psjwo1=psjwo1*sngl(ymax-ymin)*pi**3
4865 **pssalf(qt/qcdlam)**2*sqrt(qt)
4869 c-----------------------------------------------------------------------
4870 function psjwox(t,qt,xx1,xx2,s)
4871 c-----------------------------------------------------------------------
4872 double precision x,scale,upv1,dnv1,sea1,str1,chm1,gl1,
4873 *upv2,dnv2,sea2,str2,chm2,gl2
4876 call strdo1(x,scale,upv1,dnv1,sea1,str1,chm1,gl1)
4878 call strdo1(x,scale,upv2,dnv2,sea2,str2,chm2,gl2)
4880 psjwox=gl1*gl2*(psbori(s,t,0,0,1)+psbori(s,s-t,0,0,1)
4881 *+psbori(s,t,0,0,2)+psbori(s,s-t,0,0,2)+psbori(s,t,0,0,3)
4882 *+psbori(s,s-t,0,0,3))/2.
4883 *+(psbori(s,t,0,1,1)+psbori(s,s-t,0,1,1)
4884 *+psbori(s,t,0,1,2)+psbori(s,s-t,0,1,2)+psbori(s,t,0,1,3)
4885 *+psbori(s,s-t,0,1,3))*(gl2*(upv1+dnv1+4.*sea1+2.*str1+2.*chm1)+
4886 *gl1*(upv2+dnv2+4.*sea2+2.*str2+2.*chm2))
4887 *+(psbori(s,t,1,1,1)+psbori(s,s-t,1,1,1)
4888 *+psbori(s,t,1,1,2)+psbori(s,s-t,1,1,2)+psbori(s,t,1,1,3)+
4889 *psbori(s,s-t,1,1,3))/2.*
4890 *((upv1+sea1)*(upv2+sea2)+(dnv1+sea1)*(dnv2+sea2)+2.*sea1*sea2
4891 *+2.*str1*str2+2.*chm1*chm2)
4892 *+(psbori(s,t,1,-1,1)+psbori(s,s-t,1,-1,1)+psbori(s,t,1,-1,2)+
4893 *psbori(s,s-t,1,-1,2)+psbori(s,t,1,-1,3)+psbori(s,s-t,1,-1,3))*
4894 *((upv1+sea1)*sea2+sea1*(upv2+sea2)+(dnv1+sea1)*sea2+
4895 *sea1*(dnv2+sea2)+2.*str1*str2+2.*chm1*chm2)
4896 *+(psbori(s,t,1,2,1)
4897 *+psbori(s,s-t,1,2,1)+psbori(s,t,1,2,2)+psbori(s,s-t,1,2,2)
4898 *+psbori(s,t,1,2,3)+psbori(s,s-t,1,2,3))*
4899 *(upv1*dnv2+upv2*dnv1+(upv1+dnv1)*(2.*sea2+2.*str2+2.*chm2)+
4900 *(upv2+dnv2)*(2.*sea1+2.*str1+2.*chm1)+
4901 *4.*sea1*(2.*sea2+2.*str2+2.*chm2)+2.*str1*(4.*sea2+2.*chm2)+
4902 *2.*chm1*(4.*sea2+2.*str2))
4906 c------------------------------------------------------------------------
4907 subroutine pslcsh(wp1,wm1,wp2,wm2,samqt,amqpt)
4908 c-----------------------------------------------------------------------
4909 c pslcsh - sh pomeron lc momentum sharing between two strings
4910 c------------------------------------------------------------------------
4911 double precision amqt(4),yqm(4),yqm1(4),xlp(4),xlm(4),am23,sx,y2
4912 *,wp1,wp2,wm1,wm2,s,sq,psutz,yqmax,y,amjp,amjm,y1,s12,s34,x34,amqpt
4919 amqt(i)=dble(samqt(i))
4920 yqm(i)=dlog(sq/amqt(i)*psutz(s,amqt(i)**2,(amqpt-amqt(i))**2))
4922 yqmax=max(yqm(1),yqm(2))
4924 1 y=yqmax*dble(rangen())
4926 if(y.gt.yqm(j))goto 1
4928 amjp=amqt(j)*dexp(y)
4929 amjm=amqt(j)*dexp(-y)
4931 am23=amqt(3-j)+amqt(7-i)
4932 sx=(am23+amjp)*(am23+amjm)
4933 yqm1(i)=dlog(sq/amqt(i)*psutz(s,amqt(i)**2,sx))
4935 yqmax1=max(yqm1(3),yqm1(4))
4936 if(dble(rangen()).gt.yqmax1/max(yqm(3),yqm(4)))goto 1
4938 y1=yqmax1*dble(rangen())
4939 j1=int(3.5+rangen())
4940 if(y1.gt.yqm1(j1))goto 1
4942 amjp1=amqt(j1)*exp(y1)
4943 amjm1=amqt(j1)*exp(-y1)
4944 s12=(amqt(3-j)+amjp)*(amqt(3-j)+amjm)
4945 s34=(amqt(7-j1)+amjp1)*(amqt(7-j1)+amjm1)
4946 y2=dlog(sq/(amqt(3-j)+amjp)*psutz(s,s12,s34))
4948 xlp(j)=amqt(j)/sq*dexp(y+y2)
4949 xlm(j)=amqt(j)/sq*dexp(-y-y2)
4950 xlp(3-j)=amqt(3-j)/sq*dexp(y2)
4951 xlm(3-j)=amqt(3-j)/sq*dexp(-y2)
4952 x34=1.-xlm(1)-xlm(2)
4953 xlm(7-j1)=x34/(1.+amjp1/amqt(7-j1))
4954 xlm(j1)=x34-xlm(7-j1)
4955 c write (*,*)'xlc',xlp(1),xlp(2),xlm(3),xlm(4)
4956 if(dble(rangen()).gt.(xlp(1)*xlp(2)*xlm(3)*xlm(4))**(-alpqua)*
4957 *(xlp(j)*(1.d0-xlp(j))*xlm(j1)*(1.d0-xlm(j1))))goto 1
4963 c write (*,*)'wp1,wm1,wp2,wm2',wp1,wm1,wp2,wm2
4967 c------------------------------------------------------------------------
4969 c-----------------------------------------------------------------------
4970 c 4-vector squared calculation
4971 c-----------------------------------------------------------------------
4972 double precision sm2,ep(4)
4981 c------------------------------------------------------------------------
4982 subroutine psrotat(ep,s0x,c0x,s0,c0)
4983 c-----------------------------------------------------------------------
4984 c psrotat - spacial rotation to the lab. system for 4-vector ep
4985 c s0, c0 - sin and cos for the zx-rotation;
4986 c s0x, c0x - sin and cos for the xy-rotation
4987 c-----------------------------------------------------------------------
4988 dimension ep(4),ep1(3)
4991 ep1(2)=ep(2)*s0+ep(3)*c0
4992 ep1(1)=ep(2)*c0-ep(3)*s0
4995 ep(4)=ep1(2)*s0x+ep1(3)*c0x
4996 ep(3)=ep1(2)*c0x-ep1(3)*s0x
5000 cc------------------------------------------------------------------------
5001 c subroutine psrotat1(ep,s0x,c0x,s0,c0)
5002 cc-----------------------------------------------------------------------
5003 cc psrotat - spacial rotation to the lab. system for 4-vector ep
5004 cc s0, c0 - sin and cos for the zx-rotation;
5005 cc s0x, c0x - sin and cos for the xy-rotation
5006 cc-----------------------------------------------------------------------
5007 c dimension ep(4),ep1(3)
5010 c ep1(3)=-ep(3)*s0x+ep(4)*c0x
5011 c ep1(2)=ep(3)*c0x+ep(4)*s0x
5014 c ep(3)=-ep1(1)*s0+ep1(2)*c0
5015 c ep(2)=ep1(1)*c0+ep1(2)*s0
5019 c-----------------------------------------------------------------------
5021 c-----------------------------------------------------------------------
5022 c pssalf - effective qcd coupling (alpha_s/2/pi)
5023 c-----------------------------------------------------------------------
5024 include "epos.incsem"
5025 pssalf=2./(11.-naflav/1.5)/log(qq)
5029 c------------------------------------------------------------------------
5030 subroutine pstrans(ep,ey,jj)
5031 c-----------------------------------------------------------------------
5032 c pstrans - lorentz boosts according to the parameters ey ( determining
5033 c shift along the z,x,y-axis respectively (ey(1),ey(2),ey(3)))
5034 c jj=1 - inverse transformation to the lab. system;
5035 c jj=-1 - direct transformation
5036 c-----------------------------------------------------------------------
5037 dimension ey(3),ep(4)
5040 c lorentz transform to lab. system according to 1/ey(i) parameters
5042 if(ey(4-i).ne.1.)then
5043 wp=(ep(1)+ep(5-i))/ey(4-i)
5044 wm=(ep(1)-ep(5-i))*ey(4-i)
5050 c lorentz transform to lab. system according to ey(i) parameters
5053 wp=(ep(1)+ep(i+1))*ey(i)
5054 wm=(ep(1)-ep(i+1))/ey(i)
5063 c------------------------------------------------------------------------
5064 double precision function psuds(q,m)
5065 c-----------------------------------------------------------------------
5066 c psuds - spacelike sudakov formfactor
5067 c q - maximal value of the effective momentum,
5068 c m - type of parton (0 - g, 1,2, etc. - q)
5069 c-----------------------------------------------------------------------
5071 common /psar15/ sudx(40,2)
5072 include 'epos.incsem'
5073 double precision dps,qlm,ffacs,qlm0,qlmi
5078 qli=log(q/q2min)*2.+1.
5083 wi(3)=wi(2)*(wi(2)-1.)*.5
5084 wi(1)=1.-wi(2)+wi(3)
5085 wi(2)=wi(2)-2.*wi(3)
5088 dps=dps+dble(sudx(i+i1-1,j)*wi(i1))
5091 qlm0=dble(log(q2ini/qcdlam))
5092 qlm=dble(log(q/qcdlam))
5093 qlmi=qlm-qlm0 !=log(q/q2ini)
5094 psuds=(qlm*log(qlm/qlm0)-qlmi)
5096 ffacs=(11.d0-dble(naflav)/1.5d0)/12.d0
5098 psuds=(psuds-ffacs*log(qlm/qlm0)
5099 * +dps*(1.d0-dble(q2ini/q)))/ffacs
5101 psuds=(psuds-log(qlm/qlm0)*.75d0
5102 * +dps*(1.d0-dble(q2ini/q)))*4.d0/9.d0/ffacs
5111 c------------------------------------------------------------------------
5113 c-----------------------------------------------------------------------
5114 c psudx - part of the bspacelike sudakov formfactor
5115 c q - maximal value of the effective momentum,
5116 c j - type of parton (1 - g, 2 - q)
5117 c-----------------------------------------------------------------------
5118 common /ar3/ x1(7),a1(7)
5119 include 'epos.incsem'
5125 qt=.5*(q2ini+q-x1(i)*(2.*m-3.)*(q2ini-q))
5128 dps=((11.-naflav/1.5)/12.-zm**2*(1.-naflav/12.)+
5129 * (zm**3/3.-zm**4/4.)*(1.-naflav/3.))*q/qt
5133 psudx=psudx+a1(i)*dps/log(qt/qcdlam)
5140 c------------------------------------------------------------------------
5141 double precision function psutz(s,a,b)
5142 c-----------------------------------------------------------------------
5143 c psutz - kinematical function for two particle decay - light cone momen
5144 c share for the particle of mass squared a,
5145 c b - partner's mass squared,
5146 c s - two particle invariant mass
5147 c-----------------------------------------------------------------------
5148 double precision a1,b1,s1,x,dx,s,a,b
5152 x=(1.d0+(a1-b1)*(a1+b1)/s)/2.d0
5153 dx=(x-a1/s1)*(x+a1/s1)
5161 psutz=min(0.999999999d0,x)
5165 c------------------------------------------------------------------------
5167 c-----------------------------------------------------------------------
5168 c constants for numerical integration (gaussian weights)
5169 c-----------------------------------------------------------------------
5170 common /ar3/ x1(7),a1(7)
5171 common /ar4/ x4(2),a4(2)
5172 common /ar5/ x5(2),a5(2)
5173 common /ar8/ x2(4),a2
5174 common /ar9/ x9(3),a9(3)
5176 data x1/.9862838,.9284349,.8272013,.6872929,.5152486,
5178 data a1/.03511946,.08015809,.1215186,.1572032,
5179 *.1855384,.2051985,.2152639/
5180 data x2/.00960736,.0842652,.222215,.402455/
5182 data x4/ 0.339981,0.861136/
5183 data a4/ 0.652145,0.347855/
5184 data x5/.585786,3.41421/
5185 data a5/.853553,.146447/
5186 data x9/.93247,.661209,.238619/
5187 data a9/.171324,.360762,.467914/
5190 c------------------------------------------------------------------------
5191 subroutine strdo1(x,scale,upv,dnv,sea,str,chm,gl)
5192 c------------------------------------------------------------------------
5193 c :::::::::::: duke owens set 1 ::::::::::::::::::::::::::::
5194 c------------------------------------------------------------------------
5195 implicit double precision(a-h,o-z)
5197 + f(5),a(6,5),b1(3,6,5)
5198 data q0,ql1/2.d0,.2d0/
5199 data b1/3.d0,0.d0,0.d0,.419d0,.004383d0,-.007412d0,
5200 &3.46d0,.72432d0,-.065998d0,4.4d0,-4.8644d0,1.3274d0,
5202 &0.d0,0.d0,.763d0,-.23696d0,.025836d0,4.d0,.62664d0,-.019163d0,
5203 &0.d0,-.42068d0,.032809d0,6*0.d0,1.265d0,-1.1323d0,.29268d0,
5204 &0.d0,-.37162d0,-.028977d0,8.05d0,1.5877d0,-.15291d0,
5205 &0.d0,6.3059d0,-.27342d0,0.d0,-10.543d0,-3.1674d0,
5206 &0.d0,14.698d0,9.798d0,0.d0,.13479d0,-.074693d0,
5207 &-.0355d0,-.22237d0,-.057685d0,6.3494d0,3.2649d0,-.90945d0,
5208 &0.d0,-3.0331d0,1.5042d0,0.d0,17.431d0,-11.255d0,
5209 &0.d0,-17.861d0,15.571d0,1.564d0,-1.7112d0,.63751d0,
5210 &0.d0,-.94892d0,.32505d0,6.d0,1.4345d0,-1.0485d0,
5211 &9.d0,-7.1858d0,.25494d0,0.d0,-16.457d0,10.947d0,
5212 &0.d0,15.261d0,-10.085d0/
5214 s= log( log( max(q0,scale)/ql1)/ log(q0/ql1))
5217 10 a(j,i)=b1(1,j,i)+s*(b1(2,j,i)+s*b1(3,j,i))
5219 40 f(i)=a(1,i)*x**a(2,i)*(wn-x)**a(3,i)*(wn+x*
5220 & (a(4,i)+x*(a(5,i)+x*a(6,i))))
5223 50 f(i)=f(i)*utgam2(aa)/((wn+a(2,i)*a(4,i)/aa)
5224 &*utgam2(a(2,i))*utgam2(wn+a(3,i)))
5236 c------------------------------------------------------------------------
5237 function fzeroGluZZ(z,k) ! former psftild
5238 c-----------------------------------------------------------------------
5240 c fzeroGluZZComplete = fzeroGluZZ * z^(-1-dels) * gamsoft * gamhad
5242 c A = 8*pi*s0*gampar*gamtilde
5243 c integration over semihard pomeron light cone momentum share xp==u
5245 c fzeroGluZZ = (1-glusea) * engy^epszero
5246 c * int(du) u^(epszero-alppar+dels) (1-u)^alplea * (1-z/u)**betpom
5248 c z - light cone x of the gluon,
5250 c-----------------------------------------------------------------------
5251 double precision xpmin,xp
5253 common /ar3/ x1(7),a1(7)
5254 include 'epos.incsem'
5258 xpmin=xpmin**(1.-alppar+dels+epszero)
5261 xp=(.5*(1.+xpmin+(2*m-3)*x1(i)*(1.-xpmin)))**(1./
5262 * (1.-alppar+dels+epszero))
5263 fzeroGluZZ=fzeroGluZZ+a1(i)*(1.-xp)**alplea(k)*(1.-z/xp)**betpom
5267 * fzeroGluZZ*.5*(1.-xpmin)/(1.-alppar+dels+epszero)
5268 * *(1.-glusea) *engy**epszero
5272 c------------------------------------------------------------------------
5273 function fzeroSeaZZ(z,k) ! former psftile
5274 c-----------------------------------------------------------------------
5276 c fzeroSeaZZComplete = fzeroSeaZZ * z^(-1-dels) * gamsoft * gamhad
5278 c gamsoft = 8*pi*s0*gampar*gamtilde
5279 c integration over semihard pomeron light cone momentum share xp==u
5281 c fzeroSeaZZ = glusea * engy^epszero
5282 c * int(du) u^(epszero-alppar+dels) (1-u)^alplea * EsoftQZero(z/u)
5284 c z - light cone x of the quark,
5286 c-----------------------------------------------------------------------
5287 double precision xpmin,xp
5288 common /ar3/ x1(7),a1(7)
5290 include 'epos.incsem'
5294 xpmin=xpmin**(1.-alppar+dels+epszero)
5297 xp=(.5*(1.+xpmin+(2*m-3)*x1(i)*(1.-xpmin)))**(1./
5298 * (1.-alppar+dels+epszero))
5300 fzeroSeaZZ=fzeroSeaZZ+a1(i)*(1.-xp)**alplea(k)*EsoftQZero(zz)
5303 fzeroSeaZZ=fzeroSeaZZ*.5*(1.-xpmin)/(1.-alppar+dels+epszero)
5304 * *glusea *engy**epszero
5309 c########################################################################
5310 c########################################################################
5312 c########################################################################
5313 c########################################################################
5315 c-----------------------------------------------------------------------
5316 c common initialization procedure
5317 c if isetcs = 0, alpD, betD, etc ... in inirj are not used and xkappa=1
5318 c if isetcs = 1, alpD, betD, etc ... in inirj are not used but xkappa.ne.1
5319 c if isetcs = 2, alpD, betD, xkappa, etc ... in inirj are used and
5320 c cross section from calculation in inics are read.
5321 c if epos.inics doesn't exist, it produces only the calculated part of it.
5322 c if isetcs = 3, alpD, betD, xkappa, etc ... in inirj are used and
5323 c cross section from simulation in inics are read.
5324 c if epos.inics doesn't exist, it produces the calculated AND the
5325 c simulated part of it.
5326 c-----------------------------------------------------------------------
5328 include 'epos.incpar'
5329 include 'epos.incsem'
5330 logical lcalc!,lcalc2
5331 c double precision om5p,xh,yh,v3pom(4),om2p
5332 dimension gamhad0(nclha),r2had0(nclha),chad0(nclha)
5333 *,alplea0(nclha),asect1(7,4,7),asect2(7,4,7),asect3(7,7,7)
5334 *,asect4(7,7,7)!,cgam(idxD)
5335 common /psar2/ edmax,epmax
5336 common /psar4/ fhgg(11,10,8),fhqg(11,10,80)
5337 *,fhgq(11,10,80),fhqq(11,10,80),fhgg0(11,10),fhgg1(11,10,4)
5338 *,fhqg1(11,10,40),fhgg01(11),fhgg02(11),fhgg11(11,4)
5339 *,fhgg12(11,4),fhqg11(11,10,4),fhqg12(11,10,4)
5340 *,ftoint(11,14,2,2,3)
5341 common /psar7/ delx,alam3p,gam3p
5343 common /psar15/ sudx(40,2)
5344 common /psar19/ cstot(20,20,240)
5345 common /psar20/ csord(20,20,240)
5346 common /psar21/ csbor(20,160,2)
5347 common /psar22/ cstotzero(20,4,2),csborzer(20,4,2)
5348 common /psar23/ cschar(20,20,2)
5349 common /psar25/ csdsi(21,21,104)
5350 common /psar27/ csds(21,26,4),csdt(21,26,2),csdr(21,26,2)
5351 common /psar33/ asect(7,4,7),asectn(7,7,7)
5352 common /psar34/ rrr,rrrm
5353 common /psar35/ anorm,anormp
5354 common /psar41/ rrrp,rrrmp
5355 common /psar36/ alvc
5356 common /psar37/ coefom1,coefom2
5357 common /psar38/ vfro(11,14,3,2)
5358 common /psar39/ vnorm(11,14,3,2,2)
5359 c$$$ common /psar40/ coefxu1(idxD,nclha,10)
5360 c$$$ *,coefxu2(idxD,idxD,nclha,10),coefxc2(idxD,idxD,nclha,10)
5361 common /ar3/ x1(7),a1(7)
5362 common /testj/ ajeth(4),ajete(5),ajet0(7)
5363 parameter(nbkbin=40)
5364 common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
5365 common/geom/rmproj,rmtarg,bmax,bkmx
5366 character textini*38
5367 external ptfau,ptfauAA
5370 call utpri('psaini',ish,ishini,4)
5383 if(isetcs.le.1)then !for Kfit
5390 edmax=edmaxi !1.e12 defined in epos-bas
5391 epmax=epmaxi !1.e12 defined in epos-bas
5393 c fix enhanced diagrams at minimum energy = 2.5
5394 delx=1.5 !sqrt(egymin*egymin/exp(1.))
5400 c interface to 'bas'
5404 alpqua=(alppar+1.)/2.
5405 if(abs(alpqua).lt.1.e-6)call utstop('alpar should not be -1 !&')
5406 alpr=-2.+alpqua !x-exponent for remnant mass
5410 coefom0=utgam1(1.+dels-alppar)*utgam1(1.+alplea(iclpro))
5411 */utgam1(2.+alplea(iclpro)+dels-alppar)
5412 **utgam1(1.+dels-alppar)*utgam1(1.+alplea(icltar))
5413 */utgam1(2.+alplea(icltar)+dels-alppar)
5414 coefom1=1.-utgam1(1.+dels-alppar)**2*utgam1(1.+alplea(iclpro))
5415 */utgam1(1.+alplea(iclpro)+2.*(1.+dels-alppar))
5416 **utgam1(1.+dels-alppar)**2*utgam1(1.+alplea(icltar))
5417 */utgam1(1.+alplea(icltar)+2.*(1.+dels-alppar))/coefom0**2
5418 coefom2=3.*coefom1-1.
5419 *+utgam1(1.+dels-alppar)**3*utgam1(1.+alplea(iclpro))
5420 */utgam1(1.+alplea(iclpro)+3.*(1.+dels-alppar))
5421 **utgam1(1.+dels-alppar)**3*utgam1(1.+alplea(icltar))
5422 */utgam1(1.+alplea(icltar)+3.*(1.+dels-alppar))/coefom0**3
5423 if(ish.ge.4)write(ifch,*)'coefom',coefom0,coefom1,coefom2,delx
5425 c soft pomeron: abbreviations
5426 c---------------------------------------
5427 if(iappl.eq.1.or.iappl.eq.8.or.iappl.eq.9)then
5430 c---------------------------------------
5431 c auxiliary constants:
5432 c---------------------------------------
5433 stmass=.05 !string mass cutoff
5435 c---------------------------------------
5436 c parton density normalization
5437 sq=log(log(q2min/.232**2)/log(.23/.232**2))
5438 du=2.997+.753*sq-.076*sq*sq
5443 xxq=1.-xx**(1./(1.+du))
5444 qnorm=qnorm+a1(i)*(psdfh4(xxq,q2min,0.,2,1)+
5445 * psdfh4(xxq,q2min,0.,2,2))/(1.-xxq)**du
5448 qnorm=qnorm*.5/(1.+du)
5450 ckkkkk-----------------------------
5451 c rr=(1.-qnorm)/4./pi/gamhad(2)
5452 c * *utgam1(2.+betpom-dels)/utgam1(1.-dels)
5453 c * /utgam1(1.+betpom)/utgam1(1.+alplea(2))/
5454 c * utgam1(2.-alppar)*utgam1(3.+alplea(2)-alppar)
5455 c ffrr=(1.-qnorm)/4./pi/gamhad(2)
5456 c * *utgam1(2.+betpom-dels)/utgam1(1.-dels)
5457 c * /utgam1(1.+betpom)
5458 c write(6,*)'===========',ffrr
5461 * /utgam1(1.+alplea(2))/
5462 * utgam1(2.-alppar)*utgam1(3.+alplea(2)-alppar)
5464 ckkkkkkk-------------------------------
5465 if(ish.ge.4)write (ifch,*)'rr,qnorm',rr,qnorm
5468 sq=log(log(q2min/.232**2)/log(.25/.232**2))
5474 xxq=1.-xx**(1./(1.+dpi))
5475 qnorm=qnorm+a1(i)*(psdfh4(xxq,q2min,0.,1,1)+
5476 * psdfh4(xxq,q2min,0.,1,2))/(1.-xxq)**dpi
5479 qnorm=qnorm*.5/(1.+dpi)
5480 cftmp=1./(1.-qnormp)*(1.-qnorm)
5481 * *utgam1(alplea(2)+1.)/utgam1(alplea(2)+3.-alppar)
5482 * /utgam1(alplea(1)+1.)*utgam1(alplea(1)+3.-alppar)
5483 gamhad(1)=gamhad(2)*cftmp
5484 if(gamhadsi(1).lt.0.)then
5485 gamhads(1)=gamhad(1)
5487 gamhads(1)=gamhadsi(1)
5490 * write (ifch,*)'gamhad(1),gamhads(1)',gamhad(1),gamhads(1)
5492 if(gamhadsi(2).lt.0.)then
5493 gamhads(2)=gamhad(2)
5495 gamhads(2)=gamhadsi(2)
5498 * write (ifch,*)'gamhad(2),gamhads(2)',gamhad(2),gamhads(2)
5504 xxq=1.-xx**(1./(1.+dpi))
5505 qnorm=qnorm+a1(i)*(psdfh4(xxq,q2min,0.,1,1)+
5506 * psdfh4(xxq,q2min,0.,1,2))/(1.-xxq)**dpi
5509 qnorm=qnorm*.5/(1.+dpi)
5510 cftmp=1./(1.-qnormp)*(1.-qnorm)
5511 * *utgam1(alplea(2)+1.)/utgam1(alplea(2)+3.-alppar)
5512 * /utgam1(alplea(3)+1.)*utgam1(alplea(3)+3.-alppar)
5513 gamhad(3)=gamhad(2)*cftmp
5514 if(gamhadsi(3).lt.0.)then
5515 gamhads(3)=gamhad(3)
5517 gamhads(3)=gamhadsi(3)
5520 * write (ifch,*)'gamhad(3),gamhads(3)',gamhad(3),gamhads(3)
5523 gamhad(4)=gamhad(1)*(quamas/qcmass)**2
5524 if(gamhadsi(4).lt.0.)then
5525 gamhads(4)=gamhad(4)
5527 gamhads(4)=gamhadsi(4)
5530 * write (ifch,*)'gamhad(4),gamhads(4)',gamhad(4),gamhads(4)
5535 xxg=xx**(1./(1.-dels))
5536 gnorm=gnorm+a1(i)*(fzeroGluZZ(xxg,4)+fzeroSeaZZ(xxg,4))
5539 gnorm=gnorm/(1.-dels)*2.*pi*gamhad(4)*rr
5540 alvc=6./(1.-gnorm)-4.
5541 if(ish.ge.4) write (ifch,*)'rr,qnorm,gnorm,alvc',
5542 * rr,qnorm,gnorm,alvc
5544 c write (*,*)'rr-c,qnorm,gnorm,alvc',rr,qnorm,gnorm,alvc
5547 c-----------------------------------------------
5548 c tabulation of inclusive jet cross sections
5549 c--------------------------------------------------
5552 qi=q2min*exp(.5*(i-1))
5553 sudx(i,1)=psudx(qi,1)
5554 sudx(i,2)=psudx(qi,2)
5556 if(ish.ge.4)write(ifch,*)'bare cross sections ...'
5560 ccc call MakeCSTable
5562 inquire(file=fnii(1:nfnii),exist=lcalc)
5565 write(ifmt,'(3a)')'read from ',fnii(1:nfnii),' ...'
5566 open(1,file=fnii(1:nfnii),status='old')
5567 read (1,*)qcdlam0,q2min0,q2ini0,naflav0,epmax0,pt2cut0
5568 if(qcdlam0.ne.qcdlam)write(ifmt,'(a)')'initl: wrong qcdlam'
5569 if(q2min0 .ne.q2min )write(ifmt,'(a)')'initl: wrong q2min'
5570 if(q2ini0 .ne.q2ini )write(ifmt,'(a)')'initl: wrong q2ini'
5571 if(naflav0.ne.naflav)write(ifmt,'(a)')'initl: wrong naflav'
5572 if(epmax0 .ne.epmax )write(ifmt,'(a)')'initl: wrong epmax'
5573 if(pt2cut0 .ne.pt2cut )write(ifmt,'(a)')'initl: wrong pt2cut'
5574 if(qcdlam0.ne.qcdlam.or.q2min0 .ne.q2min .or.q2ini0 .ne.q2ini
5575 * .or.naflav0.ne.naflav.or.epmax0 .ne.epmax.or. pt2cut.ne.pt2cut0)
5577 write(ifmt,'(//a//)')' initl has to be reinitialized!!!'
5580 read (1,*)csbor,csord,cstot,cstotzero,csborzer
5587 write(ifmt,'(a)')'initl does not exist -> calculate tables ...'
5589 write (*,*)'Born xsection csbor'
5591 spminc=4.*q2min+qcmass**2
5592 do m=1,4 !parton type at upper end of the ladder (1...4 - g,u,d,c)
5595 sk=spmin*(epmax/2./spmin)**((k-1)/19.)
5598 sk=spminc*(epmax/2./spminc)**((k-1)/19.)
5599 p1=sk/(1.+qcmass**2/sk)
5603 qq=q2min*(qmax/q2min)**((i-1)/19.)
5604 do l=1,2 !parton type at lower end of the ladder
5605 k1=k+20*(m-1)+80*(l-1)
5607 if(m.eq.3.and.l.eq.1)then !dd~
5611 endif !born cr.-sect.
5612 csbor(i,k1,1)=log(max(1.e-30,psborn(qq,qq,qq,sk,m1,l1,0,0)))
5614 csbor(i,k1,2)=log(max(1.e-30,psborn(4.*qq,qq,qq,sk,m1,l1,1,0)))
5621 write (*,*)'ordered jet xsection csord'
5622 do m=1,4 !parton type at upper end of the ladder
5624 write (*,*)' m=',m,'/4 k=',k,'/20'
5626 sk=spmin*(epmax/2./spmin)**((k-1)/19.) !c.m. energy squared for the hard
5629 sk=spminc*(epmax/2./spminc)**((k-1)/19.)
5630 p1=sk/(1.+qcmass**2/sk)
5634 do i=1,20 !cross-sections initialization
5635 qi=q2min*(qmax/q2min)**((i-1)/19.)
5637 qq=qi*(qmax/qi)**((j-1)/19.)
5639 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
5643 do l=1,2 !parton type at lower end of the ladder
5645 if(m.eq.3.and.l.eq.1)then
5651 k1=k+20*(m-1)+60*(l-1)
5652 if(k.eq.1.or.i.eq.20.or.j.eq.20)then
5653 csord(i,j,k1)=log(max(1.e-30,psborn(qi,qq,qq,sk,m1,l1,0,0)))
5655 * log(max(1.e-30,psborn(4.*qq,qi,qq,sk,l1,m1,1,0)))
5657 csord(i,j,k1)=log(psjet1(qi,qq,qq,sk,m1,l1,0)
5658 * /(1./tmin-1./tmax)+psborn(qi,qq,qq,sk,m1,l1,0,0))
5659 csord(i,j,k1+120)=log(psjet1(qi,4.*qq,qq,sk,m1,l1,2)
5660 * /(1./tmin-1./tmax)+psborn(4.*qq,qi,qq,sk,l1,m1,1,0))
5664 if(k.eq.1.or.i.eq.20)then
5665 cschar(i,k,l)=log(max(1.e-30,psborn(q2min,qi,qq,sk,m1,l1,0,0)))
5667 cschar(i,k,l)=log(psjet1(qi,q2min,qq,sk,l1,m1,0)
5668 * /(1./tmin-1./tmax)+psborn(q2min,qi,qq,sk,m1,l1,0,0))
5677 write (ifmt,*)'tests:'
5678 write (ifmt,'(a,a)')' n-1 sk qi qj qq '
5679 * ,' born born-i ord ord-i '
5681 sk=spmin*(epmax/2./spmin)**((k-1)/19.)
5682 if(k.ge.5)sk=spmin*1.5**(k-4)
5692 qi=q2min*(qmax1/q2min)**((i-1)/3.)
5694 qj=q2min*(qmax2/q2min)**((j-1)/3.)
5702 qq=qqmin*(qqmax/qqmin)**((lq-1)/3.)
5704 tmin=2.*qq/(1.+sqrt(1.-4.*qq/sk))
5709 do m=1,1 !parton type at upper end of the ladder (1
5710 do l=1,1 !parton type at lower end of the ladder (1
5712 if(m.eq.3.and.l.eq.1)then
5717 a=psborn(qj,qi,qq,sk,l1,m1,n-1,0)*(1./tmin-1./tmax)
5718 b=psbint(qj,qi,qq,sk,l1,m1,n-1)
5719 c=psjet1(qi,qj,qq,sk,m1,l1,2*(n-1))
5720 * +psborn(qj,qi,qq,sk,l1,m1,n-1,0)*(1./tmin-1./tmax)
5721 d=psjti1(qi,qj,qq,sk,m1,l1,n-1)
5722 write (ifmt,'(i3,4f9.1,3x,4f9.4)')n-1,sk,qi,qj,qq,a,b,c,d
5731 write (*,*)'jet xsection cstot'
5733 write (*,*)'k=',k,'/20'
5734 sk=spmin*(epmax/2./spmin)**((k-1)/19.) !c.m. energy squared for the hard
5737 do i=1,20 !cross-sections initialization
5740 qi=q2min*(qmax/q2min)**((i-1)/19.)
5742 qi=q2min*(4.*qmax/q2min)**((i-1)/19.)
5746 qq=qi*(qmax/qi)**((j-1)/19.)
5748 qq=max(q2min,qi/4.)*(qmax/max(q2min,qi/4.))**
5752 tmin=2.*qq/(1.+sqrt(1.-4.*qq/sk))
5756 do m=1,3 !parton type at upper end of the ladder (1
5757 do l=1,2 !parton type at lower end of the ladder (1
5759 if(m.eq.3.and.l.eq.1)then
5764 k1=k+20*(m-1)+60*(l-1)+120*(n-1)
5765 if(k.eq.1.or.i.eq.20.or.j.eq.20)then
5766 cstot(i,j,k1)=log(max(1.e-30,psborn(qi,q2min,qq,sk,m1,l1,n-1,0)))
5769 cstot(i,j,k1)=log((psjet(qi,q2min,qq,sk,m1,l1,0)+
5770 * psjti1(qi,q2min,qq,sk,m1,l1,0)+
5771 * psjti1(q2min,qi,qq,sk,l1,m1,0)
5772 * -psbint(qi,q2min,qq,sk,m1,l1,0))/(1./tmin-1./tmax))
5774 cstot(i,j,k1)=log((psjet(qi,q2min,qq,sk,m1,l1,1)+
5775 * psjet1(qi,q2min,qq,sk,m1,l1,1)+
5776 * psjti1(q2min,qi,qq,sk,l1,m1,1))/(1./tmin-1./tmax))
5786 c total and born hard cross-sections logarithms for minimal cutoff
5787 c (q2min), interpolated in the psjti0 procedure
5789 spminc=4.*q2min+qcmass**2
5793 if(m.eq.3.and.l.eq.1)then
5800 sk=spmin*(epmax/2./spmin)**((k-1)/19.) !c.m. energy squared for the hard
5804 sk=spminc*(epmax/2./spminc)**((k-1)/19.)
5805 p1=sk/(1.+qcmass**2/sk)
5809 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
5815 k1=k+20*(m-1)+80*(l-1)
5817 * =log(max(1.e-30,psborn(q2min,q2min,qq,sk,m1,l1,0,0)))
5819 cstotzero(k,m,l)=csborzer(k,m,l)
5821 cstotzero(k,m,l)=log(psjti(q2min,qq,sk,m1,l1,0)/
5822 * (1./tmin-1./tmax))
5824 smins=2.5*q2min*(1.+sqrt(1.+4.*qcmass**2/q2min))
5826 cstotzero(k,m,l)=log(psjci(q2min,sk,l1)/(1./tmin-1./tmax))
5828 cstotzero(k,m,l)=log((psjci(q2min,sk,l1)+psjct(sk,l1))
5829 * /(1./tmin-1./tmax))
5836 write(ifmt,'(a)')'write to initl ...'
5837 open(1,file=fnii(1:nfnii),status='unknown')
5838 write (1,*)qcdlam,q2min,q2ini,naflav,epmax,pt2cut
5839 write (1,*)csbor,csord,cstot,cstotzero,csborzer,cschar
5844 if(iappl.ne.8)goto 3
5845 if(ish.ge.3)write(ifch,*)'dis cross sections ...'
5846 inquire(file=fnid(1:nfnid),exist=lcalc)
5849 write(ifmt,'(3a)')'read from ',fnid(1:nfnid),' ...'
5850 open(1,file=fnid(1:nfnid),status='old')
5851 read (1,*)qcdlam0,q2min0,q2ini0,naflav0,epmax0,edmax0
5852 if(qcdlam0.ne.qcdlam)write(ifmt,'(a)')'inidi: wrong qcdlam'
5853 if(q2min0 .ne.q2min )write(ifmt,'(a)')'inidi: wrong q2min'
5854 if(q2ini0 .ne.q2ini )write(ifmt,'(a)')'inidi: wrong q2ini'
5855 if(naflav0.ne.naflav)write(ifmt,'(a)')'inidi: wrong naflav'
5856 if(epmax0 .ne.epmax )write(ifmt,'(a)')'inidi: wrong epmax'
5857 if(edmax0 .ne.edmax )write(ifmt,'(a)')'inidi: wrong edmax'
5858 if(qcdlam0.ne.qcdlam.or.q2min0 .ne.q2min.or.q2ini0 .ne.q2ini
5859 * .or.naflav0.ne.naflav.or.epmax0 .ne.epmax
5860 * .or.edmax0 .ne.edmax)then
5861 write(ifmt,'(//a//)')' inidi has to be reinitialized!!!'
5864 read (1,*)csdsi,csds,csdt,csdr
5870 write(ifmt,'(a)')'inidi does not exist -> calculate tables ...'
5872 qq=q2min*exp(.5*(j-1)) !photon virtuality
5874 do m=1,2 !parton type at the end of the ladder
5876 s2min=4.*max(q2mass,q2min)+qq
5877 if(m.eq.2)s2min=s2min/(1.-4.*q2ini/(s2min-qq))
5879 write (*,*)'sin,j,m,k',j,m,k
5880 sk=s2min*(edmax/s2min)**(.04*(k-1)) !c.m. energy squared
5881 if(k.eq.26)sk=1.01*sk
5886 qmax=(sk-qq+sqrt((sk-qq)**2-16.*sk*q2ini))/8.
5889 do i=1,21 !cross-sections calculation
5890 qi=qmin*(qmax/qmin)**((i-1)/20.)
5892 qtq=4.*max(q2mass,qi)/(sk-qq)
5894 tmin=.5*sk*qtq/(1.+sqrt(1.-qtq))
5900 k1=k+26*(m-1)+52*(ilong-1)
5902 if(tmax.gt.1.01*tmin)then
5903 sij=psds(qi,qq,sk,m-1,ilong-1)
5904 if(sij.lt.0.)write (*,*)'qi,qq,sk,m,long,sij',
5905 * qi,qq,sk,m,ilong,sij
5906 csdsi(i,j,k1)=log(max(0.,sij)/(1./tmin-1./tmax)
5907 * +psdbor(qi,qq,sk,ilong-1))
5910 * log(max(1.e-25,psdbor(qi,qq,sk,ilong-1)))
5913 csdsi(i,j,k1)=psds(qi,qq,sk,m-1,ilong-1)
5923 qq=q2min*exp(.5*(j-1)) !photon virtuality
5924 s2min=max(4.*qq,16.*q2min) !pt2dis=qq
5928 k1=k+26*(m-1)+52*(ilong-1)
5929 csds(j,k,m+2*(ilong-1))=csdsi(1,j,k1)
5932 sk=(s2min+qq)*(edmax/(s2min+qq))**(.04*(k-1))
5933 csdt(j,k,m)=psdres(qq,sk,s2min,m-1)
5934 csdr(j,k,m)=psdrga(qq,sk-qq,s2min,m-1)
5939 write(ifmt,'(a)')'write to inidi ...'
5941 write(ifmt,'(a)')'write to inidi ...'
5942 open(1,file=fnid(1:nfnid),status='unknown')
5943 write (1,*)qcdlam,q2min,q2ini,naflav,epmax,edmax
5944 write (1,*)csdsi,csds,csdt,csdr
5948 c---------------------------------------
5949 c tabulation of semihard eikonals
5950 c---------------------------------------
5952 !!!!!!!!! if(iappl.eq.1)then
5954 if(ish.ge.4)write(ifch,*)'semihard eikonals ...'
5956 inquire(file=fnrj,exist=lcalc)
5959 write(ifmt,'(3a)')'read from ',fnrj(1:nfnrj),' ...'
5960 open(1,file=fnrj(1:nfnrj),status='old')
5961 read (1,*)alpqua0,alplea0,alppom0,slopom0,
5962 * gamhad0,r2had0,chad0,
5963 * qcdlam0,q2min0,q2ini0,betpom0,glusea0,naflav0,
5964 * factk0,pt2cut0,gamtil0
5965 if(alpqua0.ne.alpqua)write(ifmt,'(a,2f8.4)')
5966 * 'inirj: wrong alpqua',alpqua0,alpqua
5967 if(alppom0.ne.alppom)write(ifmt,'(a,2f8.4)')
5968 * 'inirj: wrong alppom',alppom0,alppom
5969 if(slopom0.ne.slopom)write(ifmt,'(a,2f8.4)')
5970 * 'inirj: wrong slopom',slopom0,slopom
5972 if(gamhad0(iii).ne.gamhad(iii))write(ifmt,'(a,i1,a,2f8.4)')
5973 * 'inirj: wrong gamhad(',iii,')',gamhad0(iii),gamhad(iii)
5975 if(r2had0(iii) .ne.r2had(iii) )write(ifmt,'(a,i1,a,2f8.4)')
5976 * 'inirj: wrong r2had(',iii,')',r2had0(iii),r2had(iii)
5977 if(chad0(iii) .ne.chad(iii) )write(ifmt,'(a,i1,a,2f8.4)')
5978 * 'inirj: wrong chad(',iii,')',chad0(iii),chad(iii)
5979 if(alplea0(iii).ne.alplea0(iii))write(ifmt,'(a,i1,a,2f8.4)')
5980 * 'inirj: wrong alplea(',iii,')',alplea0(iii),alplea(iii)
5982 if(qcdlam0.ne.qcdlam)write(ifmt,'(a,2f8.4)')
5983 * 'inirj: wrong qcdlam',qcdlam0,qcdlam
5984 if(q2min0 .ne.q2min )write(ifmt,'(a,2f8.4)')
5985 * 'inirj: wrong q2min',q2min0,q2min
5986 if(q2ini0 .ne.q2ini )write(ifmt,'(a,2f8.4)')
5987 * 'inirj: wrong q2ini',q2ini0,q2ini
5988 if(betpom0.ne.betpom)write(ifmt,'(a,2f8.4)')
5989 * 'inirj: wrong betpom',betpom0,betpom
5990 if(glusea0.ne.glusea)write(ifmt,'(a,2f8.4)')
5991 * 'inirj: wrong glusea',glusea0,glusea
5992 if(naflav0.ne.naflav)write(ifmt,'(a,2f8.4)')
5993 * 'inirj: wrong naflav',naflav0,naflav
5994 if(factk0 .ne.factk )write(ifmt,'(a,2f8.4)')
5995 * 'inirj: wrong factk', factk0,factk
5996 if(pt2cut0 .ne.pt2cut )write(ifmt,'(a,2f8.4)')
5997 * 'inirj: wrong pt2cut', pt2cut0,pt2cut
5998 if(gamtil0 .ne.gamtil )write(ifmt,'(a,2f8.4)')
5999 * 'inirj: wrong gamtil', gamtil0,gamtil
6000 if(alpqua0.ne.alpqua.or.alppom0.ne.alppom
6001 * .or.slopom0.ne.slopom.or.gamhad0(2).ne.gamhad(2)
6002 * .or.r2had0(1).ne.r2had(1).or.r2had0(2).ne.r2had(2)
6003 * .or.r2had0(3).ne.r2had(3)
6004 * .or.chad0(1).ne.chad(1).or.chad0(2).ne.chad(2)
6005 * .or.chad0(3).ne.chad(3)
6006 * .or.alplea0(1).ne.alplea(1).or.alplea0(2).ne.alplea(2)
6007 * .or.alplea0(3).ne.alplea(3)
6008 * .or.qcdlam0.ne.qcdlam.or.q2min0 .ne.q2min
6009 * .or.q2ini0 .ne.q2ini.or.gamtil0.ne.gamtil
6010 * .or.betpom0.ne.betpom.or.glusea0.ne.glusea.or.naflav0.ne.naflav
6011 * .or.factk0 .ne.factk .or.pt2cut0.ne.pt2cut)then
6012 write(ifmt,'(//a//)')' inirj has to be reinitialized!!!!'
6016 read(1,*)fhgg,fhqg,fhgq,fhqq,fhgg0,fhgg1,fhqg1
6017 * ,fhgg01,fhgg02,fhgg11,fhgg12,fhqg11,fhqg12
6018 * ,ftoint,vfro,vnorm,coefxu1,coefxu2,coefxc2
6019 read(1,*)bkbin0,iclpro10,iclpro20,icltar10,icltar20,iclegy10
6020 * ,iclegy20,egylow0,egymax0,iomega0,egyscr0,epscrw0,epscrp0
6023 if(iclpro10.ne.iclpro1)write(textini,'(a,2i8)')
6024 * 'inirj: wrong iclpro1 ',iclpro10,iclpro1
6025 if(iclpro20.ne.iclpro2)write(textini,'(a,2i8)')
6026 * 'inirj: wrong iclpro2 ',iclpro20,iclpro2
6027 if(icltar10.ne.icltar1)write(textini,'(a,2i8)')
6028 * 'inirj: wrong icltar1 ',icltar10,icltar1
6029 if(icltar20.ne.icltar2)write(textini,'(a,2i8)')
6030 * 'inirj: wrong icltar2 ',icltar20,icltar2
6031 if(iclegy10.ne.iclegy1)write(textini,'(a,2i8)')
6032 * 'inirj: wrong iclegy1 ',iclegy10,iclegy1
6033 if(iclegy20.ne.iclegy2)write(textini,'(a,2i8)')
6034 * 'inirj: wrong iclegy2 ',iclegy20,iclegy2
6035 if(egylow0.ne.egylow)write(textini,'(a,2f8.4)')
6036 * 'inirj: wrong egylow ',egylow0,egylow
6037 if(egymax0.ne.egymax)write(textini,'(a,2f8.4)')
6038 * 'inirj: wrong egymax ',egymax0,egymax
6039 if(epscrw0.ne.epscrw)write(textini,'(a,2f8.4)')
6040 * 'inirj: wrong epscrw ',epscrw0,epscrw
6041 if(epscrp0.ne.epscrp)write(textini,'(a,2f8.4)')
6042 * 'inirj: wrong epscrp ',epscrp0,epscrp
6043 if(bkbin0.ne.bkbin)write(textini,'(a,2f8.4)')
6044 * 'inirj: wrong bkbin',bkbin0,bkbin
6045 if(textini.ne.' ')then
6046 write(ifmt,'(//10x,a//10x,a//)')textini,
6047 * 'inirj has to be reinitialized!!!!'
6050 do iiipro=iclpro1,iclpro2
6051 do iiitar=icltar1,icltar2
6052 do iiiegy=iclegy1,iclegy2
6054 read(1,*)xkappafit(iiiegy,iiipro,iiitar,iiib)
6056 xkappafit(iiiegy,iiipro,iiitar,nbkbin)=1.
6058 if(xkappafit(iiiegy,iiipro,iiitar,iiib).lt.1.)then
6059 xkappafit(iiiegy,iiipro,iiitar,iiib)=max(1.,0.5*
6060 * (xkappafit(iiiegy,iiipro,iiitar,iiib-1)
6061 * +xkappafit(iiiegy,iiipro,iiitar,iiib+1)))
6065 read(1,*)alpDs(iiidf,iiiegy,iiipro,iiitar),
6066 * alpDps(iiidf,iiiegy,iiipro,iiitar),
6067 * alpDpps(iiidf,iiiegy,iiipro,iiitar),
6068 * betDs(iiidf,iiiegy,iiipro,iiitar),
6069 * betDps(iiidf,iiiegy,iiipro,iiitar),
6070 * betDpps(iiidf,iiiegy,iiipro,iiitar),
6071 * gamDs(iiidf,iiiegy,iiipro,iiitar),
6072 * delDs(iiidf,iiiegy,iiipro,iiitar)
6087 write(ifmt,'(a)')'inirj does not exist -> calculate tables ...'
6095 spminc=4.*q2min+2.*qcmass**2
6098 write(ifmt,'(a)')' tabulate om5 ...'
6101 sy=spmin*(epmax/2./spmin)**((iy-1)/10.)
6102 syc=spminc*(epmax/2./spminc)**((iy-1)/10.)
6109 fhgg01(iy)=log(om51pp(sy,1.,1.,3))
6110 fhgg02(iy)=log(om51pp(sy,1.,1.,7))
6113 do iclpro=iclpro1,iclpro2
6115 fhgg11(iy,iclpro)=-80.
6116 fhgg12(iy,iclpro)=-80.
6118 fhgg11(iy,iclpro)=log(om51pp(sy,1.,1.,4))
6119 fhgg12(iy,iclpro)=log(om51pp(sy,1.,1.,9))
6128 fhqg11(iy,ix,iclpro)=-80.
6129 fhqg12(iy,ix,iclpro)=-80.
6130 elseif(iclpro.eq.4)then
6131 fhqg11(iy,ix,iclpro)=log(om51pp(syc,1.,1.,5))
6132 fhqg12(iy,ix,iclpro)=log(om51pp(syc,1.,1.,11))
6134 fhqg11(iy,ix,iclpro)=log(om51pp(sy,xp,1.,5))
6135 fhqg12(iy,ix,iclpro)=log(om51pp(sy,xp,1.,11))
6148 fhgg0(iy,iz)=log(om51pp(sy,1.,z,6)/z)
6151 do iclpro=iclpro1,iclpro2
6153 fhgg1(iy,iz,iclpro)=-80.
6155 fhgg1(iy,iz,iclpro)=log(om51pp(sy,1.,z,8)/z)
6165 fhqg1(iy,ix,iz+10*(iclpro-1))=-80.
6166 elseif(iclpro.eq.4)then
6167 fhqg1(iy,ix,iz+10*(iclpro-1))=log(om51pp(syc,xp,z,10)/z)
6169 fhqg1(iy,ix,iz+10*(iclpro-1))=log(om51pp(sy,xp,z,10)/z)
6176 do iclpro=iclpro1,iclpro2 !hadron type (1 - pion, 2 - nucleon, 3 - kaon, 4 - charm)
6177 do icltar=icltar1,icltar2 !hadron type (2 - nucleon)
6179 sy=spmin*(epmax/2./spmin)**((iy-1)/10.)
6180 syc=spminc*(epmax/2./spminc)**((iy-1)/10.)
6184 fhgg(iy,iz,iclpro+4*(icltar-1))=-80.
6186 fhgg(iy,iz,iclpro+4*(icltar-1))=log(om51pp(sy,1.,z,0)/z)
6196 fhqg(iy,ix,iz+10*(iclpro+4*(icltar-1)-1))=-80.
6197 fhgq(iy,ix,iz+10*(iclpro+4*(icltar-1)-1))=-80.
6204 fhqg(iy,ix,iz+10*(iclpro+4*(icltar-1)-1))=
6205 * log(om51pp(syx,xp,z,1)/z)
6211 fhgq(iy,ix,iz+10*(iclpro+4*(icltar-1)-1))=
6212 * log(om51pp(syx,xp,z,2)/z)
6231 fhqq(iy,ix1,ix2+10*(iclpro+4*(icltar-1)-1))=-80.
6233 if(iclpro.ne.4.and.icltar.ne.4)then
6238 fhqq(iy,ix1,ix2+10*(iclpro+4*(icltar-1)-1))=
6239 * log(pshard(syx,xpp,xmm))
6251 write(ifmt,'(a)')' tabulate fit parameters ...'
6254 do iclpro=iclpro1,iclpro2 !hadron type (1 - pion, 2 - nucleon, 3 - kaon, 4 - charm)
6255 do icltar=icltar1,icltar2 !hadron type (2 - nucleon)
6256 do iclegy=iclegy2,iclegy1,-1
6259 do iiclegy=iclegy2,iclegy1,-1
6260 engy=egyfac**(iiclegy-1)*egylow
6270 write(ifmt,'(a)')' write to inirj ...'
6271 open(1,file=fnrj,status='unknown')
6272 write (1,*)alpqua,alplea,alppom,slopom,gamhad,r2had,chad,
6273 *qcdlam,q2min,q2ini,betpom,glusea,naflav,factk,pt2cut,gamtil
6274 write (1,*)fhgg,fhqg,fhgq,fhqq,fhgg0,fhgg1,fhqg1
6275 *,fhgg01,fhgg02,fhgg11,fhgg12,fhqg11,fhqg12
6276 *,ftoint,vfro,vnorm,coefxu1,coefxu2,coefxc2
6277 write(1,*)bkbin,iclpro1,iclpro2,icltar1,icltar2,iclegy1,iclegy2
6278 *,egylow,egymax,iomega,egyscr,epscrw,epscrp
6279 do iiipro=iclpro1,iclpro2
6280 do iiitar=icltar1,icltar2
6281 do iiiegy=iclegy1,iclegy2
6283 write(1,*)xkappafit(iiiegy,iiipro,iiitar,iiib)
6286 write(1,*)alpDs(iiidf,iiiegy,iiipro,iiitar),
6287 * alpDps(iiidf,iiiegy,iiipro,iiitar),
6288 * alpDpps(iiidf,iiiegy,iiipro,iiitar),
6289 * betDs(iiidf,iiiegy,iiipro,iiitar),
6290 * betDps(iiidf,iiiegy,iiipro,iiitar),
6291 * betDpps(iiidf,iiiegy,iiipro,iiitar),
6292 * gamDs(iiidf,iiiegy,iiipro,iiitar),
6293 * delDs(iiidf,iiiegy,iiipro,iiitar)
6311 c--------------------------------------
6312 c inelastic cross sections
6313 c---------------------------------------
6315 if(isetcs.ge.2)then !--------------------
6317 if(ish.ge.4)write(ifch,*)'cross sections ...'
6319 inquire(file=fncs,exist=lcalc)
6322 write(ifmt,'(3a)')'read from ',fncs(1:nfncs),' ...'
6323 open(1,file=fncs(1:nfncs),status='old')
6324 read (1,*)alpqua0,alplea0,alppom0,slopom0,
6325 * gamhad0,r2had0,chad0,
6326 * qcdlam0,q2min0,q2ini0,betpom0,glusea0,naflav0,
6328 if(alpqua0.ne.alpqua)write(ifmt,'(a,2f8.4)')
6329 * 'inics: wrong alpqua',alpqua0,alpqua
6330 if(alppom0.ne.alppom)write(ifmt,'(a,2f8.4)')
6331 * 'inics: wrong alppom',alppom0,alppom
6332 if(slopom0.ne.slopom)write(ifmt,'(a,2f8.4)')
6333 * 'inics: wrong slopom',slopom0,slopom
6335 if(gamhad0(iii).ne.gamhad(iii))write(ifmt,'(a,i1,a,2f8.4)')
6336 * 'inics: wrong gamhad(',iii,')',gamhad0(iii),gamhad(iii)
6338 if(r2had0(iii) .ne.r2had(iii) )write(ifmt,'(a,i1,a,2f8.4)')
6339 * 'inics: wrong r2had(',iii,')',r2had0(iii),r2had(iii)
6340 if(chad0(iii) .ne.chad(iii) )write(ifmt,'(a,i1,a,2f8.4)')
6341 * 'inics: wrong chad(',iii,')',chad0(iii),chad(iii)
6342 if(alplea0(iii).ne.alplea0(iii))write(ifmt,'(a,i1,a,2f8.4)')
6343 * 'inics: wrong alplea(',iii,')',alplea0(iii),alplea(iii)
6345 if(qcdlam0.ne.qcdlam)write(ifmt,'(a,2f8.4)')
6346 * 'inics: wrong qcdlam',qcdlam0,qcdlam
6347 if(q2min0 .ne.q2min )write(ifmt,'(a,2f8.4)')
6348 * 'inics: wrong q2min',q2min0,q2min
6349 if(q2ini0 .ne.q2ini )write(ifmt,'(a,2f8.4)')
6350 * 'inics: wrong q2ini',q2ini0,q2ini
6351 if(betpom0.ne.betpom)write(ifmt,'(a,2f8.4)')
6352 * 'inics: wrong betpom',betpom0,betpom
6353 if(glusea0.ne.glusea)write(ifmt,'(a,2f8.4)')
6354 * 'inics: wrong glusea',glusea0,glusea
6355 if(naflav0.ne.naflav)write(ifmt,'(a,2f8.4)')
6356 * 'inics: wrong naflav',naflav0,naflav
6357 if(factk0 .ne.factk )write(ifmt,'(a,2f8.4)')
6358 * 'inics: wrong factk', factk0,factk
6359 if(pt2cut0 .ne.pt2cut )write(ifmt,'(a,2f8.4)')
6360 * 'inics: wrong pt2cut', pt2cut0,pt2cut
6361 if(alpqua0.ne.alpqua.or.alppom0.ne.alppom
6362 * .or.slopom0.ne.slopom.or.gamhad0(2).ne.gamhad(2)
6363 * .or.r2had0(1).ne.r2had(1).or.r2had0(2).ne.r2had(2)
6364 * .or.r2had0(3).ne.r2had(3)
6365 * .or.chad0(1).ne.chad(1).or.chad0(2).ne.chad(2)
6366 * .or.chad0(3).ne.chad(3)
6367 * .or.alplea0(1).ne.alplea(1).or.alplea0(2).ne.alplea(2)
6368 * .or.alplea0(3).ne.alplea(3)
6369 * .or.qcdlam0.ne.qcdlam.or.q2min0 .ne.q2min
6370 * .or.q2ini0 .ne.q2ini
6371 * .or.betpom0.ne.betpom.or.glusea0.ne.glusea.or.naflav0.ne.naflav
6372 * .or.factk0 .ne.factk .or.pt2cut0.ne.pt2cut)then
6373 write(ifmt,'(//a//)')' inics has to be reinitialized!!!!'
6377 read(1,*)isetcs0,iclpro10,iclpro20,icltar10,icltar20,iclegy10
6378 * ,iclegy20,egylow0,egymax0,iomega0,egyscr0,epscrw0,epscrp0
6380 if(iclpro10.ne.iclpro1)write(ifmt,'(a,2i2)')
6381 * 'inics: wrong iclpro1',iclpro10,iclpro1
6382 if(iclpro20.ne.iclpro2)write(ifmt,'(a,2i2)')
6383 * 'inics: wrong iclpro2',iclpro20,iclpro2
6384 if(icltar10.ne.icltar1)write(ifmt,'(a,2i2)')
6385 * 'inics: wrong icltar1',icltar10,icltar1
6386 if(icltar20.ne.icltar2)write(ifmt,'(a,2i2)')
6387 * 'inics: wrong icltar2',icltar20,icltar2
6388 if(iclegy10.ne.iclegy1)write(ifmt,'(a,2i4)')
6389 * 'inics: wrong iclegy1',iclegy10,iclegy1
6390 if(iclegy20.ne.iclegy2)write(ifmt,'(a,2i4)')
6391 * 'inics: wrong iclegy2',iclegy20,iclegy2
6392 if(egylow0.ne.egylow)write(ifmt,'(a,2f8.4)')
6393 * 'inics: wrong egylow',egylow0,egylow
6394 if(egymax0.ne.egymax)write(ifmt,'(a,2f12.4)')
6395 * 'inics: wrong egymax',egymax0,egymax
6396 if(egyscr0.ne.egyscr)write(ifmt,'(a,2f8.4)')
6397 * 'inics: wrong egyscr ',egyscr0,egyscr
6398 if(epscrw0.ne.epscrw)write(ifmt,'(a,2f8.4)')
6399 * 'inics: wrong epscrw',epscrw0,epscrw
6400 if(epscrp0.ne.epscrp)write(ifmt,'(a,2f8.4)')
6401 * 'inics: wrong epscrp',epscrp0,epscrp
6402 if(isetcs0.lt.isetcs)write(ifmt,'(a,2f8.4)')
6403 * 'inics: wrong isetcs',isetcs0,isetcs
6404 if(iclpro10.ne.iclpro1.or.iclpro20.ne.iclpro2
6405 * .or.icltar10.ne.icltar1.or.icltar20.ne.icltar2
6406 * .or.iclegy10.ne.iclegy1.or.iclegy20.ne.iclegy2
6407 * .or.egylow0.ne.egylow.or.egymax0.ne.egymax
6408 * .or.egyscr0.ne.egyscr.or.epscrw0.ne.epscrw.or.isetcs0.lt.isetcs
6409 * .or.epscrp0.ne.epscrp)then
6410 write(ifmt,'(//a//)')' inics has to be reinitialized!!!!'
6414 read (1,*)asect,asect2,asectn,asect4
6415 elseif(isetcs.eq.3)then
6416 read (1,*)asect1,asect,asect3,asectn
6418 write(ifmt,'(//a//)')' Wrong isetcs in psaini !!!!'
6432 idprojinsave=idprojin
6434 idtarginsave=idtargin
6466 fctrmx=100. !to get stable pA and AA cross section, this number has to be large
6467 ifrade=0 !to save time, no fragmentation
6469 write(ifmt,'(a)')'inics does not exist -> calculate tables ...'
6475 if(iclpro.lt.iclpro1.or.iclpro.gt.iclpro2)then
6478 asect1(ie,iclpro,iia)=0.
6479 asect2(ie,iclpro,iia)=0.
6484 engy=1.5*10.**(ie-1)
6486 write(ifmt,*)' calcul. ',ie,' (',iclpro,')',engy
6491 if(matarg.eq.1)then !hadron-proton interaction
6498 rrr=rad/difnuc(matarg)
6500 anorm=1.5/pi/rrr**3/(1.+(pi/rrr)**2)/difnuc(matarg)**2
6501 gin=(ptgau(ptfau,bm,2)+ptgau1(bm,2))*10. !sig_in
6503 if(ish.ge.3)write (ifch,226)gin
6504 226 format(2x,'psaini: hadron-nucleus cross sections:'/
6506 asect1(ie,iclpro,iia)=log(gin)
6507 asect(ie,iclpro,iia)=asect1(ie,iclpro,iia)
6508 write(ifmt,*)' matarg,gin :'
6513 if(isetcssave.ge.3)then
6517 elseif(iclpro.eq.2)then
6519 elseif(iclpro.eq.3)then
6523 engy=1.5*10.**(ie-1)
6524 if(engy.le.egymin)engy=egymin
6525 if(engy.ge.egymax)engy=egymax
6526 write(ifmt,*)' simul. ',ie,' (',iclpro,')',engy
6527 write(ifch,*)' simul. ',ie,' (',iclpro,')',engy
6530 latarg=min(1,matarg/2)
6545 if(ntry.lt.10000)then
6546 c if random sign for projectile, set it here
6547 idproj=idprojin*(1-2*int(rangen()+0.5d0))
6549 if(iret.gt.0)goto 222
6555 if(a.gt.0..and.ntevt.gt.0.)then
6556 xs=float(nevent)/float(ntevt)*a*10.
6557 write(ifmt,*)' matarg,nevent,ntevt,bmax,xs :'
6558 . ,matarg,nevent,ntevt,bmax,xs
6559 write(ifch,*)' matarg,nevent,ntevt,bmax,xs :'
6560 . ,matarg,nevent,ntevt,bmax,xs
6561 asect2(ie,iclpro,iia)=log(xs)
6563 write(ifmt,*)' Problem ? ',iclpro,matarg,bmax,ntevt
6564 asect2(ie,iclpro,iia)=0.
6571 asect2(ie,iclpro,iia)=0.
6582 engy=1.5*10.**(ie-1)
6584 write(ifmt,*)' calcul. AB ',ie,engy
6588 laproj=max(1,maproj/2)
6591 latarg=max(1,matarg/2)
6593 if(matarg.eq.1.and.maproj.eq.1)then !proton-proton interaction
6601 rrr=rad/difnuc(matarg)
6603 anorm=1.5/pi/rrr**3/(1.+(pi/rrr)**2)/difnuc(matarg)**2
6604 gin=(ptgau(ptfau,bm,2)+ptgau1(bm,2))*10. !sig_in
6605 elseif(matarg.eq.1)then
6608 rrrp=radp/difnuc(maproj)
6610 anormp=1.5/pi/rrrp**3/(1.+(pi/rrrp)**2)/difnuc(maproj)**2
6611 gin=(ptgau(ptfau,bm,1)+ptgau1(bm,1))*10. !sig_in
6613 rad=radnuc(matarg)+1.
6614 radp=radnuc(maproj)+1.
6616 rrr=rad/difnuc(matarg)
6618 rrrp=radp/difnuc(maproj)
6620 anorm=1.5/pi/rrr**3/(1.+(pi/rrr)**2)/difnuc(matarg)**2
6621 anormp=1.5/pi/rrrp**3/(1.+(pi/rrrp)**2)/difnuc(maproj)**2
6622 gin=(ptgau(ptfauAA,bm,2)+ptgau2(bm))*10.
6625 if(ish.ge.3)write (ifch,227)gin
6626 227 format(2x,'psaini: nucleus-nucleus cross sections:'/
6628 asect3(ie,iia,iib)=log(gin)
6629 asectn(ie,iia,iib)=asect3(ie,iia,iib)
6630 write(ifmt,*)' maproj,matarg,gin :'
6631 * ,maproj,matarg,gin
6637 if(isetcssave.ge.3)then
6640 engy=1.5*10.**(ie-1)
6641 if(engy.le.egymin)engy=egymin
6642 if(engy.ge.egymax)engy=egymax
6643 write(ifmt,*)' AB xs ',ie,engy
6644 write(ifch,*)' AB xs ',ie,engy
6647 laproj=max(1,maproj/2)
6650 latarg=max(1,matarg/2)
6666 if(ntry.lt.10000)then
6668 if(iret.gt.0)goto 223
6674 if(a.gt.0..and.ntevt.gt.0.)then
6675 xs=float(nevent)/float(ntevt)*a*10.
6676 write(ifmt,*)' maproj,matarg,nevent,ntevt,bmax,xs :'
6677 & ,maproj,matarg,nevent,ntevt,bmax,xs
6678 write(ifch,*)' maproj,matarg,nevent,ntevt,bmax,xs :'
6679 & ,maproj,matarg,nevent,ntevt,bmax,xs
6680 asect4(ie,iia,iib)=log(xs)
6682 write(ifmt,*)' Problem ? ',maproj,matarg,bmax,ntevt
6683 asect4(ie,iia,iib)=0.
6692 asect4(ie,iia,iib)=0.
6700 idprojin=idprojinsave
6702 idtargin=idtarginsave
6726 write(ifmt,'(a)')'write to inics ...'
6727 open(1,file=fncs,status='unknown')
6728 write (1,*)alpqua,alplea,alppom,slopom,gamhad,r2had,chad,
6729 *qcdlam,q2min,q2ini,betpom,glusea,naflav,factk,pt2cut
6730 write(1,*)isetcs,iclpro1,iclpro2,icltar1,icltar2,iclegy1,iclegy2
6731 *,egylow,egymax,iomega,egyscr,epscrw,epscrp
6732 write (1,*)asect1,asect2,asect3,asect4
6741 endif !----------isetcs.ge.2-----------
6743 call utprix('psaini',ish,ishini,4)
6748 cc-----------------------------------------------------------------------
6749 c function fjetxx(jpp,je1,je2)
6750 cc-----------------------------------------------------------------------
6751 cc almost exactly psjet, just with Eqcd replaced by fparton
6753 cc gives indeed the same result as jetx
6754 cc so the integration seems correct
6755 cc-----------------------------------------------------------------------
6756 c double precision xx1,xx2,s2min,xmin,xmax,xmin1,xmax1,t,tmin
6757 c *,tmax,sh,z,qtmin,ft,fx1,fx2
6758 c common /ar3/ x1(7),a1(7)
6759 c common /ar9/ x9(3),a9(3)
6760 c include 'epos.inc'
6761 c include 'epos.incsem'
6767 c zmin=s2min/dble(s)
6770 c zmin=zmin**(-delh)
6771 c zmax=zmax**(-delh)
6774 c z=dble(.5*(zmax+zmin+(zmin-zmax)*(2*m-3)*x9(i)))**(-1./delh)
6777 c qtmin=max(dble(q2min),dble(q2ini)/(1.d0-dsqrt(z)))
6778 c tmin=max(0.d0,1.d0-4.d0*qtmin/sh)
6779 c tmin=2.d0*qtmin/(1.d0+dsqrt(tmin))
6784 c t=2.d0*tmin/(1.d0+tmin/tmax-dble(x9(i1)*(2*m1-3))
6785 c & *(1.d0-tmin/tmax))
6787 c xmax=1.d0-q2ini/qt
6788 c xmin=max(dsqrt(z),z/xmax) !xm<xp !!!
6789 c if(xmin.gt.xmax.and.ish.ge.1)write(ifmt,*)'fjetxx:xmin,xmax'
6793 c if(xmax.gt..8d0)then
6794 c xmin1=max(xmin,.8d0)
6797 c xx1=1.d0-(1.d0-xmax)*((1.d0-xmin1)/(1.d0-xmax))**
6798 c * dble(.5+x9(i2)*(m2-1.5))
6800 c fb=ffsigj(sngl(t),qt,sngl(xx1),sngl(xx2),jpp,je1,je2)
6801 c * +ffsigj(sngl(t),qt,sngl(xx2),sngl(xx1),jpp,je1,je2)
6802 c fx1=fx1+dble(a9(i2)*fb)*(1.d0/xx1-1.d0)
6803 c * *pssalf(qt/qcdlam)**2
6806 c fx1=fx1*dlog((1.d0-xmin1)/(1.d0-xmax))
6808 c if(xmin.lt..8d0)then
6809 c xmax1=min(xmax,.8d0)
6812 c xx1=xmin*(xmax1/xmin)**dble(.5+x9(i2)*(m2-1.5))
6817 c * +ffsigj(sngl(t),qt,sngl(xx1),sngl(xx2),jpp,je1,je2)
6818 c * +ffsigj(sngl(t),qt,sngl(xx2),sngl(xx1),jpp,je1,je2)
6819 c fx2=fx2+dble(a9(i2))*fb*pssalf(qt/qcdlam)**2
6822 c fx2=fx2*dlog(xmax1/xmin)
6824 c ft=ft+dble(a9(i1))*(fx1+fx2)*t**2
6827 c ft=ft*(1.d0/tmin-1.d0/tmax)
6828 c fjetxx=fjetxx+a9(i)*sngl(ft*z**(1.+delh)/sh**2)
6829 c * /z ! ffsig = xp f xm f sigma
6832 c fjetxx=fjetxx*(zmin-zmax)/delh*pi**3
6833 c ! * /2. !??????????????? kkkkkkkkk