1 c----------------------------------------------------------------------
2 c The two elementary fuctions of our approach, the profile fuction G
3 c and the Phi exponent G~, are here referred to as Gpro and Gtilde.
4 c Both functions can be written as
6 c G = \sum_type alp * xp**bet * xm**betp
8 c The parameters alp, bet, betp are calculated in GfunParK (k-mode,
9 c b is takento the one of pair k) or GfunPar (b-mode: arbitrary b) as
11 c Gpro: bet = betD' + epsGp + gamD*b**2 + epsG -alppar
12 c bet' = betD'' + epsGt + gamD*b**2 + epsG -alppar
14 c alp = alpD*f * s**(betD+gamD*b**2+epsG) * exp(-b**2/delD)
16 c Gtilde: bet~ = bet + 1
19 c alp~ = alp * gam(bet~) * gam(bet~')
20 c * gam(1+alppro) * gam(1+alptar)
21 c * gam(1+alppro+bet~) * gam(1+alptar+bet~')
22 c * (1+epsGt') * (1+epsGt')
24 c The parameters depend implicitely on s.
26 c In the program we use om1 = Gpro
27 c (they differ by a constant which is actually one)
29 c All functions related to om1 are called om1... .
31 c The inclusive Pomeron distributions are
33 c PomInc(xp,xm) = Gpro(xp,xm) * (1-xp)**alppro * (1-xm)**alptar
35 c----------------------------------------------------------------------
38 c----------------------------------------------------------------------
39 subroutine GfunParK(irea) !---MC---
40 c----------------------------------------------------------------------
41 c calculates parameters alp,bet,betp of the G functions (k-mode)
42 c and the screening exponents epsilongp(k,i), epsilongt(k,i), epsilongs(k,i)
43 c----------------------------------------------------------------------
44 c Gpro parameters written to /comtilde/atilde(,)btildep(,),btildepp(,)
45 c Gtilde parameters written to /cgtilde/atildg(,)btildgp(,),btildgpp(,)
46 c two subscripts: first=type, second=collision k
47 c Certain pieces of this routine are only done if irea is <= or = zero.
48 c----------------------------------------------------------------------
53 double precision atildg,btildgp,btildgpp
54 common/cgtilde/atildg(idxD0:idxD1,kollmx)
55 *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx)
56 double precision utgam2,coefgdp,coefgdt
58 common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
59 common /cgtilnu/ cfbetpnp,cfbetppnp,cfbetpnm,cfbetppnm,cfalpro
60 &,cfaltar,cfbpap,cfbpam,cfbppap,cfbppam
61 double precision cfbetpnp,cfbetppnp,cfbetpnm,cfbetppnm,cfalpro
62 &,cfaltar,cfbpap,cfbpam,cfbppap,cfbppam,gamv,eps
63 parameter (eps=1.d-25)
64 parameter(nxeps=20,nyeps=32)
65 common/cxeps1/w(0:nxeps,nyeps),y1(nyeps),y2(nyeps)
66 common/cxeps2/db,b1,b2
67 common/geom/rmproj,rmtarg,bmax,bkmx
74 call utpri('GfunParK',ish,ishini,5)
103 ! calculate collision number according to Glauber -----------------------
105 bglaub2=max(1.,sigine/10./pi) !10= fm^2 -> mb
110 if(r.le.bglaub)nglevt=nglevt+1
113 ! Z_parton_projectile (zparpro) and Z_parton_target (zpartar)-----------
115 if(iscreen.eq.1.or.isplit.eq.1)then
118 ctp060829 bzero=bglaub
119 b2x=epscrp*epscrp*bglaub2
120 alpfom=alpfomi*epscrw*fscra(engy/egyscr)
127 !....... targ partons seen by proj
128 if(lproj(ip).ge.1)then
129 absb=max(1.e-9,bk(k))
131 zkp=epscrw*exp(-b2/2./b2x)*fscra(engy/egyscr)
132 c epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
133 c & epscrs*min(epscrx,zkp*abs(betD(0,iclpro,icltar))))
134 c epsilongp(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar,
135 c & epscrh*min(epscrx,zkp*abs(betD(1,iclpro,icltar))))
136 if(lproj3(ip).ge.1)then
141 absb=max(1.e-9,abs(bk(kp))-bcut)
143 zkp=epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr)
145 zpartar(k)=zpartar(k)+zkp
146 c alpfom=max(alpfom,dble(zpartar(k)))
153 !...........proj partons seen by targ
154 if(ltarg(it).ge.1)then
155 absb=max(1.e-9,bk(k))
157 zkt=epscrw*exp(-b2/2./b2x)*fscra(engy/egyscr)
158 c epsilongt(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
159 c & epscrs*min(epscrx,zkt*abs(betD(0,iclpro,icltar))))
160 c epsilongt(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar,
161 c & epscrh*min(epscrx,zkt*abs(betD(1,iclpro,icltar))))
162 if(ltarg3(it).ge.1)then
167 absb=max(1.e-9,abs(bk(kt))-bcut)
169 zkt=epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr)
171 zparpro(k)=zparpro(k)+zkt
172 c alpfom=max(alpfom,dble(zparpro(k)))
190 c calculation of epsilongp epsilongt
200 ip=iproj(k) !...........projectile
201 if(lproj(ip).gt.0)then
203 epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
204 & min(epscrx,epscrs*x))
206 c & epscrs*x*abs(betDp(0,iclpro,icltar)))
207 c & epscrs*min(epscrx,x*abs(betDp(0,iclpro,icltar))))
208 epsilongp(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar,
209 & epscrh*min(epscrx,x*abs(betDp(1,iclpro,icltar))))
210 cc epsilongp(k,0)=epscrs*min(epscrx,x*abs(betD(0,iclpro,icltar)))
211 cc epsilongp(k,1)=epscrh*min(epscrx,x*abs(betD(1,iclpro,icltar)))
219 it=itarg(k) !...........target
220 if(ltarg(it).gt.0)then
222 epsilongt(k,0)=max(-betDpp(0,iclpro,icltar)-0.95+alppar,
223 & min(epscrx,epscrs*x))
225 c & epscrs*x*abs(betDpp(0,iclpro,icltar)))
226 c & epscrs*min(epscrx,x*abs(betDpp(0,iclpro,icltar))))
227 epsilongt(k,1)=max(-betDpp(1,iclpro,icltar)-0.95+alppar,
228 & epscrh*min(epscrx,x*abs(betDpp(1,iclpro,icltar))))
229 cc epsilongt(k,0)=epscrs*min(epscrx,x*abs(betD(0,iclpro,icltar)))
230 cc epsilongt(k,1)=epscrh*min(epscrx,x*abs(betD(1,iclpro,icltar)))
231 cc gammaV(k)=gammaV(k)
255 !..............alpha beta betap for Gtilde (->PhiExpo).......................
258 if(iomega.eq.2)imax=1
267 if(b.lt.(nbkbin-1)*bkbin)then
268 ibk=int(bk(k)/bkbin)+1
269 if(isetcs.gt.1.and.iclegy.lt.iclegy2)then
270 egy0=egylow*egyfac**real(iclegy-1)
271 xkappa1=xkappafit(iclegy,iclpro,icltar,ibk)
272 * +(bk(k)-bkbin*real(ibk-1))/bkbin
273 * *(xkappafit(iclegy,iclpro,icltar,ibk+1)
274 * -xkappafit(iclegy,iclpro,icltar,ibk))
275 xkappa2=xkappafit(iclegy+1,iclpro,icltar,ibk)
276 * +(bk(k)-bkbin*real(ibk-1))/bkbin
277 * *(xkappafit(iclegy+1,iclpro,icltar,ibk+1)
278 * -xkappafit(iclegy+1,iclpro,icltar,ibk))
279 xkappa=xkappa1+(xkappa2-xkappa1)/log(egyfac)
280 * *(log(engy)-log(egy0))
282 xkappa=xkappafit(iclegy,iclpro,icltar,ibk)
283 * +(bk(k)-bkbin*real(ibk-1))/bkbin
284 * *(xkappafit(iclegy,iclpro,icltar,ibk+1)
285 * -xkappafit(iclegy,iclpro,icltar,ibk))
290 gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
291 gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
300 gamb=gamD(i,iclpro,icltar)*b2+epsG
301 atildg(i,k)=dble(alpD(i,iclpro,icltar))
304 c if(i.eq.0) atildg(i,k)=atildg(i,k)
305 atildg(i,k)=atildg(i,k)
306 * *dble(xkappa*xkappa)
308 atildg(i,k)=atildg(i,k)*dble(
309 * chad(iclpro)*chad(icltar)
310 * *exp(-b2/delD(i,iclpro,icltar))
311 * *sy**(betD(i,iclpro,icltar)+gamb))
312 * *gfactorp *gfactort
315 btildgp(i,k)=dble(betDp(i,iclpro,icltar)
318 btildgpp(i,k)=dble(betDpp(i,iclpro,icltar)
322 absb=abs(bk(k))-bmxdif(iclpro,icltar)
324 atildg(i,k)=atildg(i,k)*dble(
325 * sy**betD(i,iclpro,icltar)
326 * *exp(-b2a/delD(i,iclpro,icltar)))
327 btildgp(i,k)=dble(betDp(i,iclpro,icltar)-alppar)+1.d0
328 btildgpp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar)+1.d0
330 coefgdp=utgam2(1.d0+dble(alplea(iclpro))+btildgp(i,k))
331 coefgdt=utgam2(1.d0+dble(alplea(icltar))+btildgpp(i,k))
332 atildg(i,k)=atildg(i,k)
333 * *utgam2(btildgp(i,k))*utgam2(btildgpp(i,k))
335 !...........prepare plot in xepsilon
336 if(irea.eq.0.and.i.lt.2)then
337 kk=max(1,int((bk(k)-b1)/db)+1)
338 if(i.eq.0)w(0,kk)=w(0,kk)+1
339 if(i.eq.0)w(1,kk)=w(1,kk)+epsGp
340 if(i.eq.0)w(2,kk)=w(2,kk)+epsGt
341 if(i.eq.0)w(3+i,kk)=w(3+i,kk)+0
342 if(i.eq.1)w(3+i,kk)=w(3+i,kk)+0
343 !...5-8 soft ... 9-12 semi
344 w(5+4*i,kk)=w(5+4*i,kk)
345 * +betDp(i,iclpro,icltar) !prj eff
347 w(6+4*i,kk)=w(6+4*i,kk)
348 * +betDpp(i,iclpro,icltar) !tgt eff
350 w(7+4*i,kk)=w(7+4*i,kk)
351 * +betDp(i,iclpro,icltar) !prj unscr
353 w(8+4*i,kk)=w(8+4*i,kk)
354 * +betDpp(i,iclpro,icltar) !tgt unscr
356 if(i.eq.0)w(13,kk)=w(13,kk)+zparpro(k)
357 if(i.eq.0)w(14,kk)=w(14,kk)+zpartar(k)
363 !...........................................................................
368 if(ktot.le.nyeps)then
369 w(15,ktot)=w(15,ktot)+zpav/koll
370 w(16,ktot)=w(16,ktot)+ztav/koll
371 w(17,ktot)=w(17,ktot)+1
373 n=1+float(nglevt)/(0.1*maproj*matarg)*(nyeps-1)
374 if(nglevt.ge.1.and.n.ge.1.and.n.le.nyeps)then
375 w(18,n)=w(18,n)+zpav/koll
376 w(19,n)=w(19,n)+ztav/koll
381 !........alpha beta betap for Gpro...........................................
390 if(iomega.eq.2)imax=1
398 gamb=gamD(i,iclpro,icltar)*b2+epsG
400 atilde(i,k)=dble(alpD(i,iclpro,icltar))
402 atilde(i,k)=atilde(i,k)*dble(
403 * exp(-b2/delD(i,iclpro,icltar))
404 * *sy**(betD(i,iclpro,icltar)
406 * *chad(iclpro)*chad(icltar))
409 btildep(i,k)=dble(betDp(i,iclpro,icltar)
412 btildepp(i,k)=dble(betDpp(i,iclpro,icltar)
416 absb=abs(bk(k))-bmxdif(iclpro,icltar)
418 atilde(i,k)=atilde(i,k)*dble(
419 * sy**betD(i,iclpro,icltar)
420 * *exp(-b2a/delD(i,iclpro,icltar)))
422 btildep(i,k)=dble(betDp(i,iclpro,icltar)-alppar)
423 btildepp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar)
426 if(btildep(i,k)+1.d0.lt.-eps.or.btildepp(i,k)+1.d0.lt.-eps)
427 & call utstop('Error in epos-omg in GfunPark&')
433 !...........................................................................
439 write(ifch,*)' k,b,ip,it,',k,bk(k),ip,it
440 write(ifch,*)' gammaV,epsilonGP1/2,epsilonGP1/2'
441 * ,gammaV(k),epsilongp(k,0),epsilongp(k,1)
442 * ,epsilongt(k,0),epsilongt(k,1)
443 write(ifch,*)'*******************************************'
444 write(ifch,*)" atilde,btildep,btildepp "
446 write(ifch,*)i,atilde(i,k),btildep(i,k),btildepp(i,k)
451 call utprix('GfunParK',ish,ishini,5)
456 c----------------------------------------------------------------------
457 subroutine GfunPar(m,i,b,spp,alp,bet,betp,epsp,epst,epss,gamvv)
458 c----------------------------------------------------------------------
459 c calculates parameters alp,bet,betp of the G functions for pp (b-mode)
460 c and the screening exponents epsp,epst,epss,gamvv
461 c----------------------------------------------------------------------
462 c m=1: profile function Gpro, i=0: soft i=2: diff
463 c m=2: Gtilde, i=1: semi (no screening for diff)
464 c b: impact param, spp: pp energy squared
465 c----------------------------------------------------------------------
468 include 'epos.incsem'
469 include 'epos.incpar'
470 include 'epos.incems'
472 common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
473 common /kwrite/ xkapZ
474 parameter(eps=1.e-20)
475 double precision utgam2
477 call utprj('GfunPar ',ish,ishini,6)
479 !???????????????????????????????????????????????????estimation
482 if(iclpro+icltar.eq.3)then !pi+p
483 sigtotf=15.*x**0.15+55.*x**(-1.5) !fit to data
484 sigelaf=20.*(x-1)**(-3)+6.*x**(-0.4)+0.2*alog(x)**2. !fit to data
485 elseif(iclpro+icltar.eq.5)then !K+p
486 sigtotf=12.5*x**0.15+35.*x**(-1.5) !fit to data
487 sigelaf=15.*(x-1)**(-3)+3.*x**(-0.4)+0.2*alog(x)**2 !fit to data
489 sigtotf=21*x**0.17+44.*x**(-0.8) !fit to data
490 sigelaf=30.*(x-1)**(-3)+17*x**(-0.47)+0.3*alog(x)**2 !fit to data
492 siginex=sigtotf-sigelaf
496 bglaub2=max(1.,siginex/10./pi)
498 if(ish.ge.6)write(ifch,*)'Gf in:',m,i,x,b,sigine,siginex,bglaub
500 !!!!!print*,'++++++++++++++++++',x,siginex,bglaub
505 gamb=gamD(i,iclpro,icltar)*b2
509 ctp060829 bzero=bglaub
510 absb=max(1.e-9,abs(b)-bcut)
512 b2x=epscrp*epscrp*bglaub2
513 zzp=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr)
514 zzt=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr)
517 epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar
518 & ,min(epscrx,epscrs*x))
520 c & ,epscrs*x*abs(betDp(i,iclpro,icltar)))
521 c & ,epscrs*min(epscrx,x*abs(betDp(i,iclpro,icltar))))
524 epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar
525 & ,epscrh*min(epscrx,x*abs(betDp(i,iclpro,icltar))))
532 epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar
533 & ,min(epscrx,epscrs*x))
535 c & ,epscrs*x*abs(betDpp(i,iclpro,icltar)))
536 c & ,epscrs*min(epscrx,x*abs(betDpp(i,iclpro,icltar))))
539 epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar
540 & ,epscrh*min(epscrx,x*abs(betDpp(i,iclpro,icltar))))
555 gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
556 gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
558 rho=betD(i,iclpro,icltar)
564 alp=alpD(i,iclpro,icltar)
567 * *exp(rho*log(spp)-b2/delD(i,iclpro,icltar))
568 bet=betDp(i,iclpro,icltar)
571 betp=betDpp(i,iclpro,icltar)
575 absb=abs(b)-bmxdif(iclpro,icltar)
578 * *exp(betD(i,iclpro,icltar)*log(spp)
579 * -b2a/delD(i,iclpro,icltar))
580 bet=betDp(i,iclpro,icltar)-alppar
581 betp=betDpp(i,iclpro,icltar)-alppar
584 if((bet+1.0).lt.-eps.or.(betp+1.0).lt.-eps)then
585 write(*,*)'m,i,b,spp,alp,bet,betp',m,i,b,spp,alp,bet,betp
586 call utstop('Error : beta < -1 in Gfunpar in epos-omg&')
591 c if(i.eq.0.and.b.lt.(nbkbin-1)*bkbin)then
592 if(b.lt.(nbkbin-1)*bkbin)then
594 if(isetcs.gt.1.and.iclegy.lt.iclegy2)then
595 egy0=egylow*egyfac**real(iclegy-1)
596 xkappa1=xkappafit(iclegy,iclpro,icltar,ibk)
597 * +(b-bkbin*real(ibk-1))/bkbin
598 * *(xkappafit(iclegy,iclpro,icltar,ibk+1)
599 * -xkappafit(iclegy,iclpro,icltar,ibk))
600 xkappa2=xkappafit(iclegy+1,iclpro,icltar,ibk)
601 * +(b-bkbin*real(ibk-1))/bkbin
602 * *(xkappafit(iclegy+1,iclpro,icltar,ibk+1)
603 * -xkappafit(iclegy+1,iclpro,icltar,ibk))
604 xkappa=xkappa1+(xkappa2-xkappa1)/log(egyfac)
605 * *(log(ee)-log(egy0))
607 xkappa=xkappafit(iclegy,iclpro,icltar,ibk)
608 * +(b-bkbin*real(ibk-1))/bkbin
609 * *(xkappafit(iclegy,iclpro,icltar,ibk+1)
610 * -xkappafit(iclegy,iclpro,icltar,ibk))
615 alp=alpD(i,iclpro,icltar)
624 * *exp(rho*log(spp)-b2/delD(i,iclpro,icltar))
625 * *chad(iclpro)*chad(icltar)
626 * *gfactorp *gfactort
627 bet=betDp(i,iclpro,icltar)
630 betp=betDpp(i,iclpro,icltar)
634 absb=abs(b)-bmxdif(iclpro,icltar)
637 * *exp(betD(i,iclpro,icltar)*log(spp)
638 * -b2a/delD(i,iclpro,icltar))
639 bet=betDp(i,iclpro,icltar)-alppar+1.
640 betp=betDpp(i,iclpro,icltar)-alppar+1.
642 coefgdp=utgam2(1.d0+dble(alplea(iclpro))+bet)
643 coefgdt=utgam2(1.d0+dble(alplea(icltar))+betp)
644 alp=alp*utgam2(dble(bet))*utgam2(dble(betp))/coefgdp/coefgdt
647 stop'GproPar: wrong m value. '
661 if(ish.ge.6)write(ifch,*)' GfunPar :',alp,bet,betp,epsp,epst
664 call utprjx('GfunPar ',ish,ishini,6)
667 c----------------------------------------------------------------------
668 subroutine GfomPar(b,spp)
669 c----------------------------------------------------------------------
670 c calculates parameters of the fom functions for pp (b-mode)
671 c----------------------------------------------------------------------
672 c b: impact param, spp: pp energy squared
673 c----------------------------------------------------------------------
676 include 'epos.incsem'
677 include 'epos.incpar'
678 include 'epos.incems'
679 parameter(eps=1.e-20)
681 call utprj('GfomPar ',ish,ishini,6)
683 !???????????????????????????????????????????????????estimation
686 if(iclpro+icltar.eq.3)then !pi+p
687 sigtotf=15.*x**0.15+55.*x**(-1.5) !fit to data
688 sigelaf=20.*(x-1)**(-3)+6.*x**(-0.4)+0.2*alog(x)**2. !fit to data
689 elseif(iclpro+icltar.eq.5)then !K+p
690 sigtotf=12.5*x**0.15+35.*x**(-1.5) !fit to data
691 sigelaf=15.*(x-1)**(-3)+3.*x**(-0.4)+0.2*alog(x)**2 !fit to data
693 sigtotf=21*x**0.17+44.*x**(-0.8) !fit to data
694 sigelaf=30.*(x-1)**(-3)+17*x**(-0.47)+0.3*alog(x)**2 !fit to data
696 siginex=sigtotf-sigelaf
700 bglaub2=max(1.,siginex/10./pi)
707 ctp060829 bzero=bglaub
708 absb=max(1.e-9,abs(b)-bcut)
710 b2x=epscrp*epscrp*bglaub2
711 zzp=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr)
712 zzt=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr)
719 z0=alpfomi!*epscrw*fscra(ee/egyscr)
722 zpUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
723 c zpUni=dble(4.*z0*(z1/z0)**1.5)
725 ztUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
726 c ztUni=dble(4.*z0*(z1/z0)**1.5)
732 if(ish.ge.6)write(ifch,*)' GfomPar :',zpUni,ztUni
734 call utprjx('GfomPar ',ish,ishini,6)
737 c----------------------------------------------------------------------
739 c----------------------------------------------------------------------
742 if(x2.gt.1.)fscra=log(x2)
745 c----------------------------------------------------------------------
747 c----------------------------------------------------------------------
748 include 'epos.incpar'
751 if(x2.gt.1.)fscro=sqrt(log(x2)**2+epscrv**2)
754 c----------------------------------------------------------------------
755 subroutine recalcZPtn !???????????? not updated !!!
756 c----------------------------------------------------------------------
759 include 'epos.incems'
760 stop'recalcZPtn not valid any more!!!!!!!'
761 c if(koll.eq.1.and.maproj.eq.1.and.matarg.eq.1)then
766 c zparpro(k)=max(0,npom-1)*0.2
768 c zpartar(k)=max(0,npom-1)*0.2
776 c----------------------------------------------------------------------
777 double precision function om1(xh,yp,b) !---test---
778 c----------------------------------------------------------------------
779 c om1 = G * C * gamd (C and gamd usually 1)
780 c xh - fraction of the energy squared s for the pomeron;
781 c b - impact parameter between the pomeron ends;
782 c yp - rapidity for the pomeron;
783 c----------------------------------------------------------------------
785 include 'epos.incsem'
786 include 'epos.incpar'
788 double precision Gf,xp,xm,xh,yp
795 if(iomega.eq.2)imax=1
797 call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
798 Gf=Gf+alp*real(xp)**bet*real(xm)**betp
801 * * dble(chad(iclpro)*chad(icltar))
804 c----------------------------------------------------------------------
805 double precision function om1intb(b) !---MC---
806 c----------------------------------------------------------------------
807 c om1 integrated over xp and xm for given b
808 c Calculation by analytical integration
809 c----------------------------------------------------------------------
811 include 'epos.incsem'
812 include 'epos.incpar'
813 parameter(eps=1.e-20)
817 if(iomega.eq.2)imax=1
820 call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
822 if((bet+1.0).gt.eps)then
823 cint2=cint2/(bet+1.0)
825 cint2=-cint2*log(xminDf)
827 if((betp+1.0).gt.eps)then
828 cint2=cint2/(betp+1.0)
830 cint2=-cint2*log(xminDf)
836 write(*,*) 'WARNING ! om1intb in epos-omg is <0 !!!!!'
837 write(*,*) 'WARNING ! => om1intb set to 1e-3 !!!!!'
838 write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
843 * *chad(iclpro)*chad(icltar)
848 c----------------------------------------------------------------------
849 double precision function om1intbi(b,iqq) !---MC---
850 c----------------------------------------------------------------------
851 c om1 integrated over xp and xm for given b
852 c Calculation by analytical integration of contribution iqq
853 c----------------------------------------------------------------------
855 include 'epos.incsem'
856 include 'epos.incpar'
857 parameter(eps=1.e-20)
860 call GfunPar(1,iqq,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
862 if((bet+1.0).gt.eps)then
865 cint=-cint*log(xminDf)
867 if((betp+1.0).gt.eps)then
870 cint=-cint*log(xminDf)
873 write(*,*) 'WARNING ! om1intbi in epos-omg is <0 !!!!!'
874 write(*,*) 'WARNING ! => om1intbi set to 1e-3 !!!!!'
875 write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
880 * *chad(iclpro)*chad(icltar)
885 c----------------------------------------------------------------------
886 double precision function om1intbc(b) !---MC---
887 c----------------------------------------------------------------------
888 c om1*F*F integrated over xp and xm for given b
889 c Calculation by analytical integration
890 c----------------------------------------------------------------------
892 include 'epos.incems'
893 include 'epos.incsem'
894 double precision cint,gamom,deltap,deltam
895 double precision utgam2,Fp,Fm
899 Fp=dble(ucfpro) !gamma(1+alplea)
903 if(iomega.eq.2)imax=1
908 call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
909 gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
912 cint=cint+gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0)
913 & /utgam2(2.d0+deltap+dble(alplea(iclpro)))
914 & /utgam2(2.d0+deltam+dble(alplea(icltar)))
919 if(om1intbc.lt.-1.d-10)then
920 write(*,*) 'WARNING ! om1intbc in epos-omg is <0 !!!!!'
921 write(*,*) 'WARNING ! => om1intbc set to 0. !!!!!'
928 cc----------------------------------------------------------------------
929 c double precision function om1intbci(b,iqq) !---MC---
930 cc----------------------------------------------------------------------
931 cc om1*F*F integrated over xp and xm for given b and given Pomeron type iqq
932 cc Calculation by analytical integration
933 cc----------------------------------------------------------------------
935 c include 'epos.incems'
936 c include 'epos.incsem'
937 c double precision cint,gamom,deltap,deltam
938 c double precision utgam2,Fp,Fm,eps
942 c Fp=dble(ucfpro) !gamma(1+alplea)
946 c call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
947 c gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
950 c cint=gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0)
951 c & /utgam2(2.d0+deltap+dble(alplea(iclpro)))
952 c & /utgam2(2.d0+deltam+dble(alplea(icltar)))
954 c om1intbci=cint*Fp*Fm
956 c if(om1intbci.lt.-1.d-10)then
957 c write(*,*) 'WARNING ! om1intbci in epos-omg is <0 !!!!!'
958 c write(*,*) 'WARNING ! => om1intbci set to 0. !!!!!'
965 c----------------------------------------------------------------------
966 double precision function om1intgck(k,xprem,xmrem) !---MC---
967 c----------------------------------------------------------------------
968 c om1*(xprem-xp)*(xmrem-xm) integrated over xp and xm for given k
969 c Calculation by analytical integration
970 c----------------------------------------------------------------------
972 include 'epos.incems'
973 include 'epos.incsem'
974 double precision cint,gamom,deltap,deltam,xprem,xmrem
979 if(iomega.eq.2)imax=1
983 gamom=dble(atilde(i,k))
984 deltap=dble(btildep(i,k))
985 deltam=dble(btildepp(i,k))
986 cint=cint+gamom/(deltap+1.d0)/(deltam+1.d0)
987 & /(2.d0+deltap) /(2.d0+deltam)
988 & *xprem**(deltap+2.d0)
989 & *xmrem**(deltam+2.d0)
996 c----------------------------------------------------------------------
997 double precision function om1intgc(b) !---test---
998 c----------------------------------------------------------------------
999 c om1*(1-xp)*(1-xm) integrated over xp and xm for given b
1000 c Calculation by analytical integration
1001 c----------------------------------------------------------------------
1003 include 'epos.incpar'
1004 include 'epos.incsem'
1005 double precision cint,gamom,deltap,deltam,eps
1006 parameter(eps=1.d-20)
1012 if(iomega.eq.2)imax=1
1019 call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
1020 gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
1023 if((deltap+1.d0).gt.eps)then
1024 gamom=gamom/(deltap+1.0)
1026 gamom=-gamom*log(xminDf)
1028 if((deltam+1.d0).gt.eps)then
1029 gamom=gamom/(deltam+1.0)
1031 gamom=-gamom*log(xminDf)
1033 cint=cint+gamom /(2.d0+deltap) /(2.d0+deltam)
1037 if(om1intgc.lt.eps)then
1038 write(*,*) b,deltap,deltam,gamom
1039 write(*,*) 'WARNING ! om1intgc in epos-omg is <0 !!!!!'
1040 write(*,*) 'WARNING ! => om1intgc set to 0. !!!!!'
1048 c----------------------------------------------------------------------
1049 subroutine integom1(irea)
1050 c----------------------------------------------------------------------
1051 c om1 integrated over xp and xm for all b=bk(k) if irea=0
1052 c result written to :
1053 c om1int(k) = om1intb(bk(k)) = \int om1
1054 c om1intc(k) = om1intgc(bk(k)... = \int om1*(1-xp)*(1-xm)
1055 c----------------------------------------------------------------------
1057 include 'epos.incems'
1058 include 'epos.incsem'
1059 double precision om1intb,om1intgc
1061 call utpri('intom1',ish,ishini,5)
1068 om1int(k)=om1intb(bk(k))
1069 om1intc(k)=om1intgc(bk(k))
1074 call utprix('intom1',ish,ishini,5)
1081 c----------------------------------------------------------------------
1082 double precision function om1xpk(xp,xpr1i,xmr1i,k) !---test---
1083 c----------------------------------------------------------------------
1084 c \int dxm om1*(1-xp)*(1-xm) (normalized)
1086 c----------------------------------------------------------------------
1088 include 'epos.incpar'
1089 include 'epos.incems'
1090 double precision xp,gamomx(ntymi:ntymx),cint,gamom
1091 double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps
1092 & ,xpr1,xmr1,xpr1i,xmr1i
1093 parameter(eps=1.d-20)
1096 if(xp.ge.xpr1i)return
1102 if(iomega.eq.2)imax=1
1105 deltap(i)=btildep(i,k)
1106 deltam(i)=btildepp(i,k)
1107 gamomx(i)=atilde(i,k)*xmr1**(deltam(i)+2.d0)/(2.d0+deltam(i))
1108 if((deltam(i)+1.d0).gt.eps)then
1109 gamomx(i)=gamomx(i)/(deltam(i)+1.0)
1111 gamomx(i)=-gamomx(i)*log(xminDf)
1118 gamom=gamomx(i)*xpr1**(deltap(i)+2.d0)/(2.d0+deltap(i))
1119 if((deltap(i)+1.d0).gt.eps)then
1120 gamom=gamom/(deltap(i)+1.d0)
1122 gamom=-gamom*log(xminDf)
1129 om1xpk=om1xpk+gamomx(i)*xp**deltap(i)
1140 c----------------------------------------------------------------------
1141 double precision function om1xmk(xp,xm,xpr1i,xmr1i,k) !---test---
1142 c----------------------------------------------------------------------
1143 c om1(xp,xm)*(1-xp)*(1-xm) (normalized for fixed xp)
1145 c----------------------------------------------------------------------
1147 include 'epos.incpar'
1148 include 'epos.incems'
1149 double precision xp,xm,gamomx(ntymi:ntymx),cint,gamom
1150 double precision deltam(ntymi:ntymx),eps,xpr1,xmr1,xpr1i,xmr1i
1151 parameter(eps=1.d-20)
1154 if(xp.ge.xpr1i)return
1155 if(xm.ge.xmr1i)return
1161 if(iomega.eq.2)imax=1
1164 gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xpr1-xp)
1165 deltam(i)=btildepp(i,k)
1170 gamom=gamomx(i)*xmr1**(deltam(i)+2.d0)/(2.d0+deltam(i))
1171 if((deltam(i)+1.d0).gt.eps)then
1172 gamom=gamom/(deltam(i)+1.0)
1174 gamom=-gamom*log(xminDf)
1181 om1xmk=om1xmk+gamomx(i)*xm**deltam(i)*(xmr1-xm)
1189 c----------------------------------------------------------------------
1190 double precision function om1xprk(k,xpremi,xmremi,ir) !---MC---
1191 c----------------------------------------------------------------------
1192 c Random number generated from the function om1xpk. We solve the equation
1193 c which give om1xprk by Newton-Raphson + secant method.
1195 c ir - 1 to get xp, -1 to get xm (be carrefull to inverse xpremi et xmremi
1196 c when calling with ir=-1)
1197 c----------------------------------------------------------------------
1199 include 'epos.incpar'
1200 include 'epos.incems'
1201 double precision x,x0,x1,gamomx(ntymi:ntymx),eps,f0t,f1t,f00
1202 double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,drangen,xmrem
1203 double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),xprem
1204 parameter (eps=1.d-20)
1205 double precision xpremi,xmremi,lxmin
1216 if(iomega.eq.2)imax=1
1220 deltap(i)=btildep(i,k)
1221 deltam(i)=btildepp(i,k)
1223 deltap(i)=btildepp(i,k)
1224 deltam(i)=btildep(i,k)
1226 gamomx(i)=atilde(i,k)*xmrem**(deltam(i)+2.d0)/(2.d0+deltam(i))
1227 if((deltam(i)+1.d0).gt.eps)then
1228 gamomx(i)=gamomx(i)/(deltam(i)+1.0)
1231 gamomx(i)=-gamomx(i)*lxmin
1238 if((deltap(i)+1.d0).gt.eps)then
1240 & *(xprem*exp(x0)**(1.d0+deltap(i))/(1.d0+deltap(i))
1241 & -exp(x0)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1243 & *(xprem*exp(x1)**(1.d0+deltap(i))/(1.d0+deltap(i))
1244 & -exp(x1)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1247 f0=f0+gamomx(i)*(xprem*(x0-lxmin)-exp(x0)+xminDf)
1248 f1=f1+gamomx(i)*(xprem*(x1-lxmin)-exp(x1)+xminDf)
1257 r=drangen(dble(ntry))
1260 if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1261 if(f1t*f0t.ge.eps)then
1263 write(ifmt,*)i,gamomx(i),deltap(i),deltam(i)
1265 write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k
1266 call utstop('om1xprk (2)&')
1270 if(abs(f0).le.eps) then
1274 if(abs(f1).le.eps) then
1285 if(ntry.le.1000)then
1289 if((deltap(i)+1.d0).gt.eps)then
1291 & *(xprem*exp(x)**(1.d0+deltap(i))/(1.d0+deltap(i))
1292 & -exp(x)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1293 fpx=fpx+gamomx(i)*exp(x)**deltap(i)*(xprem-exp(x))
1295 fx=fx+gamomx(i)*(xprem*(x-lxmin)-exp(x)+xminDf)
1296 fpx=fpx+gamomx(i)*(xprem/exp(x)-1.d0)
1303 if (f0*fx.lt.0.D0) then
1310 if ((xt.lt.x0.or.xt.gt.x1).and.abs(f1-f0).gt.eps) then
1311 xt=x1-f1*(x1-x0)/(f1-f0)
1316 write(ifmt,*)'Warning in om1xprk, to much try !'
1321 if(abs(x-xt).gt.deltx*0.5d0) then
1325 if(abs(x).gt.eps)then
1328 call utstop('Problem in om1xprk&')
1331 if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1342 c----------------------------------------------------------------------
1343 double precision function om1xmrk(k,xp,xpremi,xmremi,ir) !---MC---
1344 c----------------------------------------------------------------------
1345 c Random number generated from the function om1xmk. We solve the equation
1346 c which give om1xmrk by Newton-Raphson + secant method.
1348 c ir - 1 to get xm, -1 to get xp (be carrefull to inverse xpremi et xmremi
1349 c when calling with ir=-1)
1350 c----------------------------------------------------------------------
1352 include 'epos.incpar'
1353 include 'epos.incems'
1354 double precision x,x0,x1,gamomx(ntymi:ntymx),eps,xp,f0t,f1t
1355 double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,f00
1356 double precision deltam(ntymi:ntymx),drangen,xprem,xmrem
1357 double precision xpremi,xmremi
1358 parameter (eps=1.d-20)
1361 if(xp.ge.xpremi)return
1369 if(iomega.eq.2)imax=1
1373 gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xprem-xp)
1374 deltam(i)=btildepp(i,k)
1376 gamomx(i)=atilde(i,k)*xp**btildepp(i,k)*(xprem-xp)
1377 deltam(i)=btildep(i,k)
1386 if(abs(deltam(i)+1.d0).gt.eps)then
1388 & *(xmrem*exp(x0)**(1.d0+deltam(i))/(1.d0+deltam(i))
1389 & -exp(x0)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1391 & *(xmrem*exp(x1)**(1.d0+deltam(i))/(1.d0+deltam(i))
1392 & -exp(x1)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1395 f0=f0+gamomx(i)*(xmrem*(x0-lxmin)-exp(x0)+xminDf)
1396 f1=f1+gamomx(i)*(xmrem*(x1-lxmin)-exp(x1)+xminDf)
1405 r=drangen(dble(ntry))
1408 if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1409 if(f1t*f0t.ge.eps)then
1410 write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry
1411 call utstop('Error(2) in epos-omg in om1xmrk&')
1415 if(abs(f0).lt.eps) then
1419 if(abs(f1).lt.eps) then
1430 if(ntry.le.1000)then
1434 if(abs(deltam(i)+1.d0).gt.eps)then
1436 & *(xmrem*exp(x)**(1.d0+deltam(i))/(1.d0+deltam(i))
1437 & -exp(x)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1438 fpx=fpx+gamomx(i)*exp(x)**deltam(i)*(xmrem-exp(x))
1440 fx=fx+gamomx(i)*(xmrem*(x-lxmin)-exp(x)+xminDf)
1441 fpx=fpx+gamomx(i)*(xmrem/exp(x)-1.d0)
1448 if (f0*fx.lt.-eps) then
1455 if ((xt.lt.x0-eps.or.xt.gt.x1+eps).and.abs(f1-f0).gt.eps) then
1456 xt=x1-f1*(x1-x0)/(f1-f0)
1461 write(ifmt,*)'Warning in om1xmrk, to much try !'
1465 if(abs(x-xt).gt.deltx*0.5d0) then
1469 if(abs(x).gt.eps)then
1472 call utstop('Problem in om1xmrk&')
1475 if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1486 c----------------------------------------------------------------------
1487 double precision function om1xk(xh,k) !---test---
1488 c----------------------------------------------------------------------
1489 c \int dxp om1 (normalised)
1491 c----------------------------------------------------------------------
1493 include 'epos.incpar'
1494 include 'epos.incems'
1496 double precision xh,gamomx(ntymi:ntymx),cint,alpp(ntymi:ntymx)
1497 &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps
1499 parameter(eps=1.d-20)
1506 if(iomega.eq.2)imax=1
1509 gamomx(i)=atilde(i,k)
1510 deltap(i)=btildep(i,k)
1511 deltam(i)=btildepp(i,k)
1513 delta(i)=(deltap(i)+deltam(i))*0.5d0
1514 alpp(i)=deltap(i)-deltam(i)
1521 if((deltap(i)+1.d0).gt.eps)then
1522 gamom=gamom/(deltap(i)+1.d0)
1524 gamom=-gamom*log(xminDf)
1526 if((deltam(i)+1.d0).gt.eps)then
1527 gamom=gamom/(deltam(i)+1.d0)
1529 gamom=-gamom*log(xminDf)
1536 if(abs(alpp(i)).gt.eps)then
1537 om1xk=om1xk+gamomx(i)/alpp(i)*xh**deltam(i)*(1.d0-xh**alpp(i))
1539 om1xk=om1xk-gamomx(i)*xh**delta(i)*log(xh)
1548 c----------------------------------------------------------------------
1549 double precision function om1yk(xh,yp,k) !---test---
1550 c----------------------------------------------------------------------
1551 c om1 normalized for fixed xp
1552 c xh - fraction of the energy squared s for the pomeron;
1554 c----------------------------------------------------------------------
1556 include 'epos.incems'
1557 double precision xh,yp,gamomy(ntymi:ntymx),alpp(ntymi:ntymx),cint
1558 double precision deltap,deltam,eps
1559 parameter(eps=1.d-20)
1565 if(iomega.eq.2)imax=1
1568 gamomy(i)=atilde(i,k)
1570 deltam=btildepp(i,k)
1572 alpp(i)=deltap-deltam
1573 gamomy(i)=gamomy(i)*xh**((deltap+deltam)*0.5d0)
1579 if(abs(alpp(i)).gt.eps)then
1580 cint=cint-gamomy(i)/alpp(i)*xh**(alpp(i)*0.5d0)
1581 & *(1.d0-xh**(-alpp(i)))
1583 cint=cint-gamomy(i)*log(xh)
1588 if(abs(alpp(i)).gt.eps)then
1589 om1yk=om1yk+gamomy(i)*exp(alpp(i)*yp)
1591 om1yk=om1yk+gamomy(i)
1601 c----------------------------------------------------------------------
1602 double precision function om1xrk(k) !---test---
1603 c----------------------------------------------------------------------
1604 c Random number generated from the function om1xk. We solve the equation
1605 c which give om1xrk by Newton-Raphson + secant method.
1607 c----------------------------------------------------------------------
1609 include 'epos.incems'
1610 include 'epos.incpar'
1612 double precision x,x0,x1,gamomx(ntymi:ntymx),eps,prec,drangen
1613 double precision xt,fx,fpx,r,f1,f0,cint,deltx,alpp(ntymi:ntymx)
1614 &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx)
1615 parameter (eps=1.d-20)
1622 if(iomega.eq.2)imax=1
1625 gamomx(i)=atilde(i,k)
1626 deltap(i)=btildep(i,k)
1627 deltam(i)=btildepp(i,k)
1629 delta(i)=(deltap(i)+deltam(i))*0.5d0
1630 alpp(i)=deltap(i)-deltam(i)
1637 if((deltap(i)+1.d0).gt.eps)then
1638 gamom=gamom/(deltap(i)+1.d0)
1640 gamom=-gamom*log(xminDf)
1642 if((deltam(i)+1.d0).gt.eps)then
1643 gamom=gamom/(deltam(i)+1.d0)
1645 gamom=-gamom*log(xminDf)
1656 if(abs(alpp(i)).lt.eps)then
1657 if(delta(i)+1.d0.gt.eps)then
1658 f0=f0-gamomx(i)/(delta(i)+1.d0)*x0**(delta(i)+1.d0)
1659 & *(log(x0)-1.d0/(delta(i)+1.d0))
1660 f1=f1-gamomx(i)/(delta(i)+1.d0)*x1**(delta(i)+1.d0)
1661 & *(log(x1)-1.d0/(delta(i)+1.d0))
1663 f0=f0-0.5d0*gamomx(i)*(log(x0)**2-log(xminDf)**2)
1664 f1=f1-0.5d0*gamomx(i)*(log(x1)**2-log(xminDf)**2)
1667 if(abs(deltap(i)+1.d0).gt.eps
1668 & .and.abs(deltam(i)+1.d0).gt.eps)then
1669 f0=f0+gamomx(i)/alpp(i)*(x0**(deltam(i)+1.d0)/(deltam(i)+1.d0)
1670 & -x0**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1671 f1=f1+gamomx(i)/alpp(i)*(x1**(deltam(i)+1.d0)/(deltam(i)+1.d0)
1672 & -x1**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1673 elseif(abs(deltap(i)+1.d0).gt.eps)then
1674 f0=f0+gamomx(i)/alpp(i)*(log(x0/xminDf)
1675 & -x0**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1676 f1=f1+gamomx(i)/alpp(i)*(log(x1/xminDf)
1677 & -x1**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1678 elseif(abs(deltam(i)+1.d0).gt.eps)then
1679 f0=f0-gamomx(i)/alpp(i)*(log(x0/xminDf)
1680 & -x0**(deltam(i)+1.d0)/(deltam(i)+1.d0))
1681 f1=f1-gamomx(i)/alpp(i)*(log(x1/xminDf)
1682 & -x1**(deltam(i)+1.d0)/(deltam(i)+1.d0))
1690 r=drangen(dble(ntry))
1693 if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1694 if(f1t*f0t.ge.eps)then
1696 write(ifmt,*)i,gamomx(i),deltap(i),deltam(i),alpp(i),delta(i)
1698 write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k
1699 call utstop('om1xrk (1)&')
1703 c if(f1*f0.gt.eps)then
1704 c call utmsg('om1xrk')
1705 c write(ifch,*)'Poblem with x0, no root ! --> om1xrk=xminDf'
1706 c write(ifmt,*)'Poblem with x0, no root ! --> om1xrk=xminDf'
1707 c write(ifmt,*)f0,f1,cint,r
1712 if(abs(f0).lt.eps) then
1716 if(abs(f1).lt.eps) then
1730 if(ntry.le.1000)then
1734 if(abs(alpp(i)).lt.eps)then
1735 if(delta(i)+1.d0.gt.eps)then
1736 fx=fx-gamomx(i)/(delta(i)+1.d0)*x**(delta(i)+1.d0)
1737 & *(log(x)-1.d0/(delta(i)+1.d0))
1738 fpx=fpx-gamomx(i)*x**delta(i)*log(x)
1740 fx=fx-0.5d0*gamomx(i)*(log(x)**2-log(xminDf)**2)
1741 fpx=fpx-gamomx(i)*log(x)/x
1744 if(abs(deltap(i)+1.d0).gt.eps
1745 & .and.abs(deltam(i)+1.d0).gt.eps)then
1746 fx=fx+gamomx(i)/alpp(i)*(x**(deltam(i)+1.d0)/(deltam(i)+1.d0)
1747 & -x**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1748 fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i))
1749 elseif(abs(deltap(i)+1.d0).gt.eps)then
1750 fx=fx+gamomx(i)/alpp(i)*(log(x/xminDf)
1751 & -x**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1752 fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i))
1753 elseif(abs(deltam(i)+1.d0).gt.eps)then
1754 fx=fx-gamomx(i)/alpp(i)*(log(x/xminDf)
1755 & -x**(deltam(i)+1.d0)/(deltam(i)+1.d0))
1756 fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i))
1764 if (f0*fx.lt.-eps) then
1771 if ((xt.lt.x0-eps.or.xt.gt.x1+eps).and.abs(f1-f0).gt.eps) then
1772 xt=x1-f1*(x1-x0)/(f1-f0)
1777 write(ifmt,*)'Warning in om1xrk, to much try !'
1781 if(abs(x-xt).gt.deltx*0.5d0) then
1785 if(abs(x).gt.eps)then
1788 call utstop('Problem in om1xrk&')
1791 if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000)then
1802 c----------------------------------------------------------------------
1803 double precision function om1yrk(xh,k) !---test---
1804 c----------------------------------------------------------------------
1805 c Random number generated from the function om1yk(xh). We solve the
1806 c equation which give om1yrk by Newton-Raphson + secant method.
1807 c xh - fraction of the energy squared s for the pomeron;
1809 c----------------------------------------------------------------------
1811 include 'epos.incems'
1813 double precision xh,r!,y0,y1,y,gamomy(ntymi:ntymx),eps,ymin,prec,yt
1817 om1yrk=(0.5d0-r)*log(xh)
1822 c----------------------------------------------------------------------
1823 function ffom12aii(iq,je1,je2) !---test---
1824 c----------------------------------------------------------------------
1832 xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5))
1836 xp=xmin+(xmax-xmin)*(.5+tgss(ig,i1)*(m1-1.5))
1837 f=ffom12a(xp,xm,iq,iq,je1,je2)
1841 f=r1*0.5*(xmax-xmin)
1845 ffom12aii=r2*0.5*(xmax-xmin)
1848 c----------------------------------------------------------------------
1849 function ffom12ai(xp,iq1,iq2,je1,je2) !---test---
1850 c----------------------------------------------------------------------
1858 xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5))
1859 f=ffom12a(xp,xm,iq1,iq2,je1,je2)
1863 ffom12ai=r2*0.5*(xmax-xmin)
1866 c----------------------------------------------------------------------
1867 function ffom12a(xp,xm,iq1,iq2,je1,je2) !---test---
1868 c----------------------------------------------------------------------
1870 c 2*om52*F*F == PomInc
1875 c iq1 - min iq iq=2 .... val-sea
1876 c iq2 - max iq iq=3 .... sea-val
1878 c je = emission type (projectile and target side)
1879 c 0 ... no emissions
1883 c already b-averaged (\int d2b /sigine*10)
1884 c----------------------------------------------------------------------
1889 ctp060829 yp=0.5*log(xp/xm)
1893 ffom12a=ffom12a+2*om52pi(sy*xh,1.,1.,0,je1,je2)
1895 ffom12a=ffom12a+2*om52pi(sy*xh,xp,1.,1,je1,je2)
1897 ffom12a=ffom12a+2*om52pi(sy*xh,xm,1.,2,je1,je2)
1899 ffom12a=ffom12a+2*om52pi(sy*xh,xp,xm,3,je1,je2)
1903 * *alpff(iclpro)*xp**betff(1)*(1-xp)**alplea(iclpro)
1904 * *alpff(icltar)*xm**betff(2)*(1-xm)**alplea(icltar)
1908 c----------------------------------------------------------------------
1909 function ffom11a(xp,xm,iq1,iq2) !---test---
1910 c----------------------------------------------------------------------
1912 c int(db) om1ff /sigine*10
1914 c xp - xplus iq=-1 ... fit
1915 c xm - xminus iq=0 .... soft
1917 c iq1 - min iq iq=2 .... qg
1918 c iq2 - max iq iq=3 .... gq
1921 c----------------------------------------------------------------------
1923 common/geom/rmproj,rmtarg,bmax,bkmx
1929 bb=bmid*(1.+(2.*m-3)*tgss(ig,i))
1930 f=ffom11(xp,xm,bb,iq1,iq2)
1934 ffom11a=r*2.*pi*bmid /sigine*10
1938 c----------------------------------------------------------------------
1939 function ffom11(xp,xm,b,iq1,iq2) !---test---
1940 c----------------------------------------------------------------------
1942 c 2*om5*F*F == PomInc
1944 c xp - xplus iq=-1 ... fit
1945 c xm - xminus iq=0 .... soft
1946 c b - impact parameter iq=1 .... gg
1947 c iq1 - min iq iq=2 .... qg
1948 c iq2 - max iq iq=3 .... gq
1951 c----------------------------------------------------------------------
1953 double precision om51
1959 ffom11=2*om51(dble(xh),dble(yp),b,iq1,iq2)
1960 * *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
1962 else !xm integration
1970 xmm=xmin*(xmax/xmin)**(.5+tgss(ig,i)*(m-1.5))
1973 f=2*om51(dble(xh),dble(yp),b,iq1,iq2)
1974 * *(1-xmm)**alplea(icltar)*(1-xp)**alplea(iclpro)
1975 r=r+wgss(ig,i)*f*xmm
1978 ffom11=r*0.5*log(xmax/xmin)
1984 c----------------------------------------------------------------------
1985 double precision function om51(xh,yp,b,iq1,iq2) !---test---
1986 c----------------------------------------------------------------------
1987 c xh - xplus*xminus iq=-1 ... fit (om1 * 0.5)
1988 c yp - rapidity iq=0 .... soft
1989 c b - impact param iq=1 .... gg
1990 c iq1 - min iq iq=2 .... qg
1991 c iq2 - max iq iq=3 .... gq
1994 c----------------------------------------------------------------------
1996 include 'epos.incsem'
1997 include 'epos.incpar'
1998 double precision xp,xm,xh,yp,om51p,om1
2000 if(xh.le.0.d0.or.xh.gt.1.d0)return
2006 if(iq1.eq.-1.and.iq2.eq.-1)then
2007 om51=0.5d0*om1(xh,yp,b)
2008 elseif(iq1.ge.0)then
2013 call GfunPar(1,i1,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2014 om51=om51+om51p(sy*real(xh),xh,yp,b,i)
2015 * *xp**dble(epsp)*xm**dble(epst)
2016 * *dble(sy)**dble(epss)
2018 call GfunPar(1,2,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2019 om51=om51+0.5d0*alp*real(xp)**bet*real(xm)**betp
2023 stop'om5: choice of iq1 and iq2 is nonsense. '
2028 c----------------------------------------------------------------------
2029 double precision function om5s(sx,xh,yp,b,iq1,iq2) !---test---
2030 c----------------------------------------------------------------------
2032 double precision om51
2033 double precision xh,yp
2037 om5s=om51(xh,yp,b,iq1,iq2)
2041 c----------------------------------------------------------------------
2042 double precision function om5Jk(k,xh,yp,iqq) !---MC---
2043 c----------------------------------------------------------------------
2045 c xh - fraction of the energy squared s for the pomeron;
2046 c b - impact parameter between the pomeron ends;
2047 c yp - rapidity for the pomeron;
2053 c iqq=5 - diffractif
2054 c----------------------------------------------------------------------
2056 include 'epos.incems'
2057 include 'epos.incsem'
2059 double precision xh,yp,om51p
2060 double precision plc,s
2067 om5Jk=om51p(sy,xh,yp,b,iqq)
2068 if(iscreen.ne.0)then
2073 epsGp=epsilongp(k,0)
2074 epsGt=epsilongt(k,0)
2077 epsGp=epsilongp(k,1)
2078 epsGt=epsilongt(k,1)
2080 epsG=exp(-gfactor*(zparpro(k)*zpartar(k)))*epsG
2081 epsGp=exp(-gfactor*(zparpro(k)*zpartar(k)))*epsGp
2082 epsGt=exp(-gfactor*(zparpro(k)*zpartar(k)))*epsGt
2083 om5Jk=om5Jk*xp**dble(epsGp)*xm**dble(epsGt)*s**dble(epsG)
2088 om5Jk=0.5d0*atilde(2,k)*xp**btildep(2,k)*xm**btildepp(2,k)
2094 c----------------------------------------------------------------------
2095 double precision function omIgamint(b,iqq) !---test---
2096 c----------------------------------------------------------------------
2097 c - integrated chi~(b)FF/2 for cut I diagram (simple Pomeron)
2098 c b - impact parameter between the pomeron ends;
2099 c yp - rapidity for the pomeron;
2100 c iqq=0 effective one
2103 c----------------------------------------------------------------------
2105 include 'epos.incems'
2106 include 'epos.incsem'
2107 include 'epos.incpar'
2115 if(iomega.eq.2)imax=1
2118 coefp=1.+alplea(iclpro)
2119 coeft=1.+alplea(icltar)
2122 call GfunPar(1,i,b,sy,alpx,betx,betpx,epsp,epst,epss,gamv)
2125 Df=alpx*dble(utgam1(betp)*utgam1(betpp)*ucfpro
2126 * *ucftar/utgam1(betp+coefp)/utgam1(betpp+coeft))
2127 omIgamint=omIgamint+Df
2130 call utstop('Wrong iqq in omIgamint&')
2134 * *dble(chad(iclpro)*chad(icltar))
2136 omIgamint=omIgamint*0.5d0
2141 c-----------------------------------------------------------------------
2142 subroutine WomTy(w,xh,yp,k)
2143 c-----------------------------------------------------------------------
2144 c - w(ity) for group iqq of cut enhanced diagram giving
2145 c the probability of the type of the same final state.
2147 c xh - fraction of the energy squared s for the pomeron;
2148 c yp - rapidity for the pomeron;
2149 c xpr,xmr impulsion fraction of remnant
2155 c-----------------------------------------------------------------------
2157 include 'epos.incems'
2158 include 'epos.incsem'
2159 doubleprecision xh,yp,om5Jk,w(0:7)
2166 w(i)=om5Jk(k,xh,yp,i)
2169 if(iotst1.gt.0)then !????????????????????????????????????????????????
2170 corfac=float(iotst1)
2171 ww=w(0)+w(1)+w(2)+w(3)+w(4)+w(5)
2185 w(0)=w(0)/w05*(ww-whard)
2186 w(5)=w(5)/w05*(ww-whard)
2191 !write(*,'(2f11.4)')(ww-w05)/ww,(ww-w(0)-w(5))/ww
2192 endif !?????????????????????????????????????????????????????????????????
2198 c-----------------------------------------------------------------------
2199 double precision function Womegak(xp,xm,xprem,xmrem,k) !---MC---
2200 c-----------------------------------------------------------------------
2201 c - sum(omGam(xp,xm))*(1-xp)*(1-xm) for group of cut enhanced
2202 c diagram giving the same final state (without nuclear effect).
2203 c xp,xm - fraction of the loght cone momenta of the pomeron;
2205 c-----------------------------------------------------------------------
2207 include 'epos.incems'
2208 double precision xp,xm,xprem,xmrem
2213 if(iomega.eq.2)imax=1
2216 Womegak=Womegak+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k)
2220 Womegak=Womegak*(xprem-xp)*(xmrem-xm)
2226 cc----------------------------------------------------------------------
2227 c double precision function omNpcut(xp,xm,xprem,xmrem,bh,iqq) !---test---
2228 cc----------------------------------------------------------------------
2229 cc Sum of all cut diagrams
2231 cc iqq=1 exact G + diff
2232 cc----------------------------------------------------------------------
2233 c include "epos.inc"
2234 c double precision om51,xh,yp,xprem,xmrem,xp,xm!,omYcutI
2238 c if(abs(xh).gt.1.d-10)then
2239 c yp=0.5d0*log(xp/xm)
2244 c if(iqq.eq.0)omNpcut=om51(xh,yp,bh,-1,-1)
2245 c if(iqq.eq.1)omNpcut=om51(xh,yp,bh,0,5)
2247 c omNpcut=omNpcut*2.d0
2252 c----------------------------------------------------------------------
2253 double precision function omGam(xp,xm,bh) !---test---
2254 c-----------------------------------------------------------------------
2255 c Cut diagram part for calculation of probability distribution
2256 c xp,xm impulsion fraction of remnant
2257 c bh - impact parameter between the pomeron ends;
2258 c-----------------------------------------------------------------------
2260 include "epos.incems"
2261 double precision om51,xp,xm,xh,yp,eps!,omYgam
2262 parameter (eps=1.d-20)
2265 if(xp.lt.eps.or.xm.lt.eps)return
2267 if(abs(xh).gt.1.d-10)then
2273 omGam=om51(xh,yp,bh,-1,-1)
2280 c----------------------------------------------------------------------
2281 double precision function omGamk(k,xp,xm) !---MC---
2282 c-----------------------------------------------------------------------
2283 c Cut diagram part for calculation of probability distribution (for omega)
2284 c xp,xm impulsion fraction of remnant
2285 c bh - impact parameter between the pomeron ends;
2286 c-----------------------------------------------------------------------
2288 include "epos.incems"
2289 double precision xp,xm
2292 if(iomega.eq.2)imax=1
2294 omGamk=omGamk+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k)
2300 c----------------------------------------------------------------------
2301 double precision function omGamint(bh) !---test---
2302 c-----------------------------------------------------------------------
2303 c Integrated cut diagram part for calculation of probability distribution
2304 c bh - impact parameter between the pomeron ends;
2305 c-----------------------------------------------------------------------
2307 double precision omIgamint!,omYgamint
2309 omGamint=2.d0*omIgamint(bh,0)
2318 c----------------------------------------------------------------------
2320 c----------------------------------------------------------------------
2321 c constants for numerical integration (gaussian weights)
2322 c----------------------------------------------------------------------
2323 double precision dgx1,dga1
2324 common /dga20/ dgx1(10),dga1(10)
2328 & .765265211334973D-01,
2329 & .227785851141645D+00,
2330 & .373706088715420D+00,
2331 & .510867001950827D+00,
2332 & .636053680726515D+00,
2333 & .746331906460151D+00,
2334 & .839116971822219D+00,
2335 & .912234428251326D+00,
2336 & .963971927277914D+00,
2337 & .993128599185095D+00/
2339 & .152753387130726D+00,
2340 & .149172986472604D+00,
2341 & .142096109318382D+00,
2342 & .131688638449177D+00,
2343 & .118194531961518D+00,
2344 & .101930119817233D+00,
2345 & .832767415767047D-01,
2346 & .626720483341090D-01,
2347 & .406014298003871D-01,
2348 & .176140071391506D-01/
2354 c----------------------------------------------------------------------
2355 double precision function PhiExact(fj,xp,xm,s,b) !---test---
2356 c----------------------------------------------------------------------
2357 c Exact expression of the Phi function for pp collision
2358 c----------------------------------------------------------------------
2360 include 'epos.incsem'
2361 include 'epos.incems'
2362 double precision al(idxD0:idxD1),betp(idxD0:idxD1)
2364 double precision zp(idxD0:idxD1),Phitmp,betpp(idxD0:idxD1)
2366 double precision eps
2367 parameter(eps=1.d-20)
2368 dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
2370 if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in PhiExact"
2373 if(xp.gt.eps.and.xm.gt.eps.and.xp.le.1.d0+eps
2374 & .and.xm.le.1.d0+eps)then
2387 if(iomega.eq.2)imax0=1
2390 imax(i)=10+max(5,int(log10(s)))
2391 if(b.ge.1.)imax(i)=4+max(3,int(log10(sqrt(s))))
2392 imax(i)=min(30,imax(i))
2396 call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
2397 betp(i)=dble(betx)+1.d0
2398 betpp(i)=dble(betpx)+1.d0
2399 al(i)=dble(alpx*gamv)
2400 * *dble(chad(iclpro)*chad(icltar))
2406 if (ipr(0).ne.0) zp(0)=(-dble(fj)*al(0))**ipr(0)*facto(ipr(0))
2410 if (ipr(1).ne.0) zp(1)=(-dble(fj)*al(1))**ipr(1)*facto(ipr(1))
2414 if (ipr(2).ne.0) zp(2)=(-dble(fj)*al(2))**ipr(2)*facto(ipr(2))
2420 yp=yp+dble(ipr(i))*betp(i)
2421 ym=ym+dble(ipr(i))*betpp(i)
2426 z=z*xIrst(1,xp,yp,betp,ipr)
2427 z=z*xIrst(2,xm,ym,betpp,ipr)
2446 c----------------------------------------------------------------------
2447 double precision function PhiExpoK(k,xp,xm) !---MC---
2448 c----------------------------------------------------------------------
2449 c Exponential expression of the Phi function for pp collision
2451 c----------------------------------------------------------------------
2453 include 'epos.incsem'
2454 include 'epos.incems'
2456 double precision xp,xm,Phitmp,Gt1
2457 double precision atildg,btildgp,btildgpp
2458 common/cgtilde/atildg(idxD0:idxD1,kollmx)
2459 *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx)
2465 if(iomega.eq.2)imax=1
2470 Gt1=Gt1+atildg(i,k)*xp**btildgp(i,k)*xm**btildgpp(i,k)
2480 c----------------------------------------------------------------------
2481 double precision function PhiExpo(fj,xp,xm,s,b) !---MC---
2482 c----------------------------------------------------------------------
2483 c Exponential expression of the Phi function for pp collision
2485 c----------------------------------------------------------------------
2487 include 'epos.incsem'
2488 include 'epos.incems'
2489 include 'epos.incpar'
2491 parameter(nbkbin=40)
2492 common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
2493 double precision AlTi
2494 double precision BeTip,BeTipp
2495 double precision xp,xm,Phitmp,Gt1
2498 if(iomega.eq.2)imax=1
2502 call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
2506 Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp*dble(fj*xkappa**(2.5*fj-2.5))
2512 & *xp**dble(alplea(iclpro))
2513 & *xm**dble(alplea(icltar))
2518 c----------------------------------------------------------------------
2519 double precision function PhiUnit(xp,xm) !---test---
2520 c----------------------------------------------------------------------
2521 c Exponential expression of the Phi function for pp collision
2523 c----------------------------------------------------------------------
2525 include 'epos.incsem'
2526 include 'epos.incems'
2527 include 'epos.incpar'
2529 double precision AlTi
2530 double precision BeTip,BeTipp
2531 double precision xp,xm,Phitmp,Gt1
2534 if(iomega.eq.2)imax=1
2541 Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp
2542 c write(ifch,*)'Phiunit',i,xp,xm,Gt1,AlTi,BeTip,BeTipp
2548 & *xp**dble(alplea(iclpro))
2549 & *xm**dble(alplea(icltar))
2556 cc----------------------------------------------------------------------
2557 c double precision function PhiUnit(xp,xm,s,b) !---inu---
2558 cc----------------------------------------------------------------------
2559 c include 'epos.inc'
2560 c double precision xp,xm,PhiExpo,Znorm
2562 c PhiUnit=PhiExpo(1.,xp,xm,s,b)
2569 c----------------------------------------------------------------------
2570 double precision function Hrst(s,b,xp,xm) !test
2571 c----------------------------------------------------------------------
2573 include 'epos.incems'
2574 include 'epos.incsem'
2575 include 'epos.incpar'
2577 double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
2578 common/DGamUni/GbetUni( idxD0:idxD2),HbetUni( idxD0:idxD2),
2579 & GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
2580 & HalpUni(idxD0:idxD2)
2581 double precision al(idxD0:idxD2),betp(idxD0:idxD2)
2583 double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2)
2585 dimension ipr(idxD0:idxD2),imax(idxD0:idxD2)
2587 if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in Hrst"
2596 if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then
2599 if(iomega.eq.2)imax0=1
2601 if(iomega.eq.2)imax1=imax1-1
2604 imax(i)=max(2,int(log10(100.*s)/3.))
2605 c if(i.ge.2)imax(i)=imax(i)*2
2606 if(b.ge.1.5)imax(i)=2 !max(2,imax(i)/2)
2607 imax(i)=min(30,imax(i))
2609 if((zpUni*ztUni.lt.1.d-6)
2610 & .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
2613 imax(i)=1 !imax(i)/3
2621 betpp(i)=HbetpUni(i)
2626 c write(ifmt,*)'Hrst ipr0,xp,xm :',ipr0,xp,xm
2629 if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2633 if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2637 if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
2641 if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3))
2645 if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4))
2649 if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5))
2653 if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6))
2657 if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7))
2661 if (ipr(8).ne.0) zp(8)=al(8)**ipr(8)*facto(ipr(8))
2662 if (ipr(0)+ipr(1)+ipr(2)+ipr(3)+ipr(4)+ipr(5)
2663 & +ipr(6)+ipr(7)+ipr(8).ne.0) then
2668 yp=yp+dble(ipr(i))*betp(i)
2669 ym=ym+dble(ipr(i))*betpp(i)
2672 z=z*xJrst(xp,yp,GbetUni,ipr)
2673 z=z*xJrst(xm,ym,GbetpUni,ipr)
2693 c----------------------------------------------------------------------
2694 double precision function HrstI(s,b,xp,xm) !test
2695 c----------------------------------------------------------------------
2697 include 'epos.incems'
2698 include 'epos.incsem'
2699 include 'epos.incpar'
2701 double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
2702 common/DGamUni/GbetUni( idxD0:idxD2),HbetUni( idxD0:idxD2),
2703 & GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
2704 & HalpUni(idxD0:idxD2)
2705 double precision al(idxD0:idxD2),betp(idxD0:idxD2)
2707 double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2)
2709 dimension ipr(idxD0:idxD2),imax(idxD0:idxD2)
2711 if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
2721 if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then
2725 if(iomega.eq.2)imax0=1
2727 if(iomega.eq.2)imax1=imax1-1
2730 imax(i)=max(3,int(log10(s)/2.))
2731 c if(i.ge.2)imax(i)=imax(i)*2
2732 if(b.ge.1.5)imax(i)=max(2,imax(i)/2)
2733 imax(i)=min(30,imax(i))
2735 if((zpUni*ztUni.lt.1.d-6)
2736 & .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
2739 imax(i)=1 !imax(i)/3
2747 betpp(i)=HbetpUni(i)
2753 if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2757 if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2761 if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
2765 if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3))
2769 if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4))
2773 if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5))
2777 if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6))
2781 if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7))
2785 if (ipr(8).ne.0) zp(8)=al(8)**ipr(8)*facto(ipr(8))
2786 if (ipr(0)+ipr(1)+ipr(2)+ipr(3)+ipr(4)+ipr(5)
2787 & +ipr(6)+ipr(7)+ipr(8).ne.0) then
2792 yp=yp+dble(ipr(i))*betp(i)
2793 ym=ym+dble(ipr(i))*betpp(i)
2796 z=z*xJrstI(xp,yp,GbetUni,ipr)
2797 z=z*xJrstI(xm,ym,GbetpUni,ipr)
2819 cc----------------------------------------------------------------------
2820 c double precision function HrstI(s,xp,xm) !---inu---
2821 cc----------------------------------------------------------------------
2822 c include 'epos.inc'
2823 c include 'epos.incems'
2824 c include 'epos.incsem'
2825 c include 'epos.incpar'
2826 c double precision al(idxD0:idxD1),betp(idxD0:idxD1)
2827 c *,z,xJrstI!,ffacto
2828 c double precision zp(idxD0:idxD1),Htmp,betpp(idxD0:idxD1)
2830 c dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
2832 c if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
2841 c if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.lt.1.d0)then
2846 c if(iomega.eq.2)imax0=1
2848 c do i=idxDmin,imax0
2849 c imax(i)=max(5,int(log10(s)))
2850 cc if(i.ge.2)imax(i)=imax(i)*2
2851 c imax(i)=min(30,imax(i))
2854 c do i=idxDmin,imax0
2855 c betp(i)=betUni(i,1)+1.d0
2856 c betpp(i)=betpUni(i,1)+1.d0
2857 c al(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
2863 c if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2867 c if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2871 c if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
2872 c if (ipr(0)+ipr(1)+ipr(2).ne.0) then
2876 c do i=idxDmin,imax0
2877 c yp=yp+dble(ipr(i))*betp(i)
2878 c ym=ym+dble(ipr(i))*betpp(i)
2881 c z=z*xJrstI(xp,yp,betp,ipr)
2882 c z=z*xJrstI(xm,ym,betpp,ipr)
2897 cc----------------------------------------------------------------------
2898 c double precision function ffacto(n) !---test---
2899 cc----------------------------------------------------------------------
2903 c ffacto=ffacto*dble(i)
2909 c----------------------------------------------------------------------
2910 double precision function xIrst(id,x,y,bet,ipr) !---test---
2911 c----------------------------------------------------------------------
2913 double precision y,gammag,utgam2,x,bet(idxD0:idxD1)
2914 dimension ipr(idxD0:idxD1)
2916 if(id.eq.1)iclrem=iclpro
2917 if(id.eq.2)iclrem=icltar
2919 if(iomega.eq.2)imax=1
2921 xIrst=gammag(iclrem,y)*x**dble(alplea(iclrem))
2925 if(xIrst.gt.0.d0)then
2927 if(ipr(i).ne.0.and.bet(i).gt.1.d-10)
2928 & xIrst=xIrst*utgam2(bet(i))**dble(ipr(i))
2930 if (abs(y).gt.1.d-10) xIrst=xIrst*x**y
2936 c----------------------------------------------------------------------
2937 double precision function xJrst(x,y,Gbeta,ipr) !---inu---
2938 c----------------------------------------------------------------------
2941 double precision y,utgam2,x,Gbeta(idxD0:idxD2),eps,gam
2942 dimension ipr(idxD0:idxD2)
2947 if(iomega.eq.2)imax=imax-1
2951 if(gam.lt.1.d99)then
2953 if ((x-1.d0).gt.eps.or.(y-1.d0).gt.eps) then
2954 xJrst=(1.d0-x)**(y-1.d0)/gam
2956 if (ipr(i).ne.0) xJrst=xJrst*Gbeta(i)**dble(ipr(i))
2959 c write (*,*) 'Warning in xJrst, infinite value !'
2960 xJrst=(1.d0-x+eps)**(y-1.d0)/gam
2962 if (ipr(i).ne.0) xJrst=xJrst*Gbeta(i)**dble(ipr(i))
2973 c----------------------------------------------------------------------
2974 double precision function xJrstI(x,y,Gbeta,ipr) !---inu---
2975 c----------------------------------------------------------------------
2976 c Function used for the integration of H*Phi. We do the changement of
2977 c variable (1-x)=z**alpha. The power alpha can be change if necessary.
2978 c----------------------------------------------------------------------
2981 double precision y,utgam2,x,Gbeta(idxD0:idxD2),alpha,w,gam
2982 dimension ipr(idxD0:idxD2)
2985 w=alpha*(y-1.d0)+alpha-1.d0
2987 if(iomega.eq.2)imax=imax-1
2991 if(gam.lt.1.d99)then
2995 xJrstI=alpha*x**w/gam
2997 if (ipr(i).ne.0) xJrstI=xJrstI*Gbeta(i)**dble(ipr(i))
3001 write(*,*) 'x,y,bet,ipr,w',x,y,Gbeta,ipr,w
3002 stop 'Error in xJrstI in epos-omg, integration not possible'
3013 c----------------------------------------------------------------------
3014 double precision function HPhiInt(s,b) !---inu---
3015 c----------------------------------------------------------------------
3016 c Set integrated over xp and xm (x and y) H(x,y)*Phi(x,y) for a
3017 c given b by gauss method
3018 c----------------------------------------------------------------------
3021 double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
3022 common/DGamUni/GbetUni( idxD0:idxD2),HbetUni( idxD0:idxD2),
3023 & GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
3024 & HalpUni(idxD0:idxD2)
3025 double precision xhm,x,y,yhm,w,Hrst,utgam2,PhiUnit!,PhiExact
3026 c double precision zp2,zm2,HrstI,eps
3027 c common /ar3/ x1(7),a1(7)
3028 common /ar9/ x9(3),a9(3)
3039 HbetUni(i)=betUni(i,1)+1.d0
3040 HbetpUni(i)=betpUni(i,1)+1.d0
3041 GbetUni(i)=utgam2(HbetUni(i))
3042 GbetpUni(i)=utgam2(HbetpUni(i))
3043 HalpUni(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
3046 HbetUni(imax0+1+i)=betUni(i,1)+1.d0+betfom
3047 HbetUni(imax0+3+i)=betUni(i,1)+1.d0
3048 HbetUni(imax0+5+i)=betUni(i,1)+1.d0+betfom
3049 HbetpUni(imax0+1+i)=betpUni(i,1)+1.d0
3050 HbetpUni(imax0+3+i)=betpUni(i,1)+1.d0+betfom
3051 HbetpUni(imax0+5+i)=betpUni(i,1)+1.d0+betfom
3052 GbetUni(imax0+1+i)=utgam2(HbetUni(imax0+1+i))
3053 GbetUni(imax0+3+i)=utgam2(HbetUni(imax0+3+i))
3054 GbetUni(imax0+5+i)=utgam2(HbetUni(imax0+5+i))
3055 GbetpUni(imax0+1+i)=utgam2(HbetpUni(imax0+1+i))
3056 GbetpUni(imax0+3+i)=utgam2(HbetpUni(imax0+3+i))
3057 GbetpUni(imax0+5+i)=utgam2(HbetpUni(imax0+5+i))
3058 HalpUni(imax0+1+i)=ztUni*alpUni(i,1)
3059 HalpUni(imax0+3+i)=zpUni*alpUni(i,1)
3060 HalpUni(imax0+5+i)=zpUni*ztUni*alpUni(i,1)
3069 x=xhm+dble((2*m-3)*x9(i))*xhm
3070 c write(ifmt,*)'HPhiInt, xp int :',x
3074 y=yhm+dble((2*n-3)*x9(j))*yhm
3075 w=w+dble(a9(i)*a9(j))*Hrst(s,b,x,y)
3077 c & *PhiExact(1.,x,y,s,b)
3091 c x=1d0-eps+xhm+dble((2*m-3)*x1(i))*xhm
3094 c y=1d0-epsyhm+dble((2*n-3)*x1(j))*yhm
3097 c w=w+dble(a1(i)*a1(j))*HrstI(s,x,y)
3098 cc & *PhiUnit(zp2,zm2)
3099 c & *PhiExact(1.,zp2,zm2,s,b)
3105 c HPhiInt=HPhiInt+w*xhm*yhm
3112 c----------------------------------------------------------------------
3113 subroutine Kfit(iiclegy)
3114 c----------------------------------------------------------------------
3116 include "epos.incsem"
3117 double precision Znorm
3118 parameter(nbkbin=40)
3119 common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
3123 if(iiclegy.eq.-1.or.iiclegy.gt.iclegy2)then
3128 xkappafit(iiiegy,iiipro,iiitar,iiibk)=1.
3140 s=(egylow*egyfac**(iiclegy-1))**2.
3144 write(ifmt,*)"Fit xkappa ..."
3146 write(ifmt,*)"Kfit s,bkbin,iclegy,ipro,itar"
3147 * ,s,bkbin,iiclegy,iclpro,icltar
3152 if(isetcs.le.1.or.iiclegy.eq.iclegy2)then
3155 xkf=xkappafit(iiclegy+1,iclpro,icltar,1)
3162 xkappafit(iiclegy,iclpro,icltar,ib)=1.
3163 if(b.gt.3.+0.05*log(s).or.s.le.20.*q2min)then
3167 if(ib.gt.1.and.ish.ge.3)write(ifch,*)" End",delta,xkf
3168 delta=1.-real(Znorm(s,b))
3169 if(delta.le.0d0)then
3171 xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3172 delta=1.-real(Znorm(s,b))
3174 else!if(xkf.ne.1.)then
3177 if(abs(delta).lt.eps)then
3178 if(delta.lt.0d0)then
3184 elseif(ib.le.nbkbin-1)then
3186 if(delta.gt.0.d0)then
3191 xkappafit(iiclegy,iclpro,icltar,ib)=xkf2
3192 delta1=1.-real(Znorm(s,b))
3193 if(delta1.lt.0.d0)then
3199 xkf1=max(delta0,xkf2)
3212 xkfs=max(0.00001,1.-delta)
3215 if(delta.le.deltas)xkf=xkfs
3216 if(ish.ge.3)write(ifch,*)" Start",ib,b,delta,xkf,xkf0,xkf1
3217 if(xkf.eq.xkf2)delta=delta1
3224 if(n.le.nmax.and.xkf1.ne.xkf0)then
3225 if(abs(xkf-xkf2).gt.1e-6.or.abs(delta).gt.abs(deltas))then
3226 xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3227 delta=1.-real(Znorm(s,b))
3229 if(ish.ge.5)write(ifch,*)" step",ib,n,delta,xkf,delta0
3230 if(delta*delta0.ge.0.)then
3231 if(lnoch.and.abs(delta).gt.abs(delta0))goto 5
3235 if(abs(delta).gt.eps)then
3249 * write(ifmt,*)"Warning in Kfit, nmax reached : xkappafit=1."
3250 xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3256 if(ish.ge.3)write(ifch,*)" End",delta,xkf
3257 100 if(xkf.gt.1.+eps)write(ifmt,*)
3258 * "Warning in Kfit, xkappafit not yet 1"
3259 xkappafit(iiclegy,iclpro,icltar,nbkbin)=1.
3266 c----------------------------------------------------------------------
3267 double precision function Znorm(s,b) !---inu---
3268 c----------------------------------------------------------------------
3270 common /kwrite/ xkapZ
3271 double precision HPhiInt,PhiUnit!,PhiExact
3273 c write(ifmt,*)'Z calculation for (s,b) :',s,b
3275 if(iomega.eq.2)imax=1
3277 call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3278 call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3282 c write(ifch,*)'int',Znorm,' phi',PhiExact(1.,1.d0,1.d0,s,b)
3284 & +PhiUnit(1.d0,1.d0)
3285 c & +PhiExact(1.,1.d0,1.d0,s,b)
3287 !write(ifmt,*)'Z=',Znorm,xkapZ,b
3292 c------------------------------------------------------------
3293 double precision function gammag(iclrem,x) !---test---
3294 c----------------------------------------------------------------------
3296 include 'epos.incsem'
3297 double precision x,utgam2
3299 gammag=utgam2(dble(alplea(iclrem))+1.D0)
3300 & /utgam2(dble(alplea(iclrem))+1.D0+x)
3307 cc----------------------------------------------------------------------
3308 c double precision function PomNbri(iqq) !---xsigma---
3309 cc----------------------------------------------------------------------
3310 cc integral d2b om1intbci
3311 cc iqq, Pomeron type
3312 cc----------------------------------------------------------------------
3313 c include 'epos.inc'
3314 c double precision om1intbci
3315 c common/geom/rmproj,rmtarg,bmax,bkmx
3316 c common /ar3/ x1(7),a1(7)
3322 c bb=bmid*(1.+(2.*m-3)*x1(i))
3323 c PomNbri=PomNbri+dble(bb*a1(i))*om1intbci(bb,iqq)
3327 c PomNbri=PomNbri*dble(2.*pi*bmid)
3336 c####################################################################################
3337 c############# former chk #########################################################
3338 c####################################################################################
3347 cc----------------------------------------------------------------------
3348 c double precision function PomIncII(b) !---check---
3349 cc----------------------------------------------------------------------
3350 cc integral_dx_dy om1*F_remn*F_remn for a given b !---check---
3351 cc----------------------------------------------------------------------
3352 c include 'epos.inc'
3353 c include 'epos.incems'
3354 c include 'epos.incsem'
3355 c include 'epos.incpar'
3356 c double precision cint,gamom(idxD0:idxD1),deltap(idxD0:idxD1)
3357 c &,deltapp(idxD0:idxD1),utgam2
3359 cc Calculation by analytical integration (perfect but it changes
3364 c if(iomega.eq.2)imax=1
3366 c call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3367 c gamom(i)=dble(alp*gamv)*chad(iclpro)*chad(icltar)
3368 c deltap(i)=dble(bet)
3369 c deltapp(i)=dble(betp)
3371 cc Integration possible only if delta(i)>-1
3373 c if(deltap(i).le.-1.d0.or.deltapp(i).le.-1.d0)
3374 c & stop 'Error in epos-par-300 in PomIncII'
3379 c cint=cint+gamom(i)*utgam2(deltap(i)+1.d0)*utgam2(deltapp(i)+1.d0)
3380 c & *dble(ucfpro*ucftar)
3381 c & /utgam2(dble(alplea(iclpro))+deltap(i)+2.d0)
3382 c & /utgam2(dble(alplea(icltar))+deltapp(i)+2.d0)
3391 c----------------------------------------------------------------------
3392 double precision function PomIncXIExact(x) !---check---
3393 c----------------------------------------------------------------------
3394 c integral d2b PomIncXExact
3395 c----------------------------------------------------------------------
3397 double precision x,PomIncXExact
3398 common /ar3/ x1(7),a1(7)
3399 common/geom/rmproj,rmtarg,bmax,bkmx
3405 bb=bmid*(1.+(2.*m-3)*x1(i))
3406 PomIncXIExact=PomIncXIExact+dble(bb*a1(i))*PomIncXExact(x,bb)
3410 PomIncXIExact=PomIncXIExact*dble(2.*pi*bmid)
3415 c----------------------------------------------------------------------
3416 double precision function PomIncXIUnit(x) !---check---
3417 c----------------------------------------------------------------------
3418 c integral d2b PomIncXUnit
3419 c----------------------------------------------------------------------
3421 double precision x,PomIncXUnit
3422 common /ar3/ x1(7),a1(7)
3423 common/geom/rmproj,rmtarg,bmax,bkmx
3429 bb=bmid*(1.+(2.*m-3)*x1(i))
3430 PomIncXIUnit=PomIncXIUnit+dble(bb*a1(i))*PomIncXUnit(x,bb)
3434 PomIncXIUnit=PomIncXIUnit*dble(2.*pi*bmid)
3439 c----------------------------------------------------------------------
3440 double precision function PomIncPIExact(x) !---check---
3441 c----------------------------------------------------------------------
3442 c integral d2b PomIncPExact
3443 c----------------------------------------------------------------------
3445 double precision x,PomIncPExact
3446 common/geom/rmproj,rmtarg,bmax,bkmx
3447 common /ar3/ x1(7),a1(7)
3453 bb=bmid*(1.+(2.*m-3)*x1(i))
3454 PomIncPIExact=PomIncPIExact+dble(bb*a1(i))*PomIncPExact(x,bb)
3458 PomIncPIExact=PomIncPIExact*dble(2.*pi*bmid)
3463 c----------------------------------------------------------------------
3464 double precision function PomIncPIUnit(x) !---check---
3465 c----------------------------------------------------------------------
3466 c integral d2b PomIncPUnit
3467 c----------------------------------------------------------------------
3469 double precision x,PomIncPUnit
3470 common/geom/rmproj,rmtarg,bmax,bkmx
3471 common /ar3/ x1(7),a1(7)
3477 bb=bmid*(1.+(2.*m-3)*x1(i))
3478 PomIncPIUnit=PomIncPIUnit+dble(bb*a1(i))*PomIncPUnit(x,bb)
3482 PomIncPIUnit=PomIncPIUnit*dble(2.*pi*bmid)
3487 c----------------------------------------------------------------------
3488 double precision function PomIncMIExact(x) !---check---
3489 c----------------------------------------------------------------------
3490 c integral d2b PomIncMExact
3491 c----------------------------------------------------------------------
3493 double precision x,PomIncMExact
3494 common/geom/rmproj,rmtarg,bmax,bkmx
3495 common /ar3/ x1(7),a1(7)
3501 bb=bmid*(1.+(2.*m-3)*x1(i))
3502 PomIncMIExact=PomIncMIExact+dble(bb*a1(i))*PomIncMExact(x,bb)
3506 PomIncMIExact=PomIncMIExact*dble(2.*pi*bmid)
3511 c----------------------------------------------------------------------
3512 double precision function PomIncMIUnit(x) !---check---
3513 c----------------------------------------------------------------------
3514 c integral d2b PomIncMUnit
3515 c----------------------------------------------------------------------
3517 double precision x,PomIncMUnit
3518 common/geom/rmproj,rmtarg,bmax,bkmx
3519 common /ar3/ x1(7),a1(7)
3525 bb=bmid*(1.+(2.*m-3)*x1(i))
3526 PomIncMIUnit=PomIncMIUnit+dble(bb*a1(i))*PomIncMUnit(x,bb)
3530 PomIncMIUnit=PomIncMIUnit*dble(2.*pi*bmid)
3535 c----------------------------------------------------------------------
3536 double precision function PomIncMExact(xm,b) !---check---
3537 c----------------------------------------------------------------------
3538 c incluse Pomeron distribution \int dx+ { 2G F_remn F_remn }
3539 c----------------------------------------------------------------------
3541 include 'epos.incsem'
3542 include 'epos.incems'
3543 double precision AlTiP,BeTip,al,bep,bepp,xpInt,utgam2,xm
3548 if(iomega.eq.2)imax=1
3550 call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3556 xpInt=utgam2(BeTip)*dble(ucfpro)
3557 * /utgam2(1.d0+dble(alplea(iclpro))+BeTip)
3559 PomIncMExact=PomIncMExact+AlTiP*xm**bepp
3560 * *(1.d0-xm)**dble(alplea(icltar))
3566 c----------------------------------------------------------------------
3567 double precision function PomIncMUnit(xm,b) !---check---
3568 c----------------------------------------------------------------------
3569 c incluse Unitarized Pomeron distribution \int dx+
3570 c----------------------------------------------------------------------
3572 include 'epos.incsem'
3573 include 'epos.incems'
3574 double precision xh,Df,xp,xm,G2,w,xpm
3575 double precision PoInU!,Znorm
3576 common /ar3/ x1(7),a1(7)
3580 c Calculation by numeric integration :
3585 xp=xpm*(1.d0+dble((2.*m-3.)*x1(j)))
3587 ctp060829 sy=s*sngl(xh)
3590 call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3591 Df=Df+dble(alp)*xp**dble(bet)*xm**dble(betp)
3594 G2=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
3595 w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2
3601 PomIncMUnit=w!/Znorm(s,b)
3607 c----------------------------------------------------------------------
3608 double precision function PomIncPExact(xp,b) !---check---
3609 c----------------------------------------------------------------------
3610 c incluse Pomeron distribution \int dx- { 2G F_remn F_remn }
3611 c----------------------------------------------------------------------
3613 include 'epos.incsem'
3614 include 'epos.incems'
3615 double precision AlTiP,BeTipp,al,bep,bepp,xmInt,utgam2,xp
3620 if(iomega.eq.2)imax=1
3622 call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3627 xmInt=utgam2(BeTipp)*dble(ucftar)
3628 * /utgam2(1.d0+dble(alplea(icltar))+BeTipp)
3630 PomIncPExact=PomIncPExact+AlTiP*xp**bep
3631 * *(1.d0-xp)**dble(alplea(iclpro))
3637 c----------------------------------------------------------------------
3638 double precision function PomIncPUnit(xp,b) !---check---
3639 c----------------------------------------------------------------------
3640 c incluse Unitarized Pomeron distribution \int dx-
3641 c----------------------------------------------------------------------
3643 include 'epos.incsem'
3644 double precision xh,Df,xp,xm,G2,w,xmm
3645 double precision PoInU!,Znorm
3646 common /ar3/ x1(7),a1(7)
3650 if(iomega.eq.2)imax=1
3652 c Calculation by numeric integration :
3657 xm=xmm*(1.d0+dble((2.*m-3.)*x1(j)))
3659 ctp060829 sy=s*real(xh)
3662 call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3663 Df=Df+alp*real(xp)**bet*real(xm)**betp
3666 G2=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
3667 w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2
3673 PomIncPUnit=w!/Znorm(s,b)
3680 c----------------------------------------------------------------------
3681 double precision function PomIncJExact(b) !---check---
3682 c----------------------------------------------------------------------
3683 c integral of Pomeron distribution \int dy dx { 2G F_remn F_remn }
3684 c----------------------------------------------------------------------
3686 double precision allea,PomIncXExact,xh
3687 common /ar3/ x1(7),a1(7)
3689 allea=2.d0+dble(alplea(iclpro)+alplea(icltar))
3693 xh=1.d0-(.5d0+dble(x1(i)*(real(m)-1.5)))**(1.d0/allea)
3694 PomIncJExact=PomIncJExact+dble(a1(i))
3695 & *PomIncXExact(xh,b)/(1.d0-xh)**(allea-1.d0)
3698 PomIncJExact=PomIncJExact/allea/2.d0
3704 c----------------------------------------------------------------------
3705 double precision function PomIncJUnit(b) !---check---
3706 c----------------------------------------------------------------------
3707 c integral of Pomeron distribution \int dy dx { 2G F_remn F_remn }
3708 c----------------------------------------------------------------------
3710 double precision PomIncXUnit,xh,xhm
3711 common /ar3/ x1(7),a1(7)
3717 xh=xhm*(1.d0+dble(x1(i)*(2.*real(m)-3.)))
3718 PomIncJUnit=PomIncJUnit+dble(a1(i))
3719 & *PomIncXUnit(xh,b)
3722 PomIncJUnit=PomIncJUnit*xhm
3728 cc----------------------------------------------------------------------
3729 c double precision function PomIncXExact1(xh,b) !---check---
3730 cc----------------------------------------------------------------------
3731 cc incluse Pomeron distribution \int dy { 2G F_remn F_remn }
3732 cc----------------------------------------------------------------------
3733 c include 'epos.inc'
3734 c include 'epos.incsem'
3735 c double precision xh,Df,xp,xm,w,ymax
3736 c common /ar3/ x1(7),a1(7)
3739 c if(iomega.eq.2)imax=1
3741 ccctp060829 sy=s*real(xh)
3742 cc Calculation by numeric integration :
3744 c ymax=-.5d0*log(xh)
3747 c xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
3751 c call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3752 c Df=Df+alp*gamv*real(xp)**bet*real(xm)**betp
3753 c * *(1.d0-xp)**dble(alplea(iclpro))
3754 c * *(1.d0-xm)**dble(alplea(icltar))
3756 c w=w+dble(a1(j))*Df
3759 c w=w*ymax*xh**dble(-alppar)
3767 c----------------------------------------------------------------------
3768 double precision function PomIncXExact(xh,b) !---check---
3769 c----------------------------------------------------------------------
3770 c incluse Pomeron distribution \int dy { 2G F_remn F_remn }
3771 c (optimized integration but with no y dependance)
3772 c----------------------------------------------------------------------
3774 include 'epos.incsem'
3775 double precision AlTiP,bep,bepp,factor,factor1
3776 double precision xpmin,xh,xp,xm,ymax,y
3777 common /ar3/ x1(7),a1(7)
3780 if(iomega.eq.2)imax=1
3784 call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3788 PomIncXExact=PomIncXExact+AlTiP*xh**((bep+bepp)/2.d0)
3792 allea=min(alplea(iclpro),alplea(icltar))+1.
3793 xpmin=max(sqrt(xh),exp(-1.d0))
3796 xp=1.d0-(1.d0-xpmin)*(.5d0+dble(x1(i)*(real(m)-1.5)))
3797 * **(1.d0/dble(allea))
3799 factor=factor+dble(a1(i))
3800 * *((1.d0-xp)**dble(alplea(iclpro)-allea+1.)
3801 * *(1.d0-xm)**dble(alplea(icltar))
3802 * +(1.d0-xp)**dble(alplea(icltar)-allea+1.)
3803 * *(1.d0-xm)**dble(alplea(iclpro)))/xp
3806 factor=factor*(1.d0-xpmin)**dble(allea)/dble(allea)
3809 if(xpmin.gt.1.00001d0*sqrt(xh))then
3814 y=ymax*dble(x1(i)*(2*m-3))
3817 factor1=factor1+dble(a1(i))*(1.d0-xp)**dble(alplea(iclpro))
3818 * *(1.d0-xm)**dble(alplea(icltar))
3821 factor=factor+factor1*ymax
3826 PomIncXExact=PomIncXExact*factor
3832 c----------------------------------------------------------------------
3833 double precision function PomIncXUnit(xh,b) !---check---
3834 c----------------------------------------------------------------------
3835 c incluse Unitarized Pomeron distribution \int dy
3836 c----------------------------------------------------------------------
3838 include 'epos.incsem'
3839 double precision xh,Df,xp,xm,w
3840 double precision PoInU,ymax!,Znorm
3841 common /ar3/ x1(7),a1(7)
3844 if(iomega.eq.2)imax=1
3846 ctp060829 sy=s*real(xh)
3847 c Calculation by numeric integration :
3852 xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
3856 call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3857 Df=Df+alp*real(xp)**bet*real(xm)**betp
3860 Df=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
3861 w=w+dble(a1(j))*Df*PoInU(xp,xm,s,b)
3867 PomIncXUnit=w!/Znorm(s,b)
3874 c----------------------------------------------------------------------
3875 double precision function PoInU(xp,xm,s,b) !---check---
3876 c----------------------------------------------------------------------
3877 c Function : PhiU(1-xp,1-xm) + /int(H(z+ + x+,z- + x-)PhiU(z+,z-)dz+dz-)
3878 c----------------------------------------------------------------------
3880 double precision xp,xm,zp,zm,zp2,zm2,zpmin,zmmin,deltzp,deltzm
3881 double precision zpm,zmm,w,HrstI,PhiUnit,Hrst,eps!,PhiExact
3882 common /ar3/ x1(7),a1(7)
3887 if(iomega.eq.2)imax=1
3889 call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3890 call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3894 if (1.d0-xp-eps.gt.0.d0.and.1.d0-xm-eps.gt.0.d0) then
3902 zp=zpm+dble((2*m-3)*x1(i))*zpm
3905 zm=zmm+dble((2*n-3)*x1(j))*zmm
3906 w=w+dble(a1(i)*a1(j))*Hrst(s,b,zp+xp,zm+xm)
3908 c & *PhiExact(1.,zp,zm,s,b)
3925 if(abs(deltzp).gt.1.d-10.and.abs(deltzm).gt.1.d-10)then
3930 zp=zpmin+zpm+dble((2*m-3)*x1(i))*zpm
3933 zm=zmmin+zmm+dble((2*n-3)*x1(j))*zmm
3936 w=w+dble(a1(i)*a1(j))*HrstI(s,b,zp,zm)
3938 c & *PhiExact(1.,zp2,zm2,s,b)
3945 PoInU=PoInU+w*zpm*zmm
3946 & +PhiUnit(1.d0-xp,1.d0-xm)
3947 c & +PhiExact(1.,1.d0-xp,1.d0-xm,s,b)
3952 c----------------------------------------------------------------------
3953 double precision function PomIncExact(xp,xm,b) !---check---
3954 c----------------------------------------------------------------------
3955 c inclusive Pomeron distribution { 2G F_remn F_remn }
3956 c----------------------------------------------------------------------
3958 include "epos.incsem"
3959 double precision Df,xp,xm!,xh
3963 ctp060829 sy=engy**2*real(xh)
3966 if(iomega.eq.2)imax=1
3968 call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3969 Df=Df+alp*gamv*real(xp)**bet*real(xm)**betp
3971 Df=dble(chad(iclpro)*chad(icltar))
3975 * *(1.d0-xp)**dble(alplea(iclpro))
3976 * *(1.d0-xm)**dble(alplea(icltar))
3981 c----------------------------------------------------------------------
3982 double precision function PomIncUnit(xp,xm,b) !---check---
3983 c----------------------------------------------------------------------
3984 c inclusive Pomeron distribution { Sum{int{G*Phi} }
3985 c----------------------------------------------------------------------
3987 include "epos.incpar"
3988 include "epos.incsem"
3989 double precision PoInU,xp,xm,om1,xh,yp
3994 if(xm.ne.0.d0)yp=0.5d0*log(xp/xm)
3995 PomIncUnit=om1(xh,yp,b)*PoInU(xp,xm,engy*engy,b)
3996 & *(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom)
4002 cc----------------------------------------------------------------------
4003 c double precision function PomIncUnitMC(xp,xm,b) !---check---
4004 cc----------------------------------------------------------------------
4005 cc inclusive Pomeron distribution { Sum{int{G*Phi} }
4006 cc----------------------------------------------------------------------
4007 c include "epos.inc"
4008 c include "epos.incpar"
4009 c include "epos.incsem"
4010 c include "epos.incems"
4011 c parameter(mmax=20)
4012 c double precision Gtab,Phitab,xxpu(mmax),xxmu(mmax)
4013 c double precision Zm,xp,xm,pGtab,Z,omNpcut,xprem,xmrem,
4017 c if(xp.lt.1.d0.and.xm.lt.1.d0)then
4032 c Gtab=omNpcut(xp,xm,sxp,sxm,b,0)
4033 c Phitab=PhiExpo(1.,sxp,sxm,sy,b)
4039 c write(*,*)"GPhi",mtmp-1,Zm,Z
4045 c if(mod(n,1000000).eq.0)write(*,*)
4046 c & "Calculation of PomIncUnit(",mtmp,")->",n
4053 c rnau=rangen()*real(xprem-sxp)
4054 c xxpu(i)=dble(rnau)
4056 c rnau=rangen()*real(xmrem-sxm)
4057 c xxmu(i)=dble(rnau)
4060 c if(sxp.lt.xprem.and.sxm.lt.xmrem)then
4062 c Gtab=omNpcut(xxpu(i),xxmu(i),xprem-sxp,xmrem-sxm,b,0)
4065 c Zm=Zm+pGtab*PhiExpo(1.,xprem-sxp,xmrem-sxm,sy,b)
4067 c if(n.lt.nmax)goto 10
4068 c Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(mtmp)
4072 c PomIncUnitMC=Z/dble(chad(iclpro)*chad(icltar))
4079 cc----------------------------------------------------------------------
4080 c double precision function PhiMCExact(xp,xm,b) !---check---
4081 cc----------------------------------------------------------------------
4082 cc virtual emissions { Sum{int{-GFF} }
4083 cc----------------------------------------------------------------------
4084 c include "epos.inc"
4085 c include "epos.incpar"
4086 c include "epos.incsem"
4087 c include "epos.incems"
4088 c parameter(mmax=20)
4089 c double precision Gtab,xxpu(mmax),xxmu(mmax)
4090 c double precision Zm,xp,xm,pGtab,Z,om51,sxp,sxm,xh,yp
4094 c if(xp.le.1.d0.and.xm.le.1.d0)then
4105 c Z=xp**dble(alplea(iclpro))
4106 c * *xm**dble(alplea(icltar))
4111 c write(*,*)"GPhi",mtmp-1,Zm,Z/xp**dble(alplea(iclpro))
4112 c * /xm**dble(alplea(icltar))
4118 c if(mod(n,1000000).eq.0)write(*,*)
4119 c & "Calculation of PhiExact(",mtmp,")->",n
4124 c rnau=rangen()!*real(xp-sxp)
4125 c xxpu(i)=dble(rnau)
4127 c rnau=rangen()!*real(xm-sxm)
4128 c xxmu(i)=dble(rnau)
4131 c if(sxp.lt.xp.and.sxm.lt.xm)then
4133 c xh=xxpu(i)*xxmu(i)
4134 c if(abs(xh).gt.1.d-30)then
4135 c yp=0.5d0*log(xxpu(i)/xxmu(i))
4139 c Gtab=2*om51(xh,yp,b,0,4)
4140 cc * +omNpuncut(sy*real(xh),xh,yp,b,1) !om1(xh,yp,b)
4141 c pGtab=pGtab*(-Gtab)
4143 c Zm=Zm+pGtab*(xp-sxp)**dble(alplea(iclpro))
4144 c * *(xm-sxm)**dble(alplea(icltar))
4146 c if(n.lt.nmax)goto 10
4147 c Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(m)
4157 c----------------------------------------------------------------------
4158 double precision function Gammapp(sy,b,mtmp) !---check---
4159 c----------------------------------------------------------------------
4161 include "epos.incpar"
4162 include "epos.incsem"
4163 include "epos.incems"
4165 double precision Gtab,xxpu(mmax),xxmu(mmax),PhiExpo
4166 double precision Zm,xp,xm,pGtab,om1,sxp,sxm,xh,yp
4184 if(mod(n,10000).eq.0)write(*,*)
4185 & "Calculation of Gammapp(",mtmp,")->",n
4190 rnau=rangen()!*real(xp-sxp)
4193 rnau=rangen()!*real(xm-sxm)
4197 if(sxp.lt.xp.and.sxm.lt.xm)then
4200 if(abs(xh).gt.1.d-30)then
4201 yp=0.5d0*log(xxpu(i)/xxmu(i))
4208 Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4210 if(n.lt.nmax)goto 10
4211 Zm=Zm/dble(nmax)!**2.d0*(xp*xm)**dble(mtmp)
4218 cc----------------------------------------------------------------------
4219 c double precision function GammaMCnew(sy,b,mtmp) !---check---
4220 cc----------------------------------------------------------------------
4221 c include "epos.inc"
4222 c include "epos.incpar"
4223 c include "epos.incsem"
4224 c include "epos.incems"
4225 c parameter(mmax=20)
4226 c common /psar7/ delx,alam3p,gam3p
4227 c double precision Gtab,xxpu(mmax),xxmu(mmax)
4228 c double precision Zm,xp,xm,pGtab,omGam,Zmtot,
4229 c * sxp,sxm,PhiExpo!,om1,yp,xh
4250 c if(mod(n,1000000).eq.0)write(*,*)
4251 c & "Calculation of GammaMCnew(",mtmp,")->",n
4257 c xxpu(i)=dble(rnau)
4260 c xxmu(i)=dble(rnau)
4263 c if(sxp.lt.xp.and.sxm.lt.xm)then
4267 c Gtab=omGam(xxpu(i),xxmu(i),b) !om1(xh,yp,b)
4270 c Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4272 c if(n.lt.nmax)goto 10
4274 c Zmtot=Zmtot+Zm/dble(nmax)
4280 cc----------------------------------------------------------------------
4281 c double precision function GammaMC(sy,b,mtmp) !---check---
4282 cc----------------------------------------------------------------------
4283 c include "epos.inc"
4284 c include "epos.incpar"
4285 c include "epos.incsem"
4286 c include "epos.incems"
4287 c parameter(mmax=20)
4288 c double precision Gtab,xxpu(mmax),xxmu(mmax)
4289 c double precision Zm,xp,xm,pGtab,om1,
4290 c * sxp,sxm,xh,yp,PhiExpo!,omNpcut
4309 c if(mod(n,1000000).eq.0)write(*,*)
4310 c & "Calculation of GammaMC(",mtmp,")->",n
4315 c rnau=rangen()!*real(xp-sxp)
4316 c xxpu(i)=dble(rnau)
4318 c rnau=rangen()!*real(xm-sxm)
4319 c xxmu(i)=dble(rnau)
4322 c if(sxp.lt.xp.and.sxm.lt.xm)then
4324 c xh=xxpu(i)*xxmu(i)
4325 c if(abs(xh).gt.1.d-30)then
4326 c yp=0.5d0*log(xxpu(i)/xxmu(i))
4330 c Gtab=om1(xh,yp,b)!omNpcut(xxpu(i),xxmu(i),xp-sxp,xm-sxm,b,0) !om1(xh,yp,b)
4333 c Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4335 c if(n.lt.nmax)goto 10
4336 c Zm=Zm/dble(nmax)*fctrl(n-mtmp)*facto(mtmp)
4346 c----------------------------------------------------------------------
4347 double precision function GammaGauss(sy,b,mtmp) !---check---
4348 c----------------------------------------------------------------------
4350 include "epos.incpar"
4351 include "epos.incsem"
4352 include "epos.incems"
4354 common /psar7/ delx,alam3p,gam3p
4355 double precision xpmin,xmmin,Gst,zm(mmax),zp(mmax)
4356 *,xpmax,xmmax,zpmin(mmax),zmmin(mmax),zpmax(mmax)
4357 double precision xp,xm,pGtab,omGam,dzp(mmax),Gp1,Gm1,xmin,eps
4358 *,sxp,sxm,PhiExpo,zmmax(mmax),dzm(mmax),Gp2,Gm2,Gp3,Gm3,G0
4360 common /ar3/ x1(7),a1(7)
4361 c double precision dgx1,dga1
4362 c common /dga20/ dgx1(10),dga1(10)
4377 elseif(mtmp.eq.1)then
4384 elseif(mtmp.eq.2)then
4391 elseif(mtmp.eq.3)then
4399 write(*,*)"m not between 0 and 3, return ..."
4423 G0=PhiExpo(1.,sxp,sxm,sy,b)
4427 c write(*,*)'x+/-',xmmin,xmmax,xpmin,xpmax
4431 if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.1)then
4432 zmmax(1)=-log(xmmin)
4433 zmmin(1)=-log(xmmax)
4434 if(abs(xmmin-xmin).lt.eps)then
4435 zmmin(1)=-log(min(xmmax,1.d0-xmmin-xmmin))
4436 zmmax(1)=-log(max(xmmin,1.d0-xmmax-xmmax))
4438 dzm(1)=(zmmax(1)-zmmin(1))/2.d0
4442 if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.1)then
4443 zpmax(1)=-log(xpmin)
4444 zpmin(1)=-log(xpmax)
4445 if(abs(xpmin-xmin).lt.eps)then
4446 zpmin(1)=-log(min(xpmax,1.d0-xpmin-xpmin))
4447 zpmax(1)=-log(max(xpmin,1.d0-xpmax-xpmax))
4449 dzp(1)=(zpmax(1)-zpmin(1))/2.d0
4452 c write(*,*)'bornes1=',exp(-zpmax(1)),exp(-zpmin(1))
4453 c &,exp(-zmmax(1)),exp(-zmmin(1))
4458 zp(1)=zpmin(1)+dzp(1)*(1.d0+dble(2.*np1-3.)*dble(x1(jp1)))
4460 if(dzm(1).eq.0.d0)then
4466 if(dzm(1).ne.0.d0)then
4467 zm(1)=zmmin(1)+dzm(1)*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1)))
4483 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4484 Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4487 & write(*,*)'j1=',k,exp(-zp(k)),exp(-zm(k))
4488 & ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1
4491 write(*,*)'error1 ?',dzp(k),dzm(k)
4494 if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4495 if(dzm(1).ne.0.d0)then
4497 &dble(a1(jm1))*PhiExpo(1.,sxp,sxm,sy,b)
4499 c &dble(a1(jm1))*PhiExact(1.,sxp,sxm,sy,b)
4503 &dble(a1(jp1))*PhiExpo(1.,sxp,sxm,sy,b)
4505 c &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b)
4508 c write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm
4513 if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.2)then
4514 zpmin(2)=-log(min(min(xpmax,1.d0-exp(-zp(1))),
4515 & 1.d0-xpmin-exp(-zp(1))))
4516 zpmax(2)=-log(max(xpmin,1.d0-xpmax-exp(-zp(1))))
4517 if(abs(xpmax+xpmax+xpmax-3.d0*dble(1./delx)).lt.eps)then
4518 zpmin(2)=-log(xpmax)
4519 zpmax(2)=-log(xpmin)
4521 dzp(2)=(zpmax(2)-zpmin(2))/2.d0
4525 if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.2)then
4526 zmmin(2)=-log(min(min(xmmax,1.d0-exp(-zm(1))),
4527 & 1.d0-xmmin-exp(-zm(1))))
4528 zmmax(2)=-log(max(xmmin,1.d0-xmmax-exp(-zm(1))))
4529 if(abs(xmmax+xmmax+xmmax-3.d0*dble(1./delx)).lt.eps)then
4530 zmmin(2)=-log(xmmax)
4531 zmmax(2)=-log(xmmin)
4533 dzm(2)=(zmmax(2)-zmmin(2))/2.d0
4535 c write(*,*)'bornes2=',exp(-zpmax(2)),exp(-zpmin(2))
4536 c &,exp(-zmmax(2)),exp(-zmmin(2)),xpmax(2),1.d0-exp(-zp(1))
4537 c &,1.d0-xpmin(3)-exp(-zp(1)),xpmin(2),1.d0-xpmax(3)-exp(-zp(1))
4542 zp(2)=zpmin(2)+dzp(2)*(1.d0+dble(2.*np2-3.)*dble(x1(jp2)))
4544 if(dzm(2).eq.0.d0)then
4550 if(dzm(2).ne.0.d0)then
4551 zm(2)=zmmin(2)+dzm(2)*(1.d0+dble(2.*nm2-3.)*dble(x1(jm2)))
4567 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4568 Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4571 & write(*,*)'j2=',k,exp(-zp(k)),exp(-zm(k))
4572 & ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2
4575 write(*,*)'error2 ?',dzp(k),dzm(k)
4578 if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4579 if(dzm(2).ne.0.d0)then
4581 &dble(a1(jm2))*PhiExpo(1.,sxp,sxm,sy,b)
4583 c &dble(a1(jm2))*PhiExact(1.,sxp,sxm,sy,b,mk)
4587 &dble(a1(jp2))*PhiExpo(1.,sxp,sxm,sy,b)
4589 c &dble(a1(jp2))*PhiExact(1.,sxp,sxm,sy,b,mk)
4592 c write(*,*)'m=2',mtmp,Gm2,Gp2,pGtab,sxp,sxm
4597 if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.3)then
4598 zpmin(3)=-log(min(xpmax,1.d0-exp(-zp(1))-exp(-zp(2))))
4599 zpmax(3)=-log(xpmin)
4600 dzp(3)=(zpmax(3)-zpmin(3))/2.d0
4604 if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.3)then
4605 zmmin(3)=-log(min(xmmax,1.d0-exp(-zm(1))-exp(-zm(2))))
4606 zmmax(3)=-log(xmmin)
4607 dzm(3)=(zmmax(3)-zmmin(3))/2.d0
4610 c write(*,*)'bornes3=',exp(-zpmax(3)),exp(-zpmin(3))
4611 c &,exp(-zmmax(3)),exp(-zmmin(3))
4616 zp(3)=zpmin(3)+dzp(3)*(1.d0+dble(2.*np3-3.)*dble(x1(jp3)))
4618 if(dzm(3).eq.0.d0)then
4624 if(dzm(3).ne.0.d0)then
4625 zm(3)=zmmin(3)+dzm(3)*(1.d0+dble(2.*nm3-3.)*dble(x1(jm3)))
4640 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4641 Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4644 & write(*,*)'j3=',k,exp(-zp(k)),exp(-zm(k))
4645 & ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2,jp3
4648 write(*,*)'error3 ?',k,dzp(k),dzm(k)
4651 if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4652 if(dzm(3).ne.0.d0)then
4654 &*dble(a1(jm3))*PhiExpo(1.,sxp,sxm,sy,b)
4658 &*dble(a1(jp3))*PhiExpo(1.,sxp,sxm,sy,b)
4664 if(dzm(3).ne.0.d0)Gp3=Gp3+Gm3*dble(a1(jp3))*exp(-zp(3))*dzm(3)
4669 if(mtmp.gt.2.and.dzm(2).ne.0.d0)then
4670 Gm2=Gm2+Gp3*dble(a1(jm2))*exp(-zm(2))*dzp(3)
4671 elseif(mtmp.gt.2)then
4672 Gp2=Gp2+Gp3*dble(a1(jp2))*exp(-zp(2))*dzp(3)
4676 if(dzm(2).ne.0.d0)Gp2=Gp2+Gm2*dble(a1(jp2))*exp(-zp(2))*dzm(2)
4681 if(mtmp.gt.1.and.dzm(1).ne.0.d0)then
4682 Gm1=Gm1+Gp2*dble(a1(jm1))*exp(-zm(1))*dzp(2)
4683 elseif(mtmp.gt.1)then
4684 Gp1=Gp1+Gp2*dble(a1(jp1))*exp(-zp(1))*dzp(2)
4688 if(dzm(1).ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*exp(-zp(1))*dzm(1)
4694 if(mtmp.gt.0)G0=Gp1*dzp(1)
4697 GammaGauss=GammaGauss+G0
4703 c-----------------------------------------------------------------------
4704 double precision function omWi(sy,b) !---check---
4705 c-----------------------------------------------------------------------
4706 c cut enhanced diagram integrated over xp, xm, xpr,xmr
4708 c b - impact parameter between the pomeron ends;
4710 c-----------------------------------------------------------------------
4712 include "epos.incpar"
4713 include "epos.incsem"
4714 include "epos.incems"
4715 common /psar7/ delx,alam3p,gam3p
4716 double precision xpmin,xmmin,zp,zm,alpp,alpm,xjacp,xjacm
4717 *,xpmax,xmmax,zpmin,zmmin,zpmax,chg
4718 double precision xp,xm,pGtab,omGam,dzp,Gp1,Gm1,xmin,eps
4719 *,sxp,sxm,PhiExpo,zmmax,dzm!,gamp,gamm,gampp,gammp
4721 common /ar3/ x1(7),a1(7)
4722 c double precision dgx1,dga1
4723 c common /dga20/ dgx1(10),dga1(10)
4731 gamb=gamD(1,iclpro,icltar)
4732 ctp060829 gamp=dble(gamb*b2/2.-alppar)
4733 ctp060829 gamm=dble(gamb*b2/2.-alppar)
4734 ctp060829 gampp=1.d0+2.d0*gamp
4735 ctp060829 gammp=1.d0+2.d0*gamm
4762 alpp=(1.d0+2.d0*dble(gamb*b2/2.))
4768 alpp=1.d0!(1.d0+2.d0*dble(gamb*b2/2.))
4776 alpm=(1.d0+2.d0*dble(gamb*b2/2.))
4782 alpm=1.d0!(1.d0+2.d0*dble(gamb*b2/2.))
4785 c write(*,*)'x+/-',intp,intm,xmmin,xmmax,xpmin,xpmax,alpp,alpm
4789 if(abs(xmmin-xmmax).ge.eps)then
4790 if(alpm.eq.0.d0)then
4797 dzm=(zmmax-zmmin)/2.d0
4801 if(abs(xpmin-xpmax).ge.eps)then
4802 if(alpp.eq.0.d0)then
4809 dzp=(zpmax-zpmin)/2.d0
4815 if(abs(dzp).gt.eps.and.abs(dzm).gt.eps)then
4816 c write(*,*)'Ca passe ...'
4820 zp=zpmin+dzp*(1.d0+dble(2.*np1-3.)*dble(x1(jp1)))
4821 c zp=zpmin+dzp*(1.d0+dble(2.*np1-3.)*dgx1(jp1))
4822 if(alpp.eq.0.d0)then
4827 xjacp=zp**(1.d0/alpp-1.d0)/alpp
4833 zm=zmmin+dzm*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1)))
4834 c zm=zmmin+dzm*(1.d0+dble(2.*nm1-3.)*dgx1(jm1))
4835 if(alpm.eq.0.d0)then
4840 xjacm=zm**(1.d0/alpm-1.d0)/alpm
4846 if(dzp.ge.0.d0.and.dzm.ge.0.d0)then
4847 pGtab=omGam(xp,xm,b)
4849 & write(*,*)'j1=',xp,xm,xmmin,xmmax,dzp,dzm,jp1
4852 write(*,*)'error ?',dzp,dzm
4854 if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4857 &dble(a1(jm1))*PhiExpo(1.,sxp,sxm,sy,b)*xjacm
4858 c &dga1(jm1)*PhiExpo(1.,sxp,sxm,sy,b)*xjacm
4859 c &dble(a1(jm1))*PhiExact(1.,sxp,sxm,sy,b)*xjacm
4862 &dble(a1(jp1))*PhiExpo(1.,sxp,sxm,sy,b)*xjacp
4863 c &dga1(jp1)*PhiExpo(1.,sxp,sxm,sy,b)*xjacp
4864 c &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b)*xjacp
4866 c write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm
4871 if(dzm.ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*dzm*xjacp
4872 c if(dzm.ne.0.d0)Gp1=Gp1+Gm1*dga1(jp1)*dzm*xjacp
4886 c-----------------------------------------------------------------------
4887 double precision function Womint(sy,bh) !---check---
4888 c-----------------------------------------------------------------------
4889 c - chi~(xp,xm)/2. for group of cut enhanced diagram giving
4890 c the same final state integrated over xpr and xmr (with ideal G)
4891 c bh - impact parameter between the pomeron ends;
4892 c xh - fraction of the energy squared s for the pomeron;
4893 c yp - rapidity for the pomeron;
4894 c-----------------------------------------------------------------------
4896 double precision omWi
4905 c-----------------------------------------------------------------------
4906 double precision function WomGamint(bh) !---check---
4907 c-----------------------------------------------------------------------
4908 c - chi~(xp,xm)/2. for group of integrated cut enhanced diagram giving
4909 c the same final for proposal.
4910 c bh - impact parameter between the pomeron ends;
4911 c xh - fraction of the energy squared s for the pomeron;
4912 c yp - rapidity for the pomeron;
4913 c-----------------------------------------------------------------------
4915 double precision omGamint
4917 WomGamint=omGamint(bh)