ReggeGribovPartonMCInterface/src/epos-omg.f

0001 c----------------------------------------------------------------------
0002 c The two elementary fuctions of our approach, the profile fuction G
0003 c and the Phi exponent G~, are here referred to as Gpro and Gtilde.
0004 c Both functions can be written as
0005 c
0006 c               G = \sum_type  alp * xp**bet * xm**betp
0007 c
0008 c The parameters alp, bet, betp are calculated in GfunParK (k-mode,
0009 c b is takento the one of pair k) or GfunPar (b-mode: arbitrary b) as
0010 c
0011 c  Gpro:   bet  = betD'  + epsGp + gamD*b**2 + epsG -alppar
0012 c          bet' = betD'' + epsGt + gamD*b**2 + epsG -alppar
0013 c
0014 c          alp  = alpD*f * s**(betD+gamD*b**2+epsG) * exp(-b**2/delD)
0015 c
0016 c  Gtilde: bet~  = bet  + 1
0017 c          bet~' = bet' + 1
0018 c
0019 c          alp~  = alp * gam(bet~)          * gam(bet~')
0020 c                      * gam(1+alppro)      * gam(1+alptar)
0021 c                      * gam(1+alppro+bet~) * gam(1+alptar+bet~')
0022 c                      * (1+epsGt')         * (1+epsGt')
0023 c
0024 c The parameters depend implicitely on s.
0025 c
0026 c In the program we use om1 = Gpro
0027 c  (they differ by a constant which is actually one)
0028 c and om5 = om1 * 0.5
0029 c All functions related to om1 are called om1... .
0030 c
0031 c The inclusive Pomeron distributions are
0032 c
0033 c      PomInc(xp,xm) = Gpro(xp,xm) * (1-xp)**alppro * (1-xm)**alptar
0034 c
0035 c----------------------------------------------------------------------
0036
0037
0038 c----------------------------------------------------------------------
0039       subroutine GfunParK(irea)   !---MC---
0040 c----------------------------------------------------------------------
0041 c  calculates parameters alp,bet,betp of the G functions (k-mode)
0042 c  and the screening exponents epsilongp(k,i), epsilongt(k,i), epsilongs(k,i)
0043 c----------------------------------------------------------------------
0044 c  Gpro parameters written to /comtilde/atilde(,)btildep(,),btildepp(,)
0045 c Gtilde parameters written to /cgtilde/atildg(,)btildgp(,),btildgpp(,)
0046 c  two subscripts: first=type, second=collision k
0047 c Certain pieces of this routine are only done if irea is <= or = zero.
0048 c----------------------------------------------------------------------
0049       include 'epos.inc'
0050       include 'epos.incsem'
0051       include 'epos.incems'
0052       include 'epos.incpar'
0053       double precision atildg,btildgp,btildgpp
0054       common/cgtilde/atildg(idxD0:idxD1,kollmx)
0055      *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx)
0056       double precision utgam2,coefgdp,coefgdt
0057       parameter(nbkbin=40)
0058       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
0059       common /cgtilnu/ cfbetpnp,cfbetppnp,cfbetpnm,cfbetppnm,cfalpro
0060      &,cfaltar,cfbpap,cfbpam,cfbppap,cfbppam
0061       double precision cfbetpnp,cfbetppnp,cfbetpnm,cfbetppnm,cfalpro
0062      &,cfaltar,cfbpap,cfbpam,cfbppap,cfbppam,gamv,eps
0063       parameter (eps=1.d-25)
0064       parameter(nxeps=20,nyeps=32)
0065       common/cxeps1/w(0:nxeps,nyeps),y1(nyeps),y2(nyeps)
0066       common/cxeps2/db,b1,b2
0067       common/geom/rmproj,rmtarg,bmax,bkmx
0068       common/nucl3/phi,bimp
0069       common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
0070      *            ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
0071
0072       b1=0.03
0073       b2=bkmx*1.2
0074       db=(b2-b1)/nyeps
0075       call utprj('GfunParK ',ish,ishini,10)
0076
0077       cfbetpnp=0.d0
0078       cfbetppnp=0.d0
0079       cfbetpnm=0.d0
0080       cfbetppnm=0.d0
0081       cfalpro=dble(ucfpro)
0082       cfaltar=dble(ucftar)
0083       cfbpap=1.d0
0084       cfbppap=1.d0
0085       cfbpam=1.d0
0086       cfbppam=1.d0
0087       alpfom=0.
0088       sy=engy*engy
0089
0090       do k=1,koll
0091         do i=ntymi,ntymx
0092           atilde(i,k)=0.d0
0093           btildep(i,k)=0.d0
0094           btildepp(i,k)=0.d0
0095         enddo
0096         do i=idxD0,idxD1
0097           atildg(i,k)=0.d0
0098           btildgp(i,k)=0.d0
0099           btildgpp(i,k)=0.d0
0100         enddo
0101       enddo
0102
0103
0104 ! calculate collision number according to Glauber -----------------------
0105
0106       bglaub=sqrt(sigine/10./pi)
0107       nglevt=0
0108       do ko=1,koll
0109         r=bk(ko)
0110         if(r.le.bglaub)nglevt=nglevt+1
0111       enddo
0112
0113 ! Z_parton_projectile (zparpro) and Z_parton_target (zpartar)-----------
0114
0115       ztav=0.
0116       zpav=0.
0117       if(iscreen.eq.1.or.isplit.eq.1)then
0118
0119         b2x=b2xscr
0120         alpfom=alpfomi*fegypp
0121         rho0p=conrho(1,0.)
0122         rho0t=conrho(2,0.)
0123         bcut=0.
0124
0125         do k=1,koll
0126           ip=iproj(k)
0127           it=itarg(k)
0128         !....... targ partons seen by proj
0129           if(lproj(ip).gt.0)then
0130             absb=max(1.e-9,bk(k))
0131             b2=absb*absb
0132             zkp=fegypp*exp(-b2/b2x)
0133             zpartar(k)=min(zkp,epscrx)
0134             if(lproj3(ip).gt.1)then
0135               do lt=1,lproj3(ip)
0136                 kp=kproj3(ip,lt)
0137                 if(kp.ne.k)then
0138                   ikt=itarg(kp)
0139                   rtarg=sqrt(xtarg(ikt)**2+ytarg(ikt)**2+ztarg(ikt)**2)
0140                   rho=conrho(2,rtarg)/rho0t
0141                   fegyAA=epscrw*fscro(engy/egyscr,rho)
0142                   absb=max(1.e-9,abs(bk(kp))-bcut)
0143                   b2=absb*absb
0144                   zkp=fegyAA*exp(-b2/b2x)
0145                   zpartar(k)=zpartar(k)+min(zkp,epscrx)
0146                 endif
0147 c                alpfom=max(alpfom,dble(zpartar(k)))
0148               enddo
0149             endif
0150           else
0151             zpartar(k)=0.
0152           endif
0153           ztav=ztav+zpartar(k)
0154          !...........proj partons seen by targ
0155           if(ltarg(it).gt.0)then
0156             absb=max(1.e-9,bk(k))
0157             b2=absb*absb
0158             zkt=fegypp*exp(-b2/b2x)
0159             zparpro(k)=min(zkt,epscrx)
0160             if(ltarg3(it).gt.1)then
0161               do lp=1,ltarg3(it)
0162                 kt=ktarg3(it,lp)
0163                 if(kt.ne.k)then
0164                   ikp=iproj(kt)
0165                   rproj=sqrt(xproj(ikp)**2+yproj(ikp)**2+zproj(ikp)**2)
0166                   rho=conrho(1,rproj)/rho0p
0167                   fegyAA=epscrw*fscro(engy/egyscr,rho)
0168                   absb=max(1.e-9,abs(bk(kt))-bcut)
0169                   b2=absb*absb
0170                   zkt=fegyAA*exp(-b2/b2x)
0171                   zparpro(k)=zparpro(k)+min(zkt,epscrx)
0172                 endif
0173 c                alpfom=max(alpfom,dble(zparpro(k)))
0174               enddo
0175             endif
0176           else
0177             zparpro(k)=0.
0178           endif
0179           zpav=zpav+zparpro(k)
0180           xzcutpar(k)=dble(exp(-min(50.,xzcut*(zparpro(k)+zpartar(k)))))
0181         enddo
0182
0183       else                      ! no screening
0184
0185         do k=1,koll
0186           zparpro(k)=0.
0187           zpartar(k)=0.
0188           xzcutpar(k)=1d0
0189         enddo
0190
0191       endif
0192
0193 c calculation of epsilongp epsilongt
0194
0195       if(iscreen.eq.1)then
0196
0197           !ip=0
0198        do k=1,koll
0199           !ipp=ip
0200         epsG=0.
0201         epsilongs(k,0)=0.
0202         epsilongs(k,1)=0.
0203         ip=iproj(k)             !...........projectile
0204        if(lproj(ip).gt.0)then
0205           x=zparpro(k)
0206           epsilongs(k,0)=sign(abs(epscrd)*x
0207      &                       ,epscrd)
0208           epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
0209      &                        epscrs*x)
0210 c     &                       min(epscrx,epscrs*x))
0211           if(sy.ge.sfshlim)then
0212             epsilongp(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar,
0213      &                          epscrh*x)
0214 c     &                          min(epscrx,epscrh*x))
0215           else
0216             epsilongp(k,1)=epsilongp(k,0)
0217 c            epsilongs(k,1)=epsilongs(k,0)
0218           endif
0219           gammaV(k)=1.d0
0220         else
0221          epsilongp(k,0)=0.
0222          epsilongp(k,1)=0.
0223          gammaV(k)=1.d0
0224         endif
0225
0226         it=itarg(k)             !...........target
0227         if(ltarg(it).gt.0)then
0228           x=zpartar(k)
0229           epsilongs(k,1)=sign(abs(epscrd)*x
0230      &                       ,epscrd)
0231           epsilongt(k,0)=max(-betDpp(0,iclpro,icltar)-0.95+alppar,
0232      &                        epscrs*x)
0233 c     &                        min(epscrx,epscrs*x))
0234           if(sy.ge.sfshlim)then
0235             epsilongt(k,1)=max(-betDpp(1,iclpro,icltar)-0.95+alppar,
0236      &                          epscrh*x)
0237 c     &                          min(epscrx,epscrh*x))
0238           else
0239             epsilongt(k,1)=epsilongt(k,0)
0240 c            epsilongs(k,1)=epsilongs(k,0)
0241           endif
0242 cc          gammaV(k)=gammaV(k)
0243         else
0244          epsilongt(k,0)=0.
0245          epsilongt(k,1)=0.
0246          gammaV(k)=gammaV(k)
0247         endif
0248
0249        enddo
0250
0251       else                      ! no screening
0252
0253        do k=1,koll
0254         epsilongs(k,0)=0.
0255         epsilongs(k,1)=0.
0256         epsilongp(k,0)=0.
0257         epsilongp(k,1)=0.
0258         epsilongt(k,0)=0.
0259         epsilongt(k,1)=0.
0260         gammaV(k)=1.d0
0261        enddo
0262
0263       endif
0264
0265
0266 !..............alpha beta betap for Gtilde (->PhiExpo).......................
0267
0268       imax=idxD1
0269       if(iomega.eq.2)imax=1
0270
0271       do k=1,koll
0272
0273         b=bk(k)
0274         b2=bk(k)*bk(k)
0275         ip=iproj(k)
0276         it=itarg(k)
0277
0278         if(b.lt.(nbkbin-1)*bkbin)then
0279           ibk=int(bk(k)/bkbin)+1
0280           if(isetcs.gt.1.and.iclegy.lt.iclegy2)then
0281             egy0=egylow*egyfac**float(iclegy-1)
0282             xkappa1=xkappafit(iclegy,iclpro,icltar,ibk)
0283      *         +(bk(k)-bkbin*float(ibk-1))/bkbin
0284      *         *(xkappafit(iclegy,iclpro,icltar,ibk+1)
0285      *         -xkappafit(iclegy,iclpro,icltar,ibk))
0286             xkappa2=xkappafit(iclegy+1,iclpro,icltar,ibk)
0287      *         +(bk(k)-bkbin*float(ibk-1))/bkbin
0288      *         *(xkappafit(iclegy+1,iclpro,icltar,ibk+1)
0289      *         -xkappafit(iclegy+1,iclpro,icltar,ibk))
0290             xkappa=xkappa1+(xkappa2-xkappa1)/log(egyfac)
0291      *         *(log(engy)-log(egy0))
0292             xkappa=facmc*xkappa
0293           else
0294             xkappa=xkappafit(iclegy,iclpro,icltar,ibk)
0295      *         +(bk(k)-bkbin*float(ibk-1))/bkbin
0296      *         *(xkappafit(iclegy,iclpro,icltar,ibk+1)
0297      *         -xkappafit(iclegy,iclpro,icltar,ibk))
0298             xkappa=facmc*xkappa
0299           endif
0300         else
0301           xkappa=1.
0302         endif
0303         gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
0304         gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
0305
0306         do i=idxDmin,imax
0307           gamV=gammaV(k)
0308 c          if(i.lt.2)then
0309 c          if(i.eq.0)then
0310 c            epsG=epsilongs(k,i)
0311           if(i.eq.2)then
0312             if(epscrd.lt.0.)then
0313               epsG=epsilongs(k,0)+epsilongs(k,1)
0314               epsGp=0.
0315               epsGt=0.
0316             else
0317               epsG=0.
0318               epsGp=epsilongs(k,0)
0319               epsGt=epsilongs(k,1)
0320             endif
0321           else
0322             epsG=0.
0323             epsGp=0.
0324             epsGt=0.
0325           endif
0326           gamb=gamD(i,iclpro,icltar)*b2
0327           atildg(i,k)=dble(alpD(i,iclpro,icltar))
0328      *            *cfalpro*cfaltar
0329      *            *gamv
0330 c          if(i.eq.0) atildg(i,k)=atildg(i,k)
0331           atildg(i,k)=atildg(i,k)
0332      *            *dble(xkappa*xkappa)
0333           if(i.lt.2)then
0334             atildg(i,k)=atildg(i,k)*dble(
0335      *            chad(iclpro)*chad(icltar)
0336      *            *exp(-b2/delD(i,iclpro,icltar))
0337      *            *sy**(betD(i,iclpro,icltar)+gamb+epsG)
0338      *            *gfactorp *gfactort)
0339             epsGp=epsilongp(k,i)
0340             epsGt=epsilongt(k,i)
0341             btildgp(i,k)=dble(betDp(i,iclpro,icltar)
0342      *                    +epsGp
0343      *                    +gamb-alppar)+1.d0
0344             btildgpp(i,k)=dble(betDpp(i,iclpro,icltar)
0345      *                    +epsGt
0346      *                    +gamb-alppar)+1.d0
0347           else
0348             absb=abs(bk(k))-bmxdif(iclpro,icltar)
0349             b2a=absb*absb
0350             atildg(i,k)=atildg(i,k)*dble(
0351      *                  sy**(betD(i,iclpro,icltar)+epsG)
0352      *                 *exp(-b2a/delD(i,iclpro,icltar)))
0353             btildgp(i,k)=dble(betDp(i,iclpro,icltar)-alppar+epsGp)+1.d0
0354           btildgpp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar+epsGt)+1.d0
0355           endif
0356           coefgdp=utgam2(1.d0+dble(alplea(iclpro))+btildgp(i,k))
0357           coefgdt=utgam2(1.d0+dble(alplea(icltar))+btildgpp(i,k))
0358           atildg(i,k)=atildg(i,k)
0359      *            *utgam2(btildgp(i,k))*utgam2(btildgpp(i,k))
0360      *            /coefgdp/coefgdt
0361    !...........prepare plot in xepsilon
0362           if(irea.eq.0)then
0363            kk=max(1,int((bk(k)-b1)/db)+1)
0364            if(i.lt.2)then
0365              if(i.eq.0)w(0,kk)=w(0,kk)+1
0366              if(i.eq.0)w(1,kk)=w(1,kk)+epsGp
0367              if(i.eq.0)w(2,kk)=w(2,kk)+epsGt
0368 c                     w(3+i,kk)=w(3+i,kk)+abs(epsG)
0369          !...5-8 soft ... 9-12 semi
0370              w(5+4*i,kk)=w(5+4*i,kk)
0371      *              +betDp(i,iclpro,icltar)   !prj eff
0372      *              +epsGp+gamb
0373              w(6+4*i,kk)=w(6+4*i,kk)
0374      *              +betDpp(i,iclpro,icltar)  !tgt eff
0375      *              +epsGt+gamb
0376              w(7+4*i,kk)=w(7+4*i,kk)
0377      *              +betDp(i,iclpro,icltar)  !prj unscr
0378      *              +gamb
0379              w(8+4*i,kk)=w(8+4*i,kk)
0380      *              +betDpp(i,iclpro,icltar) !tgt unscr
0381      *              +gamb
0382              if(i.eq.0)w(13,kk)=w(13,kk)+zparpro(k)
0383              if(i.eq.0)w(14,kk)=w(14,kk)+zpartar(k)
0384            else
0385              if(epscrd.lt.0.)then
0386                w(3,kk)=w(3,kk)+abs(epsG)
0387              else
0388                w(3,kk)=w(3,kk)+epsGp
0389                w(4,kk)=w(4,kk)+epsGt
0390              endif
0391            endif
0392           endif
0393    !................
0394         enddo
0395       enddo
0396
0397 !...........................................................................
0398
0399            zppevt=zpav/koll
0400            zptevt=ztav/koll
0401            if(irea.eq.0)then
0402              ktot=int(bimp)+1
0403              if(ktot.le.nyeps)then
0404                w(15,ktot)=w(15,ktot)+zppevt
0405                w(16,ktot)=w(16,ktot)+zptevt
0406                w(17,ktot)=w(17,ktot)+1
0407              endif
0408              n=1+int(float(nglevt)/(0.1*maproj*matarg))*(nyeps-1)
0409              if(nglevt.ge.1.and.n.ge.1.and.n.le.nyeps)then
0410                w(18,n)=w(18,n)+zppevt
0411                w(19,n)=w(19,n)+zptevt
0412                w(20,n)=w(20,n)+1
0413              endif
0414            endif
0415
0416
0417 !........alpha beta betap for Gpro...........................................
0418
0419       if(irea.le.0)then
0420       do k=1,koll
0421         ip=iproj(k)
0422         it=itarg(k)
0423
0424         b2=bk(k)*bk(k)
0425         imax=ntymx
0426         if(iomega.eq.2)imax=1
0427         do i=ntymin,imax
0428
0429 c          if(i.lt.2)then
0430 c          if(i.eq.0)then
0431 c            epsG=epsilongs(k,i)
0432           if(i.eq.2)then
0433             if(epscrd.lt.0.)then
0434               epsG=epsilongs(k,0)+epsilongs(k,1)
0435               epsGp=0.
0436               epsGt=0.
0437             else
0438               epsG=0.
0439               epsGp=epsilongs(k,0)
0440               epsGt=epsilongs(k,1)
0441             endif
0442           else
0443             epsG=0.
0444             epsGp=0.
0445             epsGt=0.
0446           endif
0447           gamb=gamD(i,iclpro,icltar)*b2
0448
0449           atilde(i,k)=dble(alpD(i,iclpro,icltar))
0450           if(i.lt.2)then
0451           atilde(i,k)=atilde(i,k)*dble(
0452      *              exp(-b2/delD(i,iclpro,icltar))
0453      *              *sy**(betD(i,iclpro,icltar)
0454      *                    +gamb+epsG)
0455      *              *chad(iclpro)*chad(icltar))
0456             epsGp=epsilongp(k,i)
0457             epsGt=epsilongt(k,i)
0458             btildep(i,k)=dble(betDp(i,iclpro,icltar)
0459      *                    +epsGp
0460      *                    +gamb-alppar)
0461             btildepp(i,k)=dble(betDpp(i,iclpro,icltar)
0462      *                    +epsGt
0463      *                    +gamb-alppar)
0464           else
0465             absb=abs(bk(k))-bmxdif(iclpro,icltar)
0466             b2a=absb*absb
0467             atilde(i,k)=atilde(i,k)*dble(
0468      *                  sy**(betD(i,iclpro,icltar)+epsG)
0469      *                 *exp(-b2a/delD(i,iclpro,icltar)))
0470
0471             btildep(i,k)=dble(betDp(i,iclpro,icltar)-alppar+epsGp)
0472             btildepp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar+epsGt)
0473           endif
0474
0475           if(btildep(i,k)+1d0.lt.-eps.or.btildepp(i,k)+1d0.lt.-eps)then
0476             write(ifmt,*)' k,b,ip,it,gamb,alppar',k,bk(k),ip,it,gamb
0477      *                                           ,alppar
0478             write(ifmt,*)' gammaV,epsGP1/2,epsGT1/2,epsGS1/2'
0479      *           ,gammaV(k),epsilongp(k,0),epsilongp(k,1)
0480      *     ,epsilongt(k,0),epsilongt(k,1),epsilongs(k,0),epsilongs(k,1)
0481             write(ifmt,*)'*******************************************'
0482             write(ifmt,*)" atilde,btildep,btildepp "
0483             do ii=ntymin,ntymx
0484               write(ifmt,*)ii,atilde(ii,k),btildep(ii,k),btildepp(ii,k)
0485             enddo
0486             call utstop('Error in epos-omg in GfunPark&')
0487           endif
0488         enddo
0489
0490       enddo
0491       endif
0492
0493 !...........................................................................
0494
0495       if(ish.ge.10)then
0496       do k=1,koll
0497         ip=iproj(k)
0498         it=itarg(k)
0499         write(ifch,*)' k,b,ip,it,',k,bk(k),ip,it
0500         write(ifch,*)' zparpro,zpartar,xzcutpar'
0501      *      ,zparpro(k),zpartar(k),xzcutpar(k)
0502         write(ifch,*)' gammaV,epsilonGP1/2,epsilonGT1/2,epsilonGs1/2'
0503      *      ,gammaV(k),epsilongp(k,0),epsilongp(k,1)
0504      *      ,epsilongt(k,0),epsilongt(k,1),epsilongs(k,0),epsilongs(k,1)
0505         write(ifch,*)'*******************************************'
0506         write(ifch,*)" atilde,btildep,btildepp "
0507         do i=ntymin,ntymx
0508         write(ifch,*)i,atilde(i,k),btildep(i,k),btildepp(i,k)
0509         enddo
0510         write(ifch,*)" atildg,btildgp,btildgpp "
0511         do i=ntymin,ntymx
0512         write(ifch,*)i,atildg(i,k),btildgp(i,k),btildgpp(i,k)
0513         enddo
0514         call GfunPar(xpar7,xpar7,1,0,bk(k),sy,alp,bet,betp
0515      &                  ,epsp,epst,epss,gamvv)
0516         call GfunPar(xpar7,xpar7,1,1,bk(k),sy,alp,bet,betp
0517      &                  ,epsp,epst,epss,gamvv)
0518       enddo
0519       endif
0520
0521       call utprjx('GfunParK ',ish,ishini,10)
0522
0523       return
0524       end
0525
0526 c----------------------------------------------------------------------
0527       subroutine GfunPar(zzip,zzit,m,i,b,spp,alp,bet,betp,epsp,epst
0528      &                  ,epss,gamvv)
0529 c----------------------------------------------------------------------
0530 c  calculates parameters alp,bet,betp of the G functions for pp (b-mode)
0531 c  and the screening exponents epsp,epst,epss,gamvv
0532 c----------------------------------------------------------------------
0533 c  zzip:additional z component (nuclear effect projectile side)
0534 c  zzit:additional z component (nuclear effect target side)
0535 c  m=1: profile function Gpro,  i=0: soft       i=2: diff
0536 c  m=2: Gtilde,                 i=1: semi       (no screening for diff)
0537 c  b: impact param, spp: pp energy squared
0538 c----------------------------------------------------------------------
0539
0540       include 'epos.inc'
0541       include 'epos.incsem'
0542       include 'epos.incpar'
0543       include 'epos.incems'
0544       parameter(nbkbin=40)
0545       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
0546       common /kwrite/ xkapZ
0547       double precision utgam2,coefgdp,coefgdt,dalp,dbet,dbetp,eps
0548       parameter(eps=1.d-20)
0549
0550       call utprj('GfunPar ',ish,ishini,10)
0551
0552       ee=sqrt(spp)
0553 c      bglaub2=FbGlaub2(ee)
0554       rs=r2had(iclpro)+r2had(icltar)+slopom*log(spp)
0555       bglaub2=4.*.0389*rs
0556       if(ish.ge.10)write(ifch,*)'Gf in:',m,i,b,bglaub2,spp,zzip,zzit
0557      &                                    ,iclpro,icltar
0558       b2=b*b
0559       cfalpro=ucfpro
0560       cfaltar=ucftar
0561       gamb=gamD(i,iclpro,icltar)*b2
0562
0563       if(iscreen.ne.0)then
0564         absb=max(1.e-9,abs(b))
0565         b2a=absb*absb
0566         b2x=2.*epscrp*bglaub2
0567         zzp=epscrw*exp(-b2a/b2x)*fscra(ee/egyscr)
0568         zzp=min(zzp,epscrx)+zzip !saturation
0569         zzt=epscrw*exp(-b2a/b2x)*fscra(ee/egyscr)
0570         zzt=min(zzt,epscrx)+zzit !saturation
0571
0572         x=zzp
0573         epsG=0.
0574         if(i.eq.1.and.spp.ge.sfshlim)then
0575           epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar,
0576      &               epscrh*x)
0577 c     &               min(epscrx,epscrh*x))
0578         elseif(i.le.1)then
0579           epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar,
0580      &               epscrs*x)
0581 c     &               min(epscrx,epscrs*x))
0582         else
0583           if(epscrd.lt.0.)then
0584             epsG=epsG+sign(abs(epscrd)*x,epscrd)
0585             epsGp=0.
0586           else
0587             epsGp=sign(abs(epscrd)*x,epscrd)
0588             epsG=0.
0589           endif
0590         endif
0591         gamV=1.
0592         x=zzt
0593         if(i.eq.1.and.spp.ge.sfshlim)then
0594           epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar,
0595      &               epscrh*x)
0596 c     &               min(epscrx,epscrh*x))
0597         elseif(i.le.1)then
0598           epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar,
0599      &               epscrs*x)
0600 c     &               min(epscrx,epscrs*x))
0601         else
0602           if(epscrd.lt.0.)then
0603             epsG=epsG+sign(abs(epscrd)*x,epscrd)
0604             epsGt=0.
0605           else
0606             epsGt=sign(abs(epscrd)*x,epscrd)
0607             epsG=0.
0608           endif
0609         endif
0610 c        gamV=gamV
0611       else
0612         zzp=0.
0613         zzt=0.
0614         epsGp=0.
0615         epsGt=0.
0616         epsG=0.
0617         gamV=1.
0618       endif
0619
0620       gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
0621       gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub)
0622
0623       rho=betD(i,iclpro,icltar)+gamb+epsG
0624
0625
0626       if(m.eq.1)then
0627
0628         dalp=dble(alpD(i,iclpro,icltar))
0629         if(i.lt.2)then
0630           dalp=dalp
0631      *       *exp(min(50d0,dble(rho*log(spp)-b2/delD(i,iclpro,icltar))))
0632           dbet=dble(betDp(i,iclpro,icltar)
0633      *         +epsGp
0634      *         +gamb-alppar)
0635           dbetp=dble(betDpp(i,iclpro,icltar)
0636      *         +epsGt
0637      *         +gamb-alppar)
0638         else
0639           absb=abs(b)-bmxdif(iclpro,icltar)
0640           b2a=absb*absb
0641           dalp=dalp
0642      *       *exp(min(50d0,dble((betD(i,iclpro,icltar)+epsG)*log(spp)
0643      *            -b2a/delD(i,iclpro,icltar))))
0644           dbet=dble(betDp(i,iclpro,icltar)-alppar+epsGp)
0645           dbetp=dble(betDpp(i,iclpro,icltar)-alppar+epsGt)
0646         endif
0647
0648         if((dbet+1.d0).lt.-eps.or.(dbetp+1.d0).lt.-eps)then
0649           write(*,*)'m,i,b,spp,alp,bet,betp',m,i,b,spp,dalp,dbet,dbetp
0650           call utstop('Error : beta < -1 in Gfunpar in epos-omg&')
0651         endif
0652
0653       elseif(m.eq.2)then
0654         xkappa=1.
0655 c        if(i.eq.0.and.b.lt.(nbkbin-1)*bkbin)then
0656         if(b.lt.(nbkbin-1)*bkbin)then
0657           ibk=int(b/bkbin)+1
0658           if(isetcs.gt.1.and.iclegy.lt.iclegy2)then
0659             egy0=egylow*egyfac**float(iclegy-1)
0660             xkappa1=xkappafit(iclegy,iclpro,icltar,ibk)
0661      *         +(b-bkbin*float(ibk-1))/bkbin
0662      *           *(xkappafit(iclegy,iclpro,icltar,ibk+1)
0663      *            -xkappafit(iclegy,iclpro,icltar,ibk))
0664             xkappa2=xkappafit(iclegy+1,iclpro,icltar,ibk)
0665      *         +(b-bkbin*float(ibk-1))/bkbin
0666      *           *(xkappafit(iclegy+1,iclpro,icltar,ibk+1)
0667      *            -xkappafit(iclegy+1,iclpro,icltar,ibk))
0668             xkappa=xkappa1+(xkappa2-xkappa1)/log(egyfac)
0669      *         *(log(ee)-log(egy0))
0670             xkappa=facmc*xkappa
0671           else
0672             xkappa=xkappafit(iclegy,iclpro,icltar,ibk)
0673      *         +(b-bkbin*float(ibk-1))/bkbin
0674      *           *(xkappafit(iclegy,iclpro,icltar,ibk+1)
0675      *            -xkappafit(iclegy,iclpro,icltar,ibk))
0676             xkappa=facmc*xkappa
0677           endif
0678           xkapZ=xkappa
0679         endif
0680
0681         dalp=dble(alpD(i,iclpro,icltar)
0682      *        *cfalpro*cfaltar
0683      *        *gamV)
0684 c        if(i.eq.0)alp=alp
0685         dalp=dalp
0686      *        *dble(xkappa)*dble(xkappa)
0687
0688         if(i.lt.2)then
0689           dalp=dalp
0690      *     *exp(min(50d0,dble(rho*log(spp)-b2/delD(i,iclpro,icltar))))
0691      *        *dble(chad(iclpro)*chad(icltar)
0692      *        *gfactorp *gfactort)
0693           dbet=dble(betDp(i,iclpro,icltar)
0694      *        +epsGp
0695      *        +gamb-alppar+1.)
0696           dbetp=dble(betDpp(i,iclpro,icltar)
0697      *        +epsGt
0698      *        +gamb-alppar+1.)
0699         else
0700           absb=abs(b)-bmxdif(iclpro,icltar)
0701           b2a=absb*absb
0702           dalp=dalp
0703      *     *exp(min(50d0,dble((betD(i,iclpro,icltar)+epsG)*log(spp)
0704      *          -b2a/delD(i,iclpro,icltar))))
0705           dbet=dble(betDp(i,iclpro,icltar)-alppar+1.+epsGp)
0706           dbetp=dble(betDpp(i,iclpro,icltar)-alppar+1.+epsGt)
0707         endif
0708         coefgdp=utgam2(1.d0+dble(alplea(iclpro))+dbet)
0709         coefgdt=utgam2(1.d0+dble(alplea(icltar))+dbetp)
0710         dalp=dalp*utgam2(dbet)*utgam2(dbetp)/coefgdp/coefgdt
0711       else
0712
0713         stop'GproPar: wrong m value.              '
0714
0715       endif
0716
0717
0718       alp=sngl(dalp)
0719       bet=sngl(dbet)
0720       betp=sngl(dbetp)
0721       epsp=epsGp
0722       epst=epsGt
0723       epss=epsG
0724       gamvv=gamV
0725
0726       alpUni(i,m)=dalp
0727       betUni(i,m)=dbet
0728       betpUni(i,m)=dbetp
0729
0730
0731       if(ish.ge.10)write(ifch,*)'   GfunPar :',alp,bet,betp,epsp,epst
0732      &                                    ,epss,gamvv
0733
0734       call utprjx('GfunPar ',ish,ishini,10)
0735       end
0736
0737 c----------------------------------------------------------------------
0738       subroutine GfomPar(b,spp)
0739 c----------------------------------------------------------------------
0740 c  calculates parameters of the fom functions for pp (b-mode)
0741 c----------------------------------------------------------------------
0742 c  b: impact param, spp: pp energy squared
0743 c----------------------------------------------------------------------
0744
0745       include 'epos.inc'
0746       include 'epos.incsem'
0747       include 'epos.incpar'
0748       include 'epos.incems'
0749
0750       call utprj('GfomPar ',ish,ishini,6)
0751
0752       ee=sqrt(spp)
0753       rs=r2had(iclpro)+r2had(icltar)+slopom*log(spp)
0754       bglaub2=4.*.0389*rs
0755
0756       if(iscreen.ne.0)then
0757         absb=max(1.e-9,abs(b))
0758         b2a=absb*absb
0759         b2x=2.*epscrp*bglaub2
0760         zzp=(epscrw*exp(-b2a/b2x))*fscra(ee/egyscr)
0761         zzp=min(zzp,epscrx) !saturation
0762         zzt=(epscrw*exp(-b2a/b2x))*fscra(ee/egyscr)
0763         zzt=min(zzt,epscrx) !saturation
0764       else
0765         zzp=0.
0766         zzt=0.
0767       endif
0768
0769
0770       z0=alpfomi!*epscrw*fscra(ee/egyscr)
0771       if(z0.gt.0.)then
0772         z1=zzp
0773         zzpUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
0774 c      zzpUni=dble(4.*z0*(z1/z0)**1.5)
0775         z1=zzt
0776         zztUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar)))
0777 c      zztUni=dble(4.*z0*(z1/z0)**1.5)
0778       else
0779         zzpUni=0d0
0780         zztUni=0d0
0781       endif
0782
0783       if(ish.ge.6)write(ifch,*)'   GfomPar :',zzpUni,zztUni
0784
0785       call utprjx('GfomPar ',ish,ishini,6)
0786       end
0787
0788 c----------------------------------------------------------------------
0789       function fscra(x)
0790 c----------------------------------------------------------------------
0791       fscra=0
0792       x2=x*x
0793       if(x2.gt.1.)fscra=log(x2)!**2
0794       end
0795
0796 c----------------------------------------------------------------------
0797       function fscro(x,rho)
0798 c----------------------------------------------------------------------
0799       include 'epos.incpar'
0800       fscro=znurho*rho
0801       x2=x*x
0802 c      if(x2.gt.1.)fscro=sqrt(log(x2)**2+fscro**2)
0803       if(x2.gt.1.)fscro=log(x2)*(1.+fscro)
0804       end
0805
0806 c----------------------------------------------------------------------
0807       function FbGlaub2(x)
0808 c----------------------------------------------------------------------
0809 c  calculates (glauber radius)^2 from pp cross section (data fit)
0810 c(estimated if not already calculated --> not any more to have smoother xs)
0811 c----------------------------------------------------------------------
0812 c  x: pp energy
0813 c----------------------------------------------------------------------
0814
0815       include 'epos.inc'
0816
0817 c      if(sigine.eq.0.)then
0818         if(iclpro+icltar.eq.3)then !pi+p
0819           siginex=20.+0.08*log(x)**3.-0.004*log(x)**4.
0820         elseif(iclpro+icltar.eq.5)then !K+p
0821           siginex=16.+0.08*log(x)**3.-0.004*log(x)**4.
0822         else
0823           siginex=30.+0.095*log(x)**3.-0.004*log(x)**4.
0824         endif
0825 c      else
0826 c       siginex=sigine
0827 c      endif
0828       FbGlaub2=siginex/10./pi
0829
0830       return
0831       end
0832
0833 c----------------------------------------------------------------------
0834       subroutine recalcZPtn !???????????? not updated !!!
0835 c----------------------------------------------------------------------
0836
0837       include 'epos.inc'
0838       include 'epos.incems'
0839        stop'recalcZPtn not valid any more!!!!!!!'
0840 c      if(koll.eq.1.and.maproj.eq.1.and.matarg.eq.1)then
0841 c       npom=nprt(1)
0842 c       k=1
0843 c       ip=iproj(1)
0844 c       it=itarg(1)
0845 c       zparpro(k)=max(0,npom-1)*0.2
0846 c       zpartar(k)=0
0847 c       zpartar(k)=max(0,npom-1)*0.2
0848 c       ztav=zpartar(k)
0849 c       zpav=zparpro(k)
0850 c       zppevt=zpav
0851 c       zptevt=ztav
0852 c      endif
0853       end
0854
0855 c----------------------------------------------------------------------
0856       double precision function om1(xh,yp,b)   !---test---
0857 c----------------------------------------------------------------------
0858 c om1 = G * C * gamd    (C and gamd usually 1)
0859 c xh - fraction of the energy squared s for the pomeron;
0860 c b - impact parameter between the pomeron ends;
0861 c yp - rapidity for the pomeron;
0862 c----------------------------------------------------------------------
0863       include 'epos.inc'
0864       include 'epos.incsem'
0865       include 'epos.incpar'
0866
0867       double precision Gf,xp,xm,xh,yp
0868
0869       Gf=0.d0
0870       xp=sqrt(xh)*exp(yp)
0871       xm=xh/xp
0872       spp=engy**2
0873       imax=idxD1
0874       if(iomega.eq.2)imax=1
0875       do i=idxDmin,imax
0876         call Gfunpar(0.,0.,1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
0877         Gf=Gf+dble(alp)*xp**dble(bet)*xm**dble(betp)
0878       enddo
0879       om1=Gf
0880      *  * dble(chad(iclpro)*chad(icltar))
0881       end
0882
0883 c----------------------------------------------------------------------
0884       double precision function om1intb(b)   !---test---
0885 c----------------------------------------------------------------------
0886 c  om1 integrated over xp and xm for given b
0887 c  Calculation by analytical integration
0888 c----------------------------------------------------------------------
0889       include 'epos.inc'
0890       include 'epos.incsem'
0891       include 'epos.incpar'
0892       double precision cint,cint2,eps
0893       parameter(eps=1.d-20)
0894
0895       spp=engy*engy
0896       imax=idxD1
0897       if(iomega.eq.2)imax=1
0898       cint=0.d0
0899       do i=idxDmin,imax
0900         call Gfunpar(0.,0.,1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
0901         cint2=dble(gamv*alp)
0902         if((bet+1.0).gt.eps)then
0903           cint2=cint2/dble(bet+1.0)
0904         else
0905           cint2=-cint2*log(xminDf)
0906         endif
0907         if((betp+1.0).gt.eps)then
0908           cint2=cint2/dble(betp+1.0)
0909         else
0910           cint2=-cint2*log(xminDf)
0911         endif
0912         cint=cint+cint2
0913       enddo
0914
0915       if(cint.lt.-eps)then
0916         write(*,*) 'WARNING ! om1intb in epos-omg is <0 !!!!!'
0917         write(*,*) 'WARNING ! => om1intb set to 1e-3 !!!!!'
0918         write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
0919         cint=1.d-3
0920       endif
0921
0922       om1intb=cint
0923      *       *dble(chad(iclpro)*chad(icltar))
0924
0925       return
0926       end
0927
0928 c----------------------------------------------------------------------
0929       double precision function om1intbk(k)   !---MC---
0930 c----------------------------------------------------------------------
0931 c  Diffractive part of om1 integrated over xp and xm for given pair k
0932 c  Calculation by analytical integration
0933 c----------------------------------------------------------------------
0934       include 'epos.inc'
0935       include 'epos.incsem'
0936       include 'epos.incems'
0937       include 'epos.incpar'
0938       double precision cint,cint2,eps,bet,betp
0939       parameter(eps=1.d-20)
0940
0941       imax=idxD1
0942       if(iomega.eq.2)imax=1
0943       om1intbk=0d0
0944       cint=0
0945       do i=idxDmin,imax
0946         bet=btildep(i,k)
0947         betp=btildepp(i,k)
0948         cint2=atilde(i,k)
0949         if((bet+1.d0).gt.eps)then
0950           cint2=cint2/(bet+1.d0)
0951         else
0952           cint2=-cint2*log(xminDf)
0953         endif
0954         if((betp+1.d0).gt.eps)then
0955           cint2=cint2/(betp+1.d0)
0956         else
0957           cint2=-cint2*log(xminDf)
0958         endif
0959         cint=cint+cint2
0960       enddo
0961
0962       if(cint.lt.-eps)then
0963         write(*,*) 'WARNING ! om1intbk in epos-omg is <0 !!!!!'
0964         write(*,*) 'WARNING ! => om1intbk set to 1e-3 !!!!!'
0965         write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
0966         cint=1.d-3
0967       endif
0968
0969       om1intbk=cint
0970      *       *dble(chad(iclpro)*chad(icltar))
0971
0972       return
0973       end
0974
0975 c----------------------------------------------------------------------
0976       double precision function om1intbi(b,iqq)   !---MC---
0977 c----------------------------------------------------------------------
0978 c  om1 integrated over xp and xm for given b
0979 c  Calculation by analytical integration of contribution iqq
0980 c----------------------------------------------------------------------
0981       include 'epos.inc'
0982       include 'epos.incsem'
0983       include 'epos.incpar'
0984       double precision eps,cint
0985       parameter(eps=1.d-20)
0986
0987       spp=engy*engy
0988       call Gfunpar(0.,0.,1,iqq,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
0989       cint=dble(gamv*alp)
0990       if(dble(bet+1.0).gt.eps)then
0991         cint=cint/dble(bet+1.0)
0992       else
0993         cint=-cint*log(xminDf)
0994       endif
0995       if(dble(betp+1.0).gt.eps)then
0996         cint=cint/dble(betp+1.0)
0997       else
0998         cint=-cint*log(xminDf)
0999       endif
1000       if(cint.lt.-eps)then
1001         write(*,*) 'WARNING ! om1intbi in epos-omg is <0 !!!!!'
1002         write(*,*) 'WARNING ! => om1intbi set to 1e-3 !!!!!'
1003         write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!'
1004         cint=1.d-3
1005       endif
1006
1007       om1intbi=cint
1008      *       *dble(chad(iclpro)*chad(icltar))
1009
1010       return
1011       end
1012
1013 c----------------------------------------------------------------------
1014       double precision function om1intbc(b)   !---MC---
1015 c----------------------------------------------------------------------
1016 c  om1*F*F integrated over xp and xm for given b
1017 c  Calculation by analytical integration
1018 c----------------------------------------------------------------------
1019       include 'epos.inc'
1020       include 'epos.incems'
1021       include 'epos.incsem'
1022       double precision cint,gamom,deltap,deltam
1023       double precision utgam2,Fp,Fm
1024
1025       spp=engy**2
1026       om1intbc=0.d0
1027       Fp=dble(ucfpro)   !gamma(1+alplea)
1028       Fm=dble(ucftar)
1029
1030       imax=idxD1
1031       if(iomega.eq.2)imax=1
1032
1033       cint=0.d0
1034
1035       do i=idxDmin,imax
1036         call Gfunpar(0.,0.,1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
1037        gamom=dble(alp*gamv*chad(iclpro)*chad(icltar))
1038         deltap=dble(bet)
1039         deltam=dble(betp)
1040         cint=cint+gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0)
1041      &            /utgam2(2.d0+deltap+dble(alplea(iclpro)))
1042      &            /utgam2(2.d0+deltam+dble(alplea(icltar)))
1043       enddo
1044
1045       om1intbc=cint*Fp*Fm
1046
1047       if(om1intbc.lt.-1.d-10)then
1048         write(*,*) 'WARNING ! om1intbc in epos-omg is <0 !!!!!'
1049         write(*,*) 'WARNING ! => om1intbc set to 0. !!!!!'
1050         om1intbc=0.d0
1051       endif
1052
1053       return
1054       end
1055
1056 cc----------------------------------------------------------------------
1057 c      double precision function om1intbci(b,iqq)   !---MC---
1058 cc----------------------------------------------------------------------
1059 cc  om1*F*F integrated over xp and xm for given b and given Pomeron type iqq
1060 cc  Calculation by analytical integration
1061 cc----------------------------------------------------------------------
1062 c      include 'epos.inc'
1063 c      include 'epos.incems'
1064 c      include 'epos.incsem'
1065 c      double precision cint,gamom,deltap,deltam
1066 c      double precision utgam2,Fp,Fm,eps
1067 c
1068 c      spp=engy**2
1069 c      om1intbci=0.d0
1070 c      Fp=dble(ucfpro)   !gamma(1+alplea)
1071 c      Fm=dble(ucftar)
1072 c
1073 c      i=iqq
1074 c      call Gfunpar(0.,0.,1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
1075 c      gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
1076 c      deltap=dble(bet)
1077 c      deltam=dble(betp)
1078 c      cint=gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0)
1079 c     &            /utgam2(2.d0+deltap+dble(alplea(iclpro)))
1080 c     &            /utgam2(2.d0+deltam+dble(alplea(icltar)))
1081 c
1082 c      om1intbci=cint*Fp*Fm
1083 c
1084 c      if(om1intbci.lt.-1.d-10)then
1085 c        write(*,*) 'WARNING ! om1intbci in epos-omg is <0 !!!!!'
1086 c        write(*,*) 'WARNING ! => om1intbci set to 0. !!!!!'
1087 c        om1intbci=0.d0
1088 c      endif
1089 c
1090 c      return
1091 c      end
1092 c
1093 c----------------------------------------------------------------------
1094       double precision function om1intgck(k,xprem,xmrem)   !---MC---
1095 c----------------------------------------------------------------------
1096 c  om1*(xprem-xp)*(xmrem-xm) integrated over xp and xm for given k
1097 c  Calculation by analytical integration
1098 c----------------------------------------------------------------------
1099       include 'epos.inc'
1100       include 'epos.incems'
1101       include 'epos.incsem'
1102       double precision cint,gamom,deltap,deltam,xprem,xmrem
1103
1104       om1intgck=0.d0
1105
1106       imax=idxD1
1107       if(iomega.eq.2)imax=1
1108
1109       cint=0.d0
1110       do i=idxDmin,imax
1111         gamom=dble(atilde(i,k))
1112         deltap=dble(btildep(i,k))
1113         deltam=dble(btildepp(i,k))
1114         cint=cint+gamom/(deltap+1.d0)/(deltam+1.d0)
1115      &            /(2.d0+deltap) /(2.d0+deltam)
1116      &            *xprem**(deltap+2.d0)
1117      &            *xmrem**(deltam+2.d0)
1118       enddo
1119       om1intgck=cint
1120
1121       return
1122       end
1123
1124 c----------------------------------------------------------------------
1125       double precision function om1intgc(b)   !---test---
1126 c----------------------------------------------------------------------
1127 c  om1*(1-xp)*(1-xm) integrated over xp and xm for given b
1128 c  Calculation by analytical integration
1129 c----------------------------------------------------------------------
1130       include 'epos.inc'
1131       include 'epos.incpar'
1132       include 'epos.incsem'
1133       double precision cint,gamom,deltap,deltam,eps
1134       parameter(eps=1.d-20)
1135
1136       spp=engy**2
1137       om1intgc=0.d0
1138
1139       imax=idxD1
1140       if(iomega.eq.2)imax=1
1141
1142
1143
1144       cint=0.d0
1145
1146       do i=idxDmin,imax
1147         call Gfunpar(0.,0.,1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
1148         gamom=dble(alp*gamv*chad(iclpro)*chad(icltar))
1149         deltap=dble(bet)
1150         deltam=dble(betp)
1151         if((deltap+1.d0).gt.eps)then
1152           gamom=gamom/(deltap+1.d0)
1153         else
1154           gamom=-gamom*log(xminDf)
1155         endif
1156         if((deltam+1.d0).gt.eps)then
1157           gamom=gamom/(deltam+1.d0)
1158         else
1159           gamom=-gamom*log(xminDf)
1160         endif
1161         cint=cint+gamom /(2.d0+deltap) /(2.d0+deltam)
1162       enddo
1163       om1intgc=cint
1164
1165       if(om1intgc.lt.eps)then
1166         write(*,*) b,deltap,deltam,gamom
1167         write(*,*) 'WARNING ! om1intgc in epos-omg is <0 !!!!!'
1168         write(*,*) 'WARNING ! => om1intgc set to 0. !!!!!'
1169         om1intgc=0.d0
1170       endif
1171
1172       return
1173       end
1174
1175
1176 c----------------------------------------------------------------------
1177         subroutine integom1(irea)
1178 c----------------------------------------------------------------------
1179 c  om1 integrated over xp and xm for  all b=bk(k) if irea=0
1180 c  result written to :
1181 c    om1int(k)   = om1intb(bk(k))      = \int om1
1182 c    om1intc(k)  = om1intgc(bk(k)...   = \int om1*(1-xp)*(1-xm)
1183 c----------------------------------------------------------------------
1184       include 'epos.inc'
1185       include 'epos.incems'
1186       include 'epos.incsem'
1187       double precision om1intbk,om1intgck,PomIncExactk
1188
1189       if(irea.le.0)then
1190
1191
1192         do k=1,koll
1193
1194           om1int(k)=om1intbk(k) !only diffractive contribution
1195           om1intc(k)=om1intgck(k,1d0,1d0)
1196
1197         enddo
1198
1199         if(irea.eq.0.and.zrminc.gt.0..and.xzcut.gt.0.)then
1200           do k=1,koll
1201             PomInck(k)=PomIncExactk(k)
1202           enddo
1203         endif
1204
1205       endif
1206
1207
1208       return
1209       end
1210
1211
1212 c----------------------------------------------------------------------
1213       double precision function om1xpk(xp,xpr1i,k)   !---test---
1214 c----------------------------------------------------------------------
1215 c \int dxm om1*(1-xp)*(1-xm)   (normalized)
1216 c k - pair indice;
1217 c----------------------------------------------------------------------
1218       include 'epos.inc'
1219       include 'epos.incpar'
1220       include 'epos.incems'
1221       double precision xp,gamomx(ntymi:ntymx),cint,gamom
1222       double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps
1223      &                 ,xpr1,xmr1,xpr1i
1224       parameter(eps=1.d-20)
1225
1226       om1xpk=0.d0
1227       if(xp.ge.xpr1i)return
1228
1229       xpr1=1.d0
1230       xmr1=1.d0
1231       imin=ntymin
1232       imax=ntymx
1233       if(iomega.eq.2)imax=1
1234
1235       do i=imin,imax
1236           deltap(i)=btildep(i,k)
1237           deltam(i)=btildepp(i,k)
1238           gamomx(i)=atilde(i,k)*xmr1**(deltam(i)+2.d0)/(2.d0+deltam(i))
1239           if((deltam(i)+1.d0).gt.eps)then
1240             gamomx(i)=gamomx(i)/(deltam(i)+1.d0)
1241           else
1242             gamomx(i)=-gamomx(i)*log(xminDf)
1243           endif
1244       enddo
1245
1246
1247       cint=0.d0
1248       do i=imin,imax
1249           gamom=gamomx(i)*xpr1**(deltap(i)+2.d0)/(2.d0+deltap(i))
1250           if((deltap(i)+1.d0).gt.eps)then
1251             gamom=gamom/(deltap(i)+1.d0)
1252           else
1253             gamom=-gamom*log(xminDf)
1254           endif
1255           cint=cint+gamom
1256       enddo
1257
1258
1259       do i=imin,imax
1260         om1xpk=om1xpk+gamomx(i)*xp**deltap(i)
1261      &                       *(xpr1-xp)
1262
1263       enddo
1264
1265       om1xpk=om1xpk/cint
1266
1267       return
1268       end
1269
1270
1271 c----------------------------------------------------------------------
1272       double precision function om1xmk(xp,xm,xpr1i,xmr1i,k)   !---test---
1273 c----------------------------------------------------------------------
1274 c om1(xp,xm)*(1-xp)*(1-xm)   (normalized for fixed xp)
1275 c k - pair indice;
1276 c----------------------------------------------------------------------
1277       include 'epos.inc'
1278       include 'epos.incpar'
1279       include 'epos.incems'
1280       double precision xp,xm,gamomx(ntymi:ntymx),cint,gamom
1281       double precision deltam(ntymi:ntymx),eps,xpr1,xmr1,xpr1i,xmr1i
1282       parameter(eps=1.d-20)
1283
1284       om1xmk=0.d0
1285       if(xp.ge.xpr1i)return
1286       if(xm.ge.xmr1i)return
1287       xpr1=1.d0
1288       xmr1=1.d0
1289
1290       imin=ntymin
1291       imax=ntymx
1292       if(iomega.eq.2)imax=1
1293
1294       do i=imin,imax
1295           gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xpr1-xp)
1296           deltam(i)=btildepp(i,k)
1297       enddo
1298
1299       cint=0.d0
1300       do i=imin,imax
1301           gamom=gamomx(i)*xmr1**(deltam(i)+2.d0)/(2.d0+deltam(i))
1302           if((deltam(i)+1.d0).gt.eps)then
1303             gamom=gamom/(deltam(i)+1.d0)
1304           else
1305             gamom=-gamom*log(xminDf)
1306           endif
1307           cint=cint+gamom
1308       enddo
1309
1310
1311       do i=imin,imax
1312         om1xmk=om1xmk+gamomx(i)*xm**deltam(i)*(xmr1-xm)
1313       enddo
1314
1315       om1xmk=om1xmk/cint
1316
1317       return
1318       end
1319
1320 c----------------------------------------------------------------------
1321       double precision function om1xpr(atil,btilp,btilpp
1322      &                                   ,xpremi,xmremi,ir)   !---MC---
1323 c----------------------------------------------------------------------
1324 c Random number generated from the function om1xpk. We solve the equation
1325 c which give om1xprk by Newton-Raphson + secant method.
1326 c k - pair indice;
1327 c ir - 1 to get xp, -1 to get xm (be carrefull to inverse xpremi et xmremi
1328 c when calling with ir=-1)
1329 c----------------------------------------------------------------------
1330       include 'epos.inc'
1331       include 'epos.incpar'
1332       include 'epos.incems'
1333       double precision x,x0,x1,gamomx(ntymi:ntymx),eps,f0t,f1t,f00
1334       double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,drangen,xmrem
1335       double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),xprem
1336       parameter (eps=1.d-20)
1337       double precision xpremi,xmremi,xlmin,atil(ntymi:ntymx)
1338      &                ,btilp(ntymi:ntymx),btilpp(ntymi:ntymx)
1339
1340
1341       om1xpr=0.d0
1342       if(xpremi.gt.1d0.or.xmremi.gt.1d0)return
1343       x0=log(xminDf)
1344       x1=log(xpremi)
1345       xprem=xpremi
1346       xmrem=xmremi
1347       imin=ntymin
1348       imax=ntymx
1349       xlmin=1.d0
1350       xt=0d0
1351       if(iomega.eq.2)imax=1
1352
1353       do i=imin,imax
1354         if(ir.gt.0)then
1355           deltap(i)=btilp(i)
1356           deltam(i)=btilpp(i)
1357         else
1358           deltap(i)=btilpp(i)
1359           deltam(i)=btilp(i)
1360         endif
1361         gamomx(i)=atil(i)*xmrem**(deltam(i)+2.d0)/(2.d0+deltam(i))
1362         if((deltam(i)+1.d0).gt.eps)then
1363           gamomx(i)=gamomx(i)/(deltam(i)+1.d0)
1364         else
1365           xlmin=log(xminDf)
1366           gamomx(i)=-gamomx(i)*xlmin
1367         endif
1368       enddo
1369
1370       f0=0.d0
1371       f1=0.d0
1372       do i=imin,imax
1373         if((deltap(i)+1.d0).gt.eps)then
1374           f0=f0+gamomx(i)
1375      &          *(xprem*exp(x0)**(1.d0+deltap(i))/(1.d0+deltap(i))
1376      &           -exp(x0)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1377           f1=f1+gamomx(i)
1378      &          *(xprem*exp(x1)**(1.d0+deltap(i))/(1.d0+deltap(i))
1379      &           -exp(x1)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1380         else
1381           xlmin=log(xminDf)
1382           f0=f0+gamomx(i)*(xprem*(x0-xlmin)-exp(x0)+xminDf)
1383           f1=f1+gamomx(i)*(xprem*(x1-xlmin)-exp(x1)+xminDf)
1384         endif
1385       enddo
1386       f00=f0
1387       cint=f1-f00
1388       f0=-(f0-f00)/cint
1389       f1=-(f1-f00)/cint
1390       ntry=0
1391  11   ntry=ntry+1
1392       r=drangen(dble(ntry))
1393       f0t=f0+r
1394       f1t=f1+r
1395       if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1396       if(f1t*f0t.ge.eps)then
1397         do i=imin,imax
1398          write(ifmt,*)i,gamomx(i),deltap(i),deltam(i)
1399         enddo
1400         write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry
1401         call utstop('om1xpr (2)&')
1402       endif
1403       f0=f0t
1404       f1=f1t
1405       if(abs(f0).le.eps) then
1406         om1xpr=exp(x0)
1407         return
1408       endif
1409       if(abs(f1).le.eps) then
1410         om1xpr=exp(x1)
1411         return
1412       endif
1413       x=0.5d0*(x1+x0)
1414       deltx=abs(x1-x0)
1415
1416
1417       ntry=0
1418
1419  111  continue
1420       if(ntry.le.1000)then
1421       fx=0.d0
1422       fpx=0.d0
1423       do i=imin,imax
1424         if((deltap(i)+1.d0).gt.eps)then
1425           fx=fx+gamomx(i)
1426      &          *(xprem*exp(x)**(1.d0+deltap(i))/(1.d0+deltap(i))
1427      &           -exp(x)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1428           fpx=fpx+gamomx(i)*exp(x)**deltap(i)*(xprem-exp(x))
1429         else
1430           fx=fx+gamomx(i)*(xprem*(x-xlmin)-exp(x)+xminDf)
1431           fpx=fpx+gamomx(i)*(xprem/exp(x)-1.d0)
1432         endif
1433       enddo
1434       fx=-(fx-f00)/cint+r
1435       fpx=fpx/cint
1436       xt=x-fx/fpx
1437
1438       if (f0*fx.lt.0.D0) then
1439         f1=fx
1440         x1=x
1441       else
1442         f0=fx
1443         x0=x
1444       endif
1445       if ((xt.lt.x0.or.xt.gt.x1).and.abs(f1-f0).gt.eps) then
1446         xt=x1-f1*(x1-x0)/(f1-f0)
1447       endif
1448
1449        else
1450
1451         write(ifmt,*)'Warning in om1xpr, to much try !'
1452
1453       endif
1454
1455
1456       if(abs(x-xt).gt.deltx*0.5d0) then
1457         xt=(x1+x0)*0.5D0
1458       endif
1459       deltx=abs(x-xt)
1460       if(abs(x).gt.eps)then
1461         prec=abs(deltx/x)
1462       else
1463         prec=0d0
1464         call utstop('Problem in om1xpr&')
1465       endif
1466
1467       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1468          x=xt
1469          ntry=ntry+1
1470          goto 111
1471       endif
1472
1473       om1xpr=exp(x)
1474
1475       return
1476       end
1477
1478 c----------------------------------------------------------------------
1479       double precision function om1xprk(k,xpremi,xmremi,ir)   !---MC---
1480 c----------------------------------------------------------------------
1481 c Random number generated from the function om1xpk. We solve the equation
1482 c which give om1xprk by Newton-Raphson + secant method.
1483 c k - pair indice;
1484 c ir - 1 to get xp, -1 to get xm (be carrefull to inverse xpremi et xmremi
1485 c when calling with ir=-1)
1486 c----------------------------------------------------------------------
1487       include 'epos.inc'
1488       include 'epos.incpar'
1489       include 'epos.incems'
1490       double precision x,x0,x1,gamomx(ntymi:ntymx),eps,f0t,f1t,f00
1491       double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,drangen,xmrem
1492       double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),xprem
1493       parameter (eps=1.d-20)
1494       double precision xpremi,xmremi,xlmin
1495
1496
1497       om1xprk=0.d0
1498       x0=log(xmremi)
1499       x1=log(xpremi)
1500       xprem=1.d0
1501       xmrem=1.d0
1502       imin=ntymin
1503       imax=ntymx
1504       xlmin=1.d0
1505       xt=0d0
1506       if(iomega.eq.2)imax=1
1507
1508       do i=imin,imax
1509         if(ir.gt.0)then
1510           deltap(i)=btildep(i,k)
1511           deltam(i)=btildepp(i,k)
1512         else
1513           deltap(i)=btildepp(i,k)
1514           deltam(i)=btildep(i,k)
1515         endif
1516         gamomx(i)=atilde(i,k)*xmrem**(deltam(i)+2.d0)/(2.d0+deltam(i))
1517         if((deltam(i)+1.d0).gt.eps)then
1518           gamomx(i)=gamomx(i)/(deltam(i)+1.d0)
1519         else
1520           xlmin=log(xminDf)
1521           gamomx(i)=-gamomx(i)*xlmin
1522         endif
1523       enddo
1524
1525       f0=0.d0
1526       f1=0.d0
1527       do i=imin,imax
1528         if((deltap(i)+1.d0).gt.eps)then
1529           f0=f0+gamomx(i)
1530      &          *(xprem*exp(x0)**(1.d0+deltap(i))/(1.d0+deltap(i))
1531      &           -exp(x0)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1532           f1=f1+gamomx(i)
1533      &          *(xprem*exp(x1)**(1.d0+deltap(i))/(1.d0+deltap(i))
1534      &           -exp(x1)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1535         else
1536           xlmin=log(xminDf)
1537           f0=f0+gamomx(i)*(xprem*(x0-xlmin)-exp(x0)+xminDf)
1538           f1=f1+gamomx(i)*(xprem*(x1-xlmin)-exp(x1)+xminDf)
1539         endif
1540       enddo
1541       f00=f0
1542       cint=f1-f00
1543       f0=-(f0-f00)/cint
1544       f1=-(f1-f00)/cint
1545       ntry=0
1546  11   ntry=ntry+1
1547       r=drangen(dble(ntry))
1548       f0t=f0+r
1549       f1t=f1+r
1550       if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1551       if(f1t*f0t.ge.eps)then
1552         do i=imin,imax
1553          write(ifmt,*)i,gamomx(i),deltap(i),deltam(i)
1554         enddo
1555         write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k
1556         call utstop('om1xprk (2)&')
1557       endif
1558       f0=f0t
1559       f1=f1t
1560       if(abs(f0).le.eps) then
1561         om1xprk=exp(x0)
1562         return
1563       endif
1564       if(abs(f1).le.eps) then
1565         om1xprk=exp(x1)
1566         return
1567       endif
1568       x=0.5d0*(x1+x0)
1569       deltx=abs(x1-x0)
1570
1571
1572       ntry=0
1573
1574  111  continue
1575       if(ntry.le.1000)then
1576       fx=0.d0
1577       fpx=0.d0
1578       do i=imin,imax
1579         if((deltap(i)+1.d0).gt.eps)then
1580           fx=fx+gamomx(i)
1581      &          *(xprem*exp(x)**(1.d0+deltap(i))/(1.d0+deltap(i))
1582      &           -exp(x)**(2.d0+deltap(i))/(2.d0+deltap(i)))
1583           fpx=fpx+gamomx(i)*exp(x)**deltap(i)*(xprem-exp(x))
1584         else
1585           fx=fx+gamomx(i)*(xprem*(x-xlmin)-exp(x)+xminDf)
1586           fpx=fpx+gamomx(i)*(xprem/exp(x)-1.d0)
1587         endif
1588       enddo
1589       fx=-(fx-f00)/cint+r
1590       fpx=fpx/cint
1591       xt=x-fx/fpx
1592
1593       if (f0*fx.lt.0.D0) then
1594         f1=fx
1595         x1=x
1596       else
1597         f0=fx
1598         x0=x
1599       endif
1600       if ((xt.lt.x0.or.xt.gt.x1).and.abs(f1-f0).gt.eps) then
1601         xt=x1-f1*(x1-x0)/(f1-f0)
1602       endif
1603
1604        else
1605
1606         write(ifmt,*)'Warning in om1xprk, to much try !'
1607
1608       endif
1609
1610
1611       if(abs(x-xt).gt.deltx*0.5d0) then
1612         xt=(x1+x0)*0.5D0
1613       endif
1614       deltx=abs(x-xt)
1615       if(abs(x).gt.eps)then
1616         prec=abs(deltx/x)
1617       else
1618         prec=0d0
1619         call utstop('Problem in om1xprk&')
1620       endif
1621
1622       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1623          x=xt
1624          ntry=ntry+1
1625          goto 111
1626       endif
1627
1628       om1xprk=exp(x)
1629
1630       return
1631       end
1632
1633 c----------------------------------------------------------------------
1634       double precision function om1xmrk(k,xp,xpremi,xmremi,ir)   !---MC---
1635 c----------------------------------------------------------------------
1636 c Random number generated from the function om1xmk. We solve the equation
1637 c which give om1xmrk by Newton-Raphson + secant method.
1638 c k - pair indice;
1639 c ir - 1 to get xm, -1 to get xp (be carrefull to inverse xpremi et xmremi
1640 c when calling with ir=-1)
1641 c----------------------------------------------------------------------
1642       include 'epos.inc'
1643       include 'epos.incpar'
1644       include 'epos.incems'
1645       double precision x,x0,x1,gamomx(ntymi:ntymx),eps,xp,f0t,f1t
1646       double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,f00
1647       double precision deltam(ntymi:ntymx),drangen,xprem,xmrem
1648       double precision xpremi,xmremi,xlmin
1649       parameter (eps=1.d-20)
1650
1651       om1xmrk=0.d0
1652 c      if(xp.ge.xpremi)return
1653
1654       xprem=1.d0
1655       xmrem=1.d0
1656       x0=log(xpremi)
1657       x1=log(xmremi)
1658       imin=ntymin
1659       imax=ntymx
1660       if(iomega.eq.2)imax=1
1661
1662       do i=imin,imax
1663         if(ir.gt.0)then
1664           gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xprem-xp)
1665           deltam(i)=btildepp(i,k)
1666         else
1667           gamomx(i)=atilde(i,k)*xp**btildepp(i,k)*(xprem-xp)
1668           deltam(i)=btildep(i,k)
1669         endif
1670       enddo
1671
1672
1673
1674       f0=0.d0
1675       f1=0.d0
1676       xlmin=0.d0
1677       do i=imin,imax
1678         if(abs(deltam(i)+1.d0).gt.eps)then
1679           f0=f0+gamomx(i)
1680      &          *(xmrem*exp(x0)**(1.d0+deltam(i))/(1.d0+deltam(i))
1681      &           -exp(x0)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1682           f1=f1+gamomx(i)
1683      &          *(xmrem*exp(x1)**(1.d0+deltam(i))/(1.d0+deltam(i))
1684      &           -exp(x1)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1685         else
1686           xlmin=log(xminDf)
1687           f0=f0+gamomx(i)*(xmrem*(x0-xlmin)-exp(x0)+xminDf)
1688           f1=f1+gamomx(i)*(xmrem*(x1-xlmin)-exp(x1)+xminDf)
1689         endif
1690       enddo
1691       f00=f0
1692       cint=f1-f00
1693       f0=-(f0-f00)/cint
1694       f1=-(f1-f00)/cint
1695       ntry=0
1696  11   ntry=ntry+1
1697       r=drangen(dble(ntry))
1698       f0t=f0+r
1699       f1t=f1+r
1700       if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1701       if(f1t*f0t.ge.eps)then
1702         write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry
1703         call utstop('Error(2) in epos-omg in om1xmrk&')
1704       endif
1705       f0=f0t
1706       f1=f1t
1707       if(abs(f0).lt.eps) then
1708         om1xmrk=exp(x0)
1709         return
1710       endif
1711       if(abs(f1).lt.eps) then
1712         om1xmrk=exp(x1)
1713         return
1714       endif
1715       x=0.5d0*(x1+x0)
1716       deltx=abs(x1-x0)
1717
1718
1719       ntry=0
1720       xt=0d0
1721
1722  111  continue
1723       if(ntry.le.1000)then
1724       fx=0.d0
1725       fpx=0.d0
1726       do i=imin,imax
1727         if(abs(deltam(i)+1.d0).gt.eps)then
1728           fx=fx+gamomx(i)
1729      &          *(xmrem*exp(x)**(1.d0+deltam(i))/(1.d0+deltam(i))
1730      &           -exp(x)**(2.d0+deltam(i))/(2.d0+deltam(i)))
1731           fpx=fpx+gamomx(i)*exp(x)**deltam(i)*(xmrem-exp(x))
1732         else
1733           fx=fx+gamomx(i)*(xmrem*(x-xlmin)-exp(x)+xminDf)
1734           fpx=fpx+gamomx(i)*(xmrem/exp(x)-1.d0)
1735         endif
1736       enddo
1737       fx=-(fx-f00)/cint+r
1738       fpx=fpx/cint
1739       xt=x-fx/fpx
1740
1741       if (f0*fx.lt.-eps) then
1742         f1=fx
1743         x1=x
1744       else
1745         f0=fx
1746         x0=x
1747       endif
1748       if ((xt.lt.x0-eps.or.xt.gt.x1+eps).and.abs(f1-f0).gt.eps) then
1749         xt=x1-f1*(x1-x0)/(f1-f0)
1750       endif
1751
1752        else
1753
1754         write(ifmt,*)'Warning in om1xmrk, to much try !'
1755
1756       endif
1757
1758       if(abs(x-xt).gt.deltx*0.5d0) then
1759         xt=(x1+x0)*0.5D0
1760       endif
1761       deltx=abs(x-xt)
1762       if(abs(x).gt.eps)then
1763         prec=abs(deltx/x)
1764       else
1765         prec=0d0
1766         call utstop('Problem in om1xmrk&')
1767       endif
1768
1769       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then
1770          x=xt
1771          ntry=ntry+1
1772          goto 111
1773       endif
1774
1775       om1xmrk=exp(x)
1776
1777       return
1778       end
1779
1780 c----------------------------------------------------------------------
1781       double precision function om1xk(xh,k)   !---test---
1782 c----------------------------------------------------------------------
1783 c \int dxp om1   (normalised)
1784 c k - pair indice;
1785 c----------------------------------------------------------------------
1786       include 'epos.inc'
1787       include 'epos.incpar'
1788       include 'epos.incems'
1789
1790       double precision xh,gamomx(ntymi:ntymx),cint,alpp(ntymi:ntymx)
1791      &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps
1792      &,gamom
1793       parameter(eps=1.d-20)
1794
1795
1796       om1xk=0.d0
1797
1798       imin=ntymin
1799       imax=ntymx
1800       if(iomega.eq.2)imax=1
1801
1802       do i=imin,imax
1803         gamomx(i)=atilde(i,k)
1804         deltap(i)=btildep(i,k)
1805         deltam(i)=btildepp(i,k)
1806
1807         delta(i)=(deltap(i)+deltam(i))*0.5d0
1808         alpp(i)=deltap(i)-deltam(i)
1809
1810       enddo
1811
1812       cint=0.d0
1813       do i=imin,imax
1814         gamom=gamomx(i)
1815         if((deltap(i)+1.d0).gt.eps)then
1816           gamom=gamom/(deltap(i)+1.d0)
1817         else
1818           gamom=-gamom*log(xminDf)
1819         endif
1820         if((deltam(i)+1.d0).gt.eps)then
1821           gamom=gamom/(deltam(i)+1.d0)
1822         else
1823           gamom=-gamom*log(xminDf)
1824         endif
1825         cint=cint+gamom
1826       enddo
1827
1828
1829       do i=imin,imax
1830         if(abs(alpp(i)).gt.eps)then
1831         om1xk=om1xk+gamomx(i)/alpp(i)*xh**deltam(i)*(1.d0-xh**alpp(i))
1832         else
1833         om1xk=om1xk-gamomx(i)*xh**delta(i)*log(xh)
1834         endif
1835       enddo
1836
1837       om1xk=om1xk/cint
1838
1839       return
1840       end
1841
1842 c----------------------------------------------------------------------
1843       double precision function om1yk(xh,yp,k)   !---test---
1844 c----------------------------------------------------------------------
1845 c om1 normalized for fixed xp
1846 c xh - fraction of the energy squared s for the pomeron;
1847 c k - pair indice;
1848 c----------------------------------------------------------------------
1849       include 'epos.inc'
1850       include 'epos.incems'
1851       double precision xh,yp,gamomy(ntymi:ntymx),alpp(ntymi:ntymx),cint
1852       double precision deltap,deltam,eps
1853       parameter(eps=1.d-20)
1854
1855       om1yk=0.d0
1856
1857       imin=ntymin
1858       imax=ntymx
1859       if(iomega.eq.2)imax=1
1860
1861       do i=imin,imax
1862           gamomy(i)=atilde(i,k)
1863           deltap=btildep(i,k)
1864           deltam=btildepp(i,k)
1865
1866         alpp(i)=deltap-deltam
1867         gamomy(i)=gamomy(i)*xh**((deltap+deltam)*0.5d0)
1868
1869       enddo
1870
1871       cint=0.d0
1872       do i=imin,imax
1873         if(abs(alpp(i)).gt.eps)then
1874           cint=cint-gamomy(i)/alpp(i)*xh**(alpp(i)*0.5d0)
1875      &                               *(1.d0-xh**(-alpp(i)))
1876         else
1877           cint=cint-gamomy(i)*log(xh)
1878         endif
1879       enddo
1880
1881       do i=imin,imax
1882         if(abs(alpp(i)).gt.eps)then
1883         om1yk=om1yk+gamomy(i)*exp(alpp(i)*yp)
1884         else
1885         om1yk=om1yk+gamomy(i)
1886         endif
1887       enddo
1888
1889       om1yk=om1yk/cint
1890
1891
1892       return
1893       end
1894
1895 c----------------------------------------------------------------------
1896       double precision function om1xrk(k)   !---test---
1897 c----------------------------------------------------------------------
1898 c Random number generated from the function om1xk. We solve the equation
1899 c which give om1xrk by Newton-Raphson + secant method.
1900 c k - pair indice;
1901 c----------------------------------------------------------------------
1902       include 'epos.inc'
1903       include 'epos.incems'
1904       include 'epos.incpar'
1905
1906       double precision x,x0,x1,gamomx(ntymi:ntymx),eps,prec,drangen
1907       double precision xt,fx,fpx,r,f1,f0,cint,deltx,alpp(ntymi:ntymx)
1908      &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx)
1909      &,gamom,f0t,f1t
1910       parameter (eps=1.d-20)
1911
1912
1913       om1xrk=0.d0
1914
1915       imin=ntymin
1916       imax=ntymx
1917       if(iomega.eq.2)imax=1
1918
1919       do i=imin,imax
1920         gamomx(i)=atilde(i,k)
1921         deltap(i)=btildep(i,k)
1922         deltam(i)=btildepp(i,k)
1923
1924         delta(i)=(deltap(i)+deltam(i))*0.5d0
1925         alpp(i)=deltap(i)-deltam(i)
1926
1927       enddo
1928
1929       cint=0.d0
1930       do i=imin,imax
1931         gamom=gamomx(i)
1932         if((deltap(i)+1.d0).gt.eps)then
1933           gamom=gamom/(deltap(i)+1.d0)
1934         else
1935           gamom=-gamom*log(xminDf)
1936         endif
1937         if((deltam(i)+1.d0).gt.eps)then
1938           gamom=gamom/(deltam(i)+1.d0)
1939         else
1940           gamom=-gamom*log(xminDf)
1941         endif
1942         cint=cint+gamom
1943       enddo
1944
1945       x0=eps
1946       x1=1.d0
1947       f0=0.d0
1948       f1=0.d0
1949       do i=imin,imax
1950
1951         if(abs(alpp(i)).lt.eps)then
1952           if(delta(i)+1.d0.gt.eps)then
1953         f0=f0-gamomx(i)/(delta(i)+1.d0)*x0**(delta(i)+1.d0)
1954      &        *(log(x0)-1.d0/(delta(i)+1.d0))
1955         f1=f1-gamomx(i)/(delta(i)+1.d0)*x1**(delta(i)+1.d0)
1956      &        *(log(x1)-1.d0/(delta(i)+1.d0))
1957           else
1958         f0=f0-0.5d0*gamomx(i)*(log(x0)**2-log(xminDf)**2)
1959         f1=f1-0.5d0*gamomx(i)*(log(x1)**2-log(xminDf)**2)
1960           endif
1961         else
1962           if(abs(deltap(i)+1.d0).gt.eps
1963      &  .and.abs(deltam(i)+1.d0).gt.eps)then
1964         f0=f0+gamomx(i)/alpp(i)*(x0**(deltam(i)+1.d0)/(deltam(i)+1.d0)
1965      &                          -x0**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1966         f1=f1+gamomx(i)/alpp(i)*(x1**(deltam(i)+1.d0)/(deltam(i)+1.d0)
1967      &                          -x1**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1968         elseif(abs(deltap(i)+1.d0).gt.eps)then
1969         f0=f0+gamomx(i)/alpp(i)*(log(x0/xminDf)
1970      &                          -x0**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1971         f1=f1+gamomx(i)/alpp(i)*(log(x1/xminDf)
1972      &                          -x1**(deltap(i)+1.d0)/(deltap(i)+1.d0))
1973         elseif(abs(deltam(i)+1.d0).gt.eps)then
1974         f0=f0-gamomx(i)/alpp(i)*(log(x0/xminDf)
1975      &                          -x0**(deltam(i)+1.d0)/(deltam(i)+1.d0))
1976         f1=f1-gamomx(i)/alpp(i)*(log(x1/xminDf)
1977      &                          -x1**(deltam(i)+1.d0)/(deltam(i)+1.d0))
1978           endif
1979         endif
1980       enddo
1981       f0=-f0/cint
1982       f1=-f1/cint
1983       ntry=0
1984  11   ntry=ntry+1
1985       r=drangen(dble(ntry))
1986       f0t=f0+r
1987       f1t=f1+r
1988       if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11
1989       if(f1t*f0t.ge.eps)then
1990         do i=imin,imax
1991          write(ifmt,*)i,gamomx(i),deltap(i),deltam(i),alpp(i),delta(i)
1992         enddo
1993         write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k
1994         call utstop('om1xrk (1)&')
1995       endif
1996       f0=f0t
1997       f1=f1t
1998 c      if(f1*f0.gt.eps)then
1999 c        call utmsg('om1xrk')
2000 c        write(ifch,*)'Poblem with x0, no root ! --> om1xrk=xminDf'
2001 c        write(ifmt,*)'Poblem with x0, no root ! --> om1xrk=xminDf'
2002 c        write(ifmt,*)f0,f1,cint,r
2003 c        call utmsgf
2004 c        om1xrk=x0
2005 c        return
2006 c      endif
2007       if(abs(f0).lt.eps) then
2008         om1xrk=x0
2009         return
2010       endif
2011       if(abs(f1).lt.eps) then
2012         om1xrk=x1
2013         return
2014       endif
2015 c      x=(x1+x0)*0.5D0
2016       x=sqrt(x1*x0)
2017       deltx=abs(x1-x0)
2018
2019       ntry=0
2020       fx=0.d0
2021       fpx=0.d0
2022       xt=x
2023  111  continue
2024
2025       if(ntry.le.1000)then
2026       fx=0.d0
2027       fpx=0.d0
2028       do i=imin,imax
2029         if(abs(alpp(i)).lt.eps)then
2030           if(delta(i)+1.d0.gt.eps)then
2031         fx=fx-gamomx(i)/(delta(i)+1.d0)*x**(delta(i)+1.d0)
2032      &        *(log(x)-1.d0/(delta(i)+1.d0))
2033         fpx=fpx-gamomx(i)*x**delta(i)*log(x)
2034           else
2035         fx=fx-0.5d0*gamomx(i)*(log(x)**2-log(xminDf)**2)
2036         fpx=fpx-gamomx(i)*log(x)/x
2037           endif
2038         else
2039           if(abs(deltap(i)+1.d0).gt.eps
2040      &  .and.abs(deltam(i)+1.d0).gt.eps)then
2041         fx=fx+gamomx(i)/alpp(i)*(x**(deltam(i)+1.d0)/(deltam(i)+1.d0)
2042      &                          -x**(deltap(i)+1.d0)/(deltap(i)+1.d0))
2043         fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i))
2044         elseif(abs(deltap(i)+1.d0).gt.eps)then
2045         fx=fx+gamomx(i)/alpp(i)*(log(x/xminDf)
2046      &                          -x**(deltap(i)+1.d0)/(deltap(i)+1.d0))
2047         fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i))
2048         elseif(abs(deltam(i)+1.d0).gt.eps)then
2049         fx=fx-gamomx(i)/alpp(i)*(log(x/xminDf)
2050      &                          -x**(deltam(i)+1.d0)/(deltam(i)+1.d0))
2051         fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i))
2052           endif
2053         endif
2054       enddo
2055       fx=-fx/cint+r
2056       fpx=fpx/cint
2057       xt=x-fx/fpx
2058
2059       if (f0*fx.lt.-eps) then
2060         f1=fx
2061         x1=x
2062       else
2063         f0=fx
2064         x0=x
2065       endif
2066       if ((xt.lt.x0-eps.or.xt.gt.x1+eps).and.abs(f1-f0).gt.eps) then
2067         xt=x1-f1*(x1-x0)/(f1-f0)
2068       endif
2069
2070       else
2071
2072         write(ifmt,*)'Warning in om1xrk, to much try !'
2073
2074       endif
2075
2076       if(abs(x-xt).gt.deltx*0.5d0) then
2077         xt=sqrt(x1*x0)
2078       endif
2079       deltx=abs(x-xt)
2080       if(abs(x).gt.eps)then
2081         prec=deltx/x
2082       else
2083         prec=0d0
2084         call utstop('Problem in om1xrk&')
2085       endif
2086
2087       if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000)then
2088          x=xt
2089          ntry=ntry+1
2090          goto 111
2091       endif
2092
2093       om1xrk=x
2094
2095       return
2096       end
2097
2098 c----------------------------------------------------------------------
2099       double precision function om1yrk(xh)   !---test---
2100 c----------------------------------------------------------------------
2101 c Random number generated from the function om1yk(xh). We solve the
2102 c equation which give om1yrk by Newton-Raphson + secant method.
2103 c xh - fraction of the energy squared s for the pomeron;
2104 c k - pair indice;
2105 c----------------------------------------------------------------------
2106       include 'epos.inc'
2107       include 'epos.incems'
2108
2109       double precision xh,r!,y0,y1,y,gamomy(ntymi:ntymx),eps,ymin,prec,yt
2110
2111       r=dble(rangen())
2112
2113       om1yrk=(0.5d0-r)*log(xh)
2114       return
2115
2116       end
2117
2118 c----------------------------------------------------------------------
2119       function ffom12aii(iq,je1,je2)   !---test---
2120 c----------------------------------------------------------------------
2121       include 'epos.inc'
2122       ig=5
2123       xmin=0.01/engy
2124       xmax=1
2125       r2=0
2126       do i2=1,ig
2127       do m2=1,2
2128        xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5))
2129        r1=0
2130        do i1=1,ig
2131        do m1=1,2
2132         xp=xmin+(xmax-xmin)*(.5+tgss(ig,i1)*(m1-1.5))
2133         f=ffom12a(xp,xm,iq,iq,je1,je2)
2134         r1=r1+wgss(ig,i1)*f
2135        enddo
2136        enddo
2137        f=r1*0.5*(xmax-xmin)
2138        r2=r2+wgss(ig,i2)*f
2139       enddo
2140       enddo
2141       ffom12aii=r2*0.5*(xmax-xmin)
2142       end
2143
2144 c----------------------------------------------------------------------
2145       function ffom12ai(xp,iq1,iq2,je1,je2)   !---test---
2146 c----------------------------------------------------------------------
2147       include 'epos.inc'
2148       ig=5
2149       xmin=0.01/engy
2150       xmax=1
2151       r2=0
2152       do i2=1,ig
2153       do m2=1,2
2154        xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5))
2155        f=ffom12a(xp,xm,iq1,iq2,je1,je2)
2156        r2=r2+wgss(ig,i2)*f
2157       enddo
2158       enddo
2159       ffom12ai=r2*0.5*(xmax-xmin)
2160       end
2161
2162 c----------------------------------------------------------------------
2163       function ffom12a(xp,xm,iq1,iq2,je1,je2)   !---test---
2164 c----------------------------------------------------------------------
2165 c
2166 c      2*om52*F*F == PomInc
2167 c
2168 c  xp - xplus
2169 c  xm - xminus
2170 c                              iq=1 .... sea-sea
2171 c  iq1 - min iq                iq=2 .... val-sea
2172 c  iq2 - max iq                iq=3 .... sea-val
2173 c                              iq=4 .... val-val
2174 c  je = emission type (projectile and target side)
2175 c          0 ... no emissions
2176 c          1 ... emissions
2177 c       else ... all
2178 c
2179 c  already b-averaged  (\int d2b /sigine*10)
2180 c----------------------------------------------------------------------
2181       include 'epos.inc'
2182
2183       sy=engy*engy
2184       xh=xm*xp
2185 ctp060829      yp=0.5*log(xp/xm)
2186       ffom12a=0
2187       do i=iq1,iq2
2188        if(i.eq.1)then
2189          ffom12a=ffom12a+2*om52pi(sy*xh,1.,1.,0,je1,je2)
2190        elseif(i.eq.2)then
2191          ffom12a=ffom12a+2*om52pi(sy*xh,xp,1.,1,je1,je2)
2192        elseif(i.eq.3)then
2193          ffom12a=ffom12a+2*om52pi(sy*xh,xm,1.,2,je1,je2)
2194        elseif(i.eq.4)then
2195                ffom12a=ffom12a+2*om52pi(sy*xh,xp,xm,3,je1,je2)
2196        endif
2197       enddo
2198       ffom12a=ffom12a
2199      *           *alpff(iclpro)*xp**betff(1)*(1-xp)**alplea(iclpro)
2200      *           *alpff(icltar)*xm**betff(2)*(1-xm)**alplea(icltar)
2201
2202       end
2203
2204 c----------------------------------------------------------------------
2205       function ffom11a(xp,xm,iq1,iq2)   !---test---
2206 c----------------------------------------------------------------------
2207 c
2208 c      int(db) om1ff /sigine*10
2209 c
2210 c  xp - xplus                  iq=-1 ... fit
2211 c  xm - xminus                 iq=0 .... soft
2212 c                              iq=1 .... gg
2213 c  iq1 - min iq                iq=2 .... qg
2214 c  iq2 - max iq                iq=3 .... gq
2215 c                              iq=4 .... qq
2216 c                              iq=5 .... diff
2217 c----------------------------------------------------------------------
2218       include 'epos.inc'
2219       common/geom/rmproj,rmtarg,bmax,bkmx
2220       ig=5
2221       bmid=bkmx/2.
2222       r=0.d0
2223       do i=1,ig
2224         do m=1,2
2225           bb=bmid*(1.+(2.*m-3)*tgss(ig,i))
2226           f=ffom11(xp,xm,bb,iq1,iq2)
2227           r=r+bb*wgss(ig,i)*f
2228         enddo
2229       enddo
2230       ffom11a=r*2.*pi*bmid  /sigine*10
2231       return
2232       end
2233
2234 c----------------------------------------------------------------------
2235       function ffom11(xp,xm,b,iq1,iq2)   !---test---
2236 c----------------------------------------------------------------------
2237 c
2238 c       2*om5*F*F == PomInc
2239 c
2240 c  xp - xplus                  iq=-1 ... fit
2241 c  xm - xminus                 iq=0 .... soft
2242 c  b - impact parameter        iq=1 .... gg
2243 c  iq1 - min iq                iq=2 .... qg
2244 c  iq2 - max iq                iq=3 .... gq
2245 c                              iq=4 .... qq
2246 c                              iq=5 .... diff
2247 c----------------------------------------------------------------------
2248       include 'epos.inc'
2249       double precision om51
2250
2251       if(xm.ge.0.)then
2252
2253        xh=xm*xp
2254        yp=0.5*log(xp/xm)
2255        ffom11=2.*sngl(om51(dble(xh),dble(yp),b,iq1,iq2))
2256      *     *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
2257
2258       else   !xm integration
2259
2260        ig=5
2261        xmin=0.01/engy
2262        xmax=1
2263        r=0
2264        do i=1,ig
2265        do m=1,2
2266         xmm=xmin*(xmax/xmin)**(.5+tgss(ig,i)*(m-1.5))
2267         xh=xmm*xp
2268         yp=0.5*log(xp/xmm)
2269         f=2.*sngl(om51(dble(xh),dble(yp),b,iq1,iq2))
2270      *     *(1-xmm)**alplea(icltar)*(1-xp)**alplea(iclpro)
2271         r=r+wgss(ig,i)*f*xmm
2272        enddo
2273        enddo
2274        ffom11=r*0.5*log(xmax/xmin)
2275
2276       endif
2277
2278       end
2279
2280 c----------------------------------------------------------------------
2281       double precision function om51(xh,yp,b,iq1,iq2)   !---test---
2282 c----------------------------------------------------------------------
2283 c xh - xplus*xminus     iq=-1 ... fit        (om1 * 0.5)
2284 c yp - rapidity         iq=0 .... soft
2285 c b - impact param      iq=1 .... gg
2286 c iq1 - min iq          iq=2 .... qg
2287 c iq2 - max iq          iq=3 .... gq
2288 c                       iq=4 .... qq
2289 c                       iq=5 .... diff
2290 c----------------------------------------------------------------------
2291       include 'epos.inc'
2292       include 'epos.incsem'
2293       include 'epos.incpar'
2294       double precision xp,xm,xh,yp,om51p,om1
2295       om51=0.d0
2296       if(xh.le.0.d0.or.xh.gt.1.d0)return
2297
2298       sy=engy*engy
2299       xp=sqrt(xh)*exp(yp)
2300       xm=xh/xp
2301
2302       if(iq1.eq.-1.and.iq2.eq.-1)then
2303         om51=0.5d0*om1(xh,yp,b)
2304       elseif(iq1.ge.0)then
2305         om51=0.d0
2306         do i=iq1,iq2
2307           if(i.ne.5)then
2308             i1=min(i,1)
2309           call Gfunpar(0.,0.,1,i1,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2310           om51=om51+om51p(sy*sngl(xh),xh,yp,b,i)
2311      *           *xp**dble(epsp)*xm**dble(epst)
2312      *           *dble(sy)**dble(epss)
2313           else
2314           call Gfunpar(0.,0.,1,2,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2315           om51=om51+0.5d0*dble(alp)*xp**dble(bet)*xm**dble(betp)
2316           endif
2317        enddo
2318       else
2319         stop'om5: choice of iq1 and iq2 is nonsense.     '
2320       endif
2321
2322       end
2323
2324 c----------------------------------------------------------------------
2325       double precision function om5s(sx,xh,yp,b,iq1,iq2)   !---test---
2326 c----------------------------------------------------------------------
2327       include 'epos.inc'
2328       double precision om51
2329       double precision xh,yp
2330       ss=sx/xh
2331       engyx=engy
2332       engy=sqrt(ss)
2333       om5s=om51(xh,yp,b,iq1,iq2)
2334       engy=engyx
2335       end
2336
2337 c----------------------------------------------------------------------
2338       double precision function om5Jk(k,xh,yp,iqq)   !---MC---
2339 c----------------------------------------------------------------------
2340 c partial om5
2341 c xh - fraction of the energy squared s for the pomeron;
2342 c b - impact parameter between the pomeron ends;
2343 c yp - rapidity for the pomeron;
2344 c iqq=0 - soft
2345 c iqq=1 - gg
2346 c iqq=2 - qg
2347 c iqq=3 - gq
2348 c iqq=4 - qq
2349 c iqq=5 - diffractif
2350 c----------------------------------------------------------------------
2351       include 'epos.inc'
2352       include 'epos.incems'
2353       include 'epos.incsem'
2354
2355       double precision xh,yp,om51p
2356       double precision plc,s
2357       common/cems5/plc,s
2358
2359       sy=sngl(s*xh)
2360       b=bk(k)
2361       epsGp=0.
2362       epsGt=0.
2363
2364       if(iqq.ne.5)then
2365         om5Jk=om51p(sy,xh,yp,b,iqq)
2366 c Screening effect on Pomeron type set in WomTy
2367         if(iscreen.ne.0)then
2368           xp=sqrt(xh)*exp(yp)
2369           xm=xh/xp
2370           i1=min(iqq,1)
2371 c use pp screening even for nuclei (avoid to many diffractive collisions)
2372 c          call Gfunpar(0.,0.,1,i1,b,sngl(s),alp,bet,betp,epsp,epst,epss
2373 c     &                 ,gamv)
2374 c          epsGp=epsp
2375 c          epsGt=epst
2376 cc          epsG=epsilongs(k,i1)
2377           if(iqq.eq.0)then
2378             epsGp=epsilongp(k,i1)
2379             epsGt=epsilongt(k,i1)
2380           elseif(gfactor.gt.0.)then
2381             epsGp=epsilongp(k,i1)*exp(-min(50.,gfactor*zparpro(k)))
2382             epsGt=epsilongt(k,i1)*exp(-min(50.,gfactor*zpartar(k)))
2383           endif
2384           om5Jk=om5Jk*xp**dble(epsGp)*xm**dble(epsGt)!*s**dble(epsG))
2385         endif
2386       else
2387         xp=sqrt(xh)*exp(yp)
2388         xm=xh/xp
2389         om5Jk=0.5d0*atilde(2,k)*xp**btildep(2,k)*xm**btildepp(2,k)
2390       endif
2391       return
2392       end
2393
2394 c----------------------------------------------------------------------
2395       double precision function om5J(zzp,zzt,xh,yp,b,iq)   !---test---
2396 c----------------------------------------------------------------------
2397 c xh - xplus*xminus     iq=-1 ... fit        (om1 * 0.5)
2398 c yp - rapidity         iq=0 .... soft
2399 c b - impact param      iq=1 .... gg
2400 c                       iq=2 .... qg
2401 c                       iq=3 .... gq
2402 c zzp - projectile Z    iq=4 .... qq
2403 c zzt - target Z        iq=5 .... diff
2404 c----------------------------------------------------------------------
2405       include 'epos.inc'
2406       include 'epos.incsem'
2407       include 'epos.incpar'
2408       double precision xp,xm,xh,yp,om51p,om1
2409       om5J=0.d0
2410       if(xh.le.0.d0.or.xh.gt.1.d0)return
2411
2412       sy=engy*engy
2413       xp=sqrt(xh)*exp(yp)
2414       xm=xh/xp
2415       epsGp=0.
2416       epsGt=0.
2417
2418       if(iq.eq.-1)then
2419         om5J=0.5d0*om1(xh,yp,b)
2420       elseif(iq.ne.5)then
2421         i1=min(iq,1)
2422         call Gfunpar(zzp,zzt,1,i1,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2423 c        call Gfunpar(0.,0.,1,i1,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2424 c        epsG=epss
2425         if(iq.eq.0)then
2426           epsGp=epsp
2427           epsGt=epst
2428         elseif(gfactor.gt.0.)then
2429           epsGp=epsp*exp(-min(50.,gfactor*zzp))
2430           epsGt=epst*exp(-min(50.,gfactor*zzt))
2431         endif
2432         om5J=om51p(sy*sngl(xh),xh,yp,b,iq)
2433      &      *xp**dble(epsGp)*xm**dble(epsGt)!*sy**dble(epsG)
2434       else
2435         call Gfunpar(zzp,zzt,1,2,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2436 c        call Gfunpar(0.,0.,1,2,b,sy,alp,bet,betp,epsp,epst,epss,gamv)
2437       om5J=0.5d0*alp*xp**dble(bet)*xm**dble(betp)
2438       endif
2439
2440       end
2441
2442 c----------------------------------------------------------------------
2443       double precision function omIgamint(b,iqq)   !---test---
2444 c----------------------------------------------------------------------
2445 c - integrated chi~(b)FF/2 for cut I diagram (simple Pomeron)
2446 c b - impact parameter between the pomeron ends;
2447 c yp - rapidity for the pomeron;
2448 c iqq=0 effective one
2449 c iqq=1 soft
2450 c iqq=2 gg
2451 c----------------------------------------------------------------------
2452       include 'epos.inc'
2453       include 'epos.incems'
2454       include 'epos.incsem'
2455       include 'epos.incpar'
2456
2457       double precision Df
2458
2459       Df=0.d0
2460       sy=engy*engy
2461       omIgamint=0.d0
2462       imax=idxD1
2463       if(iomega.eq.2)imax=1
2464
2465       if(iqq.eq.0)then
2466         coefp=1.+alplea(iclpro)
2467         coeft=1.+alplea(icltar)
2468
2469         do i=idxDmin,imax
2470         call Gfunpar(0.,0.,1,i,b,sy,alpx,betx,betpx,epsp,epst,epss,gamv)
2471           betp=1.+betx
2472           betpp=1.+betpx
2473           Df=alpx*dble(utgam1(betp)*utgam1(betpp)*ucfpro
2474      *         *ucftar/utgam1(betp+coefp)/utgam1(betpp+coeft))
2475           omIgamint=omIgamint+Df
2476         enddo
2477       else
2478         call utstop('Wrong iqq in omIgamint&')
2479       endif
2480
2481       omIgamint=omIgamint
2482      *     *dble(chad(iclpro)*chad(icltar))
2483
2484       omIgamint=omIgamint*0.5d0
2485
2486       return
2487       end
2488
2489 c-----------------------------------------------------------------------
2490       subroutine WomTy(w,xh,yp,k)
2491 c-----------------------------------------------------------------------
2492 c - w(ity) for group iqq of cut enhanced diagram giving
2493 c the probability of the type of the same final state.
2494 c k - pair indice;
2495 c xh - fraction of the energy squared s for the pomeron;
2496 c yp - rapidity for the pomeron;
2497 c xpr,xmr impulsion fraction of remnant
2498 c    ity = 0   - soft
2499 c    ity = 1   - gg
2500 c    ity = 2   - qg
2501 c    ity = 3   - gq
2502 c    ity = 4   - qq
2503 c Modified by gfactor to increase semihard int. probability
2504 c-----------------------------------------------------------------------
2505       include 'epos.inc'
2506       include 'epos.incpar'
2507       include 'epos.incems'
2508       include 'epos.incsem'
2509       doubleprecision xh,yp,om5Jk,w(0:7),ww
2510
2511       do i=0,7
2512         w(i)=0.d0
2513       enddo
2514
2515       do i=0,5
2516         w(i)=om5Jk(k,xh,yp,i)
2517       enddo
2518       if(gfactor.lt.0..and.w(1).gt.xggfit*w(0))then      !??????????????
2519         corfac=exp(-dble(abs(gfactor)*(zparpro(k)+zpartar(k)))
2520      &                          *max(0d0,1d0-sqrt(xggfit*w(0)/w(1))))
2521         ww=w(0)+w(1)     !gg interaction probability
2522         w(0)=w(0)*corfac        !soft suppressed
2523         w(1)=ww-w(0)            !semi-hard increased
2524
2525
2526 c        corfac=float(iotst1) !??????????????????????
2527 c        ww=w(0)+w(1)+w(2)+w(3)+w(4)+w(5)
2528 c        whard=0.
2529 c        do i=1,4
2530 c          w(i)=w(i)*corfac
2531 c          whard=whard+w(i)
2532 c        enddo
2533 c        if(whard.gt.ww)then
2534 c         do i=1,4
2535 c          w(i)=w(i)/whard*ww
2536 c         enddo
2537 c         whard=ww
2538 c        endif
2539 c        w05=w(0)+w(5)
2540 c        if(whard.lt.ww)then
2541 c          w(0)=w(0)/w05*(ww-whard)
2542 c          w(5)=w(5)/w05*(ww-whard)
2543 c        else
2544 c          w(0)=0
2545 c          w(5)=0
2546 c        endif
2547         !write(*,'(2f11.4)')(ww-w05)/ww,(ww-w(0)-w(5))/ww
2548       endif
2549
2550       return
2551       end
2552
2553
2554 c-----------------------------------------------------------------------
2555       double precision function Womegak(xp,xm,xprem,xmrem,k)   !---MC---
2556 c-----------------------------------------------------------------------
2557 c - sum(omGam(xp,xm))*(1-xp)*(1-xm) for group of cut enhanced
2558 c diagram giving the same final state (without nuclear effect).
2559 c xp,xm - fraction of the loght cone momenta of the pomeron;
2560 c k - pair index
2561 c-----------------------------------------------------------------------
2562       include 'epos.inc'
2563       include 'epos.incems'
2564       double precision xp,xm,xprem,xmrem
2565
2566       Womegak=0.d0
2567
2568       imax=ntymx
2569       if(iomega.eq.2)imax=1
2570
2571       do i=ntymin,imax
2572         Womegak=Womegak+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k)
2573       enddo
2574
2575
2576       Womegak=Womegak*(xprem-xp)*(xmrem-xm)
2577
2578       return
2579       end
2580
2581
2582 cc----------------------------------------------------------------------
2583 c      double precision function omNpcut(xp,xm,xprem,xmrem,bh,iqq)   !---test---
2584 cc----------------------------------------------------------------------
2585 cc Sum of all cut diagrams
2586 cc iqq=0 ideal G
2587 cc iqq=1 exact G + diff
2588 cc----------------------------------------------------------------------
2589 c      include "epos.inc"
2590 c      double precision om51,xh,yp,xprem,xmrem,xp,xm!,omYcutI
2591 c
2592 c      omNpcut=0.d0
2593 c      xh=xp*xm
2594 c      if(abs(xh).gt.1.d-10)then
2595 c        yp=0.5d0*log(xp/xm)
2596 c      else
2597 c        yp=0.d0
2598 c      endif
2599 c
2600 c      if(iqq.eq.0)omNpcut=om51(xh,yp,bh,-1,-1)
2601 c      if(iqq.eq.1)omNpcut=om51(xh,yp,bh,0,5)
2602 c
2603 c      omNpcut=omNpcut*2.d0
2604 c
2605 c      return
2606 c      end
2607 c
2608 c----------------------------------------------------------------------
2609       double precision function omGam(xp,xm,bh)   !---test---
2610 c-----------------------------------------------------------------------
2611 c Cut diagram part for calculation of probability distribution
2612 c xp,xm impulsion fraction of remnant
2613 c bh - impact parameter between the pomeron ends;
2614 c-----------------------------------------------------------------------
2615       include "epos.inc"
2616       include "epos.incems"
2617       double precision om51,xp,xm,xh,yp,eps!,omYgam
2618       parameter (eps=1.d-20)
2619
2620       omGam=0.d0
2621       if(xp.lt.eps.or.xm.lt.eps)return
2622       xh=xp*xm
2623       if(abs(xh).gt.1.d-10)then
2624         yp=0.5d0*log(xp/xm)
2625       else
2626         yp=0.d0
2627       endif
2628
2629       omGam=om51(xh,yp,bh,-1,-1)
2630
2631       omGam=2.d0*omGam
2632
2633       return
2634       end
2635
2636 c----------------------------------------------------------------------
2637       double precision function omGamk(k,xp,xm)   !---MC---
2638 c-----------------------------------------------------------------------
2639 c Cut diagram part for calculation of probability distribution (for omega)
2640 c xp,xm impulsion fraction of remnant
2641 c bh - impact parameter between the pomeron ends;
2642 c-----------------------------------------------------------------------
2643       include "epos.inc"
2644       include "epos.incems"
2645       double precision xp,xm
2646       omGamk=0.d0
2647       imax=ntymx
2648       if(iomega.eq.2)imax=1
2649       do i=ntymin,imax
2650         omGamk=omGamk+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k)
2651       enddo
2652
2653       return
2654       end
2655
2656 c----------------------------------------------------------------------
2657       double precision function omGamint(bh)   !---test---
2658 c-----------------------------------------------------------------------
2659 c Integrated cut diagram part for calculation of probability distribution
2660 c bh - impact parameter between the pomeron ends;
2661 c-----------------------------------------------------------------------
2662       include "epos.inc"
2663       double precision omIgamint!,omYgamint
2664
2665       omGamint=2.d0*omIgamint(bh,0)
2666
2667       return
2668       end
2669
2670
2671
2672
2673
2674 c----------------------------------------------------------------------
2675       block data dgdata
2676 c----------------------------------------------------------------------
2677 c constants for numerical integration (gaussian weights)
2678 c----------------------------------------------------------------------
2679       double precision dgx1,dga1
2680       common /dga20/ dgx1(10),dga1(10)
2681
2682
2683       data dgx1/
2684      &   .765265211334973D-01,
2685      &   .227785851141645D+00,
2686      &   .373706088715420D+00,
2687      &   .510867001950827D+00,
2688      &   .636053680726515D+00,
2689      &   .746331906460151D+00,
2690      &   .839116971822219D+00,
2691      &   .912234428251326D+00,
2692      &   .963971927277914D+00,
2693      &   .993128599185095D+00/
2694       data dga1/
2695      &   .152753387130726D+00,
2696      &   .149172986472604D+00,
2697      &   .142096109318382D+00,
2698      &   .131688638449177D+00,
2699      &   .118194531961518D+00,
2700      &   .101930119817233D+00,
2701      &   .832767415767047D-01,
2702      &   .626720483341090D-01,
2703      &   .406014298003871D-01,
2704      &   .176140071391506D-01/
2705
2706       end
2707
2708
2709
2710 c----------------------------------------------------------------------
2711       double precision function Phiexact(zzip,zzit,fj,xp,xm,s,b) !---test---
2712 c----------------------------------------------------------------------
2713 c    Exact expression of the Phi function for pp collision
2714 c    zzip : additionnal component for Z (nuclear effect projectile side)
2715 c    zzit : additionnal component for Z (nuclear effect target side)
2716 c    fj   : overall factor for cross section (elastic or inelastic)
2717 c    xp,xm: momentum fraction
2718 c    s    : energy square
2719 c    b    : impact parameter
2720 c----------------------------------------------------------------------
2721       include 'epos.inc'
2722       include 'epos.incsem'
2723       include 'epos.incems'
2724       double precision al(idxD0:idxD1),betp(idxD0:idxD1)
2725      *,z,xIrst!,ffacto
2726       double precision zp(idxD0:idxD1),Phitmp,betpp(idxD0:idxD1)
2727      *,yp,ym,xm,xp
2728       double precision eps
2729       parameter(eps=1.d-20)
2730       dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
2731
2732       if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in PhiExact"
2733       Phitmp=0.d0
2734
2735       if(xp.gt.eps.and.xm.gt.eps.and.xp.le.1.d0+eps
2736      &   .and.xm.le.1.d0+eps)then
2737
2738
2739        do i=idxD0,idxD1
2740         imax(i)=0
2741         ipr(i)=0
2742         zp(i)=1.d0
2743         al(i)=0.d0
2744         betp(i)=0.d0
2745         betpp(i)=0.d0
2746        enddo
2747
2748        imax0=idxD1
2749        if(iomega.eq.2)imax0=1
2750
2751        do i=idxDmin,imax0
2752         imax(i)=10+max(5,int(log10(s)))
2753         if(b.ge.1.)imax(i)=4+max(3,int(log10(sqrt(s))))
2754         imax(i)=min(30,imax(i))
2755        enddo
2756        Phitmp=0.d0
2757        do i=idxDmin,imax0
2758          call Gfunpar(zzip,zzit,1,i,b,s,alpx,betx,betpx,epsp,epst,epss
2759      *               ,gamv)
2760          betp(i)=dble(betx)+1.d0
2761          betpp(i)=dble(betpx)+1.d0
2762          al(i)=dble(alpx*gamv)
2763      *         *dble(chad(iclpro)*chad(icltar))
2764        enddo
2765
2766        do ipr0=0,imax(0)
2767           ipr(0)=ipr0
2768           zp(0)=1.d0
2769         if (ipr(0).ne.0) zp(0)=(-dble(fj)*al(0))**ipr(0)*facto(ipr(0))
2770         do ipr1=0,imax(1)
2771            ipr(1)=ipr1
2772            zp(1)=1.d0
2773         if (ipr(1).ne.0) zp(1)=(-dble(fj)*al(1))**ipr(1)*facto(ipr(1))
2774         do ipr2=0,imax(2)
2775            ipr(2)=ipr2
2776            zp(2)=1.d0
2777         if (ipr(2).ne.0) zp(2)=(-dble(fj)*al(2))**ipr(2)*facto(ipr(2))
2778           yp=0.d0
2779           ym=0.d0
2780           z=1.d0
2781           isum=0
2782           do i=idxDmin,imax0
2783             yp=yp+dble(ipr(i))*betp(i)
2784             ym=ym+dble(ipr(i))*betpp(i)
2785             isum=isum+ipr(i)
2786             z=z*zp(i)
2787           enddo
2788
2789           z=z*xIrst(1,xp,yp,betp,ipr)
2790           z=z*xIrst(2,xm,ym,betpp,ipr)
2791
2792           Phitmp=Phitmp+z
2793
2794          enddo
2795          enddo
2796         enddo
2797
2798
2799
2800       endif
2801
2802       PhiExact=Phitmp
2803
2804
2805       return
2806       end
2807
2808
2809 c----------------------------------------------------------------------
2810       double precision function PhiExpoK(k,xp,xm)   !---MC---
2811 c----------------------------------------------------------------------
2812 c    Exponential expression of the Phi function for pp collision
2813 c    for given k
2814 c----------------------------------------------------------------------
2815       include 'epos.inc'
2816       include 'epos.incsem'
2817       include 'epos.incems'
2818
2819       double precision xp,xm,Phitmp,Gt1
2820       double precision atildg,btildgp,btildgpp
2821       common/cgtilde/atildg(idxD0:idxD1,kollmx)
2822      *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx)
2823
2824
2825       Phitmp=0.d0
2826
2827       imax=idxD1
2828       if(iomega.eq.2)imax=1
2829
2830       Phitmp=0.d0
2831       Gt1=0.d0
2832       do i=idxDmin,imax
2833        Gt1=Gt1+atildg(i,k)*xp**btildgp(i,k)*xm**btildgpp(i,k)
2834       enddo
2835
2836       Phitmp=exp(-Gt1)
2837
2838       PhiExpoK=Phitmp
2839
2840       return
2841       end
2842
2843 c----------------------------------------------------------------------
2844       double precision function PhiExpoK2(k,xp,xm)   !---xs---
2845 c----------------------------------------------------------------------
2846 c    Exponential expression of the Phi function for pp collision
2847 c    for given k without diffractive part
2848 c----------------------------------------------------------------------
2849       include 'epos.inc'
2850       include 'epos.incsem'
2851       include 'epos.incems'
2852
2853       double precision xp,xm,Phitmp,Gt1
2854       double precision atildg,btildgp,btildgpp
2855       common/cgtilde/atildg(idxD0:idxD1,kollmx)
2856      *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx)
2857
2858
2859       Phitmp=0.d0
2860
2861       imax=1
2862
2863       Phitmp=0.d0
2864       Gt1=0.d0
2865       do i=idxDmin,imax
2866        Gt1=Gt1+atildg(i,k)*xp**btildgp(i,k)*xm**btildgpp(i,k)
2867       enddo
2868
2869       Phitmp=exp(-Gt1)
2870
2871       PhiExpoK2=Phitmp
2872
2873       return
2874       end
2875
2876 c----------------------------------------------------------------------
2877       double precision function Phiexpo(zzip,zzit,fj,xp,xm,s,b)   !---MC---
2878 c----------------------------------------------------------------------
2879 c    Exponential expression of the Phi function for pp collision
2880 c    for given b
2881 c input :
2882 c    zzip : additionnal component for Z (nuclear effect projectile side)
2883 c    zzit : additionnal component for Z (nuclear effect target side)
2884 c    fj   : overall factor for cross section (elastic or inelastic)
2885 c    xp,xm: momentum fraction
2886 c    s    : energy square
2887 c    b    : impact parameter
2888 c----------------------------------------------------------------------
2889       include 'epos.inc'
2890       include 'epos.incsem'
2891       include 'epos.incems'
2892       include 'epos.incpar'
2893
2894       parameter(nbkbin=40)
2895       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
2896       double precision AlTi
2897       double precision BeTip,BeTipp
2898       double precision xp,xm,Phitmp,Gt1
2899
2900       imax=idxD1
2901       if(iomega.eq.2)imax=1
2902
2903       Gt1=0.d0
2904       do i=idxDmin,imax
2905         call Gfunpar(zzip,zzit,2,i,b,s,alpx,betx,betpx,epsp,epst,epss
2906      &              ,gamv)
2907         BeTip =dble(betx)
2908         BeTipp=dble(betpx)
2909         AlTi  =dble(alpx)
2910         Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp*dble(fj*xkappa**(fj-1.))
2911
2912       enddo
2913
2914       Phitmp=exp(-Gt1)
2915
2916       PhiExpo=Phitmp
2917      &     *xp**dble(alplea(iclpro))
2918      &     *xm**dble(alplea(icltar))
2919
2920       return
2921       end
2922
2923 c----------------------------------------------------------------------
2924       double precision function PhiUnit(xp,xm)   !---test---
2925 c----------------------------------------------------------------------
2926 c    Exponential expression of the Phi function for pp collision
2927 c    for given b
2928 c----------------------------------------------------------------------
2929       include 'epos.inc'
2930       include 'epos.incsem'
2931       include 'epos.incems'
2932       include 'epos.incpar'
2933
2934       double precision AlTi
2935       double precision BeTip,BeTipp
2936       double precision xp,xm,Phitmp,Gt1
2937
2938       imax=idxD1
2939       if(iomega.eq.2)imax=1
2940
2941       Gt1=0.d0
2942       do i=idxDmin,imax
2943         BeTip =betUni(i,2)
2944         BeTipp=betpUni(i,2)
2945         AlTi  =alpUni(i,2)
2946         Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp
2947 c        write(ifch,*)'Phiunit',i,xp,xm,Gt1,AlTi,BeTip,BeTipp
2948       enddo
2949
2950       Phitmp=exp(-Gt1)
2951
2952       PhiUnit=Phitmp
2953      &     *xp**dble(alplea(iclpro))
2954      &     *xm**dble(alplea(icltar))
2955
2956
2957       return
2958       end
2959
2960
2961 cc----------------------------------------------------------------------
2962 c      double precision function PhiUnit(xp,xm,s,b)   !---inu---
2963 cc----------------------------------------------------------------------
2964 c      include 'epos.inc'
2965 c      double precision xp,xm,PhiExpo,Znorm
2966 c
2967 c      PhiUnit=Phiexpo(0.,0.,1.,xp,xm,s,b)
2968 c     &          /Znorm(s,b)
2969 c
2970 c      return
2971 c      end
2972 c
2973
2974 c----------------------------------------------------------------------
2975       double precision function Hrst(s,b,xp,xm)   !test
2976 c----------------------------------------------------------------------
2977       include 'epos.inc'
2978       include 'epos.incems'
2979       include 'epos.incsem'
2980       include 'epos.incpar'
2981       parameter(idxD2=8)
2982       double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
2983       common/DGamUni/GbetUni(  idxD0:idxD2),HbetUni(  idxD0:idxD2),
2984      &               GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
2985      &               HalpUni(idxD0:idxD2)
2986       double precision al(idxD0:idxD2),betp(idxD0:idxD2)
2987      *,z,xJrst!,ffacto
2988       double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2)
2989      *,yp,ym,xp,xm
2990       dimension ipr(idxD0:idxD2),imax(idxD0:idxD2)
2991
2992       if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in Hrst"
2993       Htmp=0.d0
2994       do i=idxD0,idxD2
2995         imax(i)=0
2996         ipr(i)=0
2997         zp(i)=1.d0
2998         al(i)=0.d0
2999       enddo
3000
3001       if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then
3002
3003       imax0=idxD1
3004       if(iomega.eq.2)imax0=1
3005       imax1=idxD2
3006       if(iomega.eq.2)imax1=imax1-1
3007
3008       do i=idxDmin,imax1
3009         imax(i)=max(2,int(log10(100.*s)/3.))
3010 c        if(i.ge.2)imax(i)=imax(i)*2
3011         if(b.ge.1.5)imax(i)=2   !max(2,imax(i)/2)
3012         imax(i)=min(30,imax(i))
3013         if(i.gt.imax0)then
3014           if((zzpUni*zztUni.lt.1.d-6)
3015      &   .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
3016             imax(i)=0
3017           else
3018             imax(i)=1     !imax(i)/3
3019           endif
3020         endif
3021       enddo
3022
3023       Htmp=0.d0
3024         do i=idxDmin,imax1
3025           betp(i)=HbetUni(i)
3026           betpp(i)=HbetpUni(i)
3027           al(i)=HalpUni(i)
3028         enddo
3029
3030         do ipr0=0,imax(0)
3031 c          write(ifmt,*)'Hrst ipr0,xp,xm :',ipr0,xp,xm
3032            ipr(0)=ipr0
3033            zp(0)=1.d0
3034            if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
3035          do ipr1=0,imax(1)
3036             ipr(1)=ipr1
3037             zp(1)=1.d0
3038             if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
3039          do ipr2=0,imax(2)
3040             ipr(2)=ipr2
3041             zp(2)=1.d0
3042             if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
3043          do ipr3=0,imax(3)
3044             ipr(3)=ipr3
3045             zp(3)=1.d0
3046             if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3))
3047          do ipr4=0,imax(4)
3048             ipr(4)=ipr4
3049             zp(4)=1.d0
3050             if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4))
3051          do ipr5=0,imax(5)
3052             ipr(5)=ipr5
3053             zp(5)=1.d0
3054             if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5))
3055          do ipr6=0,imax(6)
3056             ipr(6)=ipr6
3057             zp(6)=1.d0
3058             if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6))
3059          do ipr7=0,imax(7)
3060             ipr(7)=ipr7
3061             zp(7)=1.d0
3062             if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7))
3063          do ipr8=0,imax(8)
3064             ipr(8)=ipr8
3065             zp(8)=1.d0
3066             if (ipr(8).ne.0) zp(8)=al(8)**ipr(8)*facto(ipr(8))
3067            if (ipr(0)+ipr(1)+ipr(2)+ipr(3)+ipr(4)+ipr(5)
3068      &        +ipr(6)+ipr(7)+ipr(8).ne.0) then
3069              yp=0.d0
3070              ym=0.d0
3071              z=1.d0
3072              do i=idxDmin,imax1
3073                yp=yp+dble(ipr(i))*betp(i)
3074                ym=ym+dble(ipr(i))*betpp(i)
3075                z=z*zp(i)
3076              enddo
3077              z=z*xJrst(xp,yp,GbetUni,ipr)
3078              z=z*xJrst(xm,ym,GbetpUni,ipr)
3079              Htmp=Htmp+z
3080            endif
3081           enddo
3082          enddo
3083        enddo
3084           enddo
3085          enddo
3086        enddo
3087           enddo
3088          enddo
3089        enddo
3090
3091       endif
3092
3093       Hrst=Htmp
3094
3095       return
3096       end
3097
3098 c----------------------------------------------------------------------
3099       double precision function HrstI(s,b,xp,xm)   !test
3100 c----------------------------------------------------------------------
3101       include 'epos.inc'
3102       include 'epos.incems'
3103       include 'epos.incsem'
3104       include 'epos.incpar'
3105       parameter(idxD2=8)
3106       double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
3107       common/DGamUni/GbetUni(  idxD0:idxD2),HbetUni(  idxD0:idxD2),
3108      &               GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
3109      &               HalpUni(idxD0:idxD2)
3110       double precision al(idxD0:idxD2),betp(idxD0:idxD2)
3111      *,z,xJrstI!,ffacto
3112       double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2)
3113      *,yp,ym,xp,xm
3114       dimension ipr(idxD0:idxD2),imax(idxD0:idxD2)
3115
3116       if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
3117       Htmp=0d0
3118       do i=idxD0,idxD2
3119         imax(i)=0
3120         ipr(i)=0
3121         zp(i)=1.d0
3122         al(i)=0.d0
3123       enddo
3124
3125
3126       if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then
3127
3128
3129       imax0=idxD1
3130       if(iomega.eq.2)imax0=1
3131       imax1=idxD2
3132       if(iomega.eq.2)imax1=imax1-1
3133
3134       do i=idxDmin,imax1
3135         imax(i)=max(3,int(log10(s)/2.))
3136 c        if(i.ge.2)imax(i)=imax(i)*2
3137         if(b.ge.1.5)imax(i)=max(2,imax(i)/2)
3138         imax(i)=min(30,imax(i))
3139         if(i.gt.imax0)then
3140           if((zzpUni*zztUni.lt.1.d-6)
3141      &   .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then
3142             imax(i)=0
3143           else
3144             imax(i)=1   !imax(i)/3
3145           endif
3146         endif
3147       enddo
3148
3149       Htmp=0.d0
3150         do i=idxDmin,imax1
3151           betp(i)=HbetUni(i)
3152           betpp(i)=HbetpUni(i)
3153           al(i)=HalpUni(i)
3154         enddo
3155         do ipr0=0,imax(0)
3156            ipr(0)=ipr0
3157            zp(0)=1.d0
3158            if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
3159          do ipr1=0,imax(1)
3160             ipr(1)=ipr1
3161             zp(1)=1.d0
3162             if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
3163          do ipr2=0,imax(2)
3164             ipr(2)=ipr2
3165             zp(2)=1.d0
3166             if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
3167          do ipr3=0,imax(3)
3168             ipr(3)=ipr3
3169             zp(3)=1.d0
3170             if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3))
3171          do ipr4=0,imax(4)
3172             ipr(4)=ipr4
3173             zp(4)=1.d0
3174             if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4))
3175          do ipr5=0,imax(5)
3176             ipr(5)=ipr5
3177             zp(5)=1.d0
3178             if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5))
3179          do ipr6=0,imax(6)
3180             ipr(6)=ipr6
3181             zp(6)=1.d0
3182             if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6))
3183          do ipr7=0,imax(7)
3184             ipr(7)=ipr7
3185             zp(7)=1.d0
3186             if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7))
3187          do ipr8=0,imax(8)
3188             ipr(8)=ipr8
3189             zp(8)=1.d0
3190             if (ipr(8).ne.0) zp(8)=al(8)**ipr(8)*facto(ipr(8))
3191            if (ipr(0)+ipr(1)+ipr(2)+ipr(3)+ipr(4)+ipr(5)
3192      &        +ipr(6)+ipr(7)+ipr(8).ne.0) then
3193              yp=0.d0
3194              ym=0.d0
3195              z=1.d0
3196              do i=idxDmin,imax1
3197                yp=yp+dble(ipr(i))*betp(i)
3198                ym=ym+dble(ipr(i))*betpp(i)
3199                z=z*zp(i)
3200              enddo
3201              z=z*xJrstI(xp,yp,GbetUni,ipr)
3202              z=z*xJrstI(xm,ym,GbetpUni,ipr)
3203              Htmp=Htmp+z
3204            endif
3205           enddo
3206          enddo
3207        enddo
3208           enddo
3209          enddo
3210        enddo
3211           enddo
3212          enddo
3213        enddo
3214
3215       endif
3216
3217       HrstI=Htmp
3218
3219       return
3220       end
3221
3222
3223
3224 cc----------------------------------------------------------------------
3225 c      double precision function HrstI(s,xp,xm)   !---inu---
3226 cc----------------------------------------------------------------------
3227 c      include 'epos.inc'
3228 c      include 'epos.incems'
3229 c      include 'epos.incsem'
3230 c      include 'epos.incpar'
3231 c      double precision al(idxD0:idxD1),betp(idxD0:idxD1)
3232 c     *,z,xJrstI!,ffacto
3233 c      double precision zp(idxD0:idxD1),Htmp,betpp(idxD0:idxD1)
3234 c     *,yp,ym,xp,xm
3235 c      dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
3236 c
3237 c      if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
3238 c      HrstI=0.d0
3239 c      do i=idxD0,idxD1
3240 c        imax(i)=0
3241 c        ipr(i)=0
3242 c        zp(i)=1.d0
3243 c        al(i)=0.d0
3244 c      enddo
3245 c
3246 c      if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.lt.1.d0)then
3247 c
3248 c      HrstI=0.d0
3249 c
3250 c      imax0=idxD1
3251 c      if(iomega.eq.2)imax0=1
3252 c
3253 c      do i=idxDmin,imax0
3254 c        imax(i)=max(5,int(log10(s)))
3255 cc        if(i.ge.2)imax(i)=imax(i)*2
3256 c        imax(i)=min(30,imax(i))
3257 c      enddo
3258 c      Htmp=0.d0
3259 c        do i=idxDmin,imax0
3260 c          betp(i)=betUni(i,1)+1.d0
3261 c          betpp(i)=betpUni(i,1)+1.d0
3262 c          al(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
3263 c        enddo
3264 c
3265 c        do ipr0=0,imax(0)
3266 c           ipr(0)=ipr0
3267 c           zp(0)=1.d0
3268 c          if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
3269 c         do ipr1=0,imax(1)
3270 c            ipr(1)=ipr1
3271 c            zp(1)=1.d0
3272 c          if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
3273 c         do ipr2=0,imax(2)
3274 c            ipr(2)=ipr2
3275 c            zp(2)=1.d0
3276 c          if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
3277 c             if (ipr(0)+ipr(1)+ipr(2).ne.0) then
3278 c             yp=0.d0
3279 c             ym=0.d0
3280 c             z=1.d0
3281 c             do i=idxDmin,imax0
3282 c               yp=yp+dble(ipr(i))*betp(i)
3283 c               ym=ym+dble(ipr(i))*betpp(i)
3284 c               z=z*zp(i)
3285 c             enddo
3286 c             z=z*xJrstI(xp,yp,betp,ipr)
3287 c             z=z*xJrstI(xm,ym,betpp,ipr)
3288 c             Htmp=Htmp+z
3289 c           endif
3290 c          enddo
3291 c         enddo
3292 c       enddo
3293 c
3294 c       HrstI=Htmp
3295 c
3296 c      endif
3297 c
3298 c      return
3299 c      end
3300 c
3301
3302 cc----------------------------------------------------------------------
3303 c        double precision function ffacto(n)   !---test---
3304 cc----------------------------------------------------------------------
3305 c
3306 c        ffacto=1.D0
3307 c        do i=1,n
3308 c          ffacto=ffacto*dble(i)
3309 c        enddo
3310 c        return
3311 c        end
3312 c
3313
3314 c----------------------------------------------------------------------
3315       double precision function xIrst(id,x,y,bet,ipr)   !---test---
3316 c----------------------------------------------------------------------
3317       include 'epos.inc'
3318       double precision y,gammag,utgam2,x,bet(idxD0:idxD1)
3319       dimension ipr(idxD0:idxD1)
3320
3321       if(id.eq.1)iclrem=iclpro
3322       if(id.eq.2)iclrem=icltar
3323       imax=idxD1
3324       if(iomega.eq.2)imax=1
3325       if(y.le.160.)then
3326        xIrst=gammag(iclrem,y)*x**dble(alplea(iclrem))
3327       else
3328        xIrst=0
3329       endif
3330       if(xIrst.gt.0.d0)then
3331         do i=idxDmin,imax
3332           if(ipr(i).ne.0.and.bet(i).gt.1.d-10)
3333      &         xIrst=xIrst*utgam2(bet(i))**dble(ipr(i))
3334         enddo
3335         if (abs(y).gt.1.d-10) xIrst=xIrst*x**y
3336       endif
3337       return
3338       end
3339
3340
3341 c----------------------------------------------------------------------
3342       double precision function xJrst(x,y,Gbeta,ipr)   !---inu---
3343 c----------------------------------------------------------------------
3344       include 'epos.inc'
3345       parameter(idxD2=8)
3346       double precision y,utgam2,x,Gbeta(idxD0:idxD2),eps,gam
3347       dimension ipr(idxD0:idxD2)
3348
3349       eps=1.d-10
3350
3351       imax=idxD2
3352       if(iomega.eq.2)imax=imax-1
3353
3354       gam=utgam2(y)
3355
3356       if(gam.lt.1.d99)then
3357
3358       if ((x-1.d0).gt.eps.or.(y-1.d0).gt.eps) then
3359                         xJrst=(1.d0-x)**(y-1.d0)/gam
3360                         do i=idxDmin,imax
3361       if (ipr(i).ne.0)   xJrst=xJrst*Gbeta(i)**dble(ipr(i))
3362                         enddo
3363           else
3364 c            write (*,*) 'Warning in xJrst, infinite value !'
3365                         xJrst=(1.d0-x+eps)**(y-1.d0)/gam
3366                         do i=idxDmin,imax
3367       if (ipr(i).ne.0)   xJrst=xJrst*Gbeta(i)**dble(ipr(i))
3368                         enddo
3369       endif
3370       else
3371         xJrst=0.d0
3372       endif
3373
3374       return
3375       end
3376
3377
3378 c----------------------------------------------------------------------
3379       double precision function xJrstI(x,y,Gbeta,ipr)   !---inu---
3380 c----------------------------------------------------------------------
3381 c Function used for the integration of H*Phi. We do the changement of
3382 c variable (1-x)=z**alpha. The power alpha can be change if necessary.
3383 c----------------------------------------------------------------------
3384       include 'epos.inc'
3385       parameter(idxD2=8)
3386       double precision y,utgam2,x,Gbeta(idxD0:idxD2),alpha,w,gam
3387       dimension ipr(idxD0:idxD2)
3388
3389       alpha=4.d0
3390       w=alpha*(y-1.d0)+alpha-1.d0
3391       imax=idxD2
3392       if(iomega.eq.2)imax=imax-1
3393
3394       gam=utgam2(y)
3395
3396       if(gam.lt.1.d99)then
3397
3398       if (w.ge.0)then
3399
3400                         xJrstI=alpha*x**w/gam
3401                         do i=idxDmin,imax
3402       if (ipr(i).ne.0)   xJrstI=xJrstI*Gbeta(i)**dble(ipr(i))
3403                         enddo
3404
3405          else
3406            write(*,*) 'x,y,bet,ipr,w',x,y,Gbeta,ipr,w
3407           stop 'Error in xJrstI in epos-omg, integration not possible'
3408        endif
3409
3410       else
3411         xJrstI=0.d0
3412       endif
3413
3414
3415       return
3416       end
3417
3418 c----------------------------------------------------------------------
3419       double precision function HPhiInt(s,b)   !---inu---
3420 c----------------------------------------------------------------------
3421 c  Set integrated over xp and xm (x and y) H(x,y)*Phi(x,y) for a
3422 c  given b by gauss method
3423 c----------------------------------------------------------------------
3424       include 'epos.inc'
3425       parameter(idxD2=8)
3426       double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
3427       common/DGamUni/GbetUni(  idxD0:idxD2),HbetUni(  idxD0:idxD2),
3428      &               GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
3429      &               HalpUni(idxD0:idxD2)
3430       double precision xhm,x,y,yhm,w,Hrst,utgam2,PhiUnit!,PhiExact
3431 c      double precision zp2,zm2,HrstI,eps
3432 c      common /ar3/  x1(7),a1(7)
3433       common /ar9/    x9(3),a9(3)
3434
3435       eps=0d0 !1.d-5
3436
3437       imax0=idxD1
3438       imax1=idxD2
3439       if(iomega.eq.2)then
3440         imax0=1
3441         imax1=imax1-1
3442       endif
3443       do i=idxDmin,imax0
3444         HbetUni(i)=betUni(i,1)+1.d0
3445         HbetpUni(i)=betpUni(i,1)+1.d0
3446         GbetUni(i)=utgam2(HbetUni(i))
3447         GbetpUni(i)=utgam2(HbetpUni(i))
3448         HalpUni(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
3449       enddo
3450       do i=0,1
3451         HbetUni(imax0+1+i)=betUni(i,1)+1.d0+betfom
3452         HbetUni(imax0+3+i)=betUni(i,1)+1.d0
3453         HbetUni(imax0+5+i)=betUni(i,1)+1.d0+betfom
3454         HbetpUni(imax0+1+i)=betpUni(i,1)+1.d0
3455         HbetpUni(imax0+3+i)=betpUni(i,1)+1.d0+betfom
3456         HbetpUni(imax0+5+i)=betpUni(i,1)+1.d0+betfom
3457         GbetUni(imax0+1+i)=utgam2(HbetUni(imax0+1+i))
3458         GbetUni(imax0+3+i)=utgam2(HbetUni(imax0+3+i))
3459         GbetUni(imax0+5+i)=utgam2(HbetUni(imax0+5+i))
3460         GbetpUni(imax0+1+i)=utgam2(HbetpUni(imax0+1+i))
3461         GbetpUni(imax0+3+i)=utgam2(HbetpUni(imax0+3+i))
3462         GbetpUni(imax0+5+i)=utgam2(HbetpUni(imax0+5+i))
3463         HalpUni(imax0+1+i)=zztUni*alpUni(i,1)
3464         HalpUni(imax0+3+i)=zzpUni*alpUni(i,1)
3465         HalpUni(imax0+5+i)=zzpUni*zztUni*alpUni(i,1)
3466       enddo
3467
3468       w=0.d0
3469       xhm=.5d0*(1d0-eps)
3470       yhm=.5d0*(1d0-eps)
3471       do m=1,2
3472         do i=1,3
3473 c        do i=1,7
3474           x=xhm+dble((2*m-3)*x9(i))*xhm
3475 c          write(ifmt,*)'HPhiInt, xp int :',x
3476           do n=1,2
3477             do j=1,3
3478 c            do j=1,7
3479               y=yhm+dble((2*n-3)*x9(j))*yhm
3480              w=w+dble(a9(i)*a9(j))*Hrst(s,b,x,y)
3481      &          *PhiUnit(x,y)
3482 c     &          *Phiexact(0.,0.,1.,x,y,s,b)
3483             enddo
3484           enddo
3485         enddo
3486       enddo
3487
3488       HPhiInt=w*xhm*yhm
3489
3490
3491 c      w=0.d0
3492 c      xhm=.5d0*eps
3493 c      yhm=.5d0*eps
3494 c      do m=1,2
3495 c        do i=1,7
3496 c          x=1d0-eps+xhm+dble((2*m-3)*x1(i))*xhm
3497 c          do n=1,2
3498 c            do j=1,7
3499 c              y=1d0-epsyhm+dble((2*n-3)*x1(j))*yhm
3500 c              zp2=1.d0-x**4
3501 c              zm2=1.d0-y**4
3502 c              w=w+dble(a1(i)*a1(j))*HrstI(s,x,y)
3503 cc     &             *PhiUnit(zp2,zm2)
3504 c     &             *Phiexact(0.,0.,1.,zp2,zm2,s,b)
3505 c            enddo
3506 c          enddo
3507 c        enddo
3508 c      enddo
3509 c
3510 c      HPhiInt=HPhiInt+w*xhm*yhm
3511
3512       return
3513       end
3514
3515
3516
3517 c----------------------------------------------------------------------
3518       subroutine Kfit(iiclegy)
3519 c----------------------------------------------------------------------
3520       include "epos.inc"
3521       include "epos.incsem"
3522       double precision Znorm
3523       parameter(nbkbin=40)
3524       common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin
3525       parameter (nmax=30)
3526       logical lnoch
3527
3528       if(iiclegy.eq.-1.or.iiclegy.gt.iclegy2)then
3529         do iiiegy=1,nclegy
3530         do iiipro=1,nclha
3531         do iiitar=1,nclha
3532         do iiibk=1,nbkbin
3533           xkappafit(iiiegy,iiipro,iiitar,iiibk)=1.
3534         enddo
3535         enddo
3536         enddo
3537         enddo
3538
3539       else
3540
3541       if(isetcs.le.1)then
3542         s=engy*engy
3543         eps=0.05
3544       else
3545         s=(egylow*egyfac**(iiclegy-1))**2.
3546         eps=0.001
3547       endif
3548
3549       write(ifmt,*)"Fit xkappa ..."
3550       if(ish.ge.2)then
3551         write(ifmt,*)"Kfit s,bkbin,iclegy,ipro,itar"
3552      *       ,s,bkbin,iiclegy,iclpro,icltar
3553       endif
3554
3555
3556       b=0.
3557       if(isetcs.le.1.or.iiclegy.eq.iclegy2)then
3558         xkf=1.
3559       else
3560         xkf=xkappafit(iiclegy+1,iclpro,icltar,1)
3561       endif
3562       xkfs=1.
3563       delta=0.
3564       deltas=0.
3565
3566
3567       do 5 ib=1,nbkbin-1
3568         b=float(ib-1)*bkbin
3569         xkappafit(iiclegy,iclpro,icltar,ib)=1.
3570         if(b.gt.3.+0.05*log(s).or.s.le.20.*q2min)then
3571           xkf=1.
3572           goto 5
3573         endif
3574         if(ib.gt.1.and.ish.ge.3)write(ifch,*)"    End",delta,xkf
3575         delta=1.-sngl(Znorm(s,b))
3576         if(delta.le.0d0)then
3577           if(xkf.ne.1.)then
3578             xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3579             delta=1.-sngl(Znorm(s,b))
3580           endif
3581         else!if(xkf.ne.1.)then
3582           goto 5
3583         endif
3584         if(abs(delta).lt.eps)then
3585           if(delta.lt.0d0)then
3586             xkfs=xkf-delta
3587             deltas=delta
3588           endif
3589           xkf=1.
3590           goto 5
3591         elseif(ib.le.nbkbin-1)then
3592
3593           if(delta.gt.0.d0)then
3594             xkf0=1.
3595             xkf1=xkf
3596             delta0=delta
3597             xkf2=xkf-delta0
3598             xkappafit(iiclegy,iclpro,icltar,ib)=xkf2
3599             delta1=1.-sngl(Znorm(s,b))
3600             if(delta1.lt.0.d0)then
3601               xkf0=xkf2
3602               xkf1=xkf
3603               delta=delta1
3604               xkf=xkf0
3605             else
3606               xkf1=max(delta0,xkf2)
3607               xkf0=0.
3608               xkf=xkf1
3609             endif
3610           else
3611             xkf0=xkf
3612             xkf1=1.-delta
3613             xkf2=xkf
3614             delta1=delta
3615           endif
3616
3617           if(ib.eq.1)then
3618             deltas=delta
3619             xkfs=max(0.00001,1.-delta)
3620           endif
3621
3622           if(delta.le.deltas)xkf=xkfs
3623           if(ish.ge.3)write(ifch,*)"    Start",ib,b,delta,xkf,xkf0,xkf1
3624           if(xkf.eq.xkf2)delta=delta1
3625
3626           n=0
3627           delta0=delta
3628           lnoch=.true.
3629  10       continue
3630           n=n+1
3631           if(n.le.nmax.and.xkf1.ne.xkf0)then
3632             if(abs(xkf-xkf2).gt.1e-6.or.abs(delta).gt.abs(deltas))then
3633               xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3634               delta=1.-sngl(Znorm(s,b))
3635             endif
3636             if(ish.ge.5)write(ifch,*)"    step",ib,n,delta,xkf,delta0
3637             if(delta*delta0.ge.0.)then
3638               if(lnoch.and.abs(delta).gt.abs(delta0))goto 5
3639             else
3640               lnoch=.false.
3641             endif
3642             if(abs(delta).gt.eps)then
3643               if(delta.gt.0.)then
3644                 xkf1=xkf
3645                 xkf=(xkf1+xkf0)*0.5
3646                 delta0=delta
3647               else
3648                 xkf0=xkf
3649                 xkf=(xkf1+xkf0)*0.5
3650                 delta0=delta
3651               endif
3652               goto 10
3653             endif
3654           else
3655             if(ish.ge.2)
3656      *      write(ifmt,*)"Warning in Kfit, nmax reached : xkappafit=1."
3657             xkappafit(iiclegy,iclpro,icltar,ib)=xkf
3658           endif
3659         endif
3660
3661  5    continue
3662
3663       if(ish.ge.3)write(ifch,*)"    End",delta,xkf
3664       if(xkf.gt.1.+eps)write(ifmt,*)
3665      *     "Warning in Kfit, xkappafit not yet 1"
3666       xkappafit(iiclegy,iclpro,icltar,nbkbin)=1.
3667
3668       endif
3669
3670       return
3671       end
3672
3673 c----------------------------------------------------------------------
3674       double precision function Znorm(s,b)   !---inu---
3675 c----------------------------------------------------------------------
3676       include 'epos.inc'
3677       common /kwrite/ xkapZ
3678       double precision HPhiInt,PhiUnit!,PhiExact
3679
3680 c      write(ifmt,*)'Z calculation for (s,b) :',s,b
3681       imax=idxD1
3682       if(iomega.eq.2)imax=1
3683       do i=idxDmin,imax
3684         call Gfunpar(0.,0.,1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3685         call Gfunpar(0.,0.,2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
3686       enddo
3687       call GfomPar(b,s)
3688       Znorm=HPhiInt(s,b)
3689 c      write(ifch,*)'int',Znorm,' phi',Phiexact(0.,0.,1.,1.d0,1.d0,s,b)
3690       Znorm=Znorm
3691      &       +PhiUnit(1.d0,1.d0)
3692 c     &       +Phiexact(0.,0.,1.,1.d0,1.d0,s,b)
3693
3694       !write(ifmt,*)'Z=',Znorm,xkapZ,b
3695       return
3696       end
3697
3698
3699 c------------------------------------------------------------
3700       double precision function gammag(iclrem,x)   !---test---
3701 c----------------------------------------------------------------------
3702       include 'epos.inc'
3703       include 'epos.incsem'
3704       double precision x,utgam2
3705
3706       gammag=utgam2(dble(alplea(iclrem))+1.D0)
3707      &       /utgam2(dble(alplea(iclrem))+1.D0+x)
3708
3709       return
3710       end
3711
3712
3713
3714 cc----------------------------------------------------------------------
3715 c        double precision function PomNbri(iqq)   !---xsigma---
3716 cc----------------------------------------------------------------------
3717 cc integral d2b om1intbci
3718 cc iqq, Pomeron type
3719 cc----------------------------------------------------------------------
3720 c      include 'epos.inc'
3721 c      double precision om1intbci
3722 c      common/geom/rmproj,rmtarg,bmax,bkmx
3723 c      common /ar3/  x1(7),a1(7)
3724 c
3725 c      bmid=bkmx/2.
3726 c      PomNbri=0.d0
3727 c      do i=1,7
3728 c        do m=1,2
3729 c          bb=bmid*(1.+(2.*m-3)*x1(i))
3730 c          PomNbri=PomNbri+dble(bb*a1(i))*om1intbci(bb,iqq)
3731 c        enddo
3732 c      enddo
3733 c
3734 c      PomNbri=PomNbri*dble(2.*pi*bmid)
3735 c
3736 c      return
3737 c      end
3738 c
3739 c
3740 c
3741
3742
3743 c####################################################################################
3744 c#############   former chk #########################################################
3745 c####################################################################################
3746
3747
3748
3749
3750
3751
3752
3753
3754 cc----------------------------------------------------------------------
3755 c      double precision function PomIncII(b)   !---check---
3756 cc----------------------------------------------------------------------
3757 cc  integral_dx_dy om1*F_remn*F_remn for a given b   !---check---
3758 cc----------------------------------------------------------------------
3759 c      include 'epos.inc'
3760 c      include 'epos.incems'
3761 c      include 'epos.incsem'
3762 c      include 'epos.incpar'
3763 c       double precision cint,gamom(idxD0:idxD1),deltap(idxD0:idxD1)
3764 c     &,deltapp(idxD0:idxD1),utgam2
3765 c
3766 cc Calculation by analytical integration (perfect but it changes
3767 cc if om1 change):
3768 c
3769 c      s=engy**2
3770 c      imax=1
3771 c      if(iomega.eq.2)imax=1
3772 c      do i=idxDmin,imax
3773 c        call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3774 c        gamom(i)=dble(alp*gamv)*chad(iclpro)*chad(icltar)
3775 c        deltap(i)=dble(bet)
3776 c        deltapp(i)=dble(betp)
3777 c
3778 cc Integration possible only if delta(i)>-1
3779 c
3780 c       if(deltap(i).le.-1.d0.or.deltapp(i).le.-1.d0)
3781 c     &       stop 'Error in epos-par-300 in PomIncII'
3782 c      enddo
3783 c
3784 c      cint=0.d0
3785 c      do i=idxDmin,imax
3786 c       cint=cint+gamom(i)*utgam2(deltap(i)+1.d0)*utgam2(deltapp(i)+1.d0)
3787 c     &            *dble(ucfpro*ucftar)
3788 c     &            /utgam2(dble(alplea(iclpro))+deltap(i)+2.d0)
3789 c     &            /utgam2(dble(alplea(icltar))+deltapp(i)+2.d0)
3790 c      enddo
3791 c
3792 c      PomIncII=cint
3793 c
3794 c      return
3795 c      end
3796 c
3797
3798 c----------------------------------------------------------------------
3799         double precision function PomIncXIExact(x)   !---check---
3800 c----------------------------------------------------------------------
3801 c integral d2b PomIncXExact
3802 c----------------------------------------------------------------------
3803       include 'epos.inc'
3804       double precision x,PomIncXExact
3805       common /ar3/  x1(7),a1(7)
3806       common/geom/rmproj,rmtarg,bmax,bkmx
3807
3808       bmid=bkmx/2.
3809       PomIncXIExact=0.d0
3810       do i=1,7
3811         do m=1,2
3812           bb=bmid*(1.+(2.*m-3)*x1(i))
3813           PomIncXIExact=PomIncXIExact+dble(bb*a1(i))*PomIncXExact(x,bb)
3814         enddo
3815       enddo
3816
3817       PomIncXIExact=PomIncXIExact*dble(2.*pi*bmid)
3818
3819       return
3820       end
3821
3822 c----------------------------------------------------------------------
3823         double precision function PomIncXIUnit(x)   !---check---
3824 c----------------------------------------------------------------------
3825 c integral d2b PomIncXUnit
3826 c----------------------------------------------------------------------
3827       include 'epos.inc'
3828       double precision x,PomIncXUnit
3829       common /ar3/  x1(7),a1(7)
3830       common/geom/rmproj,rmtarg,bmax,bkmx
3831
3832       bmid=bkmx/2.
3833       PomIncXIUnit=0.d0
3834       do i=1,7
3835         do m=1,2
3836           bb=bmid*(1.+(2.*m-3)*x1(i))
3837           PomIncXIUnit=PomIncXIUnit+dble(bb*a1(i))*PomIncXUnit(x,bb)
3838         enddo
3839       enddo
3840
3841       PomIncXIUnit=PomIncXIUnit*dble(2.*pi*bmid)
3842
3843       return
3844       end
3845
3846 c----------------------------------------------------------------------
3847         double precision function PomIncPIExact(x)   !---check---
3848 c----------------------------------------------------------------------
3849 c integral d2b PomIncPExact
3850 c----------------------------------------------------------------------
3851       include 'epos.inc'
3852       double precision x,PomIncPExact
3853       common/geom/rmproj,rmtarg,bmax,bkmx
3854       common /ar3/  x1(7),a1(7)
3855
3856       bmid=bkmx/2.
3857       PomIncPIExact=0.d0
3858       do i=1,7
3859         do m=1,2
3860           bb=bmid*(1.+(2.*m-3)*x1(i))
3861           PomIncPIExact=PomIncPIExact+dble(bb*a1(i))*PomIncPExact(x,bb)
3862         enddo
3863       enddo
3864
3865       PomIncPIExact=PomIncPIExact*dble(2.*pi*bmid)
3866
3867       return
3868       end
3869
3870 c----------------------------------------------------------------------
3871         double precision function PomIncPIUnit(x)   !---check---
3872 c----------------------------------------------------------------------
3873 c integral d2b PomIncPUnit
3874 c----------------------------------------------------------------------
3875       include 'epos.inc'
3876       double precision x,PomIncPUnit
3877       common/geom/rmproj,rmtarg,bmax,bkmx
3878       common /ar3/  x1(7),a1(7)
3879
3880       bmid=bkmx/2.
3881       PomIncPIUnit=0.d0
3882       do i=1,7
3883         do m=1,2
3884           bb=bmid*(1.+(2.*m-3)*x1(i))
3885           PomIncPIUnit=PomIncPIUnit+dble(bb*a1(i))*PomIncPUnit(x,bb)
3886         enddo
3887       enddo
3888
3889       PomIncPIUnit=PomIncPIUnit*dble(2.*pi*bmid)
3890
3891       return
3892       end
3893
3894 c----------------------------------------------------------------------
3895         double precision function PomIncMIExact(x)   !---check---
3896 c----------------------------------------------------------------------
3897 c integral d2b PomIncMExact
3898 c----------------------------------------------------------------------
3899       include 'epos.inc'
3900       double precision x,PomIncMExact
3901       common/geom/rmproj,rmtarg,bmax,bkmx
3902       common /ar3/  x1(7),a1(7)
3903
3904       bmid=bkmx/2.
3905       PomIncMIExact=0.d0
3906       do i=1,7
3907         do m=1,2
3908           bb=bmid*(1.+(2.*m-3)*x1(i))
3909           PomIncMIExact=PomIncMIExact+dble(bb*a1(i))*PomIncMExact(x,bb)
3910         enddo
3911       enddo
3912
3913       PomIncMIExact=PomIncMIExact*dble(2.*pi*bmid)
3914
3915       return
3916       end
3917
3918 c----------------------------------------------------------------------
3919         double precision function PomIncMIUnit(x)   !---check---
3920 c----------------------------------------------------------------------
3921 c integral d2b PomIncMUnit
3922 c----------------------------------------------------------------------
3923       include 'epos.inc'
3924       double precision x,PomIncMUnit
3925       common/geom/rmproj,rmtarg,bmax,bkmx
3926       common /ar3/  x1(7),a1(7)
3927
3928       bmid=bkmx/2.
3929       PomIncMIUnit=0.d0
3930       do i=1,7
3931         do m=1,2
3932           bb=bmid*(1.+(2.*m-3)*x1(i))
3933           PomIncMIUnit=PomIncMIUnit+dble(bb*a1(i))*PomIncMUnit(x,bb)
3934         enddo
3935       enddo
3936
3937       PomIncMIUnit=PomIncMIUnit*dble(2.*pi*bmid)
3938
3939       return
3940       end
3941
3942 c----------------------------------------------------------------------
3943         double precision function PomIncMExact(xm,b)   !---check---
3944 c----------------------------------------------------------------------
3945 c incluse Pomeron distribution \int dx+ { 2G F_remn F_remn }
3946 c----------------------------------------------------------------------
3947       include 'epos.inc'
3948       include 'epos.incsem'
3949       include 'epos.incems'
3950       double precision AlTiP,BeTip,al,bep,bepp,xpInt,utgam2,xm
3951
3952       s=engy**2
3953       PomIncMExact=0.d0
3954       imax=1
3955       if(iomega.eq.2)imax=1
3956       do i=idxDmin,imax
3957         call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3958         bep =dble(bet)
3959         bepp=dble(betp)
3960         al  =dble(alp*gamv)
3961
3962         BeTip=bep+1.d0
3963         xpInt=utgam2(BeTip)*dble(ucfpro)
3964      *                    /utgam2(1.d0+dble(alplea(iclpro))+BeTip)
3965         AlTiP=al*xpInt
3966         PomIncMExact=PomIncMExact+AlTiP*xm**bepp
3967      *                            *(1.d0-xm)**dble(alplea(icltar))
3968       enddo
3969
3970       return
3971       end
3972
3973 c----------------------------------------------------------------------
3974       double precision function PomIncMUnit(xm,b)   !---check---
3975 c----------------------------------------------------------------------
3976 c incluse  Unitarized Pomeron distribution  \int dx+
3977 c----------------------------------------------------------------------
3978       include 'epos.inc'
3979       include 'epos.incsem'
3980       include 'epos.incems'
3981       double precision Df,xp,xm,G2,w,xpm
3982       double precision PoInU!,Znorm
3983       common /ar3/  x1(7),a1(7)
3984
3985       s=engy**2
3986
3987 c Calculation by numeric integration :
3988       w=0.d0
3989       xpm=.5d0
3990       do m=1,2
3991         do j=1,7
3992           xp=xpm*(1.d0+dble((2.*m-3.)*x1(j)))
3993           Df=0.d0
3994           do i=idxDmin,idxD1
3995             call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3996             Df=Df+dble(alp)*xp**dble(bet)*xm**dble(betp)
3997           enddo
3998           call GfomPar(b,s)
3999           G2=Df*(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
4000           w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2
4001         enddo
4002       enddo
4003       w=w*xpm
4004
4005
4006       PomIncMUnit=w!/Znorm(s,b)
4007
4008       return
4009       end
4010
4011
4012 c----------------------------------------------------------------------
4013       double precision function PomIncPExact(xp,b)   !---check---
4014 c----------------------------------------------------------------------
4015 c incluse Pomeron distribution  \int dx- { 2G F_remn F_remn }
4016 c----------------------------------------------------------------------
4017       include 'epos.inc'
4018       include 'epos.incsem'
4019       include 'epos.incems'
4020       double precision AlTiP,BeTipp,al,bep,bepp,xmInt,utgam2,xp
4021
4022       s=engy**2
4023       PomIncPExact=0.d0
4024       imax=1
4025       if(iomega.eq.2)imax=1
4026       do i=idxDmin,imax
4027         call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
4028         bep=dble(bet)
4029         bepp=dble(betp)
4030         al=dble(alp*gamv)
4031         BeTipp=bepp+1.d0
4032         xmInt=utgam2(BeTipp)*dble(ucftar)
4033      *                    /utgam2(1.d0+dble(alplea(icltar))+BeTipp)
4034         AlTiP=al*xmInt
4035         PomIncPExact=PomIncPExact+AlTiP*xp**bep
4036      *                            *(1.d0-xp)**dble(alplea(iclpro))
4037       enddo
4038
4039       return
4040       end
4041
4042 c----------------------------------------------------------------------
4043       double precision function PomIncPUnit(xp,b)   !---check---
4044 c----------------------------------------------------------------------
4045 c incluse  Unitarized Pomeron distribution  \int dx-
4046 c----------------------------------------------------------------------
4047       include 'epos.inc'
4048       include 'epos.incsem'
4049       double precision Df,xp,xm,G2,w,xmm
4050       double precision PoInU!,Znorm
4051       common /ar3/  x1(7),a1(7)
4052
4053       s=engy**2
4054       imax=1
4055       if(iomega.eq.2)imax=1
4056
4057 c Calculation by numeric integration :
4058       w=0.d0
4059       xmm=.5d0
4060       do m=1,2
4061         do j=1,7
4062           xm=xmm*(1.d0+dble((2.*m-3.)*x1(j)))
4063           Df=0.d0
4064           do i=idxDmin,imax
4065             call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
4066             Df=Df+alp*xp**dble(bet)*xm**dble(betp)
4067           enddo
4068           call GfomPar(b,s)
4069           G2=Df*(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
4070           w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2
4071         enddo
4072       enddo
4073       w=w*xmm
4074
4075
4076       PomIncPUnit=w!/Znorm(s,b)
4077
4078       return
4079       end
4080
4081
4082 c----------------------------------------------------------------------
4083         double precision function PomIncJExact(b)   !---check---
4084 c----------------------------------------------------------------------
4085 c integral of Pomeron distribution  \int dy dx { 2G F_remn F_remn }
4086 c----------------------------------------------------------------------
4087       include 'epos.inc'
4088       double precision allea,PomIncXExact,xh
4089       common /ar3/  x1(7),a1(7)
4090
4091       allea=2.d0+dble(alplea(iclpro)+alplea(icltar))
4092       PomIncJExact=0.d0
4093       do i=1,7
4094         do m=1,2
4095           xh=1.d0-(.5d0+dble(x1(i)*(float(m)-1.5)))**(1.d0/allea)
4096           PomIncJExact=PomIncJExact+dble(a1(i))
4097      &       *PomIncXExact(xh,b)/(1.d0-xh)**(allea-1.d0)
4098         enddo
4099       enddo
4100       PomIncJExact=PomIncJExact/allea/2.d0
4101
4102       return
4103       end
4104
4105
4106 c----------------------------------------------------------------------
4107         double precision function PomIncExactk(k)   !---MC---
4108 c----------------------------------------------------------------------
4109 c integral of Pomeron distribution  \int dy dx { 2G F_remn F_remn }
4110 c----------------------------------------------------------------------
4111       include 'epos.inc'
4112       include 'epos.incems'
4113       double precision allea,PomIncXExactk,xh
4114       common /ar3/  x1(7),a1(7)
4115
4116       allea=2.d0+dble(alplea(iclpro)+alplea(icltar))
4117       PomIncExactk=0.d0
4118       do i=1,7
4119         do m=1,2
4120           xh=1.d0-(.5d0+dble(x1(i)*(float(m)-1.5)))**(1.d0/allea)
4121           PomIncExactk=PomIncExactk+dble(a1(i))
4122      &       *PomIncXExactk(k,xh)/(1.d0-xh)**(allea-1.d0)
4123         enddo
4124       enddo
4125       PomIncExactk=PomIncExactk/allea/2.d0
4126
4127       return
4128       end
4129
4130
4131 c----------------------------------------------------------------------
4132         double precision function PomIncJUnit(b)   !---check---
4133 c----------------------------------------------------------------------
4134 c integral of Pomeron distribution  \int dy dx { 2G F_remn F_remn }
4135 c----------------------------------------------------------------------
4136       include 'epos.inc'
4137       double precision PomIncXUnit,xh,xhm
4138       common /ar3/  x1(7),a1(7)
4139
4140       PomIncJUnit=0.d0
4141       xhm=.5d0
4142       do i=1,7
4143         do m=1,2
4144           xh=xhm*(1.d0+dble(x1(i)*(2.*float(m)-3.)))
4145           PomIncJUnit=PomIncJUnit+dble(a1(i))
4146      &                                *PomIncXUnit(xh,b)
4147         enddo
4148       enddo
4149       PomIncJUnit=PomIncJUnit*xhm
4150
4151       return
4152       end
4153
4154
4155 c----------------------------------------------------------------------
4156       double precision function PomIncXExactk(k,xh)   !---MC---
4157 c----------------------------------------------------------------------
4158 c incluse Pomeron distribution  \int dy { 2G F_remn F_remn }
4159 c----------------------------------------------------------------------
4160       include 'epos.inc'
4161       include 'epos.incsem'
4162       include 'epos.incems'
4163       double precision xh,Df,xp,xm,w,ymax,bet,betp,alp
4164       common /ar3/  x1(7),a1(7)
4165
4166       imax=1
4167       if(iomega.eq.2)imax=1
4168 c Calculation by numeric integration :
4169       w=0.d0
4170       ymax=-.5d0*log(xh)
4171       do m=1,2
4172         do j=1,7
4173           xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
4174           xm=xh/xp
4175           Df=0.d0
4176           do i=idxDmin,imax
4177             bet=btildep(i,k)+dble(alppar)
4178             betp=btildepp(i,k)+dble(alppar)
4179             alp=atilde(i,k)
4180             Df=Df+alp*xp**bet*xm**betp
4181      *      *(1.d0-xp)**dble(alplea(iclpro))
4182      *      *(1.d0-xm)**dble(alplea(icltar))
4183           enddo
4184           w=w+dble(a1(j))*Df
4185         enddo
4186       enddo
4187       w=w*ymax*xh**dble(-alppar)
4188
4189
4190       PomIncXExactk=w
4191
4192       return
4193       end
4194
4195 c----------------------------------------------------------------------
4196         double precision function PomIncXExact(xh,b)   !---check---
4197 c----------------------------------------------------------------------
4198 c incluse Pomeron distribution  \int dy { 2G F_remn F_remn }
4199 c (optimized integration but with no y dependance)
4200 c----------------------------------------------------------------------
4201       include 'epos.inc'
4202       include 'epos.incsem'
4203       double precision AlTiP,bep,bepp,factor,factor1
4204       double precision xpmin,xh,xp,xm,ymax,y
4205       common /ar3/  x1(7),a1(7)
4206
4207       imax=1
4208       if(iomega.eq.2)imax=1
4209       s=engy**2
4210       PomIncXExact=0.d0
4211       do i=idxDmin,imax
4212         call Gfunpar(0.,0.,1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
4213         bep  =betx
4214         bepp =betpx
4215         AlTiP=alpx*gamv
4216         PomIncXExact=PomIncXExact+AlTiP*xh**((bep+bepp)/2.d0)
4217       enddo
4218
4219       factor=0.d0
4220       allea=min(alplea(iclpro),alplea(icltar))+1.
4221       xpmin=max(sqrt(xh),exp(-1.d0))
4222       do i=1,7
4223         do m=1,2
4224           xp=1.d0-(1.d0-xpmin)*(.5d0+dble(x1(i)*(float(m)-1.5)))
4225      *                                         **(1.d0/dble(allea))
4226           xm=xh/xp
4227           factor=factor+dble(a1(i))
4228      *        *((1.d0-xp)**dble(alplea(iclpro)-allea+1.)
4229      *        *(1.d0-xm)**dble(alplea(icltar))
4230      *        +(1.d0-xp)**dble(alplea(icltar)-allea+1.)
4231      *        *(1.d0-xm)**dble(alplea(iclpro)))/xp
4232         enddo
4233       enddo
4234       factor=factor*(1.d0-xpmin)**dble(allea)/dble(allea)
4235
4236
4237       if(xpmin.gt.1.00001d0*sqrt(xh))then
4238         ymax=-log(xh)-2.d0
4239         factor1=0.d0
4240         do i=1,7
4241           do m=1,2
4242             y=ymax*dble(x1(i)*(2*m-3))
4243             xp=sqrt(xh*exp(y))
4244             xm=xh/xp
4245             factor1=factor1+dble(a1(i))*(1.d0-xp)**dble(alplea(iclpro))
4246      *                                 *(1.d0-xm)**dble(alplea(icltar))
4247           enddo
4248         enddo
4249         factor=factor+factor1*ymax
4250       endif
4251
4252       factor=factor/2.d0
4253
4254       PomIncXExact=PomIncXExact*factor
4255
4256       return
4257       end
4258
4259
4260 c----------------------------------------------------------------------
4261       double precision function PomIncXUnit(xh,b)   !---check---
4262 c----------------------------------------------------------------------
4263 c incluse  Unitarized Pomeron distribution  \int dy
4264 c----------------------------------------------------------------------
4265       include 'epos.inc'
4266       include 'epos.incsem'
4267       double precision xh,Df,xp,xm,w
4268       double precision PoInU,ymax!,Znorm
4269       common /ar3/  x1(7),a1(7)
4270
4271       imax=1
4272       if(iomega.eq.2)imax=1
4273       s=engy**2
4274 ctp060829      sy=s*sngl(xh)
4275 c Calculation by numeric integration :
4276       w=0.d0
4277       ymax=-.5d0*log(xh)
4278       do m=1,2
4279         do j=1,7
4280           xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
4281           xm=xh/xp
4282           Df=0.d0
4283           do i=idxDmin,imax
4284             call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
4285             Df=Df+alp*xp**dble(bet)*xm**dble(betp)
4286           enddo
4287           call GfomPar(b,s)
4288           Df=Df*(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
4289           w=w+dble(a1(j))*Df*PoInU(xp,xm,s,b)
4290         enddo
4291       enddo
4292       w=w*ymax
4293
4294
4295       PomIncXUnit=w!/Znorm(s,b)
4296
4297       return
4298       end
4299
4300
4301
4302 c----------------------------------------------------------------------
4303       double precision function PoInU(xp,xm,s,b)   !---check---
4304 c----------------------------------------------------------------------
4305 c Function : PhiU(1-xp,1-xm) + /int(H(z+ + x+,z- + x-)PhiU(z+,z-)dz+dz-)
4306 c----------------------------------------------------------------------
4307       include 'epos.inc'
4308       double precision xp,xm,zp,zm,zp2,zm2,zpmin,zmmin,deltzp,deltzm
4309       double precision zpm,zmm,w,HrstI,PhiUnit,Hrst,eps!,PhiExact
4310       common /ar3/  x1(7),a1(7)
4311
4312       eps=1.d-5
4313
4314       imax=idxD1
4315       if(iomega.eq.2)imax=1
4316       do i=idxDmin,imax
4317         call Gfunpar(0.,0.,1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
4318         call Gfunpar(0.,0.,2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv)
4319       enddo
4320       call GfomPar(b,s)
4321
4322       if (1.d0-xp-eps.gt.0.d0.and.1.d0-xm-eps.gt.0.d0) then
4323         w=0.d0
4324         zpmin=1.d0-xp-eps
4325         zmmin=1.d0-xm-eps
4326         zpm=.5d0*zpmin
4327         zmm=.5d0*zmmin
4328         do m=1,2
4329           do i=1,7
4330             zp=zpm+dble((2*m-3)*x1(i))*zpm
4331             do n=1,2
4332               do j=1,7
4333                 zm=zmm+dble((2*n-3)*x1(j))*zmm
4334             w=w+dble(a1(i)*a1(j))*Hrst(s,b,zp+xp,zm+xm)
4335      &                              *PhiUnit(zp,zm)
4336 c     &                             *Phiexact(0.,0.,1.,zp,zm,s,b)
4337              enddo
4338            enddo
4339          enddo
4340        enddo
4341        PoInU=w*zpm*zmm
4342        deltzp=eps
4343        deltzm=eps
4344       else
4345         PoInU=0.d0
4346         zpmin=0.d0
4347         zmmin=0.d0
4348         deltzp=1.d0-xp
4349         deltzm=1.d0-xm
4350       endif
4351
4352       w=0.d0
4353       zpm=0d0
4354       zmm=0d0
4355       if(abs(deltzp).gt.1.d-10.and.abs(deltzm).gt.1.d-10)then
4356       zpm=.5d0*deltzp
4357       zmm=.5d0*deltzm
4358       do m=1,2
4359         do i=1,7
4360           zp=zpmin+zpm+dble((2*m-3)*x1(i))*zpm
4361           do n=1,2
4362             do j=1,7
4363               zm=zmmin+zmm+dble((2*n-3)*x1(j))*zmm
4364               zp2=1.d0-xp-zp**2
4365               zm2=1.d0-xm-zm**2
4366               w=w+dble(a1(i)*a1(j))*HrstI(s,b,zp,zm)
4367      &             *PhiUnit(zp2,zm2)
4368 c     &             *Phiexact(0.,0.,1.,zp2,zm2,s,b)
4369             enddo
4370           enddo
4371         enddo
4372       enddo
4373       endif
4374
4375       PoInU=PoInU+w*zpm*zmm
4376      &           +PhiUnit(1.d0-xp,1.d0-xm)
4377 c     &           +Phiexact(0.,0.,1.,1.d0-xp,1.d0-xm,s,b)
4378
4379       return
4380       end
4381
4382 c----------------------------------------------------------------------
4383         double precision function PomIncExact(xp,xm,b)   !---check---
4384 c----------------------------------------------------------------------
4385 c inclusive Pomeron distribution  { 2G F_remn F_remn }
4386 c----------------------------------------------------------------------
4387       include "epos.inc"
4388       include "epos.incsem"
4389       double precision Df,xp,xm
4390
4391       Df=0.d0
4392       s=engy**2
4393       imax=1
4394       if(iomega.eq.2)imax=1
4395       do i=idxDmin,imax
4396         call Gfunpar(0.,0.,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
4397         Df=Df+alp*gamv*xp**dble(bet)*xm**dble(betp)
4398       enddo
4399       Df=dble(chad(iclpro)*chad(icltar))
4400      *  *Df
4401
4402       PomIncExact=Df
4403      *            *(1.d0-xp)**dble(alplea(iclpro))
4404      *            *(1.d0-xm)**dble(alplea(icltar))
4405
4406       return
4407       end
4408
4409 c----------------------------------------------------------------------
4410         double precision function PomIncUnit(xp,xm,b)   !---check---
4411 c----------------------------------------------------------------------
4412 c inclusive Pomeron distribution  { Sum{int{G*Phi} }
4413 c----------------------------------------------------------------------
4414       include "epos.inc"
4415       include "epos.incpar"
4416       include "epos.incsem"
4417       double precision PoInU,xp,xm,om1,xh,yp
4418
4419
4420       xh=xp*xm
4421       yp=0.d0
4422       if(xm.ne.0.d0)yp=0.5d0*log(xp/xm)
4423       PomIncUnit=om1(xh,yp,b)*PoInU(xp,xm,engy*engy,b)
4424      &          *(1.d0+zztUni*xp**betfom)*(1.d0+zzpUni*xm**betfom)
4425
4426       return
4427       end
4428
4429
4430 cc----------------------------------------------------------------------
4431 c        double precision function PomIncUnitMC(xp,xm,b)   !---check---
4432 cc----------------------------------------------------------------------
4433 cc inclusive Pomeron distribution  { Sum{int{G*Phi} }
4434 cc----------------------------------------------------------------------
4435 c      include "epos.inc"
4436 c      include "epos.incpar"
4437 c      include "epos.incsem"
4438 c      include "epos.incems"
4439 c      parameter(mmax=20)
4440 c      double precision Gtab,Phitab,xxpu(mmax),xxmu(mmax)
4441 c      double precision Zm,xp,xm,pGtab,Z,omNpcut,xprem,xmrem,
4442 c     *                 sxp,sxm,PhiExpo
4443 c
4444 c      PomIncUnitMC=0.d0
4445 c      if(xp.lt.1.d0.and.xm.lt.1.d0)then
4446 c      m=10
4447 c
4448 c      sy=engy*engy
4449 c      nmcint=2000
4450 c      nmax=nmcint
4451 c      do i=1,mmax
4452 c        xxpu(i)=0.d0
4453 c        xxmu(i)=0.d0
4454 c      enddo
4455 c      xprem=1.d0
4456 c      xmrem=1.d0
4457 c      sxp=xprem-xp
4458 c      sxm=xmrem-xm
4459 c
4460 c      Gtab=omNpcut(xp,xm,sxp,sxm,b,0)
4461 c      Phitab=Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
4462 c      Z=Gtab*Phitab
4463 c      Zm=0.d0
4464 c
4465 c      do mtmp=2,m
4466 c
4467 c        write(*,*)"GPhi",mtmp-1,Zm,Z
4468 c        Zm=0.d0
4469 c        n=0
4470 c
4471 c 10     continue
4472 c          n=n+1
4473 c          if(mod(n,1000000).eq.0)write(*,*)
4474 c     &              "Calculation of PomIncUnit(",mtmp,")->",n
4475 c          xxpu(1)=xp
4476 c          xxmu(1)=xm
4477 c          sxp=xxpu(1)
4478 c          sxm=xxmu(1)
4479 c          pGtab=1.d0
4480 c          do i=2,mtmp
4481 c            rnau=rangen()*sngl(xprem-sxp)
4482 c            xxpu(i)=dble(rnau)
4483 c            sxp=sxp+xxpu(i)
4484 c            rnau=rangen()*sngl(xmrem-sxm)
4485 c            xxmu(i)=dble(rnau)
4486 c            sxm=sxm+xxmu(i)
4487 c          enddo
4488 c          if(sxp.lt.xprem.and.sxm.lt.xmrem)then
4489 c            do i=1,mtmp
4490 c              Gtab=omNpcut(xxpu(i),xxmu(i),xprem-sxp,xmrem-sxm,b,0)
4491 c              pGtab=pGtab*Gtab
4492 c            enddo
4493 c          Zm=Zm+pGtab*Phiexpo(0.,0.,1.,xprem-sxp,xmrem-sxm,sy,b)
4494 c          endif
4495 c        if(n.lt.nmax)goto 10
4496 c        Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(mtmp)
4497 c        Z=Z+Zm
4498 c      enddo
4499 c
4500 c      PomIncUnitMC=Z/dble(chad(iclpro)*chad(icltar))
4501 c      endif
4502 c
4503 c      return
4504 c      end
4505 c
4506 c
4507 cc----------------------------------------------------------------------
4508 c        double precision function PhiMCExact(xp,xm,b)   !---check---
4509 cc----------------------------------------------------------------------
4510 cc virtual emissions  { Sum{int{-GFF} }
4511 cc----------------------------------------------------------------------
4512 c      include "epos.inc"
4513 c      include "epos.incpar"
4514 c      include "epos.incsem"
4515 c      include "epos.incems"
4516 c      parameter(mmax=20)
4517 c      double precision Gtab,xxpu(mmax),xxmu(mmax)
4518 c      double precision Zm,xp,xm,pGtab,Z,om51,sxp,sxm,xh,yp
4519 cc     *                 ,omNpuncut
4520 c
4521 c      PhiMCExact=0.d0
4522 c      if(xp.le.1.d0.and.xm.le.1.d0)then
4523 c      m=6
4524 c
4525 c      sy=engy*engy
4526 c      nmcint=50000
4527 c      nmax=nmcint
4528 c      do i=1,mmax
4529 c        xxpu(i)=0.d0
4530 c        xxmu(i)=0.d0
4531 c      enddo
4532 c
4533 c      Z=xp**dble(alplea(iclpro))
4534 c     * *xm**dble(alplea(icltar))
4535 c      Zm=0.d0
4536 c
4537 c      do mtmp=1,m
4538 c
4539 c        write(*,*)"GPhi",mtmp-1,Zm,Z/xp**dble(alplea(iclpro))
4540 c     *                              /xm**dble(alplea(icltar))
4541 c        Zm=0.d0
4542 c        n=0
4543 c
4544 c 10     continue
4545 c          n=n+1
4546 c          if(mod(n,1000000).eq.0)write(*,*)
4547 c     &              "Calculation of Phiexact(0.,0.,",mtmp,")->",n
4548 c          sxp=0.d0
4549 c          sxm=0.d0
4550 c          pGtab=1.d0
4551 c          do i=1,mtmp
4552 c            rnau=rangen()!*sngl(xp-sxp)
4553 c            xxpu(i)=dble(rnau)
4554 c            sxp=sxp+xxpu(i)
4555 c            rnau=rangen()!*sngl(xm-sxm)
4556 c            xxmu(i)=dble(rnau)
4557 c            sxm=sxm+xxmu(i)
4558 c          enddo
4559 c          if(sxp.lt.xp.and.sxm.lt.xm)then
4560 c            do i=1,mtmp
4561 c              xh=xxpu(i)*xxmu(i)
4562 c              if(abs(xh).gt.1.d-30)then
4563 c                yp=0.5d0*log(xxpu(i)/xxmu(i))
4564 c              else
4565 c                yp=0.d0
4566 c              endif
4567 c              Gtab=2*om51(xh,yp,b,0,4)
4568 cc     *            +omNpuncut(sy*sngl(xh),xh,yp,b,1) !om1(xh,yp,b)
4569 c              pGtab=pGtab*(-Gtab)
4570 c            enddo
4571 c            Zm=Zm+pGtab*(xp-sxp)**dble(alplea(iclpro))
4572 c     *           *(xm-sxm)**dble(alplea(icltar))
4573 c          endif
4574 c        if(n.lt.nmax)goto 10
4575 c        Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(m)
4576 c        Z=Z+Zm
4577 c      enddo
4578 c
4579 c      PhiMCExact=Z
4580 c      endif
4581 c
4582 c      return
4583 c      end
4584
4585 c----------------------------------------------------------------------
4586         double precision function Gammapp(sy,b,mtmp)   !---check---
4587 c----------------------------------------------------------------------
4588       include "epos.inc"
4589       include "epos.incpar"
4590       include "epos.incsem"
4591       include "epos.incems"
4592       parameter(mmax=20)
4593       double precision Gtab,xxpu(mmax),xxmu(mmax),PhiExpo
4594       double precision Zm,xp,xm,pGtab,om1,sxp,sxm,xh,yp
4595
4596       Gammapp=0.d0
4597
4598       xp=1.d0
4599       xm=1.d0
4600       nmcint=20000
4601       nmax=nmcint
4602       do i=1,mmax
4603         xxpu(i)=0.d0
4604         xxmu(i)=0.d0
4605       enddo
4606       Zm=0.d0
4607
4608         n=0
4609
4610  10     continue
4611           n=n+1
4612           if(mod(n,10000).eq.0)write(*,*)
4613      &              "Calculation of Gammapp(",mtmp,")->",n
4614           sxp=0.d0
4615           sxm=0.d0
4616           pGtab=1.d0
4617           do i=1,mtmp
4618             rnau=rangen()!*sngl(xp-sxp)
4619             xxpu(i)=dble(rnau)
4620             sxp=sxp+xxpu(i)
4621             rnau=rangen()!*sngl(xm-sxm)
4622             xxmu(i)=dble(rnau)
4623             sxm=sxm+xxmu(i)
4624           enddo
4625           if(sxp.lt.xp.and.sxm.lt.xm)then
4626             do i=1,mtmp
4627               xh=xxpu(i)*xxmu(i)
4628               if(abs(xh).gt.1.d-30)then
4629                 yp=0.5d0*log(xxpu(i)/xxmu(i))
4630               else
4631                 yp=0.d0
4632               endif
4633               Gtab=om1(xh,yp,b)
4634               pGtab=pGtab*Gtab
4635             enddo
4636             Zm=Zm+pGtab*Phiexpo(0.,0.,1.,xp-sxp,xm-sxm,sy,b)
4637           endif
4638         if(n.lt.nmax)goto 10
4639         Zm=Zm/dble(nmax)!**2.d0*(xp*xm)**dble(mtmp)
4640
4641       Gammapp=Zm
4642
4643       return
4644       end
4645
4646 cc----------------------------------------------------------------------
4647 c        double precision function GammaMCnew(sy,b,mtmp)   !---check---
4648 cc----------------------------------------------------------------------
4649 c      include "epos.inc"
4650 c      include "epos.incpar"
4651 c      include "epos.incsem"
4652 c      include "epos.incems"
4653 c      parameter(mmax=20)
4654 c      common /psar7/ delx,alam3p,gam3p
4655 c      double precision Gtab,xxpu(mmax),xxmu(mmax)
4656 c      double precision Zm,xp,xm,pGtab,omGam,Zmtot,
4657 c     *                 sxp,sxm,PhiExpo!,om1,yp,xh
4658 c
4659 c      GammaMCnew=0.d0
4660 c      Zmtot=0.d0
4661 c      xp=1.d0
4662 c      xm=1.d0
4663 c      nmcint=1000
4664 c      nmax=nmcint
4665 c
4666 c
4667 c      do i=1,mmax
4668 c        xxpu(i)=0.d0
4669 c        xxmu(i)=0.d0
4670 c      enddo
4671 c
4672 c      Zm=0.d0
4673 c
4674 c        n=0
4675 c
4676 c 10     continue
4677 c          n=n+1
4678 c          if(mod(n,1000000).eq.0)write(*,*)
4679 c     &              "Calculation of GammaMCnew(",mtmp,")->",n
4680 c          sxp=0.d0
4681 c          sxm=0.d0
4682 c          pGtab=1.d0
4683 c          do i=1,mtmp
4684 c            rnau=rangen()
4685 c            xxpu(i)=dble(rnau)
4686 c            sxp=sxp+xxpu(i)
4687 c            rnau=rangen()
4688 c            xxmu(i)=dble(rnau)
4689 c            sxm=sxm+xxmu(i)
4690 c          enddo
4691 c          if(sxp.lt.xp.and.sxm.lt.xm)then
4692 c            i=0
4693 c            do k=1,mtmp
4694 c                i=i+1
4695 c                Gtab=omGam(xxpu(i),xxmu(i),b) !om1(xh,yp,b)
4696 c                pGtab=pGtab*Gtab
4697 c              enddo
4698 c            Zm=Zm+pGtab*Phiexpo(0.,0.,1.,xp-sxp,xm-sxm,sy,b)
4699 c          endif
4700 c        if(n.lt.nmax)goto 10
4701 c
4702 c          Zmtot=Zmtot+Zm/dble(nmax)
4703 c
4704 c      GammaMCnew=Zmtot
4705 c
4706 c      return
4707 c      end
4708 cc----------------------------------------------------------------------
4709 c      double precision function GammaMC(sy,b,mtmp)   !---check---
4710 cc----------------------------------------------------------------------
4711 c      include "epos.inc"
4712 c      include "epos.incpar"
4713 c      include "epos.incsem"
4714 c      include "epos.incems"
4715 c      parameter(mmax=20)
4716 c      double precision Gtab,xxpu(mmax),xxmu(mmax)
4717 c      double precision Zm,xp,xm,pGtab,om1,
4718 c     *                 sxp,sxm,xh,yp,PhiExpo!,omNpcut
4719 c
4720 c      GammaMC=0.d0
4721 c
4722 c      xp=1.d0
4723 c      xm=1.d0
4724 c      nmcint=50000
4725 c      nmax=nmcint
4726 c      do i=1,mmax
4727 c        xxpu(i)=0.d0
4728 c        xxmu(i)=0.d0
4729 c      enddo
4730 c
4731 c      Zm=0.d0
4732 c
4733 c        n=0
4734 c
4735 c 10     continue
4736 c          n=n+1
4737 c          if(mod(n,1000000).eq.0)write(*,*)
4738 c     &              "Calculation of GammaMC(",mtmp,")->",n
4739 c          sxp=0.d0
4740 c          sxm=0.d0
4741 c          pGtab=1.d0
4742 c          do i=1,mtmp
4743 c            rnau=rangen()!*sngl(xp-sxp)
4744 c            xxpu(i)=dble(rnau)
4745 c            sxp=sxp+xxpu(i)
4746 c            rnau=rangen()!*sngl(xm-sxm)
4747 c            xxmu(i)=dble(rnau)
4748 c            sxm=sxm+xxmu(i)
4749 c          enddo
4750 c          if(sxp.lt.xp.and.sxm.lt.xm)then
4751 c            do i=1,mtmp
4752 c              xh=xxpu(i)*xxmu(i)
4753 c              if(abs(xh).gt.1.d-30)then
4754 c                yp=0.5d0*log(xxpu(i)/xxmu(i))
4755 c              else
4756 c                yp=0.d0
4757 c              endif
4758 c              Gtab=om1(xh,yp,b)!omNpcut(xxpu(i),xxmu(i),xp-sxp,xm-sxm,b,0) !om1(xh,yp,b)
4759 c              pGtab=pGtab*Gtab
4760 c            enddo
4761 c            Zm=Zm+pGtab*Phiexpo(0.,0.,1.,xp-sxp,xm-sxm,sy,b)
4762 c          endif
4763 c        if(n.lt.nmax)goto 10
4764 c        Zm=Zm/dble(nmax)*fctrl(n-mtmp)*facto(mtmp)
4765 c
4766 c      GammaMC=Zm
4767 c
4768 c      return
4769 c      end
4770 c
4771 c
4772 c
4773
4774 c----------------------------------------------------------------------
4775         double precision function GammaGauss(sy,b,mtmp)   !---check---
4776 c----------------------------------------------------------------------
4777       include "epos.inc"
4778       include "epos.incpar"
4779       include "epos.incsem"
4780       include "epos.incems"
4781       parameter(mmax=3)
4782       common /psar7/ delx,alam3p,gam3p
4783       double precision xpmin,xmmin,Gst,zm(mmax),zp(mmax)
4784      *,xpmax,xmmax,zpmin(mmax),zmmin(mmax),zpmax(mmax)
4785       double precision xp,xm,pGtab,omGam,dzp(mmax),Gp1,Gm1,xmin,eps
4786      *,sxp,sxm,PhiExpo,zmmax(mmax),dzm(mmax),Gp2,Gm2,Gp3,Gm3,G0
4787 c     *,PhiExact
4788       common /ar3/  x1(7),a1(7)
4789 c      double precision dgx1,dga1
4790 c      common /dga20/ dgx1(10),dga1(10)
4791
4792       GammaGauss=0.d0
4793       xp=1.d0
4794       xm=1.d0
4795       xmin=1.d-13
4796       eps=1.d-15
4797
4798       if(mtmp.eq.0)then
4799         nmax1=0
4800         jmax1=0
4801         nmax2=0
4802         jmax2=0
4803         nmax3=0
4804         jmax3=0
4805       elseif(mtmp.eq.1)then
4806         nmax1=2
4807         jmax1=7
4808         nmax2=0
4809         jmax2=0
4810         nmax3=0
4811         jmax3=0
4812       elseif(mtmp.eq.2)then
4813         nmax1=2
4814         jmax1=7
4815         nmax2=2
4816         jmax2=7
4817         nmax3=0
4818         jmax3=0
4819       elseif(mtmp.eq.3)then
4820         nmax1=2
4821         jmax1=7
4822         nmax2=2
4823         jmax2=7
4824         nmax3=2
4825         jmax3=7
4826       else
4827         write(*,*)"m not between 0 and 3, return ..."
4828         return
4829       endif
4830
4831         xpmin=xmin
4832         xmmin=xmin
4833         xpmax=1.d0
4834         xmmax=1.d0
4835       do i=1,mmax
4836         zp(i)=0.d0
4837         zm(i)=0.d0
4838         dzp(i)=0.d0
4839         dzm(i)=0.d0
4840         zmmin(i)=0.d0
4841         zpmax(i)=0.d0
4842         zpmin(i)=0.d0
4843         zmmax(i)=0.d0
4844       enddo
4845
4846         G0=1.d0
4847
4848         if(mtmp.eq.0)then
4849           sxp=xp
4850           sxm=xm
4851           G0=Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
4852         endif
4853
4854
4855 c        write(*,*)'x+/-',xmmin,xmmax,xpmin,xpmax
4856
4857
4858         dzm(1)=0.d0
4859         if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.1)then
4860         zmmax(1)=-log(xmmin)
4861         zmmin(1)=-log(xmmax)
4862         if(abs(xmmin-xmin).lt.eps)then
4863           zmmin(1)=-log(min(xmmax,1.d0-xmmin-xmmin))
4864           zmmax(1)=-log(max(xmmin,1.d0-xmmax-xmmax))
4865         endif
4866         dzm(1)=(zmmax(1)-zmmin(1))/2.d0
4867         endif
4868
4869         dzp(1)=0.d0
4870         if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.1)then
4871         zpmax(1)=-log(xpmin)
4872         zpmin(1)=-log(xpmax)
4873         if(abs(xpmin-xmin).lt.eps)then
4874           zpmin(1)=-log(min(xpmax,1.d0-xpmin-xpmin))
4875           zpmax(1)=-log(max(xpmin,1.d0-xpmax-xpmax))
4876         endif
4877         dzp(1)=(zpmax(1)-zpmin(1))/2.d0
4878         endif
4879
4880 c        write(*,*)'bornes1=',exp(-zpmax(1)),exp(-zpmin(1))
4881 c     &,exp(-zmmax(1)),exp(-zmmin(1))
4882
4883         Gp1=0.d0
4884         do np1=1,nmax1
4885         do jp1=1,jmax1
4886           zp(1)=zpmin(1)+dzp(1)*(1.d0+dble(2.*np1-3.)*dble(x1(jp1)))
4887           Gm1=0.d0
4888           if(dzm(1).eq.0.d0)then
4889             nmax1=1
4890             jmax1=1
4891           endif
4892           do nm1=1,nmax1
4893             do jm1=1,jmax1
4894               if(dzm(1).ne.0.d0)then
4895               zm(1)=zmmin(1)+dzm(1)*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1)))
4896               else
4897               zm(1)=zp(1)
4898               endif
4899
4900               if(mtmp.eq.1)then
4901               sxp=xp
4902               sxm=xm
4903               do i=1,mtmp
4904                 sxp=sxp-exp(-zp(i))
4905                 sxm=sxm-exp(-zm(i))
4906               enddo
4907               pGtab=1.d0
4908               k=0
4909               do l=1,mtmp
4910                 k=k+1
4911                 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4912                   Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4913                   pGtab=pGtab*Gst
4914                   if(Gst.eq.0.d0)
4915      &                write(*,*)'j1=',k,exp(-zp(k)),exp(-zm(k))
4916      &     ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1
4917                 else
4918                   pGtab=0.d0
4919                   write(*,*)'error1 ?',dzp(k),dzm(k)
4920                 endif
4921               enddo
4922               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
4923                 if(dzm(1).ne.0.d0)then
4924                   Gm1=Gm1+pGtab*
4925      &dble(a1(jm1))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
4926      &*exp(-zm(1))
4927 c     &dble(a1(jm1))*Phiexact(0.,0.,1.,sxp,sxm,sy,b)
4928 c     &*exp(-zm(1))
4929                   else
4930                   Gp1=Gp1+pGtab*
4931      &dble(a1(jp1))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
4932      &*exp(-zp(1))
4933 c     &dble(a1(jp1))*Phiexact(0.,0.,1.,sxp,sxm,sy,b)
4934 c     &*exp(-zp(1))
4935                 endif
4936 c          write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm
4937               endif
4938               endif
4939
4940         dzp(2)=0.d0
4941         if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.2)then
4942         zpmin(2)=-log(min(min(xpmax,1.d0-exp(-zp(1))),
4943      &                     1.d0-xpmin-exp(-zp(1))))
4944         zpmax(2)=-log(max(xpmin,1.d0-xpmax-exp(-zp(1))))
4945       if(abs(xpmax+xpmax+xpmax-3.d0*dble(1./delx)).lt.eps)then
4946           zpmin(2)=-log(xpmax)
4947           zpmax(2)=-log(xpmin)
4948         endif
4949           dzp(2)=(zpmax(2)-zpmin(2))/2.d0
4950         endif
4951
4952         dzm(2)=0.d0
4953         if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.2)then
4954             zmmin(2)=-log(min(min(xmmax,1.d0-exp(-zm(1))),
4955      &                     1.d0-xmmin-exp(-zm(1))))
4956             zmmax(2)=-log(max(xmmin,1.d0-xmmax-exp(-zm(1))))
4957       if(abs(xmmax+xmmax+xmmax-3.d0*dble(1./delx)).lt.eps)then
4958             zmmin(2)=-log(xmmax)
4959             zmmax(2)=-log(xmmin)
4960           endif
4961           dzm(2)=(zmmax(2)-zmmin(2))/2.d0
4962         endif
4963 c        write(*,*)'bornes2=',exp(-zpmax(2)),exp(-zpmin(2))
4964 c     &,exp(-zmmax(2)),exp(-zmmin(2)),xpmax(2),1.d0-exp(-zp(1))
4965 c     &,1.d0-xpmin(3)-exp(-zp(1)),xpmin(2),1.d0-xpmax(3)-exp(-zp(1))
4966
4967         Gp2=0.d0
4968         do np2=1,nmax2
4969         do jp2=1,jmax2
4970           zp(2)=zpmin(2)+dzp(2)*(1.d0+dble(2.*np2-3.)*dble(x1(jp2)))
4971           Gm2=0.d0
4972           if(dzm(2).eq.0.d0)then
4973             nmax2=1
4974             jmax2=1
4975           endif
4976           do nm2=1,nmax2
4977             do jm2=1,jmax2
4978               if(dzm(2).ne.0.d0)then
4979               zm(2)=zmmin(2)+dzm(2)*(1.d0+dble(2.*nm2-3.)*dble(x1(jm2)))
4980               else
4981               zm(2)=zp(2)
4982               endif
4983
4984               if(mtmp.eq.2)then
4985               sxp=xp
4986               sxm=xm
4987               do i=1,mtmp
4988                 sxp=sxp-exp(-zp(i))
4989                 sxm=sxm-exp(-zm(i))
4990               enddo
4991               pGtab=1.d0
4992               k=0
4993               do l=1,mtmp
4994                 k=k+1
4995                 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
4996                Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
4997                   pGtab=pGtab*Gst
4998                   if(Gst.eq.0.d0)
4999      &                write(*,*)'j2=',k,exp(-zp(k)),exp(-zm(k))
5000      &     ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2
5001                 else
5002                   pGtab=0.d0
5003                   write(*,*)'error2 ?',dzp(k),dzm(k)
5004                 endif
5005               enddo
5006               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
5007                 if(dzm(2).ne.0.d0)then
5008                   Gm2=Gm2+pGtab*
5009      &dble(a1(jm2))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
5010      &*exp(-zm(2))
5011 c     &dble(a1(jm2))*Phiexact(0.,0.,1.,sxp,sxm,sy,b,mk)
5012 c     &*exp(-zm(2))
5013                   else
5014                   Gp2=Gp2+pGtab*
5015      &dble(a1(jp2))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
5016      &*exp(-zp(2))
5017 c     &dble(a1(jp2))*Phiexact(0.,0.,1.,sxp,sxm,sy,b,mk)
5018 c     &*exp(-zp(2))
5019                 endif
5020 c          write(*,*)'m=2',mtmp,Gm2,Gp2,pGtab,sxp,sxm
5021               endif
5022               endif
5023
5024         dzp(3)=0.d0
5025         if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.3)then
5026         zpmin(3)=-log(min(xpmax,1.d0-exp(-zp(1))-exp(-zp(2))))
5027         zpmax(3)=-log(xpmin)
5028         dzp(3)=(zpmax(3)-zpmin(3))/2.d0
5029         endif
5030
5031         dzm(3)=0.d0
5032         if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.3)then
5033         zmmin(3)=-log(min(xmmax,1.d0-exp(-zm(1))-exp(-zm(2))))
5034         zmmax(3)=-log(xmmin)
5035         dzm(3)=(zmmax(3)-zmmin(3))/2.d0
5036         endif
5037
5038 c        write(*,*)'bornes3=',exp(-zpmax(3)),exp(-zpmin(3))
5039 c     &,exp(-zmmax(3)),exp(-zmmin(3))
5040
5041         Gp3=0.d0
5042         do np3=1,nmax3
5043         do jp3=1,jmax3
5044           zp(3)=zpmin(3)+dzp(3)*(1.d0+dble(2.*np3-3.)*dble(x1(jp3)))
5045           Gm3=0.d0
5046           if(dzm(3).eq.0.d0)then
5047             nmax3=1
5048             jmax3=1
5049           endif
5050           do nm3=1,nmax3
5051             do jm3=1,jmax3
5052               if(dzm(3).ne.0.d0)then
5053               zm(3)=zmmin(3)+dzm(3)*(1.d0+dble(2.*nm3-3.)*dble(x1(jm3)))
5054               else
5055               zm(3)=zp(3)
5056               endif
5057
5058               sxp=xp
5059               sxm=xm
5060               do i=1,mtmp
5061                 sxp=sxp-exp(-zp(i))
5062                 sxm=sxm-exp(-zm(i))
5063               enddo
5064               pGtab=1.d0
5065               k=0
5066               do l=1,mtmp
5067                 k=k+1
5068                 if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then
5069                Gst=omGam(exp(-zp(k)),exp(-zm(k)),b)
5070                   pGtab=pGtab*Gst
5071                   if(Gst.eq.0.d0)
5072      &                write(*,*)'j3=',k,exp(-zp(k)),exp(-zm(k))
5073      &   ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2,jp3
5074                 else
5075                   pGtab=0.d0
5076                   write(*,*)'error3 ?',k,dzp(k),dzm(k)
5077                 endif
5078               enddo
5079               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
5080                 if(dzm(3).ne.0.d0)then
5081                   Gm3=Gm3+pGtab
5082      &*dble(a1(jm3))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
5083      &*exp(-zm(3))
5084                   else
5085                   Gp3=Gp3+pGtab
5086      &*dble(a1(jp3))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)
5087      &*exp(-zp(3))
5088                 endif
5089               endif
5090             enddo
5091           enddo
5092          if(dzm(3).ne.0.d0)Gp3=Gp3+Gm3*dble(a1(jp3))*exp(-zp(3))*dzm(3)
5093          nmax3=2
5094          jmax3=7
5095         enddo
5096       enddo
5097               if(mtmp.gt.2.and.dzm(2).ne.0.d0)then
5098                 Gm2=Gm2+Gp3*dble(a1(jm2))*exp(-zm(2))*dzp(3)
5099               elseif(mtmp.gt.2)then
5100                 Gp2=Gp2+Gp3*dble(a1(jp2))*exp(-zp(2))*dzp(3)
5101               endif
5102             enddo
5103           enddo
5104          if(dzm(2).ne.0.d0)Gp2=Gp2+Gm2*dble(a1(jp2))*exp(-zp(2))*dzm(2)
5105          nmax2=2
5106          jmax2=7
5107         enddo
5108       enddo
5109               if(mtmp.gt.1.and.dzm(1).ne.0.d0)then
5110                 Gm1=Gm1+Gp2*dble(a1(jm1))*exp(-zm(1))*dzp(2)
5111               elseif(mtmp.gt.1)then
5112                 Gp1=Gp1+Gp2*dble(a1(jp1))*exp(-zp(1))*dzp(2)
5113               endif
5114             enddo
5115           enddo
5116          if(dzm(1).ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*exp(-zp(1))*dzm(1)
5117          nmax1=2
5118          jmax1=7
5119         enddo
5120       enddo
5121
5122       if(mtmp.gt.0)G0=Gp1*dzp(1)
5123       write(*,*)"int:",G0
5124
5125       GammaGauss=GammaGauss+G0
5126
5127       return
5128
5129       end
5130
5131 c-----------------------------------------------------------------------
5132       double precision function omWi(sy,b)   !---check---
5133 c-----------------------------------------------------------------------
5134 c cut enhanced diagram integrated over xp, xm, xpr,xmr
5135 c (with ideal G)
5136 c b - impact parameter between the pomeron ends;
5137 c sy- total energy
5138 c-----------------------------------------------------------------------
5139       include "epos.inc"
5140       include "epos.incpar"
5141       include "epos.incsem"
5142       include "epos.incems"
5143       common /psar7/ delx,alam3p,gam3p
5144       double precision xpmin,xmmin,zp,zm,alpp,alpm,xjacp,xjacm
5145      *,xpmax,xmmax,zpmin,zmmin,zpmax,chg
5146       double precision xp,xm,pGtab,omGam,dzp,Gp1,Gm1,xmin,eps
5147      *,sxp,sxm,PhiExpo,zmmax,dzm!,gamp,gamm,gampp,gammp
5148 c     *,PhiExact
5149       common /ar3/  x1(7),a1(7)
5150 c      double precision dgx1,dga1
5151 c      common /dga20/ dgx1(10),dga1(10)
5152
5153       omWi=0.d0
5154
5155         xmin=1.d-30
5156         eps=1.d-15
5157         chg=1.d0/dble(delx)
5158         b2=b*b
5159         gamb=gamD(1,iclpro,icltar)
5160 ctp060829        gamp=dble(gamb*b2/2.-alppar)
5161 ctp060829        gamm=dble(gamb*b2/2.-alppar)
5162 ctp060829        gampp=1.d0+2.d0*gamp
5163 ctp060829        gammp=1.d0+2.d0*gamm
5164
5165         nmax=2
5166         jmax=7
5167
5168         xpmin=xmin
5169         xmmin=xmin
5170         xpmax=1.d0
5171         xmmax=1.d0
5172         zpmin=0.d0
5173         zmmin=0.d0
5174         zpmax=0.d0
5175         zmmax=0.d0
5176         zp=0.d0
5177         zm=0.d0
5178         dzp=0.d0
5179         dzm=0.d0
5180
5181
5182
5183         do intp=1,2
5184         do intm=1,2
5185
5186           if(intp.eq.1)then
5187
5188           xpmin=xmin
5189           xpmax=chg
5190           alpp=(1.d0+2.d0*dble(gamb*b2/2.))
5191
5192           else
5193
5194           xpmin=chg
5195           xpmax=1.d0
5196           alpp=1.d0!(1.d0+2.d0*dble(gamb*b2/2.))
5197
5198           endif
5199
5200           if(intm.eq.1)then
5201
5202           xmmin=xmin
5203           xmmax=chg
5204           alpm=(1.d0+2.d0*dble(gamb*b2/2.))
5205
5206           else
5207
5208           xmmin=chg
5209           xmmax=1.d0
5210           alpm=1.d0!(1.d0+2.d0*dble(gamb*b2/2.))
5211
5212           endif
5213 c        write(*,*)'x+/-',intp,intm,xmmin,xmmax,xpmin,xpmax,alpp,alpm
5214
5215
5216         dzm=0.d0
5217         if(abs(xmmin-xmmax).ge.eps)then
5218           if(alpm.eq.0.d0)then
5219             zmmax=-log(xmmin)
5220             zmmin=-log(xmmax)
5221           else
5222             zmmin=xmmin**alpm
5223             zmmax=xmmax**alpm
5224           endif
5225           dzm=(zmmax-zmmin)/2.d0
5226         endif
5227
5228         dzp=0.d0
5229         if(abs(xpmin-xpmax).ge.eps)then
5230           if(alpp.eq.0.d0)then
5231             zpmax=-log(xpmin)
5232             zpmin=-log(xpmax)
5233           else
5234             zpmin=xpmin**alpp
5235             zpmax=xpmax**alpp
5236           endif
5237           dzp=(zpmax-zpmin)/2.d0
5238         endif
5239
5240
5241         Gp1=0.d0
5242
5243         if(abs(dzp).gt.eps.and.abs(dzm).gt.eps)then
5244 c        write(*,*)'Ca passe ...'
5245
5246         do np1=1,nmax
5247         do jp1=1,jmax
5248           zp=zpmin+dzp*(1.d0+dble(2.*np1-3.)*dble(x1(jp1)))
5249 c          zp=zpmin+dzp*(1.d0+dble(2.*np1-3.)*dgx1(jp1))
5250           if(alpp.eq.0.d0)then
5251             xp=exp(-zp)
5252             xjacp=xp
5253           else
5254             xp=zp**(1.d0/alpp)
5255             xjacp=zp**(1.d0/alpp-1.d0)/alpp
5256           endif
5257
5258           Gm1=0.d0
5259           do nm1=1,nmax
5260             do jm1=1,jmax
5261                 zm=zmmin+dzm*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1)))
5262 c                zm=zmmin+dzm*(1.d0+dble(2.*nm1-3.)*dgx1(jm1))
5263                 if(alpm.eq.0.d0)then
5264                   xm=exp(-zm)
5265                   xjacm=xm
5266                 else
5267                   xm=zm**(1.d0/alpm)
5268                   xjacm=zm**(1.d0/alpm-1.d0)/alpm
5269                 endif
5270
5271               sxp=1.d0-xp
5272               sxm=1.d0-xm
5273               pGtab=1.d0
5274               if(dzp.ge.0.d0.and.dzm.ge.0.d0)then
5275                 pGtab=omGam(xp,xm,b)
5276                 if(pGtab.eq.0.d0)
5277      &          write(*,*)'j1=',xp,xm,xmmin,xmmax,dzp,dzm,jp1
5278               else
5279                 pGtab=0.d0
5280                 write(*,*)'error ?',dzp,dzm
5281               endif
5282               if(sxp.gt.0.d0.and.sxm.gt.0.d0)then
5283                 if(dzm.ne.0.d0)then
5284                   Gm1=Gm1+pGtab*
5285      &dble(a1(jm1))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)*xjacm
5286 c     &dga1(jm1)*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)*xjacm
5287 c     &dble(a1(jm1))*Phiexact(0.,0.,1.,sxp,sxm,sy,b)*xjacm
5288                   else
5289                   Gp1=Gp1+pGtab*
5290      &dble(a1(jp1))*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)*xjacp
5291 c     &dga1(jp1)*Phiexpo(0.,0.,1.,sxp,sxm,sy,b)*xjacp
5292 c     &dble(a1(jp1))*Phiexact(0.,0.,1.,sxp,sxm,sy,b)*xjacp
5293                 endif
5294 c          write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm
5295               endif
5296
5297             enddo
5298           enddo
5299           if(dzm.ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*dzm*xjacp
5300 c          if(dzm.ne.0.d0)Gp1=Gp1+Gm1*dga1(jp1)*dzm*xjacp
5301         enddo
5302         enddo
5303         endif
5304
5305         omWi=omWi+Gp1*dzp
5306
5307         enddo
5308         enddo
5309
5310
5311       return
5312       end
5313
5314 c-----------------------------------------------------------------------
5315       double precision function Womint(sy,bh)   !---check---
5316 c-----------------------------------------------------------------------
5317 c - chi~(xp,xm)/2. for group of cut enhanced diagram giving
5318 c the same final state integrated over xpr and xmr (with ideal G)
5319 c bh - impact parameter between the pomeron ends;
5320 c xh - fraction of the energy squared s for the pomeron;
5321 c yp - rapidity for the pomeron;
5322 c-----------------------------------------------------------------------
5323       include 'epos.inc'
5324       double precision omWi
5325
5326       Womint=omWi(sy,bh)
5327
5328       return
5329       end
5330
5331
5332
5333 c-----------------------------------------------------------------------
5334       double precision function WomGamint(bh)   !---check---
5335 c-----------------------------------------------------------------------
5336 c - chi~(xp,xm)/2. for group of integrated cut enhanced diagram giving
5337 c the same final for proposal.
5338 c bh - impact parameter between the pomeron ends;
5339 c xh - fraction of the energy squared s for the pomeron;
5340 c yp - rapidity for the pomeron;
5341 c-----------------------------------------------------------------------
5342       include 'epos.inc'
5343       double precision omGamint
5344
5345       WomGamint=omGamint(bh)
5346
5347       return
5348       end
5349