EPOS/epos167/epos-omg-160.f

   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
   5 c
   6 c               G = \sum_type  alp * xp**bet * xm**betp
   7 c
   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
  10 c
  11 c  Gpro:   bet  = betD'  + epsGp + gamD*b**2 + epsG -alppar
  12 c          bet' = betD'' + epsGt + gamD*b**2 + epsG -alppar
  13 c
  14 c          alp  = alpD*f * s**(betD+gamD*b**2+epsG) * exp(-b**2/delD)
  15 c
  16 c  Gtilde: bet~  = bet  + 1
  17 c          bet~' = bet' + 1
  18 c
  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')
  23 c
  24 c The parameters depend implicitely on s.
  25 c
  26 c In the program we use om1 = Gpro
  27 c  (they differ by a constant which is actually one)
  28 c and om5 = om1 * 0.5
  29 c All functions related to om1 are called om1... .
  30 c
  31 c The inclusive Pomeron distributions are
  32 c
  33 c      PomInc(xp,xm) = Gpro(xp,xm) * (1-xp)**alppro * (1-xm)**alptar
  34 c
  35 c----------------------------------------------------------------------
  36
  37
  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----------------------------------------------------------------------
  49       include 'epos.inc'
  50       include 'epos.incsem'
  51       include 'epos.incems'
  52       include 'epos.incpar'
  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
  57       parameter(nbkbin=40)
  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
  68       common/nucl3/phi,bimp
  69
  70       b1=0.03
  71       b2=bkmx*1.2
  72       db=(b2-b1)/nyeps
  73
  74       call utpri('GfunParK',ish,ishini,5)
  75
  76       cfbetpnp=0.d0
  77       cfbetppnp=0.d0
  78       cfbetpnm=0.d0
  79       cfbetppnm=0.d0
  80       cfalpro=dble(ucfpro)
  81       cfaltar=dble(ucftar)
  82       cfbpap=1.d0
  83       cfbppap=1.d0
  84       cfbpam=1.d0
  85       cfbppam=1.d0
  86       alpfom=0.
  87       sy=engy*engy
  88
  89       do k=1,koll
  90         do i=ntymi,ntymx
  91           atilde(i,k)=0.d0
  92           btildep(i,k)=0.d0
  93           btildepp(i,k)=0.d0
  94         enddo
  95         do i=idxD0,idxD1
  96           atildg(i,k)=0.d0
  97           btildgp(i,k)=0.d0
  98           btildgpp(i,k)=0.d0
  99         enddo
 100       enddo
 101
 102
 103 ! calculate collision number according to Glauber -----------------------
 104
 105       bglaub2=max(1.,sigine/10./pi)        !10= fm^2 -> mb
 106       bglaub=sqrt(bglaub2)
 107       nglevt=0
 108       do ko=1,koll
 109         r=bk(ko)
 110         if(r.le.bglaub)nglevt=nglevt+1
 111       enddo
 112
 113 ! Z_parton_projectile (zparpro) and Z_parton_target (zpartar)-----------
 114
 115       if(iscreen.eq.1.or.isplit.eq.1)then
 116
 117         bcut=0. !zbcut
 118 ctp060829        bzero=bglaub
 119         b2x=epscrp*epscrp*bglaub2
 120         alpfom=alpfomi*epscrw*fscra(engy/egyscr)
 121         ztav=0.
 122         zpav=0.
 123
 124         do k=1,koll
 125           ip=iproj(k)
 126           it=itarg(k)
 127         !....... targ partons seen by proj
 128           if(lproj(ip).ge.1)then
 129             absb=max(1.e-9,bk(k))
 130             b2=absb*absb
 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
 137               zpartar(k)=zkp
 138               do lt=1,lproj3(ip)
 139                 kp=kproj3(ip,lt)
 140                 if(kp.ne.k)then
 141                   absb=max(1.e-9,abs(bk(kp))-bcut)
 142                   b2=absb*absb
 143                   zkp=epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr)
 144                 endif
 145                 zpartar(k)=zpartar(k)+zkp
 146 c                alpfom=max(alpfom,dble(zpartar(k)))
 147               enddo
 148             endif
 149           else
 150             zpartar(k)=0.
 151           endif
 152           ztav=ztav+zpartar(k)
 153          !...........proj partons seen by targ
 154           if(ltarg(it).ge.1)then
 155             absb=max(1.e-9,bk(k))
 156             b2=absb*absb
 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
 163               zparpro(k)=zkt
 164               do lp=1,ltarg3(it)
 165                 kt=ktarg3(it,lp)
 166                 if(kt.ne.k)then
 167                   absb=max(1.e-9,abs(bk(kt))-bcut)
 168                   b2=absb*absb
 169                   zkt=epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr)
 170                 endif
 171                 zparpro(k)=zparpro(k)+zkt
 172 c                alpfom=max(alpfom,dble(zparpro(k)))
 173               enddo
 174             endif
 175           else
 176             zparpro(k)=0.
 177           endif
 178           zpav=zpav+zparpro(k)
 179         enddo
 180
 181       else                      ! no screening
 182
 183         do k=1,koll
 184           zparpro(k)=0.
 185           zpartar(k)=0.
 186         enddo
 187
 188       endif
 189
 190 c calculation of epsilongp epsilongt
 191
 192       if(iscreen.eq.1)then
 193
 194           !ip=0
 195        do k=1,koll
 196           !ipp=ip
 197         epsG=0.
 198         epsilongs(k,0)=epsG
 199         epsilongs(k,1)=epsG
 200         ip=iproj(k)             !...........projectile
 201         if(lproj(ip).gt.0)then
 202           x=zparpro(k)
 203           epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
 204      &                        min(epscrx,epscrs*x))
 205 c     &                        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)))
 212           gammaV(k)=1.d0
 213           else
 214          epsilongp(k,0)=0.
 215          epsilongp(k,1)=0.
 216          gammaV(k)=1.d0
 217         endif
 218
 219         it=itarg(k)             !...........target
 220         if(ltarg(it).gt.0)then
 221           x=zpartar(k)
 222           epsilongt(k,0)=max(-betDpp(0,iclpro,icltar)-0.95+alppar,
 223      &                        min(epscrx,epscrs*x))
 224 c     &                        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)
 232         else
 233          epsilongt(k,0)=0.
 234          epsilongt(k,1)=0.
 235          gammaV(k)=gammaV(k)
 236         endif
 237
 238        enddo
 239
 240       else                      ! no screening
 241
 242        do k=1,koll
 243         epsilongs(k,0)=0.
 244         epsilongs(k,1)=0.
 245         epsilongp(k,0)=0.
 246         epsilongp(k,1)=0.
 247         epsilongt(k,0)=0.
 248         epsilongt(k,1)=0.
 249         gammaV(k)=1.d0
 250        enddo
 251
 252       endif
 253
 254
 255 !..............alpha beta betap for Gtilde (->PhiExpo).......................
 256
 257       imax=idxD1
 258       if(iomega.eq.2)imax=1
 259
 260       do k=1,koll
 261
 262         b=bk(k)
 263         b2=bk(k)*bk(k)
 264         ip=iproj(k)
 265         it=itarg(k)
 266
 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))
 281           else
 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))
 286           endif
 287         else
 288           xkappa=1.
 289         endif
 290         gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
 291         gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
 292
 293         do i=idxDmin,imax
 294           gamV=gammaV(k)
 295           if(i.lt.2)then
 296             epsG=epsilongs(k,i)
 297           else
 298             epsG=0.
 299           endif
 300           gamb=gamD(i,iclpro,icltar)*b2+epsG
 301           atildg(i,k)=dble(alpD(i,iclpro,icltar))
 302      *            *cfalpro*cfaltar
 303      *            *gamv
 304 c          if(i.eq.0) atildg(i,k)=atildg(i,k)
 305           atildg(i,k)=atildg(i,k)
 306      *            *dble(xkappa*xkappa)
 307           if(i.lt.2)then
 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
 313             epsGp=epsilongp(k,i)
 314             epsGt=epsilongt(k,i)
 315             btildgp(i,k)=dble(betDp(i,iclpro,icltar)
 316      *                    +epsGp
 317      *                    +gamb-alppar)+1.d0
 318             btildgpp(i,k)=dble(betDpp(i,iclpro,icltar)
 319      *                    +epsGt
 320      *                    +gamb-alppar)+1.d0
 321           else
 322             absb=abs(bk(k))-bmxdif(iclpro,icltar)
 323             b2a=absb*absb
 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
 329           endif
 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))
 334      *            /coefgdp/coefgdt
 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
 346      *              +epsGp+gamb
 347            w(6+4*i,kk)=w(6+4*i,kk)
 348      *              +betDpp(i,iclpro,icltar)  !tgt eff
 349      *              +epsGt+gamb
 350            w(7+4*i,kk)=w(7+4*i,kk)
 351      *              +betDp(i,iclpro,icltar)  !prj unscr
 352      *              +gamb
 353            w(8+4*i,kk)=w(8+4*i,kk)
 354      *              +betDpp(i,iclpro,icltar) !tgt unscr
 355      *              +gamb
 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)
 358           endif
 359    !................
 360         enddo
 361       enddo
 362
 363 !...........................................................................
 364
 365            zppevt=zpav/koll
 366            zptevt=ztav/koll
 367            ktot=int(bimp)+1
 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
 372            endif
 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
 377             w(20,n)=w(20,n)+1
 378            endif
 379
 380
 381 !........alpha beta betap for Gpro...........................................
 382
 383       if(irea.le.0)then
 384       do k=1,koll
 385         ip=iproj(k)
 386         it=itarg(k)
 387
 388         b2=bk(k)*bk(k)
 389         imax=ntymx
 390         if(iomega.eq.2)imax=1
 391         do i=ntymin,imax
 392
 393           if(i.lt.2)then
 394             epsG=epsilongs(k,i)
 395           else
 396             epsG=0.
 397           endif
 398           gamb=gamD(i,iclpro,icltar)*b2+epsG
 399
 400           atilde(i,k)=dble(alpD(i,iclpro,icltar))
 401           if(i.lt.2)then
 402           atilde(i,k)=atilde(i,k)*dble(
 403      *              exp(-b2/delD(i,iclpro,icltar))
 404      *              *sy**(betD(i,iclpro,icltar)
 405      *                    +gamb)
 406      *              *chad(iclpro)*chad(icltar))
 407             epsGp=epsilongp(k,i)
 408             epsGt=epsilongt(k,i)
 409             btildep(i,k)=dble(betDp(i,iclpro,icltar)
 410      *                    +epsGp
 411      *                    +gamb-alppar)
 412             btildepp(i,k)=dble(betDpp(i,iclpro,icltar)
 413      *                    +epsGt
 414      *                    +gamb-alppar)
 415           else
 416             absb=abs(bk(k))-bmxdif(iclpro,icltar)
 417             b2a=absb*absb
 418             atilde(i,k)=atilde(i,k)*dble(
 419      *                  sy**betD(i,iclpro,icltar)
 420      *                 *exp(-b2a/delD(i,iclpro,icltar)))
 421
 422             btildep(i,k)=dble(betDp(i,iclpro,icltar)-alppar)
 423             btildepp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar)
 424           endif
 425
 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&')
 428         enddo
 429
 430       enddo
 431       endif
 432
 433 !...........................................................................
 434
 435       if(ish.ge.8)then
 436       do k=1,koll
 437         ip=iproj(k)
 438         it=itarg(k)
 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 "
 445         do i=ntymin,ntymx
 446         write(ifch,*)i,atilde(i,k),btildep(i,k),btildepp(i,k)
 447         enddo
 448       enddo
 449       endif
 450
 451       call utprix('GfunParK',ish,ishini,5)
 452
 453       return
 454       end
 455
 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----------------------------------------------------------------------
 466
 467       include 'epos.inc'
 468       include 'epos.incsem'
 469       include 'epos.incpar'
 470       include 'epos.incems'
 471       parameter(nbkbin=40)
 472       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
 473       common /kwrite/ xkapZ
 474       parameter(eps=1.e-20)
 475       double precision utgam2
 476
 477       call utprj('GfunPar ',ish,ishini,6)
 478
 479    !???????????????????????????????????????????????????estimation
 480       x=sqrt(spp)
 481       if(sigine.eq.0.)then
 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
 488         else
 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
 491         endif
 492         siginex=sigtotf-sigelaf
 493       else
 494        siginex=sigine
 495       endif
 496       bglaub2=max(1.,siginex/10./pi)
 497       bglaub=sqrt(bglaub2)
 498       if(ish.ge.6)write(ifch,*)'Gf in:',m,i,x,b,sigine,siginex,bglaub
 499      &                                 ,iclpro,icltar
 500       !!!!!print*,'++++++++++++++++++',x,siginex,bglaub
 501       b2=b*b
 502       ee=sqrt(spp)
 503       cfalpro=ucfpro
 504       cfaltar=ucftar
 505       gamb=gamD(i,iclpro,icltar)*b2
 506
 507       if(iscreen.ne.0)then
 508         bcut=0.!zbcut
 509 ctp060829        bzero=bglaub
 510         absb=max(1.e-9,abs(b)-bcut)
 511         b2a=absb*absb
 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)
 515         x=zzp
 516         if(i.eq.0)then
 517           epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar
 518      &             ,min(epscrx,epscrs*x))
 519 c     &             ,epscrs*x)
 520 c     &             ,epscrs*x*abs(betDp(i,iclpro,icltar)))
 521 c     &             ,epscrs*min(epscrx,x*abs(betDp(i,iclpro,icltar))))
 522 c          epsG=epsG+epsGp
 523         elseif(i.eq.1)then
 524           epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar
 525      &             ,epscrh*min(epscrx,x*abs(betDp(i,iclpro,icltar))))
 526         else
 527           epsGp=0.
 528         endif
 529         gamV=1.
 530          x=zzt
 531         if(i.eq.0)then
 532           epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar
 533      &             ,min(epscrx,epscrs*x))
 534 c     &             ,epscrs*x)
 535 c     &             ,epscrs*x*abs(betDpp(i,iclpro,icltar)))
 536 c     &             ,epscrs*min(epscrx,x*abs(betDpp(i,iclpro,icltar))))
 537 c          epsG=epsG+epsGt
 538         elseif(i.eq.1)then
 539           epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar
 540      &             ,epscrh*min(epscrx,x*abs(betDpp(i,iclpro,icltar))))
 541         else
 542           epsGt=0.
 543         endif
 544 c        gamV=gamV
 545         epsG=0.
 546       else
 547         zzp=0.
 548         zzt=0.
 549         epsGp=0.
 550         epsGt=0.
 551         epsG=0.
 552         gamV=1.
 553       endif
 554
 555       gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
 556       gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
 557
 558       rho=betD(i,iclpro,icltar)
 559      *   +gamb+epsG
 560
 561
 562       if(m.eq.1)then
 563
 564         alp=alpD(i,iclpro,icltar)
 565         if(i.lt.2)then
 566           alp=alp
 567      *       *exp(rho*log(spp)-b2/delD(i,iclpro,icltar))
 568           bet=betDp(i,iclpro,icltar)
 569      *         +epsGp
 570      *         +gamb-alppar
 571           betp=betDpp(i,iclpro,icltar)
 572      *         +epsGt
 573      *         +gamb-alppar
 574         else
 575           absb=abs(b)-bmxdif(iclpro,icltar)
 576           b2a=absb*absb
 577           alp=alp
 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
 582         endif
 583
 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&')
 587         endif
 588
 589       elseif(m.eq.2)then
 590         xkappa=1.
 591 c        if(i.eq.0.and.b.lt.(nbkbin-1)*bkbin)then
 592         if(b.lt.(nbkbin-1)*bkbin)then
 593           ibk=int(b/bkbin)+1
 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))
 606           else
 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))
 611           endif
 612           xkapZ=xkappa
 613         endif
 614
 615         alp=alpD(i,iclpro,icltar)
 616      *        *cfalpro*cfaltar
 617      *        *gamV
 618 c        if(i.eq.0)alp=alp
 619         alp=alp
 620      *        *xkappa*xkappa
 621
 622         if(i.lt.2)then
 623           alp=alp
 624      *     *exp(rho*log(spp)-b2/delD(i,iclpro,icltar))
 625      *        *chad(iclpro)*chad(icltar)
 626      *        *gfactorp *gfactort
 627           bet=betDp(i,iclpro,icltar)
 628      *        +epsGp
 629      *        +gamb-alppar+1.
 630           betp=betDpp(i,iclpro,icltar)
 631      *        +epsGt
 632      *        +gamb-alppar+1.
 633         else
 634           absb=abs(b)-bmxdif(iclpro,icltar)
 635           b2a=absb*absb
 636           alp=alp
 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.
 641         endif
 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
 645       else
 646
 647         stop'GproPar: wrong m value.              '
 648
 649       endif
 650
 651       epsp=epsGp
 652       epst=epsGt
 653       epss=epsG
 654       gamvv=gamV
 655
 656       alpUni(i,m)=alp
 657       betUni(i,m)=bet
 658       betpUni(i,m)=betp
 659
 660
 661       if(ish.ge.6)write(ifch,*)'   GfunPar :',alp,bet,betp,epsp,epst
 662      &                                    ,epss,gamvv
 663
 664       call utprjx('GfunPar ',ish,ishini,6)
 665       end
 666
 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----------------------------------------------------------------------
 674
 675       include 'epos.inc'
 676       include 'epos.incsem'
 677       include 'epos.incpar'
 678       include 'epos.incems'
 679       parameter(eps=1.e-20)
 680
 681       call utprj('GfomPar ',ish,ishini,6)
 682
 683    !???????????????????????????????????????????????????estimation
 684       x=sqrt(spp)
 685       if(sigine.eq.0.)then
 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
 692         else
 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
 695         endif
 696         siginex=sigtotf-sigelaf
 697       else
 698        siginex=sigine
 699       endif
 700       bglaub2=max(1.,siginex/10./pi)
 701       bglaub=sqrt(bglaub2)
 702       b2=b*b
 703       ee=sqrt(spp)
 704
 705       if(iscreen.ne.0)then
 706         bcut=0.!zbcut
 707 ctp060829        bzero=bglaub
 708         absb=max(1.e-9,abs(b)-bcut)
 709         b2a=absb*absb
 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)
 713       else
 714         zzp=0.
 715         zzt=0.
 716       endif
 717
 718
 719       z0=alpfomi!*epscrw*fscra(ee/egyscr)
 720       if(z0.gt.0.)then
 721         z1=zzp
 722         zpUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
 723 c      zpUni=dble(4.*z0*(z1/z0)**1.5)
 724         z1=zzt
 725         ztUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
 726 c      ztUni=dble(4.*z0*(z1/z0)**1.5)
 727       else
 728         zpUni=0d0
 729         ztUni=0d0
 730       endif
 731
 732       if(ish.ge.6)write(ifch,*)'   GfomPar :',zpUni,ztUni
 733
 734       call utprjx('GfomPar ',ish,ishini,6)
 735       end
 736
 737 c----------------------------------------------------------------------
 738       function fscra(x)
 739 c----------------------------------------------------------------------
 740       fscra=0
 741       x2=x*x
 742       if(x2.gt.1.)fscra=log(x2)
 743       end
 744
 745 c----------------------------------------------------------------------
 746       function fscro(x)
 747 c----------------------------------------------------------------------
 748       include 'epos.incpar'
 749       fscro=epscrv
 750       x2=x*x
 751       if(x2.gt.1.)fscro=sqrt(log(x2)**2+epscrv**2)
 752       end
 753
 754 c----------------------------------------------------------------------
 755       subroutine recalcZPtn !???????????? not updated !!!
 756 c----------------------------------------------------------------------
 757
 758       include 'epos.inc'
 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
 762 c       npom=nprt(1)
 763 c       k=1
 764 c       ip=iproj(1)
 765 c       it=itarg(1)
 766 c       zparpro(k)=max(0,npom-1)*0.2
 767 c       zpartar(k)=0
 768 c       zpartar(k)=max(0,npom-1)*0.2
 769 c       ztav=zpartar(k)
 770 c       zpav=zparpro(k)
 771 c       zppevt=zpav
 772 c       zptevt=ztav
 773 c      endif
 774       end
 775
 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----------------------------------------------------------------------
 784       include 'epos.inc'
 785       include 'epos.incsem'
 786       include 'epos.incpar'
 787
 788       double precision Gf,xp,xm,xh,yp
 789
 790       Gf=0.d0
 791       xp=sqrt(xh)*exp(yp)
 792       xm=xh/xp
 793       spp=engy**2
 794       imax=idxD1
 795       if(iomega.eq.2)imax=1
 796       do i=idxDmin,imax
 797         call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
 798         Gf=Gf+alp*real(xp)**bet*real(xm)**betp
 799       enddo
 800       om1=Gf
 801      *  * dble(chad(iclpro)*chad(icltar))
 802       end
 803
 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----------------------------------------------------------------------
 810       include 'epos.inc'
 811       include 'epos.incsem'
 812       include 'epos.incpar'
 813       parameter(eps=1.e-20)
 814
 815       spp=engy*engy
 816       imax=idxD1
 817       if(iomega.eq.2)imax=1
 818       cint=0
 819       do i=idxDmin,imax
 820         call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
 821         cint2=gamv*alp
 822         if((bet+1.0).gt.eps)then
 823           cint2=cint2/(bet+1.0)
 824         else
 825           cint2=-cint2*log(xminDf)
 826         endif
 827         if((betp+1.0).gt.eps)then
 828           cint2=cint2/(betp+1.0)
 829         else
 830           cint2=-cint2*log(xminDf)
 831         endif
 832         cint=cint+cint2
 833       enddo
 834
 835       if(cint.lt.-eps)then
 836         write(*,*) 'WARNING ! om1intb in epos-omg is <0 !!!!!'
 837         write(*,*) 'WARNING ! => om1intb set to 1e-3 !!!!!'
 838         write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
 839         cint=1.e-3
 840       endif
 841
 842       om1intb=cint
 843      *       *chad(iclpro)*chad(icltar)
 844
 845       return
 846       end
 847
 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----------------------------------------------------------------------
 854       include 'epos.inc'
 855       include 'epos.incsem'
 856       include 'epos.incpar'
 857       parameter(eps=1.e-20)
 858
 859       spp=engy*engy
 860       call GfunPar(1,iqq,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
 861       cint=gamv*alp
 862       if((bet+1.0).gt.eps)then
 863         cint=cint/(bet+1.0)
 864       else
 865         cint=-cint*log(xminDf)
 866       endif
 867       if((betp+1.0).gt.eps)then
 868         cint=cint/(betp+1.0)
 869       else
 870         cint=-cint*log(xminDf)
 871       endif
 872       if(cint.lt.-eps)then
 873         write(*,*) 'WARNING ! om1intbi in epos-omg is <0 !!!!!'
 874         write(*,*) 'WARNING ! => om1intbi set to 1e-3 !!!!!'
 875         write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
 876         cint=1.e-3
 877       endif
 878
 879       om1intbi=cint
 880      *       *chad(iclpro)*chad(icltar)
 881
 882       return
 883       end
 884
 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----------------------------------------------------------------------
 891       include 'epos.inc'
 892       include 'epos.incems'
 893       include 'epos.incsem'
 894       double precision cint,gamom,deltap,deltam
 895       double precision utgam2,Fp,Fm
 896
 897       spp=engy**2
 898       om1intbc=0.d0
 899       Fp=dble(ucfpro)   !gamma(1+alplea)
 900       Fm=dble(ucftar)
 901
 902       imax=idxD1
 903       if(iomega.eq.2)imax=1
 904
 905       cint=0.d0
 906
 907       do i=idxDmin,imax
 908         call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
 909        gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
 910         deltap=dble(bet)
 911         deltam=dble(betp)
 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)))
 915       enddo
 916
 917       om1intbc=cint*Fp*Fm
 918
 919       if(om1intbc.lt.-1.d-10)then
 920         write(*,*) 'WARNING ! om1intbc in epos-omg is <0 !!!!!'
 921         write(*,*) 'WARNING ! => om1intbc set to 0. !!!!!'
 922         om1intbc=0.d0
 923       endif
 924
 925       return
 926       end
 927
 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----------------------------------------------------------------------
 934 c      include 'epos.inc'
 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
 939 c
 940 c      spp=engy**2
 941 c      om1intbci=0.d0
 942 c      Fp=dble(ucfpro)   !gamma(1+alplea)
 943 c      Fm=dble(ucftar)
 944 c
 945 c      i=iqq
 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)
 948 c      deltap=dble(bet)
 949 c      deltam=dble(betp)
 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)))
 953 c
 954 c      om1intbci=cint*Fp*Fm
 955 c
 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. !!!!!'
 959 c        om1intbci=0.d0
 960 c      endif
 961 c
 962 c      return
 963 c      end
 964 c
 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----------------------------------------------------------------------
 971       include 'epos.inc'
 972       include 'epos.incems'
 973       include 'epos.incsem'
 974       double precision cint,gamom,deltap,deltam,xprem,xmrem
 975
 976       om1intgck=0.d0
 977
 978       imax=idxD1
 979       if(iomega.eq.2)imax=1
 980
 981       cint=0.d0
 982       do i=idxDmin,imax
 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)
 990       enddo
 991       om1intgck=cint
 992
 993       return
 994       end
 995
 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----------------------------------------------------------------------
1002       include 'epos.inc'
1003       include 'epos.incpar'
1004       include 'epos.incsem'
1005       double precision cint,gamom,deltap,deltam,eps
1006       parameter(eps=1.d-20)
1007
1008       spp=engy**2
1009       om1intgc=0.d0
1010
1011       imax=idxD1
1012       if(iomega.eq.2)imax=1
1013
1014
1015
1016       cint=0.d0
1017
1018       do i=idxDmin,imax
1019         call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
1020         gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
1021         deltap=dble(bet)
1022         deltam=dble(betp)
1023         if((deltap+1.d0).gt.eps)then
1024           gamom=gamom/(deltap+1.0)
1025         else
1026           gamom=-gamom*log(xminDf)
1027         endif
1028         if((deltam+1.d0).gt.eps)then
1029           gamom=gamom/(deltam+1.0)
1030         else
1031           gamom=-gamom*log(xminDf)
1032         endif
1033         cint=cint+gamom /(2.d0+deltap) /(2.d0+deltam)
1034       enddo
1035       om1intgc=cint
1036
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. !!!!!'
1041         om1intgc=0.d0
1042       endif
1043
1044       return
1045       end
1046
1047
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----------------------------------------------------------------------
1056       include 'epos.inc'
1057       include 'epos.incems'
1058       include 'epos.incsem'
1059       double precision om1intb,om1intgc
1060
1061       call utpri('intom1',ish,ishini,5)
1062
1063       if(irea.eq.0)then
1064
1065
1066       do k=1,koll
1067
1068         om1int(k)=om1intb(bk(k))
1069         om1intc(k)=om1intgc(bk(k))
1070
1071       enddo
1072
1073       endif
1074       call utprix('intom1',ish,ishini,5)
1075
1076
1077       return
1078       end
1079
1080
1081 c----------------------------------------------------------------------
1082       double precision function om1xpk(xp,xpr1i,xmr1i,k)   !---test---
1083 c----------------------------------------------------------------------
1084 c \int dxm om1*(1-xp)*(1-xm)   (normalized)
1085 c k - pair indice;
1086 c----------------------------------------------------------------------
1087       include 'epos.inc'
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)
1094
1095       om1xpk=0.d0
1096       if(xp.ge.xpr1i)return
1097
1098       xpr1=1.d0
1099       xmr1=1.d0
1100       imin=ntymin
1101       imax=ntymx
1102       if(iomega.eq.2)imax=1
1103
1104       do i=imin,imax
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)
1110           else
1111             gamomx(i)=-gamomx(i)*log(xminDf)
1112           endif
1113       enddo
1114
1115
1116       cint=0.d0
1117       do i=imin,imax
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)
1121           else
1122             gamom=-gamom*log(xminDf)
1123           endif
1124           cint=cint+gamom
1125       enddo
1126
1127
1128       do i=imin,imax
1129         om1xpk=om1xpk+gamomx(i)*xp**deltap(i)
1130      &                       *(xpr1-xp)
1131
1132       enddo
1133
1134       om1xpk=om1xpk/cint
1135
1136       return
1137       end
1138
1139
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)
1144 c k - pair indice;
1145 c----------------------------------------------------------------------
1146       include 'epos.inc'
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)
1152
1153       om1xmk=0.d0
1154       if(xp.ge.xpr1i)return
1155       if(xm.ge.xmr1i)return
1156       xpr1=1.d0
1157       xmr1=1.d0
1158
1159       imin=ntymin
1160       imax=ntymx
1161       if(iomega.eq.2)imax=1
1162
1163       do i=imin,imax
1164           gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xpr1-xp)
1165           deltam(i)=btildepp(i,k)
1166       enddo
1167
1168       cint=0.d0
1169       do i=imin,imax
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)
1173           else
1174             gamom=-gamom*log(xminDf)
1175           endif
1176           cint=cint+gamom
1177       enddo
1178
1179
1180       do i=imin,imax
1181         om1xmk=om1xmk+gamomx(i)*xm**deltam(i)*(xmr1-xm)
1182       enddo
1183
1184       om1xmk=om1xmk/cint
1185
1186       return
1187       end
1188
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.
1194 c k - pair indice;
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----------------------------------------------------------------------
1198       include 'epos.inc'
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
1206
1207
1208       om1xprk=0.d0
1209       x0=log(xminDf)
1210       x1=log(xpremi)
1211       xprem=1.d0
1212       xmrem=1.d0
1213       imin=ntymin
1214       imax=ntymx
1215       lxmin=1.d0
1216       if(iomega.eq.2)imax=1
1217
1218       do i=imin,imax
1219         if(ir.gt.0)then
1220           deltap(i)=btildep(i,k)
1221           deltam(i)=btildepp(i,k)
1222         else
1223           deltap(i)=btildepp(i,k)
1224           deltam(i)=btildep(i,k)
1225         endif
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)
1229         else
1230           lxmin=log(xminDf)
1231           gamomx(i)=-gamomx(i)*lxmin
1232         endif
1233       enddo
1234
1235       f0=0.d0
1236       f1=0.d0
1237       do i=imin,imax
1238         if((deltap(i)+1.d0).gt.eps)then
1239           f0=f0+gamomx(i)
1240      &          *(xprem*exp(x0)**(1.d0+deltap(i))/(1.d0+deltap(i))
1241      &           -exp(x0)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1242           f1=f1+gamomx(i)
1243      &          *(xprem*exp(x1)**(1.d0+deltap(i))/(1.d0+deltap(i))
1244      &           -exp(x1)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1245         else
1246           lxmin=log(xminDf)
1247           f0=f0+gamomx(i)*(xprem*(x0-lxmin)-exp(x0)+xminDf)
1248           f1=f1+gamomx(i)*(xprem*(x1-lxmin)-exp(x1)+xminDf)
1249         endif
1250       enddo
1251       f00=f0
1252       cint=f1-f00
1253       f0=-(f0-f00)/cint
1254       f1=-(f1-f00)/cint
1255       ntry=0
1256  11   ntry=ntry+1
1257       r=drangen(dble(ntry))
1258       f0t=f0+r
1259       f1t=f1+r
1260       if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1261       if(f1t*f0t.ge.eps)then
1262         do i=imin,imax
1263          write(ifmt,*)i,gamomx(i),deltap(i),deltam(i)
1264         enddo
1265         write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k
1266         call utstop('om1xprk (2)&')
1267       endif
1268       f0=f0t
1269       f1=f1t
1270       if(abs(f0).le.eps) then
1271         om1xprk=exp(x0)
1272         return
1273       endif
1274       if(abs(f1).le.eps) then
1275         om1xprk=exp(x1)
1276         return
1277       endif
1278       x=0.5d0*(x1+x0)
1279       deltx=abs(x1-x0)
1280
1281
1282       ntry=0
1283
1284  111  continue
1285       if(ntry.le.1000)then
1286       fx=0.d0
1287       fpx=0.d0
1288       do i=imin,imax
1289         if((deltap(i)+1.d0).gt.eps)then
1290           fx=fx+gamomx(i)
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))
1294         else
1295           fx=fx+gamomx(i)*(xprem*(x-lxmin)-exp(x)+xminDf)
1296           fpx=fpx+gamomx(i)*(xprem/exp(x)-1.d0)
1297         endif
1298       enddo
1299       fx=-(fx-f00)/cint+r
1300       fpx=fpx/cint
1301       xt=x-fx/fpx
1302
1303       if (f0*fx.lt.0.D0) then
1304         f1=fx
1305         x1=x
1306       else
1307         f0=fx
1308         x0=x
1309       endif
1310       if ((xt.lt.x0.or.xt.gt.x1).and.abs(f1-f0).gt.eps) then
1311         xt=x1-f1*(x1-x0)/(f1-f0)
1312       endif
1313
1314        else
1315
1316         write(ifmt,*)'Warning in om1xprk, to much try !'
1317
1318       endif
1319
1320
1321       if(abs(x-xt).gt.deltx*0.5d0) then
1322         xt=(x1+x0)*0.5D0
1323       endif
1324       deltx=abs(x-xt)
1325       if(abs(x).gt.eps)then
1326         prec=abs(deltx/x)
1327       else
1328         call utstop('Problem in om1xprk&')
1329       endif
1330
1331       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1332          x=xt
1333          ntry=ntry+1
1334          goto 111
1335       endif
1336
1337       om1xprk=exp(x)
1338
1339       return
1340       end
1341
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.
1347 c k - pair indice;
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----------------------------------------------------------------------
1351       include 'epos.inc'
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)
1359
1360       om1xmrk=0.d0
1361       if(xp.ge.xpremi)return
1362
1363       xprem=1.d0
1364       xmrem=1.d0
1365       x0=log(xminDf)
1366       x1=log(xmremi)
1367       imin=ntymin
1368       imax=ntymx
1369       if(iomega.eq.2)imax=1
1370
1371       do i=imin,imax
1372         if(ir.gt.0)then
1373           gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xprem-xp)
1374           deltam(i)=btildepp(i,k)
1375         else
1376           gamomx(i)=atilde(i,k)*xp**btildepp(i,k)*(xprem-xp)
1377           deltam(i)=btildep(i,k)
1378         endif
1379       enddo
1380
1381
1382
1383       f0=0.d0
1384       f1=0.d0
1385       do i=imin,imax
1386         if(abs(deltam(i)+1.d0).gt.eps)then
1387           f0=f0+gamomx(i)
1388      &          *(xmrem*exp(x0)**(1.d0+deltam(i))/(1.d0+deltam(i))
1389      &           -exp(x0)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1390           f1=f1+gamomx(i)
1391      &          *(xmrem*exp(x1)**(1.d0+deltam(i))/(1.d0+deltam(i))
1392      &           -exp(x1)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1393         else
1394           lxmin=log(xminDf)
1395           f0=f0+gamomx(i)*(xmrem*(x0-lxmin)-exp(x0)+xminDf)
1396           f1=f1+gamomx(i)*(xmrem*(x1-lxmin)-exp(x1)+xminDf)
1397         endif
1398       enddo
1399       f00=f0
1400       cint=f1-f00
1401       f0=-(f0-f00)/cint
1402       f1=-(f1-f00)/cint
1403       ntry=0
1404  11   ntry=ntry+1
1405       r=drangen(dble(ntry))
1406       f0t=f0+r
1407       f1t=f1+r
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&')
1412       endif
1413       f0=f0t
1414       f1=f1t
1415       if(abs(f0).lt.eps) then
1416         om1xmrk=exp(x0)
1417         return
1418       endif
1419       if(abs(f1).lt.eps) then
1420         om1xmrk=exp(x1)
1421         return
1422       endif
1423       x=0.5d0*(x1+x0)
1424       deltx=abs(x1-x0)
1425
1426
1427       ntry=0
1428
1429  111  continue
1430       if(ntry.le.1000)then
1431       fx=0.d0
1432       fpx=0.d0
1433       do i=imin,imax
1434         if(abs(deltam(i)+1.d0).gt.eps)then
1435           fx=fx+gamomx(i)
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))
1439         else
1440           fx=fx+gamomx(i)*(xmrem*(x-lxmin)-exp(x)+xminDf)
1441           fpx=fpx+gamomx(i)*(xmrem/exp(x)-1.d0)
1442         endif
1443       enddo
1444       fx=-(fx-f00)/cint+r
1445       fpx=fpx/cint
1446       xt=x-fx/fpx
1447
1448       if (f0*fx.lt.-eps) then
1449         f1=fx
1450         x1=x
1451       else
1452         f0=fx
1453         x0=x
1454       endif
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)
1457       endif
1458
1459        else
1460
1461         write(ifmt,*)'Warning in om1xmrk, to much try !'
1462
1463       endif
1464
1465       if(abs(x-xt).gt.deltx*0.5d0) then
1466         xt=(x1+x0)*0.5D0
1467       endif
1468       deltx=abs(x-xt)
1469       if(abs(x).gt.eps)then
1470         prec=abs(deltx/x)
1471       else
1472         call utstop('Problem in om1xmrk&')
1473       endif
1474
1475       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1476          x=xt
1477          ntry=ntry+1
1478          goto 111
1479       endif
1480
1481       om1xmrk=exp(x)
1482
1483       return
1484       end
1485
1486 c----------------------------------------------------------------------
1487       double precision function om1xk(xh,k)   !---test---
1488 c----------------------------------------------------------------------
1489 c \int dxp om1   (normalised)
1490 c k - pair indice;
1491 c----------------------------------------------------------------------
1492       include 'epos.inc'
1493       include 'epos.incpar'
1494       include 'epos.incems'
1495
1496       double precision xh,gamomx(ntymi:ntymx),cint,alpp(ntymi:ntymx)
1497      &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps
1498      &,gamom
1499       parameter(eps=1.d-20)
1500
1501
1502       om1xk=0.d0
1503
1504       imin=ntymin
1505       imax=ntymx
1506       if(iomega.eq.2)imax=1
1507
1508       do i=imin,imax
1509         gamomx(i)=atilde(i,k)
1510         deltap(i)=btildep(i,k)
1511         deltam(i)=btildepp(i,k)
1512
1513         delta(i)=(deltap(i)+deltam(i))*0.5d0
1514         alpp(i)=deltap(i)-deltam(i)
1515
1516       enddo
1517
1518       cint=0.d0
1519       do i=imin,imax
1520         gamom=gamomx(i)
1521         if((deltap(i)+1.d0).gt.eps)then
1522           gamom=gamom/(deltap(i)+1.d0)
1523         else
1524           gamom=-gamom*log(xminDf)
1525         endif
1526         if((deltam(i)+1.d0).gt.eps)then
1527           gamom=gamom/(deltam(i)+1.d0)
1528         else
1529           gamom=-gamom*log(xminDf)
1530         endif
1531         cint=cint+gamom
1532       enddo
1533
1534
1535       do i=imin,imax
1536         if(abs(alpp(i)).gt.eps)then
1537         om1xk=om1xk+gamomx(i)/alpp(i)*xh**deltam(i)*(1.d0-xh**alpp(i))
1538         else
1539         om1xk=om1xk-gamomx(i)*xh**delta(i)*log(xh)
1540         endif
1541       enddo
1542
1543       om1xk=om1xk/cint
1544
1545       return
1546       end
1547
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;
1553 c k - pair indice;
1554 c----------------------------------------------------------------------
1555       include 'epos.inc'
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)
1560
1561       om1yk=0.d0
1562
1563       imin=ntymin
1564       imax=ntymx
1565       if(iomega.eq.2)imax=1
1566
1567       do i=imin,imax
1568           gamomy(i)=atilde(i,k)
1569           deltap=btildep(i,k)
1570           deltam=btildepp(i,k)
1571
1572         alpp(i)=deltap-deltam
1573         gamomy(i)=gamomy(i)*xh**((deltap+deltam)*0.5d0)
1574
1575       enddo
1576
1577       cint=0.d0
1578       do i=imin,imax
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)))
1582         else
1583           cint=cint-gamomy(i)*log(xh)
1584         endif
1585       enddo
1586
1587       do i=imin,imax
1588         if(abs(alpp(i)).gt.eps)then
1589         om1yk=om1yk+gamomy(i)*exp(alpp(i)*yp)
1590         else
1591         om1yk=om1yk+gamomy(i)
1592         endif
1593       enddo
1594
1595       om1yk=om1yk/cint
1596
1597
1598       return
1599       end
1600
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.
1606 c k - pair indice;
1607 c----------------------------------------------------------------------
1608       include 'epos.inc'
1609       include 'epos.incems'
1610       include 'epos.incpar'
1611
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)
1616
1617
1618       om1xrk=0.d0
1619
1620       imin=ntymin
1621       imax=ntymx
1622       if(iomega.eq.2)imax=1
1623
1624       do i=imin,imax
1625         gamomx(i)=atilde(i,k)
1626         deltap(i)=btildep(i,k)
1627         deltam(i)=btildepp(i,k)
1628
1629         delta(i)=(deltap(i)+deltam(i))*0.5d0
1630         alpp(i)=deltap(i)-deltam(i)
1631
1632       enddo
1633
1634       cint=0.d0
1635       do i=imin,imax
1636         gamom=gamomx(i)
1637         if((deltap(i)+1.d0).gt.eps)then
1638           gamom=gamom/(deltap(i)+1.d0)
1639         else
1640           gamom=-gamom*log(xminDf)
1641         endif
1642         if((deltam(i)+1.d0).gt.eps)then
1643           gamom=gamom/(deltam(i)+1.d0)
1644         else
1645           gamom=-gamom*log(xminDf)
1646         endif
1647         cint=cint+gamom
1648       enddo
1649
1650       x0=eps
1651       x1=1.d0
1652       f0=0.d0
1653       f1=0.d0
1654       do i=imin,imax
1655
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))
1662           else
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)
1665           endif
1666         else
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))
1683           endif
1684         endif
1685       enddo
1686       f0=-f0/cint
1687       f1=-f1/cint
1688       ntry=0
1689  11   ntry=ntry+1
1690       r=drangen(dble(ntry))
1691       f0t=f0+r
1692       f1t=f1+r
1693       if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1694       if(f1t*f0t.ge.eps)then
1695         do i=imin,imax
1696          write(ifmt,*)i,gamomx(i),deltap(i),deltam(i),alpp(i),delta(i)
1697         enddo
1698         write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k
1699         call utstop('om1xrk (1)&')
1700       endif
1701       f0=f0t
1702       f1=f1t
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
1708 c        call utmsgf
1709 c        om1xrk=x0
1710 c        return
1711 c      endif
1712       if(abs(f0).lt.eps) then
1713         om1xrk=x0
1714         return
1715       endif
1716       if(abs(f1).lt.eps) then
1717         om1xrk=x1
1718         return
1719       endif
1720 c      x=(x1+x0)*0.5D0
1721       x=sqrt(x1*x0)
1722       deltx=abs(x1-x0)
1723
1724       ntry=0
1725       fx=0.d0
1726       fpx=0.d0
1727       xt=x
1728  111  continue
1729
1730       if(ntry.le.1000)then
1731       fx=0.d0
1732       fpx=0.d0
1733       do i=imin,imax
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)
1739           else
1740         fx=fx-0.5d0*gamomx(i)*(log(x)**2-log(xminDf)**2)
1741         fpx=fpx-gamomx(i)*log(x)/x
1742           endif
1743         else
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))
1757           endif
1758         endif
1759       enddo
1760       fx=-fx/cint+r
1761       fpx=fpx/cint
1762       xt=x-fx/fpx
1763
1764       if (f0*fx.lt.-eps) then
1765         f1=fx
1766         x1=x
1767       else
1768         f0=fx
1769         x0=x
1770       endif
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)
1773       endif
1774
1775       else
1776
1777         write(ifmt,*)'Warning in om1xrk, to much try !'
1778
1779       endif
1780
1781       if(abs(x-xt).gt.deltx*0.5d0) then
1782         xt=sqrt(x1*x0)
1783       endif
1784       deltx=abs(x-xt)
1785       if(abs(x).gt.eps)then
1786         prec=deltx/x
1787       else
1788         call utstop('Problem in om1xrk&')
1789       endif
1790
1791       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000)then
1792          x=xt
1793          ntry=ntry+1
1794          goto 111
1795       endif
1796
1797       om1xrk=x
1798
1799       return
1800       end
1801
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;
1808 c k - pair indice;
1809 c----------------------------------------------------------------------
1810       include 'epos.inc'
1811       include 'epos.incems'
1812
1813       double precision xh,r!,y0,y1,y,gamomy(ntymi:ntymx),eps,ymin,prec,yt
1814
1815       r=dble(rangen())
1816
1817       om1yrk=(0.5d0-r)*log(xh)
1818       return
1819
1820       end
1821
1822 c----------------------------------------------------------------------
1823       function ffom12aii(iq,je1,je2)   !---test---
1824 c----------------------------------------------------------------------
1825       include 'epos.inc'
1826       ig=5
1827       xmin=0.01/engy
1828       xmax=1
1829       r2=0
1830       do i2=1,ig
1831       do m2=1,2
1832        xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5))
1833        r1=0
1834        do i1=1,ig
1835        do m1=1,2
1836         xp=xmin+(xmax-xmin)*(.5+tgss(ig,i1)*(m1-1.5))
1837         f=ffom12a(xp,xm,iq,iq,je1,je2)
1838         r1=r1+wgss(ig,i1)*f
1839        enddo
1840        enddo
1841        f=r1*0.5*(xmax-xmin)
1842        r2=r2+wgss(ig,i2)*f
1843       enddo
1844       enddo
1845       ffom12aii=r2*0.5*(xmax-xmin)
1846       end
1847
1848 c----------------------------------------------------------------------
1849       function ffom12ai(xp,iq1,iq2,je1,je2)   !---test---
1850 c----------------------------------------------------------------------
1851       include 'epos.inc'
1852       ig=5
1853       xmin=0.01/engy
1854       xmax=1
1855       r2=0
1856       do i2=1,ig
1857       do m2=1,2
1858        xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5))
1859        f=ffom12a(xp,xm,iq1,iq2,je1,je2)
1860        r2=r2+wgss(ig,i2)*f
1861       enddo
1862       enddo
1863       ffom12ai=r2*0.5*(xmax-xmin)
1864       end
1865
1866 c----------------------------------------------------------------------
1867       function ffom12a(xp,xm,iq1,iq2,je1,je2)   !---test---
1868 c----------------------------------------------------------------------
1869 c
1870 c      2*om52*F*F == PomInc
1871 c
1872 c  xp - xplus
1873 c  xm - xminus
1874 c                              iq=1 .... sea-sea
1875 c  iq1 - min iq                iq=2 .... val-sea
1876 c  iq2 - max iq                iq=3 .... sea-val
1877 c                              iq=4 .... val-val
1878 c  je = emission type (projectile and target side)
1879 c          0 ... no emissions
1880 c          1 ... emissions
1881 c       else ... all
1882 c
1883 c  already b-averaged  (\int d2b /sigine*10)
1884 c----------------------------------------------------------------------
1885       include 'epos.inc'
1886
1887       sy=engy*engy
1888       xh=xm*xp
1889 ctp060829      yp=0.5*log(xp/xm)
1890       ffom12a=0
1891       do i=iq1,iq2
1892        if(i.eq.1)then
1893          ffom12a=ffom12a+2*om52pi(sy*xh,1.,1.,0,je1,je2)
1894        elseif(i.eq.2)then
1895          ffom12a=ffom12a+2*om52pi(sy*xh,xp,1.,1,je1,je2)
1896        elseif(i.eq.3)then
1897          ffom12a=ffom12a+2*om52pi(sy*xh,xm,1.,2,je1,je2)
1898        elseif(i.eq.4)then
1899                ffom12a=ffom12a+2*om52pi(sy*xh,xp,xm,3,je1,je2)
1900        endif
1901       enddo
1902       ffom12a=ffom12a
1903      *           *alpff(iclpro)*xp**betff(1)*(1-xp)**alplea(iclpro)
1904      *           *alpff(icltar)*xm**betff(2)*(1-xm)**alplea(icltar)
1905
1906       end
1907
1908 c----------------------------------------------------------------------
1909       function ffom11a(xp,xm,iq1,iq2)   !---test---
1910 c----------------------------------------------------------------------
1911 c
1912 c      int(db) om1ff /sigine*10
1913 c
1914 c  xp - xplus                  iq=-1 ... fit
1915 c  xm - xminus                 iq=0 .... soft
1916 c                              iq=1 .... gg
1917 c  iq1 - min iq                iq=2 .... qg
1918 c  iq2 - max iq                iq=3 .... gq
1919 c                              iq=4 .... qq
1920 c                              iq=5 .... diff
1921 c----------------------------------------------------------------------
1922       include 'epos.inc'
1923       common/geom/rmproj,rmtarg,bmax,bkmx
1924       ig=5
1925       bmid=bkmx/2.
1926       r=0.d0
1927       do i=1,ig
1928         do m=1,2
1929           bb=bmid*(1.+(2.*m-3)*tgss(ig,i))
1930           f=ffom11(xp,xm,bb,iq1,iq2)
1931           r=r+bb*wgss(ig,i)*f
1932         enddo
1933       enddo
1934       ffom11a=r*2.*pi*bmid  /sigine*10
1935       return
1936       end
1937
1938 c----------------------------------------------------------------------
1939       function ffom11(xp,xm,b,iq1,iq2)   !---test---
1940 c----------------------------------------------------------------------
1941 c
1942 c       2*om5*F*F == PomInc
1943 c
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
1949 c                              iq=4 .... qq
1950 c                              iq=5 .... diff
1951 c----------------------------------------------------------------------
1952       include 'epos.inc'
1953       double precision om51
1954
1955       if(xm.ge.0.)then
1956
1957        xh=xm*xp
1958        yp=0.5*log(xp/xm)
1959        ffom11=2*om51(dble(xh),dble(yp),b,iq1,iq2)
1960      *     *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
1961
1962       else   !xm integration
1963
1964        ig=5
1965        xmin=0.01/engy
1966        xmax=1
1967        r=0
1968        do i=1,ig
1969        do m=1,2
1970         xmm=xmin*(xmax/xmin)**(.5+tgss(ig,i)*(m-1.5))
1971         xh=xmm*xp
1972         yp=0.5*log(xp/xmm)
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
1976        enddo
1977        enddo
1978        ffom11=r*0.5*log(xmax/xmin)
1979
1980       endif
1981
1982       end
1983
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
1992 c                       iq=4 .... qq
1993 c                       iq=5 .... diff
1994 c----------------------------------------------------------------------
1995       include 'epos.inc'
1996       include 'epos.incsem'
1997       include 'epos.incpar'
1998       double precision xp,xm,xh,yp,om51p,om1
1999       om51=0.d0
2000       if(xh.le.0.d0.or.xh.gt.1.d0)return
2001
2002       sy=engy*engy
2003       xp=sqrt(xh)*exp(yp)
2004       xm=xh/xp
2005
2006       if(iq1.eq.-1.and.iq2.eq.-1)then
2007         om51=0.5d0*om1(xh,yp,b)
2008       elseif(iq1.ge.0)then
2009         om51=0.d0
2010         do i=iq1,iq2
2011           if(i.ne.5)then
2012             i1=min(i,1)
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)
2017           else
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
2020           endif
2021        enddo
2022       else
2023         stop'om5: choice of iq1 and iq2 is nonsense.     '
2024       endif
2025
2026       end
2027
2028 c----------------------------------------------------------------------
2029       double precision function om5s(sx,xh,yp,b,iq1,iq2)   !---test---
2030 c----------------------------------------------------------------------
2031       include 'epos.inc'
2032       double precision om51
2033       double precision xh,yp
2034       ss=sx/xh
2035       engyx=engy
2036       engy=sqrt(ss)
2037       om5s=om51(xh,yp,b,iq1,iq2)
2038       engy=engyx
2039       end
2040
2041 c----------------------------------------------------------------------
2042       double precision function om5Jk(k,xh,yp,iqq)   !---MC---
2043 c----------------------------------------------------------------------
2044 c partial om5
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;
2048 c iqq=0 - soft
2049 c iqq=1 - gg
2050 c iqq=2 - qg
2051 c iqq=3 - gq
2052 c iqq=4 - qq
2053 c iqq=5 - diffractif
2054 c----------------------------------------------------------------------
2055       include 'epos.inc'
2056       include 'epos.incems'
2057       include 'epos.incsem'
2058
2059       double precision xh,yp,om51p
2060       double precision plc,s
2061       common/cems5/plc,s
2062
2063       sy=real(s*xh)
2064       b=bk(k)
2065
2066       if(iqq.ne.5)then
2067         om5Jk=om51p(sy,xh,yp,b,iqq)
2068         if(iscreen.ne.0)then
2069           xp=sqrt(xh)*exp(yp)
2070           xm=xh/xp
2071           if(iqq.eq.0)then
2072             epsG=epsilongs(k,0)
2073             epsGp=epsilongp(k,0)
2074             epsGt=epsilongt(k,0)
2075           else
2076             epsG=epsilongs(k,1)
2077             epsGp=epsilongp(k,1)
2078             epsGt=epsilongt(k,1)
2079           endif
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)
2084         endif
2085       else
2086         xp=sqrt(xh)*exp(yp)
2087         xm=xh/xp
2088         om5Jk=0.5d0*atilde(2,k)*xp**btildep(2,k)*xm**btildepp(2,k)
2089       endif
2090
2091       return
2092       end
2093
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
2101 c iqq=1 soft
2102 c iqq=2 gg
2103 c----------------------------------------------------------------------
2104       include 'epos.inc'
2105       include 'epos.incems'
2106       include 'epos.incsem'
2107       include 'epos.incpar'
2108
2109       double precision Df
2110
2111       Df=0.d0
2112       sy=engy*engy
2113       omIgamint=0.d0
2114       imax=idxD1
2115       if(iomega.eq.2)imax=1
2116
2117       if(iqq.eq.0)then
2118         coefp=1.+alplea(iclpro)
2119         coeft=1.+alplea(icltar)
2120
2121         do i=idxDmin,imax
2122           call GfunPar(1,i,b,sy,alpx,betx,betpx,epsp,epst,epss,gamv)
2123           betp=1.+betx
2124           betpp=1.+betpx
2125           Df=alpx*dble(utgam1(betp)*utgam1(betpp)*ucfpro
2126      *         *ucftar/utgam1(betp+coefp)/utgam1(betpp+coeft))
2127           omIgamint=omIgamint+Df
2128         enddo
2129       else
2130         call utstop('Wrong iqq in omIgamint&')
2131       endif
2132
2133       omIgamint=omIgamint
2134      *     *dble(chad(iclpro)*chad(icltar))
2135
2136       omIgamint=omIgamint*0.5d0
2137
2138       return
2139       end
2140
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.
2146 c k - pair indice;
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
2150 c    ity = 0   - soft
2151 c    ity = 1   - gg
2152 c    ity = 2   - qg
2153 c    ity = 3   - gq
2154 c    ity = 4   - qq
2155 c-----------------------------------------------------------------------
2156       include 'epos.inc'
2157       include 'epos.incems'
2158       include 'epos.incsem'
2159       doubleprecision xh,yp,om5Jk,w(0:7)
2160
2161       do i=0,7
2162         w(i)=0.d0
2163       enddo
2164
2165       do i=0,5
2166         w(i)=om5Jk(k,xh,yp,i)
2167       enddo
2168
2169       if(iotst1.gt.0)then !????????????????????????????????????????????????
2170         corfac=float(iotst1)
2171         ww=w(0)+w(1)+w(2)+w(3)+w(4)+w(5)
2172         whard=0.
2173         do i=1,4
2174           w(i)=w(i)*corfac
2175           whard=whard+w(i)
2176         enddo
2177         if(whard.gt.ww)then
2178          do i=1,4
2179           w(i)=w(i)/whard*ww
2180          enddo
2181          whard=ww
2182         endif
2183         w05=w(0)+w(5)
2184         if(whard.lt.ww)then
2185           w(0)=w(0)/w05*(ww-whard)
2186           w(5)=w(5)/w05*(ww-whard)
2187         else
2188           w(0)=0
2189           w(5)=0
2190         endif
2191         !write(*,'(2f11.4)')(ww-w05)/ww,(ww-w(0)-w(5))/ww
2192       endif !?????????????????????????????????????????????????????????????????
2193
2194       return
2195       end
2196
2197
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;
2204 c k - pair index
2205 c-----------------------------------------------------------------------
2206       include 'epos.inc'
2207       include 'epos.incems'
2208       double precision xp,xm,xprem,xmrem
2209
2210       Womegak=0.d0
2211
2212       imax=ntymx
2213       if(iomega.eq.2)imax=1
2214
2215       do i=ntymin,imax
2216         Womegak=Womegak+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k)
2217       enddo
2218
2219
2220       Womegak=Womegak*(xprem-xp)*(xmrem-xm)
2221
2222       return
2223       end
2224
2225
2226 cc----------------------------------------------------------------------
2227 c      double precision function omNpcut(xp,xm,xprem,xmrem,bh,iqq)   !---test---
2228 cc----------------------------------------------------------------------
2229 cc Sum of all cut diagrams
2230 cc iqq=0 ideal G
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
2235 c
2236 c      omNpcut=0.d0
2237 c      xh=xp*xm
2238 c      if(abs(xh).gt.1.d-10)then
2239 c        yp=0.5d0*log(xp/xm)
2240 c      else
2241 c        yp=0.d0
2242 c      endif
2243 c
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)
2246 c
2247 c      omNpcut=omNpcut*2.d0
2248 c
2249 c      return
2250 c      end
2251 c
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-----------------------------------------------------------------------
2259       include "epos.inc"
2260       include "epos.incems"
2261       double precision om51,xp,xm,xh,yp,eps!,omYgam
2262       parameter (eps=1.d-20)
2263
2264       omGam=0.d0
2265       if(xp.lt.eps.or.xm.lt.eps)return
2266       xh=xp*xm
2267       if(abs(xh).gt.1.d-10)then
2268         yp=0.5d0*log(xp/xm)
2269       else
2270         yp=0.d0
2271       endif
2272
2273       omGam=om51(xh,yp,bh,-1,-1)
2274
2275       omGam=2.d0*omGam
2276
2277       return
2278       end
2279
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-----------------------------------------------------------------------
2287       include "epos.inc"
2288       include "epos.incems"
2289       double precision xp,xm
2290       omGamk=0.d0
2291       imax=ntymx
2292       if(iomega.eq.2)imax=1
2293       do i=ntymin,imax
2294         omGamk=omGamk+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k)
2295       enddo
2296
2297       return
2298       end
2299
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-----------------------------------------------------------------------
2306       include "epos.inc"
2307       double precision omIgamint!,omYgamint
2308
2309       omGamint=2.d0*omIgamint(bh,0)
2310
2311       return
2312       end
2313
2314
2315
2316
2317
2318 c----------------------------------------------------------------------
2319       block data dgdata
2320 c----------------------------------------------------------------------
2321 c constants for numerical integration (gaussian weights)
2322 c----------------------------------------------------------------------
2323       double precision dgx1,dga1
2324       common /dga20/ dgx1(10),dga1(10)
2325
2326
2327       data dgx1/
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/
2338       data dga1/
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/
2349
2350       end
2351
2352
2353
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----------------------------------------------------------------------
2359       include 'epos.inc'
2360       include 'epos.incsem'
2361       include 'epos.incems'
2362       double precision al(idxD0:idxD1),betp(idxD0:idxD1)
2363      *,z,xIrst!,ffacto
2364       double precision zp(idxD0:idxD1),Phitmp,betpp(idxD0:idxD1)
2365      *,yp,ym,xm,xp
2366       double precision eps
2367       parameter(eps=1.d-20)
2368       dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
2369
2370       if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in PhiExact"
2371       PhiExact=0.d0
2372
2373       if(xp.gt.eps.and.xm.gt.eps.and.xp.le.1.d0+eps
2374      &   .and.xm.le.1.d0+eps)then
2375
2376
2377        do i=idxD0,idxD1
2378         imax(i)=0
2379         ipr(i)=0
2380         zp(i)=1.d0
2381         al(i)=0.d0
2382         betp(i)=0.d0
2383         betpp(i)=0.d0
2384        enddo
2385
2386        imax0=idxD1
2387        if(iomega.eq.2)imax0=1
2388
2389        do i=idxDmin,imax0
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))
2393        enddo
2394        Phitmp=0.d0
2395        do i=idxDmin,imax0
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))
2401        enddo
2402
2403        do ipr0=0,imax(0)
2404           ipr(0)=ipr0
2405           zp(0)=1.d0
2406         if (ipr(0).ne.0) zp(0)=(-dble(fj)*al(0))**ipr(0)*facto(ipr(0))
2407         do ipr1=0,imax(1)
2408            ipr(1)=ipr1
2409            zp(1)=1.d0
2410         if (ipr(1).ne.0) zp(1)=(-dble(fj)*al(1))**ipr(1)*facto(ipr(1))
2411         do ipr2=0,imax(2)
2412            ipr(2)=ipr2
2413            zp(2)=1.d0
2414         if (ipr(2).ne.0) zp(2)=(-dble(fj)*al(2))**ipr(2)*facto(ipr(2))
2415           yp=0.d0
2416           ym=0.d0
2417           z=1.d0
2418           isum=0
2419           do i=idxDmin,imax0
2420             yp=yp+dble(ipr(i))*betp(i)
2421             ym=ym+dble(ipr(i))*betpp(i)
2422             isum=isum+ipr(i)
2423             z=z*zp(i)
2424           enddo
2425
2426           z=z*xIrst(1,xp,yp,betp,ipr)
2427           z=z*xIrst(2,xm,ym,betpp,ipr)
2428
2429           Phitmp=Phitmp+z
2430
2431          enddo
2432          enddo
2433         enddo
2434
2435
2436
2437       endif
2438
2439       PhiExact=Phitmp
2440
2441
2442       return
2443       end
2444
2445
2446 c----------------------------------------------------------------------
2447       double precision function PhiExpoK(k,xp,xm)   !---MC---
2448 c----------------------------------------------------------------------
2449 c    Exponential expression of the Phi function for pp collision
2450 c    for given k
2451 c----------------------------------------------------------------------
2452       include 'epos.inc'
2453       include 'epos.incsem'
2454       include 'epos.incems'
2455
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)
2460
2461
2462       Phitmp=0.d0
2463
2464       imax=idxD1
2465       if(iomega.eq.2)imax=1
2466
2467       Phitmp=0.d0
2468       Gt1=0.d0
2469       do i=idxDmin,imax
2470        Gt1=Gt1+atildg(i,k)*xp**btildgp(i,k)*xm**btildgpp(i,k)
2471       enddo
2472
2473       Phitmp=exp(-Gt1)
2474
2475       PhiExpoK=Phitmp
2476
2477       return
2478       end
2479
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
2484 c    for given b
2485 c----------------------------------------------------------------------
2486       include 'epos.inc'
2487       include 'epos.incsem'
2488       include 'epos.incems'
2489       include 'epos.incpar'
2490
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
2496
2497       imax=idxD1
2498       if(iomega.eq.2)imax=1
2499
2500       Gt1=0.d0
2501       do i=idxDmin,imax
2502         call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
2503         BeTip =dble(betx)
2504         BeTipp=dble(betpx)
2505         AlTi  =dble(alpx)
2506         Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp*dble(fj*xkappa**(2.5*fj-2.5))
2507       enddo
2508
2509       Phitmp=exp(-Gt1)
2510
2511       PhiExpo=Phitmp
2512      &     *xp**dble(alplea(iclpro))
2513      &     *xm**dble(alplea(icltar))
2514
2515       return
2516       end
2517
2518 c----------------------------------------------------------------------
2519       double precision function PhiUnit(xp,xm)   !---test---
2520 c----------------------------------------------------------------------
2521 c    Exponential expression of the Phi function for pp collision
2522 c    for given b
2523 c----------------------------------------------------------------------
2524       include 'epos.inc'
2525       include 'epos.incsem'
2526       include 'epos.incems'
2527       include 'epos.incpar'
2528
2529       double precision AlTi
2530       double precision BeTip,BeTipp
2531       double precision xp,xm,Phitmp,Gt1
2532
2533       imax=idxD1
2534       if(iomega.eq.2)imax=1
2535
2536       Gt1=0.d0
2537       do i=idxDmin,imax
2538         BeTip =betUni(i,2)
2539         BeTipp=betpUni(i,2)
2540         AlTi  =alpUni(i,2)
2541         Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp
2542 c        write(ifch,*)'Phiunit',i,xp,xm,Gt1,AlTi,BeTip,BeTipp
2543       enddo
2544
2545       Phitmp=exp(-Gt1)
2546
2547       PhiUnit=Phitmp
2548      &     *xp**dble(alplea(iclpro))
2549      &     *xm**dble(alplea(icltar))
2550
2551
2552       return
2553       end
2554
2555
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
2561 c
2562 c      PhiUnit=PhiExpo(1.,xp,xm,s,b)
2563 c     &          /Znorm(s,b)
2564 c
2565 c      return
2566 c      end
2567 c
2568
2569 c----------------------------------------------------------------------
2570       double precision function Hrst(s,b,xp,xm)   !test
2571 c----------------------------------------------------------------------
2572       include 'epos.inc'
2573       include 'epos.incems'
2574       include 'epos.incsem'
2575       include 'epos.incpar'
2576       parameter(idxD2=8)
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)
2582      *,z,xJrst!,ffacto
2583       double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2)
2584      *,yp,ym,xp,xm
2585       dimension ipr(idxD0:idxD2),imax(idxD0:idxD2)
2586
2587       if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in Hrst"
2588       Hrst=0.d0
2589       do i=idxD0,idxD2
2590         imax(i)=0
2591         ipr(i)=0
2592         zp(i)=1.d0
2593         al(i)=0.d0
2594       enddo
2595
2596       if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then
2597
2598       imax0=idxD1
2599       if(iomega.eq.2)imax0=1
2600       imax1=idxD2
2601       if(iomega.eq.2)imax1=imax1-1
2602
2603       do i=idxDmin,imax1
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))
2608         if(i.gt.imax0)then
2609           if((zpUni*ztUni.lt.1.d-6)
2610      &   .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
2611             imax(i)=0
2612           else
2613             imax(i)=1     !imax(i)/3
2614           endif
2615         endif
2616       enddo
2617
2618       Htmp=0.d0
2619         do i=idxDmin,imax1
2620           betp(i)=HbetUni(i)
2621           betpp(i)=HbetpUni(i)
2622           al(i)=HalpUni(i)
2623         enddo
2624
2625         do ipr0=0,imax(0)
2626 c          write(ifmt,*)'Hrst ipr0,xp,xm :',ipr0,xp,xm
2627            ipr(0)=ipr0
2628            zp(0)=1.d0
2629            if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2630          do ipr1=0,imax(1)
2631             ipr(1)=ipr1
2632             zp(1)=1.d0
2633             if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2634          do ipr2=0,imax(2)
2635             ipr(2)=ipr2
2636             zp(2)=1.d0
2637             if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
2638          do ipr3=0,imax(3)
2639             ipr(3)=ipr3
2640             zp(3)=1.d0
2641             if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3))
2642          do ipr4=0,imax(4)
2643             ipr(4)=ipr4
2644             zp(4)=1.d0
2645             if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4))
2646          do ipr5=0,imax(5)
2647             ipr(5)=ipr5
2648             zp(5)=1.d0
2649             if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5))
2650          do ipr6=0,imax(6)
2651             ipr(6)=ipr6
2652             zp(6)=1.d0
2653             if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6))
2654          do ipr7=0,imax(7)
2655             ipr(7)=ipr7
2656             zp(7)=1.d0
2657             if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7))
2658          do ipr8=0,imax(8)
2659             ipr(8)=ipr8
2660             zp(8)=1.d0
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
2664              yp=0.d0
2665              ym=0.d0
2666              z=1.d0
2667              do i=idxDmin,imax1
2668                yp=yp+dble(ipr(i))*betp(i)
2669                ym=ym+dble(ipr(i))*betpp(i)
2670                z=z*zp(i)
2671              enddo
2672              z=z*xJrst(xp,yp,GbetUni,ipr)
2673              z=z*xJrst(xm,ym,GbetpUni,ipr)
2674              Htmp=Htmp+z
2675            endif
2676           enddo
2677          enddo
2678        enddo
2679           enddo
2680          enddo
2681        enddo
2682           enddo
2683          enddo
2684        enddo
2685
2686       endif
2687
2688       Hrst=Htmp
2689
2690       return
2691       end
2692
2693 c----------------------------------------------------------------------
2694       double precision function HrstI(s,b,xp,xm)   !test
2695 c----------------------------------------------------------------------
2696       include 'epos.inc'
2697       include 'epos.incems'
2698       include 'epos.incsem'
2699       include 'epos.incpar'
2700       parameter(idxD2=8)
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)
2706      *,z,xJrstI!,ffacto
2707       double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2)
2708      *,yp,ym,xp,xm
2709       dimension ipr(idxD0:idxD2),imax(idxD0:idxD2)
2710
2711       if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
2712       HrstI=0.d0
2713       do i=idxD0,idxD2
2714         imax(i)=0
2715         ipr(i)=0
2716         zp(i)=1.d0
2717         al(i)=0.d0
2718       enddo
2719
2720
2721       if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then
2722
2723
2724       imax0=idxD1
2725       if(iomega.eq.2)imax0=1
2726       imax1=idxD2
2727       if(iomega.eq.2)imax1=imax1-1
2728
2729       do i=idxDmin,imax1
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))
2734         if(i.gt.imax0)then
2735           if((zpUni*ztUni.lt.1.d-6)
2736      &   .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
2737             imax(i)=0
2738           else
2739             imax(i)=1   !imax(i)/3
2740           endif
2741         endif
2742       enddo
2743
2744       Htmp=0.d0
2745         do i=idxDmin,imax1
2746           betp(i)=HbetUni(i)
2747           betpp(i)=HbetpUni(i)
2748           al(i)=HalpUni(i)
2749         enddo
2750         do ipr0=0,imax(0)
2751            ipr(0)=ipr0
2752            zp(0)=1.d0
2753            if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2754          do ipr1=0,imax(1)
2755             ipr(1)=ipr1
2756             zp(1)=1.d0
2757             if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2758          do ipr2=0,imax(2)
2759             ipr(2)=ipr2
2760             zp(2)=1.d0
2761             if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
2762          do ipr3=0,imax(3)
2763             ipr(3)=ipr3
2764             zp(3)=1.d0
2765             if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3))
2766          do ipr4=0,imax(4)
2767             ipr(4)=ipr4
2768             zp(4)=1.d0
2769             if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4))
2770          do ipr5=0,imax(5)
2771             ipr(5)=ipr5
2772             zp(5)=1.d0
2773             if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5))
2774          do ipr6=0,imax(6)
2775             ipr(6)=ipr6
2776             zp(6)=1.d0
2777             if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6))
2778          do ipr7=0,imax(7)
2779             ipr(7)=ipr7
2780             zp(7)=1.d0
2781             if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7))
2782          do ipr8=0,imax(8)
2783             ipr(8)=ipr8
2784             zp(8)=1.d0
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
2788              yp=0.d0
2789              ym=0.d0
2790              z=1.d0
2791              do i=idxDmin,imax1
2792                yp=yp+dble(ipr(i))*betp(i)
2793                ym=ym+dble(ipr(i))*betpp(i)
2794                z=z*zp(i)
2795              enddo
2796              z=z*xJrstI(xp,yp,GbetUni,ipr)
2797              z=z*xJrstI(xm,ym,GbetpUni,ipr)
2798              Htmp=Htmp+z
2799            endif
2800           enddo
2801          enddo
2802        enddo
2803           enddo
2804          enddo
2805        enddo
2806           enddo
2807          enddo
2808        enddo
2809
2810       endif
2811
2812       HrstI=Htmp
2813
2814       return
2815       end
2816
2817
2818
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)
2829 c     *,yp,ym,xp,xm
2830 c      dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
2831 c
2832 c      if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
2833 c      HrstI=0.d0
2834 c      do i=idxD0,idxD1
2835 c        imax(i)=0
2836 c        ipr(i)=0
2837 c        zp(i)=1.d0
2838 c        al(i)=0.d0
2839 c      enddo
2840 c
2841 c      if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.lt.1.d0)then
2842 c
2843 c      HrstI=0.d0
2844 c
2845 c      imax0=idxD1
2846 c      if(iomega.eq.2)imax0=1
2847 c
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))
2852 c      enddo
2853 c      Htmp=0.d0
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))
2858 c        enddo
2859 c
2860 c        do ipr0=0,imax(0)
2861 c           ipr(0)=ipr0
2862 c           zp(0)=1.d0
2863 c          if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2864 c         do ipr1=0,imax(1)
2865 c            ipr(1)=ipr1
2866 c            zp(1)=1.d0
2867 c          if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2868 c         do ipr2=0,imax(2)
2869 c            ipr(2)=ipr2
2870 c            zp(2)=1.d0
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
2873 c             yp=0.d0
2874 c             ym=0.d0
2875 c             z=1.d0
2876 c             do i=idxDmin,imax0
2877 c               yp=yp+dble(ipr(i))*betp(i)
2878 c               ym=ym+dble(ipr(i))*betpp(i)
2879 c               z=z*zp(i)
2880 c             enddo
2881 c             z=z*xJrstI(xp,yp,betp,ipr)
2882 c             z=z*xJrstI(xm,ym,betpp,ipr)
2883 c             Htmp=Htmp+z
2884 c           endif
2885 c          enddo
2886 c         enddo
2887 c       enddo
2888 c
2889 c       HrstI=Htmp
2890 c
2891 c      endif
2892 c
2893 c      return
2894 c      end
2895 c
2896
2897 cc----------------------------------------------------------------------
2898 c        double precision function ffacto(n)   !---test---
2899 cc----------------------------------------------------------------------
2900 c
2901 c        ffacto=1.D0
2902 c        do i=1,n
2903 c          ffacto=ffacto*dble(i)
2904 c        enddo
2905 c        return
2906 c        end
2907 c
2908
2909 c----------------------------------------------------------------------
2910       double precision function xIrst(id,x,y,bet,ipr)   !---test---
2911 c----------------------------------------------------------------------
2912       include 'epos.inc'
2913       double precision y,gammag,utgam2,x,bet(idxD0:idxD1)
2914       dimension ipr(idxD0:idxD1)
2915
2916       if(id.eq.1)iclrem=iclpro
2917       if(id.eq.2)iclrem=icltar
2918       imax=idxD1
2919       if(iomega.eq.2)imax=1
2920       if(y.le.160.)then
2921        xIrst=gammag(iclrem,y)*x**dble(alplea(iclrem))
2922       else
2923        xIrst=0
2924       endif
2925       if(xIrst.gt.0.d0)then
2926         do i=idxDmin,imax
2927           if(ipr(i).ne.0.and.bet(i).gt.1.d-10)
2928      &         xIrst=xIrst*utgam2(bet(i))**dble(ipr(i))
2929         enddo
2930         if (abs(y).gt.1.d-10) xIrst=xIrst*x**y
2931       endif
2932       return
2933       end
2934
2935
2936 c----------------------------------------------------------------------
2937       double precision function xJrst(x,y,Gbeta,ipr)   !---inu---
2938 c----------------------------------------------------------------------
2939       include 'epos.inc'
2940       parameter(idxD2=8)
2941       double precision y,utgam2,x,Gbeta(idxD0:idxD2),eps,gam
2942       dimension ipr(idxD0:idxD2)
2943
2944       eps=1.d-10
2945
2946       imax=idxD2
2947       if(iomega.eq.2)imax=imax-1
2948
2949       gam=utgam2(y)
2950
2951       if(gam.lt.1.d99)then
2952
2953       if ((x-1.d0).gt.eps.or.(y-1.d0).gt.eps) then
2954                         xJrst=(1.d0-x)**(y-1.d0)/gam
2955                         do i=idxDmin,imax
2956       if (ipr(i).ne.0)   xJrst=xJrst*Gbeta(i)**dble(ipr(i))
2957                         enddo
2958           else
2959 c            write (*,*) 'Warning in xJrst, infinite value !'
2960                         xJrst=(1.d0-x+eps)**(y-1.d0)/gam
2961                         do i=idxDmin,imax
2962       if (ipr(i).ne.0)   xJrst=xJrst*Gbeta(i)**dble(ipr(i))
2963                         enddo
2964       endif
2965       else
2966         xJrst=0.d0
2967       endif
2968
2969       return
2970       end
2971
2972
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----------------------------------------------------------------------
2979       include 'epos.inc'
2980       parameter(idxD2=8)
2981       double precision y,utgam2,x,Gbeta(idxD0:idxD2),alpha,w,gam
2982       dimension ipr(idxD0:idxD2)
2983
2984       alpha=4.d0
2985       w=alpha*(y-1.d0)+alpha-1.d0
2986       imax=idxD2
2987       if(iomega.eq.2)imax=imax-1
2988
2989       gam=utgam2(y)
2990
2991       if(gam.lt.1.d99)then
2992
2993       if (w.ge.0)then
2994
2995                         xJrstI=alpha*x**w/gam
2996                         do i=idxDmin,imax
2997       if (ipr(i).ne.0)   xJrstI=xJrstI*Gbeta(i)**dble(ipr(i))
2998                         enddo
2999
3000          else
3001            write(*,*) 'x,y,bet,ipr,w',x,y,Gbeta,ipr,w
3002           stop 'Error in xJrstI in epos-omg, integration not possible'
3003        endif
3004
3005       else
3006         xJrstI=0.d0
3007       endif
3008
3009
3010       return
3011       end
3012
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----------------------------------------------------------------------
3019       include 'epos.inc'
3020       parameter(idxD2=8)
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)
3029
3030       eps=0d0 !1.d-5
3031
3032       imax0=idxD1
3033       imax1=idxD2
3034       if(iomega.eq.2)then
3035         imax0=1
3036         imax1=imax1-1
3037       endif
3038       do i=idxDmin,imax0
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))
3044       enddo
3045       do i=0,1
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)
3061       enddo
3062
3063       w=0.d0
3064       xhm=.5d0*(1d0-eps)
3065       yhm=.5d0*(1d0-eps)
3066       do m=1,2
3067         do i=1,3
3068 c        do i=1,7
3069           x=xhm+dble((2*m-3)*x9(i))*xhm
3070 c          write(ifmt,*)'HPhiInt, xp int :',x
3071           do n=1,2
3072             do j=1,3
3073 c            do j=1,7
3074               y=yhm+dble((2*n-3)*x9(j))*yhm
3075              w=w+dble(a9(i)*a9(j))*Hrst(s,b,x,y)
3076      &          *PhiUnit(x,y)
3077 c     &          *PhiExact(1.,x,y,s,b)
3078             enddo
3079           enddo
3080         enddo
3081       enddo
3082
3083       HPhiInt=w*xhm*yhm
3084
3085
3086 c      w=0.d0
3087 c      xhm=.5d0*eps
3088 c      yhm=.5d0*eps
3089 c      do m=1,2
3090 c        do i=1,7
3091 c          x=1d0-eps+xhm+dble((2*m-3)*x1(i))*xhm
3092 c          do n=1,2
3093 c            do j=1,7
3094 c              y=1d0-epsyhm+dble((2*n-3)*x1(j))*yhm
3095 c              zp2=1.d0-x**4
3096 c              zm2=1.d0-y**4
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)
3100 c            enddo
3101 c          enddo
3102 c        enddo
3103 c      enddo
3104 c
3105 c      HPhiInt=HPhiInt+w*xhm*yhm
3106
3107       return
3108       end
3109
3110
3111
3112 c----------------------------------------------------------------------
3113       subroutine Kfit(iiclegy)
3114 c----------------------------------------------------------------------
3115       include "epos.inc"
3116       include "epos.incsem"
3117       double precision Znorm
3118       parameter(nbkbin=40)
3119       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
3120       parameter (nmax=30)
3121       logical lnoch
3122
3123       if(iiclegy.eq.-1.or.iiclegy.gt.iclegy2)then
3124         do iiiegy=1,nclegy
3125         do iiipro=1,nclha
3126         do iiitar=1,nclha
3127         do iiibk=1,nbkbin
3128           xkappafit(iiiegy,iiipro,iiitar,iiibk)=1.
3129         enddo
3130         enddo
3131         enddo
3132         enddo
3133
3134       else
3135
3136       if(isetcs.le.1)then
3137         s=engy*engy
3138         eps=0.05
3139       else
3140         s=(egylow*egyfac**(iiclegy-1))**2.
3141         eps=0.001
3142       endif
3143
3144       write(ifmt,*)"Fit xkappa ..."
3145       if(ish.ge.2)then
3146         write(ifmt,*)"Kfit s,bkbin,iclegy,ipro,itar"
3147      *       ,s,bkbin,iiclegy,iclpro,icltar
3148       endif
3149
3150
3151       b=0.
3152       if(isetcs.le.1.or.iiclegy.eq.iclegy2)then
3153         xkf=1.
3154       else
3155         xkf=xkappafit(iiclegy+1,iclpro,icltar,1)
3156       endif
3157       delta=0.
3158
3159
3160       do 5 ib=1,nbkbin-1
3161         b=real(ib-1)*bkbin
3162         xkappafit(iiclegy,iclpro,icltar,ib)=1.
3163         if(b.gt.3.+0.05*log(s).or.s.le.20.*q2min)then
3164           xkf=1.
3165           goto 5
3166         endif
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
3170           if(xkf.ne.1.)then
3171             xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3172             delta=1.-real(Znorm(s,b))
3173           endif
3174         else!if(xkf.ne.1.)then
3175           goto 5
3176         endif
3177         if(abs(delta).lt.eps)then
3178           if(delta.lt.0d0)then
3179             xkfs=xkf-delta
3180             deltas=delta
3181           endif
3182           xkf=1.
3183           goto 5
3184         elseif(ib.le.nbkbin-1)then
3185
3186           if(delta.gt.0.d0)then
3187             xkf0=1.
3188             xkf1=xkf
3189             delta0=delta
3190             xkf2=xkf-delta0
3191             xkappafit(iiclegy,iclpro,icltar,ib)=xkf2
3192             delta1=1.-real(Znorm(s,b))
3193             if(delta1.lt.0.d0)then
3194               xkf0=xkf2
3195               xkf1=xkf
3196               delta=delta1
3197               xkf=xkf0
3198             else
3199               xkf1=max(delta0,xkf2)
3200               xkf0=0.
3201               xkf=xkf1
3202             endif
3203           else
3204             xkf0=xkf
3205             xkf1=1.-delta
3206             xkf2=xkf
3207             delta1=delta
3208           endif
3209
3210           if(ib.eq.1)then
3211             deltas=delta
3212             xkfs=max(0.00001,1.-delta)
3213           endif
3214
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
3218
3219           n=0
3220           delta0=delta
3221           lnoch=.true.
3222  10       continue
3223           n=n+1
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))
3228             endif
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
3232             else
3233               lnoch=.false.
3234             endif
3235             if(abs(delta).gt.eps)then
3236               if(delta.gt.0.)then
3237                 xkf1=xkf
3238                 xkf=(xkf1+xkf0)*0.5
3239                 delta0=delta
3240               else
3241                 xkf0=xkf
3242                 xkf=(xkf1+xkf0)*0.5
3243                 delta0=delta
3244               endif
3245               goto 10
3246             endif
3247           else
3248             if(ish.ge.2)
3249      *      write(ifmt,*)"Warning in Kfit, nmax reached : xkappafit=1."
3250             xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3251           endif
3252         endif
3253
3254  5    continue
3255
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.
3260
3261       endif
3262
3263       return
3264       end
3265
3266 c----------------------------------------------------------------------
3267       double precision function Znorm(s,b)   !---inu---
3268 c----------------------------------------------------------------------
3269       include 'epos.inc'
3270       common /kwrite/ xkapZ
3271       double precision HPhiInt,PhiUnit!,PhiExact
3272
3273 c      write(ifmt,*)'Z calculation for (s,b) :',s,b
3274       imax=idxD1
3275       if(iomega.eq.2)imax=1
3276       do i=idxDmin,imax
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)
3279       enddo
3280       call GfomPar(b,s)
3281       Znorm=HPhiInt(s,b)
3282 c      write(ifch,*)'int',Znorm,' phi',PhiExact(1.,1.d0,1.d0,s,b)
3283       Znorm=Znorm
3284      &       +PhiUnit(1.d0,1.d0)
3285 c     &       +PhiExact(1.,1.d0,1.d0,s,b)
3286
3287       !write(ifmt,*)'Z=',Znorm,xkapZ,b
3288       return
3289       end
3290
3291
3292 c------------------------------------------------------------
3293       double precision function gammag(iclrem,x)   !---test---
3294 c----------------------------------------------------------------------
3295       include 'epos.inc'
3296       include 'epos.incsem'
3297       double precision x,utgam2
3298
3299       gammag=utgam2(dble(alplea(iclrem))+1.D0)
3300      &       /utgam2(dble(alplea(iclrem))+1.D0+x)
3301
3302       return
3303       end
3304
3305
3306
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)
3317 c
3318 c      bmid=bkmx/2.
3319 c      PomNbri=0.d0
3320 c      do i=1,7
3321 c        do m=1,2
3322 c          bb=bmid*(1.+(2.*m-3)*x1(i))
3323 c          PomNbri=PomNbri+dble(bb*a1(i))*om1intbci(bb,iqq)
3324 c        enddo
3325 c      enddo
3326 c
3327 c      PomNbri=PomNbri*dble(2.*pi*bmid)
3328 c
3329 c      return
3330 c      end
3331 c
3332 c
3333 c
3334
3335
3336 c####################################################################################
3337 c#############   former chk #########################################################
3338 c####################################################################################
3339
3340
3341
3342
3343
3344
3345
3346
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
3358 c
3359 cc Calculation by analytical integration (perfect but it changes
3360 cc if om1 change):
3361 c
3362 c      s=engy**2
3363 c      imax=1
3364 c      if(iomega.eq.2)imax=1
3365 c      do i=idxDmin,imax
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)
3370 c
3371 cc Integration possible only if delta(i)>-1
3372 c
3373 c       if(deltap(i).le.-1.d0.or.deltapp(i).le.-1.d0)
3374 c     &       stop 'Error in epos-par-300 in PomIncII'
3375 c      enddo
3376 c
3377 c      cint=0.d0
3378 c      do i=idxDmin,imax
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)
3383 c      enddo
3384 c
3385 c      PomIncII=cint
3386 c
3387 c      return
3388 c      end
3389 c
3390
3391 c----------------------------------------------------------------------
3392         double precision function PomIncXIExact(x)   !---check---
3393 c----------------------------------------------------------------------
3394 c integral d2b PomIncXExact
3395 c----------------------------------------------------------------------
3396       include 'epos.inc'
3397       double precision x,PomIncXExact
3398       common /ar3/  x1(7),a1(7)
3399       common/geom/rmproj,rmtarg,bmax,bkmx
3400
3401       bmid=bkmx/2.
3402       PomIncXIExact=0.d0
3403       do i=1,7
3404         do m=1,2
3405           bb=bmid*(1.+(2.*m-3)*x1(i))
3406           PomIncXIExact=PomIncXIExact+dble(bb*a1(i))*PomIncXExact(x,bb)
3407         enddo
3408       enddo
3409
3410       PomIncXIExact=PomIncXIExact*dble(2.*pi*bmid)
3411
3412       return
3413       end
3414
3415 c----------------------------------------------------------------------
3416         double precision function PomIncXIUnit(x)   !---check---
3417 c----------------------------------------------------------------------
3418 c integral d2b PomIncXUnit
3419 c----------------------------------------------------------------------
3420       include 'epos.inc'
3421       double precision x,PomIncXUnit
3422       common /ar3/  x1(7),a1(7)
3423       common/geom/rmproj,rmtarg,bmax,bkmx
3424
3425       bmid=bkmx/2.
3426       PomIncXIUnit=0.d0
3427       do i=1,7
3428         do m=1,2
3429           bb=bmid*(1.+(2.*m-3)*x1(i))
3430           PomIncXIUnit=PomIncXIUnit+dble(bb*a1(i))*PomIncXUnit(x,bb)
3431         enddo
3432       enddo
3433
3434       PomIncXIUnit=PomIncXIUnit*dble(2.*pi*bmid)
3435
3436       return
3437       end
3438
3439 c----------------------------------------------------------------------
3440         double precision function PomIncPIExact(x)   !---check---
3441 c----------------------------------------------------------------------
3442 c integral d2b PomIncPExact
3443 c----------------------------------------------------------------------
3444       include 'epos.inc'
3445       double precision x,PomIncPExact
3446       common/geom/rmproj,rmtarg,bmax,bkmx
3447       common /ar3/  x1(7),a1(7)
3448
3449       bmid=bkmx/2.
3450       PomIncPIExact=0.d0
3451       do i=1,7
3452         do m=1,2
3453           bb=bmid*(1.+(2.*m-3)*x1(i))
3454           PomIncPIExact=PomIncPIExact+dble(bb*a1(i))*PomIncPExact(x,bb)
3455         enddo
3456       enddo
3457
3458       PomIncPIExact=PomIncPIExact*dble(2.*pi*bmid)
3459
3460       return
3461       end
3462
3463 c----------------------------------------------------------------------
3464         double precision function PomIncPIUnit(x)   !---check---
3465 c----------------------------------------------------------------------
3466 c integral d2b PomIncPUnit
3467 c----------------------------------------------------------------------
3468       include 'epos.inc'
3469       double precision x,PomIncPUnit
3470       common/geom/rmproj,rmtarg,bmax,bkmx
3471       common /ar3/  x1(7),a1(7)
3472
3473       bmid=bkmx/2.
3474       PomIncPIUnit=0.d0
3475       do i=1,7
3476         do m=1,2
3477           bb=bmid*(1.+(2.*m-3)*x1(i))
3478           PomIncPIUnit=PomIncPIUnit+dble(bb*a1(i))*PomIncPUnit(x,bb)
3479         enddo
3480       enddo
3481
3482       PomIncPIUnit=PomIncPIUnit*dble(2.*pi*bmid)
3483
3484       return
3485       end
3486
3487 c----------------------------------------------------------------------
3488         double precision function PomIncMIExact(x)   !---check---
3489 c----------------------------------------------------------------------
3490 c integral d2b PomIncMExact
3491 c----------------------------------------------------------------------
3492       include 'epos.inc'
3493       double precision x,PomIncMExact
3494       common/geom/rmproj,rmtarg,bmax,bkmx
3495       common /ar3/  x1(7),a1(7)
3496
3497       bmid=bkmx/2.
3498       PomIncMIExact=0.d0
3499       do i=1,7
3500         do m=1,2
3501           bb=bmid*(1.+(2.*m-3)*x1(i))
3502           PomIncMIExact=PomIncMIExact+dble(bb*a1(i))*PomIncMExact(x,bb)
3503         enddo
3504       enddo
3505
3506       PomIncMIExact=PomIncMIExact*dble(2.*pi*bmid)
3507
3508       return
3509       end
3510
3511 c----------------------------------------------------------------------
3512         double precision function PomIncMIUnit(x)   !---check---
3513 c----------------------------------------------------------------------
3514 c integral d2b PomIncMUnit
3515 c----------------------------------------------------------------------
3516       include 'epos.inc'
3517       double precision x,PomIncMUnit
3518       common/geom/rmproj,rmtarg,bmax,bkmx
3519       common /ar3/  x1(7),a1(7)
3520
3521       bmid=bkmx/2.
3522       PomIncMIUnit=0.d0
3523       do i=1,7
3524         do m=1,2
3525           bb=bmid*(1.+(2.*m-3)*x1(i))
3526           PomIncMIUnit=PomIncMIUnit+dble(bb*a1(i))*PomIncMUnit(x,bb)
3527         enddo
3528       enddo
3529
3530       PomIncMIUnit=PomIncMIUnit*dble(2.*pi*bmid)
3531
3532       return
3533       end
3534
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----------------------------------------------------------------------
3540       include 'epos.inc'
3541       include 'epos.incsem'
3542       include 'epos.incems'
3543       double precision AlTiP,BeTip,al,bep,bepp,xpInt,utgam2,xm
3544
3545       s=engy**2
3546       PomIncMExact=0.d0
3547       imax=1
3548       if(iomega.eq.2)imax=1
3549       do i=idxDmin,imax
3550         call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3551         bep =dble(bet)
3552         bepp=dble(betp)
3553         al  =dble(alp*gamv)
3554
3555         BeTip=bep+1.d0
3556         xpInt=utgam2(BeTip)*dble(ucfpro)
3557      *                    /utgam2(1.d0+dble(alplea(iclpro))+BeTip)
3558         AlTiP=al*xpInt
3559         PomIncMExact=PomIncMExact+AlTiP*xm**bepp
3560      *                            *(1.d0-xm)**dble(alplea(icltar))
3561       enddo
3562
3563       return
3564       end
3565
3566 c----------------------------------------------------------------------
3567       double precision function PomIncMUnit(xm,b)   !---check---
3568 c----------------------------------------------------------------------
3569 c incluse  Unitarized Pomeron distribution  \int dx+
3570 c----------------------------------------------------------------------
3571       include 'epos.inc'
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)
3577
3578       s=engy**2
3579
3580 c Calculation by numeric integration :
3581       w=0.d0
3582       xpm=.5d0
3583       do m=1,2
3584         do j=1,7
3585           xp=xpm*(1.d0+dble((2.*m-3.)*x1(j)))
3586           xh=xp*xm
3587 ctp060829          sy=s*sngl(xh)
3588           Df=0.d0
3589           do i=idxDmin,idxD1
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)
3592           enddo
3593           call GfomPar(b,s)
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
3596         enddo
3597       enddo
3598       w=w*xpm
3599
3600
3601       PomIncMUnit=w!/Znorm(s,b)
3602
3603       return
3604       end
3605
3606
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----------------------------------------------------------------------
3612       include 'epos.inc'
3613       include 'epos.incsem'
3614       include 'epos.incems'
3615       double precision AlTiP,BeTipp,al,bep,bepp,xmInt,utgam2,xp
3616
3617       s=engy**2
3618       PomIncPExact=0.d0
3619       imax=1
3620       if(iomega.eq.2)imax=1
3621       do i=idxDmin,imax
3622         call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3623         bep=dble(bet)
3624         bepp=dble(betp)
3625         al=dble(alp*gamv)
3626         BeTipp=bepp+1.d0
3627         xmInt=utgam2(BeTipp)*dble(ucftar)
3628      *                    /utgam2(1.d0+dble(alplea(icltar))+BeTipp)
3629         AlTiP=al*xmInt
3630         PomIncPExact=PomIncPExact+AlTiP*xp**bep
3631      *                            *(1.d0-xp)**dble(alplea(iclpro))
3632       enddo
3633
3634       return
3635       end
3636
3637 c----------------------------------------------------------------------
3638       double precision function PomIncPUnit(xp,b)   !---check---
3639 c----------------------------------------------------------------------
3640 c incluse  Unitarized Pomeron distribution  \int dx-
3641 c----------------------------------------------------------------------
3642       include 'epos.inc'
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)
3647
3648       s=engy**2
3649       imax=1
3650       if(iomega.eq.2)imax=1
3651
3652 c Calculation by numeric integration :
3653       w=0.d0
3654       xmm=.5d0
3655       do m=1,2
3656         do j=1,7
3657           xm=xmm*(1.d0+dble((2.*m-3.)*x1(j)))
3658           xh=xp*xm
3659 ctp060829          sy=s*real(xh)
3660           Df=0.d0
3661           do i=idxDmin,imax
3662             call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3663             Df=Df+alp*real(xp)**bet*real(xm)**betp
3664           enddo
3665           call GfomPar(b,s)
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
3668         enddo
3669       enddo
3670       w=w*xmm
3671
3672
3673       PomIncPUnit=w!/Znorm(s,b)
3674
3675       return
3676       end
3677
3678
3679
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----------------------------------------------------------------------
3685       include 'epos.inc'
3686       double precision allea,PomIncXExact,xh
3687       common /ar3/  x1(7),a1(7)
3688
3689       allea=2.d0+dble(alplea(iclpro)+alplea(icltar))
3690       PomIncJExact=0.d0
3691       do i=1,7
3692         do m=1,2
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)
3696         enddo
3697       enddo
3698       PomIncJExact=PomIncJExact/allea/2.d0
3699
3700       return
3701       end
3702
3703
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----------------------------------------------------------------------
3709       include 'epos.inc'
3710       double precision PomIncXUnit,xh,xhm
3711       common /ar3/  x1(7),a1(7)
3712
3713       PomIncJUnit=0.d0
3714       xhm=.5d0
3715       do i=1,7
3716         do m=1,2
3717           xh=xhm*(1.d0+dble(x1(i)*(2.*real(m)-3.)))
3718           PomIncJUnit=PomIncJUnit+dble(a1(i))
3719      &                                *PomIncXUnit(xh,b)
3720         enddo
3721       enddo
3722       PomIncJUnit=PomIncJUnit*xhm
3723
3724       return
3725       end
3726
3727
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)
3737 c
3738 c      imax=1
3739 c      if(iomega.eq.2)imax=1
3740 c      s=engy**2
3741 ccctp060829      sy=s*real(xh)
3742 cc Calculation by numeric integration :
3743 c      w=0.d0
3744 c      ymax=-.5d0*log(xh)
3745 c      do m=1,2
3746 c        do j=1,7
3747 c          xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
3748 c          xm=xh/xp
3749 c          Df=0.d0
3750 c          do i=idxDmin,imax
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))
3755 c          enddo
3756 c          w=w+dble(a1(j))*Df
3757 c        enddo
3758 c      enddo
3759 c      w=w*ymax*xh**dble(-alppar)
3760 c
3761 c
3762 c      PomIncXExact1=w
3763 c
3764 c      return
3765 c      end
3766
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----------------------------------------------------------------------
3773       include 'epos.inc'
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)
3778
3779       imax=1
3780       if(iomega.eq.2)imax=1
3781       s=engy**2
3782       PomIncXExact=0.d0
3783       do i=idxDmin,imax
3784         call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3785         bep  =betx
3786         bepp =betpx
3787         AlTiP=alpx*gamv
3788         PomIncXExact=PomIncXExact+AlTiP*xh**((bep+bepp)/2.d0)
3789       enddo
3790
3791       factor=0.d0
3792       allea=min(alplea(iclpro),alplea(icltar))+1.
3793       xpmin=max(sqrt(xh),exp(-1.d0))
3794       do i=1,7
3795         do m=1,2
3796           xp=1.d0-(1.d0-xpmin)*(.5d0+dble(x1(i)*(real(m)-1.5)))
3797      *                                         **(1.d0/dble(allea))
3798           xm=xh/xp
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
3804         enddo
3805       enddo
3806       factor=factor*(1.d0-xpmin)**dble(allea)/dble(allea)
3807
3808
3809       if(xpmin.gt.1.00001d0*sqrt(xh))then
3810         ymax=-log(xh)-2.d0
3811         factor1=0.d0
3812         do i=1,7
3813           do m=1,2
3814             y=ymax*dble(x1(i)*(2*m-3))
3815             xp=sqrt(xh*exp(y))
3816             xm=xh/xp
3817             factor1=factor1+dble(a1(i))*(1.d0-xp)**dble(alplea(iclpro))
3818      *                                 *(1.d0-xm)**dble(alplea(icltar))
3819           enddo
3820         enddo
3821         factor=factor+factor1*ymax
3822       endif
3823
3824       factor=factor/2.d0
3825
3826       PomIncXExact=PomIncXExact*factor
3827
3828       return
3829       end
3830
3831
3832 c----------------------------------------------------------------------
3833       double precision function PomIncXUnit(xh,b)   !---check---
3834 c----------------------------------------------------------------------
3835 c incluse  Unitarized Pomeron distribution  \int dy
3836 c----------------------------------------------------------------------
3837       include 'epos.inc'
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)
3842
3843       imax=1
3844       if(iomega.eq.2)imax=1
3845       s=engy**2
3846 ctp060829      sy=s*real(xh)
3847 c Calculation by numeric integration :
3848       w=0.d0
3849       ymax=-.5d0*log(xh)
3850       do m=1,2
3851         do j=1,7
3852           xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
3853           xm=xh/xp
3854           Df=0.d0
3855           do i=idxDmin,imax
3856             call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3857             Df=Df+alp*real(xp)**bet*real(xm)**betp
3858           enddo
3859           call GfomPar(b,s)
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)
3862         enddo
3863       enddo
3864       w=w*ymax
3865
3866
3867       PomIncXUnit=w!/Znorm(s,b)
3868
3869       return
3870       end
3871
3872
3873
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----------------------------------------------------------------------
3879       include 'epos.inc'
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)
3883
3884       eps=1.d-5
3885
3886       imax=idxD1
3887       if(iomega.eq.2)imax=1
3888       do i=idxDmin,imax
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)
3891       enddo
3892       call GfomPar(b,s)
3893
3894       if (1.d0-xp-eps.gt.0.d0.and.1.d0-xm-eps.gt.0.d0) then
3895         w=0.d0
3896         zpmin=1.d0-xp-eps
3897         zmmin=1.d0-xm-eps
3898         zpm=.5d0*zpmin
3899         zmm=.5d0*zmmin
3900         do m=1,2
3901           do i=1,7
3902             zp=zpm+dble((2*m-3)*x1(i))*zpm
3903             do n=1,2
3904               do j=1,7
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)
3907      &                              *PhiUnit(zp,zm)
3908 c     &                             *PhiExact(1.,zp,zm,s,b)
3909              enddo
3910            enddo
3911          enddo
3912        enddo
3913        PoInU=w*zpm*zmm
3914        deltzp=eps
3915        deltzm=eps
3916       else
3917         PoInU=0.d0
3918         zpmin=0.d0
3919         zmmin=0.d0
3920         deltzp=1.d0-xp
3921         deltzm=1.d0-xm
3922       endif
3923
3924       w=0.d0
3925       if(abs(deltzp).gt.1.d-10.and.abs(deltzm).gt.1.d-10)then
3926       zpm=.5d0*deltzp
3927       zmm=.5d0*deltzm
3928       do m=1,2
3929         do i=1,7
3930           zp=zpmin+zpm+dble((2*m-3)*x1(i))*zpm
3931           do n=1,2
3932             do j=1,7
3933               zm=zmmin+zmm+dble((2*n-3)*x1(j))*zmm
3934               zp2=1.d0-xp-zp**2
3935               zm2=1.d0-xm-zm**2
3936               w=w+dble(a1(i)*a1(j))*HrstI(s,b,zp,zm)
3937      &             *PhiUnit(zp2,zm2)
3938 c     &             *PhiExact(1.,zp2,zm2,s,b)
3939             enddo
3940           enddo
3941         enddo
3942       enddo
3943       endif
3944
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)
3948
3949       return
3950       end
3951
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----------------------------------------------------------------------
3957       include "epos.inc"
3958       include "epos.incsem"
3959       double precision Df,xp,xm!,xh
3960
3961       Df=0.d0
3962       xh=xp*xm
3963 ctp060829      sy=engy**2*real(xh)
3964       s=engy**2
3965       imax=1
3966       if(iomega.eq.2)imax=1
3967       do i=idxDmin,imax
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
3970       enddo
3971       Df=dble(chad(iclpro)*chad(icltar))
3972      *  *Df
3973
3974       PomIncExact=Df
3975      *            *(1.d0-xp)**dble(alplea(iclpro))
3976      *            *(1.d0-xm)**dble(alplea(icltar))
3977
3978       return
3979       end
3980
3981 c----------------------------------------------------------------------
3982         double precision function PomIncUnit(xp,xm,b)   !---check---
3983 c----------------------------------------------------------------------
3984 c inclusive Pomeron distribution  { Sum{int{G*Phi} }
3985 c----------------------------------------------------------------------
3986       include "epos.inc"
3987       include "epos.incpar"
3988       include "epos.incsem"
3989       double precision PoInU,xp,xm,om1,xh,yp
3990
3991
3992       xh=xp*xm
3993       yp=0.d0
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)
3997
3998       return
3999       end
4000
4001
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,
4014 c     *                 sxp,sxm,PhiExpo
4015 c
4016 c      PomIncUnitMC=0.d0
4017 c      if(xp.lt.1.d0.and.xm.lt.1.d0)then
4018 c      m=10
4019 c
4020 c      sy=engy*engy
4021 c      nmcint=2000
4022 c      nmax=nmcint
4023 c      do i=1,mmax
4024 c        xxpu(i)=0.d0
4025 c        xxmu(i)=0.d0
4026 c      enddo
4027 c      xprem=1.d0
4028 c      xmrem=1.d0
4029 c      sxp=xprem-xp
4030 c      sxm=xmrem-xm
4031 c
4032 c      Gtab=omNpcut(xp,xm,sxp,sxm,b,0)
4033 c      Phitab=PhiExpo(1.,sxp,sxm,sy,b)
4034 c      Z=Gtab*Phitab
4035 c      Zm=0.d0
4036 c
4037 c      do mtmp=2,m
4038 c
4039 c        write(*,*)"GPhi",mtmp-1,Zm,Z
4040 c        Zm=0.d0
4041 c        n=0
4042 c
4043 c 10     continue
4044 c          n=n+1
4045 c          if(mod(n,1000000).eq.0)write(*,*)
4046 c     &              "Calculation of PomIncUnit(",mtmp,")->",n
4047 c          xxpu(1)=xp
4048 c          xxmu(1)=xm
4049 c          sxp=xxpu(1)
4050 c          sxm=xxmu(1)
4051 c          pGtab=1.d0
4052 c          do i=2,mtmp
4053 c            rnau=rangen()*real(xprem-sxp)
4054 c            xxpu(i)=dble(rnau)
4055 c            sxp=sxp+xxpu(i)
4056 c            rnau=rangen()*real(xmrem-sxm)
4057 c            xxmu(i)=dble(rnau)
4058 c            sxm=sxm+xxmu(i)
4059 c          enddo
4060 c          if(sxp.lt.xprem.and.sxm.lt.xmrem)then
4061 c            do i=1,mtmp
4062 c              Gtab=omNpcut(xxpu(i),xxmu(i),xprem-sxp,xmrem-sxm,b,0)
4063 c              pGtab=pGtab*Gtab
4064 c            enddo
4065 c          Zm=Zm+pGtab*PhiExpo(1.,xprem-sxp,xmrem-sxm,sy,b)
4066 c          endif
4067 c        if(n.lt.nmax)goto 10
4068 c        Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(mtmp)
4069 c        Z=Z+Zm
4070 c      enddo
4071 c
4072 c      PomIncUnitMC=Z/dble(chad(iclpro)*chad(icltar))
4073 c      endif
4074 c
4075 c      return
4076 c      end
4077 c
4078 c
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
4091 cc     *                 ,omNpuncut
4092 c
4093 c      PhiMCExact=0.d0
4094 c      if(xp.le.1.d0.and.xm.le.1.d0)then
4095 c      m=6
4096 c
4097 c      sy=engy*engy
4098 c      nmcint=50000
4099 c      nmax=nmcint
4100 c      do i=1,mmax
4101 c        xxpu(i)=0.d0
4102 c        xxmu(i)=0.d0
4103 c      enddo
4104 c
4105 c      Z=xp**dble(alplea(iclpro))
4106 c     * *xm**dble(alplea(icltar))
4107 c      Zm=0.d0
4108 c
4109 c      do mtmp=1,m
4110 c
4111 c        write(*,*)"GPhi",mtmp-1,Zm,Z/xp**dble(alplea(iclpro))
4112 c     *                              /xm**dble(alplea(icltar))
4113 c        Zm=0.d0
4114 c        n=0
4115 c
4116 c 10     continue
4117 c          n=n+1
4118 c          if(mod(n,1000000).eq.0)write(*,*)
4119 c     &              "Calculation of PhiExact(",mtmp,")->",n
4120 c          sxp=0.d0
4121 c          sxm=0.d0
4122 c          pGtab=1.d0
4123 c          do i=1,mtmp
4124 c            rnau=rangen()!*real(xp-sxp)
4125 c            xxpu(i)=dble(rnau)
4126 c            sxp=sxp+xxpu(i)
4127 c            rnau=rangen()!*real(xm-sxm)
4128 c            xxmu(i)=dble(rnau)
4129 c            sxm=sxm+xxmu(i)
4130 c          enddo
4131 c          if(sxp.lt.xp.and.sxm.lt.xm)then
4132 c            do i=1,mtmp
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))
4136 c              else
4137 c                yp=0.d0
4138 c              endif
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)
4142 c            enddo
4143 c            Zm=Zm+pGtab*(xp-sxp)**dble(alplea(iclpro))
4144 c     *           *(xm-sxm)**dble(alplea(icltar))
4145 c          endif
4146 c        if(n.lt.nmax)goto 10
4147 c        Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(m)
4148 c        Z=Z+Zm
4149 c      enddo
4150 c
4151 c      PhiMCExact=Z
4152 c      endif
4153 c
4154 c      return
4155 c      end
4156
4157 c----------------------------------------------------------------------
4158         double precision function Gammapp(sy,b,mtmp)   !---check---
4159 c----------------------------------------------------------------------
4160       include "epos.inc"
4161       include "epos.incpar"
4162       include "epos.incsem"
4163       include "epos.incems"
4164       parameter(mmax=20)
4165       double precision Gtab,xxpu(mmax),xxmu(mmax),PhiExpo
4166       double precision Zm,xp,xm,pGtab,om1,sxp,sxm,xh,yp
4167
4168       Gammapp=0.d0
4169
4170       xp=1.d0
4171       xm=1.d0
4172       nmcint=20000
4173       nmax=nmcint
4174       do i=1,mmax
4175         xxpu(i)=0.d0
4176         xxmu(i)=0.d0
4177       enddo
4178       Zm=0.d0
4179
4180         n=0
4181
4182  10     continue
4183           n=n+1
4184           if(mod(n,10000).eq.0)write(*,*)
4185      &              "Calculation of Gammapp(",mtmp,")->",n
4186           sxp=0.d0
4187           sxm=0.d0
4188           pGtab=1.d0
4189           do i=1,mtmp
4190             rnau=rangen()!*real(xp-sxp)
4191             xxpu(i)=dble(rnau)
4192             sxp=sxp+xxpu(i)
4193             rnau=rangen()!*real(xm-sxm)
4194             xxmu(i)=dble(rnau)
4195             sxm=sxm+xxmu(i)
4196           enddo
4197           if(sxp.lt.xp.and.sxm.lt.xm)then
4198             do i=1,mtmp
4199               xh=xxpu(i)*xxmu(i)
4200               if(abs(xh).gt.1.d-30)then
4201                 yp=0.5d0*log(xxpu(i)/xxmu(i))
4202               else
4203                 yp=0.d0
4204               endif
4205               Gtab=om1(xh,yp,b)
4206               pGtab=pGtab*Gtab
4207             enddo
4208             Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4209           endif
4210         if(n.lt.nmax)goto 10
4211         Zm=Zm/dble(nmax)!**2.d0*(xp*xm)**dble(mtmp)
4212
4213       Gammapp=Zm
4214
4215       return
4216       end
4217
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
4230 c
4231 c      GammaMCnew=0.d0
4232 c      Zmtot=0.d0
4233 c      xp=1.d0
4234 c      xm=1.d0
4235 c      nmcint=1000
4236 c      nmax=nmcint
4237 c
4238 c
4239 c      do i=1,mmax
4240 c        xxpu(i)=0.d0
4241 c        xxmu(i)=0.d0
4242 c      enddo
4243 c
4244 c      Zm=0.d0
4245 c
4246 c        n=0
4247 c
4248 c 10     continue
4249 c          n=n+1
4250 c          if(mod(n,1000000).eq.0)write(*,*)
4251 c     &              "Calculation of GammaMCnew(",mtmp,")->",n
4252 c          sxp=0.d0
4253 c          sxm=0.d0
4254 c          pGtab=1.d0
4255 c          do i=1,mtmp
4256 c            rnau=rangen()
4257 c            xxpu(i)=dble(rnau)
4258 c            sxp=sxp+xxpu(i)
4259 c            rnau=rangen()
4260 c            xxmu(i)=dble(rnau)
4261 c            sxm=sxm+xxmu(i)
4262 c          enddo
4263 c          if(sxp.lt.xp.and.sxm.lt.xm)then
4264 c            i=0
4265 c            do k=1,mtmp
4266 c                i=i+1
4267 c                Gtab=omGam(xxpu(i),xxmu(i),b) !om1(xh,yp,b)
4268 c                pGtab=pGtab*Gtab
4269 c              enddo
4270 c            Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4271 c          endif
4272 c        if(n.lt.nmax)goto 10
4273 c
4274 c          Zmtot=Zmtot+Zm/dble(nmax)
4275 c
4276 c      GammaMCnew=Zmtot
4277 c
4278 c      return
4279 c      end
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
4291 c
4292 c      GammaMC=0.d0
4293 c
4294 c      xp=1.d0
4295 c      xm=1.d0
4296 c      nmcint=50000
4297 c      nmax=nmcint
4298 c      do i=1,mmax
4299 c        xxpu(i)=0.d0
4300 c        xxmu(i)=0.d0
4301 c      enddo
4302 c
4303 c      Zm=0.d0
4304 c
4305 c        n=0
4306 c
4307 c 10     continue
4308 c          n=n+1
4309 c          if(mod(n,1000000).eq.0)write(*,*)
4310 c     &              "Calculation of GammaMC(",mtmp,")->",n
4311 c          sxp=0.d0
4312 c          sxm=0.d0
4313 c          pGtab=1.d0
4314 c          do i=1,mtmp
4315 c            rnau=rangen()!*real(xp-sxp)
4316 c            xxpu(i)=dble(rnau)
4317 c            sxp=sxp+xxpu(i)
4318 c            rnau=rangen()!*real(xm-sxm)
4319 c            xxmu(i)=dble(rnau)
4320 c            sxm=sxm+xxmu(i)
4321 c          enddo
4322 c          if(sxp.lt.xp.and.sxm.lt.xm)then
4323 c            do i=1,mtmp
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))
4327 c              else
4328 c                yp=0.d0
4329 c              endif
4330 c              Gtab=om1(xh,yp,b)!omNpcut(xxpu(i),xxmu(i),xp-sxp,xm-sxm,b,0) !om1(xh,yp,b)
4331 c              pGtab=pGtab*Gtab
4332 c            enddo
4333 c            Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4334 c          endif
4335 c        if(n.lt.nmax)goto 10
4336 c        Zm=Zm/dble(nmax)*fctrl(n-mtmp)*facto(mtmp)
4337 c
4338 c      GammaMC=Zm
4339 c
4340 c      return
4341 c      end
4342 c
4343 c
4344 c
4345
4346 c----------------------------------------------------------------------
4347         double precision function GammaGauss(sy,b,mtmp)   !---check---
4348 c----------------------------------------------------------------------
4349       include "epos.inc"
4350       include "epos.incpar"
4351       include "epos.incsem"
4352       include "epos.incems"
4353       parameter(mmax=3)
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
4359 c     *,PhiExact
4360       common /ar3/  x1(7),a1(7)
4361 c      double precision dgx1,dga1
4362 c      common /dga20/ dgx1(10),dga1(10)
4363
4364       GammaGauss=0.d0
4365       xp=1.d0
4366       xm=1.d0
4367       xmin=1.d-13
4368       eps=1.d-15
4369
4370       if(mtmp.eq.0)then
4371         nmax1=0
4372         jmax1=0
4373         nmax2=0
4374         jmax2=0
4375         nmax3=0
4376         jmax3=0
4377       elseif(mtmp.eq.1)then
4378         nmax1=2
4379         jmax1=7
4380         nmax2=0
4381         jmax2=0
4382         nmax3=0
4383         jmax3=0
4384       elseif(mtmp.eq.2)then
4385         nmax1=2
4386         jmax1=7
4387         nmax2=2
4388         jmax2=7
4389         nmax3=0
4390         jmax3=0
4391       elseif(mtmp.eq.3)then
4392         nmax1=2
4393         jmax1=7
4394         nmax2=2
4395         jmax2=7
4396         nmax3=2
4397         jmax3=7
4398       else
4399         write(*,*)"m not between 0 and 3, return ..."
4400         return
4401       endif
4402
4403         xpmin=xmin
4404         xmmin=xmin
4405         xpmax=1.d0
4406         xmmax=1.d0
4407       do i=1,mmax
4408         zp(i)=0.d0
4409         zm(i)=0.d0
4410         dzp(i)=0.d0
4411         dzm(i)=0.d0
4412         zmmin(i)=0.d0
4413         zpmax(i)=0.d0
4414         zpmin(i)=0.d0
4415         zmmax(i)=0.d0
4416       enddo
4417
4418         G0=1.d0
4419
4420         if(mtmp.eq.0)then
4421           sxp=xp
4422           sxm=xm
4423           G0=PhiExpo(1.,sxp,sxm,sy,b)
4424         endif
4425
4426
4427 c        write(*,*)'x+/-',xmmin,xmmax,xpmin,xpmax
4428
4429
4430         dzm(1)=0.d0
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))
4437         endif
4438         dzm(1)=(zmmax(1)-zmmin(1))/2.d0
4439         endif
4440
4441         dzp(1)=0.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))
4448         endif
4449         dzp(1)=(zpmax(1)-zpmin(1))/2.d0
4450         endif
4451
4452 c        write(*,*)'bornes1=',exp(-zpmax(1)),exp(-zpmin(1))
4453 c     &,exp(-zmmax(1)),exp(-zmmin(1))
4454
4455         Gp1=0.d0
4456         do np1=1,nmax1
4457         do jp1=1,jmax1
4458           zp(1)=zpmin(1)+dzp(1)*(1.d0+dble(2.*np1-3.)*dble(x1(jp1)))
4459           Gm1=0.d0
4460           if(dzm(1).eq.0.d0)then
4461             nmax1=1
4462             jmax1=1
4463           endif
4464           do nm1=1,nmax1
4465             do jm1=1,jmax1
4466               if(dzm(1).ne.0.d0)then
4467               zm(1)=zmmin(1)+dzm(1)*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1)))
4468               else
4469               zm(1)=zp(1)
4470               endif
4471
4472               if(mtmp.eq.1)then
4473               sxp=xp
4474               sxm=xm
4475               do i=1,mtmp
4476                 sxp=sxp-exp(-zp(i))
4477                 sxm=sxm-exp(-zm(i))
4478               enddo
4479               pGtab=1.d0
4480               k=0
4481               do l=1,mtmp
4482                 k=k+1
4483                 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4484                   Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4485                   pGtab=pGtab*Gst
4486                   if(Gst.eq.0.d0)
4487      &                write(*,*)'j1=',k,exp(-zp(k)),exp(-zm(k))
4488      &     ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1
4489                 else
4490                   pGtab=0.d0
4491                   write(*,*)'error1 ?',dzp(k),dzm(k)
4492                 endif
4493               enddo
4494               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4495                 if(dzm(1).ne.0.d0)then
4496                   Gm1=Gm1+pGtab*
4497      &dble(a1(jm1))*PhiExpo(1.,sxp,sxm,sy,b)
4498      &*exp(-zm(1))
4499 c     &dble(a1(jm1))*PhiExact(1.,sxp,sxm,sy,b)
4500 c     &*exp(-zm(1))
4501                   else
4502                   Gp1=Gp1+pGtab*
4503      &dble(a1(jp1))*PhiExpo(1.,sxp,sxm,sy,b)
4504      &*exp(-zp(1))
4505 c     &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b)
4506 c     &*exp(-zp(1))
4507                 endif
4508 c          write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm
4509               endif
4510               endif
4511
4512         dzp(2)=0.d0
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)
4520         endif
4521           dzp(2)=(zpmax(2)-zpmin(2))/2.d0
4522         endif
4523
4524         dzm(2)=0.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)
4532           endif
4533           dzm(2)=(zmmax(2)-zmmin(2))/2.d0
4534         endif
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))
4538
4539         Gp2=0.d0
4540         do np2=1,nmax2
4541         do jp2=1,jmax2
4542           zp(2)=zpmin(2)+dzp(2)*(1.d0+dble(2.*np2-3.)*dble(x1(jp2)))
4543           Gm2=0.d0
4544           if(dzm(2).eq.0.d0)then
4545             nmax2=1
4546             jmax2=1
4547           endif
4548           do nm2=1,nmax2
4549             do jm2=1,jmax2
4550               if(dzm(2).ne.0.d0)then
4551               zm(2)=zmmin(2)+dzm(2)*(1.d0+dble(2.*nm2-3.)*dble(x1(jm2)))
4552               else
4553               zm(2)=zp(2)
4554               endif
4555
4556               if(mtmp.eq.2)then
4557               sxp=xp
4558               sxm=xm
4559               do i=1,mtmp
4560                 sxp=sxp-exp(-zp(i))
4561                 sxm=sxm-exp(-zm(i))
4562               enddo
4563               pGtab=1.d0
4564               k=0
4565               do l=1,mtmp
4566                 k=k+1
4567                 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4568                Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4569                   pGtab=pGtab*Gst
4570                   if(Gst.eq.0.d0)
4571      &                write(*,*)'j2=',k,exp(-zp(k)),exp(-zm(k))
4572      &     ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2
4573                 else
4574                   pGtab=0.d0
4575                   write(*,*)'error2 ?',dzp(k),dzm(k)
4576                 endif
4577               enddo
4578               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4579                 if(dzm(2).ne.0.d0)then
4580                   Gm2=Gm2+pGtab*
4581      &dble(a1(jm2))*PhiExpo(1.,sxp,sxm,sy,b)
4582      &*exp(-zm(2))
4583 c     &dble(a1(jm2))*PhiExact(1.,sxp,sxm,sy,b,mk)
4584 c     &*exp(-zm(2))
4585                   else
4586                   Gp2=Gp2+pGtab*
4587      &dble(a1(jp2))*PhiExpo(1.,sxp,sxm,sy,b)
4588      &*exp(-zp(2))
4589 c     &dble(a1(jp2))*PhiExact(1.,sxp,sxm,sy,b,mk)
4590 c     &*exp(-zp(2))
4591                 endif
4592 c          write(*,*)'m=2',mtmp,Gm2,Gp2,pGtab,sxp,sxm
4593               endif
4594               endif
4595
4596         dzp(3)=0.d0
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
4601         endif
4602
4603         dzm(3)=0.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
4608         endif
4609
4610 c        write(*,*)'bornes3=',exp(-zpmax(3)),exp(-zpmin(3))
4611 c     &,exp(-zmmax(3)),exp(-zmmin(3))
4612
4613         Gp3=0.d0
4614         do np3=1,nmax3
4615         do jp3=1,jmax3
4616           zp(3)=zpmin(3)+dzp(3)*(1.d0+dble(2.*np3-3.)*dble(x1(jp3)))
4617           Gm3=0.d0
4618           if(dzm(3).eq.0.d0)then
4619             nmax3=1
4620             jmax3=1
4621           endif
4622           do nm3=1,nmax3
4623             do jm3=1,jmax3
4624               if(dzm(3).ne.0.d0)then
4625               zm(3)=zmmin(3)+dzm(3)*(1.d0+dble(2.*nm3-3.)*dble(x1(jm3)))
4626               else
4627               zm(3)=zp(3)
4628               endif
4629
4630               sxp=xp
4631               sxm=xm
4632               do i=1,mtmp
4633                 sxp=sxp-exp(-zp(i))
4634                 sxm=sxm-exp(-zm(i))
4635               enddo
4636               pGtab=1.d0
4637               k=0
4638               do l=1,mtmp
4639                 k=k+1
4640                 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4641                Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4642                   pGtab=pGtab*Gst
4643                   if(Gst.eq.0.d0)
4644      &                write(*,*)'j3=',k,exp(-zp(k)),exp(-zm(k))
4645      &   ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2,jp3
4646                 else
4647                   pGtab=0.d0
4648                   write(*,*)'error3 ?',k,dzp(k),dzm(k)
4649                 endif
4650               enddo
4651               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4652                 if(dzm(3).ne.0.d0)then
4653                   Gm3=Gm3+pGtab
4654      &*dble(a1(jm3))*PhiExpo(1.,sxp,sxm,sy,b)
4655      &*exp(-zm(3))
4656                   else
4657                   Gp3=Gp3+pGtab
4658      &*dble(a1(jp3))*PhiExpo(1.,sxp,sxm,sy,b)
4659      &*exp(-zp(3))
4660                 endif
4661               endif
4662             enddo
4663           enddo
4664          if(dzm(3).ne.0.d0)Gp3=Gp3+Gm3*dble(a1(jp3))*exp(-zp(3))*dzm(3)
4665          nmax3=2
4666          jmax3=7
4667         enddo
4668       enddo
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)
4673               endif
4674             enddo
4675           enddo
4676          if(dzm(2).ne.0.d0)Gp2=Gp2+Gm2*dble(a1(jp2))*exp(-zp(2))*dzm(2)
4677          nmax2=2
4678          jmax2=7
4679         enddo
4680       enddo
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)
4685               endif
4686             enddo
4687           enddo
4688          if(dzm(1).ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*exp(-zp(1))*dzm(1)
4689          nmax1=2
4690          jmax1=7
4691         enddo
4692       enddo
4693
4694       if(mtmp.gt.0)G0=Gp1*dzp(1)
4695       write(*,*)"int:",G0
4696
4697       GammaGauss=GammaGauss+G0
4698
4699       return
4700
4701       end
4702
4703 c-----------------------------------------------------------------------
4704       double precision function omWi(sy,b)   !---check---
4705 c-----------------------------------------------------------------------
4706 c cut enhanced diagram integrated over xp, xm, xpr,xmr
4707 c (with ideal G)
4708 c b - impact parameter between the pomeron ends;
4709 c sy- total energy
4710 c-----------------------------------------------------------------------
4711       include "epos.inc"
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
4720 c     *,PhiExact
4721       common /ar3/  x1(7),a1(7)
4722 c      double precision dgx1,dga1
4723 c      common /dga20/ dgx1(10),dga1(10)
4724
4725       omWi=0.d0
4726
4727         xmin=1.d-30
4728         eps=1.d-15
4729         chg=1.d0/dble(delx)
4730         b2=b*b
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
4736
4737         nmax=2
4738         jmax=7
4739
4740         xpmin=xmin
4741         xmmin=xmin
4742         xpmax=1.d0
4743         xmmax=1.d0
4744         zpmin=0.d0
4745         zmmin=0.d0
4746         zpmax=0.d0
4747         zmmax=0.d0
4748         zp=0.d0
4749         zm=0.d0
4750         dzp=0.d0
4751         dzm=0.d0
4752
4753
4754
4755         do intp=1,2
4756         do intm=1,2
4757
4758           if(intp.eq.1)then
4759
4760           xpmin=xmin
4761           xpmax=chg
4762           alpp=(1.d0+2.d0*dble(gamb*b2/2.))
4763
4764           else
4765
4766           xpmin=chg
4767           xpmax=1.d0
4768           alpp=1.d0!(1.d0+2.d0*dble(gamb*b2/2.))
4769
4770           endif
4771
4772           if(intm.eq.1)then
4773
4774           xmmin=xmin
4775           xmmax=chg
4776           alpm=(1.d0+2.d0*dble(gamb*b2/2.))
4777
4778           else
4779
4780           xmmin=chg
4781           xmmax=1.d0
4782           alpm=1.d0!(1.d0+2.d0*dble(gamb*b2/2.))
4783
4784           endif
4785 c        write(*,*)'x+/-',intp,intm,xmmin,xmmax,xpmin,xpmax,alpp,alpm
4786
4787
4788         dzm=0.d0
4789         if(abs(xmmin-xmmax).ge.eps)then
4790           if(alpm.eq.0.d0)then
4791             zmmax=-log(xmmin)
4792             zmmin=-log(xmmax)
4793           else
4794             zmmin=xmmin**alpm
4795             zmmax=xmmax**alpm
4796           endif
4797           dzm=(zmmax-zmmin)/2.d0
4798         endif
4799
4800         dzp=0.d0
4801         if(abs(xpmin-xpmax).ge.eps)then
4802           if(alpp.eq.0.d0)then
4803             zpmax=-log(xpmin)
4804             zpmin=-log(xpmax)
4805           else
4806             zpmin=xpmin**alpp
4807             zpmax=xpmax**alpp
4808           endif
4809           dzp=(zpmax-zpmin)/2.d0
4810         endif
4811
4812
4813         Gp1=0.d0
4814
4815         if(abs(dzp).gt.eps.and.abs(dzm).gt.eps)then
4816 c        write(*,*)'Ca passe ...'
4817
4818         do np1=1,nmax
4819         do jp1=1,jmax
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
4823             xp=exp(-zp)
4824             xjacp=xp
4825           else
4826             xp=zp**(1.d0/alpp)
4827             xjacp=zp**(1.d0/alpp-1.d0)/alpp
4828           endif
4829
4830           Gm1=0.d0
4831           do nm1=1,nmax
4832             do jm1=1,jmax
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
4836                   xm=exp(-zm)
4837                   xjacm=xm
4838                 else
4839                   xm=zm**(1.d0/alpm)
4840                   xjacm=zm**(1.d0/alpm-1.d0)/alpm
4841                 endif
4842
4843               sxp=1.d0-xp
4844               sxm=1.d0-xm
4845               pGtab=1.d0
4846               if(dzp.ge.0.d0.and.dzm.ge.0.d0)then
4847                 pGtab=omGam(xp,xm,b)
4848                 if(pGtab.eq.0.d0)
4849      &          write(*,*)'j1=',xp,xm,xmmin,xmmax,dzp,dzm,jp1
4850               else
4851                 pGtab=0.d0
4852                 write(*,*)'error ?',dzp,dzm
4853               endif
4854               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4855                 if(dzm.ne.0.d0)then
4856                   Gm1=Gm1+pGtab*
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
4860                   else
4861                   Gp1=Gp1+pGtab*
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
4865                 endif
4866 c          write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm
4867               endif
4868
4869             enddo
4870           enddo
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
4873         enddo
4874         enddo
4875         endif
4876
4877         omWi=omWi+Gp1*dzp
4878
4879         enddo
4880         enddo
4881
4882
4883       return
4884       end
4885
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-----------------------------------------------------------------------
4895       include 'epos.inc'
4896       double precision omWi
4897
4898       Womint=omWi(sy,bh)
4899
4900       return
4901       end
4902
4903
4904
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-----------------------------------------------------------------------
4914       include 'epos.inc'
4915       double precision omGamint
4916
4917       WomGamint=omGamint(bh)
4918
4919       return
4920       end
4921