ReggeGribovPartonMCInterface/src/epos-xpr.f

0001 c----------------------------------------------------------------------
0002       double precision function xDfit(zz,i1,i2,s,xp,xm,b)
0003 c----------------------------------------------------------------------
0004
0005       include 'epos.inc'
0006       include 'epos.incpar'
0007       include 'epos.incsem'
0008
0009
0010
0011       xDfit=0.d0
0012       do i=max(0,i1),i2
0013         call GfunPar(zz,zz,1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
0014         if(i1.ge.0)then
0015           corp=alppar-epsp
0016           cort=alppar-epst
0017           cors=-epss
0018         else
0019           corp=alppar
0020           cort=alppar
0021           cors=0.
0022         endif
0023 c        write(ifch,*)'xdfit',i,zz,b,s,alp,bet,betp,epsp,epst,epss
0024         xDfit=xDfit+dble(alp*xp**(bet+corp)*xm**(betp+cort)*s**cors)
0025       enddo
0026       return
0027       end
0028
0029
0030 c----------------------------------------------------------------------
0031       subroutine xFitD1
0032 c----------------------------------------------------------------------
0033
0034       include 'epos.inc'
0035       include 'epos.incsem'
0036       include 'epos.incpar'
0037       double precision x,y,Dsoftshval,om51p,xminr,tmp,xtmp,xDfit
0038       character chenergy*12
0039
0040       nptg=50                   !number of point for the graphs
0041       iii=nint(xpar1)
0042       biniDf=xpar2                   !value of biniDf (impact parameter)
0043       y=dble(xpar3)                    !value of y (rapidity)
0044       xtmp=xmaxDf
0045       xmaxDf=dexp(-2.d0*y)
0046       zz=xpar7
0047
0048       chenergy='E=          '
0049       if (engy.ge.10000.) then
0050         write(chenergy(4:8),'(I5)')int(engy)
0051         ke=10
0052       elseif (engy.ge.1000.) then
0053         write(chenergy(4:7),'(I4)')int(engy)
0054         ke=9
0055       elseif (engy.ge.100.) then
0056         write(chenergy(4:6),'(I3)')int(engy)
0057         ke=8
0058       elseif (engy.ge.10.) then
0059         write(chenergy(4:5),'(I2)')int(engy)
0060         ke=7
0061       else
0062         write(chenergy(4:4),'(I1)')int(engy)
0063         ke=6
0064       endif
0065       chenergy(ke:ke+2)='GeV'
0066
0067       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
0068
0069       if(iii/10.eq.1)then !...................................................
0070
0071       write(ifhi,'(a)')'!----------------------------------------------'
0072       write(ifhi,'(a)')'!     D exact all      (blue)      '
0073       write(ifhi,'(a)')'!----------------------------------------------'
0074
0075       write(ifhi,'(a)')'openhisto name DExact-'//chenergy(4:ke-2)
0076       write(ifhi,'(a)')       'htyp lbu'
0077       write(ifhi,'(a)')       'xmod log ymod log'
0078       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0079       write(ifhi,'(a)')       'yrange .01 auto'
0080       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0081       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0082       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact" '
0083       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0084         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0085         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0086       write(ifhi,'(a)')       'array 2'
0087
0088       do i=0,nptg
0089         x=xminr
0090         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0091         tmp=Dsoftshval(real(x)*smaxDf,x,y,biniDf,0)
0092         write(ifhi,*) x,tmp
0093       enddo
0094
0095       write(ifhi,'(a)')    '  endarray'
0096       write(ifhi,'(a)')    'closehisto plot 0-'
0097
0098
0099       write(ifhi,'(a)')'!----------------------------------------------'
0100       write(ifhi,'(a)')'!     D exact soft      (red dot)      '
0101       write(ifhi,'(a)')'!----------------------------------------------'
0102
0103       write(ifhi,'(a)')'openhisto name DExactSoft-'//chenergy(4:ke-2)
0104       write(ifhi,'(a)')       'htyp lra'
0105       write(ifhi,'(a)')       'xmod log ymod log'
0106       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0107       write(ifhi,'(a)')       'yrange .01 auto'
0108       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0109       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0110       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0111       if (xpar8.eq.1.) then
0112         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0113         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0114       endif
0115       write(ifhi,'(a)')       'array 2'
0116
0117       do i=0,nptg
0118         x=xminr
0119         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0120         write(ifhi,*) x,2.d0*om51p(real(x)*smaxDf,x,y,biniDf,0)
0121      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
0122       enddo
0123
0124       write(ifhi,'(a)')    '  endarray'
0125       write(ifhi,'(a)')    'closehisto plot 0-'
0126
0127
0128       write(ifhi,'(a)')'!---------------------------------------------'
0129       write(ifhi,'(a)')'!     D exact sea-sea      (yellow-dot)     '
0130       write(ifhi,'(a)')'!---------------------------------------------'
0131
0132       write(ifhi,'(a)')'openhisto name DExactSemi-'//chenergy(4:ke-2)
0133       write(ifhi,'(a)')       'htyp lya'
0134       write(ifhi,'(a)')       'xmod log ymod log'
0135       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0136       write(ifhi,'(a)')       'yrange .01 auto'
0137       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0138       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0139       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0140       if (xpar8.eq.1.) then
0141         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0142         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0143       endif
0144       write(ifhi,'(a)')       'array 2'
0145
0146       do i=0,nptg
0147         x=xminr
0148         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0149         tmp=2.d0*om51p(real(x)*smaxDf,x,y,biniDf,1)
0150      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
0151         write(ifhi,*) x,tmp
0152       enddo
0153
0154       write(ifhi,'(a)')    '  endarray'
0155       write(ifhi,'(a)')    'closehisto plot 0-'
0156
0157       write(ifhi,'(a)')'!---------------------------------------------'
0158       write(ifhi,'(a)')'!     D exact semi      (blue dot)          '
0159       write(ifhi,'(a)')'!---------------------------------------------'
0160
0161       write(ifhi,'(a)')'openhisto name DExactVal-'//chenergy(4:ke-2)
0162       write(ifhi,'(a)')       'htyp lba'
0163       write(ifhi,'(a)')       'xmod log ymod log'
0164       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0165       write(ifhi,'(a)')       'yrange .01 auto'
0166       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0167       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0168       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0169       if (xpar8.eq.1.) then
0170         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0171         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0172       endif
0173       write(ifhi,'(a)')       'array 2'
0174
0175       do i=0,nptg
0176         x=xminr
0177         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0178         tmp=2.d0*(om51p(real(x)*smaxDf,x,y,biniDf,2)
0179      &       +om51p(real(x)*smaxDf,x,y,biniDf,1)
0180      &       +om51p(real(x)*smaxDf,x,y,biniDf,3)
0181      &       +om51p(real(x)*smaxDf,x,y,biniDf,4))
0182      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
0183         write(ifhi,*) x,tmp
0184       enddo
0185
0186       write(ifhi,'(a)')    '  endarray'
0187
0188       if(iii.eq.11)then
0189         write(ifhi,'(a)')    'closehisto plot 0-'
0190       else
0191         write(ifhi,'(a)')    'closehisto plot 0'
0192       endif
0193
0194       endif !................................................................
0195       if(mod(iii,10).eq.1)then !.............................................
0196
0197       write(ifhi,'(a)')'!---------------------------------------------'
0198       write(ifhi,'(a)')'!     D exact all      (blue)              '
0199       write(ifhi,'(a)')'!---------------------------------------------'
0200
0201       write(ifhi,'(a)')'openhisto name DExact-'//chenergy(4:ke-2)
0202       write(ifhi,'(a)')       'htyp lbu'
0203       write(ifhi,'(a)')       'xmod log ymod log'
0204       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0205       write(ifhi,'(a)')       'yrange .01 auto'
0206       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0207       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0208       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact+fit" '
0209       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0210       if (xpar8.eq.1.) then
0211         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0212         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0213       endif
0214       write(ifhi,'(a)')       'array 2'
0215
0216       do i=0,nptg
0217         x=xminr
0218         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0219         tmp=Dsoftshval(real(x)*smaxDf,x,y,biniDf,0)
0220         write(ifhi,*) x,tmp
0221       enddo
0222
0223       write(ifhi,'(a)')    '  endarray'
0224       write(ifhi,'(a)')    'closehisto plot 0-'
0225
0226       write(ifhi,'(a)')'!---------------------------------------------'
0227       write(ifhi,'(a)')'!     D exact all      (green)              '
0228       write(ifhi,'(a)')'!---------------------------------------------'
0229
0230       write(ifhi,'(a)')'openhisto name DExact-f-'//chenergy(4:ke-2)
0231       write(ifhi,'(a)')       'htyp lga'
0232       write(ifhi,'(a)')       'xmod log ymod log'
0233       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0234       write(ifhi,'(a)')       'yrange .01 auto'
0235       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0236       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0237       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact+fit" '
0238       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0239       if (xpar8.eq.1.) then
0240         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0241         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0242       endif
0243       write(ifhi,'(a)')       'array 2'
0244
0245       do i=0,nptg
0246         x=xminr
0247         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0248         tmp=Dsoftshval(real(x)*smaxDf,x,y,biniDf,-1)
0249         write(ifhi,*) x,tmp
0250       enddo
0251
0252       write(ifhi,'(a)')    '  endarray'
0253       write(ifhi,'(a)')    'closehisto plot 0-'
0254
0255       write(ifhi,'(a)')'!---------------------------------------------'
0256       write(ifhi,'(a)')'!     fit soft      (red dot)             '
0257       write(ifhi,'(a)')'!---------------------------------------------'
0258
0259       write(ifhi,'(a)')       'openhisto'
0260       write(ifhi,'(a)')       'htyp lro'
0261       write(ifhi,'(a)')       'xmod log ymod log'
0262       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0263       write(ifhi,'(a)')       'yrange .01 auto'
0264       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0265       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0266       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0267       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0268       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0269       write(ifhi,'(a)')       'array 2'
0270
0271
0272       do i=0,nptg
0273         x=xminr
0274         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0275         xp=sqrt(real(x))*exp(real(y))
0276         xm=sqrt(real(x))*exp(-real(y))
0277         tmp=xDfit(zz,0,0,smaxDf,xp,xm,biniDf)
0278         write(ifhi,*) x,tmp
0279       enddo
0280       write(ifhi,'(a)')    '  endarray'
0281       write(ifhi,'(a)')    'closehisto plot 0-'
0282
0283       write(ifhi,'(a)')'!---------------------------------------------'
0284       write(ifhi,'(a)')'!     fit semi      (blue dot)              '
0285       write(ifhi,'(a)')'!---------------------------------------------'
0286
0287       write(ifhi,'(a)')       'openhisto'
0288       write(ifhi,'(a)')       'htyp lbo'
0289       write(ifhi,'(a)')       'xmod log ymod log'
0290       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0291       write(ifhi,'(a)')       'yrange .01 auto'
0292       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0293       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,,s,b)" '
0294       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0295       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0296       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0297       write(ifhi,'(a)')       'array 2'
0298
0299
0300       do i=0,nptg
0301         x=xminr
0302         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0303         xp=sqrt(real(x))*exp(real(y))
0304         xm=sqrt(real(x))*exp(-real(y))
0305         tmp=xDfit(zz,1,1,smaxDf,xp,xm,biniDf)
0306         write(ifhi,*) x,tmp
0307       enddo
0308       write(ifhi,'(a)')    '  endarray'
0309       write(ifhi,'(a)')    'closehisto plot 0-'
0310
0311       write(ifhi,'(a)')'!---------------------------------------------'
0312       write(ifhi,'(a)')'!     fit all      (red)      '
0313       write(ifhi,'(a)')'!---------------------------------------------'
0314
0315       write(ifhi,'(a)')       'openhisto'
0316       write(ifhi,'(a)')       'htyp lru'
0317       write(ifhi,'(a)')       'xmod log ymod log'
0318       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0319       write(ifhi,'(a)')       'yrange .01 auto'
0320       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0321       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,,s,b)" '
0322       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
0323       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
0324       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
0325       write(ifhi,'(a)')       'array 2'
0326
0327       do i=0,nptg
0328         x=xminr
0329         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0330         xp=sqrt(real(x))*exp(real(y))
0331         xm=sqrt(real(x))*exp(-real(y))
0332         tmp=xDfit(zz,0,1,smaxDf,xp,xm,biniDf)
0333         write(ifhi,*) x,tmp
0334       enddo
0335       write(ifhi,'(a)')    '  endarray'
0336       write(ifhi,'(a)')    'closehisto plot 0'
0337
0338       endif !...............................................................
0339
0340       xmaxDf=xtmp
0341
0342       end
0343
0344 c----------------------------------------------------------------------
0345       subroutine xFitD2
0346 c----------------------------------------------------------------------
0347
0348       include 'epos.inc'
0349       include 'epos.incsem'
0350       include 'epos.incpar'
0351
0352       double precision x,om51p,xDfit,z(0:200),xminr,y,xtmp,om,om5,om51
0353 c     & ,omYuncut,omNpuncut
0354       character chenergy*12,chf*3,texte*15,textb*17,texty*17
0355
0356       nptg=30                  !number of point for the graphs
0357       biniDf=xpar2                 !value of biniDf (impact parameter)
0358       y=dble(xpar3)                 !value of y (rapidity)
0359       jj1=nint(xpar4)
0360       jj2=nint(xpar5)
0361       if(jj1.ne.1.and.jj1.ne.2)jj1=3
0362       if(jj2.ne.1.and.jj2.ne.2)jj2=3
0363       zz=xpar7
0364       xtmp=xmaxDf
0365       xmaxDf=dexp(-2.d0*y)
0366
0367       chenergy='E=          '
0368       if (engy.ge.10000.) then
0369         write(chenergy(4:8),'(I5)')int(engy)
0370         ke=10
0371       elseif (engy.ge.1000.) then
0372         write(chenergy(4:7),'(I4)')int(engy)
0373         ke=9
0374       elseif (engy.ge.100.) then
0375         write(chenergy(4:6),'(I3)')int(engy)
0376         ke=8
0377       elseif (engy.ge.10.) then
0378         write(chenergy(4:5),'(I2)')int(engy)
0379         ke=7
0380       else
0381         write(chenergy(4:4),'(I1)')int(engy)
0382         ke=6
0383       endif
0384       chenergy(ke:ke+2)='GeV'
0385
0386       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
0387
0388          do jj=jj1,jj2
0389
0390       if(jj.eq.1)chf='  D'
0391       if(jj.eq.2)chf='  G'
0392       if(jj.eq.3)chf='FFG'
0393       texte='text 0.05 0.9 "'
0394       textb='text 0.05 0.8 "b='
0395       texty='text 0.05 0.7 "y='
0396       if(jj.eq.2)texty='text 0.15 0.7 "y='
0397       if(jj.eq.3)texte='text 0.05 0.3 "'
0398       if(jj.eq.3)textb='text 0.05 0.2 "b='
0399       if(jj.eq.3)texty='text 0.05 0.1 "y='
0400
0401       write(ifhi,'(a)')'!---------------------------------------------'
0402       write(ifhi,'(a)')'!     '//chf//' exact all      (green)         '
0403       write(ifhi,'(a)')'!---------------------------------------------'
0404
0405       write(ifhi,'(a)')'openhisto name '//chf//'ExaI-'//chenergy(4:ke-2)
0406       write(ifhi,'(a)')       'htyp lga'
0407       write(ifhi,'(a)')       'xmod log ymod log'
0408       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0409       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0410       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
0411       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact+fit" '
0412       write(ifhi,'(3a)')        texte,chenergy,'"'
0413       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
0414       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
0415       write(ifhi,'(a)')       'array 2'
0416       do i=0,nptg
0417         x=xminr
0418         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0419         xp=sqrt(real(x))*exp(real(y))
0420         xm=sqrt(real(x))*exp(-real(y))
0421         om=0
0422         do j=0,4
0423          om=om+om51p(engy**2*real(x),x,y,biniDf,j)
0424         enddo
0425         om=2.d0*om
0426         if(jj.eq.1)om=om/(x**dble(-alppar))
0427         if(jj.eq.3)om=om
0428      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
0429         write(ifhi,*) x,om
0430       enddo
0431       write(ifhi,'(a)')    '  endarray'
0432       write(ifhi,'(a)')    'closehisto plot 0-'
0433
0434       write(ifhi,'(a)')'!---------------------------------------------'
0435       write(ifhi,'(a)')'!     '//chf//' exact all +diff (blue)        '
0436       write(ifhi,'(a)')'!---------------------------------------------'
0437
0438       write(ifhi,'(a)')'openhisto name '//chf//'ExaD-'//chenergy(4:ke-2)
0439       write(ifhi,'(a)')       'htyp lbu'
0440       write(ifhi,'(a)')       'xmod log ymod log'
0441       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0442       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0443       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
0444       write(ifhi,'(3a)')        texte,chenergy,'"'
0445       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
0446       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
0447       write(ifhi,'(a)')       'array 2'
0448       do i=0,nptg
0449         x=xminr
0450         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0451         xp=sqrt(real(x))*exp(real(y))
0452         xm=sqrt(real(x))*exp(-real(y))
0453         om=0
0454         do j=0,4
0455          om=om+om51p(engy**2*real(x),x,y,biniDf,j)
0456         enddo
0457         om5=om51(x,y,biniDf,5,5)
0458         om=2.d0*(om+om5)
0459         if(jj.eq.1)om=om/(x**dble(-alppar))
0460         if(jj.eq.3)om=om
0461      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
0462         write(ifhi,*) x,om
0463       enddo
0464       write(ifhi,'(a)')    '  endarray'
0465       write(ifhi,'(a)')    'closehisto plot 0-'
0466
0467       write(ifhi,'(a)')'!---------------------------------------------'
0468       write(ifhi,'(a)')'!     '//chf//' param all      (red)          '
0469       write(ifhi,'(a)')'!---------------------------------------------'
0470
0471       write(ifhi,'(a)')'openhisto name '//chf//'Par-'//chenergy(4:ke-2)
0472       write(ifhi,'(a)')       'htyp lru'
0473         write(ifhi,'(a)')       'xmod log ymod log'
0474       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0475       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0476       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
0477       write(ifhi,'(3a)')        texte,chenergy,'"'
0478       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
0479       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
0480       write(ifhi,'(a)')       'array 2'
0481       imax=idxD1
0482       if(iomega.eq.2)imax=1
0483       do i=0,nptg
0484         x=xminr
0485         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0486         xp=sqrt(real(x))*exp(real(y))
0487         xm=sqrt(real(x))*exp(-real(y))
0488           z(i)=xDfit(zz,0,imax,engy**2,xp,xm,biniDf)
0489          if(jj.ge.2)z(i)=z(i)*(x**dble(-alppar))
0490          if(jj.eq.3)z(i)=z(i)
0491      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
0492          write(ifhi,*) x,z(i)
0493        enddo
0494       write(ifhi,'(a)')    '  endarray'
0495       write(ifhi,'(a)')    'closehisto plot 0-'
0496
0497       write(ifhi,'(a)')'!---------------------------------------------'
0498       write(ifhi,'(a)')'!     '//chf//' param all      (yellow)       '
0499       write(ifhi,'(a)')'!---------------------------------------------'
0500
0501       write(ifhi,'(a)')'openhisto name '//chf//'Scr-'//chenergy(4:ke-2)
0502       write(ifhi,'(a)')       'htyp lyi'
0503         write(ifhi,'(a)')       'xmod log ymod log'
0504       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
0505       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
0506       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
0507       write(ifhi,'(3a)')        texte,chenergy,'"'
0508       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
0509       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
0510       write(ifhi,'(a)')       'array 2'
0511       imax=idxD1
0512       if(iomega.eq.2)imax=1
0513       do i=0,nptg
0514         x=xminr
0515         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
0516         xp=sqrt(real(x))*exp(real(y))
0517         xm=sqrt(real(x))*exp(-real(y))
0518           z(i)=xDfit(zz,-1,imax,engy**2,xp,xm,biniDf)
0519          if(jj.ge.2)z(i)=z(i)*(x**dble(-alppar))
0520          if(jj.eq.3)z(i)=z(i)
0521      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
0522          write(ifhi,*) x,z(i)
0523        enddo
0524       write(ifhi,'(a)')    '  endarray'
0525       write(ifhi,'(a)')    'closehisto plot 0'
0526
0527          enddo
0528
0529       xmaxDf=xtmp
0530
0531       end
0532
0533 c----------------------------------------------------------------------
0534       subroutine xbExaD
0535 c----------------------------------------------------------------------
0536
0537       include 'epos.inc'
0538       include 'epos.incsem'
0539       include 'epos.incpar'
0540
0541       double precision x,y,Dsoftshval,om51p,z,xDfit!,omNpuncut
0542
0543
0544       nptg=50     !number of point for the graphs
0545       bmax=xpar2
0546       bmax=max(0.1,bmax)
0547              !value max of b (impact parameter)
0548       y=dble(xpar3)                    !value of y (rapidity)
0549       x=dble(xpar4)
0550       zz=xpar7
0551       if(x.eq.0.d0)x=1.d0
0552
0553       if (engy.ge.10.) then
0554        if (engy.ge.100.) then
0555         if (engy.ge.1000.) then
0556          if (engy.ge.10000.) then
0557       write(ifhi,'(a,I5)')  'openhisto name DExactb-',int(engy)
0558          else
0559       write(ifhi,'(a,I4)')  'openhisto name DExactb-',int(engy)
0560         endif
0561         else
0562       write(ifhi,'(a,I3)')  'openhisto name DExactb-',int(engy)
0563        endif
0564        else
0565       write(ifhi,'(a,I2)')  'openhisto name DExactb-',int(engy)
0566       endif
0567       else
0568       write(ifhi,'(a,I1)')  'openhisto name DExactb-',int(engy)
0569       endif
0570       write(ifhi,'(a)')       'htyp lbu'
0571       write(ifhi,'(a)')       'xmod lin ymod lin'
0572       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0573       write(ifhi,'(a)')       'yrange .00001 auto'
0574 c      write(ifhi,'(a)')       'yrange auto auto'
0575       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0576       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0577       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0578       if (xpar8.eq.1.) then
0579       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0580       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0581       endif
0582       write(ifhi,'(a)')       'array 2'
0583
0584       do i=0,nptg
0585         b=-bmax+2.*real(i)/real(nptg)*bmax
0586         z=Dsoftshval(real(x)*smaxDf,x,y,b,0)
0587         write(ifhi,*) b,z
0588       enddo
0589
0590       write(ifhi,'(a)')    '  endarray'
0591       write(ifhi,'(a)')    'closehisto plot 0-'
0592
0593 c**********************************************************************
0594
0595       if (engy.ge.10.) then
0596        if (engy.ge.100.) then
0597         if (engy.ge.1000.) then
0598          if (engy.ge.10000.) then
0599       write(ifhi,'(a,I5)')       'openhisto name DParamb-',int(engy)
0600          else
0601       write(ifhi,'(a,I4)')       'openhisto name DParamb-',int(engy)
0602         endif
0603         else
0604       write(ifhi,'(a,I3)')       'openhisto name DParamb-',int(engy)
0605        endif
0606        else
0607       write(ifhi,'(a,I2)')       'openhisto name DParamb-',int(engy)
0608       endif
0609       else
0610       write(ifhi,'(a,I1)')       'openhisto name DParamb-',int(engy)
0611       endif
0612       write(ifhi,'(a)')       'htyp lrd'
0613       write(ifhi,'(a)')       'xmod lin ymod lin'
0614       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0615 c      write(ifhi,'(a)')       'yrange -.01 auto'
0616       write(ifhi,'(a)')       'yrange auto auto'
0617       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0618       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0619       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0620       if (xpar8.eq.1.) then
0621       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0622       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0623       endif
0624       write(ifhi,'(a)')       'array 2'
0625
0626       imax=idxD1
0627       if(iomega.eq.2)imax=1
0628       do i=0,nptg
0629         b=-bmax+2.*real(i)/real(nptg)*bmax
0630           z=xDfit(zz,0,imax,smaxDf,
0631      &       real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),b)
0632         write(ifhi,*) b,z
0633       enddo
0634
0635       write(ifhi,'(a)')    '  endarray'
0636       write(ifhi,'(a)')    'closehisto plot 0-'
0637
0638
0639 c**********************************************************************
0640
0641       if (engy.ge.10.) then
0642        if (engy.ge.100.) then
0643         if (engy.ge.1000.) then
0644          if (engy.ge.10000.) then
0645       write(ifhi,'(a,I5)')  'openhisto name DExactSoftb-',int(engy)
0646          else
0647       write(ifhi,'(a,I4)')  'openhisto name DExactSoftb-',int(engy)
0648         endif
0649         else
0650       write(ifhi,'(a,I3)')  'openhisto name DExactSoftb-',int(engy)
0651        endif
0652        else
0653       write(ifhi,'(a,I2)')  'openhisto name DExactSoftb-',int(engy)
0654       endif
0655       else
0656       write(ifhi,'(a,I1)')  'openhisto name DExactSoftb-',int(engy)
0657       endif
0658       write(ifhi,'(a)')       'htyp pfc'
0659       write(ifhi,'(a)')       'xmod lin ymod lin'
0660       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0661 c      write(ifhi,'(a)')       'yrange -.01 auto'
0662       write(ifhi,'(a)')       'yrange auto auto'
0663       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0664       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0665       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0666       if (xpar8.eq.1.) then
0667       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0668       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0669       endif
0670       write(ifhi,'(a)')       'array 2'
0671
0672       do i=0,nptg
0673         b=-bmax+2.*real(i)/real(nptg)*bmax
0674         z=2.d0*om51p(real(x)*smaxDf,x,y,b,0)
0675      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
0676         write(ifhi,*) b,z
0677       enddo
0678
0679       write(ifhi,'(a)')    '  endarray'
0680       write(ifhi,'(a)')    'closehisto plot 0-'
0681
0682
0683 c**********************************************************************
0684
0685       if (engy.ge.10.) then
0686        if (engy.ge.100.) then
0687         if (engy.ge.1000.) then
0688          if (engy.ge.10000.) then
0689       write(ifhi,'(a,I5)')  'openhisto name DExactSemib-',int(engy)
0690          else
0691       write(ifhi,'(a,I4)')  'openhisto name DExactSemib-',int(engy)
0692         endif
0693         else
0694       write(ifhi,'(a,I3)')  'openhisto name DExactSemib-',int(engy)
0695        endif
0696        else
0697       write(ifhi,'(a,I2)')  'openhisto name DExactSemib-',int(engy)
0698       endif
0699       else
0700       write(ifhi,'(a,I1)')  'openhisto name DExactSemib-',int(engy)
0701       endif
0702       write(ifhi,'(a)')       'htyp pft'
0703       write(ifhi,'(a)')       'xmod lin ymod lin'
0704       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0705 c      write(ifhi,'(a)')       'yrange -.01 auto'
0706       write(ifhi,'(a)')       'yrange auto auto'
0707       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0708       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0709       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0710       if (xpar8.eq.1.) then
0711       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0712       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0713       endif
0714       write(ifhi,'(a)')       'array 2'
0715
0716       do i=0,nptg
0717         b=-bmax+2.*real(i)/real(nptg)*bmax
0718         z=2.d0*om51p(real(x)*smaxDf,x,y,b,1)
0719      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
0720         write(ifhi,*) b,z
0721       enddo
0722
0723       write(ifhi,'(a)')    '  endarray'
0724       write(ifhi,'(a)')    'closehisto plot 0-'
0725
0726
0727 c**********************************************************************
0728
0729       if (engy.ge.10.) then
0730        if (engy.ge.100.) then
0731         if (engy.ge.1000.) then
0732          if (engy.ge.10000.) then
0733       write(ifhi,'(a,I5)')  'openhisto name DExactValb-',int(engy)
0734          else
0735       write(ifhi,'(a,I4)')  'openhisto name DExactValb-',int(engy)
0736         endif
0737         else
0738       write(ifhi,'(a,I3)')  'openhisto name DExactValb-',int(engy)
0739        endif
0740        else
0741       write(ifhi,'(a,I2)')  'openhisto name DExactValb-',int(engy)
0742       endif
0743       else
0744       write(ifhi,'(a,I1)')  'openhisto name DExactValb-',int(engy)
0745       endif
0746       write(ifhi,'(a)')       'htyp pfs'
0747       write(ifhi,'(a)')       'xmod lin ymod lin'
0748       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0749 c      write(ifhi,'(a)')       'yrange -.01 auto'
0750       write(ifhi,'(a)')       'yrange auto auto'
0751       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0752       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0753       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0754       if (xpar8.eq.1.) then
0755       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0756       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0757       endif
0758       write(ifhi,'(a)')       'array 2'
0759
0760       do i=0,nptg
0761         b=-bmax+2.*real(i)/real(nptg)*bmax
0762         z=2.d0*(om51p(real(x)*smaxDf,x,y,b,2)+
0763      &  om51p(real(x)*smaxDf,x,y,b,3)+om51p(real(x)*smaxDf,x,y,b,4))
0764      &           /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
0765         write(ifhi,*) b,z
0766       enddo
0767
0768       write(ifhi,'(a)')    '  endarray'
0769
0770       write(ifhi,'(a)')    'closehisto plot 0'
0771
0772
0773       end
0774
0775
0776 c----------------------------------------------------------------------
0777       subroutine xbnExaD
0778 c----------------------------------------------------------------------
0779
0780       include 'epos.inc'
0781       include 'epos.incpar'
0782
0783       double precision x,y,Dsoftshval,z,xDfit,Dint
0784
0785
0786
0787       nptg=50     !number of point for the graphs
0788       bmax=xpar2
0789       bmax=max(0.1,bmax)
0790       y=dble(xpar2)                   !value of y (rapidity)
0791       x=dble(xpar3)
0792       zz=xpar7
0793
0794       if (engy.ge.10.) then
0795        if (engy.ge.100.) then
0796         if (engy.ge.1000.) then
0797          if (engy.ge.10000.) then
0798       write(ifhi,'(a,I5)')  'openhisto name DExactbn-',int(engy)
0799          else
0800       write(ifhi,'(a,I4)')  'openhisto name DExactbn-',int(engy)
0801         endif
0802         else
0803       write(ifhi,'(a,I3)')  'openhisto name DExactbn-',int(engy)
0804        endif
0805        else
0806       write(ifhi,'(a,I2)')  'openhisto name DExactbn-',int(engy)
0807       endif
0808       else
0809       write(ifhi,'(a,I1)')  'openhisto name DExactbn-',int(engy)
0810       endif
0811       write(ifhi,'(a)')       'htyp lbu'
0812       write(ifhi,'(a)')       'xmod lin ymod lin'
0813       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0814 c      write(ifhi,'(a)')       'yrange -.01 auto'
0815       write(ifhi,'(a)')       'yrange auto auto'
0816       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0817       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0818       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0819       if (xpar8.eq.1.) then
0820       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0821       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0822       endif
0823       write(ifhi,'(a)')       'array 2'
0824
0825       Dint=Dsoftshval(real(x)*smaxDf,x,y,0.,0)
0826       do i=0,nptg
0827         b=-bmax+2.*real(i)/real(nptg)*bmax
0828         z=Dsoftshval(real(x)*smaxDf,x,y,b,0)/Dint
0829         write(ifhi,*) b,z
0830       enddo
0831
0832       write(ifhi,'(a)')    '  endarray'
0833       write(ifhi,'(a)')    'closehisto plot 0-'
0834
0835 c**********************************************************************
0836
0837       if (engy.ge.10.) then
0838        if (engy.ge.100.) then
0839         if (engy.ge.1000.) then
0840          if (engy.ge.10000.) then
0841       write(ifhi,'(a,I5)')       'openhisto name DParambn-',int(engy)
0842          else
0843       write(ifhi,'(a,I4)')       'openhisto name DParambn-',int(engy)
0844         endif
0845         else
0846       write(ifhi,'(a,I3)')       'openhisto name DParambn-',int(engy)
0847        endif
0848        else
0849       write(ifhi,'(a,I2)')       'openhisto name DParambn-',int(engy)
0850       endif
0851       else
0852       write(ifhi,'(a,I1)')       'openhisto name DParambn-',int(engy)
0853       endif
0854       write(ifhi,'(a)')       'htyp lrd'
0855       write(ifhi,'(a)')       'xmod lin ymod lin'
0856       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0857 c      write(ifhi,'(a)')       'yrange -.01 auto'
0858       write(ifhi,'(a)')       'yrange auto auto'
0859       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0860       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0861       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0862       if (xpar8.eq.1.) then
0863       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0864       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0865       endif
0866       write(ifhi,'(a)')       'array 2'
0867
0868       imax=idxD1
0869       if(iomega.eq.2)imax=1
0870
0871
0872       Dint=xDfit(zz,0,imax,engy**2,
0873      &     real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),0.)
0874
0875       do i=0,nptg
0876         b=-bmax+2.*real(i)/real(nptg)*bmax
0877         z=xDfit(zz,0,imax,engy**2,
0878      &       real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),b)
0879
0880         write(ifhi,*) b,z/Dint
0881       enddo
0882
0883       write(ifhi,'(a)')    '  endarray'
0884       write(ifhi,'(a)')    'closehisto plot 0-'
0885
0886
0887 c**********************************************************************
0888
0889       if (engy.ge.10.) then
0890        if (engy.ge.100.) then
0891         if (engy.ge.1000.) then
0892          if (engy.ge.10000.) then
0893       write(ifhi,'(a,I5)')  'openhisto name DEfitb-',int(engy)
0894          else
0895       write(ifhi,'(a,I4)')  'openhisto name DEfitb-',int(engy)
0896         endif
0897         else
0898       write(ifhi,'(a,I3)')  'openhisto name DEfitb-',int(engy)
0899        endif
0900        else
0901       write(ifhi,'(a,I2)')  'openhisto name DEfitb-',int(engy)
0902       endif
0903       else
0904       write(ifhi,'(a,I1)')  'openhisto name DEfitb-',int(engy)
0905       endif
0906       write(ifhi,'(a)')       'htyp pfc'
0907       write(ifhi,'(a)')       'xmod lin ymod lin'
0908       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0909 c      write(ifhi,'(a)')       'yrange -.01 auto'
0910       write(ifhi,'(a)')       'yrange auto auto'
0911       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0912       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0913       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0914       if (xpar8.eq.1.) then
0915       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0916       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0917       endif
0918       write(ifhi,'(a)')       'array 2'
0919
0920       sig2=sigma2(x,2)
0921       if(sig2.le.0.) sig2=1.e+10
0922       do i=0,nptg
0923         b=-bmax+2.*real(i)/real(nptg)*bmax
0924         write(ifhi,*) b,exp(-b**2/sig2)
0925       enddo
0926
0927       write(ifhi,'(a)')    '  endarray'
0928       write(ifhi,'(a)')    'closehisto plot 0-'
0929
0930
0931 c**********************************************************************
0932
0933       if (engy.ge.10.) then
0934        if (engy.ge.100.) then
0935         if (engy.ge.1000.) then
0936          if (engy.ge.10000.) then
0937       write(ifhi,'(a,I5)')  'openhisto name DEfitbnSoft-',int(engy)
0938          else
0939       write(ifhi,'(a,I4)')  'openhisto name DEfitbnSoft-',int(engy)
0940         endif
0941         else
0942       write(ifhi,'(a,I3)')  'openhisto name DEfitbnSoft-',int(engy)
0943        endif
0944        else
0945       write(ifhi,'(a,I2)')  'openhisto name DEfitbnSoft-',int(engy)
0946       endif
0947       else
0948       write(ifhi,'(a,I1)')  'openhisto name DEfitbnSoft-',int(engy)
0949       endif
0950       write(ifhi,'(a)')       'htyp pft'
0951       write(ifhi,'(a)')       'xmod lin ymod lin'
0952       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0953 c      write(ifhi,'(a)')       'yrange -.01 auto'
0954       write(ifhi,'(a)')       'yrange auto auto'
0955       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
0956       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
0957       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
0958       if (xpar8.eq.1.) then
0959       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
0960       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
0961       endif
0962       write(ifhi,'(a)')       'array 2'
0963
0964       sig2=sigma2(x,0)
0965       if(sig2.le.0.) sig2=1.e+10
0966       do i=0,nptg
0967         b=-bmax+2.*real(i)/real(nptg)*bmax
0968         write(ifhi,*) b,exp(-b**2/sig2)
0969       enddo
0970
0971       write(ifhi,'(a)')    '  endarray'
0972       write(ifhi,'(a)')    'closehisto plot 0-'
0973
0974
0975 c**********************************************************************
0976
0977
0978       if (engy.ge.10.) then
0979        if (engy.ge.100.) then
0980         if (engy.ge.1000.) then
0981          if (engy.ge.10000.) then
0982       write(ifhi,'(a,I5)')  'openhisto name DEfitbnSh-',int(engy)
0983          else
0984       write(ifhi,'(a,I4)')  'openhisto name DEfitbnSh-',int(engy)
0985         endif
0986         else
0987       write(ifhi,'(a,I3)')  'openhisto name DEfitbnSh-',int(engy)
0988        endif
0989        else
0990       write(ifhi,'(a,I2)')  'openhisto name DEfitbnSh-',int(engy)
0991       endif
0992       else
0993       write(ifhi,'(a,I1)')  'openhisto name DEfitbnSh-',int(engy)
0994       endif
0995       write(ifhi,'(a)')       'htyp poc'
0996       write(ifhi,'(a)')       'xmod lin ymod lin'
0997       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
0998 c      write(ifhi,'(a)')       'yrange -.01 auto'
0999       write(ifhi,'(a)')       'yrange auto auto'
1000       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1001       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1002       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1003       if (xpar8.eq.1.) then
1004       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1005       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1006       endif
1007       write(ifhi,'(a)')       'array 2'
1008
1009       sig2=sigma2(x,1)
1010       if(sig2.le.0.) sig2=1.e+10
1011       do i=0,nptg
1012         b=-bmax+2.*real(i)/real(nptg)*bmax
1013         write(ifhi,*) b,exp(-b**2/sig2)
1014       enddo
1015
1016       write(ifhi,'(a)')    '  endarray'
1017
1018       write(ifhi,'(a)')    'closehisto plot 0'
1019
1020       end
1021
1022 c----------------------------------------------------------------------
1023       subroutine xbnParD
1024 c----------------------------------------------------------------------
1025
1026       include 'epos.inc'
1027       include 'epos.incpar'
1028
1029       double precision x,y,Dsoftshval,z,xDfit
1030
1031
1032
1033       nptg=50     !number of point for the graphs
1034       bmax=xpar2                   !value max of b (impact parameter)
1035       bmax=max(0.1,bmax)
1036       y=dble(xpar3)              !value of y (rapidity)
1037       x=dble(xpar4)
1038       zz=xpar7
1039
1040       if (engy.ge.10.) then
1041        if (engy.ge.100.) then
1042         if (engy.ge.1000.) then
1043          if (engy.ge.10000.) then
1044       write(ifhi,'(a,I5)')  'openhisto name DExactbn-',int(engy)
1045          else
1046       write(ifhi,'(a,I4)')  'openhisto name DExactbn-',int(engy)
1047         endif
1048         else
1049       write(ifhi,'(a,I3)')  'openhisto name DExactbn-',int(engy)
1050        endif
1051        else
1052       write(ifhi,'(a,I2)')  'openhisto name DExactbn-',int(engy)
1053       endif
1054       else
1055       write(ifhi,'(a,I1)')  'openhisto name DExactbn-',int(engy)
1056       endif
1057       write(ifhi,'(a)')       'htyp lbu'
1058       write(ifhi,'(a)')       'xmod lin ymod lin'
1059       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1060 c      write(ifhi,'(a)')       'yrange -.01 auto'
1061       write(ifhi,'(a)')       'yrange auto auto'
1062       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1063       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1064       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1065       if (xpar8.eq.1.) then
1066       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1067       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1068       endif
1069       write(ifhi,'(a)')       'array 2'
1070
1071       Dint=Dsoftshval(real(x)*smaxDf,x,y,0.,0)
1072       do i=0,nptg
1073         b=-bmax+2.*real(i)/real(nptg)*bmax
1074         z=Dsoftshval(real(x)*smaxDf,x,y,b,0)/Dint
1075         write(ifhi,*) b,z
1076       enddo
1077
1078       write(ifhi,'(a)')    '  endarray'
1079       write(ifhi,'(a)')    'closehisto plot 0-'
1080
1081 c**********************************************************************
1082
1083       if (engy.ge.10.) then
1084        if (engy.ge.100.) then
1085         if (engy.ge.1000.) then
1086          if (engy.ge.10000.) then
1087       write(ifhi,'(a,I5)')       'openhisto name DParambn-',int(engy)
1088          else
1089       write(ifhi,'(a,I4)')       'openhisto name DParambn-',int(engy)
1090         endif
1091         else
1092       write(ifhi,'(a,I3)')       'openhisto name DParambn-',int(engy)
1093        endif
1094        else
1095       write(ifhi,'(a,I2)')       'openhisto name DParambn-',int(engy)
1096       endif
1097       else
1098       write(ifhi,'(a,I1)')       'openhisto name DParambn-',int(engy)
1099       endif
1100       write(ifhi,'(a)')       'htyp lrd'
1101       write(ifhi,'(a)')       'xmod lin ymod lin'
1102       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1103 c      write(ifhi,'(a)')       'yrange -.01 auto'
1104       write(ifhi,'(a)')       'yrange auto auto'
1105       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1106       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1107       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1108       if (xpar8.eq.1.) then
1109       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1110       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1111       endif
1112       write(ifhi,'(a)')       'array 2'
1113
1114       Dint=xDfit(zz,0,1,smaxDf,
1115      &     real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),0.)
1116
1117       imax=idxD1
1118       if(iomega.eq.2)imax=1
1119
1120       do i=0,nptg
1121         b=-bmax+2.*real(i)/real(nptg)*bmax
1122         z=xDfit(zz,0,imax,smaxDf,
1123      &       real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),b)
1124
1125         write(ifhi,*) b,z/Dint
1126       enddo
1127
1128       write(ifhi,'(a)')    '  endarray'
1129       write(ifhi,'(a)')    'closehisto plot 0-'
1130
1131
1132 c**********************************************************************
1133
1134       if (engy.ge.10.) then
1135        if (engy.ge.100.) then
1136         if (engy.ge.1000.) then
1137          if (engy.ge.10000.) then
1138       write(ifhi,'(a,I5)')  'openhisto name DEfitb-',int(engy)
1139          else
1140       write(ifhi,'(a,I4)')  'openhisto name DEfitb-',int(engy)
1141         endif
1142         else
1143       write(ifhi,'(a,I3)')  'openhisto name DEfitb-',int(engy)
1144        endif
1145        else
1146       write(ifhi,'(a,I2)')  'openhisto name DEfitb-',int(engy)
1147       endif
1148       else
1149       write(ifhi,'(a,I1)')  'openhisto name DEfitb-',int(engy)
1150       endif
1151       write(ifhi,'(a)')       'htyp pfc'
1152       write(ifhi,'(a)')       'xmod lin ymod lin'
1153       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1154 c      write(ifhi,'(a)')       'yrange -.01 auto'
1155       write(ifhi,'(a)')       'yrange auto auto'
1156       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1157       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1158       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1159       if (xpar8.eq.1.) then
1160       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1161       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1162       endif
1163       write(ifhi,'(a)')       'array 2'
1164
1165       sig2=sigma2(x,2)
1166       if(sig2.le.0.) sig2=1.e+10
1167       do i=0,nptg
1168         b=-bmax+2.*real(i)/real(nptg)*bmax
1169         write(ifhi,*) b,exp(-b**2/sig2)
1170       enddo
1171
1172       write(ifhi,'(a)')    '  endarray'
1173       write(ifhi,'(a)')    'closehisto plot 0-'
1174
1175
1176 c**********************************************************************
1177
1178       if (engy.ge.10.) then
1179        if (engy.ge.100.) then
1180         if (engy.ge.1000.) then
1181          if (engy.ge.10000.) then
1182       write(ifhi,'(a,I5)')  'openhisto name DEintb-',int(engy)
1183          else
1184       write(ifhi,'(a,I4)')  'openhisto name DEintb-',int(engy)
1185         endif
1186         else
1187       write(ifhi,'(a,I3)')  'openhisto name DEintb-',int(engy)
1188        endif
1189        else
1190       write(ifhi,'(a,I2)')  'openhisto name DEintb-',int(engy)
1191       endif
1192       else
1193       write(ifhi,'(a,I1)')  'openhisto name DEintb-',int(engy)
1194       endif
1195       write(ifhi,'(a)')       'htyp pft'
1196       write(ifhi,'(a)')       'xmod lin ymod lin'
1197       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1198 c      write(ifhi,'(a)')       'yrange -.01 auto'
1199       write(ifhi,'(a)')       'yrange auto auto'
1200       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1201       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1202       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1203       if (xpar8.eq.1.) then
1204       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1205       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1206       endif
1207       write(ifhi,'(a)')       'array 2'
1208
1209       sig2=sigma1i(x)
1210       do i=0,nptg
1211         b=-bmax+2.*real(i)/real(nptg)*bmax
1212         write(ifhi,*) b,exp(-b**2/sig2)
1213       enddo
1214
1215       write(ifhi,'(a)')    '  endarray'
1216       write(ifhi,'(a)')    'closehisto plot 0-'
1217
1218
1219 c**********************************************************************
1220
1221
1222       if (engy.ge.10.) then
1223        if (engy.ge.100.) then
1224         if (engy.ge.1000.) then
1225          if (engy.ge.10000.) then
1226       write(ifhi,'(a,I5)')  'openhisto name DPfitb-',int(engy)
1227          else
1228       write(ifhi,'(a,I4)')  'openhisto name DPfitb-',int(engy)
1229         endif
1230         else
1231       write(ifhi,'(a,I3)')  'openhisto name DPfitb-',int(engy)
1232        endif
1233        else
1234       write(ifhi,'(a,I2)')  'openhisto name DPfitb-',int(engy)
1235       endif
1236       else
1237       write(ifhi,'(a,I1)')  'openhisto name DPfitb-',int(engy)
1238       endif
1239       write(ifhi,'(a)')       'htyp poc'
1240       write(ifhi,'(a)')       'xmod lin ymod lin'
1241       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1242 c      write(ifhi,'(a)')       'yrange -.01 auto'
1243       write(ifhi,'(a)')       'yrange auto auto'
1244       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1245       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1246       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1247       if (xpar8.eq.1.) then
1248       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1249       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1250       endif
1251       write(ifhi,'(a)')       'array 2'
1252
1253       sig2=xsigmafit(x)
1254       do i=0,nptg
1255         b=-bmax+2.*real(i)/real(nptg)*bmax
1256         write(ifhi,*) b,exp(-b**2/sig2)
1257       enddo
1258
1259       write(ifhi,'(a)')    '  endarray'
1260       write(ifhi,'(a)')    'closehisto plot 0'
1261
1262
1263
1264       end
1265
1266 c----------------------------------------------------------------------
1267       subroutine xbParD
1268 c----------------------------------------------------------------------
1269
1270       include 'epos.inc'
1271       include 'epos.incsem'
1272       include 'epos.incpar'
1273
1274       double precision x,om5s,xDfit,z(0:200),y
1275 c     & ,omYuncut,omNpuncut,omYcut
1276
1277       nptg=50                  !number of point for the graphs
1278       x=dble(xpar4)                 !value of biniDf (impact parameter)
1279       y=dble(xpar3)                 !value of y (rapidity)
1280       if(x.gt.dexp(-2.d0*y))x=x*dexp(-2.d0*y)
1281       zz=xpar7
1282
1283       bmax=xpar2
1284       bmax=max(0.1,bmax)
1285       t=1.
1286 c      iqqN=0
1287 c      iqq=int(xpar7)
1288       iqq1=-1
1289       iqq2=-1
1290
1291       if (engy.ge.10.) then
1292        if (engy.ge.100.) then
1293         if (engy.ge.1000.) then
1294          if (engy.ge.10000.) then
1295       write(ifhi,'(a,I5)')       'openhisto name DbParam-',int(engy)
1296          else
1297       write(ifhi,'(a,I4)')       'openhisto name DbParam-',int(engy)
1298         endif
1299         else
1300       write(ifhi,'(a,I3)')       'openhisto name DbParam-',int(engy)
1301        endif
1302        else
1303       write(ifhi,'(a,I2)')       'openhisto name DbParam-',int(engy)
1304       endif
1305       else
1306       write(ifhi,'(a,I1)')       'openhisto name DbParam-',int(engy)
1307       endif
1308       write(ifhi,'(a)')       'htyp lbu'
1309       if(xpar5.eq.0.)then
1310         write(ifhi,'(a)')       'xmod lin ymod lin'
1311       else
1312         write(ifhi,'(a)')       'xmod lin ymod log'
1313       endif
1314       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1315       write(ifhi,'(a)')       'yrange auto auto'
1316       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1317       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1318       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1319       if (xpar8.eq.1.) then
1320       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1321       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1322       endif
1323       write(ifhi,'(a)')       'array 2'
1324
1325       xp=sqrt(real(x))*exp(real(y))
1326       xm=sqrt(real(x))*exp(-real(y))
1327       imax=idxD1
1328       if(iomega.eq.2)imax=1
1329
1330       if(xpar6.eq.1.)t=real(xDfit(zz,0,imax,smaxDf,xp,xm,0.))
1331       if(abs(t).lt.1.d-8)t=1.
1332       do i=0,nptg
1333         b=-bmax+2.*real(i)/real(nptg)*bmax
1334          z(i)=xDfit(zz,0,imax,smaxDf,xp,xm,b)/dble(t)
1335          write(ifhi,*) b,z(i)
1336        enddo
1337
1338       write(ifhi,'(a)')    '  endarray'
1339
1340       write(ifhi,'(a)')    'closehisto plot 0-'
1341
1342
1343 c**********************************************************************
1344
1345       if (engy.ge.10.) then
1346         if (engy.ge.100.) then
1347           if (engy.ge.1000.) then
1348             if (engy.ge.10000.) then
1349            write(ifhi,'(a,I5)')  'openhisto name DbParamI-',int(engy)
1350             else
1351            write(ifhi,'(a,I4)')  'openhisto name DbParamI-',int(engy)
1352             endif
1353           else
1354            write(ifhi,'(a,I3)')  'openhisto name DbParamI-',int(engy)
1355           endif
1356         else
1357           write(ifhi,'(a,I2)')  'openhisto name DbParamI-',int(engy)
1358         endif
1359       else
1360         write(ifhi,'(a,I1)')  'openhisto name DbParamI-',int(engy)
1361       endif
1362       write(ifhi,'(a)')       'htyp lga'
1363       write(ifhi,'(a)')       'xmod lin ymod lin'
1364       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1365       write(ifhi,'(a)')       'yrange .01 auto'
1366       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1367       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1368       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1369       if (xpar8.eq.1.) then
1370         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1371         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1372       endif
1373       write(ifhi,'(a)')       'array 2'
1374
1375       if(xpar6.eq.1.)t=(alpDs(idxD0,iclegy,iclpro,icltar)
1376      &   *(smaxDf*real(x))**betDs(idxD0,iclegy,iclpro,icltar)
1377      &                +alpDs(1,iclegy,iclpro,icltar)
1378      &   *(smaxDf*real(x))**betDs(1,iclegy,iclpro,icltar))
1379      &   /real(idxD0+1)
1380       if(abs(t).lt.1.d-8)t=1.
1381       do i=0,nptg
1382         b=-bmax+2.*real(i)/real(nptg)*bmax
1383         tmp=(alpDs(idxD0,iclegy,iclpro,icltar)/t
1384      &   *(smaxDf*real(x))**(betDs(idxD0,iclegy,iclpro,icltar)
1385      &           +gamDs(idxD0,iclegy,iclpro,icltar)*b**2.)
1386      &        *exp(-b**2./delDs(idxD0,iclegy,iclpro,icltar))
1387      &                 +alpDs(1,iclegy,iclpro,icltar)/t
1388      &   *(smaxDf*real(x))**(betDs(1,iclegy,iclpro,icltar)
1389      &           +gamDs(1,iclegy,iclpro,icltar)*b**2.)
1390      &        *exp(-b**2./delDs(1,iclegy,iclpro,icltar)))
1391      &   /real(idxD0+1)
1392         write(ifhi,*) b,tmp
1393       enddo
1394
1395       write(ifhi,'(a)')    '  endarray'
1396       write(ifhi,'(a)')    'closehisto plot 0-'
1397
1398 c**********************************************************************
1399
1400       if (engy.ge.10.) then
1401         if (engy.ge.100.) then
1402           if (engy.ge.1000.) then
1403             if (engy.ge.10000.) then
1404            write(ifhi,'(a,I5)')  'openhisto name DbExaI-',int(engy)
1405             else
1406            write(ifhi,'(a,I4)')  'openhisto name DbExaI-',int(engy)
1407             endif
1408           else
1409            write(ifhi,'(a,I3)')  'openhisto name DbExaI-',int(engy)
1410           endif
1411         else
1412           write(ifhi,'(a,I2)')  'openhisto name DbExaI-',int(engy)
1413         endif
1414       else
1415         write(ifhi,'(a,I1)')  'openhisto name DbExaI-',int(engy)
1416       endif
1417       write(ifhi,'(a)')       'htyp poc'
1418       write(ifhi,'(a)')       'xmod lin ymod lin'
1419       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1420       write(ifhi,'(a)')       'yrange .01 auto'
1421       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1422       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1423       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1424       if (xpar8.eq.1.) then
1425         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1426         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1427       endif
1428       write(ifhi,'(a)')       'array 2'
1429
1430       if(xpar6.eq.1.)t=2*real(om5s(real(x)*smaxDf,x,0.d0,0.,iqq1,iqq2)
1431      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar))))
1432       if(abs(t).lt.1.d-8)t=1.
1433       nptg=nptg/2
1434       do i=0,nptg
1435         b=-bmax+2.*real(i)/real(nptg)*bmax
1436         tmp=2*real(om5s(real(x)*smaxDf,x,y,b,iqq1,iqq2)
1437      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar))))/t
1438         write(ifhi,*) b,tmp
1439       enddo
1440
1441       write(ifhi,'(a)')    '  endarray'
1442
1443       write(ifhi,'(a)')    'closehisto plot 0'
1444
1445       end
1446
1447 c----------------------------------------------------------------------
1448       subroutine xGexaJ
1449 c----------------------------------------------------------------------
1450
1451       include 'epos.inc'
1452       include 'epos.incems'
1453       include 'epos.incsem'
1454       include 'epos.incpar'
1455
1456       parameter(nptg=50) !number of point for the graphs
1457       double precision x(0:nptg),omJ1(0:nptg),xtmp,w(0:nptg),corfa
1458       double precision z(0:nptg),om5J!,omYuncut,omYuncutJ
1459       double precision xminr,y,t,omJ2(0:nptg),omJ3(0:nptg),omJ4(0:nptg)
1460      &,omJ5(0:nptg)!,omJ6(0:nptg),omJ7(0:nptg)
1461
1462
1463       kollsave=koll        !koll modified in zzfz
1464       koll=1
1465       biniDf=xpar2                 !value of biniDf (impact parameter)
1466       y=dble(xpar3)
1467       t=1.d0
1468       xtmp=xmaxDf
1469       xmaxDf=exp(-2.d0*y)
1470
1471       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
1472
1473       nnnmax=1
1474       if(iscreen.ne.0)nnnmax=2
1475       iscreensave=iscreen
1476
1477       do nnn=1,nnnmax
1478
1479         if(nnn.eq.2)then
1480           iscreen=0
1481           corfa=0d0
1482           zzp=0.
1483           zzt=0.
1484         elseif(iscreen.ne.0)then
1485           call zzfz(zzp,zzt,kollth,biniDf)
1486           koll=kollth
1487           rs=r2had(iclpro)+r2had(icltar)+slopom*log(engy*engy)
1488           rpom=4.*.0389*rs
1489           b2x=epscrp*rpom
1490           zzini=epscrw*fscra(engy/egyscr)*exp(-biniDf*biniDf/2./b2x)
1491           zzini=min(zzini,epscrx) !saturation
1492           zzp=zzini+zzp
1493           zzt=zzini+zzt
1494           if(gfactor.lt.0.)then
1495             corfa=dble(abs(gfactor)*(zzp+zzt))
1496           else
1497             corfa=0d0
1498           endif
1499         else
1500           corfa=0d0
1501           zzp=0.
1502           zzt=0.
1503         endif
1504
1505
1506       if (engy.ge.10.) then
1507        if (engy.ge.100.) then
1508         if (engy.ge.1000.) then
1509          if (engy.ge.10000.) then
1510       write(ifhi,'(a,I5)')  'openhisto name GIexa-',int(engy)
1511          else
1512       write(ifhi,'(a,I4)')  'openhisto name GIexa-',int(engy)
1513         endif
1514         else
1515       write(ifhi,'(a,I3)')  'openhisto name GIexa-',int(engy)
1516        endif
1517        else
1518       write(ifhi,'(a,I2)')  'openhisto name GIexa-',int(engy)
1519       endif
1520       else
1521       write(ifhi,'(a,I1)')  'openhisto name GIexa-',int(engy)
1522       endif
1523       write(ifhi,'(a)')       'htyp lru'
1524       if(xpar5.eq.0.)then
1525         write(ifhi,'(a)')       'xmod log ymod log'
1526       else
1527         write(ifhi,'(a)')       'xmod log ymod lin'
1528       endif
1529       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1530       write(ifhi,'(a)')       'yrange auto auto'
1531       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1532       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1533       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1534       if (xpar8.eq.1.) then
1535       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1536       endif
1537       write(ifhi,'(a)')       'array 2'
1538
1539       x(0)=xminr
1540       do i=0,nptg
1541         if (i.ne.0) x(i)=x(0)*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1542         z(i)=0.d0
1543         omJ1(i)=om5J(zzp,zzt,x(i),y,biniDf,0)
1544         omJ2(i)=om5J(zzp,zzt,x(i),y,biniDf,1)
1545         omJ3(i)=om5J(zzp,zzt,x(i),y,biniDf,2)
1546         omJ3(i)=omJ3(i)+om5J(zz,zz,x(i),y,biniDf,3)
1547         omJ4(i)=om5J(zzp,zzt,x(i),y,biniDf,4)
1548         omJ5(i)=om5J(zzp,zzt,x(i),y,biniDf,5)
1549         z(i)=omJ1(i)+omJ2(i)+omJ3(i)+omJ4(i)+omJ5(i)
1550         w(i)=om5J(zzp,zzt,x(i),y,biniDf,-1)
1551         if(xpar6.eq.0)t=z(i)/w(i)
1552         if(xpar6.eq.1)t=z(i)
1553         write(ifhi,*) x(i),z(i)/t
1554       enddo
1555
1556       write(ifhi,'(a)')    '  endarray'
1557       write(ifhi,'(a)')    'closehisto plot 0-'
1558
1559
1560 c*************************************************************************
1561
1562       if (engy.ge.10.) then
1563        if (engy.ge.100.) then
1564         if (engy.ge.1000.) then
1565          if (engy.ge.10000.) then
1566       write(ifhi,'(a,I5)')       'openhisto name GIsoft-',int(engy)
1567          else
1568       write(ifhi,'(a,I4)')       'openhisto name GIsoft-',int(engy)
1569         endif
1570         else
1571       write(ifhi,'(a,I3)')       'openhisto name GIsoft-',int(engy)
1572        endif
1573        else
1574       write(ifhi,'(a,I2)')       'openhisto name GIsoft-',int(engy)
1575       endif
1576       else
1577       write(ifhi,'(a,I1)')       'openhisto name GIsoft-',int(engy)
1578       endif
1579       write(ifhi,'(a)')       'htyp lba'
1580       write(ifhi,'(a)')       'xmod log ymod log'
1581       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1582       write(ifhi,'(a)')       'yrange auto auto'
1583       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1584       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1585       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1586       if (xpar8.eq.1.) then
1587       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1588       endif
1589       write(ifhi,'(a)')       'array 2'
1590
1591       do i=0,nptg
1592         if(xpar6.eq.0)t=z(i)/w(i)
1593         if(xpar6.eq.1)t=z(i)
1594         gfa=1.
1595       if(omJ2(i).gt.0)gfa=exp(-max(0d0,corfa
1596      &                   *(1d0-sqrt(dble(xggfit)*omJ1(i)/omJ2(i)))))
1597         write(ifhi,*) x(i),omJ1(i)*gfa/t
1598       enddo
1599
1600       write(ifhi,'(a)')    '  endarray'
1601       write(ifhi,'(a)')    'closehisto plot 0-'
1602
1603 c*************************************************************************
1604
1605       if (engy.ge.10.) then
1606        if (engy.ge.100.) then
1607         if (engy.ge.1000.) then
1608          if (engy.ge.10000.) then
1609       write(ifhi,'(a,I5)')       'openhisto name GIgg-',int(engy)
1610          else
1611       write(ifhi,'(a,I4)')       'openhisto name GIgg-',int(engy)
1612         endif
1613         else
1614       write(ifhi,'(a,I3)')       'openhisto name GIgg-',int(engy)
1615        endif
1616        else
1617       write(ifhi,'(a,I2)')       'openhisto name GIgg-',int(engy)
1618       endif
1619       else
1620       write(ifhi,'(a,I1)')       'openhisto name GIgg-',int(engy)
1621       endif
1622       write(ifhi,'(a)')       'htyp lra'
1623       write(ifhi,'(a)')       'xmod log ymod log'
1624       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1625       write(ifhi,'(a)')       'yrange auto auto'
1626       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1627       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1628       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1629       if (xpar8.eq.1.) then
1630       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1631       endif
1632       write(ifhi,'(a)')       'array 2'
1633
1634       do i=0,nptg
1635         if(xpar6.eq.0)t=z(i)/w(i)
1636         if(xpar6.eq.1)t=z(i)
1637         gfa=1.
1638         if(omJ2(i).gt.0)gfa=exp(-max(0d0,corfa
1639      &                     *(1d0-sqrt(dble(xggfit)*omJ1(i)/omJ2(i)))))
1640         write(ifhi,*) x(i),(omJ2(i)+(1.-gfa)*omJ1(i))/t
1641       enddo
1642
1643       write(ifhi,'(a)')    '  endarray'
1644       write(ifhi,'(a)')    'closehisto plot 0-'
1645
1646 c*************************************************************************
1647
1648       if (engy.ge.10.) then
1649        if (engy.ge.100.) then
1650         if (engy.ge.1000.) then
1651          if (engy.ge.10000.) then
1652       write(ifhi,'(a,I5)')       'openhisto name GIgq-',int(engy)
1653          else
1654       write(ifhi,'(a,I4)')       'openhisto name GIgq-',int(engy)
1655         endif
1656         else
1657       write(ifhi,'(a,I3)')       'openhisto name GIgq-',int(engy)
1658        endif
1659        else
1660       write(ifhi,'(a,I2)')       'openhisto name GIgq-',int(engy)
1661       endif
1662       else
1663       write(ifhi,'(a,I1)')       'openhisto name GIgq-',int(engy)
1664       endif
1665       write(ifhi,'(a)')       'htyp lga'
1666       write(ifhi,'(a)')       'xmod log ymod log'
1667       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1668       write(ifhi,'(a)')       'yrange auto auto'
1669       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1670       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1671       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1672       if (xpar8.eq.1.) then
1673       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1674       endif
1675       write(ifhi,'(a)')       'array 2'
1676
1677       do i=0,nptg
1678         if(xpar6.eq.0)t=z(i)/w(i)
1679         if(xpar6.eq.1)t=z(i)
1680         write(ifhi,*) x(i),omJ3(i)/t
1681       enddo
1682
1683       write(ifhi,'(a)')    '  endarray'
1684       write(ifhi,'(a)')    'closehisto plot 0-'
1685
1686 c*************************************************************************
1687
1688       if (engy.ge.10.) then
1689        if (engy.ge.100.) then
1690         if (engy.ge.1000.) then
1691          if (engy.ge.10000.) then
1692       write(ifhi,'(a,I5)')       'openhisto name GIqq-',int(engy)
1693          else
1694       write(ifhi,'(a,I4)')       'openhisto name GIqq-',int(engy)
1695         endif
1696         else
1697       write(ifhi,'(a,I3)')       'openhisto name GIqq-',int(engy)
1698        endif
1699        else
1700       write(ifhi,'(a,I2)')       'openhisto name GIqq-',int(engy)
1701       endif
1702       else
1703       write(ifhi,'(a,I1)')       'openhisto name GIqq-',int(engy)
1704       endif
1705       write(ifhi,'(a)')       'htyp lya'
1706       write(ifhi,'(a)')       'xmod log ymod log'
1707       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1708       write(ifhi,'(a)')       'yrange auto auto'
1709       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1710       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1711       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1712       if (xpar8.eq.1.) then
1713       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1714       endif
1715       write(ifhi,'(a)')       'array 2'
1716
1717       do i=0,nptg
1718         if(xpar6.eq.0)t=z(i)/w(i)
1719         if(xpar6.eq.1)t=z(i)
1720         write(ifhi,*) x(i),omJ4(i)/t
1721       enddo
1722
1723       write(ifhi,'(a)')    '  endarray'
1724       write(ifhi,'(a)')    'closehisto plot 0-'
1725
1726 c*************************************************************************
1727
1728       if (engy.ge.10.) then
1729        if (engy.ge.100.) then
1730         if (engy.ge.1000.) then
1731          if (engy.ge.10000.) then
1732       write(ifhi,'(a,I5)')       'openhisto name GIdif-',int(engy)
1733          else
1734       write(ifhi,'(a,I4)')       'openhisto name GIdif-',int(engy)
1735         endif
1736         else
1737       write(ifhi,'(a,I3)')       'openhisto name GIdif-',int(engy)
1738        endif
1739        else
1740       write(ifhi,'(a,I2)')       'openhisto name GIdif-',int(engy)
1741       endif
1742       else
1743       write(ifhi,'(a,I1)')       'openhisto name GIdif-',int(engy)
1744       endif
1745       write(ifhi,'(a)')       'htyp lfa'
1746       write(ifhi,'(a)')       'xmod log ymod log'
1747       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1748       write(ifhi,'(a)')       'yrange auto auto'
1749       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1750       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1751       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1752       if (xpar8.eq.1.) then
1753       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1754       endif
1755       write(ifhi,'(a)')       'array 2'
1756
1757       do i=0,nptg
1758         if(xpar6.eq.0)t=z(i)/w(i)
1759         if(xpar6.eq.1)t=z(i)
1760         write(ifhi,*) x(i),omJ5(i)/t
1761       enddo
1762
1763       write(ifhi,'(a)')    '  endarray'
1764       write(ifhi,'(a)')    'closehisto plot 0-'
1765
1766 c*************************************************************************
1767
1768       if (engy.ge.10.) then
1769        if (engy.ge.100.) then
1770         if (engy.ge.1000.) then
1771          if (engy.ge.10000.) then
1772       write(ifhi,'(a,I5)')       'openhisto name GItot-',int(engy)
1773          else
1774       write(ifhi,'(a,I4)')       'openhisto name GItot-',int(engy)
1775         endif
1776         else
1777       write(ifhi,'(a,I3)')       'openhisto name GItot-',int(engy)
1778        endif
1779        else
1780       write(ifhi,'(a,I2)')       'openhisto name GItot-',int(engy)
1781       endif
1782       else
1783       write(ifhi,'(a,I1)')       'openhisto name GItot-',int(engy)
1784       endif
1785       write(ifhi,'(a)')       'htyp pfc'
1786       write(ifhi,'(a)')       'xmod log ymod log'
1787       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1788       write(ifhi,'(a)')       'yrange auto auto'
1789       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1790       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1791       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1792       if (xpar8.eq.1.) then
1793       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1794       endif
1795       write(ifhi,'(a)')       'array 2'
1796
1797       do i=0,nptg
1798         if(xpar6.eq.0)t=z(i)/w(i)
1799         if(xpar6.eq.1)t=z(i)
1800         write(ifhi,*) x(i),(omJ1(i)+omJ2(i)+omJ3(i)+omJ4(i)+omJ5(i))
1801      &                  /t
1802       enddo
1803
1804       write(ifhi,'(a)')    '  endarray'
1805
1806       write(ifhi,'(a)')    'closehisto plot 0'
1807
1808         if(nnn.eq.2)iscreen=iscreensave
1809
1810       enddo
1811
1812       xmaxDf=xtmp
1813       koll=kollsave
1814
1815       end
1816
1817
1818
1819 c----------------------------------------------------------------------
1820       subroutine xsParD
1821 c----------------------------------------------------------------------
1822
1823       include 'epos.inc'
1824       include 'epos.incsem'
1825       include 'epos.incpar'
1826
1827       double precision x,xminr,y,om5s,xtmp,z
1828 c     & ,t,omYuncut
1829
1830       nptg=50                  !number of point for the graphs
1831       biniDf=xpar2                 !value of biniDf (impact parameter)
1832       y=dble(xpar3)                 !value of y (rapidity)
1833       xtmp=xmaxDf
1834       xmaxDf=dexp(-2.d0*y)
1835       call Class('xsParD     ')
1836       iqq1=-1
1837       iqq2=-1
1838 c      iqqN=0
1839
1840       xminr=dble(egylow/engy)**2.d0  !value of xminr for plotting the function
1841
1842 c**********************************************************************
1843
1844       if (engy.ge.10.) then
1845         if (engy.ge.100.) then
1846           if (engy.ge.1000.) then
1847             if (engy.ge.10000.) then
1848            write(ifhi,'(a,I5)')  'openhisto name DParamI-',int(engy)
1849             else
1850            write(ifhi,'(a,I4)')  'openhisto name DParamI-',int(engy)
1851             endif
1852           else
1853            write(ifhi,'(a,I3)')  'openhisto name DParamI-',int(engy)
1854           endif
1855         else
1856           write(ifhi,'(a,I2)')  'openhisto name DParamI-',int(engy)
1857         endif
1858       else
1859         write(ifhi,'(a,I1)')  'openhisto name DParamI-',int(engy)
1860       endif
1861       write(ifhi,'(a)')       'htyp lga'
1862       if(xpar5.eq.0.)then
1863         write(ifhi,'(a)')       'xmod log ymod log'
1864       else
1865         write(ifhi,'(a)')       'xmod log ymod lin'
1866       endif
1867       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1868       write(ifhi,'(a)')       'yrange auto auto'
1869       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1870       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1871       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1872       if (xpar8.eq.1.) then
1873         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1874       endif
1875       write(ifhi,'(a)')       'array 2'
1876
1877       do i=0,nptg
1878         x=xminr
1879         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1880         egy=sqrt(real(x))*engy
1881         iclegy=1+int((log(egy)-log(egylow))/log(egyfac))
1882         write(ifhi,*) x,(alpDs(idxD0,iclegy,iclpro,icltar)
1883      &   *(smaxDf*real(x))**(betDs(idxD0,iclegy,iclpro,icltar)
1884      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
1885      &     *exp(-biniDf**2./delDs(idxD0,iclegy,iclpro,icltar))
1886      &                 +alpDs(1,iclegy,iclpro,icltar)
1887      &   *(smaxDf*real(x))**(betDs(1,iclegy,iclpro,icltar)
1888      &           +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
1889      &        *exp(-biniDf**2./delDs(1,iclegy,iclpro,icltar)))
1890      &   /real(idxD0+1)
1891       enddo
1892
1893       write(ifhi,'(a)')    '  endarray'
1894       write(ifhi,'(a)')    'closehisto plot 0-'
1895
1896 c**********************************************************************
1897
1898       if (engy.ge.10.) then
1899         if (engy.ge.100.) then
1900           if (engy.ge.1000.) then
1901             if (engy.ge.10000.) then
1902            write(ifhi,'(a,I5)')  'openhisto name DExaI-',int(engy)
1903             else
1904            write(ifhi,'(a,I4)')  'openhisto name DExaI-',int(engy)
1905             endif
1906           else
1907            write(ifhi,'(a,I3)')  'openhisto name DExaI-',int(engy)
1908           endif
1909         else
1910           write(ifhi,'(a,I2)')  'openhisto name DExaI-',int(engy)
1911         endif
1912       else
1913         write(ifhi,'(a,I1)')  'openhisto name DExaI-',int(engy)
1914       endif
1915       write(ifhi,'(a)')       'htyp poc'
1916       write(ifhi,'(a)')       'xmod log ymod log'
1917       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1918       write(ifhi,'(a)')       'yrange .01 auto'
1919       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1920       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1921       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1922       if (xpar8.eq.1.) then
1923         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1924       endif
1925       write(ifhi,'(a)')       'array 2'
1926
1927       nptg=nptg/2
1928       do i=0,nptg
1929         x=xminr
1930         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1931         egy=sqrt(real(x))*engy
1932         iclegy=1+int((log(egy)-log(egylow))/log(egyfac))
1933         z=2d0*om5s(real(x)*smaxDf,1.d0,0.d0,biniDf,iqq1,iqq2)
1934      &              /dble(chad(iclpro)*chad(icltar))
1935       write(ifhi,*) x,z
1936       enddo
1937
1938       write(ifhi,'(a)')    '  endarray'
1939       write(ifhi,'(a)')    'closehisto plot 0'
1940
1941       xmaxDf=xtmp
1942
1943       end
1944
1945 c----------------------------------------------------------------------
1946       subroutine xyParD
1947 c----------------------------------------------------------------------
1948
1949       include 'epos.inc'
1950       include 'epos.incsem'
1951       include 'epos.incpar'
1952
1953       double precision x,om5s,xDfit,z(0:200),ymax,y,t
1954 c     & ,omYuncut,omNpuncut
1955
1956       nptg=50                  !number of point for the graphs
1957       biniDf=xpar2                 !value of biniDf (impact parameter)
1958
1959       x=dble(xpar4)
1960       if(x.le.1.d-20)x=1.d0/dble(engy)
1961       ymax=-.5d0*dlog(x)
1962       zz=xpar7
1963       iqq1=-1
1964       iqq2=-1
1965 c      iqqN=0
1966
1967
1968       if (engy.ge.10.) then
1969        if (engy.ge.100.) then
1970         if (engy.ge.1000.) then
1971          if (engy.ge.10000.) then
1972       write(ifhi,'(a,I5)')       'openhisto name DyParam-',int(engy)
1973          else
1974       write(ifhi,'(a,I4)')       'openhisto name DyParam-',int(engy)
1975         endif
1976         else
1977       write(ifhi,'(a,I3)')       'openhisto name DyParam-',int(engy)
1978        endif
1979        else
1980       write(ifhi,'(a,I2)')       'openhisto name DyParam-',int(engy)
1981       endif
1982       else
1983       write(ifhi,'(a,I1)')       'openhisto name DyParam-',int(engy)
1984       endif
1985       write(ifhi,'(a)')       'htyp lbu'
1986       if(xpar5.eq.0.)then
1987         write(ifhi,'(a)')       'xmod lin ymod lin'
1988       else
1989         write(ifhi,'(a)')       'xmod lin ymod log'
1990       endif
1991       write(ifhi,'(a,2e11.3)')'xrange ',-ymax-1.d0,ymax+1.d0
1992       write(ifhi,'(a)')       'yrange auto auto'
1993       write(ifhi,'(a)')      'text 0 0 "xaxis y"'
1994       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1995       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1996       if (xpar8.eq.1.) then
1997       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1998       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "x=',x,'"'
1999       endif
2000       write(ifhi,'(a)')       'array 2'
2001
2002       imax=idxD1
2003       if(iomega.eq.2)imax=1
2004
2005       do i=0,nptg
2006         y=-ymax+(ymax+ymax)*(dble(i)/dble(nptg))
2007         xp=sqrt(real(x))*exp(real(y))
2008         xm=sqrt(real(x))*exp(-real(y))
2009          z(i)=xDfit(zz,0,imax,engy**2,xp,xm,biniDf)
2010          write(ifhi,*) y,z(i)
2011        enddo
2012
2013       write(ifhi,'(a)')    '  endarray'
2014
2015       write(ifhi,'(a)')    'closehisto plot 0-'
2016
2017
2018 c**********************************************************************
2019
2020       if (engy.ge.10.) then
2021         if (engy.ge.100.) then
2022           if (engy.ge.1000.) then
2023             if (engy.ge.10000.) then
2024            write(ifhi,'(a,I5)')  'openhisto name DyParamI-',int(engy)
2025             else
2026            write(ifhi,'(a,I4)')  'openhisto name DyParamI-',int(engy)
2027             endif
2028           else
2029            write(ifhi,'(a,I3)')  'openhisto name DyParamI-',int(engy)
2030           endif
2031         else
2032           write(ifhi,'(a,I2)')  'openhisto name DyParamI-',int(engy)
2033         endif
2034       else
2035         write(ifhi,'(a,I1)')  'openhisto name DyParamI-',int(engy)
2036       endif
2037       write(ifhi,'(a)')       'htyp lga'
2038       write(ifhi,'(a)')       'xmod lin ymod lin'
2039       write(ifhi,'(a,2e11.3)')'xrange ',-ymax-1.d0,ymax+1.d0
2040       write(ifhi,'(a)')       'yrange auto auto'
2041       write(ifhi,'(a)')      'text 0 0 "xaxis y"'
2042       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
2043       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
2044       if (xpar8.eq.1.) then
2045         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
2046         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "x=',x,'"'
2047       endif
2048       write(ifhi,'(a)')       'array 2'
2049
2050       do i=0,nptg
2051         y=-ymax+(ymax+ymax)*(dble(i)/dble(nptg))
2052         xp=sqrt(real(x))*exp(real(y))
2053         xm=sqrt(real(x))*exp(-real(y))
2054         write(ifhi,*) y,(alpDs(idxD0,iclegy,iclpro,icltar)
2055      &     *smaxDf**(betDs(idxD0,iclegy,iclpro,icltar)
2056      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
2057      &     *xp**(betDps(idxD0,iclegy,iclpro,icltar)
2058      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
2059      &     *xm**(betDpps(idxD0,iclegy,iclpro,icltar)
2060      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
2061      &     *exp(-biniDf**2./delDs(idxD0,iclegy,iclpro,icltar))
2062      &                  +alpDs(1,iclegy,iclpro,icltar)
2063      &     *smaxDf**(betDs(1,iclegy,iclpro,icltar)
2064      &            +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
2065      &     *xp**(betDps(1,iclegy,iclpro,icltar)
2066      &            +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
2067      &     *xm**(betDpps(1,iclegy,iclpro,icltar)
2068      &            +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
2069      &        *exp(-biniDf**2./delDs(1,iclegy,iclpro,icltar)))
2070      &     /real(idxD0+1)
2071       enddo
2072
2073       write(ifhi,'(a)')    '  endarray'
2074       write(ifhi,'(a)')    'closehisto plot 0-'
2075
2076 c**********************************************************************
2077
2078       if (engy.ge.10.) then
2079         if (engy.ge.100.) then
2080           if (engy.ge.1000.) then
2081             if (engy.ge.10000.) then
2082            write(ifhi,'(a,I5)')  'openhisto name DyExaI-',int(engy)
2083             else
2084            write(ifhi,'(a,I4)')  'openhisto name DyExaI-',int(engy)
2085             endif
2086           else
2087            write(ifhi,'(a,I3)')  'openhisto name DyExaI-',int(engy)
2088           endif
2089         else
2090           write(ifhi,'(a,I2)')  'openhisto name DyExaI-',int(engy)
2091         endif
2092       else
2093         write(ifhi,'(a,I1)')  'openhisto name DyExaI-',int(engy)
2094       endif
2095       write(ifhi,'(a)')       'htyp poc'
2096       write(ifhi,'(a)')       'xmod lin ymod lin'
2097       write(ifhi,'(a,2e11.3)')'xrange ',-ymax-1.d0,ymax+1.d0
2098       write(ifhi,'(a)')       'yrange auto auto'
2099       write(ifhi,'(a)')      'text 0 0 "xaxis y"'
2100       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
2101       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
2102       if (xpar8.eq.1.) then
2103         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
2104         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "x=',x,'"'
2105       endif
2106       write(ifhi,'(a)')       'array 2'
2107
2108       nptg=nptg/2
2109       do i=0,nptg
2110         y=-ymax+(ymax+ymax)*(dble(i)/dble(nptg))
2111         t=2*om5s(real(x)*smaxDf,x,y,biniDf,iqq1,iqq2)
2112      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
2113         write(ifhi,*) y,t
2114       enddo
2115
2116       write(ifhi,'(a)')    '  endarray'
2117       write(ifhi,'(a)')    'closehisto plot 0'
2118
2119       end
2120
2121 c----------------------------------------------------------------------
2122       subroutine xParSigma
2123 c----------------------------------------------------------------------
2124
2125       include 'epos.inc'
2126       include 'epos.incpar'
2127
2128       double precision x,w(0:200),z(0:200),xminr,t
2129
2130       nptg=20                  !number of point for the graphs
2131
2132       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
2133
2134       if (engy.ge.10.) then
2135        if (engy.ge.100.) then
2136         if (engy.ge.1000.) then
2137          if (engy.ge.10000.) then
2138       write(ifhi,'(a,I5)')    'openhisto name SigmaReel-',int(engy)
2139          else
2140       write(ifhi,'(a,I4)')    'openhisto name SigmaReel-',int(engy)
2141         endif
2142         else
2143       write(ifhi,'(a,I3)')    'openhisto name SigmaReel-',int(engy)
2144        endif
2145        else
2146       write(ifhi,'(a,I2)')    'openhisto name SigmaReel-',int(engy)
2147       endif
2148       else
2149       write(ifhi,'(a,I1)')    'openhisto name SigmaReel-',int(engy)
2150       endif
2151       write(ifhi,'(a)')       'htyp lru'
2152       write(ifhi,'(a)')       'xmod log ymod lin'
2153       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2154       write(ifhi,'(a)')'yrange auto auto'
2155       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2156       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2157       if (xpar8.eq.1.) then
2158       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2159       endif
2160       write(ifhi,'(a)')       'array 2'
2161
2162       do i=0,nptg
2163         X=xminr
2164         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2165         z(i)=dble(sigma2(X,2))
2166         if(z(i).le.0.) z(i)=0.d0
2167         write(ifhi,'(2e14.6)') X,z(i)
2168       enddo
2169
2170       write(ifhi,'(a)')    '  endarray'
2171       write(ifhi,'(a)')    'closehisto plot 0-'
2172
2173
2174 c**********************************************************************
2175
2176       if (engy.ge.10.) then
2177        if (engy.ge.100.) then
2178         if (engy.ge.1000.) then
2179          if (engy.ge.10000.) then
2180       write(ifhi,'(a,I5)')    'openhisto name SigmaParam-',int(engy)
2181          else
2182       write(ifhi,'(a,I4)')    'openhisto name SigmaParam-',int(engy)
2183         endif
2184         else
2185       write(ifhi,'(a,I3)')    'openhisto name SigmaParam-',int(engy)
2186        endif
2187        else
2188       write(ifhi,'(a,I2)')    'openhisto name SigmaParam-',int(engy)
2189       endif
2190       else
2191       write(ifhi,'(a,I1)')    'openhisto name SigmaParam-',int(engy)
2192       endif
2193       write(ifhi,'(a)')       'htyp pfc'
2194       write(ifhi,'(a)')       'xmod log ymod lin'
2195       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2196       write(ifhi,'(a)')'yrange auto auto'
2197       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2198       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!?param!(X)"'
2199       if (xpar8.eq.1.) then
2200       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2201         endif
2202       write(ifhi,'(a)')       'array 2'
2203
2204       do i=0,nptg
2205         X=xminr
2206         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2207         w(i)=dble(xsigmafit(X))
2208         write(ifhi,'(2e14.6)') X,w(i)
2209       enddo
2210
2211       write(ifhi,'(a)')    '  endarray'
2212       write(ifhi,'(a)')    'closehisto plot 0-'
2213
2214
2215 c**********************************************************************
2216
2217       if (engy.ge.10.) then
2218        if (engy.ge.100.) then
2219         if (engy.ge.1000.) then
2220          if (engy.ge.10000.) then
2221       write(ifhi,'(a,I5)')    'openhisto name SigmaInt-',int(engy)
2222          else
2223       write(ifhi,'(a,I4)')    'openhisto name SigmaInt-',int(engy)
2224         endif
2225         else
2226       write(ifhi,'(a,I3)')    'openhisto name SigmaInt-',int(engy)
2227        endif
2228        else
2229       write(ifhi,'(a,I2)')    'openhisto name SigmaInt-',int(engy)
2230       endif
2231       else
2232       write(ifhi,'(a,I1)')    'openhisto name SigmaInt-',int(engy)
2233       endif
2234       write(ifhi,'(a)')       'htyp lba'
2235       write(ifhi,'(a)')       'xmod log ymod lin'
2236       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2237       write(ifhi,'(a)')'yrange auto auto'
2238       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2239       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!?Int!(X)"'
2240       if (xpar8.eq.1.) then
2241       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2242       endif
2243       write(ifhi,'(a)')       'array 2'
2244
2245       do i=0,nptg
2246         X=xminr
2247         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2248         t=dble(sigma1i(X))
2249         write(ifhi,'(2e14.6)') X,t
2250       enddo
2251
2252       write(ifhi,'(a)')    '  endarray'
2253
2254       if(xpar8.eq.1)then
2255       write(ifhi,'(a)')    'closehisto plot 0-'
2256
2257
2258       if (engy.ge.10.) then
2259        if (engy.ge.100.) then
2260         if (engy.ge.1000.) then
2261          if (engy.ge.10000.) then
2262       write(ifhi,'(a,I5)')    'openhisto name SigmaReelSoft-',int(engy)
2263          else
2264       write(ifhi,'(a,I4)')    'openhisto name SigmaReelSoft-',int(engy)
2265         endif
2266         else
2267       write(ifhi,'(a,I3)')    'openhisto name SigmaReelSoft-',int(engy)
2268        endif
2269        else
2270       write(ifhi,'(a,I2)')    'openhisto name SigmaReelSoft-',int(engy)
2271       endif
2272       else
2273       write(ifhi,'(a,I1)')    'openhisto name SigmaReelSoft-',int(engy)
2274       endif
2275       write(ifhi,'(a)')       'htyp pft'
2276       write(ifhi,'(a)')       'xmod log ymod lin'
2277       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2278       write(ifhi,'(a)')'yrange auto auto'
2279       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2280       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2281       if (xpar8.eq.1.) then
2282       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2283       endif
2284       write(ifhi,'(a)')       'array 2'
2285
2286       do i=0,nptg
2287         X=xminr
2288         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2289         t=dble(sigma2(X,0))
2290         if(t.gt.0.) write(ifhi,'(2e14.6)') X,t
2291       enddo
2292
2293       write(ifhi,'(a)')    '  endarray'
2294       write(ifhi,'(a)')    'closehisto plot 0-'
2295
2296
2297       if (engy.ge.10.) then
2298        if (engy.ge.100.) then
2299         if (engy.ge.1000.) then
2300          if (engy.ge.10000.) then
2301       write(ifhi,'(a,I5)')    'openhisto name SigmaReelSh-',int(engy)
2302          else
2303       write(ifhi,'(a,I4)')    'openhisto name SigmaReelSh-',int(engy)
2304         endif
2305         else
2306       write(ifhi,'(a,I3)')    'openhisto name SigmaReelSh-',int(engy)
2307        endif
2308        else
2309       write(ifhi,'(a,I2)')    'openhisto name SigmaReelSh-',int(engy)
2310       endif
2311       else
2312       write(ifhi,'(a,I1)')    'openhisto name SigmaReelSh-',int(engy)
2313       endif
2314       write(ifhi,'(a)')       'htyp pot'
2315       write(ifhi,'(a)')       'xmod log ymod lin'
2316       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2317       write(ifhi,'(a)')'yrange auto auto'
2318       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2319       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2320       if (xpar8.eq.1.) then
2321       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2322       endif
2323       write(ifhi,'(a)')       'array 2'
2324
2325       do i=0,nptg
2326         X=xminr
2327         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2328         t=dble(sigma2(X,1))
2329         if(t.gt.0.)write(ifhi,'(2e14.6)') X,t
2330       enddo
2331
2332       write(ifhi,'(a)')    '  endarray'
2333
2334       write(ifhi,'(a)')    'closehisto plot 0-'
2335
2336 c**********************************************************************
2337       write(ifhi,'(a)')    'openhisto'
2338       write(ifhi,'(a)')       'htyp lya'
2339       write(ifhi,'(a)')       'xmod log ymod lin'
2340       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2341       write(ifhi,'(a)')'yrange auto auto'
2342       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2343       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2344       if (xpar8.eq.1.) then
2345       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2346       endif
2347       write(ifhi,'(a)')       'array 2'
2348
2349       do i=0,nptg
2350         X=xminr
2351         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2352         t=1.d0/dble(-gamD(0,iclpro,icltar)*log(X*smaxDf)
2353      &     +1./delD(0,iclpro,icltar))
2354         if(t.gt.0.)write(ifhi,'(2e14.6)') X,t
2355       enddo
2356
2357       write(ifhi,'(a)')    '  endarray'
2358
2359
2360       write(ifhi,'(a)')    'closehisto plot 0-'
2361
2362 c**********************************************************************
2363       write(ifhi,'(a)')    'openhisto'
2364       write(ifhi,'(a)')       'htyp lyo'
2365       write(ifhi,'(a)')       'xmod log ymod lin'
2366       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2367       write(ifhi,'(a)')'yrange auto auto'
2368       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2369       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2370       if (xpar8.eq.1.) then
2371       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2372       endif
2373       write(ifhi,'(a)')       'array 2'
2374
2375       do i=0,nptg
2376         X=xminr
2377         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2378         t=1.d0/dble(-gamD(1,iclpro,icltar)*log(X*smaxDf)
2379      *     +1./delD(1,iclpro,icltar))
2380         if(t.gt.0.)write(ifhi,'(2e14.6)') X,t
2381       enddo
2382
2383       write(ifhi,'(a)')    '  endarray'
2384
2385       endif
2386
2387
2388       write(ifhi,'(a)')    'closehisto plot 0'
2389
2390 c**********************************************************************
2391 c**********************************************************************
2392
2393       if (engy.ge.10.) then
2394        if (engy.ge.100.) then
2395         if (engy.ge.1000.) then
2396          if (engy.ge.10000.) then
2397       write(ifhi,'(a,I5)')    'openhisto name SigmaDiff-',int(engy)
2398          else
2399       write(ifhi,'(a,I4)')    'openhisto name SigmaDiff-',int(engy)
2400         endif
2401         else
2402       write(ifhi,'(a,I3)')    'openhisto name SigmaDiff-',int(engy)
2403        endif
2404        else
2405       write(ifhi,'(a,I2)')    'openhisto name SigmaDiff-',int(engy)
2406       endif
2407       else
2408       write(ifhi,'(a,I1)')    'openhisto name SigmaDiff-',int(engy)
2409       endif
2410       write(ifhi,'(a)')       'htyp lin'
2411       write(ifhi,'(a)')       'xmod log ymod lin'
2412       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2413       write(ifhi,'(a,2e11.3)')'yrange auto auto'
2414       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2415       write(ifhi,'(a)')    'text 0 0 "yaxis [D][s]/[s]"'
2416       if (xpar8.eq.1.) then
2417       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.9 "Emax=',engy,' GeV"'
2418       endif
2419       write(ifhi,'(a)')       'array 2'
2420
2421       do i=0,nptg
2422         X=xminr
2423         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2424         t=0.d0
2425         if(w(i).gt.0.d0)t=(z(i)-w(i))/w(i)
2426         if(abs(t).gt.0.15d0) t=dsign(0.15d0,t)
2427         write(ifhi,'(2e14.6)') X,t
2428       enddo
2429
2430       write(ifhi,'(a)')    '  endarray'
2431       write(ifhi,'(a)')    'closehisto plot 0'
2432
2433       end
2434
2435 c----------------------------------------------------------------------
2436       subroutine xParGauss
2437 c----------------------------------------------------------------------
2438
2439       include 'epos.inc'
2440       include 'epos.incpar'
2441
2442       double precision x,om5s,xDfit,y,enh,t!,omNpuncut
2443
2444       nptg=50                  !number of point for the graphs
2445       x=dble(xpar4)                  !value of x (energy)
2446       y=dble(xpar2)                  !value of rapidity
2447       zz=xpar7
2448
2449       if (engy.ge.10.) then
2450        if (engy.ge.100.) then
2451         if (engy.ge.1000.) then
2452          if (engy.ge.10000.) then
2453       write(ifhi,'(a,I5)')    'openhisto name GaussExact-',int(engy)
2454            else
2455       write(ifhi,'(a,I4)')    'openhisto name GaussExact-',int(engy)
2456         endif
2457         else
2458       write(ifhi,'(a,I3)')    'openhisto name GaussExact-',int(engy)
2459        endif
2460        else
2461       write(ifhi,'(a,I2)')    'openhisto name GaussExact-',int(engy)
2462       endif
2463       else
2464       write(ifhi,'(a,I1)')    'openhisto name GaussExact-',int(engy)
2465       endif
2466       write(ifhi,'(a)')       'htyp lru'
2467       write(ifhi,'(a)')       'xmod lin ymod lin'
2468       write(ifhi,'(a,2e11.3)')'xrange',0.,bmaxDf
2469       write(ifhi,'(a)')'yrange auto auto'
2470       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2471       write(ifhi,'(a)')    'text 0 0 "yaxis b*D(x+,x-,s,b)"'
2472 c      if (xpar8.eq.1.) then
2473       write(ifhi,'(a,e7.2,a)')  'text 0.6 0.9 "E=',engy,' GeV"'
2474       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.8 "x=',x,'"'
2475       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.7 "y=',y,'"'
2476 c      endif
2477       write(ifhi,'(a)')       'array 2'
2478
2479       do i=0,nptg
2480         b=bmaxDf*(real(i)/real(nptg))
2481         enh=0.d0
2482         t=dble(b)*(2*om5s(real(x)*smaxDf,x,y,b,-1,-1)
2483      &                +enh)
2484      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
2485        write(ifhi,*) b,t
2486       enddo
2487
2488       write(ifhi,'(a)')    '  endarray'
2489       write(ifhi,'(a)')    'closehisto plot 0-'
2490
2491
2492
2493 c**********************************************************************
2494
2495       if (engy.ge.10.) then
2496        if (engy.ge.100.) then
2497         if (engy.ge.1000.) then
2498          if (engy.ge.10000.) then
2499       write(ifhi,'(a,I5)')    'openhisto name GaussParam-',int(engy)
2500          else
2501       write(ifhi,'(a,I4)')    'openhisto name GaussParam-',int(engy)
2502         endif
2503         else
2504       write(ifhi,'(a,I3)')    'openhisto name GaussParam-',int(engy)
2505        endif
2506        else
2507       write(ifhi,'(a,I2)')    'openhisto name GaussParam-',int(engy)
2508       endif
2509       else
2510       write(ifhi,'(a,I1)')    'openhisto name GaussParam-',int(engy)
2511       endif
2512       write(ifhi,'(a)')       'htyp pfc'
2513       write(ifhi,'(a)')       'xmod log ymod lin'
2514       write(ifhi,'(a,2e11.3)')'xrange',0.,bmaxDf
2515       write(ifhi,'(a)')'yrange auto auto'
2516       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2517       write(ifhi,'(a)')    'text 0 0 "yaxis b*D(x+,x-,s,b)"'
2518 c      if (xpar8.eq.1.) then
2519       write(ifhi,'(a,e7.2,a)')  'text 0.6 0.9 "E=',engy,' GeV"'
2520       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.8 "x=',x,'"'
2521       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.7 "y=',y,'"'
2522 c      endif
2523       write(ifhi,'(a)')       'array 2'
2524
2525       imax=idxD1
2526       if(iomega.eq.2)imax=1
2527
2528       do i=0,nptg
2529         b=bmaxDf*(real(i)/real(nptg))
2530         xp=sqrt(real(x))*exp(real(y))
2531         xm=sqrt(real(x))*exp(-real(y))
2532         t=xDfit(zz,0,imax,engy**2,xp,xm,b)
2533         write(ifhi,'(2e14.6)') b,dble(b)*t
2534       enddo
2535
2536       write(ifhi,'(a)')    '  endarray'
2537       write(ifhi,'(a)')    'closehisto plot 0'
2538
2539       end
2540
2541
2542
2543 c----------------------------------------------------------------------
2544       subroutine xParOmega1
2545 c----------------------------------------------------------------------
2546
2547       include 'epos.inc'
2548       include 'epos.incpar'
2549
2550       double precision x,w(0:200),z(0:200)
2551       double precision yp,om1,xminr,Dsoftshval,t
2552
2553       nptg=50                 !number of point for the graphs
2554       biniDf=xpar2               !value of biniDf (impact parameter)
2555       yp=xpar3                   !value of yp (rapidity)
2556       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
2557
2558       if (engy.ge.10.) then
2559        if (engy.ge.100.) then
2560         if (engy.ge.1000.) then
2561          if (engy.ge.10000.) then
2562       write(ifhi,'(a,I5)')    'openhisto name Om1Exact-',int(engy)
2563          else
2564       write(ifhi,'(a,I4)')    'openhisto name Om1Exact-',int(engy)
2565         endif
2566         else
2567       write(ifhi,'(a,I3)')    'openhisto name Om1Exact-',int(engy)
2568        endif
2569        else
2570       write(ifhi,'(a,I2)')    'openhisto name Om1Exact-',int(engy)
2571       endif
2572       else
2573       write(ifhi,'(a,I1)')    'openhisto name Om1Exact-',int(engy)
2574       endif
2575       write(ifhi,'(a)')       'htyp lru'
2576       write(ifhi,'(a)')       'xmod log ymod log'
2577       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2578       write(ifhi,'(a)')       'yrange auto auto'
2579       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2580       write(ifhi,'(a)')    'text 0 0 "yaxis [h](x,0,b)"'
2581       if (xpar8.eq.1.) then
2582       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2583       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
2584       endif
2585       write(ifhi,'(a)')       'array 2'
2586
2587       do i=0,nptg
2588         x=xminr
2589         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2590         w(i)=Dsoftshval(real(x)*engy**2,x,0.d0,biniDf,0)
2591      &       *(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
2592         write(ifhi,*) x,w(i)
2593       enddo
2594
2595       write(ifhi,'(a)')    '  endarray'
2596       write(ifhi,'(a)')    'closehisto plot 0-'
2597
2598
2599 c**********************************************************************
2600
2601       if (engy.ge.10.) then
2602        if (engy.ge.100.) then
2603         if (engy.ge.1000.) then
2604          if (engy.ge.10000.) then
2605       write(ifhi,'(a,I5)')    'openhisto name om5param-',int(engy)
2606          else
2607       write(ifhi,'(a,I4)')    'openhisto name om5param-',int(engy)
2608         endif
2609         else
2610       write(ifhi,'(a,I3)')    'openhisto name om5param-',int(engy)
2611        endif
2612        else
2613       write(ifhi,'(a,I2)')    'openhisto name om5param-',int(engy)
2614       endif
2615       else
2616       write(ifhi,'(a,I1)')    'openhisto name om5param-',int(engy)
2617       endif
2618       write(ifhi,'(a)')       'htyp pfc'
2619       write(ifhi,'(a)')       'xmod log ymod log'
2620       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2621       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2622       write(ifhi,'(a)')    'yrange auto auto'
2623       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2624       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1!(x,0,b)"'
2625       if (xpar8.eq.1.) then
2626       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2627       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
2628       endif
2629       write(ifhi,'(a)')       'array 2'
2630
2631       do i=0,nptg
2632         x=xminr
2633         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2634         z(i)=om1(x,yp,biniDf)
2635         write(ifhi,*) x,z(i)
2636       enddo
2637
2638       write(ifhi,'(a)')    '  endarray'
2639       write(ifhi,'(a)')    'closehisto plot 0'
2640
2641
2642 c**********************************************************************
2643 c**********************************************************************
2644
2645       if (engy.ge.10.) then
2646        if (engy.ge.100.) then
2647         if (engy.ge.1000.) then
2648          if (engy.ge.10000.) then
2649       write(ifhi,'(a,I5)')    'openhisto name Om1Diff-',int(engy)
2650          else
2651       write(ifhi,'(a,I4)')    'openhisto name Om1Diff-',int(engy)
2652         endif
2653         else
2654       write(ifhi,'(a,I3)')    'openhisto name Om1Diff-',int(engy)
2655        endif
2656        else
2657       write(ifhi,'(a,I2)')    'openhisto name Om1Diff-',int(engy)
2658       endif
2659       else
2660       write(ifhi,'(a,I1)')    'openhisto name Om1Diff-',int(engy)
2661       endif
2662       write(ifhi,'(a)')       'htyp lru'
2663       write(ifhi,'(a)')       'xmod log ymod lin'
2664       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2665       write(ifhi,'(a,2e11.3)')'yrange auto auto'
2666       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
2667       write(ifhi,'(a)')    'text 0 0 "yaxis ([w]?1!-[h])/[h]"'
2668       if (xpar8.eq.1.) then
2669       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2670       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
2671       endif
2672       write(ifhi,'(a)')       'array 2'
2673
2674       do i=0,nptg
2675         x=xminr
2676         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2677          t=z(i)!(z(i)-w(i))
2678 c         if(abs(w(i)).gt.0.)t=t/w(i)
2679 c         if(abs(t).gt.0.15d0) t=dsign(0.15d0,t)
2680          write(ifhi,*) x,t
2681       enddo
2682
2683       write(ifhi,'(a)')    '  endarray'
2684       write(ifhi,'(a)')    'closehisto plot 0'
2685
2686
2687       end
2688
2689 c----------------------------------------------------------------------
2690       subroutine xEpsilon(iii)
2691 c----------------------------------------------------------------------
2692 c iii:  modus (0,1,2)
2693 c----------------------------------------------------------------------
2694       include 'epos.inc'
2695       include 'epos.incsem'
2696       include 'epos.incems'
2697       include 'epos.incpar'
2698
2699       parameter(nxeps=20,nyeps=32)
2700       common/cxeps1/w(0:nxeps,nyeps),y1(nyeps),y2(nyeps)
2701       common/cxeps2/db,b1,b2
2702       common/geom/rmproj,rmtarg,bmax,bkmx
2703       character ch*3
2704       common /psar34/ rrr,rrrm
2705       common /psar41/ rrrp,rrrmp
2706       external pprzz,pttzz
2707
2708       b1=0.03
2709       b2=bkmx*1.2
2710       db=(b2-b1)/nyeps
2711
2712         if(iii.eq.0)then
2713
2714       do j=0,nxeps
2715        do k=1,nyeps
2716         w(j,k)=0
2717        enddo
2718       enddo
2719
2720         elseif(iii.eq.2)then
2721
2722
2723       do nj=1,14
2724        y1(nj)=1e+20
2725        y2(nj)=1e-20
2726        do k=1,nyeps
2727         if(w(0,k).ne.0)then
2728          y1(nj)=min(y1(nj),w(nj,k)/w(0,k))
2729          y2(nj)=max(y2(nj),w(nj,k)/w(0,k))
2730         endif
2731        enddo
2732        if(y1(nj).ge.0)then
2733          y1(nj)=max(y1(nj)*.2,1e-4)
2734        else
2735          y1(nj)=y1(nj)*2.
2736        endif
2737        y2(nj)=min(y2(nj)*5,1e4)
2738       enddo
2739       y2(13)=max(y2(13),y2(14))
2740       y2(14)=max(y2(13),y2(14))
2741       y2(1)=max(y2(1),y2(2))
2742       y2(2)=max(y2(1),y2(2))
2743       y2(5)=max(y2(5),y2(6))
2744       y2(6)=max(y2(5),y2(6))
2745       y2(7)=y2(5)
2746       y2(8)=y2(5)
2747       y2(9)=max(y2(9),y2(10))
2748       y2(10)=max(y2(9),y2(10))
2749       y2(11)=y2(9)
2750       y2(12)=y2(9)
2751       do nj=1,14
2752        if(nj.le.9)write(ifhi,'(a,i1)')'openhisto name xEps',nj
2753        if(nj.gt.9)write(ifhi,'(a,i2)')'openhisto name xEps',nj
2754        ch='lin'
2755        if(nj.eq.7.or.nj.eq.11)ch='lyo'
2756        if(nj.eq.8.or.nj.eq.12)ch='lgo'
2757        write(ifhi,'(a)')     'htyp '//ch//' xmod lin'
2758        if(y1(nj).ge.0.)then
2759          write(ifhi,'(a)')     'ymod log'
2760        else
2761          write(ifhi,'(a)')     'ymod lin'
2762        endif
2763        write(ifhi,'(a,e9.2)')'xrange 0 ',b2
2764        if(nj.eq.1.or.nj.eq.3.or.nj.eq.5.or.nj.eq.9)then
2765        write(ifhi,'(a,2e9.2)')     'yrange ',min(y1(nj),y1(nj+1))
2766      *                                      ,max(y2(nj),y2(nj+1))
2767        else
2768        write(ifhi,'(a,2e9.2)')     'yrange ',y1(nj),y2(nj)
2769        endif
2770        write(ifhi,'(a)')     'text 0 0 "xaxis b"'
2771        if(nj.eq.1) write(ifhi,'(a)')'txt "yaxis [e]?GP/T!(b)"'
2772        if(nj.eq.1) write(ifhi,'(a)')'txt "title soft pro   soft tar"'
2773        if(nj.eq.3) write(ifhi,'(a)')'txt "yaxis [e]?G!(b)"'
2774        if(nj.eq.3) write(ifhi,'(a)')'txt "title diff"'
2775 c       if(nj.eq.3) write(ifhi,'(a)')'txt "title soft   semi"'
2776        if(nj.eq.5) write(ifhi,'(a)')'txt "yaxis [b]?eff!(b)"'
2777        if(nj.eq.5) write(ifhi,'(a)')'txt "title soft pro   soft tar"'
2778        if(nj.eq.9) write(ifhi,'(a)')'txt "yaxis [b]?eff!(b)"'
2779        if(nj.eq.9) write(ifhi,'(a)')'txt "title semi pro   semi tar"'
2780        if(nj.eq.13)write(ifhi,'(a)')'txt "yaxis Z?P/T!"'
2781        write(ifhi,'(a)')       'array 2'
2782        do k=1,nyeps
2783         b=b1+(k-0.5)*db
2784         y=0
2785         if(w(0,k).ne.0)y=w(nj,k)/w(0,k)
2786         write(ifhi,'(2e11.3)')b,y
2787        enddo
2788        write(ifhi,'(a)')    '  endarray'
2789        if(nj.eq.2.or.nj.eq.4.or.nj.eq.8.or.nj.eq.12.or.nj.eq.14
2790      &    .or.nj.eq.16)then
2791         write(ifhi,'(a)')    'closehisto plot 0'
2792        else
2793         write(ifhi,'(a)')    'closehisto plot 0-'
2794        endif
2795       enddo
2796       !----15-16-17---
2797       write(ifhi,'(a)') 'openhisto name xEps15'
2798       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2799       write(ifhi,'(a)') 'xrange 0 10'
2800       write(ifhi,'(a)') 'text 0 0 "xaxis b?0!"'
2801       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(b?0!)"'
2802       write(ifhi,'(a)') 'array 2'
2803       do k=1,10
2804        b=(k-0.5)
2805        y=0
2806        if(w(17,k).ne.0)y=w(15,k)/w(17,k)
2807        write(ifhi,'(2e11.3)')b,y
2808       enddo
2809       write(ifhi,'(a)')  '  endarray'
2810       write(ifhi,'(a)')  'closehisto plot 0-'
2811
2812       write(ifhi,'(a)') 'openhisto name xEps16'
2813       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2814       write(ifhi,'(a)') 'xrange 0 10'
2815       write(ifhi,'(a)') 'text 0 0 "xaxis b?0!"'
2816       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(b?0!)"'
2817       write(ifhi,'(a)') 'array 2'
2818       do k=1,10
2819        b=(k-0.5)
2820        y=0
2821        if(w(17,k).ne.0)y=w(16,k)/w(17,k)
2822        write(ifhi,'(2e11.3)')b,y
2823       enddo
2824       write(ifhi,'(a)')  '  endarray'
2825       write(ifhi,'(a)')  'closehisto plot 0'
2826       !----18-19-20---
2827       kk=2
2828       do k=3,32
2829         if(w(18,k).ne.0)kk=k
2830       enddo
2831       xmx=(kk-1)/31.*0.1*maproj*matarg
2832       write(ifhi,'(a)') 'openhisto name xEps18'
2833       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2834       write(ifhi,'(a,f10.2)') 'xrange 0 ',xmx
2835       write(ifhi,'(a)') 'text 0 0 "xaxis n?Gl!"'
2836       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(n?Gl!)"'
2837       write(ifhi,'(a)') 'array 2'
2838       do k=1,32
2839        x=(k-1.)*0.1*maproj*matarg/(nyeps-1.)
2840        y=0
2841        if(w(20,k).ne.0)y=w(18,k)/w(20,k)
2842        write(ifhi,'(2e11.3)')x,y
2843       enddo
2844       write(ifhi,'(a)')  '  endarray'
2845       write(ifhi,'(a)')  'closehisto plot 0-'
2846
2847       write(ifhi,'(a)') 'openhisto name xEps19'
2848       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2849       write(ifhi,'(a,f10.2)') 'xrange 0 ',xmx
2850       write(ifhi,'(a)') 'text 0 0 "xaxis n?Gl!"'
2851       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(n?Gl!)"'
2852       write(ifhi,'(a)') 'array 2'
2853       do k=1,32
2854        x=(k-1.)*0.1*maproj*matarg/(nyeps-1.)
2855        y=0
2856        if(w(20,k).ne.0)y=w(19,k)/w(20,k)
2857        write(ifhi,'(2e11.3)')x,y
2858       enddo
2859       write(ifhi,'(a)')  '  endarray'
2860       write(ifhi,'(a)')  'closehisto plot 0'
2861
2862         endif
2863
2864       end
2865
2866 c----------------------------------------------------------------------
2867       subroutine xZnucTheo
2868 c----------------------------------------------------------------------
2869 c Theoretical mean Z as a function of impact parameter between proj and targ
2870 c----------------------------------------------------------------------
2871       include 'epos.inc'
2872       include 'epos.incsem'
2873       include 'epos.incems'
2874       include 'epos.incpar'
2875
2876       common/geom/rmproj,rmtarg,bmax,bkmx
2877       common /psar34/ rrr,rrrm
2878       common /psar41/ rrrp,rrrmp
2879       external pprzz,pttzz
2880
2881
2882       rs=r2had(iclpro)+r2had(icltar)+slopom*log(engy**2)
2883       bglaub2=4.*.0389*rs
2884       b2x=epscrp*bglaub2
2885       zzini=epscrw*fscra(engy/egyscr)
2886
2887       if(maproj.gt.1)then
2888         rrrp=radnuc(maproj)/difnuc(maproj)
2889         rrrmp=rrrp+log(9.)
2890       endif
2891       if(matarg.gt.1)then
2892         rrr=radnuc(matarg)/difnuc(matarg)
2893         rrrm=rrr+log(9.)
2894       endif
2895
2896       write(ifhi,'(a)') 'openhisto name xZTheo'
2897       write(ifhi,'(a)') 'htyp lyi xmod lin ymod lin'
2898       write(ifhi,'(a)') 'xrange 0 10'
2899       write(ifhi,'(a)') '- text 0 0 "xaxis b?0!"'
2900       write(ifhi,'(a)') '+ txt "yaxis Z?P!(b?0!)"'
2901       write(ifhi,'(a)') '+ txt "yaxis Z?T!(b?0!)"'
2902       write(ifhi,'(a)') 'array -3'
2903       do k=1,10
2904        b=(k-0.5)
2905        call zzfz(zzp,zzt,kollth,b)
2906        zz=min(epscrw,zzini*exp(-b*b/2./b2x))
2907        write(ifhi,'(2e11.3)')b,zz+zzp,zz+zzt
2908       enddo
2909       write(ifhi,'(a)')  '  endarray'
2910       write(ifhi,'(a)')  'closehisto'
2911       write(ifhi,'(a)')  ' plot -htyp lyi xZTheo+1-',
2912      &                   ' plot -htyp lga xZTheo+2'
2913
2914       end
2915
2916 c----------------------------------------------------------------------
2917       subroutine xParOmegaN
2918 c----------------------------------------------------------------------
2919 c xpar1=engy
2920 c xpar2=b
2921 c xpar4=xremnant
2922 c xpar5: 0=log scale (x dep of om) 1=lin scale (b dep of om)
2923 c----------------------------------------------------------------------
2924
2925       include 'epos.inc'
2926       include 'epos.incpar'
2927       include 'epos.incsem'
2928       include 'epos.incems'
2929       double precision x,w(0:200),z(0:200),xminr,t,ghh
2930      *,xprem,omGam,Womint,Gammapp,WomGamint,omGamk
2931 c     *,yp,SomY,omYuncut,y,xtmp
2932
2933       nptg=30                  !number of point for the graphs
2934       biniDf=xpar2               !value of biniDf (impact parameter)
2935       xprem=dble(xpar4)            !value of x remnant
2936       bmax=3.
2937       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
2938
2939 c**********************************************************************
2940
2941       do i=0,3
2942         b=bmax*(real(i)/real(3))
2943         z(i)=1.d0
2944 c        if(xpar5.eq.0.)z(i)=Gammapp(engy**2.,b,1)
2945       enddo
2946
2947       write(ifhi,'(a)')    'openhisto name Womint-1'
2948       write(ifhi,'(a)')       'htyp lru'
2949       if(xpar5.eq.0.)then
2950         write(ifhi,'(a)')       'xmod lin ymod log'
2951       else
2952         write(ifhi,'(a)')       'xmod lin ymod lin'
2953       endif
2954       write(ifhi,'(a,2e11.3)')'xrange',0.,bmax
2955       write(ifhi,'(a)')    'yrange auto auto'
2956       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2957       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?int!(s,b)"'
2958       if (xpar8.eq.1.) then
2959       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2960       endif
2961       write(ifhi,'(a)')       'array 2'
2962
2963       ierr=0
2964       ghh=0.d0
2965       do i=0,3
2966         b=bmax*(real(i)/real(3))
2967         w(i)=Womint(engy**2.,b)
2968         if(b.eq.biniDf)then
2969           write(*,*)'Womint(',b,',1)=',w(i)
2970           ghh=ghh+w(i)
2971         endif
2972         if(w(i).lt.0.d0.and.ierr.eq.0.and.xpar5.eq.0.)then
2973           write(*,*)'Warning Womint(1)<0 =',w(i)
2974           w(i)=-w(i)
2975           ierr=1
2976           elseif(w(i).lt.0.d0.and.xpar5.eq.0.)then
2977             w(i)=-w(i)
2978           elseif(w(i).ge.0.d0.and.ierr.eq.1.and.xpar5.eq.0.)then
2979             ierr=0
2980             write(*,*)'Warning Womint(1)>0 =',w(i)
2981         endif
2982         write(ifhi,*) b,w(i)/z(i)
2983       enddo
2984
2985
2986       write(ifhi,'(a)')    '  endarray'
2987       write(ifhi,'(a)')    'closehisto plot 0-'
2988
2989 c**********************************************************************
2990
2991       write(ifhi,'(a)')    'openhisto'
2992       write(ifhi,'(a)')       'htyp pfc'
2993       if(xpar5.eq.0.)then
2994         write(ifhi,'(a)')       'xmod lin ymod log'
2995       else
2996         write(ifhi,'(a)')       'xmod lin ymod lin'
2997       endif
2998       write(ifhi,'(a,2e11.3)')'xrange',0.,bmax
2999       write(ifhi,'(a)')    'yrange auto auto'
3000       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
3001       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?int!(s,b)"'
3002       if (xpar8.eq.1.) then
3003       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3004       endif
3005       write(ifhi,'(a)')       'array 2'
3006
3007       do i=0,nptg
3008         b=bmax*(real(i)/real(nptg))
3009         z(i)=1.d0
3010         if(xpar5.eq.0.)then
3011           z(i)=z(i)+dabs(WomGamint(b))
3012         endif
3013       enddo
3014       ierr=0
3015       do i=0,nptg
3016         b=bmax*(real(i)/real(nptg))
3017         w(i)=WomGamint(b)
3018         if(w(i).lt.0.d0.and.ierr.eq.0.and.xpar5.eq.0.)then
3019           write(*,*)'Warning WomGamint(1)<0 =',w(i)
3020           w(i)=-w(i)
3021           ierr=1
3022           elseif(w(i).lt.0.d0.and.xpar5.eq.0.)then
3023             w(i)=-w(i)
3024           elseif(w(i).ge.0.d0.and.ierr.eq.1.and.xpar5.eq.0.)then
3025             ierr=0
3026             write(*,*)'Warning WomGamint(1)>0 =',w(i)
3027         endif
3028         write(ifhi,*) b,w(i)/z(i)
3029       enddo
3030
3031
3032       write(ifhi,'(a)')    '  endarray'
3033       write(ifhi,'(a)')    'closehisto plot 0'
3034       t=Gammapp(engy**2.,biniDf,1)
3035       write(*,*)'--> gamma(',biniDf,')=',ghh,t
3036
3037 c**********************************************************************
3038 c**********************************************************************
3039       do k=1,koll
3040         bk(k)=biniDf
3041       enddo
3042       call GfunPark(0)
3043       call integom1(0)
3044
3045       if (engy.ge.10.) then
3046        if (engy.ge.100.) then
3047         if (engy.ge.1000.) then
3048          if (engy.ge.10000.) then
3049       write(ifhi,'(a,I5)')    'openhisto name xOmNG-',int(engy)
3050          else
3051       write(ifhi,'(a,I4)')    'openhisto name xOmNG-',int(engy)
3052         endif
3053         else
3054       write(ifhi,'(a,I3)')    'openhisto name xOmNG-',int(engy)
3055        endif
3056        else
3057       write(ifhi,'(a,I2)')    'openhisto name xOmNG-',int(engy)
3058       endif
3059       else
3060       write(ifhi,'(a,I1)')    'openhisto name xOmNG-',int(engy)
3061       endif
3062       write(ifhi,'(a)')       'htyp lru'
3063       write(ifhi,'(a)')       'xmod log ymod log'
3064       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3065       write(ifhi,'(a)')    'yrange auto auto'
3066       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3067       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?NPi!(x+rem,x-rem,b)"'
3068       if (xpar8.eq.1.) then
3069       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3070       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3071       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.85 "x+?rem!=',xprem,'"'
3072       endif
3073       write(ifhi,'(a)')       'array 2'
3074
3075       do i=0,nptg
3076         x=xminr
3077         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3078         t=dabs(omGam(xprem,x,biniDf))
3079         write(ifhi,*) x,t
3080       enddo
3081
3082       write(ifhi,'(a)')    '  endarray'
3083       write(ifhi,'(a)')    'closehisto plot 0-'
3084
3085       if (engy.ge.10.) then
3086        if (engy.ge.100.) then
3087         if (engy.ge.1000.) then
3088          if (engy.ge.10000.) then
3089       write(ifhi,'(a,I5)')    'openhisto name xOmNG-',int(engy)
3090          else
3091       write(ifhi,'(a,I4)')    'openhisto name xOmNG-',int(engy)
3092         endif
3093         else
3094       write(ifhi,'(a,I3)')    'openhisto name xOmNG-',int(engy)
3095        endif
3096        else
3097       write(ifhi,'(a,I2)')    'openhisto name xOmNG-',int(engy)
3098       endif
3099       else
3100       write(ifhi,'(a,I1)')    'openhisto name xOmNG-',int(engy)
3101       endif
3102       write(ifhi,'(a)')       'htyp pfc'
3103       write(ifhi,'(a)')       'xmod log ymod log'
3104       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3105       write(ifhi,'(a)')    'yrange auto auto'
3106       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3107       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?NPi!(x+rem,x-rem,b)"'
3108       if (xpar8.eq.1.) then
3109       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3110       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3111       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.85 "x+?rem!=',xprem,'"'
3112       endif
3113       write(ifhi,'(a)')       'array 2'
3114
3115       do i=0,nptg
3116         x=xminr
3117         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3118         t=dabs(omGamk(1,xprem,x))
3119         write(ifhi,*) x,t
3120       enddo
3121
3122       write(ifhi,'(a)')    '  endarray'
3123
3124       write(ifhi,'(a)')    'closehisto plot 0'
3125
3126       end
3127
3128 c----------------------------------------------------------------------
3129       subroutine xParGampp
3130 c----------------------------------------------------------------------
3131
3132       include 'epos.inc'
3133       include 'epos.incpar'
3134
3135       double precision Gammapp,GammaGauss,sg,sgmc,Znorm,Zn,t
3136      *                 ,w(0:200),z(0:200)!,GammaMC
3137
3138       nptg=2                  !number of point for the graphs
3139       biniDf=xpar2              !value of biniDf (impact parameter)
3140
3141
3142 c**************************************************
3143
3144       if (biniDf.lt.1.) then
3145         k=int(10.*biniDf)
3146         write(ifhi,'(a,I1)')  'openhisto name Gamma-b0.',k
3147       else
3148         write(ifhi,'(a,f3.1)')'openhisto name Gamma-b',biniDf
3149       endif
3150       write(ifhi,'(a)')       'htyp lru'
3151       write(ifhi,'(a)')       'xmod lin ymod lin'
3152       write(ifhi,'(a,2e11.3)')'xrange',0.,real(nptg)
3153       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3154       write(ifhi,'(a)')    'text 0 0 "xaxis m"'
3155       write(ifhi,'(a)')    'text 0 0 "yaxis [g]?h1h2!"'
3156       if (xpar8.eq.1.) then
3157         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3158         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3159       endif
3160       write(ifhi,'(a)')       'array 2'
3161       sg=0.d0
3162       do i=0,nptg
3163         w(i)=Gammapp(engy**2.,biniDf,i)
3164         sg=sg+w(i)
3165         write(ifhi,*) i,w(i)
3166       write(*,*) 'G12',i,w(i)
3167       enddo
3168
3169       write(ifhi,'(a)')    '  endarray'
3170       write(ifhi,'(a)')    'closehisto plot 0-'
3171
3172       if (biniDf.lt.1.) then
3173         k=int(10.*biniDf)
3174         write(ifhi,'(a,I1)')  'openhisto name GammaMC-b0.',k
3175       else
3176         write(ifhi,'(a,f3.1)')'openhisto name GammaMC-b',biniDf
3177       endif
3178       write(ifhi,'(a)')       'htyp pfc'
3179       write(ifhi,'(a)')       'xmod lin ymod lin'
3180       write(ifhi,'(a,2e11.3)')'xrange',0.,real(nptg)
3181       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3182       write(ifhi,'(a)')    'text 0 0 "xaxis m"'
3183       write(ifhi,'(a)')    'text 0 0 "yaxis [g]?h1h2!"'
3184       if (xpar8.eq.1.) then
3185         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3186         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3187       endif
3188       write(ifhi,'(a)')       'array 2'
3189
3190       sgmc=0.d0
3191       do i=0,nptg
3192         z(i)=GammaGauss(engy**2.,biniDf,i)
3193         sgmc=sgmc+z(i)
3194         write(ifhi,*) i,z(i)
3195         write(*,*) 'G12gauss',i,z(i)
3196       enddo
3197
3198       write(ifhi,'(a)')    '  endarray'
3199       write(ifhi,'(a)')    'closehisto plot 0'
3200
3201 c**********************************************************************
3202 c**********************************************************************
3203       Zn=Znorm(engy**2,biniDf)
3204
3205       write(ifhi,'(a)')    'openhisto name GammaDiff'
3206       write(ifhi,'(a)')       'htyp lru'
3207       write(ifhi,'(a)')       'xmod lin ymod lin'
3208       write(ifhi,'(a,2e11.3)')'xrange',0.,real(nptg)
3209       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3210       write(ifhi,'(a)')    'text 0 0 "xaxis m"'
3211       write(ifhi,'(a)')    'text 0 0 "yaxis (G-GMC)/G"'
3212       if (xpar8.eq.1.) then
3213       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3214       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
3215       endif
3216       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "[S]?Guncut!=',sg,'"'
3217       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "[S]?Gcut!=',sgmc,'"'
3218       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.7 "Z=',Zn,'"'
3219       write(ifhi,'(a)')       'array 2'
3220
3221       do i=0,nptg
3222         if(w(i).ne.0d0) t=(z(i)-w(i))/w(i)
3223 c         if(abs(t).gt.0.5d0) t=dsign(0.5d0,t)
3224          write(ifhi,*) i,t
3225       enddo
3226
3227       write(ifhi,'(a)')    '  endarray'
3228       write(ifhi,'(a)')    'closehisto plot 0'
3229
3230       end
3231
3232 c----------------------------------------------------------------------
3233       subroutine xParPomInc
3234 c----------------------------------------------------------------------
3235
3236       include 'epos.inc'
3237       include 'epos.incpar'
3238
3239       double precision x,PomIncExact,PomIncUnit,xminr,xm,t
3240      *                 ,w(0:200),z(0:200)
3241
3242       nptg=10                  !number of point for the graphs
3243       biniDf=xpar2              !value of biniDf (impact parameter)
3244       xm=dble(xpar4)              !value of biniDf (impact parameter)
3245       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3246
3247
3248 c********************* red = PomIncXExact *****************************
3249
3250       if (biniDf.lt.1.) then
3251         k=int(10.*biniDf)
3252         write(ifhi,'(a,I1)')  'openhisto name PomIncExact-b0.',k
3253       else
3254         write(ifhi,'(a,f3.1)')'openhisto name PomIncExact-b',biniDf
3255       endif
3256       write(ifhi,'(a)')       'htyp lru'
3257       write(ifhi,'(a)')       'xmod log ymod lin'
3258       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3259       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3260       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3261       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+/dx-"'
3262       if (xpar8.eq.1.) then
3263         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3264         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3265         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x-=',xm,'"'
3266       endif
3267       write(ifhi,'(a)')       'array 2'
3268
3269       do i=0,nptg
3270         x=xminr
3271         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3272 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3273         w(i)=PomIncExact(dsqrt(x),dsqrt(x),biniDf)
3274         write(ifhi,*) x,w(i)
3275       write(*,*) 'Xe',i,w(i)
3276       enddo
3277
3278       write(ifhi,'(a)')    '  endarray'
3279       write(ifhi,'(a)')    'closehisto plot 0-'
3280
3281
3282
3283 c************************* dot = PomIncXUnit **************************
3284 c      nptg=50     !number of point for the graphs
3285
3286       if (biniDf.lt.1.) then
3287         k=int(10.*biniDf)
3288         write(ifhi,'(a,I1)')  'openhisto name PomIncUnit-b0.',k
3289       else
3290         write(ifhi,'(a,f3.1)')'openhisto name PomIncUnit-b',biniDf
3291       endif
3292       write(ifhi,'(a)')       'htyp pfc'
3293       write(ifhi,'(a)')       'xmod log ymod lin'
3294       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3295       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3296       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3297       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+/dx-"'
3298       if (xpar8.eq.1.) then
3299         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3300         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3301         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x-=',xm,' fm"'
3302       endif
3303       write(ifhi,'(a)')       'array 2'
3304
3305       do i=0,nptg
3306         x=xminr
3307         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3308 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3309         z(i)=PomIncUnit(dsqrt(x),dsqrt(x),biniDf)
3310         write(ifhi,*) x,z(i)
3311         write(*,*) 'Xu',i,z(i)
3312       enddo
3313
3314       write(ifhi,'(a)')    '  endarray'
3315       write(ifhi,'(a)')    'closehisto plot 0'
3316
3317 c**********************************************************************
3318 c**********************************************************************
3319
3320       write(ifhi,'(a)')    'openhisto name PomIncDiff'
3321       write(ifhi,'(a)')       'htyp lru'
3322       write(ifhi,'(a)')       'xmod log ymod lin'
3323       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3324       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3325       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3326       write(ifhi,'(a)')    'text 0 0 "yaxis ([w]?5!-G)/G"'
3327       if (xpar8.eq.1.) then
3328       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3329       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3330       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x-=',xm,'"'
3331       endif
3332       write(ifhi,'(a)')       'array 2'
3333
3334       do i=0,nptg
3335         x=xminr
3336         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3337 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3338         t=0.d0
3339         if(w(i).ne.0d0) t=(z(i)-w(i))/w(i)
3340 c         if(abs(t).gt.0.5d0) t=dsign(0.5d0,t)
3341          write(ifhi,*) x,t
3342       enddo
3343
3344       write(ifhi,'(a)')    '  endarray'
3345       write(ifhi,'(a)')    'closehisto plot 0'
3346
3347       end
3348
3349 c----------------------------------------------------------------------
3350       subroutine xParPomIncX
3351 c----------------------------------------------------------------------
3352
3353       include 'epos.inc'
3354       include 'epos.incpar'
3355
3356       double precision x,PomIncXExact,PomIncXUnit,xminr,y
3357
3358       nptg=20                   !number of point for the graphs
3359       biniDf=xpar2              !value of biniDf (impact parameter)
3360       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3361
3362
3363 c********************* red = PomIncXExact *****************************
3364
3365       if (biniDf.lt.1.) then
3366         k=int(10.*biniDf)
3367         write(ifhi,'(a,I1)')  'openhisto name PomIncXExact-b0.',k
3368       else
3369         write(ifhi,'(a,f3.1)')'openhisto name PomIncXExact-b',biniDf
3370       endif
3371       write(ifhi,'(a)')       'htyp lru'
3372       write(ifhi,'(a)')       'xmod log ymod lin'
3373       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3374       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3375       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3376       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx(x,b)"'
3377       if (xpar8.eq.1.) then
3378         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3379         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3380       endif
3381       write(ifhi,'(a)')       'array 2'
3382
3383       do i=0,nptg
3384         x=xminr
3385         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3386 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3387         y=PomIncXExact(x,biniDf)
3388         write(ifhi,*) x,y
3389 c      write(*,*) 'Xe',i
3390       enddo
3391
3392       write(ifhi,'(a)')    '  endarray'
3393       write(ifhi,'(a)')    'closehisto plot 0-'
3394
3395
3396
3397 c************************* dot = PomIncXUnit **************************
3398 c      nptg=50     !number of point for the graphs
3399
3400       if (biniDf.lt.1.) then
3401         k=int(10.*biniDf)
3402         write(ifhi,'(a,I1)')  'openhisto name PomIncXUnit-b0.',k
3403       else
3404         write(ifhi,'(a,f3.1)')'openhisto name PomIncXUnit-b',biniDf
3405       endif
3406       write(ifhi,'(a)')       'htyp pfc'
3407       write(ifhi,'(a)')       'xmod log ymod lin'
3408       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3409       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3410       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3411       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx(x,b)"'
3412       if (xpar8.eq.1.) then
3413         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3414         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3415       endif
3416       write(ifhi,'(a)')       'array 2'
3417
3418       do i=0,nptg
3419         x=xminr
3420         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3421 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3422         write(ifhi,*) x,PomIncXUnit(x,biniDf)
3423         write(*,*) 'Xu',i
3424       enddo
3425
3426       write(ifhi,'(a)')    '  endarray'
3427       write(ifhi,'(a)')    'closehisto plot 0'
3428
3429       end
3430
3431 c----------------------------------------------------------------------
3432       subroutine xParPomIncXI
3433 c----------------------------------------------------------------------
3434
3435       include 'epos.inc'
3436       include 'epos.incpar'
3437
3438       double precision x,xminr
3439       double precision PomIncXIExact,PomIncXIUnit
3440
3441       nptg=20                   !number of point for the graphs
3442       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3443
3444 c*********************red = PomIncXIExact *****************************
3445
3446       write(ifhi,'(a)')       'openhisto name PomIncXIExact'
3447       write(ifhi,'(a)')       'htyp lru'
3448       write(ifhi,'(a)')       'xmod log ymod log'
3449       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3450       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3451       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx(x)"'
3452       if (xpar8.eq.1.) then
3453         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3454       endif
3455       write(ifhi,'(a)')       'array 2'
3456
3457       do i=0,nptg
3458         x=xminr
3459         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3460 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3461         write(ifhi,*) x,PomIncXIExact(x)
3462 c.......write(*,*) 'XIe',i
3463       enddo
3464
3465       write(ifhi,'(a)')    '  endarray'
3466       write(ifhi,'(a)')    'closehisto plot 0-'
3467
3468 c***************************dot = PomIncXIUnit ************************
3469 c.....nptg=50     !number of point for the graphs
3470
3471       write(ifhi,'(a)')       'openhisto name PomIncXIUnit'
3472       write(ifhi,'(a)')       'htyp pfc'
3473       write(ifhi,'(a)')       'xmod log ymod log'
3474       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3475       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3476       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx(x)"'
3477       if (xpar8.eq.1.) then
3478         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3479       endif
3480       write(ifhi,'(a)')       'array 2'
3481
3482       do i=0,nptg
3483         x=xminr
3484         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3485 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3486         write(ifhi,*) x,PomIncXIUnit(x)
3487         write(*,*) 'XIu',i
3488       enddo
3489
3490       write(ifhi,'(a)')    '  endarray'
3491       write(ifhi,'(a)')    'closehisto plot 0'
3492
3493       end
3494
3495 c----------------------------------------------------------------------
3496       subroutine xParPomIncP
3497 c----------------------------------------------------------------------
3498
3499       include 'epos.inc'
3500       include 'epos.incpar'
3501
3502       double precision x,PomIncPUnit,xminr
3503       double precision PomIncPExact
3504
3505       nptg=30                  !number of point for the graphs
3506       biniDf=xpar1              !value of biniDf (impact parameter)
3507       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3508
3509 c*********************red = PomIncPExact *****************************
3510
3511       if (biniDf.lt.1.) then
3512         k=int(10.*biniDf)
3513         write(ifhi,'(a,I1)')  'openhisto name PomIncPExact-b0.',k
3514       else
3515         write(ifhi,'(a,f3.1)')'openhisto name PomIncPExact-b',biniDf
3516       endif
3517       write(ifhi,'(a)')       'htyp lru'
3518       write(ifhi,'(a)')       'xmod log ymod log'
3519       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3520       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3521       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+(x+,b)"'
3522       if (xpar8.eq.1.) then
3523         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3524         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3525       endif
3526       write(ifhi,'(a)')       'array 2'
3527
3528       do i=0,nptg
3529         x=xminr
3530         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3531         write(ifhi,*) x,PomIncPExact(x,biniDf)
3532       enddo
3533
3534       write(ifhi,'(a)')    '  endarray'
3535       write(ifhi,'(a)')    'closehisto plot 0-'
3536
3537
3538 c**************************dot = PomIncPUnit **************************
3539 c.....nptg=50     !number of point for the graphs
3540
3541       if (biniDf.lt.1.) then
3542         k=int(10.*biniDf)
3543         write(ifhi,'(a,I1)')  'openhisto name PomIncPUnit-b0.',k
3544       else
3545         write(ifhi,'(a,f3.1)')'openhisto name PomIncPUnit-b',biniDf
3546       endif
3547       write(ifhi,'(a)')       'htyp pfc'
3548       write(ifhi,'(a)')       'xmod log ymod log'
3549       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3550       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3551       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+(x+,b)"'
3552       if (xpar8.eq.1.) then
3553         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3554         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3555       endif
3556       write(ifhi,'(a)')       'array 2'
3557
3558       do i=0,nptg
3559         x=xminr
3560         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3561         write(*,*) 'Pu',i
3562         write(ifhi,*) x,PomIncPUnit(x,biniDf)
3563       enddo
3564
3565       write(ifhi,'(a)')    '  endarray'
3566       write(ifhi,'(a)')    'closehisto plot 0'
3567
3568       end
3569
3570 c----------------------------------------------------------------------
3571       subroutine xParPomIncPI
3572 c----------------------------------------------------------------------
3573
3574       include 'epos.inc'
3575       include 'epos.incpar'
3576
3577       double precision x,xminr
3578       double precision PomIncPIExact,PomIncPIUnit
3579
3580       nptg=50                  !number of point for the graphs
3581       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3582
3583 c*********************red = PomIncPIExact *****************************
3584 c.....nptg=100     !number of point for the graphs
3585
3586       write(ifhi,'(a)')       'openhisto name PomIncPIExact'
3587       write(ifhi,'(a)')       'htyp lru'
3588       write(ifhi,'(a)')       'xmod log ymod log'
3589       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3590       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3591       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx+(x+)"'
3592       if (xpar8.eq.1.) then
3593         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3594       endif
3595       write(ifhi,'(a)')       'array 2'
3596
3597       do i=0,nptg
3598         x=xminr
3599         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3600         write(ifhi,*) x,PomIncPIExact(x)
3601       enddo
3602
3603       write(ifhi,'(a)')    '  endarray'
3604       write(ifhi,'(a)')    'closehisto plot 0-'
3605
3606 c***************************dot = PomIncPIUnit ************************
3607 c.....nptg=10     !number of point for the graphs
3608
3609       write(ifhi,'(a)')       'openhisto name PomIncPIUnit'
3610       write(ifhi,'(a)')       'htyp pfc'
3611       write(ifhi,'(a)')       'xmod log ymod log'
3612       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3613       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3614       write(ifhi,'(a)')    'text 0 0 "yaxis n?Pom!/dx+(x+)"'
3615       if (xpar8.eq.1.) then
3616         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3617       endif
3618       write(ifhi,'(a)')       'array 2'
3619
3620       do i=0,nptg
3621         x=xminr
3622         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3623         write(ifhi,*) x,PomIncPIUnit(x)
3624         write(*,*) 'PIu',i
3625       enddo
3626
3627       write(ifhi,'(a)')    '  endarray'
3628       write(ifhi,'(a)')    'closehisto plot 0'
3629
3630       end
3631
3632 c----------------------------------------------------------------------
3633       subroutine xParPomIncM
3634 c----------------------------------------------------------------------
3635
3636       include 'epos.inc'
3637       include 'epos.incpar'
3638
3639       double precision x,xminr,PomIncMUnit,PomIncMExact
3640
3641       nptg=50                  !number of point for the graphs
3642       biniDf=xpar1              !value of biniDf (impact parameter)
3643       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3644
3645 c**********************red = PomIncMExact *****************************
3646
3647       if (biniDf.lt.1.) then
3648         k=int(10.*biniDf)
3649         write(ifhi,'(a,I1)')  'openhisto name PomIncMExact-b0.',k
3650       else
3651         write(ifhi,'(a,f3.1)')'openhisto name PomIncMExact-b',biniDf
3652       endif
3653       write(ifhi,'(a)')       'htyp lru'
3654       write(ifhi,'(a)')       'xmod log ymod log'
3655       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3656       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3657       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx-(x-,b)"'
3658       if (xpar8.eq.1.) then
3659         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3660         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3661       endif
3662       write(ifhi,'(a)')       'array 2'
3663
3664       do i=0,nptg
3665         x=xminr
3666         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3667         write(ifhi,*) x,PomIncMExact(x,biniDf)
3668       enddo
3669
3670       write(ifhi,'(a)')    '  endarray'
3671       write(ifhi,'(a)')    'closehisto plot 0-'
3672
3673
3674 c**************************dot = PomIncMUnit **************************
3675 c.....nptg=100     !number of point for the graphs
3676
3677       if (biniDf.lt.1.) then
3678         k=int(10.*biniDf)
3679         write(ifhi,'(a,I1)')  'openhisto name PomIncMUnit-b0.',k
3680       else
3681         write(ifhi,'(a,f3.1)')'openhisto name PomIncMUnit-b',biniDf
3682       endif
3683       write(ifhi,'(a)')       'htyp pfc'
3684       write(ifhi,'(a)')       'xmod log ymod log'
3685       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3686       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3687       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx-(x-,b)"'
3688       if (xpar8.eq.1.) then
3689         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3690         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3691       endif
3692       write(ifhi,'(a)')       'array 2'
3693
3694       do i=0,nptg
3695         x=xminr
3696         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3697         write(ifhi,*) x,PomIncMUnit(x,biniDf)
3698         write(*,*) 'Mu',i
3699       enddo
3700
3701       write(ifhi,'(a)')    '  endarray'
3702       write(ifhi,'(a)')    'closehisto plot 0'
3703
3704       end
3705
3706 c----------------------------------------------------------------------
3707       subroutine xParPomIncMI
3708 c----------------------------------------------------------------------
3709
3710       include 'epos.inc'
3711       include 'epos.incpar'
3712
3713       double precision x,xminr
3714       double precision PomIncMIExact,PomIncMIUnit
3715
3716       nptg=30                 !number of point for the graphs
3717       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3718
3719 c*********************red = PomIncMIExact *****************************
3720 c.....nptg=100     !number of point for the graphs
3721
3722       write(ifhi,'(a)')       'openhisto name PomIncMIExact'
3723       write(ifhi,'(a)')       'htyp lru'
3724       write(ifhi,'(a)')       'xmod log ymod log'
3725       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3726       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3727       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx-(x-)"'
3728       if (xpar8.eq.1.) then
3729         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3730       endif
3731       write(ifhi,'(a)')       'array 2'
3732
3733       do i=0,nptg
3734         x=xminr
3735         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3736         write(ifhi,*) x,PomIncMIExact(x)
3737       enddo
3738
3739       write(ifhi,'(a)')    '  endarray'
3740       write(ifhi,'(a)')    'closehisto plot 0-'
3741
3742 c***************************dot = PomIncMIUnit ************************
3743 c.....nptg=100     !number of point for the graphs
3744
3745       write(ifhi,'(a)')       'openhisto name PomIncMIUnit'
3746       write(ifhi,'(a)')       'htyp pfc'
3747       write(ifhi,'(a)')       'xmod log ymod log'
3748       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3749       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3750       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx-(x-)"'
3751       if (xpar8.eq.1.) then
3752         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3753       endif
3754       write(ifhi,'(a)')       'array 2'
3755
3756       do i=0,nptg
3757         x=xminr
3758         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3759         write(ifhi,*) x,PomIncMIUnit(x)
3760         write(*,*) 'MIu',i
3761       enddo
3762
3763       write(ifhi,'(a)')    '  endarray'
3764       write(ifhi,'(a)')    'closehisto plot 0'
3765
3766       end
3767
3768 c----------------------------------------------------------------------
3769       subroutine xParPomIncJ
3770 c----------------------------------------------------------------------
3771
3772       include 'epos.inc'
3773       include 'epos.incpar'
3774       common/geom/rmproj,rmtarg,bmax,bkmx
3775       parameter(nbib=48)
3776
3777       double precision PomIncJExact,PomIncJUnit
3778
3779
3780       b1=0
3781       b2=bkmx*1.2
3782       db=(b2-b1)/nbib
3783
3784 c*************************red = PomIncJExact **************************
3785
3786       write(ifhi,'(a)')       'openhisto name PomIncJExact'
3787       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
3788       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
3789       write(ifhi,'(a)')    'text 0 0 "yaxis  n?Pom!(b)"'
3790       if (xpar8.eq.1.) then
3791         write(ifhi,'(a,e7.2,a)')'text 0.5 0.8 "s=',engy**2,' GeV^2!"'
3792       endif
3793       write(ifhi,'(a)')       'array 2'
3794       do k=1,nbib
3795         b=b1+(k-0.5)*db
3796         write(ifhi,*)b,PomIncJExact(b)
3797       enddo
3798       write(ifhi,'(a)')    '  endarray'
3799       write(ifhi,'(a)')    'closehisto plot 0-'
3800
3801 c****************************dot = PomIncJUnit ***********************
3802
3803       write(ifhi,'(a)')       'openhisto name PomIncJUnit'
3804       write(ifhi,'(a)')       'htyp pfc xmod lin ymod lin'
3805       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
3806       write(ifhi,'(a)')    'text 0 0 "yaxis  n?Pom!(b)"'
3807       if (xpar8.eq.1.) then
3808         write(ifhi,'(a,e7.2,a)')'text 0.5 0.8 "s=',engy**2,' GeV^2!"'
3809       endif
3810       write(ifhi,'(a)')       'array 2'
3811       do k=1,nbib
3812         b=b1+(k-0.5)*db
3813         write(ifhi,*)b,PomIncJUnit(b)
3814         write(*,*) 'Ju',k
3815       enddo
3816       write(ifhi,'(a)')    '  endarray'
3817       write(ifhi,'(a)')    'closehisto plot 0'
3818
3819
3820       end
3821
3822 c----------------------------------------------------------------------
3823       subroutine xParPhi1
3824 c----------------------------------------------------------------------
3825
3826       include 'epos.inc'
3827       include 'epos.incems'
3828       include 'epos.incsem'
3829       include 'epos.incpar'
3830       double precision x,xminr,y
3831       double precision PhiExpo
3832       double precision PhiExact
3833
3834       nptg=30                  !number of point for the graphs
3835       biniDf=xpar2              !value of biniDf (impact parameter)
3836       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3837       zz=xpar7
3838
3839 c************************* red = PhiMExact ***************************
3840
3841       write(ifhi,'(a)')       'openhisto name Phi1Exact'
3842       write(ifhi,'(a)')       'htyp lru'
3843       write(ifhi,'(a)')      'xmod log ymod lin'
3844       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3845       write(ifhi,'(a)')       'yrange auto auto'
3846       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3847       write(ifhi,'(a)') 'txt "yaxis  [F](x)/x^[a]!"'
3848       write(ifhi,'(a,i4,a,f4.1,a)')
3849      * 'txt  "title E=',nint(engy),' b=',biniDf,'"'
3850       write(ifhi,'(a)')       'array 2'
3851
3852       do i=0,nptg
3853         x=xminr
3854         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3855 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3856         y=0.d0
3857         if(engy**2..lt.5.e06)
3858      &  y=Phiexact(zz,zz,.5,dsqrt(x),dsqrt(x),engy**2,biniDf)
3859      &       *dsqrt(x)**dble(-alplea(iclpro))
3860      &       *dsqrt(x)**dble(-alplea(icltar))
3861         write(ifhi,*)x,y
3862       enddo
3863
3864       write(ifhi,'(a)')    '  endarray'
3865       write(ifhi,'(a)')    'closehisto plot 0-'
3866
3867 c********************** blue = PhiMExpo ******************************
3868
3869       write(ifhi,'(a)')       'openhisto name Phi1Expo'
3870       write(ifhi,'(a)')       'htyp lba'
3871       write(ifhi,'(a)')       'xmod log ymod lin'
3872       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3873       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3874        write(ifhi,'(a)')       'array 2'
3875
3876       do i=0,nptg
3877         x=xminr
3878         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3879 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3880         y=Phiexpo(zz,zz,.5,dsqrt(x),dsqrt(x),engy**2,biniDf)
3881      &       *dsqrt(x)**dble(-alplea(iclpro))
3882      &       *dsqrt(x)**dble(-alplea(icltar))
3883         write(ifhi,*) x,y
3884       enddo
3885
3886       write(ifhi,'(a)')    '  endarray'
3887       write(ifhi,'(a)')    'closehisto plot 0'
3888
3889       end
3890
3891
3892 cc----------------------------------------------------------------------
3893 c      subroutine xParPhi2
3894 cc----------------------------------------------------------------------
3895 c
3896 c      include 'epos.inc'
3897 c      include 'epos.incems'
3898 c      include 'epos.incsem'
3899 c      include 'epos.incpar'
3900 c      double precision x,xminr,xm,y,u(0:100),v(0:100),w(0:100)!,z
3901 c      double precision PhiExpo,omGam,PhiExpoK
3902 c      double precision PhiExact
3903 c
3904 c      nptg=30                  !number of point for the graphs
3905 c      biniDf=xpar2              !value of biniDf (impact parameter)
3906 c      xminr=1.d-3 !/dble(engy**2)  !value of xminr for plotting the function
3907 c      xm=dble(xpar4)
3908 c
3909 cc************************* yellow = PhiMExact ***************************
3910 c
3911 c      write(ifhi,'(a)')       'openhisto name Phi1Exact'
3912 c      write(ifhi,'(a)')       'htyp lru'
3913 c      write(ifhi,'(a)')      'xmod lin ymod lin'
3914 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3915 c      write(ifhi,'(a)')       'yrange auto auto'
3916 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3917 c      write(ifhi,'(a)') 'txt "yaxis  [F](x)/x^[a]!"'
3918 c      write(ifhi,'(a,i4,a,f4.1,a)')
3919 c     * 'txt  "title E=',nint(engy),' b=',biniDf,'"'
3920 c      write(ifhi,'(a)')       'array 2'
3921 c
3922 c      do i=0,nptg
3923 c       ! x=xminr
3924 c       ! if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3925 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3926 c        y=0.d0
3927 c        if(engy**2..lt.5.e06)
3928 c     &  y=Phiexact(0.,0.,1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
3929 c    ! &       *dsqrt(x)**dble(-alplea(iclpro))
3930 c    ! &       *dsqrt(x)**dble(-alplea(icltar))
3931 c        write(ifhi,*)x,y
3932 c      enddo
3933 c
3934 c      write(ifhi,'(a)')    '  endarray'
3935 c      write(ifhi,'(a)')    'closehisto plot 0-'
3936 c
3937 cc********************** blue = PhiMExpo ******************************
3938 c
3939 c      write(ifhi,'(a)')       'openhisto name Phi1Expo'
3940 c      write(ifhi,'(a)')       'htyp lbu'
3941 c      write(ifhi,'(a)')       'xmod lin ymod lin'
3942 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3943 c      write(ifhi,'(a,2e11.3)')'yrange auto auto'
3944 c       write(ifhi,'(a)')       'array 2'
3945 c
3946 c      do i=0,nptg
3947 c       ! x=xminr
3948 c       ! if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3949 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3950 c        y=Phiexpo(0.,0.,1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
3951 c  !   &       *dsqrt(x)**dble(-alplea(iclpro))
3952 c  !   &       *dsqrt(x)**dble(-alplea(icltar))
3953 c        write(ifhi,*) x,y
3954 c      enddo
3955 c
3956 c      write(ifhi,'(a)')    '  endarray'
3957 c      write(ifhi,'(a)')    'closehisto plot 0-'
3958 c
3959 cc**********************************************************************
3960 cc**********************************************************************
3961 c      do k=1,koll
3962 c        bk(k)=biniDf
3963 c      enddo
3964 c      call GfunPark(0)
3965 c      call integom1(0)
3966 c
3967 cc********************* points = PhiExpoK*********************************
3968 c
3969 c      if (biniDf.lt.1.) then
3970 c        k=int(10.*biniDf)
3971 c        write(ifhi,'(a,I1)')  'openhisto name PhiExpok-b0.',k
3972 c      else
3973 c        write(ifhi,'(a,f3.1)')  'openhisto name PhiExpok-b',biniDf
3974 c      endif
3975 c      write(ifhi,'(a)')       'htyp pfc'
3976 c      write(ifhi,'(a)')       'xmod log ymod lin'
3977 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3978 c      write(ifhi,'(a,2e11.3)')'yrange auto auto'
3979 c      write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3980 c      write(ifhi,'(a)') 'text 0 0.1 "yaxis  [F](x+,x-)/x^[a]?remn!!"'
3981 c      if (xpar8.eq.1.) then
3982 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3983 c        write(ifhi,'(a,f5.2,a)')'text 0.5 0.9 "b=',biniDf,' fm"'
3984 c      endif
3985 c      write(ifhi,'(a)')       'array 2'
3986 c
3987 c      do i=0,nptg
3988 c        x=xminr
3989 c        if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3990 c        y=PhiExpoK(1,dsqrt(x),dsqrt(x))
3991 c        write(ifhi,*) x,y
3992 c      enddo
3993 c
3994 c      write(ifhi,'(a)')    '  endarray'
3995 c      write(ifhi,'(a)')    'closehisto plot 0'
3996 c
3997 cc************************* red = PhiMExact*omGam ************************
3998 c
3999 c      write(ifhi,'(a)')       'openhisto name GPhiExact'
4000 c      write(ifhi,'(a)')       'htyp lru'
4001 c      if(xpar5.eq.0.)then
4002 c        write(ifhi,'(a)')       'xmod lin ymod lin'
4003 c      else
4004 c        write(ifhi,'(a)')       'xmod log ymod log'
4005 c      endif
4006 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4007 c      write(ifhi,'(a)')       'yrange auto auto'
4008 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
4009 c      write(ifhi,'(a,a)') 'text 0 0.1 "yaxis  G(x+,x-)*[F]'
4010 c     *,'(x+,x-)/x^[a]?remn!!"'
4011 c      if (xpar8.eq.1.) then
4012 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.2 "s=',engy**2,' GeV^2!"'
4013 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "b=',biniDf,' fm"'
4014 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "x-=',xm,'"'
4015 c      endif
4016 c      write(ifhi,'(a)')       'array 2'
4017 c
4018 c      do i=0,nptg
4019 c        x=xminr
4020 c        if(xpar5.ne.0.)then
4021 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4022 c        else
4023 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4024 c        endif
4025 cc        z=1.d0-dsqrt(x)
4026 c        v(i)=0.d0
4027 c        if(engy**2..lt.5.e06)
4028 c     *  v(i)=Phiexact(0.,0.,1.,1.d0-x,1.d0-xm,engy**2,biniDf)
4029 cc     *  v(i)=Phiexact(0.,0.,1.,z,z,engy**2,biniDf)
4030 c        u(i)=omGam(x,xm,biniDf)
4031 cc        u(i)=omGam(dsqrt(x),dsqrt(x),biniDf)
4032 c        y=u(i)*v(i)
4033 c        if(xpar5.ne.0.)y=dabs(y)
4034 c        write(ifhi,*)x,y
4035 c      enddo
4036 c
4037 c      write(ifhi,'(a)')    '  endarray'
4038 c      write(ifhi,'(a)')    'closehisto plot 0-'
4039 c
4040 cc************************* red = PhiMExpo*omGam ************************
4041 c
4042 c      write(ifhi,'(a)')       'openhisto name GPhiExpo'
4043 c      write(ifhi,'(a)')       'htyp lba'
4044 c      if(xpar5.eq.0.)then
4045 c        write(ifhi,'(a)')       'xmod lin ymod lin'
4046 c      else
4047 c        write(ifhi,'(a)')       'xmod log ymod log'
4048 c      endif
4049 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4050 c      write(ifhi,'(a)')       'yrange auto auto'
4051 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
4052 c      write(ifhi,'(a,a)')
4053 c     * 'text 0 0.1 "yaxis  G(x+,x-)*[F]?'
4054 c     * ,'(1-x+,1-x-)/x^[a]?remn!!"'
4055 c      if (xpar8.eq.1.) then
4056 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4057 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.8 "b=',biniDf,' fm"'
4058 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.7 "x-=',xm,'"'
4059 c      endif
4060 c      write(ifhi,'(a)')       'array 2'
4061 c
4062 c      do i=0,nptg
4063 c        x=xminr
4064 c        if(xpar5.ne.0.)then
4065 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4066 c        else
4067 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4068 c        endif
4069 cc        z=1.d0-dsqrt(x)
4070 cc        w(i)=Phiexpo(0.,0.,1.,z,z,engy**2,biniDf)
4071 c        w(i)=Phiexpo(0.,0.,1.,1.d0-x,1.d0-xm,engy**2,biniDf)
4072 c        y=u(i)*w(i)
4073 c        if(xpar5.ne.0.)y=dabs(y)
4074 c        write(ifhi,*)x,y
4075 c      enddo
4076 c
4077 c      write(ifhi,'(a)')    '  endarray'
4078 c      if(xpar5.ne.0.)then
4079 c      write(ifhi,'(a)')    'closehisto plot 0-'
4080 c
4081 cc************************* green = omGam ************************
4082 c
4083 c      write(ifhi,'(a)')       'openhisto name GM'
4084 c      write(ifhi,'(a)')       'htyp lgo'
4085 c      write(ifhi,'(a)')       'array 2'
4086 c
4087 c      do i=0,nptg
4088 c        x=xminr
4089 c        if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4090 c        write(ifhi,*)x,dabs(u(i))
4091 c      enddo
4092 c
4093 c      write(ifhi,'(a)')    '  endarray'
4094 c      write(ifhi,'(a)')    'closehisto plot 0-'
4095 c
4096 cc************************* circle = PhiMExact  ************************
4097 c
4098 c      write(ifhi,'(a)')       'openhisto name PhiExact'
4099 c      write(ifhi,'(a)')       'htyp poc'
4100 c      write(ifhi,'(a)')       'array 2'
4101 c
4102 c      do i=0,nptg
4103 c        x=xminr
4104 c        if(xpar5.ne.0.)then
4105 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4106 c        else
4107 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4108 c        endif
4109 c        y=v(i)
4110 c        if(xpar5.ne.0.)y=dabs(y)
4111 c        write(ifhi,*)x,y
4112 c      enddo
4113 c
4114 c      write(ifhi,'(a)')    '  endarray'
4115 c      write(ifhi,'(a)')    'closehisto plot 0-'
4116 c
4117 cc************************* triangle = PhiMExpo ************************
4118 c
4119 c      write(ifhi,'(a)')       'openhisto name PhiExpo'
4120 c      write(ifhi,'(a)')       'htyp pot'
4121 c      write(ifhi,'(a)')       'array 2'
4122 c
4123 c      do i=0,nptg
4124 c        x=xminr
4125 c        if(xpar5.ne.0.)then
4126 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4127 c        else
4128 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4129 c        endif
4130 c        y=w(i)
4131 c        if(xpar5.ne.0.)y=dabs(y)
4132 c        write(ifhi,*)x,y
4133 c      enddo
4134 c
4135 c      write(ifhi,'(a)')    '  endarray'
4136 c      endif
4137 c      write(ifhi,'(a)')    'closehisto plot 0'
4138 c
4139 c      end
4140 c
4141
4142 c----------------------------------------------------------------------
4143       subroutine xParPhi
4144 c----------------------------------------------------------------------
4145
4146       include 'epos.inc'
4147       include 'epos.incems'
4148       include 'epos.incsem'
4149       include 'epos.incpar'
4150       double precision x,xminr,y,z(0:200)!,Zn,Znorm
4151       double precision PhiExpo,PhiExact,PhiExpoK
4152
4153
4154       nptg=10                  !number of point for the graphs
4155       biniDf=xpar2              !value of biniDf (impact parameter)
4156       xminr=max(1.d-6,1.d0/dble(engy**2))  !value of xminr for plotting the function
4157       zz=xpar7
4158
4159 c********************** full-red = PhiExact ***************************
4160
4161       if (biniDf.lt.1.) then
4162         k=int(10.*biniDf)
4163         write(ifhi,'(a,I1)')  'openhisto name PhiExact-b0.',k
4164       else
4165         write(ifhi,'(a,f3.1)')  'openhisto name PhiExact-b',biniDf
4166       endif
4167       write(ifhi,'(a)')       'htyp lru'
4168       write(ifhi,'(a)')       'xmod lin ymod lin'
4169       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4170       write(ifhi,'(a)')       'yrange auto auto'
4171       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
4172       write(ifhi,'(a)') 'text 0 0 "yaxis  [F](x)/x^[a]"'
4173       write(ifhi,'(a,i4,a,f4.1,a)')
4174      * 'txt  "title E=',nint(engy),' b=',biniDf,'"'
4175       write(ifhi,'(a)')       'array 2'
4176
4177       do i=0,nptg
4178         !x=xminr
4179         !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4180         x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4181         y=Phiexact(zz,zz,1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
4182      &       *dsqrt(x)**dble(-alplea(iclpro))
4183      &       *dsqrt(x)**dble(-alplea(icltar))
4184         write(ifhi,*)x,y
4185       enddo
4186
4187       write(ifhi,'(a)')    '  endarray'
4188       write(ifhi,'(a)')    'closehisto plot 0-'
4189
4190 c******************** blue = PhiExpo ***************************
4191
4192       if (biniDf.lt.1.) then
4193         k=int(10.*biniDf)
4194         write(ifhi,'(a,I1)')  'openhisto name PhiExpo-b0.',k
4195       else
4196         write(ifhi,'(a,f3.1)')  'openhisto name PhiExpo-b',biniDf
4197       endif
4198       write(ifhi,'(a)')       'htyp lbu'
4199       write(ifhi,'(a)')       'xmod lin ymod lin'
4200       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4201       write(ifhi,'(a,2e11.3)')'yrange auto auto'
4202       write(ifhi,'(a)')       'array 2'
4203
4204       do i=0,nptg
4205         !x=xminr
4206         !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4207         x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4208         z(i)=Phiexpo(zz,zz,1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
4209      &       *dsqrt(x)**dble(-alplea(iclpro))
4210      &       *dsqrt(x)**dble(-alplea(icltar))
4211         write(ifhi,*) x,z(i)
4212       enddo
4213
4214       write(ifhi,'(a)')    '  endarray'
4215 c      write(ifhi,'(a)')    'closehisto plot 0-'
4216 c
4217 cc*********************yellow = PhiUnit*********************************
4218 c
4219 c      if (biniDf.lt.1.) then
4220 c        k=int(10.*biniDf)
4221 c        write(ifhi,'(a,I1)')  'openhisto name PhiUnit-b0.',k
4222 c      else
4223 c        write(ifhi,'(a,f3.1)')  'openhisto name PhiUnit-b',biniDf
4224 c      endif
4225 c      write(ifhi,'(a)')       'htyp lyu'
4226 c      write(ifhi,'(a)')       'xmod lin ymod lin'
4227 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4228 c      write(ifhi,'(a,2e11.3)')'yrange auto auto'
4229 c      write(ifhi,'(a)')       'array 2'
4230 c
4231 c      Zn=Znorm(engy**2,biniDf)
4232 c
4233 c      do i=0,nptg
4234 c        !x=xminr
4235 c        !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4236 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4237 c        write(ifhi,*) x,z(i)/Zn
4238 c      enddo
4239 c
4240 c      write(ifhi,'(a)')    '  endarray'
4241
4242       if(koll.ge.1)then
4243
4244       write(ifhi,'(a)')    'closehisto plot 0-'
4245
4246 c**********************************************************************
4247 c**********************************************************************
4248       do k=1,koll
4249         bk(k)=biniDf
4250       enddo
4251       call GfunPark(0)
4252       call integom1(0)
4253
4254
4255 c*********************green = PhiExpoK*********************************
4256
4257       if (biniDf.lt.1.) then
4258         k=int(10.*biniDf)
4259         write(ifhi,'(a,I1)')  'openhisto name PhiExpok-b0.',k
4260       else
4261         write(ifhi,'(a,f3.1)')  'openhisto name PhiExpok-b',biniDf
4262       endif
4263       write(ifhi,'(a)')       'htyp lga'
4264       write(ifhi,'(a)')       'xmod lin ymod lin'
4265       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4266       write(ifhi,'(a,2e11.3)')'yrange auto auto'
4267       write(ifhi,'(a)')       'array 2'
4268
4269       do i=0,nptg
4270         !x=xminr
4271         !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4272         x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4273         z(i)=PhiExpoK(1,dsqrt(x),dsqrt(x))
4274         write(ifhi,*) x,z(i)
4275       enddo
4276
4277       write(ifhi,'(a)')    '  endarray'
4278       endif
4279       write(ifhi,'(a)')    'closehisto plot 0'
4280
4281 c
4282       end
4283
4284 c----------------------------------------------------------------------
4285       subroutine xParH
4286 c----------------------------------------------------------------------
4287
4288       include 'epos.inc'
4289       include 'epos.incpar'
4290       parameter(idxD2=8)
4291       double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
4292       common/DGamUni/GbetUni(  idxD0:idxD2),HbetUni(  idxD0:idxD2),
4293      &               GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
4294      &               HalpUni(idxD0:idxD2)
4295       double precision x,xminr,y,xm,utgam2
4296       double precision Hrst
4297
4298       nptg=20                  !number of point for the graphs
4299       biniDf=xpar2              !value of biniDf (impact parameter)
4300       xm=dble(xpar4)            !value of xminus
4301 c.....xminr=0.d0   !value of xminr for plotting the function
4302       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
4303
4304       imax=idxD1
4305       if(iomega.eq.2)imax=1
4306       do i=idxDmin,imax
4307         call Gfunpar(0.,0.,1,i,biniDf,smaxDf,alpx,betx,betpx,epsp,epst
4308      &               ,epss,gamv)
4309       enddo
4310       call GfomPar(biniDf,smaxDf)
4311       imax0=idxD1
4312       if(iomega.eq.2)imax0=1
4313       imax1=idxD2
4314       if(iomega.eq.2)imax1=imax1-1
4315       do i=idxDmin,imax0
4316         GbetUni(i)=utgam2(betUni(i,1)+1.d0)
4317         GbetpUni(i)=utgam2(betpUni(i,1)+1.d0)
4318         HbetUni(i)=utgam2(GbetUni(i))
4319         HbetpUni(i)=utgam2(GbetpUni(i))
4320         HalpUni(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
4321       enddo
4322       do i=0,1
4323         HbetUni(imax0+1+i)=betUni(i,1)+1.d0+betfom
4324         HbetUni(imax0+3+i)=betUni(i,1)+1.d0
4325         HbetUni(imax0+5+i)=betUni(i,1)+1.d0+betfom
4326         HbetpUni(imax0+1+i)=betpUni(i,1)+1.d0
4327         HbetpUni(imax0+3+i)=betpUni(i,1)+1.d0+betfom
4328         HbetpUni(imax0+5+i)=betpUni(i,1)+1.d0+betfom
4329         GbetUni(imax0+1+i)=utgam2(HbetUni(imax0+1+i))
4330         GbetUni(imax0+3+i)=utgam2(HbetUni(imax0+3+i))
4331         GbetUni(imax0+5+i)=utgam2(HbetUni(imax0+5+i))
4332         GbetpUni(imax0+1+i)=utgam2(HbetpUni(imax0+1+i))
4333         GbetpUni(imax0+3+i)=utgam2(HbetpUni(imax0+3+i))
4334         GbetpUni(imax0+5+i)=utgam2(HbetpUni(imax0+5+i))
4335         HalpUni(imax0+1+i)=zztUni*alpUni(i,1)
4336         HalpUni(imax0+3+i)=zzpUni*alpUni(i,1)
4337         HalpUni(imax0+5+i)=zzpUni*zztUni*alpUni(i,1)
4338       enddo
4339
4340 c***********************  red = Hrst  *********************************
4341
4342       write(ifhi,'(a)')       'openhisto name Hrst'
4343       write(ifhi,'(a)')       'htyp lru'
4344       write(ifhi,'(a)')       'xmod log ymod lin'
4345       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4346       write(ifhi,'(a)')       'yrange auto auto'
4347       write(ifhi,'(a)')
4348      *     'text 0 0 "yaxis  H?2!(x+,x-)"'
4349       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
4350       write(ifhi,'(a)')
4351      *     'text 0 0 "yaxis  H?2!(x+,x-)"'
4352       if (xpar8.eq.1.) then
4353         write(ifhi,'(a,e7.2,a)')'text 0.1 0.2 "s=',engy**2,' GeV^2!"'
4354         write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "b=',biniDf,' fm"'
4355         write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "x-=',xm,'"'
4356       endif
4357       write(ifhi,'(a)')       'array 2'
4358
4359       do i=0,nptg
4360         x=xminr
4361         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4362 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4363         y=Hrst(smaxDf,biniDf,dsqrt(x),dsqrt(x))
4364         write(ifhi,*)x,y
4365       enddo
4366
4367       write(ifhi,'(a)')    '  endarray'
4368       write(ifhi,'(a)')    'closehisto plot 0'
4369
4370       end
4371
4372
4373
4374
4375 cc----------------------------------------------------------------------
4376 c      subroutine xParHPhiIntnew
4377 cc----------------------------------------------------------------------
4378 c
4379 c      include 'epos.inc'
4380 c      include 'epos.incsem'
4381 c      include 'epos.incems'
4382 c      include 'epos.incpar'
4383 c      double precision x,xminr,xm,y
4384 c      double precision PhiExact,omGam
4385 cc      double precision PhiExpo
4386 c
4387 c      nptg=30                  !number of point for the graphs
4388 c      biniDf=xpar2              !value of biniDf (impact parameter)
4389 c      xm=dble(xpar4)            !value of xminus
4390 cc.....xminr=0.d0   !value of xminr for plotting the function
4391 c      xminr=1.d-3 !/dble(engy**2)  !value of xminr for plotting the function
4392 cc************************* black = PhiExact ***************************
4393 c
4394 c      write(ifhi,'(a)')       'openhisto name Phi1Exact'
4395 c      write(ifhi,'(a)')       'htyp lru'
4396 c      if(xpar5.eq.0.)then
4397 c        write(ifhi,'(a)')       'xmod lin ymod lin'
4398 c      else
4399 c        write(ifhi,'(a)')       'xmod log ymod lin'
4400 c      endif
4401 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4402 c      write(ifhi,'(a)')       'yrange auto auto'
4403 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
4404 c      write(ifhi,'(a)')
4405 c     * 'text 0 0.1 "yaxis  [F]?(x+,x-)/x^[a]?remn!!"'
4406 c      if (xpar8.eq.1.) then
4407 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.2 "s=',engy**2,' GeV^2!"'
4408 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "b=',biniDf,' fm"'
4409 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "x-=',xm,'"'
4410 c      endif
4411 c      write(ifhi,'(a)')       'array 2'
4412 c
4413 c      do i=0,nptg
4414 c        x=xminr
4415 c        if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4416 cc.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4417 c      y=Phiexact(0.,0.,1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
4418 c     &       *omGam(dsqrt(x),dsqrt(x),biniDf)
4419 c        write(ifhi,*)x,y
4420 c      enddo
4421 c
4422 c      write(ifhi,'(a)')    '  endarray'
4423 c      write(ifhi,'(a)')    'closehisto plot 0'
4424 c
4425 c      end
4426 c
4427 c----------------------------------------------------------------------
4428       subroutine xParHPhiInt
4429 c----------------------------------------------------------------------
4430
4431       include 'epos.inc'
4432       include 'epos.incpar'
4433       double precision y,HPhiInt
4434       common/geom/rmproj,rmtarg,bmax,bkmx
4435       parameter(nbib=32)
4436
4437
4438 c************************ dotted = gauss integration ******************
4439
4440       b1=0
4441       b2=max(abs(bkmx),3.)*1.2
4442       db=(b2-b1)/nbib
4443
4444       write(ifhi,'(a)')       'openhisto name HPhiExpoInt'
4445       write(ifhi,'(a)')       'htyp pfc xmod lin ymod lin'
4446       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4447       write(ifhi,'(a)')    'text 0 0 "yaxis  Int(H[F]?pp!)(s,b)"'
4448       if (xpar8.eq.1.) then
4449         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4450       endif
4451       write(ifhi,'(a)')       'array 2'
4452       imax=idxD1
4453       if(iomega.eq.2)imax=1
4454       do k=1,nbib
4455         b=b1+(k-0.5)*db
4456         do i=idxDmin,imax
4457           call Gfunpar(0.,0.,1,i,b,smaxDf,alpx,betx,betpx,epsp,epst,epss
4458      &                ,gamv)
4459           call Gfunpar(0.,0.,2,i,b,smaxDf,alpx,betx,betpx,epsp,epst,epss
4460      &                ,gamv)
4461         enddo
4462         y=HPhiInt(smaxDf,b)
4463         write(ifhi,*)b,y
4464       enddo
4465       write(ifhi,'(a)')    '  endarray'
4466       write(ifhi,'(a)')    'closehisto plot 0'
4467
4468       end
4469
4470 c----------------------------------------------------------------------
4471       subroutine xParZ
4472 c----------------------------------------------------------------------
4473       include 'epos.inc'
4474       include 'epos.incpar'
4475       double precision Znorm,y
4476       common/geom/rmproj,rmtarg,bmax,bkmx
4477       parameter(nbib=12)
4478
4479       b1=0
4480       b2=max(abs(bkmx),3.)*1.2
4481       db=(b2-b1)/nbib
4482
4483 c************************full-red = Znorm *****************************
4484
4485       write(ifhi,'(a)')       'openhisto name Znorm'
4486       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4487       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4488       write(ifhi,'(a)')    'text 0 0 "yaxis  Z(s,b) "'
4489       if (xpar8.eq.1.) then
4490         write(ifhi,'(a,e7.2,a)')'text 0.1 0.1 "s=',engy**2,' GeV^2!"'
4491       endif
4492       write(ifhi,'(a)')       'array 2'
4493 c        y=Znorm(engy**2,xpar2)
4494       do k=1,nbib
4495         b=b1+(k-0.5)*db
4496         y=Znorm(smaxDf,b)
4497         write(ifhi,*)b,y
4498       enddo
4499       write(ifhi,'(a)')    '  endarray'
4500       write(ifhi,'(a)')    'closehisto plot 0-'
4501
4502
4503 c**********************************************************************
4504
4505       write(ifhi,'(a)')       'openhisto name un'
4506       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4507       write(ifhi,'(a)')       'array 2'
4508       do k=1,nbib
4509         b=b1+(k-0.5)*db
4510         y=1
4511         write(ifhi,*)b,y
4512       enddo
4513       write(ifhi,'(a)')    '  endarray'
4514       write(ifhi,'(a)')    'closehisto plot 0'
4515
4516
4517       end
4518
4519 c----------------------------------------------------------------------
4520       subroutine xParPro
4521 c----------------------------------------------------------------------
4522       include 'epos.inc'
4523       include 'epos.incsem'
4524       include 'epos.incpar'
4525       double precision PhiExact,PhiExpo,y,om1intb,om1intbc,om1intgc
4526      &,om1intbi!,PhiUnit
4527       common/geom/rmproj,rmtarg,bmax,bkmx
4528       parameter(nbib=12)
4529
4530       b1=0
4531       b2=max(abs(bkmx),3.)*1.2
4532       db=(b2-b1)/nbib
4533       zz=xpar7
4534
4535 c********************* full-red = 1-PhiExact **************************
4536
4537       write(ifhi,'(a)')       'openhisto name 1-PhiExact'
4538       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4539       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4540       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4541       if (xpar8.eq.1.) then
4542         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4543       endif
4544       write(ifhi,'(a)')       'array 2'
4545       do k=1,nbib
4546         b=b1+(k-0.5)*db
4547         y=1.d0-Phiexact(zz,zz,1.,1.d0,1.d0,engy**2,b)
4548         write(ifhi,*)b,y
4549       enddo
4550       write(ifhi,'(a)')    '  endarray'
4551       write(ifhi,'(a)')    'closehisto plot 0-'
4552
4553 c************************** blue-dashed = 1-PhiExpo *******************
4554
4555       write(ifhi,'(a)')       'openhisto name 1-PhiExpo'
4556       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4557       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4558       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4559       if (xpar8.eq.1.) then
4560         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4561       endif
4562       write(ifhi,'(a)')       'array 2'
4563       do k=1,nbib
4564         b=b1+(k-0.5)*db
4565         y=1.d0-Phiexpo(zz,zz,1.,1.d0,1.d0,engy**2,b)
4566         write(ifhi,*)b,y
4567       enddo
4568       write(ifhi,'(a)')    '  endarray'
4569       write(ifhi,'(a)')    'closehisto plot 0-'
4570
4571
4572 c**********************************************************************
4573
4574       write(ifhi,'(a)')       'openhisto'
4575       write(ifhi,'(a)')       'htyp lga xmod lin ymod lin'
4576       write(ifhi,'(a)')       'array 2'
4577       do k=1,nbib
4578         b=b1+(k-0.5)*db
4579         y=1
4580         write(ifhi,*)b,y
4581       enddo
4582       write(ifhi,'(a)')    '  endarray'
4583       write(ifhi,'(a)')    'closehisto plot 0'
4584
4585 c****************************** red = om1intbc ********************
4586
4587       write(ifhi,'(a)')       'openhisto name om1intbc'
4588       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4589       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4590       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1bc!(b) "'
4591       if (xpar8.eq.1.) then
4592         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4593       endif
4594       write(ifhi,'(a)')       'array 2'
4595
4596       do k=1,nbib
4597         b=b1+(k-0.5)*db
4598         y=om1intbc(b)
4599         write(ifhi,*)b,y
4600       enddo
4601       write(ifhi,'(a)')    '  endarray'
4602       write(ifhi,'(a)')    'closehisto plot 0-'
4603
4604
4605 c************************* blue dashed =  om1intb ********************
4606
4607       write(ifhi,'(a)')       'openhisto name om1intb'
4608       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4609       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4610       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1b!(b) "'
4611       if (xpar8.eq.1.) then
4612         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4613       endif
4614       write(ifhi,'(a)')       'array 2'
4615
4616       do k=1,nbib
4617         b=b1+(k-0.5)*db
4618         y=om1intb(b)
4619         write(ifhi,*)b,y
4620       enddo
4621       write(ifhi,'(a)')    '  endarray'
4622       write(ifhi,'(a)')    'closehisto plot 0-'
4623
4624 c****************************** green dot = om1intgc ********************
4625
4626       write(ifhi,'(a)')       'openhisto name om1intgc'
4627       write(ifhi,'(a)')       'htyp lgo xmod lin ymod lin'
4628       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4629       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1gc!(b) "'
4630       if (xpar8.eq.1.) then
4631         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4632       endif
4633       write(ifhi,'(a)')       'array 2'
4634
4635       do k=1,nbib
4636         b=b1+(k-0.5)*db
4637         y=om1intgc(b)
4638         write(ifhi,*)b,y
4639       enddo
4640       write(ifhi,'(a)')    '  endarray'
4641       write(ifhi,'(a)')    'closehisto plot 0'
4642
4643 c****************************** red = om1intbi(0) ********************
4644
4645       write(ifhi,'(a)')       'openhisto name om1intbc'
4646       write(ifhi,'(a)')       'htyp lru xmod lin ymod log'
4647       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4648       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1bc!(b) "'
4649       if (xpar8.eq.1.) then
4650         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4651       endif
4652       write(ifhi,'(a)')       'array 2'
4653
4654       do k=1,nbib
4655         b=b1+(k-0.5)*db
4656         y=om1intbi(b,0)
4657         write(ifhi,*)b,y
4658       enddo
4659       write(ifhi,'(a)')    '  endarray'
4660       write(ifhi,'(a)')    'closehisto plot 0-'
4661
4662
4663 c************************* blue dashed =  om1intbi(1) ********************
4664
4665       write(ifhi,'(a)')       'openhisto name om1intb'
4666       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4667       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4668       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1b!(b) "'
4669       if (xpar8.eq.1.) then
4670         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4671       endif
4672       write(ifhi,'(a)')       'array 2'
4673
4674       do k=1,nbib
4675         b=b1+(k-0.5)*db
4676         y=om1intbi(b,1)
4677         write(ifhi,*)b,y
4678       enddo
4679       write(ifhi,'(a)')    '  endarray'
4680       write(ifhi,'(a)')    'closehisto plot 0-'
4681
4682 c****************************** green dot = om1intbi(2) ********************
4683
4684       write(ifhi,'(a)')       'openhisto name om1intgc'
4685       write(ifhi,'(a)')       'htyp lgo xmod lin ymod lin'
4686       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4687       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1gc!(b) "'
4688       if (xpar8.eq.1.) then
4689         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4690       endif
4691       write(ifhi,'(a)')       'array 2'
4692
4693       do k=1,nbib
4694         b=b1+(k-0.5)*db
4695         y=om1intbi(b,2)
4696         write(ifhi,*)b,y
4697       enddo
4698       write(ifhi,'(a)')    '  endarray'
4699       write(ifhi,'(a)')    'closehisto plot 0'
4700
4701
4702
4703       end
4704
4705 c----------------------------------------------------------------------
4706       subroutine xParPro1
4707 c----------------------------------------------------------------------
4708       include 'epos.inc'
4709       include 'epos.incsem'
4710       include 'epos.incpar'
4711       double precision PhiExact,PhiExpo,y
4712       common/geom/rmproj,rmtarg,bmax,bkmx
4713       parameter(nbib=12)
4714
4715       b1=0
4716       b2=max(abs(bkmx),3.)*1.2
4717       db=(b2-b1)/nbib
4718       zz=xpar7
4719
4720 c********************* full-red = 1-PhiExact **************************
4721
4722       write(ifhi,'(a)')       'openhisto name 1-PhiExact'
4723       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4724       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4725       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4726       if (xpar8.eq.1.) then
4727         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4728       endif
4729       write(ifhi,'(a)')       'array 2'
4730       do k=1,nbib
4731         b=b1+(k-0.5)*db
4732         y=1.d0-Phiexact(zz,zz,.5,1.d0,1.d0,engy**2,b)
4733         write(ifhi,*)b,y
4734       enddo
4735       write(ifhi,'(a)')    '  endarray'
4736       write(ifhi,'(a)')    'closehisto plot 0-'
4737
4738 c************************** blue-dashed = 1-PhiExpo *******************
4739
4740       write(ifhi,'(a)')       'openhisto name 1-PhiExpo'
4741       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4742       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4743       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4744       if (xpar8.eq.1.) then
4745         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4746       endif
4747       write(ifhi,'(a)')       'array 2'
4748       do k=1,nbib
4749         b=b1+(k-0.5)*db
4750         y=1.d0-Phiexpo(zz,zz,.5,1.d0,1.d0,engy**2,b)
4751         write(ifhi,*)b,y
4752       enddo
4753       write(ifhi,'(a)')    '  endarray'
4754       write(ifhi,'(a)')    'closehisto plot 0'
4755
4756
4757
4758
4759       end
4760
4761 c----------------------------------------------------------------------
4762       subroutine xParGam
4763 c----------------------------------------------------------------------
4764
4765       include 'epos.inc'
4766       include 'epos.incsem'
4767       include 'epos.incpar'
4768       dimension bet(idxD0:idxD1)
4769       double precision utgam2,xgammag2!,xrem
4770       dimension ip(idxD0:idxD1),imax(idxD0:idxD1)
4771
4772       nptg=50                  !number of point for the graphs
4773       b=xpar2
4774       gamp=xpar6
4775       zmax=6.
4776 c      xrem=dble(xpar4)
4777       if(idxD0.ne.0.or.idxD1.ne.2) stop "Check xPargam"
4778
4779       do i=idxD0,idxD1
4780         imax(i)=4
4781         bet(i)=0
4782       enddo
4783       nmax=idxD1
4784       imax(idxD0)=int(zmax)
4785       imax(1)=imax(idxD0)
4786       imax(2)=imax(idxD0)
4787
4788       do i=idxD0,nmax
4789         gam=gamD(i,iclpro,icltar)*b**2
4790         bet(i)=gam+betDp(i,iclpro,icltar)-alppar+1.
4791       enddo
4792       write(ifhi,'(a)')       'openhisto name gExact'
4793       write(ifhi,'(a)')       'htyp pfs'
4794       write(ifhi,'(a)')       'xmod lin ymod log'
4795       write(ifhi,'(a,2e11.3)')'xrange',0.,zmax
4796       write(ifhi,'(a)')       'yrange auto auto'
4797 c      write(ifhi,'(a)')'yrange 1.e-10 2'
4798       write(ifhi,'(a)')    'text 0 0 "xaxis z"'
4799       write(ifhi,'(a)')    'text 0 0 "yaxis g(z) "'
4800       write(ifhi,'(a)')       'array 2'
4801
4802       do ip0=0,imax(0)
4803         ip(0)=ip0
4804         do ip1=0,imax(1)
4805           ip(1)=ip1
4806           do ip2=0,imax(2)
4807             ip(2)=ip2
4808           t=0.
4809           do i=idxD0,nmax
4810             t=t+real(ip(i))*bet(i)
4811           enddo
4812           write(ifhi,'(2e14.6)')t,utgam2(dble(alplea(2))+1.D0)
4813      &       /utgam2(dble(alplea(2))+1.D0+dble(t))
4814         enddo
4815       enddo
4816       enddo
4817
4818       write(ifhi,'(a)')    '  endarray'
4819       write(ifhi,'(a)')    'closehisto plot 0-'
4820
4821
4822 c**********************************************************************
4823
4824       write(ifhi,'(a)')       'openhisto name gExpo'
4825       write(ifhi,'(a)')       'htyp lbu'
4826       write(ifhi,'(a)')       'xmod lin ymod log'
4827       write(ifhi,'(a,2e11.3)')'xrange',0.,zmax
4828       write(ifhi,'(a)')'yrange auto auto'
4829       write(ifhi,'(a)')    'text 0 0 "xaxis z"'
4830       write(ifhi,'(a)')    'text 0 0 "yaxis g(z)"'
4831       write(ifhi,'(a)')       'array 2'
4832
4833       do i=0,nptg
4834         t=zmax*(real(i)/real(nptg))
4835         write(ifhi,'(2e14.6)') t,dexp(-dble(t))
4836       enddo
4837
4838       write(ifhi,'(a)')    '  endarray'
4839       write(ifhi,'(a)')    'closehisto plot 0-'
4840
4841
4842 c**********************************************************************
4843
4844
4845
4846       write(ifhi,'(a)')       'openhisto name gPower'
4847       write(ifhi,'(a)')       'htyp poc'
4848       write(ifhi,'(a)')       'xmod log ymod log'
4849       write(ifhi,'(a,2e11.3)')'xrange',0.,zmax
4850       write(ifhi,'(a)')       'yrange auto auto'
4851       write(ifhi,'(a)')    'text 0 0 "xaxis t"'
4852       write(ifhi,'(2a)')    'text 0 0 ',
4853      &     '"yaxis [P][G](1+[a]?L!)/[G](1+[a]?L!+[b])"'
4854       write(ifhi,'(a)')       'array 2'
4855
4856       do ip0=0,imax(0)
4857         ip(0)=ip0
4858         do ip1=0,imax(1)
4859           ip(1)=ip1
4860           do ip2=0,imax(2)
4861             ip(2)=ip2
4862 c$$$            do ip3=0,imax(3)
4863 c$$$              ip(3)=ip3
4864               t=0.
4865               do i=idxD0,nmax
4866                 t=t+real(ip(i))*bet(i)
4867               enddo
4868               write(ifhi,'(2e14.6)')t,xgammag2(iclpro,bet,ip,gamp)
4869 c$$$            enddo
4870           enddo
4871         enddo
4872       enddo
4873
4874       write(ifhi,'(a)')    '  endarray'
4875       write(ifhi,'(a)')    'closehisto plot 0'
4876
4877       end
4878
4879 c----------------------------------------------------------------------
4880       subroutine xParOmega1xy
4881 c----------------------------------------------------------------------
4882 c xpar2=b
4883 c xpar4=xh
4884 c xpar3=y
4885 c xpar7 : nucl coef
4886 c xpar8 : 2=xp/xm instead of xh/yp
4887 c----------------------------------------------------------------------
4888       include 'epos.incems'
4889       include 'epos.incsem'
4890       include 'epos.inc'
4891       include 'epos.incpar'
4892     !  double precision om1x,om1y
4893       double precision x,ranhis(0:51),y,ymax,xh
4894      &,om1xpk,om1xmk,t,om1xk,om1yk,xpr1,xmr1!,xp,xm,xmin,xmax
4895       common /psar7/ delx,alam3p,gam3p
4896
4897       nptg=50                  !number of point for the graphs
4898       biniDf=xpar2              !value of biniDf (impact paramter)
4899       xh=0.99d0
4900       xpr1=1.d0
4901       xmr1=1.d0
4902       if(xpar4.lt.1..and.xpar4.gt.0.)then
4903         xh=dble(xpar4**2)          !value of x
4904         xpr1=dble(xpar4)
4905         xmr1=dble(xpar4)
4906       endif
4907       do i=0,51
4908         ranhis(i)=0.d0
4909       enddo
4910 c$$$      xp=dsqrt(xh)*dble(exp(xpar3)) !y=xpar3
4911 c$$$      xm=1.d0
4912 c$$$      if(xp.ne.0.d0)xm=xh/xp
4913
4914
4915
4916       if(koll.eq.0)then
4917
4918       call xhistomran1(ranhis,biniDf)
4919
4920       stop'om1x not defined'
4921
4922
4923       write(ifhi,'(a)')       'openhisto name Om1x'
4924       write(ifhi,'(a)')       'htyp lin'
4925       write(ifhi,'(a)')       'xmod log ymod log'
4926       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4927       write(ifhi,'(a)')'yrange auto auto'
4928       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
4929       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x!"'
4930       if (xpar8.eq.1.) then
4931         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4932         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4933       endif
4934       write(ifhi,'(a)')       'array 2'
4935
4936       do i=0,nptg
4937         x=xminDf
4938         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
4939 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
4940         write(ifhi,*) x,   0   !om1x(x,biniDf)
4941       enddo
4942
4943       write(ifhi,'(a)')    '  endarray'
4944       write(ifhi,'(a)')    'closehisto plot 0-'
4945
4946       write(ifhi,'(a)')       'openhisto name Om1xRan'
4947       write(ifhi,'(a)')       'htyp his'
4948       write(ifhi,'(a)')       'xmod log ymod log'
4949       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4950       write(ifhi,'(a)')'yrange auto auto'
4951       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
4952       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x! random"'
4953       if (xpar8.eq.1.) then
4954         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4955         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4956       endif
4957       write(ifhi,'(a)')       'array 2'
4958
4959       do i=0,51
4960         x=xminDf
4961         if (i.ne.0) x=x*(xmaxDf/xminDf)**((dble(i)+.5d0)/51.d0)
4962 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
4963         write(ifhi,*) x,ranhis(i)
4964       enddo
4965
4966       write(ifhi,'(a)')    '  endarray'
4967       write(ifhi,'(a)')    'closehisto plot 0'
4968
4969       x=xh
4970       ymax=-.5D0*dlog(x)
4971       do i=0,51
4972         ranhis(i)=0.d0
4973       enddo
4974       call xhistomran2(ranhis,x,biniDf)
4975
4976       stop'om1y not defined'
4977
4978       write(ifhi,'(a)')       'openhisto name Om1y'
4979       write(ifhi,'(a)')       'htyp lin'
4980       write(ifhi,'(a)')       'xmod lin ymod lin'
4981       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
4982       write(ifhi,'(a)')'yrange auto auto'
4983       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
4984       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y!"'
4985       if (xpar8.eq.1.) then
4986         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4987         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4988       endif
4989       write(ifhi,'(a)')       'array 2'
4990
4991       do i=0,nptg
4992         y=-ymax+(2.d0*ymax)*(dble(i)/dble(nptg))
4993         write(ifhi,*) y, 0  !om1y(x,y,biniDf)
4994       enddo
4995
4996       write(ifhi,'(a)')    '  endarray'
4997       write(ifhi,'(a)')    'closehisto plot 0-'
4998
4999       write(ifhi,'(a)')       'openhisto name Om1yRan'
5000       write(ifhi,'(a)')       'htyp his'
5001       write(ifhi,'(a)')       'xmod lin ymod lin'
5002       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
5003       write(ifhi,'(a)')'yrange auto auto'
5004       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
5005       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y! random"'
5006       if (xpar8.eq.1.) then
5007         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5008         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5009       endif
5010       write(ifhi,'(a)')       'array 2'
5011
5012       do i=0,50
5013         y=-ymax+(2.d0*ymax)*(dble(i)/50.d0)
5014         write(ifhi,*) y,ranhis(i)
5015       enddo
5016
5017       write(ifhi,'(a)')    '  endarray'
5018       write(ifhi,'(a)')    'closehisto plot 0'
5019
5020
5021 c**********************************************************************
5022
5023       else
5024
5025       do k=1,koll
5026         bk(k)=biniDf
5027       enddo
5028       call GfunPark(0)
5029       call integom1(0)
5030
5031       if(xpar8.eq.2.)then
5032
5033         call xhistomran8(ranhis,xpr1,xmr1)
5034
5035
5036       write(ifhi,'(a)')       'openhisto name Om1xp'
5037       write(ifhi,'(a)')       'htyp lin'
5038       write(ifhi,'(a)')       'xmod log ymod log'
5039       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5040       write(ifhi,'(a)')'yrange auto auto'
5041       write(ifhi,'(a)')    'text 0 0 "xaxis X+"'
5042       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x+!"'
5043       if (xpar8.eq.1.) then
5044         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5045         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5046       endif
5047       write(ifhi,'(a)')       'array 2'
5048
5049       do i=0,nptg
5050         x=xminDf
5051         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5052 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
5053         t=om1xpk(x,xpr1,1)
5054         write(ifhi,*) x,t
5055       enddo
5056
5057       write(ifhi,'(a)')    '  endarray'
5058       write(ifhi,'(a)')    'closehisto plot 0-'
5059
5060       write(ifhi,'(a)')       'openhisto name Om1xpRan'
5061       write(ifhi,'(a)')       'htyp his'
5062       write(ifhi,'(a)')       'xmod log ymod log'
5063       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5064       write(ifhi,'(a)')'yrange auto auto'
5065       write(ifhi,'(a)')    'text 0 0 "xaxis X+"'
5066       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x+! random"'
5067       if (xpar8.eq.1.) then
5068         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5069         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5070       endif
5071       write(ifhi,'(a)')       'array 2'
5072
5073       do i=0,51
5074         x=xminDf
5075         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5076 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
5077         write(ifhi,*) x,ranhis(i)
5078       enddo
5079
5080       write(ifhi,'(a)')    '  endarray'
5081       write(ifhi,'(a)')    'closehisto plot 0'
5082
5083       do i=0,51
5084         ranhis(i)=0.d0
5085       enddo
5086
5087       call xhistomran9(ranhis,xh,xpr1,xmr1)
5088
5089
5090       write(ifhi,'(a)')       'openhisto name Om1xm'
5091       write(ifhi,'(a)')       'htyp lin'
5092       write(ifhi,'(a)')       'xmod log ymod log'
5093       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5094       write(ifhi,'(a)')'yrange auto auto'
5095       write(ifhi,'(a)')    'text 0 0 "xaxis X-"'
5096       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x-!"'
5097       if (xpar8.eq.1.) then
5098         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5099         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5100         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x+=',xh,'"'
5101       endif
5102       write(ifhi,'(a)')       'array 2'
5103
5104       do i=0,nptg
5105         x=xminDf
5106         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5107 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
5108         t=om1xmk(xh,x,xpr1,xmr1,1)
5109         write(ifhi,*) x,t
5110       enddo
5111
5112       write(ifhi,'(a)')    '  endarray'
5113       write(ifhi,'(a)')    'closehisto plot 0-'
5114
5115       write(ifhi,'(a)')       'openhisto name Om1xmRan'
5116       write(ifhi,'(a)')       'htyp his'
5117       write(ifhi,'(a)')       'xmod log ymod log'
5118       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5119       write(ifhi,'(a)')'yrange auto auto'
5120       write(ifhi,'(a)')    'text 0 0 "xaxis X-"'
5121       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x-! random"'
5122       if (xpar8.eq.1.) then
5123         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5124         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5125         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x+=',xh,'"'
5126       endif
5127       write(ifhi,'(a)')       'array 2'
5128
5129       do i=0,51
5130         x=xminDf
5131         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5132 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
5133         write(ifhi,*) x,ranhis(i)
5134       enddo
5135
5136       write(ifhi,'(a)')    '  endarray'
5137       write(ifhi,'(a)')    'closehisto plot 0'
5138
5139 c**********************************************************************
5140       else
5141
5142       call xhistomran10(ranhis)
5143
5144
5145       write(ifhi,'(a)')       'openhisto name Om1xk'
5146       write(ifhi,'(a)')       'htyp lin'
5147       write(ifhi,'(a)')       'xmod log ymod log'
5148       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5149       write(ifhi,'(a)')'yrange auto auto'
5150       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
5151       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x!"'
5152       if (xpar8.eq.1.) then
5153         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5154         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5155       endif
5156       write(ifhi,'(a)')       'array 2'
5157
5158       do i=0,nptg
5159         x=xminDf
5160         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5161 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
5162         t=om1xk(x,1)
5163         write(ifhi,*) x,t
5164       enddo
5165
5166       write(ifhi,'(a)')    '  endarray'
5167       write(ifhi,'(a)')    'closehisto plot 0-'
5168
5169       write(ifhi,'(a)')       'openhisto name Om1xpRan'
5170       write(ifhi,'(a)')       'htyp his'
5171       write(ifhi,'(a)')       'xmod log ymod log'
5172       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5173       write(ifhi,'(a)')'yrange auto auto'
5174       write(ifhi,'(a)')    'text 0 0 "xaxis X+"'
5175       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x+! random"'
5176       if (xpar8.eq.1.) then
5177         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5178         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5179       endif
5180       write(ifhi,'(a)')       'array 2'
5181
5182       do i=0,51
5183         x=xminDf
5184         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5185 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
5186         write(ifhi,*) x,ranhis(i)
5187       enddo
5188
5189       write(ifhi,'(a)')    '  endarray'
5190       write(ifhi,'(a)')    'closehisto plot 0'
5191
5192       ymax=-.5D0*dlog(xh)
5193       do i=0,51
5194         ranhis(i)=0.d0
5195       enddo
5196
5197       call xhistomran11(ranhis,xh)
5198
5199
5200
5201       write(ifhi,'(a)')       'openhisto name Om1yk'
5202       write(ifhi,'(a)')       'htyp lin'
5203       write(ifhi,'(a)')       'xmod lin ymod lin'
5204       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
5205       write(ifhi,'(a)')'yrange auto auto'
5206       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
5207       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y!"'
5208       if (xpar8.eq.1.) then
5209         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5210         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5211       endif
5212       write(ifhi,'(a)')       'array 2'
5213
5214       do i=0,nptg
5215         y=-ymax+(2.d0*ymax)*(dble(i)/dble(nptg))
5216         t=om1yk(xh,y,1)
5217         write(ifhi,*) y,t
5218       enddo
5219
5220       write(ifhi,'(a)')    '  endarray'
5221       write(ifhi,'(a)')    'closehisto plot 0-'
5222
5223       write(ifhi,'(a)')       'openhisto name Om1yRan'
5224       write(ifhi,'(a)')       'htyp his'
5225       write(ifhi,'(a)')       'xmod lin ymod lin'
5226       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
5227       write(ifhi,'(a)')'yrange auto auto'
5228       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
5229       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y! random"'
5230       if (xpar8.eq.1.) then
5231         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5232         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5233       endif
5234       write(ifhi,'(a)')       'array 2'
5235
5236       do i=0,50
5237         y=-ymax+(2.d0*ymax)*(dble(i)/50.d0)
5238         write(ifhi,*) y,ranhis(i)
5239       enddo
5240
5241       write(ifhi,'(a)')    '  endarray'
5242       write(ifhi,'(a)')    'closehisto plot 0'
5243
5244
5245       endif
5246
5247       endif
5248
5249       return
5250       end
5251
5252 c----------------------------------------------------------------------
5253       subroutine xRanPt
5254 c----------------------------------------------------------------------
5255 c xpar2=xcut
5256 c----------------------------------------------------------------------
5257       include 'epos.inc'
5258       parameter (nptg1=501)      !number of point for the graphs
5259       double precision ranhis(0:nptg1)
5260       common /cranpt/conv
5261
5262       nptg=nptg1-1
5263       xcut=xpar2              !value of biniDf (impact paramter)
5264       xfact=xpar3
5265       xadd=xpar4
5266       xmax=10.
5267       conv=10./float(nptg)
5268       if(xcut.le.0.)xcut=float(nptg)
5269       if(xfact.le.0.)xfact=1.
5270       do i=0,nptg1
5271         ranhis(i)=0.d0
5272       enddo
5273 c$$$      xp=dsqrt(xh)*dble(exp(xpar3)) !y=xpar3
5274 c$$$      xm=1.d0
5275 c$$$      if(xp.ne.0.d0)xm=xh/xp
5276
5277       if(xpar1.ge.1.)then
5278
5279       call xranptg(ranhis,xcut,xfact,xadd)
5280
5281
5282
5283       write(ifhi,'(a)')       'openhisto name ranpt'
5284       write(ifhi,'(a)')       'htyp lin'
5285       write(ifhi,'(a)')       'xmod lin ymod log'
5286       write(ifhi,'(a,2e11.3)')'xrange',0.,min(xmax,xfact*xcut+xadd+1.)
5287       write(ifhi,'(a)')'yrange auto auto'
5288       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5289       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5290       write(ifhi,'(a)')       'array 2'
5291
5292
5293       do i=0,nptg
5294         x=float(i)*conv
5295         write(ifhi,*) x,ranhis(i)
5296       enddo
5297
5298       write(ifhi,'(a)')    '  endarray'
5299       write(ifhi,'(a)')    'closehisto plot 0-'
5300
5301       endif
5302
5303       if(xpar1.ge.2.)then
5304
5305       call xranpte(ranhis,xcut,xfact,xadd)
5306
5307
5308
5309       write(ifhi,'(a)')       'openhisto name ranpt'
5310       write(ifhi,'(a)')       'htyp lin'
5311       write(ifhi,'(a)')       'xmod lin ymod log'
5312       write(ifhi,'(a,2e11.3)')'xrange',0.,min(xmax,xfact*xcut+xadd+1.)
5313       write(ifhi,'(a)')'yrange auto auto'
5314       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5315       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5316       write(ifhi,'(a)')       'array 2'
5317
5318
5319       do i=0,nptg
5320         x=float(i)*conv
5321         write(ifhi,*) x,ranhis(i)
5322       enddo
5323
5324       write(ifhi,'(a)')    '  endarray'
5325       write(ifhi,'(a)')    'closehisto plot 0-'
5326
5327
5328       endif
5329
5330       if(xpar1.ge.3.)then
5331
5332       call xranpts(ranhis,xcut,xfact,xadd)
5333
5334
5335
5336       write(ifhi,'(a)')       'openhisto name ranpt'
5337       write(ifhi,'(a)')       'htyp lin'
5338       write(ifhi,'(a)')       'xmod lin ymod log'
5339       write(ifhi,'(a,2e11.3)')'xrange',0.,min(xmax,xfact*xcut+xadd+1.)
5340       write(ifhi,'(a)')'yrange auto auto'
5341       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5342       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5343       write(ifhi,'(a)')       'array 2'
5344
5345
5346       do i=0,nptg
5347         x=float(i)*conv
5348         write(ifhi,*) x,ranhis(i)
5349       enddo
5350
5351       write(ifhi,'(a)')    '  endarray'
5352       write(ifhi,'(a)')    'closehisto plot 0-'
5353
5354       endif
5355
5356       call xranptc(ranhis,xcut,xfact,xadd)
5357
5358
5359
5360       write(ifhi,'(a)')       'openhisto name ranpt'
5361       write(ifhi,'(a)')       'htyp lin'
5362       write(ifhi,'(a)')       'xmod lin ymod log'
5363       if(xpar1.ge.1)then
5364         write(ifhi,'(a,2e11.3)')'xrange',0.,min(xmax,xfact*xcut+xadd+1.)
5365       else
5366         write(ifhi,'(a,2e11.3)')'xrange',0.,1.
5367       endif
5368       write(ifhi,'(a)')'yrange auto auto'
5369       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5370       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5371       write(ifhi,'(a)')       'array 2'
5372
5373
5374       do i=0,nptg
5375         x=float(i)*conv
5376         write(ifhi,*) x,ranhis(i)
5377       enddo
5378
5379       write(ifhi,'(a)')    '  endarray'
5380
5381
5382       return
5383       end
5384
5385 c----------------------------------------------------------------------
5386
5387       double precision function xgammag2(iclrem,bet,ip,gamp)
5388
5389 c----------------------------------------------------------------------
5390       include 'epos.inc'
5391       include 'epos.incsem'
5392       double precision utgam2
5393       dimension bet(idxD0:idxD1),ip(idxD0:idxD1)
5394
5395       xgammag2=1.d0
5396
5397       imax=idxD1
5398
5399       do i=idxD0,imax
5400       if(ip(i).ne.0) xgammag2=xgammag2
5401      &   *(utgam2(dble(alplea(iclrem))+1.d0+dble(gamp))
5402      &   /(max(0.d0,dble(int(gamp+0.5))+1))
5403      &   /utgam2(dble(alplea(iclrem)+bet(i)+gamp)+1.D0))
5404      &                                          **dble(ip(i))
5405       enddo
5406
5407       return
5408       end
5409
5410 c----------------------------------------------------------------------
5411
5412       function xsigmafit(x)
5413
5414 c----------------------------------------------------------------------
5415
5416       include 'epos.inc'
5417       include 'epos.incpar'
5418       double precision x,xDfit,sfsh,varifit,range,sig2
5419       double precision bf(maxdataDf),Db(maxdataDf)
5420       external varifit
5421
5422
5423       sig2=bmaxDf/2.
5424       range=sig2
5425       xp=real(dsqrt(x))
5426       xm=xp
5427       zz=xpar7
5428
5429
5430       sfsh=xDfit(zz,0,1,smaxDf,xp,xm,0.)
5431       if(dabs(sfsh).ge.1.d-5)then
5432       do i=0,nptf-1
5433         bf(i+1)=dble(-bmaxDf+real(i)*2.*bmaxDf/real(nptf-1))
5434         Db(i+1)=xDfit(zz,0,1,smaxDf,xp,xm,real(bf(i+1)))/sfsh
5435       enddo
5436
5437 c.....Fit of D(X,b) between -bmaxDf and bmaxDf
5438       call minfit(varifit,bf,Db,nptf,sig2,range)
5439
5440       xsigmafit=real(sig2)
5441       else
5442       xsigmafit=0.
5443       endif
5444
5445
5446       return
5447       end
5448
5449
5450 c----------------------------------------------------------------------
5451
5452       subroutine xhistomran1(histo,b)
5453
5454 c----------------------------------------------------------------------
5455 c.....Make Histogram of om1xr
5456 c----------------------------------------------------------------------
5457
5458       include 'epos.incpar'
5459
5460       double precision histo(0:51),x,x1  !,om1xr
5461       integer*4 n
5462
5463
5464       n=100000
5465       do i=0,51
5466         histo(i)=0.d0
5467       enddo
5468       do 111 j=1,n
5469         if(mod(j,10000).eq.0)write(*,*)"x1",j,b
5470         x=0  !om1xr(b)
5471         stop'om1xr(b) not defined'
5472         if(x.lt.xminDf)goto 111
5473 c.........Exponential
5474           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5475 c.........Linear
5476 c.........k=int(x*50.d0)
5477         histo(k)=histo(k)+1.d0
5478  111  continue
5479       do i=0,51
5480
5481 c.......Exponential
5482
5483         x=xminDf
5484         x1=xminDf
5485         x=x**(1.d0-dble(i)/51.d0)
5486         x1=x1**(1.d0-dble(i+1)/51.d0)
5487         if(i.eq.51)then
5488           x1=1.d0
5489           x=0.d0
5490         endif
5491         histo(i)=histo(i)/dble(n)/(x1-x)
5492
5493 c.......Linear
5494 c        histo(i)=histo(i)/dble(n)*51.d0
5495       enddo
5496
5497       return
5498       end
5499
5500
5501 c----------------------------------------------------------------------
5502
5503       subroutine xhistomran2(histo,xh,b)
5504
5505 c----------------------------------------------------------------------
5506 c.....Make Histogram of om1yr
5507 c----------------------------------------------------------------------
5508
5509       double precision histo(0:51),x,xh,dx,ymax   !,om1yr
5510       integer*4 n
5511
5512       ymax=-.5D0*dlog(xh)
5513       dx=ymax/25.d0
5514
5515       n=100000
5516       do i=0,50
5517         histo(i)=0.d0
5518       enddo
5519       do j=1,n
5520         if(mod(j,10000).eq.0)write(*,*)"y1",j,b
5521         x= 0 !  om1yr(xh,b)
5522         stop'om1yr(xh,b) not defined'
5523         k=int((x/ymax+1.d0)*25.d0)
5524 c.......write(*,*)x,k
5525         histo(k)=histo(k)+1.d0
5526       enddo
5527       do i=0,50
5528         histo(i)=histo(i)/dble(n)/dx
5529       enddo
5530
5531       return
5532       end
5533
5534
5535 cc----------------------------------------------------------------------
5536 c
5537 c      subroutine xhistomran6(histo,bx,by,bmax,del)
5538 c
5539 cc----------------------------------------------------------------------
5540 cc.....Make Histogram of b1 (impact parameter of vertex in Y and X)
5541 cc----------------------------------------------------------------------
5542 c
5543 c      double precision histo(0:51),dx
5544 c      integer*4 n
5545 c
5546 c      n=100000
5547 c      dx=dble(bmax)/50.d0
5548 c      do i=0,50
5549 c        histo(i)=0.d0
5550 c      enddo
5551 c      do j=1,n
5552 c        if(mod(j,10000).eq.0)write(*,*)"b1",j
5553 c        z=rangen()
5554 c        zp=rangen()
5555 c        bb1x=(bx+sqrt(-del*log(z))*cos(2.*3.14*zp))/2.
5556 c        bb1y=(by+sqrt(-del*log(z))*sin(2.*3.14*zp))/2.
5557 c        x=sqrt((bx-bb1x)*(bx-bb1x)+(by-bb1y)*(by-bb1y))
5558 c        k=int(x/bmax*50.)
5559 c        if(k.le.50)then
5560 c          histo(k)=histo(k)+1.d0
5561 c        else
5562 c          histo(51)=histo(51)+1.d0
5563 c        endif
5564 c      enddo
5565 c      do i=0,50
5566 c        histo(i)=histo(i)/dble(n)/dx
5567 c      enddo
5568 c
5569 c      return
5570 c      end
5571 c
5572
5573 c----------------------------------------------------------------------
5574
5575       subroutine xhistomran8(histo,xpr,xmr)
5576
5577 c----------------------------------------------------------------------
5578 c.....Make Histogram of om1xprk
5579 c----------------------------------------------------------------------
5580
5581       include 'epos.incpar'
5582
5583       double precision histo(0:51),x,x1,om1xprk,xpr,xmr
5584       integer*4 n
5585
5586
5587       n=100000
5588       do i=0,51
5589         histo(i)=0.d0
5590       enddo
5591       do 111 j=1,n
5592         if(mod(j,10000).eq.0)write(*,*)"x+",j,xmr
5593           x=om1xprk(1,xpr,xminDf,1)
5594         if(x.lt.xminDf)goto 111
5595 c.........Exponential
5596           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5597 c.........Linear
5598 c.........k=int(x*50.d0)
5599         histo(k)=histo(k)+1.d0
5600  111  continue
5601       do i=0,51
5602
5603 c.......Exponential
5604
5605         x=xminDf
5606         x1=xminDf
5607         x=x**(1.d0-dble(i)/51.d0)
5608         x1=x1**(1.d0-dble(i+1)/51.d0)
5609         if(i.eq.51)then
5610           x1=1.d0
5611           x=0.d0
5612         endif
5613         histo(i)=histo(i)/dble(n)/(x1-x)*xpr
5614
5615 c.......Linear
5616 c        histo(i)=histo(i)/dble(n)*51.d0
5617       enddo
5618
5619       return
5620       end
5621
5622 c----------------------------------------------------------------------
5623
5624       subroutine xhistomran9(histo,xp,xpr,xmr)
5625
5626 c----------------------------------------------------------------------
5627 c.....Make Histogram of om1xmrk
5628 c----------------------------------------------------------------------
5629
5630       include 'epos.incpar'
5631
5632       double precision histo(0:51),x,x1,om1xmrk,xp,xpr,xmr
5633       integer*4 n
5634
5635
5636       n=100000
5637       do i=0,51
5638         histo(i)=0.d0
5639       enddo
5640       do 111 j=1,n
5641         if(mod(j,10000).eq.0)write(*,*)"x-",j
5642           x=om1xmrk(1,xp,xpr,xminDf,1)
5643         if(x.lt.xminDf)goto 111
5644 c.........Exponential
5645           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5646 c.........Linear
5647 c.........k=int(x*50.d0)
5648         histo(k)=histo(k)+1.d0
5649  111  continue
5650       do i=0,51
5651
5652 c.......Exponential
5653
5654         x=xminDf
5655         x1=xminDf
5656         x=x**(1.d0-dble(i)/51.d0)
5657         x1=x1**(1.d0-dble(i+1)/51.d0)
5658         if(i.eq.51)then
5659           x1=1.d0
5660           x=0.d0
5661         endif
5662         histo(i)=histo(i)/dble(n)/(x1-x)*xmr
5663
5664 c.......Linear
5665 c        histo(i)=histo(i)/dble(n)*51.d0
5666       enddo
5667
5668       return
5669       end
5670
5671 c----------------------------------------------------------------------
5672
5673       subroutine xhistomran10(histo)
5674
5675 c----------------------------------------------------------------------
5676 c.....Make Histogram of om1xrk
5677 c----------------------------------------------------------------------
5678
5679       include 'epos.incpar'
5680
5681       double precision histo(0:51),x,x1,om1xrk
5682       integer*4 n
5683
5684
5685       n=100000
5686       do i=0,51
5687         histo(i)=0.d0
5688       enddo
5689       do 111 j=1,n
5690         if(mod(j,10000).eq.0)write(*,*)"xk",j
5691         x=om1xrk(1)
5692         if(x.lt.xminDf)goto 111
5693 c.........Exponential
5694           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5695 c.........Linear
5696 c.........k=int(x*50.d0)
5697         histo(k)=histo(k)+1.d0
5698  111  continue
5699       do i=0,51
5700
5701 c.......Exponential
5702
5703         x=xminDf
5704         x1=xminDf
5705         x=x**(1.d0-dble(i)/51.d0)
5706         x1=x1**(1.d0-dble(i+1)/51.d0)
5707         if(i.eq.51)then
5708           x1=1.d0
5709           x=0.d0
5710         endif
5711         histo(i)=histo(i)/dble(n)/(x1-x)
5712
5713 c.......Linear
5714 c        histo(i)=histo(i)/dble(n)*51.d0
5715       enddo
5716
5717       return
5718       end
5719
5720
5721 c----------------------------------------------------------------------
5722
5723       subroutine xhistomran11(histo,xh)
5724
5725 c----------------------------------------------------------------------
5726 c.....Make Histogram of om1yrk
5727 c----------------------------------------------------------------------
5728
5729       double precision histo(0:51),x,xh,dx,om1yrk,ymax
5730       integer*4 n
5731
5732       ymax=-.5D0*dlog(xh)
5733       dx=ymax/25.d0
5734
5735       n=100000
5736       do i=0,50
5737         histo(i)=0.d0
5738       enddo
5739       do j=1,n
5740         if(mod(j,10000).eq.0)write(*,*)"yk",j
5741         x=om1yrk(xh)
5742         k=int((x/ymax+1.d0)*25.d0)
5743 c.......write(*,*)x,k
5744         histo(k)=histo(k)+1.d0
5745       enddo
5746       do i=0,50
5747         histo(i)=histo(i)/dble(n)/dx
5748       enddo
5749
5750       return
5751       end
5752
5753 c----------------------------------------------------------------------
5754
5755       subroutine xranptg(histo,xcut,xfact,xadd)
5756
5757 c----------------------------------------------------------------------
5758 c.....Make Histogram of random distribution
5759 c----------------------------------------------------------------------
5760
5761       include 'epos.incpar'
5762
5763       parameter (nptg1=501)      !number of point for the graphs
5764       common /cranpt/conv
5765       double precision histo(0:nptg1)
5766       integer*4 n
5767
5768
5769       n=100000
5770       do i=0,nptg1
5771         histo(i)=0.d0
5772       enddo
5773       do j=1,n
5774         if(mod(j,10000).eq.0)write(*,*)"ptg",j
5775 c .........exp(-x**2)
5776  12   x=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
5777
5778       if(xcut.gt.0.)then
5779         if(rangen().lt.x/xcut)goto 12
5780       endif
5781       x=x*xfact+xadd
5782 c.........Exponential
5783 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5784 c.........Linear
5785         k=int(x/conv)
5786         k=min(k,nptg1)
5787         histo(k)=histo(k)+1.d0
5788       enddo
5789       do i=0,nptg1
5790
5791 c.......Exponential
5792
5793 c        x=xminDf
5794 c        x1=xminDf
5795 c        x=x**(1.d0-dble(i)/51.d0)
5796 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5797 c        if(i.eq.51)then
5798 c          x1=1.d0
5799 c          x=0.d0
5800 c        endif
5801 c        histo(i)=histo(i)/dble(n)/(x1-x)
5802
5803 c.......Linear
5804         histo(i)=histo(i)/dble(n)*float(nptg1)
5805       enddo
5806
5807       return
5808       end
5809
5810 c----------------------------------------------------------------------
5811
5812       subroutine xranpte(histo,xcut,xfact,xadd)
5813
5814 c----------------------------------------------------------------------
5815 c.....Make Histogram of random distribution
5816 c----------------------------------------------------------------------
5817
5818       include 'epos.incpar'
5819
5820       parameter (nptg1=501)      !number of point for the graphs
5821       common /cranpt/conv
5822       double precision histo(0:nptg1)
5823       integer*4 n
5824
5825
5826       n=100000
5827       do i=0,nptg1
5828         histo(i)=0.d0
5829       enddo
5830       do j=1,n
5831         if(mod(j,10000).eq.0)write(*,*)"pte",j
5832 c .........exp(-x)
5833   12  xmx=50
5834       r=2.
5835       x=0.
5836       do while (r.gt.1.)
5837   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
5838         if(x.eq.0.)goto11
5839         r=rangen()  /  ( exp(-x)*(1+x**2) )
5840       enddo
5841       x=x/2.
5842
5843       if(xcut.gt.0.)then
5844         if(rangen().lt.x/xcut)goto 12
5845       endif
5846       x=x*xfact+xadd
5847 c.........Exponential
5848 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5849 c.........Linear
5850         k=int(x/conv)
5851         k=min(k,nptg1)
5852         histo(k)=histo(k)+1.d0
5853       enddo
5854       do i=0,nptg1
5855
5856 c.......Exponential
5857
5858 c        x=xminDf
5859 c        x1=xminDf
5860 c        x=x**(1.d0-dble(i)/51.d0)
5861 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5862 c        if(i.eq.51)then
5863 c          x1=1.d0
5864 c          x=0.d0
5865 c        endif
5866 c        histo(i)=histo(i)/dble(n)/(x1-x)
5867
5868 c.......Linear
5869         histo(i)=histo(i)/dble(n)*float(nptg1)
5870       enddo
5871
5872       return
5873       end
5874
5875 c----------------------------------------------------------------------
5876
5877       subroutine xranpts(histo,xcut,xfact,xadd)
5878
5879 c----------------------------------------------------------------------
5880 c.....Make Histogram of random distribution
5881 c----------------------------------------------------------------------
5882
5883       include 'epos.incpar'
5884
5885       parameter (nptg1=501)      !number of point for the graphs
5886       common /cranpt/conv
5887       double precision histo(0:nptg1)
5888       integer*4 n
5889
5890
5891       n=100000
5892       do i=0,nptg1
5893         histo(i)=0.d0
5894       enddo
5895       do j=1,n
5896         if(mod(j,10000).eq.0)write(*,*)"pts",j
5897 c .........exp(-sqrt(x))
5898  12   xmx=500
5899       r=2.
5900       x=0.
5901       do while (r.gt.1.)
5902         x=sqrt(exp(rangen()*log(1+xmx**2))-1)
5903         r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
5904       enddo
5905       x=x/20.
5906
5907       if(xcut.gt.0.)then
5908         if(rangen().lt.x/xcut)goto 12
5909       endif
5910       x=x*xfact+xadd
5911 c.........Exponential
5912 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5913 c.........Linear
5914         k=int(x/conv)
5915         k=min(k,nptg1)
5916         histo(k)=histo(k)+1.d0
5917       enddo
5918       do i=0,nptg1
5919
5920 c.......Exponential
5921
5922 c        x=xminDf
5923 c        x1=xminDf
5924 c        x=x**(1.d0-dble(i)/51.d0)
5925 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5926 c        if(i.eq.51)then
5927 c          x1=1.d0
5928 c          x=0.d0
5929 c        endif
5930 c        histo(i)=histo(i)/dble(n)/(x1-x)
5931
5932 c.......Linear
5933         histo(i)=histo(i)/dble(n)*float(nptg1)
5934       enddo
5935
5936       return
5937       end
5938
5939 c----------------------------------------------------------------------
5940
5941       subroutine xranptc(histo,xcut,xfact,xadd)
5942
5943 c----------------------------------------------------------------------
5944 c.....Make Histogram of random distribution
5945 c----------------------------------------------------------------------
5946
5947       include 'epos.incpar'
5948
5949       parameter (nptg1=501)      !number of point for the graphs
5950       common /cranpt/conv
5951       double precision histo(0:nptg1)
5952       integer*4 n
5953
5954
5955       n=100000
5956       do i=0,nptg1
5957         histo(i)=0.d0
5958       enddo
5959       do j=1,n
5960         if(mod(j,10000).eq.0)write(*,*)"ptc",j
5961
5962         x=ranptcut(xcut)*xfact+xadd
5963 c.........Exponential
5964 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5965 c.........Linear
5966         k=int(x/conv)
5967         k=min(k,nptg1)
5968         histo(k)=histo(k)+1.d0
5969       enddo
5970       do i=0,nptg1
5971
5972 c.......Exponential
5973
5974 c        x=xminDf
5975 c        x1=xminDf
5976 c        x=x**(1.d0-dble(i)/51.d0)
5977 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5978 c        if(i.eq.51)then
5979 c          x1=1.d0
5980 c          x=0.d0
5981 c        endif
5982 c        histo(i)=histo(i)/dble(n)/(x1-x)
5983
5984 c.......Linear
5985         histo(i)=histo(i)/dble(n)*float(nptg1)
5986       enddo
5987
5988       return
5989       end
5990
5991
5992