EPOS/epos167/epos-xpr-165.f

   1 c----------------------------------------------------------------------
   2       double precision function xDfit(i1,i2,s,xp,xm,b)
   3 c----------------------------------------------------------------------
   4
   5       include 'epos.inc'
   6       include 'epos.incsem'
   7
   8
   9
  10       xDfit=0.d0
  11       do i=i1,i2
  12         call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
  13         if(i.lt.2)then
  14           corp=alppar-epsp
  15           cort=alppar-epst
  16         else
  17           corp=alppar
  18           cort=alppar
  19         endif
  20         xDfit=xDfit+dble(alp*xp**(bet+corp)*xm**(betp+cort)*s**(-epss))
  21       enddo
  22       return
  23       end
  24
  25
  26 c----------------------------------------------------------------------
  27       subroutine xFitD1
  28 c----------------------------------------------------------------------
  29
  30       include 'epos.inc'
  31       include 'epos.incsem'
  32       include 'epos.incpar'
  33       double precision x,y,Dsoftshval,om51p,xminr,tmp,xtmp,xDfit
  34       character chenergy*12
  35
  36       nptg=50                   !number of point for the graphs
  37       iii=nint(xpar1)
  38       biniDf=xpar2                   !value of biniDf (impact parameter)
  39       y=dble(xpar3)                    !value of y (rapidity)
  40       xtmp=xmaxDf
  41       xmaxDf=dexp(-2.d0*y)
  42
  43       chenergy='E=          '
  44       if (engy.ge.10000.) then
  45         write(chenergy(4:8),'(I5)')int(engy)
  46         ke=10
  47       elseif (engy.ge.1000.) then
  48         write(chenergy(4:7),'(I4)')int(engy)
  49         ke=9
  50       elseif (engy.ge.100.) then
  51         write(chenergy(4:6),'(I3)')int(engy)
  52         ke=8
  53       elseif (engy.ge.10.) then
  54         write(chenergy(4:5),'(I2)')int(engy)
  55         ke=7
  56       else
  57         write(chenergy(4:4),'(I1)')int(engy)
  58         ke=6
  59       endif
  60       chenergy(ke:ke+2)='GeV'
  61
  62       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
  63
  64       if(iii/10.eq.1)then !...................................................
  65
  66       write(ifhi,'(a)')'!----------------------------------------------'
  67       write(ifhi,'(a)')'!     D exact all      (blue)      '
  68       write(ifhi,'(a)')'!----------------------------------------------'
  69
  70       write(ifhi,'(a)')'openhisto name DExact-'//chenergy(4:ke-2)
  71       write(ifhi,'(a)')       'htyp lbu'
  72       write(ifhi,'(a)')       'xmod log ymod log'
  73       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
  74       write(ifhi,'(a)')       'yrange .01 auto'
  75       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
  76       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
  77       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact" '
  78       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
  79         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
  80         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
  81       write(ifhi,'(a)')       'array 2'
  82
  83       do i=0,nptg
  84         x=xminr
  85         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
  86         tmp=Dsoftshval(real(x)*smaxDf,x,y,biniDf,0)
  87         write(ifhi,*) x,tmp
  88       enddo
  89
  90       write(ifhi,'(a)')    '  endarray'
  91       write(ifhi,'(a)')    'closehisto plot 0-'
  92
  93
  94       write(ifhi,'(a)')'!----------------------------------------------'
  95       write(ifhi,'(a)')'!     D exact soft      (red dot)      '
  96       write(ifhi,'(a)')'!----------------------------------------------'
  97
  98       write(ifhi,'(a)')'openhisto name DExactSoft-'//chenergy(4:ke-2)
  99       write(ifhi,'(a)')       'htyp lro'
 100       write(ifhi,'(a)')       'xmod log ymod log'
 101       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 102       write(ifhi,'(a)')       'yrange .01 auto'
 103       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 104       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 105       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 106       if (xpar8.eq.1.) then
 107         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 108         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 109       endif
 110       write(ifhi,'(a)')       'array 2'
 111
 112       do i=0,nptg
 113         x=xminr
 114         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 115         write(ifhi,*) x,2.d0*om51p(real(x)*smaxDf,x,y,biniDf,0)
 116      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
 117       enddo
 118
 119       write(ifhi,'(a)')    '  endarray'
 120       write(ifhi,'(a)')    'closehisto plot 0-'
 121
 122
 123       write(ifhi,'(a)')'!---------------------------------------------'
 124       write(ifhi,'(a)')'!     D exact sea-sea      (yellow-dot)     '
 125       write(ifhi,'(a)')'!---------------------------------------------'
 126
 127       write(ifhi,'(a)')'openhisto name DExactSemi-'//chenergy(4:ke-2)
 128       write(ifhi,'(a)')       'htyp lyo'
 129       write(ifhi,'(a)')       'xmod log ymod log'
 130       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 131       write(ifhi,'(a)')       'yrange .01 auto'
 132       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 133       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 134       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 135       if (xpar8.eq.1.) then
 136         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 137         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 138       endif
 139       write(ifhi,'(a)')       'array 2'
 140
 141       do i=0,nptg
 142         x=xminr
 143         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 144         write(ifhi,*) x,2.d0*om51p(real(x)*smaxDf,x,y,biniDf,1)
 145      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
 146       enddo
 147
 148       write(ifhi,'(a)')    '  endarray'
 149       write(ifhi,'(a)')    'closehisto plot 0-'
 150
 151       write(ifhi,'(a)')'!---------------------------------------------'
 152       write(ifhi,'(a)')'!     D exact semi      (blue dot)          '
 153       write(ifhi,'(a)')'!---------------------------------------------'
 154
 155       write(ifhi,'(a)')'openhisto name DExactVal-'//chenergy(4:ke-2)
 156       write(ifhi,'(a)')       'htyp lbo'
 157       write(ifhi,'(a)')       'xmod log ymod log'
 158       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 159       write(ifhi,'(a)')       'yrange .01 auto'
 160       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 161       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 162       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 163       if (xpar8.eq.1.) then
 164         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 165         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 166       endif
 167       write(ifhi,'(a)')       'array 2'
 168
 169       do i=0,nptg
 170         x=xminr
 171         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 172         write(ifhi,*) x,2.d0*(om51p(real(x)*smaxDf,x,y,biniDf,2)
 173      &       +om51p(real(x)*smaxDf,x,y,biniDf,1)
 174      &       +om51p(real(x)*smaxDf,x,y,biniDf,3)
 175      &       +om51p(real(x)*smaxDf,x,y,biniDf,4))
 176      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
 177       enddo
 178
 179       write(ifhi,'(a)')    '  endarray'
 180       write(ifhi,'(a)')    'closehisto plot 0'
 181
 182       endif !................................................................
 183       if(mod(iii,10).eq.1)then !.............................................
 184
 185       write(ifhi,'(a)')'!---------------------------------------------'
 186       write(ifhi,'(a)')'!     D exact all      (blue)              '
 187       write(ifhi,'(a)')'!---------------------------------------------'
 188
 189       write(ifhi,'(a)')'openhisto name DExact-'//chenergy(4:ke-2)
 190       write(ifhi,'(a)')       'htyp lbu'
 191       write(ifhi,'(a)')       'xmod log ymod log'
 192       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 193       write(ifhi,'(a)')       'yrange .01 auto'
 194       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 195       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 196       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact+fit" '
 197       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 198       if (xpar8.eq.1.) then
 199         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 200         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 201       endif
 202       write(ifhi,'(a)')       'array 2'
 203
 204       do i=0,nptg
 205         x=xminr
 206         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 207         tmp=Dsoftshval(real(x)*smaxDf,x,y,biniDf,0)
 208         write(ifhi,*) x,tmp
 209       enddo
 210
 211       write(ifhi,'(a)')    '  endarray'
 212       write(ifhi,'(a)')    'closehisto plot 0-'
 213
 214       write(ifhi,'(a)')'!---------------------------------------------'
 215       write(ifhi,'(a)')'!     D exact all      (green)              '
 216       write(ifhi,'(a)')'!---------------------------------------------'
 217
 218       write(ifhi,'(a)')'openhisto name DExact-f-'//chenergy(4:ke-2)
 219       write(ifhi,'(a)')       'htyp lga'
 220       write(ifhi,'(a)')       'xmod log ymod log'
 221       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 222       write(ifhi,'(a)')       'yrange .01 auto'
 223       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 224       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 225       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact+fit" '
 226       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 227       if (xpar8.eq.1.) then
 228         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 229         write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 230       endif
 231       write(ifhi,'(a)')       'array 2'
 232
 233       do i=0,nptg
 234         x=xminr
 235         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 236         tmp=Dsoftshval(real(x)*smaxDf,x,y,biniDf,-1)
 237         write(ifhi,*) x,tmp
 238       enddo
 239
 240       write(ifhi,'(a)')    '  endarray'
 241       write(ifhi,'(a)')    'closehisto plot 0-'
 242
 243       write(ifhi,'(a)')'!---------------------------------------------'
 244       write(ifhi,'(a)')'!     fit soft      (red dot)             '
 245       write(ifhi,'(a)')'!---------------------------------------------'
 246
 247       write(ifhi,'(a)')       'openhisto'
 248       write(ifhi,'(a)')       'htyp lro'
 249       write(ifhi,'(a)')       'xmod log ymod log'
 250       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 251       write(ifhi,'(a)')       'yrange .01 auto'
 252       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 253       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 254       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 255       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 256       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 257       write(ifhi,'(a)')       'array 2'
 258
 259
 260       do i=0,nptg
 261         x=xminr
 262         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 263         xp=sqrt(real(x))*exp(real(y))
 264         xm=sqrt(real(x))*exp(-real(y))
 265         write(ifhi,*) x,xDfit(0,0,smaxDf,xp,xm,biniDf)
 266       enddo
 267       write(ifhi,'(a)')    '  endarray'
 268       write(ifhi,'(a)')    'closehisto plot 0-'
 269
 270       write(ifhi,'(a)')'!---------------------------------------------'
 271       write(ifhi,'(a)')'!     fit semi      (blue dot)              '
 272       write(ifhi,'(a)')'!---------------------------------------------'
 273
 274       write(ifhi,'(a)')       'openhisto'
 275       write(ifhi,'(a)')       'htyp lbo'
 276       write(ifhi,'(a)')       'xmod log ymod log'
 277       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 278       write(ifhi,'(a)')       'yrange .01 auto'
 279       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 280       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,,s,b)" '
 281       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 282       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 283       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 284       write(ifhi,'(a)')       'array 2'
 285
 286
 287       do i=0,nptg
 288         x=xminr
 289         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 290         xp=sqrt(real(x))*exp(real(y))
 291         xm=sqrt(real(x))*exp(-real(y))
 292         write(ifhi,*) x,xDfit(1,1,smaxDf,xp,xm,biniDf)
 293       enddo
 294       write(ifhi,'(a)')    '  endarray'
 295       write(ifhi,'(a)')    'closehisto plot 0-'
 296
 297       write(ifhi,'(a)')'!---------------------------------------------'
 298       write(ifhi,'(a)')'!     fit all      (red)      '
 299       write(ifhi,'(a)')'!---------------------------------------------'
 300
 301       write(ifhi,'(a)')       'openhisto'
 302       write(ifhi,'(a)')       'htyp lru'
 303       write(ifhi,'(a)')       'xmod log ymod log'
 304       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 305       write(ifhi,'(a)')       'yrange .01 auto'
 306       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 307       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,,s,b)" '
 308       write(ifhi,'(3a)')          'text 0.05 0.9 "',chenergy,'"'
 309       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.8 "b=',biniDf,' fm"'
 310       write(ifhi,'(a,f5.2,a)')  'text 0.05 0.7 "y=',y,'"'
 311       write(ifhi,'(a)')       'array 2'
 312
 313       do i=0,nptg
 314         x=xminr
 315         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 316         xp=sqrt(real(x))*exp(real(y))
 317         xm=sqrt(real(x))*exp(-real(y))
 318         write(ifhi,*) x,xDfit(0,1,smaxDf,xp,xm,biniDf)
 319       enddo
 320       write(ifhi,'(a)')    '  endarray'
 321       write(ifhi,'(a)')    'closehisto plot 0'
 322
 323       endif !...............................................................
 324
 325       xmaxDf=xtmp
 326
 327       end
 328
 329 c----------------------------------------------------------------------
 330       subroutine xFitD2
 331 c----------------------------------------------------------------------
 332
 333       include 'epos.inc'
 334       include 'epos.incsem'
 335       include 'epos.incpar'
 336
 337       double precision x,om51p,xDfit,z(0:200),xminr,y,xtmp,om,om5,om51
 338 c     & ,omYuncut,omNpuncut
 339       character chenergy*12,chf*3,texte*15,textb*17,texty*17
 340
 341       nptg=30                  !number of point for the graphs
 342       biniDf=xpar2                 !value of biniDf (impact parameter)
 343       y=dble(xpar3)                 !value of y (rapidity)
 344       jj1=nint(xpar4)
 345       jj2=nint(xpar5)
 346       if(jj1.ne.1.and.jj1.ne.2)jj1=3
 347       if(jj2.ne.1.and.jj2.ne.2)jj2=3
 348       xtmp=xmaxDf
 349       xmaxDf=dexp(-2.d0*y)
 350
 351       chenergy='E=          '
 352       if (engy.ge.10000.) then
 353         write(chenergy(4:8),'(I5)')int(engy)
 354         ke=10
 355       elseif (engy.ge.1000.) then
 356         write(chenergy(4:7),'(I4)')int(engy)
 357         ke=9
 358       elseif (engy.ge.100.) then
 359         write(chenergy(4:6),'(I3)')int(engy)
 360         ke=8
 361       elseif (engy.ge.10.) then
 362         write(chenergy(4:5),'(I2)')int(engy)
 363         ke=7
 364       else
 365         write(chenergy(4:4),'(I1)')int(engy)
 366         ke=6
 367       endif
 368       chenergy(ke:ke+2)='GeV'
 369
 370       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
 371
 372          do jj=jj1,jj2
 373
 374       if(jj.eq.1)chf='  D'
 375       if(jj.eq.2)chf='  G'
 376       if(jj.eq.3)chf='FFG'
 377       texte='text 0.05 0.9 "'
 378       textb='text 0.05 0.8 "b='
 379       texty='text 0.05 0.7 "y='
 380       if(jj.eq.2)texty='text 0.15 0.7 "y='
 381       if(jj.eq.3)texte='text 0.05 0.3 "'
 382       if(jj.eq.3)textb='text 0.05 0.2 "b='
 383       if(jj.eq.3)texty='text 0.05 0.1 "y='
 384
 385       write(ifhi,'(a)')'!---------------------------------------------'
 386       write(ifhi,'(a)')'!     '//chf//' exact all      (green)         '
 387       write(ifhi,'(a)')'!---------------------------------------------'
 388
 389       write(ifhi,'(a)')'openhisto name '//chf//'ExaI-'//chenergy(4:ke-2)
 390       write(ifhi,'(a)')       'htyp lga'
 391       write(ifhi,'(a)')       'xmod log ymod log'
 392       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 393       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 394       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
 395       write(ifhi,'(a,a)')    'text 0.65 0.9 "exact+fit" '
 396       write(ifhi,'(3a)')        texte,chenergy,'"'
 397       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
 398       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
 399       write(ifhi,'(a)')       'array 2'
 400       do i=0,nptg
 401         x=xminr
 402         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 403         xp=sqrt(real(x))*exp(real(y))
 404         xm=sqrt(real(x))*exp(-real(y))
 405         om=0
 406         do j=0,4
 407          om=om+om51p(engy**2*real(x),x,y,biniDf,j)
 408         enddo
 409         om=2.d0*om
 410         if(jj.eq.1)om=om/(x**dble(-alppar))
 411         if(jj.eq.3)om=om
 412      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
 413         write(ifhi,*) x,om
 414       enddo
 415       write(ifhi,'(a)')    '  endarray'
 416       write(ifhi,'(a)')    'closehisto plot 0-'
 417
 418       write(ifhi,'(a)')'!---------------------------------------------'
 419       write(ifhi,'(a)')'!     '//chf//' exact all +diff (blue)        '
 420       write(ifhi,'(a)')'!---------------------------------------------'
 421
 422       write(ifhi,'(a)')'openhisto name '//chf//'ExaD-'//chenergy(4:ke-2)
 423       write(ifhi,'(a)')       'htyp lbu'
 424       write(ifhi,'(a)')       'xmod log ymod log'
 425       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 426       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 427       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
 428       write(ifhi,'(3a)')        texte,chenergy,'"'
 429       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
 430       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
 431       write(ifhi,'(a)')       'array 2'
 432       do i=0,nptg
 433         x=xminr
 434         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 435         xp=sqrt(real(x))*exp(real(y))
 436         xm=sqrt(real(x))*exp(-real(y))
 437         om=0
 438         do j=0,4
 439          om=om+om51p(engy**2*real(x),x,y,biniDf,j)
 440         enddo
 441         om5=om51(x,y,biniDf,5,5)
 442         om=2.d0*(om+om5)
 443         if(jj.eq.1)om=om/(x**dble(-alppar))
 444         if(jj.eq.3)om=om
 445      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
 446         write(ifhi,*) x,om
 447       enddo
 448       write(ifhi,'(a)')    '  endarray'
 449       write(ifhi,'(a)')    'closehisto plot 0-'
 450
 451       write(ifhi,'(a)')'!---------------------------------------------'
 452       write(ifhi,'(a)')'!     '//chf//' param all      (red)          '
 453       write(ifhi,'(a)')'!---------------------------------------------'
 454
 455       write(ifhi,'(a)')'openhisto name '//chf//'Par-'//chenergy(4:ke-2)
 456       write(ifhi,'(a)')       'htyp lru'
 457         write(ifhi,'(a)')       'xmod log ymod log'
 458       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
 459       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
 460       write(ifhi,'(a)')      'text 0 0 "yaxis '//chf//'(x+,x-,s,b)" '
 461       write(ifhi,'(3a)')        texte,chenergy,'"'
 462       write(ifhi,'(a,f5.2,a)')  textb,biniDf,' fm"'
 463       write(ifhi,'(a,f5.2,a)')  texty,y,'"'
 464       write(ifhi,'(a)')       'array 2'
 465       imax=idxD1
 466       if(iomega.eq.2)imax=1
 467       do i=0,nptg
 468         x=xminr
 469         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
 470         xp=sqrt(real(x))*exp(real(y))
 471         xm=sqrt(real(x))*exp(-real(y))
 472           z(i)=xDfit(0,imax,engy**2,xp,xm,biniDf)
 473          if(jj.ge.2)z(i)=z(i)*(x**dble(-alppar))
 474          if(jj.eq.3)z(i)=z(i)
 475      &               *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro)
 476          write(ifhi,*) x,z(i)
 477        enddo
 478       write(ifhi,'(a)')    '  endarray'
 479       write(ifhi,'(a)')    'closehisto plot 0'
 480
 481          enddo
 482
 483       xmaxDf=xtmp
 484
 485       end
 486
 487 c----------------------------------------------------------------------
 488       subroutine xbExaD
 489 c----------------------------------------------------------------------
 490
 491       include 'epos.inc'
 492       include 'epos.incsem'
 493       include 'epos.incpar'
 494
 495       double precision x,y,Dsoftshval,om51p,z,xDfit!,omNpuncut
 496
 497
 498       nptg=50     !number of point for the graphs
 499       bmax=xpar2
 500       bmax=max(0.1,bmax)
 501              !value max of b (impact parameter)
 502       y=dble(xpar3)                    !value of y (rapidity)
 503       x=dble(xpar4)
 504       if(x.eq.0.d0)x=1.d0
 505
 506       if (engy.ge.10.) then
 507        if (engy.ge.100.) then
 508         if (engy.ge.1000.) then
 509          if (engy.ge.10000.) then
 510       write(ifhi,'(a,I5)')  'openhisto name DExactb-',int(engy)
 511          else
 512       write(ifhi,'(a,I4)')  'openhisto name DExactb-',int(engy)
 513         endif
 514         else
 515       write(ifhi,'(a,I3)')  'openhisto name DExactb-',int(engy)
 516        endif
 517        else
 518       write(ifhi,'(a,I2)')  'openhisto name DExactb-',int(engy)
 519       endif
 520       else
 521       write(ifhi,'(a,I1)')  'openhisto name DExactb-',int(engy)
 522       endif
 523       write(ifhi,'(a)')       'htyp lbu'
 524       write(ifhi,'(a)')       'xmod lin ymod lin'
 525       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 526       write(ifhi,'(a)')       'yrange .00001 auto'
 527 c      write(ifhi,'(a)')       'yrange auto auto'
 528       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 529       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 530       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 531       if (xpar8.eq.1.) then
 532       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 533       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 534       endif
 535       write(ifhi,'(a)')       'array 2'
 536
 537       do i=0,nptg
 538         b=-bmax+2.*real(i)/real(nptg)*bmax
 539         write(ifhi,*) b,Dsoftshval(real(x)*smaxDf,x,y,b,0)
 540       enddo
 541
 542       write(ifhi,'(a)')    '  endarray'
 543       write(ifhi,'(a)')    'closehisto plot 0-'
 544
 545 c**********************************************************************
 546
 547       if (engy.ge.10.) then
 548        if (engy.ge.100.) then
 549         if (engy.ge.1000.) then
 550          if (engy.ge.10000.) then
 551       write(ifhi,'(a,I5)')       'openhisto name DParamb-',int(engy)
 552          else
 553       write(ifhi,'(a,I4)')       'openhisto name DParamb-',int(engy)
 554         endif
 555         else
 556       write(ifhi,'(a,I3)')       'openhisto name DParamb-',int(engy)
 557        endif
 558        else
 559       write(ifhi,'(a,I2)')       'openhisto name DParamb-',int(engy)
 560       endif
 561       else
 562       write(ifhi,'(a,I1)')       'openhisto name DParamb-',int(engy)
 563       endif
 564       write(ifhi,'(a)')       'htyp lrd'
 565       write(ifhi,'(a)')       'xmod lin ymod lin'
 566       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 567 c      write(ifhi,'(a)')       'yrange -.01 auto'
 568       write(ifhi,'(a)')       'yrange auto auto'
 569       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 570       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 571       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 572       if (xpar8.eq.1.) then
 573       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 574       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 575       endif
 576       write(ifhi,'(a)')       'array 2'
 577
 578       imax=idxD1
 579       if(iomega.eq.2)imax=1
 580       do i=0,nptg
 581         b=-bmax+2.*real(i)/real(nptg)*bmax
 582           z=xDfit(0,imax,smaxDf,
 583      &       real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),b)
 584         write(ifhi,*) b,z
 585       enddo
 586
 587       write(ifhi,'(a)')    '  endarray'
 588       write(ifhi,'(a)')    'closehisto plot 0-'
 589
 590
 591 c**********************************************************************
 592
 593       if (engy.ge.10.) then
 594        if (engy.ge.100.) then
 595         if (engy.ge.1000.) then
 596          if (engy.ge.10000.) then
 597       write(ifhi,'(a,I5)')  'openhisto name DExactSoftb-',int(engy)
 598          else
 599       write(ifhi,'(a,I4)')  'openhisto name DExactSoftb-',int(engy)
 600         endif
 601         else
 602       write(ifhi,'(a,I3)')  'openhisto name DExactSoftb-',int(engy)
 603        endif
 604        else
 605       write(ifhi,'(a,I2)')  'openhisto name DExactSoftb-',int(engy)
 606       endif
 607       else
 608       write(ifhi,'(a,I1)')  'openhisto name DExactSoftb-',int(engy)
 609       endif
 610       write(ifhi,'(a)')       'htyp pfc'
 611       write(ifhi,'(a)')       'xmod lin ymod lin'
 612       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 613 c      write(ifhi,'(a)')       'yrange -.01 auto'
 614       write(ifhi,'(a)')       'yrange auto auto'
 615       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 616       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 617       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 618       if (xpar8.eq.1.) then
 619       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 620       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 621       endif
 622       write(ifhi,'(a)')       'array 2'
 623
 624       do i=0,nptg
 625         b=-bmax+2.*real(i)/real(nptg)*bmax
 626         write(ifhi,*) b,2.d0*om51p(real(x)*smaxDf,x,y,b,0)
 627      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
 628       enddo
 629
 630       write(ifhi,'(a)')    '  endarray'
 631       write(ifhi,'(a)')    'closehisto plot 0-'
 632
 633
 634 c**********************************************************************
 635
 636       if (engy.ge.10.) then
 637        if (engy.ge.100.) then
 638         if (engy.ge.1000.) then
 639          if (engy.ge.10000.) then
 640       write(ifhi,'(a,I5)')  'openhisto name DExactSemib-',int(engy)
 641          else
 642       write(ifhi,'(a,I4)')  'openhisto name DExactSemib-',int(engy)
 643         endif
 644         else
 645       write(ifhi,'(a,I3)')  'openhisto name DExactSemib-',int(engy)
 646        endif
 647        else
 648       write(ifhi,'(a,I2)')  'openhisto name DExactSemib-',int(engy)
 649       endif
 650       else
 651       write(ifhi,'(a,I1)')  'openhisto name DExactSemib-',int(engy)
 652       endif
 653       write(ifhi,'(a)')       'htyp pft'
 654       write(ifhi,'(a)')       'xmod lin ymod lin'
 655       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 656 c      write(ifhi,'(a)')       'yrange -.01 auto'
 657       write(ifhi,'(a)')       'yrange auto auto'
 658       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 659       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 660       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 661       if (xpar8.eq.1.) then
 662       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 663       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 664       endif
 665       write(ifhi,'(a)')       'array 2'
 666
 667       do i=0,nptg
 668         b=-bmax+2.*real(i)/real(nptg)*bmax
 669         write(ifhi,*) b,2.d0*om51p(real(x)*smaxDf,x,y,b,1)
 670      &       /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
 671       enddo
 672
 673       write(ifhi,'(a)')    '  endarray'
 674       write(ifhi,'(a)')    'closehisto plot 0-'
 675
 676
 677 c**********************************************************************
 678
 679       if (engy.ge.10.) then
 680        if (engy.ge.100.) then
 681         if (engy.ge.1000.) then
 682          if (engy.ge.10000.) then
 683       write(ifhi,'(a,I5)')  'openhisto name DExactValb-',int(engy)
 684          else
 685       write(ifhi,'(a,I4)')  'openhisto name DExactValb-',int(engy)
 686         endif
 687         else
 688       write(ifhi,'(a,I3)')  'openhisto name DExactValb-',int(engy)
 689        endif
 690        else
 691       write(ifhi,'(a,I2)')  'openhisto name DExactValb-',int(engy)
 692       endif
 693       else
 694       write(ifhi,'(a,I1)')  'openhisto name DExactValb-',int(engy)
 695       endif
 696       write(ifhi,'(a)')       'htyp pfs'
 697       write(ifhi,'(a)')       'xmod lin ymod lin'
 698       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 699 c      write(ifhi,'(a)')       'yrange -.01 auto'
 700       write(ifhi,'(a)')       'yrange auto auto'
 701       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 702       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 703       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 704       if (xpar8.eq.1.) then
 705       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 706       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 707       endif
 708       write(ifhi,'(a)')       'array 2'
 709
 710       do i=0,nptg
 711         b=-bmax+2.*real(i)/real(nptg)*bmax
 712         write(ifhi,*) b,2.d0*(om51p(real(x)*smaxDf,x,y,b,2)+
 713      &  om51p(real(x)*smaxDf,x,y,b,3)+om51p(real(x)*smaxDf,x,y,b,4))
 714      &           /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
 715       enddo
 716
 717       write(ifhi,'(a)')    '  endarray'
 718
 719       write(ifhi,'(a)')    'closehisto plot 0'
 720
 721
 722       end
 723
 724
 725 c----------------------------------------------------------------------
 726       subroutine xbnExaD
 727 c----------------------------------------------------------------------
 728
 729       include 'epos.inc'
 730       include 'epos.incpar'
 731
 732       double precision x,y,Dsoftshval,z,xDfit,Dint
 733
 734
 735
 736       nptg=50     !number of point for the graphs
 737       bmax=xpar2
 738       bmax=max(0.1,bmax)
 739       y=dble(xpar2)                   !value of y (rapidity)
 740       x=dble(xpar3)
 741
 742       if (engy.ge.10.) then
 743        if (engy.ge.100.) then
 744         if (engy.ge.1000.) then
 745          if (engy.ge.10000.) then
 746       write(ifhi,'(a,I5)')  'openhisto name DExactbn-',int(engy)
 747          else
 748       write(ifhi,'(a,I4)')  'openhisto name DExactbn-',int(engy)
 749         endif
 750         else
 751       write(ifhi,'(a,I3)')  'openhisto name DExactbn-',int(engy)
 752        endif
 753        else
 754       write(ifhi,'(a,I2)')  'openhisto name DExactbn-',int(engy)
 755       endif
 756       else
 757       write(ifhi,'(a,I1)')  'openhisto name DExactbn-',int(engy)
 758       endif
 759       write(ifhi,'(a)')       'htyp lbu'
 760       write(ifhi,'(a)')       'xmod lin ymod lin'
 761       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 762 c      write(ifhi,'(a)')       'yrange -.01 auto'
 763       write(ifhi,'(a)')       'yrange auto auto'
 764       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 765       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 766       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 767       if (xpar8.eq.1.) then
 768       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 769       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 770       endif
 771       write(ifhi,'(a)')       'array 2'
 772
 773       Dint=Dsoftshval(real(x)*smaxDf,x,y,0.,0)
 774       do i=0,nptg
 775         b=-bmax+2.*real(i)/real(nptg)*bmax
 776         write(ifhi,*) b,Dsoftshval(real(x)*smaxDf,x,y,b,0)/Dint
 777       enddo
 778
 779       write(ifhi,'(a)')    '  endarray'
 780       write(ifhi,'(a)')    'closehisto plot 0-'
 781
 782 c**********************************************************************
 783
 784       if (engy.ge.10.) then
 785        if (engy.ge.100.) then
 786         if (engy.ge.1000.) then
 787          if (engy.ge.10000.) then
 788       write(ifhi,'(a,I5)')       'openhisto name DParambn-',int(engy)
 789          else
 790       write(ifhi,'(a,I4)')       'openhisto name DParambn-',int(engy)
 791         endif
 792         else
 793       write(ifhi,'(a,I3)')       'openhisto name DParambn-',int(engy)
 794        endif
 795        else
 796       write(ifhi,'(a,I2)')       'openhisto name DParambn-',int(engy)
 797       endif
 798       else
 799       write(ifhi,'(a,I1)')       'openhisto name DParambn-',int(engy)
 800       endif
 801       write(ifhi,'(a)')       'htyp lrd'
 802       write(ifhi,'(a)')       'xmod lin ymod lin'
 803       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 804 c      write(ifhi,'(a)')       'yrange -.01 auto'
 805       write(ifhi,'(a)')       'yrange auto auto'
 806       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 807       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 808       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 809       if (xpar8.eq.1.) then
 810       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 811       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 812       endif
 813       write(ifhi,'(a)')       'array 2'
 814
 815       imax=idxD1
 816       if(iomega.eq.2)imax=1
 817
 818
 819       Dint=xDfit(0,imax,engy**2,
 820      &     real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),0.)
 821
 822       do i=0,nptg
 823         b=-bmax+2.*real(i)/real(nptg)*bmax
 824         z=xDfit(0,imax,engy**2,
 825      &       real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),b)
 826
 827         write(ifhi,*) b,z/Dint
 828       enddo
 829
 830       write(ifhi,'(a)')    '  endarray'
 831       write(ifhi,'(a)')    'closehisto plot 0-'
 832
 833
 834 c**********************************************************************
 835
 836       if (engy.ge.10.) then
 837        if (engy.ge.100.) then
 838         if (engy.ge.1000.) then
 839          if (engy.ge.10000.) then
 840       write(ifhi,'(a,I5)')  'openhisto name DEfitb-',int(engy)
 841          else
 842       write(ifhi,'(a,I4)')  'openhisto name DEfitb-',int(engy)
 843         endif
 844         else
 845       write(ifhi,'(a,I3)')  'openhisto name DEfitb-',int(engy)
 846        endif
 847        else
 848       write(ifhi,'(a,I2)')  'openhisto name DEfitb-',int(engy)
 849       endif
 850       else
 851       write(ifhi,'(a,I1)')  'openhisto name DEfitb-',int(engy)
 852       endif
 853       write(ifhi,'(a)')       'htyp pfc'
 854       write(ifhi,'(a)')       'xmod lin ymod lin'
 855       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 856 c      write(ifhi,'(a)')       'yrange -.01 auto'
 857       write(ifhi,'(a)')       'yrange auto auto'
 858       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 859       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 860       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 861       if (xpar8.eq.1.) then
 862       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 863       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 864       endif
 865       write(ifhi,'(a)')       'array 2'
 866
 867       sig2=sigma2(x,2)
 868       if(sig2.le.0.) sig2=1.e+10
 869       do i=0,nptg
 870         b=-bmax+2.*real(i)/real(nptg)*bmax
 871         write(ifhi,*) b,exp(-b**2/sig2)
 872       enddo
 873
 874       write(ifhi,'(a)')    '  endarray'
 875       write(ifhi,'(a)')    'closehisto plot 0-'
 876
 877
 878 c**********************************************************************
 879
 880       if (engy.ge.10.) then
 881        if (engy.ge.100.) then
 882         if (engy.ge.1000.) then
 883          if (engy.ge.10000.) then
 884       write(ifhi,'(a,I5)')  'openhisto name DEfitbnSoft-',int(engy)
 885          else
 886       write(ifhi,'(a,I4)')  'openhisto name DEfitbnSoft-',int(engy)
 887         endif
 888         else
 889       write(ifhi,'(a,I3)')  'openhisto name DEfitbnSoft-',int(engy)
 890        endif
 891        else
 892       write(ifhi,'(a,I2)')  'openhisto name DEfitbnSoft-',int(engy)
 893       endif
 894       else
 895       write(ifhi,'(a,I1)')  'openhisto name DEfitbnSoft-',int(engy)
 896       endif
 897       write(ifhi,'(a)')       'htyp pft'
 898       write(ifhi,'(a)')       'xmod lin ymod lin'
 899       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 900 c      write(ifhi,'(a)')       'yrange -.01 auto'
 901       write(ifhi,'(a)')       'yrange auto auto'
 902       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 903       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 904       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 905       if (xpar8.eq.1.) then
 906       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 907       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 908       endif
 909       write(ifhi,'(a)')       'array 2'
 910
 911       sig2=sigma2(x,0)
 912       if(sig2.le.0.) sig2=1.e+10
 913       do i=0,nptg
 914         b=-bmax+2.*real(i)/real(nptg)*bmax
 915         write(ifhi,*) b,exp(-b**2/sig2)
 916       enddo
 917
 918       write(ifhi,'(a)')    '  endarray'
 919       write(ifhi,'(a)')    'closehisto plot 0-'
 920
 921
 922 c**********************************************************************
 923
 924
 925       if (engy.ge.10.) then
 926        if (engy.ge.100.) then
 927         if (engy.ge.1000.) then
 928          if (engy.ge.10000.) then
 929       write(ifhi,'(a,I5)')  'openhisto name DEfitbnSh-',int(engy)
 930          else
 931       write(ifhi,'(a,I4)')  'openhisto name DEfitbnSh-',int(engy)
 932         endif
 933         else
 934       write(ifhi,'(a,I3)')  'openhisto name DEfitbnSh-',int(engy)
 935        endif
 936        else
 937       write(ifhi,'(a,I2)')  'openhisto name DEfitbnSh-',int(engy)
 938       endif
 939       else
 940       write(ifhi,'(a,I1)')  'openhisto name DEfitbnSh-',int(engy)
 941       endif
 942       write(ifhi,'(a)')       'htyp poc'
 943       write(ifhi,'(a)')       'xmod lin ymod lin'
 944       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
 945 c      write(ifhi,'(a)')       'yrange -.01 auto'
 946       write(ifhi,'(a)')       'yrange auto auto'
 947       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
 948       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
 949       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
 950       if (xpar8.eq.1.) then
 951       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
 952       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
 953       endif
 954       write(ifhi,'(a)')       'array 2'
 955
 956       sig2=sigma2(x,1)
 957       if(sig2.le.0.) sig2=1.e+10
 958       do i=0,nptg
 959         b=-bmax+2.*real(i)/real(nptg)*bmax
 960         write(ifhi,*) b,exp(-b**2/sig2)
 961       enddo
 962
 963       write(ifhi,'(a)')    '  endarray'
 964
 965       write(ifhi,'(a)')    'closehisto plot 0'
 966
 967       end
 968
 969 c----------------------------------------------------------------------
 970       subroutine xbnParD
 971 c----------------------------------------------------------------------
 972
 973       include 'epos.inc'
 974       include 'epos.incpar'
 975
 976       double precision x,y,Dsoftshval,z,xDfit
 977
 978
 979
 980       nptg=50     !number of point for the graphs
 981       bmax=xpar2                   !value max of b (impact parameter)
 982       bmax=max(0.1,bmax)
 983       y=dble(xpar3)              !value of y (rapidity)
 984       x=dble(xpar4)
 985
 986       if (engy.ge.10.) then
 987        if (engy.ge.100.) then
 988         if (engy.ge.1000.) then
 989          if (engy.ge.10000.) then
 990       write(ifhi,'(a,I5)')  'openhisto name DExactbn-',int(engy)
 991          else
 992       write(ifhi,'(a,I4)')  'openhisto name DExactbn-',int(engy)
 993         endif
 994         else
 995       write(ifhi,'(a,I3)')  'openhisto name DExactbn-',int(engy)
 996        endif
 997        else
 998       write(ifhi,'(a,I2)')  'openhisto name DExactbn-',int(engy)
 999       endif
1000       else
1001       write(ifhi,'(a,I1)')  'openhisto name DExactbn-',int(engy)
1002       endif
1003       write(ifhi,'(a)')       'htyp lbu'
1004       write(ifhi,'(a)')       'xmod lin ymod lin'
1005       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1006 c      write(ifhi,'(a)')       'yrange -.01 auto'
1007       write(ifhi,'(a)')       'yrange auto auto'
1008       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1009       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1010       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1011       if (xpar8.eq.1.) then
1012       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1013       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1014       endif
1015       write(ifhi,'(a)')       'array 2'
1016
1017       Dint=Dsoftshval(real(x)*smaxDf,x,y,0.,0)
1018       do i=0,nptg
1019         b=-bmax+2.*real(i)/real(nptg)*bmax
1020         write(ifhi,*) b,Dsoftshval(real(x)*smaxDf,x,y,b,0)/Dint
1021       enddo
1022
1023       write(ifhi,'(a)')    '  endarray'
1024       write(ifhi,'(a)')    'closehisto plot 0-'
1025
1026 c**********************************************************************
1027
1028       if (engy.ge.10.) then
1029        if (engy.ge.100.) then
1030         if (engy.ge.1000.) then
1031          if (engy.ge.10000.) then
1032       write(ifhi,'(a,I5)')       'openhisto name DParambn-',int(engy)
1033          else
1034       write(ifhi,'(a,I4)')       'openhisto name DParambn-',int(engy)
1035         endif
1036         else
1037       write(ifhi,'(a,I3)')       'openhisto name DParambn-',int(engy)
1038        endif
1039        else
1040       write(ifhi,'(a,I2)')       'openhisto name DParambn-',int(engy)
1041       endif
1042       else
1043       write(ifhi,'(a,I1)')       'openhisto name DParambn-',int(engy)
1044       endif
1045       write(ifhi,'(a)')       'htyp lrd'
1046       write(ifhi,'(a)')       'xmod lin ymod lin'
1047       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1048 c      write(ifhi,'(a)')       'yrange -.01 auto'
1049       write(ifhi,'(a)')       'yrange auto auto'
1050       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1051       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1052       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1053       if (xpar8.eq.1.) then
1054       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1055       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1056       endif
1057       write(ifhi,'(a)')       'array 2'
1058
1059       Dint=xDfit(0,1,smaxDf,
1060      &     real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),0.)
1061
1062       imax=idxD1
1063       if(iomega.eq.2)imax=1
1064
1065       do i=0,nptg
1066         b=-bmax+2.*real(i)/real(nptg)*bmax
1067         z=xDfit(0,imax,smaxDf,
1068      &       real(dsqrt(x)*dexp(y)),real(dsqrt(x)*dexp(-y)),b)
1069
1070         write(ifhi,*) b,z/Dint
1071       enddo
1072
1073       write(ifhi,'(a)')    '  endarray'
1074       write(ifhi,'(a)')    'closehisto plot 0-'
1075
1076
1077 c**********************************************************************
1078
1079       if (engy.ge.10.) then
1080        if (engy.ge.100.) then
1081         if (engy.ge.1000.) then
1082          if (engy.ge.10000.) then
1083       write(ifhi,'(a,I5)')  'openhisto name DEfitb-',int(engy)
1084          else
1085       write(ifhi,'(a,I4)')  'openhisto name DEfitb-',int(engy)
1086         endif
1087         else
1088       write(ifhi,'(a,I3)')  'openhisto name DEfitb-',int(engy)
1089        endif
1090        else
1091       write(ifhi,'(a,I2)')  'openhisto name DEfitb-',int(engy)
1092       endif
1093       else
1094       write(ifhi,'(a,I1)')  'openhisto name DEfitb-',int(engy)
1095       endif
1096       write(ifhi,'(a)')       'htyp pfc'
1097       write(ifhi,'(a)')       'xmod lin ymod lin'
1098       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1099 c      write(ifhi,'(a)')       'yrange -.01 auto'
1100       write(ifhi,'(a)')       'yrange auto auto'
1101       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1102       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1103       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1104       if (xpar8.eq.1.) then
1105       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1106       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1107       endif
1108       write(ifhi,'(a)')       'array 2'
1109
1110       sig2=sigma2(x,2)
1111       if(sig2.le.0.) sig2=1.e+10
1112       do i=0,nptg
1113         b=-bmax+2.*real(i)/real(nptg)*bmax
1114         write(ifhi,*) b,exp(-b**2/sig2)
1115       enddo
1116
1117       write(ifhi,'(a)')    '  endarray'
1118       write(ifhi,'(a)')    'closehisto plot 0-'
1119
1120
1121 c**********************************************************************
1122
1123       if (engy.ge.10.) then
1124        if (engy.ge.100.) then
1125         if (engy.ge.1000.) then
1126          if (engy.ge.10000.) then
1127       write(ifhi,'(a,I5)')  'openhisto name DEintb-',int(engy)
1128          else
1129       write(ifhi,'(a,I4)')  'openhisto name DEintb-',int(engy)
1130         endif
1131         else
1132       write(ifhi,'(a,I3)')  'openhisto name DEintb-',int(engy)
1133        endif
1134        else
1135       write(ifhi,'(a,I2)')  'openhisto name DEintb-',int(engy)
1136       endif
1137       else
1138       write(ifhi,'(a,I1)')  'openhisto name DEintb-',int(engy)
1139       endif
1140       write(ifhi,'(a)')       'htyp pft'
1141       write(ifhi,'(a)')       'xmod lin ymod lin'
1142       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1143 c      write(ifhi,'(a)')       'yrange -.01 auto'
1144       write(ifhi,'(a)')       'yrange auto auto'
1145       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1146       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1147       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1148       if (xpar8.eq.1.) then
1149       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1150       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1151       endif
1152       write(ifhi,'(a)')       'array 2'
1153
1154       sig2=sigma1i(x)
1155       do i=0,nptg
1156         b=-bmax+2.*real(i)/real(nptg)*bmax
1157         write(ifhi,*) b,exp(-b**2/sig2)
1158       enddo
1159
1160       write(ifhi,'(a)')    '  endarray'
1161       write(ifhi,'(a)')    'closehisto plot 0-'
1162
1163
1164 c**********************************************************************
1165
1166
1167       if (engy.ge.10.) then
1168        if (engy.ge.100.) then
1169         if (engy.ge.1000.) then
1170          if (engy.ge.10000.) then
1171       write(ifhi,'(a,I5)')  'openhisto name DPfitb-',int(engy)
1172          else
1173       write(ifhi,'(a,I4)')  'openhisto name DPfitb-',int(engy)
1174         endif
1175         else
1176       write(ifhi,'(a,I3)')  'openhisto name DPfitb-',int(engy)
1177        endif
1178        else
1179       write(ifhi,'(a,I2)')  'openhisto name DPfitb-',int(engy)
1180       endif
1181       else
1182       write(ifhi,'(a,I1)')  'openhisto name DPfitb-',int(engy)
1183       endif
1184       write(ifhi,'(a)')       'htyp poc'
1185       write(ifhi,'(a)')       'xmod lin ymod lin'
1186       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1187 c      write(ifhi,'(a)')       'yrange -.01 auto'
1188       write(ifhi,'(a)')       'yrange auto auto'
1189       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1190       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1191       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1192       if (xpar8.eq.1.) then
1193       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1194       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1195       endif
1196       write(ifhi,'(a)')       'array 2'
1197
1198       sig2=xsigmafit(x)
1199       do i=0,nptg
1200         b=-bmax+2.*real(i)/real(nptg)*bmax
1201         write(ifhi,*) b,exp(-b**2/sig2)
1202       enddo
1203
1204       write(ifhi,'(a)')    '  endarray'
1205       write(ifhi,'(a)')    'closehisto plot 0'
1206
1207
1208
1209       end
1210
1211 c----------------------------------------------------------------------
1212       subroutine xbParD
1213 c----------------------------------------------------------------------
1214
1215       include 'epos.inc'
1216       include 'epos.incsem'
1217       include 'epos.incpar'
1218
1219       double precision x,om5s,xDfit,z(0:200),y
1220 c     & ,omYuncut,omNpuncut,omYcut
1221
1222       nptg=50                  !number of point for the graphs
1223       x=dble(xpar4)                 !value of biniDf (impact parameter)
1224       y=dble(xpar3)                 !value of y (rapidity)
1225       if(x.gt.dexp(-2.d0*y))x=x*dexp(-2.d0*y)
1226
1227       bmax=xpar2
1228       bmax=max(0.1,bmax)
1229       t=1.
1230 c      iqqN=0
1231 c      iqq=int(xpar7)
1232       iqq1=-1
1233       iqq2=-1
1234
1235       if (engy.ge.10.) then
1236        if (engy.ge.100.) then
1237         if (engy.ge.1000.) then
1238          if (engy.ge.10000.) then
1239       write(ifhi,'(a,I5)')       'openhisto name DbParam-',int(engy)
1240          else
1241       write(ifhi,'(a,I4)')       'openhisto name DbParam-',int(engy)
1242         endif
1243         else
1244       write(ifhi,'(a,I3)')       'openhisto name DbParam-',int(engy)
1245        endif
1246        else
1247       write(ifhi,'(a,I2)')       'openhisto name DbParam-',int(engy)
1248       endif
1249       else
1250       write(ifhi,'(a,I1)')       'openhisto name DbParam-',int(engy)
1251       endif
1252       write(ifhi,'(a)')       'htyp lbu'
1253       if(xpar5.eq.0.)then
1254         write(ifhi,'(a)')       'xmod lin ymod lin'
1255       else
1256         write(ifhi,'(a)')       'xmod lin ymod log'
1257       endif
1258       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1259       write(ifhi,'(a)')       'yrange auto auto'
1260       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1261       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1262       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1263       if (xpar8.eq.1.) then
1264       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1265       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1266       endif
1267       write(ifhi,'(a)')       'array 2'
1268
1269       xp=sqrt(real(x))*exp(real(y))
1270       xm=sqrt(real(x))*exp(-real(y))
1271       imax=idxD1
1272       if(iomega.eq.2)imax=1
1273
1274       if(xpar6.eq.1.)t=real(xDfit(0,imax,smaxDf,xp,xm,0.))
1275       if(abs(t).lt.1.d-8)t=1.
1276       do i=0,nptg
1277         b=-bmax+2.*real(i)/real(nptg)*bmax
1278          z(i)=xDfit(0,imax,smaxDf,xp,xm,b)/dble(t)
1279          write(ifhi,*) b,z(i)
1280        enddo
1281
1282       write(ifhi,'(a)')    '  endarray'
1283
1284       write(ifhi,'(a)')    'closehisto plot 0-'
1285
1286
1287 c**********************************************************************
1288
1289       if (engy.ge.10.) then
1290         if (engy.ge.100.) then
1291           if (engy.ge.1000.) then
1292             if (engy.ge.10000.) then
1293            write(ifhi,'(a,I5)')  'openhisto name DbParamI-',int(engy)
1294             else
1295            write(ifhi,'(a,I4)')  'openhisto name DbParamI-',int(engy)
1296             endif
1297           else
1298            write(ifhi,'(a,I3)')  'openhisto name DbParamI-',int(engy)
1299           endif
1300         else
1301           write(ifhi,'(a,I2)')  'openhisto name DbParamI-',int(engy)
1302         endif
1303       else
1304         write(ifhi,'(a,I1)')  'openhisto name DbParamI-',int(engy)
1305       endif
1306       write(ifhi,'(a)')       'htyp lga'
1307       write(ifhi,'(a)')       'xmod lin ymod lin'
1308       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1309       write(ifhi,'(a)')       'yrange .01 auto'
1310       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1311       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1312       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1313       if (xpar8.eq.1.) then
1314         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1315         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1316       endif
1317       write(ifhi,'(a)')       'array 2'
1318
1319       if(xpar6.eq.1.)t=(alpDs(idxD0,iclegy,iclpro,icltar)
1320      &   *(smaxDf*real(x))**betDs(idxD0,iclegy,iclpro,icltar)
1321      &                +alpDs(1,iclegy,iclpro,icltar)
1322      &   *(smaxDf*real(x))**betDs(1,iclegy,iclpro,icltar))
1323      &   /real(idxD0+1)
1324       if(abs(t).lt.1.d-8)t=1.
1325       do i=0,nptg
1326         b=-bmax+2.*real(i)/real(nptg)*bmax
1327         write(ifhi,*) b,(alpDs(idxD0,iclegy,iclpro,icltar)/t
1328      &   *(smaxDf*real(x))**(betDs(idxD0,iclegy,iclpro,icltar)
1329      &           +gamDs(idxD0,iclegy,iclpro,icltar)*b**2.)
1330      &        *exp(-b**2./delDs(idxD0,iclegy,iclpro,icltar))
1331      &                 +alpDs(1,iclegy,iclpro,icltar)/t
1332      &   *(smaxDf*real(x))**(betDs(1,iclegy,iclpro,icltar)
1333      &           +gamDs(1,iclegy,iclpro,icltar)*b**2.)
1334      &        *exp(-b**2./delDs(1,iclegy,iclpro,icltar)))
1335      &   /real(idxD0+1)
1336       enddo
1337
1338       write(ifhi,'(a)')    '  endarray'
1339       write(ifhi,'(a)')    'closehisto plot 0-'
1340
1341 c**********************************************************************
1342
1343       if (engy.ge.10.) then
1344         if (engy.ge.100.) then
1345           if (engy.ge.1000.) then
1346             if (engy.ge.10000.) then
1347            write(ifhi,'(a,I5)')  'openhisto name DbExaI-',int(engy)
1348             else
1349            write(ifhi,'(a,I4)')  'openhisto name DbExaI-',int(engy)
1350             endif
1351           else
1352            write(ifhi,'(a,I3)')  'openhisto name DbExaI-',int(engy)
1353           endif
1354         else
1355           write(ifhi,'(a,I2)')  'openhisto name DbExaI-',int(engy)
1356         endif
1357       else
1358         write(ifhi,'(a,I1)')  'openhisto name DbExaI-',int(engy)
1359       endif
1360       write(ifhi,'(a)')       'htyp poc'
1361       write(ifhi,'(a)')       'xmod lin ymod lin'
1362       write(ifhi,'(a,2e11.3)')'xrange',-bmax,bmax
1363       write(ifhi,'(a)')       'yrange .01 auto'
1364       write(ifhi,'(a)')      'text 0 0 "xaxis b"'
1365       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1366       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1367       if (xpar8.eq.1.) then
1368         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "x=',x,'"'
1369         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "y=',y,'"'
1370       endif
1371       write(ifhi,'(a)')       'array 2'
1372
1373       if(xpar6.eq.1.)t=2*real(om5s(real(x)*smaxDf,x,0.d0,0.,iqq1,iqq2)
1374      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar))))
1375       if(abs(t).lt.1.d-8)t=1.
1376       nptg=nptg/2
1377       do i=0,nptg
1378         b=-bmax+2.*real(i)/real(nptg)*bmax
1379         write(ifhi,*) b,2*real(om5s(real(x)*smaxDf,x,y,b,iqq1,iqq2)
1380      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar))))/t
1381       enddo
1382
1383       write(ifhi,'(a)')    '  endarray'
1384
1385       write(ifhi,'(a)')    'closehisto plot 0'
1386
1387       end
1388
1389 c----------------------------------------------------------------------
1390       subroutine xGexaJ
1391 c----------------------------------------------------------------------
1392
1393       include 'epos.inc'
1394       include 'epos.incsem'
1395       include 'epos.incpar'
1396
1397       parameter(nptg=50) !number of point for the graphs
1398       double precision x,omJ1(0:nptg),xtmp,w(0:nptg)
1399       double precision z(0:nptg),om51!,omYuncut,omYuncutJ
1400       double precision xminr,y,t,omJ2(0:nptg),omJ3(0:nptg),omJ4(0:nptg)
1401      &,omJ5(0:nptg)!,omJ6(0:nptg),omJ7(0:nptg)
1402
1403
1404       biniDf=xpar2                 !value of biniDf (impact parameter)
1405       y=dble(xpar3)
1406       t=1.d0
1407       xtmp=xmaxDf
1408       xmaxDf=dexp(-2.d0*y)
1409
1410       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
1411
1412       nnnmax=1
1413       if(iscreen.ne.0)nnnmax=2
1414
1415       do nnn=1,nnnmax
1416
1417         if(nnn.eq.2)then
1418           iscreensave=iscreen
1419           iscreen=0
1420         endif
1421
1422       if (engy.ge.10.) then
1423        if (engy.ge.100.) then
1424         if (engy.ge.1000.) then
1425          if (engy.ge.10000.) then
1426       write(ifhi,'(a,I5)')  'openhisto name GIexa-',int(engy)
1427          else
1428       write(ifhi,'(a,I4)')  'openhisto name GIexa-',int(engy)
1429         endif
1430         else
1431       write(ifhi,'(a,I3)')  'openhisto name GIexa-',int(engy)
1432        endif
1433        else
1434       write(ifhi,'(a,I2)')  'openhisto name GIexa-',int(engy)
1435       endif
1436       else
1437       write(ifhi,'(a,I1)')  'openhisto name GIexa-',int(engy)
1438       endif
1439       write(ifhi,'(a)')       'htyp lru'
1440       if(xpar5.eq.0.)then
1441         write(ifhi,'(a)')       'xmod log ymod log'
1442       else
1443         write(ifhi,'(a)')       'xmod log ymod lin'
1444       endif
1445       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1446       write(ifhi,'(a)')       'yrange auto auto'
1447       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1448       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1449       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1450       if (xpar8.eq.1.) then
1451       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1452       endif
1453       write(ifhi,'(a)')       'array 2'
1454
1455       do i=0,nptg
1456         x=xminr
1457         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1458         z(i)=0.d0
1459         do j=0,5
1460           z(i)=z(i)+om51(x,y,biniDf,j,j)
1461         enddo
1462         w(i)=om51(x,y,biniDf,-1,-1)
1463         if(xpar6.eq.0)t=z(i)/w(i)
1464         if(xpar6.eq.1)t=z(i)
1465         write(ifhi,*) x,z(i)/t
1466       enddo
1467
1468       write(ifhi,'(a)')    '  endarray'
1469       write(ifhi,'(a)')    'closehisto plot 0-'
1470
1471
1472 c*************************************************************************
1473
1474       if (engy.ge.10.) then
1475        if (engy.ge.100.) then
1476         if (engy.ge.1000.) then
1477          if (engy.ge.10000.) then
1478       write(ifhi,'(a,I5)')       'openhisto name GIsoft-',int(engy)
1479          else
1480       write(ifhi,'(a,I4)')       'openhisto name GIsoft-',int(engy)
1481         endif
1482         else
1483       write(ifhi,'(a,I3)')       'openhisto name GIsoft-',int(engy)
1484        endif
1485        else
1486       write(ifhi,'(a,I2)')       'openhisto name GIsoft-',int(engy)
1487       endif
1488       else
1489       write(ifhi,'(a,I1)')       'openhisto name GIsoft-',int(engy)
1490       endif
1491       write(ifhi,'(a)')       'htyp lba'
1492       write(ifhi,'(a)')       'xmod log ymod log'
1493       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1494       write(ifhi,'(a)')       'yrange auto auto'
1495       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1496       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1497       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1498       if (xpar8.eq.1.) then
1499       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1500       endif
1501       write(ifhi,'(a)')       'array 2'
1502
1503       do i=0,nptg
1504         x=xminr
1505         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1506         omJ1(i)=om51(x,y,biniDf,0,0)
1507         if(xpar6.eq.0)t=z(i)/w(i)
1508         if(xpar6.eq.1)t=z(i)
1509         write(ifhi,*) x,omJ1(i)/t
1510       enddo
1511
1512       write(ifhi,'(a)')    '  endarray'
1513       write(ifhi,'(a)')    'closehisto plot 0-'
1514
1515 c*************************************************************************
1516
1517       if (engy.ge.10.) then
1518        if (engy.ge.100.) then
1519         if (engy.ge.1000.) then
1520          if (engy.ge.10000.) then
1521       write(ifhi,'(a,I5)')       'openhisto name GIgg-',int(engy)
1522          else
1523       write(ifhi,'(a,I4)')       'openhisto name GIgg-',int(engy)
1524         endif
1525         else
1526       write(ifhi,'(a,I3)')       'openhisto name GIgg-',int(engy)
1527        endif
1528        else
1529       write(ifhi,'(a,I2)')       'openhisto name GIgg-',int(engy)
1530       endif
1531       else
1532       write(ifhi,'(a,I1)')       'openhisto name GIgg-',int(engy)
1533       endif
1534       write(ifhi,'(a)')       'htyp lra'
1535       write(ifhi,'(a)')       'xmod log ymod log'
1536       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1537       write(ifhi,'(a)')       'yrange auto auto'
1538       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1539       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1540       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1541       if (xpar8.eq.1.) then
1542       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1543       endif
1544       write(ifhi,'(a)')       'array 2'
1545
1546       do i=0,nptg
1547         x=xminr
1548         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1549         omJ2(i)=om51(x,y,biniDf,1,1)
1550         if(xpar6.eq.0)t=z(i)/w(i)
1551         if(xpar6.eq.1)t=z(i)
1552         write(ifhi,*) x,omJ2(i)/t
1553       enddo
1554
1555       write(ifhi,'(a)')    '  endarray'
1556       write(ifhi,'(a)')    'closehisto plot 0-'
1557
1558 c*************************************************************************
1559
1560       if (engy.ge.10.) then
1561        if (engy.ge.100.) then
1562         if (engy.ge.1000.) then
1563          if (engy.ge.10000.) then
1564       write(ifhi,'(a,I5)')       'openhisto name GIgq-',int(engy)
1565          else
1566       write(ifhi,'(a,I4)')       'openhisto name GIgq-',int(engy)
1567         endif
1568         else
1569       write(ifhi,'(a,I3)')       'openhisto name GIgq-',int(engy)
1570        endif
1571        else
1572       write(ifhi,'(a,I2)')       'openhisto name GIgq-',int(engy)
1573       endif
1574       else
1575       write(ifhi,'(a,I1)')       'openhisto name GIgq-',int(engy)
1576       endif
1577       write(ifhi,'(a)')       'htyp lga'
1578       write(ifhi,'(a)')       'xmod log ymod log'
1579       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1580       write(ifhi,'(a)')       'yrange auto auto'
1581       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1582       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1583       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1584       if (xpar8.eq.1.) then
1585       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1586       endif
1587       write(ifhi,'(a)')       'array 2'
1588
1589       do i=0,nptg
1590         x=xminr
1591         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1592         omJ3(i)=om51(x,y,biniDf,2,3)
1593         if(xpar6.eq.0)t=z(i)/w(i)
1594         if(xpar6.eq.1)t=z(i)
1595         write(ifhi,*) x,omJ3(i)/t
1596       enddo
1597
1598       write(ifhi,'(a)')    '  endarray'
1599       write(ifhi,'(a)')    'closehisto plot 0-'
1600
1601 c*************************************************************************
1602
1603       if (engy.ge.10.) then
1604        if (engy.ge.100.) then
1605         if (engy.ge.1000.) then
1606          if (engy.ge.10000.) then
1607       write(ifhi,'(a,I5)')       'openhisto name GIqq-',int(engy)
1608          else
1609       write(ifhi,'(a,I4)')       'openhisto name GIqq-',int(engy)
1610         endif
1611         else
1612       write(ifhi,'(a,I3)')       'openhisto name GIqq-',int(engy)
1613        endif
1614        else
1615       write(ifhi,'(a,I2)')       'openhisto name GIqq-',int(engy)
1616       endif
1617       else
1618       write(ifhi,'(a,I1)')       'openhisto name GIqq-',int(engy)
1619       endif
1620       write(ifhi,'(a)')       'htyp lya'
1621       write(ifhi,'(a)')       'xmod log ymod log'
1622       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1623       write(ifhi,'(a)')       'yrange auto auto'
1624       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1625       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1626       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1627       if (xpar8.eq.1.) then
1628       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1629       endif
1630       write(ifhi,'(a)')       'array 2'
1631
1632       do i=0,nptg
1633         x=xminr
1634         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1635         omJ4(i)=om51(x,y,biniDf,4,4)
1636         if(xpar6.eq.0)t=z(i)/w(i)
1637         if(xpar6.eq.1)t=z(i)
1638         write(ifhi,*) x,omJ4(i)/t
1639       enddo
1640
1641       write(ifhi,'(a)')    '  endarray'
1642       write(ifhi,'(a)')    'closehisto plot 0-'
1643
1644 c*************************************************************************
1645
1646       if (engy.ge.10.) then
1647        if (engy.ge.100.) then
1648         if (engy.ge.1000.) then
1649          if (engy.ge.10000.) then
1650       write(ifhi,'(a,I5)')       'openhisto name GIdif-',int(engy)
1651          else
1652       write(ifhi,'(a,I4)')       'openhisto name GIdif-',int(engy)
1653         endif
1654         else
1655       write(ifhi,'(a,I3)')       'openhisto name GIdif-',int(engy)
1656        endif
1657        else
1658       write(ifhi,'(a,I2)')       'openhisto name GIdif-',int(engy)
1659       endif
1660       else
1661       write(ifhi,'(a,I1)')       'openhisto name GIdif-',int(engy)
1662       endif
1663       write(ifhi,'(a)')       'htyp lfa'
1664       write(ifhi,'(a)')       'xmod log ymod log'
1665       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1666       write(ifhi,'(a)')       'yrange auto auto'
1667       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1668       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1669       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1670       if (xpar8.eq.1.) then
1671       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1672       endif
1673       write(ifhi,'(a)')       'array 2'
1674
1675       do i=0,nptg
1676         x=xminr
1677         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1678         omJ5(i)=om51(x,y,biniDf,5,5)
1679         if(xpar6.eq.0)t=z(i)/w(i)
1680         if(xpar6.eq.1)t=z(i)
1681         write(ifhi,*) x,omJ5(i)/t
1682       enddo
1683
1684       write(ifhi,'(a)')    '  endarray'
1685       write(ifhi,'(a)')    'closehisto plot 0-'
1686
1687 c*************************************************************************
1688
1689       if (engy.ge.10.) then
1690        if (engy.ge.100.) then
1691         if (engy.ge.1000.) then
1692          if (engy.ge.10000.) then
1693       write(ifhi,'(a,I5)')       'openhisto name GItot-',int(engy)
1694          else
1695       write(ifhi,'(a,I4)')       'openhisto name GItot-',int(engy)
1696         endif
1697         else
1698       write(ifhi,'(a,I3)')       'openhisto name GItot-',int(engy)
1699        endif
1700        else
1701       write(ifhi,'(a,I2)')       'openhisto name GItot-',int(engy)
1702       endif
1703       else
1704       write(ifhi,'(a,I1)')       'openhisto name GItot-',int(engy)
1705       endif
1706       write(ifhi,'(a)')       'htyp pfc'
1707       write(ifhi,'(a)')       'xmod log ymod log'
1708       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1709       write(ifhi,'(a)')       'yrange auto auto'
1710       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1711       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,s,b)" '
1712       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1713       if (xpar8.eq.1.) then
1714       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1715       endif
1716       write(ifhi,'(a)')       'array 2'
1717
1718       do i=0,nptg
1719         x=xminr
1720         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1721         if(xpar6.eq.0)t=z(i)/w(i)
1722         if(xpar6.eq.1)t=z(i)
1723         write(ifhi,*) x,(omJ1(i)+omJ2(i)+omJ3(i)+omJ4(i)+omJ5(i))
1724      &                  /t
1725       enddo
1726
1727       write(ifhi,'(a)')    '  endarray'
1728
1729       write(ifhi,'(a)')    'closehisto plot 0'
1730
1731         if(nnn.eq.2)iscreen=iscreensave
1732
1733       enddo
1734
1735       xmaxDf=xtmp
1736
1737       end
1738
1739
1740
1741 c----------------------------------------------------------------------
1742       subroutine xsParD
1743 c----------------------------------------------------------------------
1744
1745       include 'epos.inc'
1746       include 'epos.incsem'
1747       include 'epos.incpar'
1748
1749       double precision x,xminr,y,om5s,xtmp
1750 c     & ,t,omYuncut
1751
1752       nptg=50                  !number of point for the graphs
1753       biniDf=xpar2                 !value of biniDf (impact parameter)
1754       y=dble(xpar3)                 !value of y (rapidity)
1755       xtmp=xmaxDf
1756       xmaxDf=dexp(-2.d0*y)
1757       call Class('xsParD     ')
1758       iqq1=-1
1759       iqq2=-1
1760 c      iqqN=0
1761
1762       xminr=dble(egylow/engy)**2.d0  !value of xminr for plotting the function
1763
1764 c**********************************************************************
1765
1766       if (engy.ge.10.) then
1767         if (engy.ge.100.) then
1768           if (engy.ge.1000.) then
1769             if (engy.ge.10000.) then
1770            write(ifhi,'(a,I5)')  'openhisto name DParamI-',int(engy)
1771             else
1772            write(ifhi,'(a,I4)')  'openhisto name DParamI-',int(engy)
1773             endif
1774           else
1775            write(ifhi,'(a,I3)')  'openhisto name DParamI-',int(engy)
1776           endif
1777         else
1778           write(ifhi,'(a,I2)')  'openhisto name DParamI-',int(engy)
1779         endif
1780       else
1781         write(ifhi,'(a,I1)')  'openhisto name DParamI-',int(engy)
1782       endif
1783       write(ifhi,'(a)')       'htyp lga'
1784       if(xpar5.eq.0.)then
1785         write(ifhi,'(a)')       'xmod log ymod log'
1786       else
1787         write(ifhi,'(a)')       'xmod log ymod lin'
1788       endif
1789       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1790       write(ifhi,'(a)')       'yrange auto auto'
1791       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1792       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1793       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1794       if (xpar8.eq.1.) then
1795         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1796       endif
1797       write(ifhi,'(a)')       'array 2'
1798
1799       do i=0,nptg
1800         x=xminr
1801         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1802         egy=sqrt(real(x))*engy
1803         iclegy=1+int((log(egy)-log(egylow))/log(egyfac))
1804         write(ifhi,*) x,(alpDs(idxD0,iclegy,iclpro,icltar)
1805      &   *(smaxDf*real(x))**(betDs(idxD0,iclegy,iclpro,icltar)
1806      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
1807      &     *exp(-biniDf**2./delDs(idxD0,iclegy,iclpro,icltar))
1808      &                 +alpDs(1,iclegy,iclpro,icltar)
1809      &   *(smaxDf*real(x))**(betDs(1,iclegy,iclpro,icltar)
1810      &           +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
1811      &        *exp(-biniDf**2./delDs(1,iclegy,iclpro,icltar)))
1812      &   /real(idxD0+1)
1813       enddo
1814
1815       write(ifhi,'(a)')    '  endarray'
1816       write(ifhi,'(a)')    'closehisto plot 0-'
1817
1818 c**********************************************************************
1819
1820       if (engy.ge.10.) then
1821         if (engy.ge.100.) then
1822           if (engy.ge.1000.) then
1823             if (engy.ge.10000.) then
1824            write(ifhi,'(a,I5)')  'openhisto name DExaI-',int(engy)
1825             else
1826            write(ifhi,'(a,I4)')  'openhisto name DExaI-',int(engy)
1827             endif
1828           else
1829            write(ifhi,'(a,I3)')  'openhisto name DExaI-',int(engy)
1830           endif
1831         else
1832           write(ifhi,'(a,I2)')  'openhisto name DExaI-',int(engy)
1833         endif
1834       else
1835         write(ifhi,'(a,I1)')  'openhisto name DExaI-',int(engy)
1836       endif
1837       write(ifhi,'(a)')       'htyp poc'
1838       write(ifhi,'(a)')       'xmod log ymod log'
1839       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
1840       write(ifhi,'(a)')       'yrange .01 auto'
1841       write(ifhi,'(a)')      'text 0 0 "xaxis x"'
1842       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1843       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1844       if (xpar8.eq.1.) then
1845         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1846       endif
1847       write(ifhi,'(a)')       'array 2'
1848
1849       nptg=nptg/2
1850       do i=0,nptg
1851         x=xminr
1852         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
1853         egy=sqrt(real(x))*engy
1854         iclegy=1+int((log(egy)-log(egylow))/log(egyfac))
1855       write(ifhi,*) x,2*om5s(real(x)*smaxDf,1.d0,0.d0,biniDf,iqq1,iqq2)
1856      &              /dble(chad(iclpro)*chad(icltar))
1857       enddo
1858
1859       write(ifhi,'(a)')    '  endarray'
1860       write(ifhi,'(a)')    'closehisto plot 0'
1861
1862       xmaxDf=xtmp
1863
1864       end
1865
1866 c----------------------------------------------------------------------
1867       subroutine xyParD
1868 c----------------------------------------------------------------------
1869
1870       include 'epos.inc'
1871       include 'epos.incsem'
1872       include 'epos.incpar'
1873
1874       double precision x,om5s,xDfit,z(0:200),ymax,y
1875 c     & ,omYuncut,omNpuncut
1876
1877       nptg=50                  !number of point for the graphs
1878       biniDf=xpar2                 !value of biniDf (impact parameter)
1879
1880       x=dble(xpar4)
1881       if(x.le.1.d-20)x=1.d0/dble(engy)
1882       ymax=-.5d0*dlog(x)
1883       iqq1=-1
1884       iqq2=-1
1885 c      iqqN=0
1886
1887
1888       if (engy.ge.10.) then
1889        if (engy.ge.100.) then
1890         if (engy.ge.1000.) then
1891          if (engy.ge.10000.) then
1892       write(ifhi,'(a,I5)')       'openhisto name DyParam-',int(engy)
1893          else
1894       write(ifhi,'(a,I4)')       'openhisto name DyParam-',int(engy)
1895         endif
1896         else
1897       write(ifhi,'(a,I3)')       'openhisto name DyParam-',int(engy)
1898        endif
1899        else
1900       write(ifhi,'(a,I2)')       'openhisto name DyParam-',int(engy)
1901       endif
1902       else
1903       write(ifhi,'(a,I1)')       'openhisto name DyParam-',int(engy)
1904       endif
1905       write(ifhi,'(a)')       'htyp lbu'
1906       if(xpar5.eq.0.)then
1907         write(ifhi,'(a)')       'xmod lin ymod lin'
1908       else
1909         write(ifhi,'(a)')       'xmod lin ymod log'
1910       endif
1911       write(ifhi,'(a,2e11.3)')'xrange ',-ymax-1.d0,ymax+1.d0
1912       write(ifhi,'(a)')       'yrange auto auto'
1913       write(ifhi,'(a)')      'text 0 0 "xaxis y"'
1914       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1915       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1916       if (xpar8.eq.1.) then
1917       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1918       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "x=',x,'"'
1919       endif
1920       write(ifhi,'(a)')       'array 2'
1921
1922       imax=idxD1
1923       if(iomega.eq.2)imax=1
1924
1925       do i=0,nptg
1926         y=-ymax+(ymax+ymax)*(dble(i)/dble(nptg))
1927         xp=sqrt(real(x))*exp(real(y))
1928         xm=sqrt(real(x))*exp(-real(y))
1929          z(i)=xDfit(0,imax,engy**2,xp,xm,biniDf)
1930          write(ifhi,*) y,z(i)
1931        enddo
1932
1933       write(ifhi,'(a)')    '  endarray'
1934
1935       write(ifhi,'(a)')    'closehisto plot 0-'
1936
1937
1938 c**********************************************************************
1939
1940       if (engy.ge.10.) then
1941         if (engy.ge.100.) then
1942           if (engy.ge.1000.) then
1943             if (engy.ge.10000.) then
1944            write(ifhi,'(a,I5)')  'openhisto name DyParamI-',int(engy)
1945             else
1946            write(ifhi,'(a,I4)')  'openhisto name DyParamI-',int(engy)
1947             endif
1948           else
1949            write(ifhi,'(a,I3)')  'openhisto name DyParamI-',int(engy)
1950           endif
1951         else
1952           write(ifhi,'(a,I2)')  'openhisto name DyParamI-',int(engy)
1953         endif
1954       else
1955         write(ifhi,'(a,I1)')  'openhisto name DyParamI-',int(engy)
1956       endif
1957       write(ifhi,'(a)')       'htyp lga'
1958       write(ifhi,'(a)')       'xmod lin ymod lin'
1959       write(ifhi,'(a,2e11.3)')'xrange ',-ymax-1.d0,ymax+1.d0
1960       write(ifhi,'(a)')       'yrange auto auto'
1961       write(ifhi,'(a)')      'text 0 0 "xaxis y"'
1962       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
1963       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
1964       if (xpar8.eq.1.) then
1965         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
1966         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "x=',x,'"'
1967       endif
1968       write(ifhi,'(a)')       'array 2'
1969
1970       do i=0,nptg
1971         y=-ymax+(ymax+ymax)*(dble(i)/dble(nptg))
1972         xp=sqrt(real(x))*exp(real(y))
1973         xm=sqrt(real(x))*exp(-real(y))
1974         write(ifhi,*) y,(alpDs(idxD0,iclegy,iclpro,icltar)
1975      &     *smaxDf**(betDs(idxD0,iclegy,iclpro,icltar)
1976      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
1977      &     *xp**(betDps(idxD0,iclegy,iclpro,icltar)
1978      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
1979      &     *xm**(betDpps(idxD0,iclegy,iclpro,icltar)
1980      &          +gamDs(idxD0,iclegy,iclpro,icltar)*biniDf**2.)
1981      &     *exp(-biniDf**2./delDs(idxD0,iclegy,iclpro,icltar))
1982      &                  +alpDs(1,iclegy,iclpro,icltar)
1983      &     *smaxDf**(betDs(1,iclegy,iclpro,icltar)
1984      &            +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
1985      &     *xp**(betDps(1,iclegy,iclpro,icltar)
1986      &            +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
1987      &     *xm**(betDpps(1,iclegy,iclpro,icltar)
1988      &            +gamDs(1,iclegy,iclpro,icltar)*biniDf**2.)
1989      &        *exp(-biniDf**2./delDs(1,iclegy,iclpro,icltar)))
1990      &     /real(idxD0+1)
1991       enddo
1992
1993       write(ifhi,'(a)')    '  endarray'
1994       write(ifhi,'(a)')    'closehisto plot 0-'
1995
1996 c**********************************************************************
1997
1998       if (engy.ge.10.) then
1999         if (engy.ge.100.) then
2000           if (engy.ge.1000.) then
2001             if (engy.ge.10000.) then
2002            write(ifhi,'(a,I5)')  'openhisto name DyExaI-',int(engy)
2003             else
2004            write(ifhi,'(a,I4)')  'openhisto name DyExaI-',int(engy)
2005             endif
2006           else
2007            write(ifhi,'(a,I3)')  'openhisto name DyExaI-',int(engy)
2008           endif
2009         else
2010           write(ifhi,'(a,I2)')  'openhisto name DyExaI-',int(engy)
2011         endif
2012       else
2013         write(ifhi,'(a,I1)')  'openhisto name DyExaI-',int(engy)
2014       endif
2015       write(ifhi,'(a)')       'htyp poc'
2016       write(ifhi,'(a)')       'xmod lin ymod lin'
2017       write(ifhi,'(a,2e11.3)')'xrange ',-ymax-1.d0,ymax+1.d0
2018       write(ifhi,'(a)')       'yrange auto auto'
2019       write(ifhi,'(a)')      'text 0 0 "xaxis y"'
2020       write(ifhi,'(a,a)')    'text 0 0 "yaxis D(x+,x-,s,b)" '
2021       write(ifhi,'(a,e8.2,a)')  'text 0.1 0.9 "E=',engy,' GeV"'
2022       if (xpar8.eq.1.) then
2023         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
2024         write(ifhi,'(a,f5.2,a)')  'text 0.1 0.7 "x=',x,'"'
2025       endif
2026       write(ifhi,'(a)')       'array 2'
2027
2028       nptg=nptg/2
2029       do i=0,nptg
2030         y=-ymax+(ymax+ymax)*(dble(i)/dble(nptg))
2031         write(ifhi,*) y,2*om5s(real(x)*smaxDf,x,y,biniDf,iqq1,iqq2)
2032      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
2033       enddo
2034
2035       write(ifhi,'(a)')    '  endarray'
2036       write(ifhi,'(a)')    'closehisto plot 0'
2037
2038       end
2039
2040 c----------------------------------------------------------------------
2041       subroutine xParSigma
2042 c----------------------------------------------------------------------
2043
2044       include 'epos.inc'
2045       include 'epos.incpar'
2046
2047       double precision x,w(0:200),z(0:200),xminr,t
2048
2049       nptg=20                  !number of point for the graphs
2050
2051       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
2052
2053       if (engy.ge.10.) then
2054        if (engy.ge.100.) then
2055         if (engy.ge.1000.) then
2056          if (engy.ge.10000.) then
2057       write(ifhi,'(a,I5)')    'openhisto name SigmaReel-',int(engy)
2058          else
2059       write(ifhi,'(a,I4)')    'openhisto name SigmaReel-',int(engy)
2060         endif
2061         else
2062       write(ifhi,'(a,I3)')    'openhisto name SigmaReel-',int(engy)
2063        endif
2064        else
2065       write(ifhi,'(a,I2)')    'openhisto name SigmaReel-',int(engy)
2066       endif
2067       else
2068       write(ifhi,'(a,I1)')    'openhisto name SigmaReel-',int(engy)
2069       endif
2070       write(ifhi,'(a)')       'htyp lru'
2071       write(ifhi,'(a)')       'xmod log ymod lin'
2072       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2073       write(ifhi,'(a)')'yrange auto auto'
2074       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2075       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2076       if (xpar8.eq.1.) then
2077       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2078       endif
2079       write(ifhi,'(a)')       'array 2'
2080
2081       do i=0,nptg
2082         X=xminr
2083         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2084         z(i)=dble(sigma2(X,2))
2085         if(z(i).le.0.) z(i)=0.d0
2086         write(ifhi,'(2e14.6)') X,z(i)
2087       enddo
2088
2089       write(ifhi,'(a)')    '  endarray'
2090       write(ifhi,'(a)')    'closehisto plot 0-'
2091
2092
2093 c**********************************************************************
2094
2095       if (engy.ge.10.) then
2096        if (engy.ge.100.) then
2097         if (engy.ge.1000.) then
2098          if (engy.ge.10000.) then
2099       write(ifhi,'(a,I5)')    'openhisto name SigmaParam-',int(engy)
2100          else
2101       write(ifhi,'(a,I4)')    'openhisto name SigmaParam-',int(engy)
2102         endif
2103         else
2104       write(ifhi,'(a,I3)')    'openhisto name SigmaParam-',int(engy)
2105        endif
2106        else
2107       write(ifhi,'(a,I2)')    'openhisto name SigmaParam-',int(engy)
2108       endif
2109       else
2110       write(ifhi,'(a,I1)')    'openhisto name SigmaParam-',int(engy)
2111       endif
2112       write(ifhi,'(a)')       'htyp pfc'
2113       write(ifhi,'(a)')       'xmod log ymod lin'
2114       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2115       write(ifhi,'(a)')'yrange auto auto'
2116       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2117       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!?param!(X)"'
2118       if (xpar8.eq.1.) then
2119       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2120         endif
2121       write(ifhi,'(a)')       'array 2'
2122
2123       do i=0,nptg
2124         X=xminr
2125         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2126         w(i)=dble(xsigmafit(X))
2127         write(ifhi,'(2e14.6)') X,w(i)
2128       enddo
2129
2130       write(ifhi,'(a)')    '  endarray'
2131       write(ifhi,'(a)')    'closehisto plot 0-'
2132
2133
2134 c**********************************************************************
2135
2136       if (engy.ge.10.) then
2137        if (engy.ge.100.) then
2138         if (engy.ge.1000.) then
2139          if (engy.ge.10000.) then
2140       write(ifhi,'(a,I5)')    'openhisto name SigmaInt-',int(engy)
2141          else
2142       write(ifhi,'(a,I4)')    'openhisto name SigmaInt-',int(engy)
2143         endif
2144         else
2145       write(ifhi,'(a,I3)')    'openhisto name SigmaInt-',int(engy)
2146        endif
2147        else
2148       write(ifhi,'(a,I2)')    'openhisto name SigmaInt-',int(engy)
2149       endif
2150       else
2151       write(ifhi,'(a,I1)')    'openhisto name SigmaInt-',int(engy)
2152       endif
2153       write(ifhi,'(a)')       'htyp lba'
2154       write(ifhi,'(a)')       'xmod log ymod lin'
2155       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2156       write(ifhi,'(a)')'yrange auto auto'
2157       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2158       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!?Int!(X)"'
2159       if (xpar8.eq.1.) then
2160       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2161       endif
2162       write(ifhi,'(a)')       'array 2'
2163
2164       do i=0,nptg
2165         X=xminr
2166         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2167         t=dble(sigma1i(X))
2168         write(ifhi,'(2e14.6)') X,t
2169       enddo
2170
2171       write(ifhi,'(a)')    '  endarray'
2172
2173       if(xpar8.eq.1)then
2174       write(ifhi,'(a)')    'closehisto plot 0-'
2175
2176
2177       if (engy.ge.10.) then
2178        if (engy.ge.100.) then
2179         if (engy.ge.1000.) then
2180          if (engy.ge.10000.) then
2181       write(ifhi,'(a,I5)')    'openhisto name SigmaReelSoft-',int(engy)
2182          else
2183       write(ifhi,'(a,I4)')    'openhisto name SigmaReelSoft-',int(engy)
2184         endif
2185         else
2186       write(ifhi,'(a,I3)')    'openhisto name SigmaReelSoft-',int(engy)
2187        endif
2188        else
2189       write(ifhi,'(a,I2)')    'openhisto name SigmaReelSoft-',int(engy)
2190       endif
2191       else
2192       write(ifhi,'(a,I1)')    'openhisto name SigmaReelSoft-',int(engy)
2193       endif
2194       write(ifhi,'(a)')       'htyp pft'
2195       write(ifhi,'(a)')       'xmod log ymod lin'
2196       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2197       write(ifhi,'(a)')'yrange auto auto'
2198       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2199       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2200       if (xpar8.eq.1.) then
2201       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2202       endif
2203       write(ifhi,'(a)')       'array 2'
2204
2205       do i=0,nptg
2206         X=xminr
2207         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2208         t=dble(sigma2(X,0))
2209         if(t.gt.0.) write(ifhi,'(2e14.6)') X,t
2210       enddo
2211
2212       write(ifhi,'(a)')    '  endarray'
2213       write(ifhi,'(a)')    'closehisto plot 0-'
2214
2215
2216       if (engy.ge.10.) then
2217        if (engy.ge.100.) then
2218         if (engy.ge.1000.) then
2219          if (engy.ge.10000.) then
2220       write(ifhi,'(a,I5)')    'openhisto name SigmaReelSh-',int(engy)
2221          else
2222       write(ifhi,'(a,I4)')    'openhisto name SigmaReelSh-',int(engy)
2223         endif
2224         else
2225       write(ifhi,'(a,I3)')    'openhisto name SigmaReelSh-',int(engy)
2226        endif
2227        else
2228       write(ifhi,'(a,I2)')    'openhisto name SigmaReelSh-',int(engy)
2229       endif
2230       else
2231       write(ifhi,'(a,I1)')    'openhisto name SigmaReelSh-',int(engy)
2232       endif
2233       write(ifhi,'(a)')       'htyp pot'
2234       write(ifhi,'(a)')       'xmod log ymod lin'
2235       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2236       write(ifhi,'(a)')'yrange auto auto'
2237       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2238       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2239       if (xpar8.eq.1.) then
2240       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2241       endif
2242       write(ifhi,'(a)')       'array 2'
2243
2244       do i=0,nptg
2245         X=xminr
2246         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2247         t=dble(sigma2(X,1))
2248         if(t.gt.0.)write(ifhi,'(2e14.6)') X,t
2249       enddo
2250
2251       write(ifhi,'(a)')    '  endarray'
2252
2253       write(ifhi,'(a)')    'closehisto plot 0-'
2254
2255 c**********************************************************************
2256       write(ifhi,'(a)')    'openhisto'
2257       write(ifhi,'(a)')       'htyp lya'
2258       write(ifhi,'(a)')       'xmod log ymod lin'
2259       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2260       write(ifhi,'(a)')'yrange auto auto'
2261       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2262       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2263       if (xpar8.eq.1.) then
2264       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2265       endif
2266       write(ifhi,'(a)')       'array 2'
2267
2268       do i=0,nptg
2269         X=xminr
2270         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2271         t=1.d0/dble(-gamD(0,iclpro,icltar)*log(X*smaxDf)
2272      &     +1./delD(0,iclpro,icltar))
2273         if(t.gt.0.)write(ifhi,'(2e14.6)') X,t
2274       enddo
2275
2276       write(ifhi,'(a)')    '  endarray'
2277
2278
2279       write(ifhi,'(a)')    'closehisto plot 0-'
2280
2281 c**********************************************************************
2282       write(ifhi,'(a)')    'openhisto'
2283       write(ifhi,'(a)')       'htyp lyo'
2284       write(ifhi,'(a)')       'xmod log ymod lin'
2285       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2286       write(ifhi,'(a)')'yrange auto auto'
2287       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2288       write(ifhi,'(a)')    'text 0 0 "yaxis [s]^2!(X)"'
2289       if (xpar8.eq.1.) then
2290       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.75 "Emax=',engy,' GeV"'
2291       endif
2292       write(ifhi,'(a)')       'array 2'
2293
2294       do i=0,nptg
2295         X=xminr
2296         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2297         t=1.d0/dble(-gamD(1,iclpro,icltar)*log(X*smaxDf)
2298      *     +1./delD(1,iclpro,icltar))
2299         if(t.gt.0.)write(ifhi,'(2e14.6)') X,t
2300       enddo
2301
2302       write(ifhi,'(a)')    '  endarray'
2303
2304       endif
2305
2306
2307       write(ifhi,'(a)')    'closehisto plot 0'
2308
2309 c**********************************************************************
2310 c**********************************************************************
2311
2312       if (engy.ge.10.) then
2313        if (engy.ge.100.) then
2314         if (engy.ge.1000.) then
2315          if (engy.ge.10000.) then
2316       write(ifhi,'(a,I5)')    'openhisto name SigmaDiff-',int(engy)
2317          else
2318       write(ifhi,'(a,I4)')    'openhisto name SigmaDiff-',int(engy)
2319         endif
2320         else
2321       write(ifhi,'(a,I3)')    'openhisto name SigmaDiff-',int(engy)
2322        endif
2323        else
2324       write(ifhi,'(a,I2)')    'openhisto name SigmaDiff-',int(engy)
2325       endif
2326       else
2327       write(ifhi,'(a,I1)')    'openhisto name SigmaDiff-',int(engy)
2328       endif
2329       write(ifhi,'(a)')       'htyp lin'
2330       write(ifhi,'(a)')       'xmod log ymod lin'
2331       write(ifhi,'(a,2e11.3)')'xrange',xminr,1.
2332       write(ifhi,'(a,2e11.3)')'yrange auto auto'
2333       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2334       write(ifhi,'(a)')    'text 0 0 "yaxis [D][s]/[s]"'
2335       if (xpar8.eq.1.) then
2336       write(ifhi,'(a,e7.2,a)')  'text 0.05 0.9 "Emax=',engy,' GeV"'
2337       endif
2338       write(ifhi,'(a)')       'array 2'
2339
2340       do i=0,nptg
2341         X=xminr
2342         if (i.ne.0) X=X*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2343         t=0.d0
2344         if(w(i).gt.0.d0)t=(z(i)-w(i))/w(i)
2345         if(abs(t).gt.0.15d0) t=dsign(0.15d0,t)
2346         write(ifhi,'(2e14.6)') X,t
2347       enddo
2348
2349       write(ifhi,'(a)')    '  endarray'
2350       write(ifhi,'(a)')    'closehisto plot 0'
2351
2352       end
2353
2354 c----------------------------------------------------------------------
2355       subroutine xParGauss
2356 c----------------------------------------------------------------------
2357
2358       include 'epos.inc'
2359       include 'epos.incpar'
2360
2361       double precision x,om5s,xDfit,y,enh!,omNpuncut
2362
2363       nptg=50                  !number of point for the graphs
2364       x=dble(xpar4)                  !value of x (energy)
2365       y=dble(xpar2)                  !value of rapidity
2366
2367       if (engy.ge.10.) then
2368        if (engy.ge.100.) then
2369         if (engy.ge.1000.) then
2370          if (engy.ge.10000.) then
2371       write(ifhi,'(a,I5)')    'openhisto name GaussExact-',int(engy)
2372            else
2373       write(ifhi,'(a,I4)')    'openhisto name GaussExact-',int(engy)
2374         endif
2375         else
2376       write(ifhi,'(a,I3)')    'openhisto name GaussExact-',int(engy)
2377        endif
2378        else
2379       write(ifhi,'(a,I2)')    'openhisto name GaussExact-',int(engy)
2380       endif
2381       else
2382       write(ifhi,'(a,I1)')    'openhisto name GaussExact-',int(engy)
2383       endif
2384       write(ifhi,'(a)')       'htyp lru'
2385       write(ifhi,'(a)')       'xmod lin ymod lin'
2386       write(ifhi,'(a,2e11.3)')'xrange',0.,bmaxDf
2387       write(ifhi,'(a)')'yrange auto auto'
2388       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2389       write(ifhi,'(a)')    'text 0 0 "yaxis b*D(x+,x-,s,b)"'
2390 c      if (xpar8.eq.1.) then
2391       write(ifhi,'(a,e7.2,a)')  'text 0.6 0.9 "E=',engy,' GeV"'
2392       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.8 "x=',x,'"'
2393       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.7 "y=',y,'"'
2394 c      endif
2395       write(ifhi,'(a)')       'array 2'
2396
2397       do i=0,nptg
2398         b=bmaxDf*(real(i)/real(nptg))
2399         enh=0.d0
2400        write(ifhi,*) b,dble(b)*(2*om5s(real(x)*smaxDf,x,y,b,-1,-1)
2401      &                +enh)
2402      &         /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
2403       enddo
2404
2405       write(ifhi,'(a)')    '  endarray'
2406       write(ifhi,'(a)')    'closehisto plot 0-'
2407
2408
2409
2410 c**********************************************************************
2411
2412       if (engy.ge.10.) then
2413        if (engy.ge.100.) then
2414         if (engy.ge.1000.) then
2415          if (engy.ge.10000.) then
2416       write(ifhi,'(a,I5)')    'openhisto name GaussParam-',int(engy)
2417          else
2418       write(ifhi,'(a,I4)')    'openhisto name GaussParam-',int(engy)
2419         endif
2420         else
2421       write(ifhi,'(a,I3)')    'openhisto name GaussParam-',int(engy)
2422        endif
2423        else
2424       write(ifhi,'(a,I2)')    'openhisto name GaussParam-',int(engy)
2425       endif
2426       else
2427       write(ifhi,'(a,I1)')    'openhisto name GaussParam-',int(engy)
2428       endif
2429       write(ifhi,'(a)')       'htyp pfc'
2430       write(ifhi,'(a)')       'xmod log ymod lin'
2431       write(ifhi,'(a,2e11.3)')'xrange',0.,bmaxDf
2432       write(ifhi,'(a)')'yrange auto auto'
2433       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2434       write(ifhi,'(a)')    'text 0 0 "yaxis b*D(x+,x-,s,b)"'
2435 c      if (xpar8.eq.1.) then
2436       write(ifhi,'(a,e7.2,a)')  'text 0.6 0.9 "E=',engy,' GeV"'
2437       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.8 "x=',x,'"'
2438       write(ifhi,'(a,f5.2,a)')  'text 0.6 0.7 "y=',y,'"'
2439 c      endif
2440       write(ifhi,'(a)')       'array 2'
2441
2442       imax=idxD1
2443       if(iomega.eq.2)imax=1
2444
2445       do i=0,nptg
2446         b=bmaxDf*(real(i)/real(nptg))
2447         xp=sqrt(real(x))*exp(real(y))
2448         xm=sqrt(real(x))*exp(-real(y))
2449         write(ifhi,'(2e14.6)') b,dble(b)*
2450      &       xDfit(0,imax,engy**2,xp,xm,b)
2451       enddo
2452
2453       write(ifhi,'(a)')    '  endarray'
2454       write(ifhi,'(a)')    'closehisto plot 0'
2455
2456       end
2457
2458
2459
2460 c----------------------------------------------------------------------
2461       subroutine xParOmega1
2462 c----------------------------------------------------------------------
2463
2464       include 'epos.inc'
2465       include 'epos.incpar'
2466
2467       double precision x,w(0:200),z(0:200)
2468       double precision yp,om1,xminr,Dsoftshval,t
2469
2470       nptg=50                 !number of point for the graphs
2471       biniDf=xpar2               !value of biniDf (impact parameter)
2472       yp=xpar3                   !value of yp (rapidity)
2473       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
2474
2475       if (engy.ge.10.) then
2476        if (engy.ge.100.) then
2477         if (engy.ge.1000.) then
2478          if (engy.ge.10000.) then
2479       write(ifhi,'(a,I5)')    'openhisto name Om1Exact-',int(engy)
2480          else
2481       write(ifhi,'(a,I4)')    'openhisto name Om1Exact-',int(engy)
2482         endif
2483         else
2484       write(ifhi,'(a,I3)')    'openhisto name Om1Exact-',int(engy)
2485        endif
2486        else
2487       write(ifhi,'(a,I2)')    'openhisto name Om1Exact-',int(engy)
2488       endif
2489       else
2490       write(ifhi,'(a,I1)')    'openhisto name Om1Exact-',int(engy)
2491       endif
2492       write(ifhi,'(a)')       'htyp lru'
2493       write(ifhi,'(a)')       'xmod log ymod log'
2494       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2495       write(ifhi,'(a)')       'yrange auto auto'
2496       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2497       write(ifhi,'(a)')    'text 0 0 "yaxis [h](x,0,b)"'
2498       if (xpar8.eq.1.) then
2499       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2500       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
2501       endif
2502       write(ifhi,'(a)')       'array 2'
2503
2504       do i=0,nptg
2505         x=xminr
2506         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2507         w(i)=Dsoftshval(real(x)*engy**2,x,0.d0,biniDf,0)
2508      &       *(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
2509         write(ifhi,*) x,w(i)
2510       enddo
2511
2512       write(ifhi,'(a)')    '  endarray'
2513       write(ifhi,'(a)')    'closehisto plot 0-'
2514
2515
2516 c**********************************************************************
2517
2518       if (engy.ge.10.) then
2519        if (engy.ge.100.) then
2520         if (engy.ge.1000.) then
2521          if (engy.ge.10000.) then
2522       write(ifhi,'(a,I5)')    'openhisto name om5param-',int(engy)
2523          else
2524       write(ifhi,'(a,I4)')    'openhisto name om5param-',int(engy)
2525         endif
2526         else
2527       write(ifhi,'(a,I3)')    'openhisto name om5param-',int(engy)
2528        endif
2529        else
2530       write(ifhi,'(a,I2)')    'openhisto name om5param-',int(engy)
2531       endif
2532       else
2533       write(ifhi,'(a,I1)')    'openhisto name om5param-',int(engy)
2534       endif
2535       write(ifhi,'(a)')       'htyp pfc'
2536       write(ifhi,'(a)')       'xmod log ymod log'
2537       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2538       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2539       write(ifhi,'(a)')    'yrange auto auto'
2540       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
2541       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1!(x,0,b)"'
2542       if (xpar8.eq.1.) then
2543       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2544       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
2545       endif
2546       write(ifhi,'(a)')       'array 2'
2547
2548       do i=0,nptg
2549         x=xminr
2550         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2551         z(i)=om1(x,yp,biniDf)
2552         write(ifhi,*) x,z(i)
2553       enddo
2554
2555       write(ifhi,'(a)')    '  endarray'
2556       write(ifhi,'(a)')    'closehisto plot 0'
2557
2558
2559 c**********************************************************************
2560 c**********************************************************************
2561
2562       if (engy.ge.10.) then
2563        if (engy.ge.100.) then
2564         if (engy.ge.1000.) then
2565          if (engy.ge.10000.) then
2566       write(ifhi,'(a,I5)')    'openhisto name Om1Diff-',int(engy)
2567          else
2568       write(ifhi,'(a,I4)')    'openhisto name Om1Diff-',int(engy)
2569         endif
2570         else
2571       write(ifhi,'(a,I3)')    'openhisto name Om1Diff-',int(engy)
2572        endif
2573        else
2574       write(ifhi,'(a,I2)')    'openhisto name Om1Diff-',int(engy)
2575       endif
2576       else
2577       write(ifhi,'(a,I1)')    'openhisto name Om1Diff-',int(engy)
2578       endif
2579       write(ifhi,'(a)')       'htyp lru'
2580       write(ifhi,'(a)')       'xmod log ymod lin'
2581       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2582       write(ifhi,'(a,2e11.3)')'yrange auto auto'
2583       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
2584       write(ifhi,'(a)')    'text 0 0 "yaxis ([w]?1!-[h])/[h]"'
2585       if (xpar8.eq.1.) then
2586       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2587       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
2588       endif
2589       write(ifhi,'(a)')       'array 2'
2590
2591       do i=0,nptg
2592         x=xminr
2593         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2594          t=z(i)!(z(i)-w(i))
2595 c         if(abs(w(i)).gt.0.)t=t/w(i)
2596 c         if(abs(t).gt.0.15d0) t=dsign(0.15d0,t)
2597          write(ifhi,*) x,t
2598       enddo
2599
2600       write(ifhi,'(a)')    '  endarray'
2601       write(ifhi,'(a)')    'closehisto plot 0'
2602
2603
2604       end
2605
2606 c----------------------------------------------------------------------
2607       subroutine xEpsilon(iii)
2608 c----------------------------------------------------------------------
2609 c iii:  modus (0,1,2)
2610 c----------------------------------------------------------------------
2611       include 'epos.inc'
2612       include 'epos.incsem'
2613       include 'epos.incems'
2614       include 'epos.incpar'
2615
2616       parameter(nxeps=20,nyeps=32)
2617       common/cxeps1/w(0:nxeps,nyeps),y1(nyeps),y2(nyeps)
2618       common/cxeps2/db,b1,b2
2619       common/geom/rmproj,rmtarg,bmax,bkmx
2620       character ch*3
2621
2622       b1=0.03
2623       b2=bkmx*1.2
2624       db=(b2-b1)/nyeps
2625
2626         if(iii.eq.0)then
2627
2628       do j=0,nxeps
2629        do k=1,nyeps
2630         w(j,k)=0
2631        enddo
2632       enddo
2633
2634         elseif(iii.eq.2)then
2635
2636       do nj=1,14
2637        y1(nj)=1e+20
2638        y2(nj)=1e-20
2639        do k=1,nyeps
2640         if(w(0,k).ne.0)then
2641          y1(nj)=min(y1(nj),w(nj,k)/w(0,k))
2642          y2(nj)=max(y2(nj),w(nj,k)/w(0,k))
2643         endif
2644        enddo
2645        y1(nj)=max(y1(nj)*.2,1e-4)
2646        y2(nj)=min(y2(nj)*5,1e4)
2647       enddo
2648       y2(13)=max(y2(13),y2(14))
2649       y2(14)=max(y2(13),y2(14))
2650       y2(1)=max(y2(1),y2(2))
2651       y2(2)=max(y2(1),y2(2))
2652       y2(5)=max(y2(5),y2(6))
2653       y2(6)=max(y2(5),y2(6))
2654       y2(7)=y2(5)
2655       y2(8)=y2(5)
2656       y2(9)=max(y2(9),y2(10))
2657       y2(10)=max(y2(9),y2(10))
2658       y2(11)=y2(9)
2659       y2(12)=y2(9)
2660       do nj=1,14
2661        if(nj.le.9)write(ifhi,'(a,i1)')'openhisto name xEps',nj
2662        if(nj.gt.9)write(ifhi,'(a,i2)')'openhisto name xEps',nj
2663        ch='lin'
2664        if(nj.eq.7.or.nj.eq.11)ch='lyo'
2665        if(nj.eq.8.or.nj.eq.12)ch='lgo'
2666        write(ifhi,'(a)')     'htyp '//ch//' xmod lin ymod log'
2667        write(ifhi,'(a,e9.2)')'xrange 0 ',b2
2668        if(nj.eq.1.or.nj.eq.3.or.nj.eq.5.or.nj.eq.9)then
2669        write(ifhi,'(a,2e9.2)')     'yrange ',min(y1(nj),y1(nj+1))
2670      *                                      ,max(y2(nj),y2(nj+1))
2671        else
2672        write(ifhi,'(a,2e9.2)')     'yrange ',y1(nj),y2(nj)
2673        endif
2674        write(ifhi,'(a)')     'text 0 0 "xaxis b"'
2675        if(nj.eq.1) write(ifhi,'(a)')'txt "yaxis [e]?GP/T!(b)"'
2676        if(nj.eq.1) write(ifhi,'(a)')'txt "title soft pro   soft tar"'
2677        if(nj.eq.3) write(ifhi,'(a)')'txt "yaxis [e]?G!(b)"'
2678        if(nj.eq.3) write(ifhi,'(a)')'txt "title soft   semi"'
2679        if(nj.eq.5) write(ifhi,'(a)')'txt "yaxis [b]?eff!(b)"'
2680        if(nj.eq.5) write(ifhi,'(a)')'txt "title soft pro   soft tar"'
2681        if(nj.eq.9) write(ifhi,'(a)')'txt "yaxis [b]?eff!(b)"'
2682        if(nj.eq.9) write(ifhi,'(a)')'txt "title semi pro   semi tar"'
2683        if(nj.eq.13)write(ifhi,'(a)')'txt "yaxis Z?P/T!"'
2684        write(ifhi,'(a)')       'array 2'
2685        do k=1,nyeps
2686         b=b1+(k-0.5)*db
2687         y=0
2688         if(w(0,k).ne.0)y=w(nj,k)/w(0,k)
2689         write(ifhi,'(2e11.3)')b,y
2690        enddo
2691        write(ifhi,'(a)')    '  endarray'
2692        if(nj.eq.2.or.nj.eq.4.or.nj.eq.8.or.nj.eq.12.or.nj.eq.14
2693      &    .or.nj.eq.16)then
2694         write(ifhi,'(a)')    'closehisto plot 0'
2695        else
2696         write(ifhi,'(a)')    'closehisto plot 0-'
2697        endif
2698       enddo
2699       !----15-16-17---
2700       write(ifhi,'(a)') 'openhisto name xEps15'
2701       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2702       write(ifhi,'(a)') 'xrange 0 10'
2703       write(ifhi,'(a)') 'text 0 0 "xaxis b?0!"'
2704       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(b?0!)"'
2705       write(ifhi,'(a)') 'array 2'
2706       do k=1,10
2707        b=(k-0.5)
2708        y=0
2709        if(w(17,k).ne.0)y=w(15,k)/w(17,k)
2710        write(ifhi,'(2e11.3)')b,y
2711       enddo
2712       write(ifhi,'(a)')  '  endarray'
2713       write(ifhi,'(a)')  'closehisto plot 0-'
2714
2715       write(ifhi,'(a)') 'openhisto name xEps16'
2716       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2717       write(ifhi,'(a)') 'xrange 0 10'
2718       write(ifhi,'(a)') 'text 0 0 "xaxis b?0!"'
2719       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(b?0!)"'
2720       write(ifhi,'(a)') 'array 2'
2721       do k=1,10
2722        b=(k-0.5)
2723        y=0
2724        if(w(17,k).ne.0)y=w(16,k)/w(17,k)
2725        write(ifhi,'(2e11.3)')b,y
2726       enddo
2727       write(ifhi,'(a)')  '  endarray'
2728       write(ifhi,'(a)')  'closehisto plot 0'
2729       !----18-19-20---
2730       kk=2
2731       do k=3,32
2732         if(w(18,k).ne.0)kk=k
2733       enddo
2734       xmx=(kk-1)/31.*0.1*maproj*matarg
2735       write(ifhi,'(a)') 'openhisto name xEps18'
2736       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2737       write(ifhi,'(a,f10.2)') 'xrange 0 ',xmx
2738       write(ifhi,'(a)') 'text 0 0 "xaxis n?Gl!"'
2739       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(n?Gl!)"'
2740       write(ifhi,'(a)') 'array 2'
2741       do k=1,32
2742        x=(k-1.)*0.1*maproj*matarg/(nyeps-1.)
2743        y=0
2744        if(w(20,k).ne.0)y=w(18,k)/w(20,k)
2745        write(ifhi,'(2e11.3)')x,y
2746       enddo
2747       write(ifhi,'(a)')  '  endarray'
2748       write(ifhi,'(a)')  'closehisto plot 0-'
2749
2750       write(ifhi,'(a)') 'openhisto name xEps19'
2751       write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
2752       write(ifhi,'(a,f10.2)') 'xrange 0 ',xmx
2753       write(ifhi,'(a)') 'text 0 0 "xaxis n?Gl!"'
2754       write(ifhi,'(a)') 'txt "yaxis Z?P/T!(n?Gl!)"'
2755       write(ifhi,'(a)') 'array 2'
2756       do k=1,32
2757        x=(k-1.)*0.1*maproj*matarg/(nyeps-1.)
2758        y=0
2759        if(w(20,k).ne.0)y=w(19,k)/w(20,k)
2760        write(ifhi,'(2e11.3)')x,y
2761       enddo
2762       write(ifhi,'(a)')  '  endarray'
2763       write(ifhi,'(a)')  'closehisto plot 0'
2764
2765         endif
2766
2767       end
2768
2769 c----------------------------------------------------------------------
2770       subroutine xParOmegaN
2771 c----------------------------------------------------------------------
2772 c xpar1=engy
2773 c xpar2=b
2774 c xpar4=xremnant
2775 c xpar5: 0=log scale (x dep of om) 1=lin scale (b dep of om)
2776 c----------------------------------------------------------------------
2777
2778       include 'epos.inc'
2779       include 'epos.incpar'
2780       include 'epos.incsem'
2781       include 'epos.incems'
2782       double precision x,w(0:200),z(0:200),xminr,t,ghh
2783      *,xprem,omGam,Womint,Gammapp,WomGamint,omGamk
2784 c     *,yp,SomY,omYuncut,y,xtmp
2785
2786       nptg=30                  !number of point for the graphs
2787       biniDf=xpar2               !value of biniDf (impact parameter)
2788       xprem=dble(xpar4)            !value of x remnant
2789       bmax=3.
2790       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
2791
2792 c**********************************************************************
2793
2794       do i=0,3
2795         b=bmax*(real(i)/real(3))
2796         z(i)=1.d0
2797 c        if(xpar5.eq.0.)z(i)=Gammapp(engy**2.,b,1)
2798       enddo
2799
2800       write(ifhi,'(a)')    'openhisto name Womint-1'
2801       write(ifhi,'(a)')       'htyp lru'
2802       if(xpar5.eq.0.)then
2803         write(ifhi,'(a)')       'xmod lin ymod log'
2804       else
2805         write(ifhi,'(a)')       'xmod lin ymod lin'
2806       endif
2807       write(ifhi,'(a,2e11.3)')'xrange',0.,bmax
2808       write(ifhi,'(a)')    'yrange auto auto'
2809       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2810       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?int!(s,b)"'
2811       if (xpar8.eq.1.) then
2812       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2813       endif
2814       write(ifhi,'(a)')       'array 2'
2815
2816       ierr=0
2817       ghh=0.d0
2818       do i=0,3
2819         b=bmax*(real(i)/real(3))
2820         w(i)=Womint(engy**2.,b)
2821         if(b.eq.biniDf)then
2822           write(*,*)'Womint(',b,',1)=',w(i)
2823           ghh=ghh+w(i)
2824         endif
2825         if(w(i).lt.0.d0.and.ierr.eq.0.and.xpar5.eq.0.)then
2826           write(*,*)'Warning Womint(1)<0 =',w(i)
2827           w(i)=-w(i)
2828           ierr=1
2829           elseif(w(i).lt.0.d0.and.xpar5.eq.0.)then
2830             w(i)=-w(i)
2831           elseif(w(i).ge.0.d0.and.ierr.eq.1.and.xpar5.eq.0.)then
2832             ierr=0
2833             write(*,*)'Warning Womint(1)>0 =',w(i)
2834         endif
2835         write(ifhi,*) b,w(i)/z(i)
2836       enddo
2837
2838
2839       write(ifhi,'(a)')    '  endarray'
2840       write(ifhi,'(a)')    'closehisto plot 0-'
2841
2842 c**********************************************************************
2843
2844       write(ifhi,'(a)')    'openhisto'
2845       write(ifhi,'(a)')       'htyp pfc'
2846       if(xpar5.eq.0.)then
2847         write(ifhi,'(a)')       'xmod lin ymod log'
2848       else
2849         write(ifhi,'(a)')       'xmod lin ymod lin'
2850       endif
2851       write(ifhi,'(a,2e11.3)')'xrange',0.,bmax
2852       write(ifhi,'(a)')    'yrange auto auto'
2853       write(ifhi,'(a)')    'text 0 0 "xaxis b"'
2854       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?int!(s,b)"'
2855       if (xpar8.eq.1.) then
2856       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2857       endif
2858       write(ifhi,'(a)')       'array 2'
2859
2860       do i=0,nptg
2861         b=bmax*(real(i)/real(nptg))
2862         z(i)=1.d0
2863         if(xpar5.eq.0.)then
2864           z(i)=z(i)+dabs(WomGamint(b))
2865         endif
2866       enddo
2867       ierr=0
2868       do i=0,nptg
2869         b=bmax*(real(i)/real(nptg))
2870         w(i)=WomGamint(b)
2871         if(w(i).lt.0.d0.and.ierr.eq.0.and.xpar5.eq.0.)then
2872           write(*,*)'Warning WomGamint(1)<0 =',w(i)
2873           w(i)=-w(i)
2874           ierr=1
2875           elseif(w(i).lt.0.d0.and.xpar5.eq.0.)then
2876             w(i)=-w(i)
2877           elseif(w(i).ge.0.d0.and.ierr.eq.1.and.xpar5.eq.0.)then
2878             ierr=0
2879             write(*,*)'Warning WomGamint(1)>0 =',w(i)
2880         endif
2881         write(ifhi,*) b,w(i)/z(i)
2882       enddo
2883
2884
2885       write(ifhi,'(a)')    '  endarray'
2886       write(ifhi,'(a)')    'closehisto plot 0'
2887       t=Gammapp(engy**2.,biniDf,1)
2888       write(*,*)'--> gamma(',biniDf,')=',ghh,t
2889
2890 c**********************************************************************
2891 c**********************************************************************
2892       do k=1,koll
2893         bk(k)=biniDf
2894       enddo
2895       call GfunPark(0)
2896       call integom1(0)
2897
2898       if (engy.ge.10.) then
2899        if (engy.ge.100.) then
2900         if (engy.ge.1000.) then
2901          if (engy.ge.10000.) then
2902       write(ifhi,'(a,I5)')    'openhisto name xOmNG-',int(engy)
2903          else
2904       write(ifhi,'(a,I4)')    'openhisto name xOmNG-',int(engy)
2905         endif
2906         else
2907       write(ifhi,'(a,I3)')    'openhisto name xOmNG-',int(engy)
2908        endif
2909        else
2910       write(ifhi,'(a,I2)')    'openhisto name xOmNG-',int(engy)
2911       endif
2912       else
2913       write(ifhi,'(a,I1)')    'openhisto name xOmNG-',int(engy)
2914       endif
2915       write(ifhi,'(a)')       'htyp lru'
2916       write(ifhi,'(a)')       'xmod log ymod log'
2917       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2918       write(ifhi,'(a)')    'yrange auto auto'
2919       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
2920       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?NPi!(x+rem,x-rem,b)"'
2921       if (xpar8.eq.1.) then
2922       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2923       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
2924       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.85 "x+?rem!=',xprem,'"'
2925       endif
2926       write(ifhi,'(a)')       'array 2'
2927
2928       do i=0,nptg
2929         x=xminr
2930         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2931         t=dabs(omGam(xprem,x,biniDf))
2932         write(ifhi,*) x,t
2933       enddo
2934
2935       write(ifhi,'(a)')    '  endarray'
2936       write(ifhi,'(a)')    'closehisto plot 0-'
2937
2938       if (engy.ge.10.) then
2939        if (engy.ge.100.) then
2940         if (engy.ge.1000.) then
2941          if (engy.ge.10000.) then
2942       write(ifhi,'(a,I5)')    'openhisto name xOmNG-',int(engy)
2943          else
2944       write(ifhi,'(a,I4)')    'openhisto name xOmNG-',int(engy)
2945         endif
2946         else
2947       write(ifhi,'(a,I3)')    'openhisto name xOmNG-',int(engy)
2948        endif
2949        else
2950       write(ifhi,'(a,I2)')    'openhisto name xOmNG-',int(engy)
2951       endif
2952       else
2953       write(ifhi,'(a,I1)')    'openhisto name xOmNG-',int(engy)
2954       endif
2955       write(ifhi,'(a)')       'htyp pfc'
2956       write(ifhi,'(a)')       'xmod log ymod log'
2957       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
2958       write(ifhi,'(a)')    'yrange auto auto'
2959       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
2960       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?NPi!(x+rem,x-rem,b)"'
2961       if (xpar8.eq.1.) then
2962       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
2963       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
2964       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.85 "x+?rem!=',xprem,'"'
2965       endif
2966       write(ifhi,'(a)')       'array 2'
2967
2968       do i=0,nptg
2969         x=xminr
2970         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
2971         t=dabs(omGamk(1,xprem,x))
2972         write(ifhi,*) x,t
2973       enddo
2974
2975       write(ifhi,'(a)')    '  endarray'
2976
2977       write(ifhi,'(a)')    'closehisto plot 0'
2978
2979       end
2980
2981 c----------------------------------------------------------------------
2982       subroutine xParGampp
2983 c----------------------------------------------------------------------
2984
2985       include 'epos.inc'
2986       include 'epos.incpar'
2987
2988       double precision Gammapp,GammaGauss,sg,sgmc,Znorm,Zn,t
2989      *                 ,w(0:200),z(0:200)!,GammaMC
2990
2991       nptg=2                  !number of point for the graphs
2992       biniDf=xpar2              !value of biniDf (impact parameter)
2993
2994
2995 c**************************************************
2996
2997       if (biniDf.lt.1.) then
2998         k=int(10.*biniDf)
2999         write(ifhi,'(a,I1)')  'openhisto name Gamma-b0.',k
3000       else
3001         write(ifhi,'(a,f3.1)')'openhisto name Gamma-b',biniDf
3002       endif
3003       write(ifhi,'(a)')       'htyp lru'
3004       write(ifhi,'(a)')       'xmod lin ymod lin'
3005       write(ifhi,'(a,2e11.3)')'xrange',0.,real(nptg)
3006       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3007       write(ifhi,'(a)')    'text 0 0 "xaxis m"'
3008       write(ifhi,'(a)')    'text 0 0 "yaxis [g]?h1h2!"'
3009       if (xpar8.eq.1.) then
3010         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3011         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3012       endif
3013       write(ifhi,'(a)')       'array 2'
3014       sg=0.d0
3015       do i=0,nptg
3016         w(i)=Gammapp(engy**2.,biniDf,i)
3017         sg=sg+w(i)
3018         write(ifhi,*) i,w(i)
3019       write(*,*) 'G12',i,w(i)
3020       enddo
3021
3022       write(ifhi,'(a)')    '  endarray'
3023       write(ifhi,'(a)')    'closehisto plot 0-'
3024
3025       if (biniDf.lt.1.) then
3026         k=int(10.*biniDf)
3027         write(ifhi,'(a,I1)')  'openhisto name GammaMC-b0.',k
3028       else
3029         write(ifhi,'(a,f3.1)')'openhisto name GammaMC-b',biniDf
3030       endif
3031       write(ifhi,'(a)')       'htyp pfc'
3032       write(ifhi,'(a)')       'xmod lin ymod lin'
3033       write(ifhi,'(a,2e11.3)')'xrange',0.,real(nptg)
3034       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3035       write(ifhi,'(a)')    'text 0 0 "xaxis m"'
3036       write(ifhi,'(a)')    'text 0 0 "yaxis [g]?h1h2!"'
3037       if (xpar8.eq.1.) then
3038         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3039         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3040       endif
3041       write(ifhi,'(a)')       'array 2'
3042
3043       sgmc=0.d0
3044       do i=0,nptg
3045         z(i)=GammaGauss(engy**2.,biniDf,i)
3046         sgmc=sgmc+z(i)
3047         write(ifhi,*) i,z(i)
3048         write(*,*) 'G12gauss',i,z(i)
3049       enddo
3050
3051       write(ifhi,'(a)')    '  endarray'
3052       write(ifhi,'(a)')    'closehisto plot 0'
3053
3054 c**********************************************************************
3055 c**********************************************************************
3056       Zn=Znorm(engy**2,biniDf)
3057
3058       write(ifhi,'(a)')    'openhisto name GammaDiff'
3059       write(ifhi,'(a)')       'htyp lru'
3060       write(ifhi,'(a)')       'xmod lin ymod lin'
3061       write(ifhi,'(a,2e11.3)')'xrange',0.,real(nptg)
3062       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3063       write(ifhi,'(a)')    'text 0 0 "xaxis m"'
3064       write(ifhi,'(a)')    'text 0 0 "yaxis (G-GMC)/G"'
3065       if (xpar8.eq.1.) then
3066       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3067       write(ifhi,'(a,f5.2,a)')  'text 0.1 0.8 "b=',biniDf,' fm"'
3068       endif
3069       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "[S]?Guncut!=',sg,'"'
3070       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "[S]?Gcut!=',sgmc,'"'
3071       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.7 "Z=',Zn,'"'
3072       write(ifhi,'(a)')       'array 2'
3073
3074       do i=0,nptg
3075         if(w(i).ne.0d0) t=(z(i)-w(i))/w(i)
3076 c         if(abs(t).gt.0.5d0) t=dsign(0.5d0,t)
3077          write(ifhi,*) i,t
3078       enddo
3079
3080       write(ifhi,'(a)')    '  endarray'
3081       write(ifhi,'(a)')    'closehisto plot 0'
3082
3083       end
3084
3085 c----------------------------------------------------------------------
3086       subroutine xParPomInc
3087 c----------------------------------------------------------------------
3088
3089       include 'epos.inc'
3090       include 'epos.incpar'
3091
3092       double precision x,PomIncExact,PomIncUnit,xminr,xm,t
3093      *                 ,w(0:200),z(0:200)
3094
3095       nptg=10                  !number of point for the graphs
3096       biniDf=xpar2              !value of biniDf (impact parameter)
3097       xm=dble(xpar4)              !value of biniDf (impact parameter)
3098       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3099
3100
3101 c********************* red = PomIncXExact *****************************
3102
3103       if (biniDf.lt.1.) then
3104         k=int(10.*biniDf)
3105         write(ifhi,'(a,I1)')  'openhisto name PomIncExact-b0.',k
3106       else
3107         write(ifhi,'(a,f3.1)')'openhisto name PomIncExact-b',biniDf
3108       endif
3109       write(ifhi,'(a)')       'htyp lru'
3110       write(ifhi,'(a)')       'xmod log ymod lin'
3111       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3112       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3113       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3114       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+/dx-"'
3115       if (xpar8.eq.1.) then
3116         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3117         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3118         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x-=',xm,'"'
3119       endif
3120       write(ifhi,'(a)')       'array 2'
3121
3122       do i=0,nptg
3123         x=xminr
3124         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3125 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3126         w(i)=PomIncExact(dsqrt(x),dsqrt(x),biniDf)
3127         write(ifhi,*) x,w(i)
3128       write(*,*) 'Xe',i,w(i)
3129       enddo
3130
3131       write(ifhi,'(a)')    '  endarray'
3132       write(ifhi,'(a)')    'closehisto plot 0-'
3133
3134
3135
3136 c************************* dot = PomIncXUnit **************************
3137 c      nptg=50     !number of point for the graphs
3138
3139       if (biniDf.lt.1.) then
3140         k=int(10.*biniDf)
3141         write(ifhi,'(a,I1)')  'openhisto name PomIncUnit-b0.',k
3142       else
3143         write(ifhi,'(a,f3.1)')'openhisto name PomIncUnit-b',biniDf
3144       endif
3145       write(ifhi,'(a)')       'htyp pfc'
3146       write(ifhi,'(a)')       'xmod log ymod lin'
3147       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3148       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3149       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3150       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+/dx-"'
3151       if (xpar8.eq.1.) then
3152         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3153         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3154         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x-=',xm,' fm"'
3155       endif
3156       write(ifhi,'(a)')       'array 2'
3157
3158       do i=0,nptg
3159         x=xminr
3160         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3161 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3162         z(i)=PomIncUnit(dsqrt(x),dsqrt(x),biniDf)
3163         write(ifhi,*) x,z(i)
3164         write(*,*) 'Xu',i,z(i)
3165       enddo
3166
3167       write(ifhi,'(a)')    '  endarray'
3168       write(ifhi,'(a)')    'closehisto plot 0'
3169
3170 c**********************************************************************
3171 c**********************************************************************
3172
3173       write(ifhi,'(a)')    'openhisto name PomIncDiff'
3174       write(ifhi,'(a)')       'htyp lru'
3175       write(ifhi,'(a)')       'xmod log ymod lin'
3176       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3177       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3178       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3179       write(ifhi,'(a)')    'text 0 0 "yaxis ([w]?5!-G)/G"'
3180       if (xpar8.eq.1.) then
3181       write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3182       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3183       write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x-=',xm,'"'
3184       endif
3185       write(ifhi,'(a)')       'array 2'
3186
3187       do i=0,nptg
3188         x=xminr
3189         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3190 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3191         t=0.d0
3192         if(w(i).ne.0d0) t=(z(i)-w(i))/w(i)
3193 c         if(abs(t).gt.0.5d0) t=dsign(0.5d0,t)
3194          write(ifhi,*) x,t
3195       enddo
3196
3197       write(ifhi,'(a)')    '  endarray'
3198       write(ifhi,'(a)')    'closehisto plot 0'
3199
3200       end
3201
3202 c----------------------------------------------------------------------
3203       subroutine xParPomIncX
3204 c----------------------------------------------------------------------
3205
3206       include 'epos.inc'
3207       include 'epos.incpar'
3208
3209       double precision x,PomIncXExact,PomIncXUnit,xminr,y
3210
3211       nptg=20                   !number of point for the graphs
3212       biniDf=xpar2              !value of biniDf (impact parameter)
3213       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3214
3215
3216 c********************* red = PomIncXExact *****************************
3217
3218       if (biniDf.lt.1.) then
3219         k=int(10.*biniDf)
3220         write(ifhi,'(a,I1)')  'openhisto name PomIncXExact-b0.',k
3221       else
3222         write(ifhi,'(a,f3.1)')'openhisto name PomIncXExact-b',biniDf
3223       endif
3224       write(ifhi,'(a)')       'htyp lru'
3225       write(ifhi,'(a)')       'xmod log ymod lin'
3226       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3227       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3228       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3229       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx(x,b)"'
3230       if (xpar8.eq.1.) then
3231         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3232         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3233       endif
3234       write(ifhi,'(a)')       'array 2'
3235
3236       do i=0,nptg
3237         x=xminr
3238         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3239 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3240         y=PomIncXExact(x,biniDf)
3241         write(ifhi,*) x,y
3242 c      write(*,*) 'Xe',i
3243       enddo
3244
3245       write(ifhi,'(a)')    '  endarray'
3246       write(ifhi,'(a)')    'closehisto plot 0-'
3247
3248
3249
3250 c************************* dot = PomIncXUnit **************************
3251 c      nptg=50     !number of point for the graphs
3252
3253       if (biniDf.lt.1.) then
3254         k=int(10.*biniDf)
3255         write(ifhi,'(a,I1)')  'openhisto name PomIncXUnit-b0.',k
3256       else
3257         write(ifhi,'(a,f3.1)')'openhisto name PomIncXUnit-b',biniDf
3258       endif
3259       write(ifhi,'(a)')       'htyp pfc'
3260       write(ifhi,'(a)')       'xmod log ymod lin'
3261       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3262       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3263       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3264       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx(x,b)"'
3265       if (xpar8.eq.1.) then
3266         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3267         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3268       endif
3269       write(ifhi,'(a)')       'array 2'
3270
3271       do i=0,nptg
3272         x=xminr
3273         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3274 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3275         write(ifhi,*) x,PomIncXUnit(x,biniDf)
3276         write(*,*) 'Xu',i
3277       enddo
3278
3279       write(ifhi,'(a)')    '  endarray'
3280       write(ifhi,'(a)')    'closehisto plot 0'
3281
3282       end
3283
3284 c----------------------------------------------------------------------
3285       subroutine xParPomIncXI
3286 c----------------------------------------------------------------------
3287
3288       include 'epos.inc'
3289       include 'epos.incpar'
3290
3291       double precision x,xminr
3292       double precision PomIncXIExact,PomIncXIUnit
3293
3294       nptg=20                   !number of point for the graphs
3295       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3296
3297 c*********************red = PomIncXIExact *****************************
3298
3299       write(ifhi,'(a)')       'openhisto name PomIncXIExact'
3300       write(ifhi,'(a)')       'htyp lru'
3301       write(ifhi,'(a)')       'xmod log ymod log'
3302       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3303       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3304       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx(x)"'
3305       if (xpar8.eq.1.) then
3306         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3307       endif
3308       write(ifhi,'(a)')       'array 2'
3309
3310       do i=0,nptg
3311         x=xminr
3312         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3313 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3314         write(ifhi,*) x,PomIncXIExact(x)
3315 c.......write(*,*) 'XIe',i
3316       enddo
3317
3318       write(ifhi,'(a)')    '  endarray'
3319       write(ifhi,'(a)')    'closehisto plot 0-'
3320
3321 c***************************dot = PomIncXIUnit ************************
3322 c.....nptg=50     !number of point for the graphs
3323
3324       write(ifhi,'(a)')       'openhisto name PomIncXIUnit'
3325       write(ifhi,'(a)')       'htyp pfc'
3326       write(ifhi,'(a)')       'xmod log ymod log'
3327       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3328       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3329       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx(x)"'
3330       if (xpar8.eq.1.) then
3331         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3332       endif
3333       write(ifhi,'(a)')       'array 2'
3334
3335       do i=0,nptg
3336         x=xminr
3337         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3338 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3339         write(ifhi,*) x,PomIncXIUnit(x)
3340         write(*,*) 'XIu',i
3341       enddo
3342
3343       write(ifhi,'(a)')    '  endarray'
3344       write(ifhi,'(a)')    'closehisto plot 0'
3345
3346       end
3347
3348 c----------------------------------------------------------------------
3349       subroutine xParPomIncP
3350 c----------------------------------------------------------------------
3351
3352       include 'epos.inc'
3353       include 'epos.incpar'
3354
3355       double precision x,PomIncPUnit,xminr
3356       double precision PomIncPExact
3357
3358       nptg=30                  !number of point for the graphs
3359       biniDf=xpar1              !value of biniDf (impact parameter)
3360       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3361
3362 c*********************red = PomIncPExact *****************************
3363
3364       if (biniDf.lt.1.) then
3365         k=int(10.*biniDf)
3366         write(ifhi,'(a,I1)')  'openhisto name PomIncPExact-b0.',k
3367       else
3368         write(ifhi,'(a,f3.1)')'openhisto name PomIncPExact-b',biniDf
3369       endif
3370       write(ifhi,'(a)')       'htyp lru'
3371       write(ifhi,'(a)')       'xmod log ymod log'
3372       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3373       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3374       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+(x+,b)"'
3375       if (xpar8.eq.1.) then
3376         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3377         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3378       endif
3379       write(ifhi,'(a)')       'array 2'
3380
3381       do i=0,nptg
3382         x=xminr
3383         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3384         write(ifhi,*) x,PomIncPExact(x,biniDf)
3385       enddo
3386
3387       write(ifhi,'(a)')    '  endarray'
3388       write(ifhi,'(a)')    'closehisto plot 0-'
3389
3390
3391 c**************************dot = PomIncPUnit **************************
3392 c.....nptg=50     !number of point for the graphs
3393
3394       if (biniDf.lt.1.) then
3395         k=int(10.*biniDf)
3396         write(ifhi,'(a,I1)')  'openhisto name PomIncPUnit-b0.',k
3397       else
3398         write(ifhi,'(a,f3.1)')'openhisto name PomIncPUnit-b',biniDf
3399       endif
3400       write(ifhi,'(a)')       'htyp pfc'
3401       write(ifhi,'(a)')       'xmod log ymod log'
3402       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3403       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3404       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx+(x+,b)"'
3405       if (xpar8.eq.1.) then
3406         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3407         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3408       endif
3409       write(ifhi,'(a)')       'array 2'
3410
3411       do i=0,nptg
3412         x=xminr
3413         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3414         write(*,*) 'Pu',i
3415         write(ifhi,*) x,PomIncPUnit(x,biniDf)
3416       enddo
3417
3418       write(ifhi,'(a)')    '  endarray'
3419       write(ifhi,'(a)')    'closehisto plot 0'
3420
3421       end
3422
3423 c----------------------------------------------------------------------
3424       subroutine xParPomIncPI
3425 c----------------------------------------------------------------------
3426
3427       include 'epos.inc'
3428       include 'epos.incpar'
3429
3430       double precision x,xminr
3431       double precision PomIncPIExact,PomIncPIUnit
3432
3433       nptg=50                  !number of point for the graphs
3434       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3435
3436 c*********************red = PomIncPIExact *****************************
3437 c.....nptg=100     !number of point for the graphs
3438
3439       write(ifhi,'(a)')       'openhisto name PomIncPIExact'
3440       write(ifhi,'(a)')       'htyp lru'
3441       write(ifhi,'(a)')       'xmod log ymod log'
3442       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3443       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3444       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx+(x+)"'
3445       if (xpar8.eq.1.) then
3446         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3447       endif
3448       write(ifhi,'(a)')       'array 2'
3449
3450       do i=0,nptg
3451         x=xminr
3452         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3453         write(ifhi,*) x,PomIncPIExact(x)
3454       enddo
3455
3456       write(ifhi,'(a)')    '  endarray'
3457       write(ifhi,'(a)')    'closehisto plot 0-'
3458
3459 c***************************dot = PomIncPIUnit ************************
3460 c.....nptg=10     !number of point for the graphs
3461
3462       write(ifhi,'(a)')       'openhisto name PomIncPIUnit'
3463       write(ifhi,'(a)')       'htyp pfc'
3464       write(ifhi,'(a)')       'xmod log ymod log'
3465       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3466       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3467       write(ifhi,'(a)')    'text 0 0 "yaxis n?Pom!/dx+(x+)"'
3468       if (xpar8.eq.1.) then
3469         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3470       endif
3471       write(ifhi,'(a)')       'array 2'
3472
3473       do i=0,nptg
3474         x=xminr
3475         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3476         write(ifhi,*) x,PomIncPIUnit(x)
3477         write(*,*) 'PIu',i
3478       enddo
3479
3480       write(ifhi,'(a)')    '  endarray'
3481       write(ifhi,'(a)')    'closehisto plot 0'
3482
3483       end
3484
3485 c----------------------------------------------------------------------
3486       subroutine xParPomIncM
3487 c----------------------------------------------------------------------
3488
3489       include 'epos.inc'
3490       include 'epos.incpar'
3491
3492       double precision x,xminr,PomIncMUnit,PomIncMExact
3493
3494       nptg=50                  !number of point for the graphs
3495       biniDf=xpar1              !value of biniDf (impact parameter)
3496       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3497
3498 c**********************red = PomIncMExact *****************************
3499
3500       if (biniDf.lt.1.) then
3501         k=int(10.*biniDf)
3502         write(ifhi,'(a,I1)')  'openhisto name PomIncMExact-b0.',k
3503       else
3504         write(ifhi,'(a,f3.1)')'openhisto name PomIncMExact-b',biniDf
3505       endif
3506       write(ifhi,'(a)')       'htyp lru'
3507       write(ifhi,'(a)')       'xmod log ymod log'
3508       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3509       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3510       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx-(x-,b)"'
3511       if (xpar8.eq.1.) then
3512         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3513         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3514       endif
3515       write(ifhi,'(a)')       'array 2'
3516
3517       do i=0,nptg
3518         x=xminr
3519         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3520         write(ifhi,*) x,PomIncMExact(x,biniDf)
3521       enddo
3522
3523       write(ifhi,'(a)')    '  endarray'
3524       write(ifhi,'(a)')    'closehisto plot 0-'
3525
3526
3527 c**************************dot = PomIncMUnit **************************
3528 c.....nptg=100     !number of point for the graphs
3529
3530       if (biniDf.lt.1.) then
3531         k=int(10.*biniDf)
3532         write(ifhi,'(a,I1)')  'openhisto name PomIncMUnit-b0.',k
3533       else
3534         write(ifhi,'(a,f3.1)')'openhisto name PomIncMUnit-b',biniDf
3535       endif
3536       write(ifhi,'(a)')       'htyp pfc'
3537       write(ifhi,'(a)')       'xmod log ymod log'
3538       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3539       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3540       write(ifhi,'(a)')    'text 0 0 "yaxis dn?Pom!/dx-(x-,b)"'
3541       if (xpar8.eq.1.) then
3542         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3543         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
3544       endif
3545       write(ifhi,'(a)')       'array 2'
3546
3547       do i=0,nptg
3548         x=xminr
3549         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3550         write(ifhi,*) x,PomIncMUnit(x,biniDf)
3551         write(*,*) 'Mu',i
3552       enddo
3553
3554       write(ifhi,'(a)')    '  endarray'
3555       write(ifhi,'(a)')    'closehisto plot 0'
3556
3557       end
3558
3559 c----------------------------------------------------------------------
3560       subroutine xParPomIncMI
3561 c----------------------------------------------------------------------
3562
3563       include 'epos.inc'
3564       include 'epos.incpar'
3565
3566       double precision x,xminr
3567       double precision PomIncMIExact,PomIncMIUnit
3568
3569       nptg=30                 !number of point for the graphs
3570       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3571
3572 c*********************red = PomIncMIExact *****************************
3573 c.....nptg=100     !number of point for the graphs
3574
3575       write(ifhi,'(a)')       'openhisto name PomIncMIExact'
3576       write(ifhi,'(a)')       'htyp lru'
3577       write(ifhi,'(a)')       'xmod log ymod log'
3578       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3579       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3580       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx-(x-)"'
3581       if (xpar8.eq.1.) then
3582         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3583       endif
3584       write(ifhi,'(a)')       'array 2'
3585
3586       do i=0,nptg
3587         x=xminr
3588         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3589         write(ifhi,*) x,PomIncMIExact(x)
3590       enddo
3591
3592       write(ifhi,'(a)')    '  endarray'
3593       write(ifhi,'(a)')    'closehisto plot 0-'
3594
3595 c***************************dot = PomIncMIUnit ************************
3596 c.....nptg=100     !number of point for the graphs
3597
3598       write(ifhi,'(a)')       'openhisto name PomIncMIUnit'
3599       write(ifhi,'(a)')       'htyp pfc'
3600       write(ifhi,'(a)')       'xmod log ymod log'
3601       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3602       write(ifhi,'(a)')    'text 0 0 "xaxis x-"'
3603       write(ifhi,'(a)')    'text 0 0 "yaxis d[s]?Pom!/dx-(x-)"'
3604       if (xpar8.eq.1.) then
3605         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3606       endif
3607       write(ifhi,'(a)')       'array 2'
3608
3609       do i=0,nptg
3610         x=xminr
3611         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3612         write(ifhi,*) x,PomIncMIUnit(x)
3613         write(*,*) 'MIu',i
3614       enddo
3615
3616       write(ifhi,'(a)')    '  endarray'
3617       write(ifhi,'(a)')    'closehisto plot 0'
3618
3619       end
3620
3621 c----------------------------------------------------------------------
3622       subroutine xParPomIncJ
3623 c----------------------------------------------------------------------
3624
3625       include 'epos.inc'
3626       include 'epos.incpar'
3627       common/geom/rmproj,rmtarg,bmax,bkmx
3628       parameter(nbib=48)
3629
3630       double precision PomIncJExact,PomIncJUnit
3631
3632
3633       b1=0
3634       b2=bkmx*1.2
3635       db=(b2-b1)/nbib
3636
3637 c*************************red = PomIncJExact **************************
3638
3639       write(ifhi,'(a)')       'openhisto name PomIncJExact'
3640       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
3641       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
3642       write(ifhi,'(a)')    'text 0 0 "yaxis  n?Pom!(b)"'
3643       if (xpar8.eq.1.) then
3644         write(ifhi,'(a,e7.2,a)')'text 0.5 0.8 "s=',engy**2,' GeV^2!"'
3645       endif
3646       write(ifhi,'(a)')       'array 2'
3647       do k=1,nbib
3648         b=b1+(k-0.5)*db
3649         write(ifhi,*)b,PomIncJExact(b)
3650       enddo
3651       write(ifhi,'(a)')    '  endarray'
3652       write(ifhi,'(a)')    'closehisto plot 0-'
3653
3654 c****************************dot = PomIncJUnit ***********************
3655
3656       write(ifhi,'(a)')       'openhisto name PomIncJUnit'
3657       write(ifhi,'(a)')       'htyp pfc xmod lin ymod lin'
3658       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
3659       write(ifhi,'(a)')    'text 0 0 "yaxis  n?Pom!(b)"'
3660       if (xpar8.eq.1.) then
3661         write(ifhi,'(a,e7.2,a)')'text 0.5 0.8 "s=',engy**2,' GeV^2!"'
3662       endif
3663       write(ifhi,'(a)')       'array 2'
3664       do k=1,nbib
3665         b=b1+(k-0.5)*db
3666         write(ifhi,*)b,PomIncJUnit(b)
3667         write(*,*) 'Ju',k
3668       enddo
3669       write(ifhi,'(a)')    '  endarray'
3670       write(ifhi,'(a)')    'closehisto plot 0'
3671
3672
3673       end
3674
3675 c----------------------------------------------------------------------
3676       subroutine xParPhi1
3677 c----------------------------------------------------------------------
3678
3679       include 'epos.inc'
3680       include 'epos.incems'
3681       include 'epos.incsem'
3682       include 'epos.incpar'
3683       double precision x,xminr,y
3684       double precision PhiExpo
3685       double precision PhiExact
3686
3687       nptg=30                  !number of point for the graphs
3688       biniDf=xpar2              !value of biniDf (impact parameter)
3689       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
3690
3691 c************************* red = PhiMExact ***************************
3692
3693       write(ifhi,'(a)')       'openhisto name Phi1Exact'
3694       write(ifhi,'(a)')       'htyp lru'
3695       write(ifhi,'(a)')      'xmod log ymod lin'
3696       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3697       write(ifhi,'(a)')       'yrange auto auto'
3698       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3699       write(ifhi,'(a)') 'txt "yaxis  [F](x)/x^[a]!"'
3700       write(ifhi,'(a,i4,a,f4.1,a)')
3701      * 'txt  "title E=',nint(engy),' b=',biniDf,'"'
3702       write(ifhi,'(a)')       'array 2'
3703
3704       do i=0,nptg
3705         x=xminr
3706         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3707 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3708         y=0.d0
3709         if(engy**2..lt.5.e06)
3710      &  y=PhiExact(.5,dsqrt(x),dsqrt(x),engy**2,biniDf)
3711      &       *dsqrt(x)**dble(-alplea(iclpro))
3712      &       *dsqrt(x)**dble(-alplea(icltar))
3713         write(ifhi,*)x,y
3714       enddo
3715
3716       write(ifhi,'(a)')    '  endarray'
3717       write(ifhi,'(a)')    'closehisto plot 0-'
3718
3719 c********************** blue = PhiMExpo ******************************
3720
3721       write(ifhi,'(a)')       'openhisto name Phi1Expo'
3722       write(ifhi,'(a)')       'htyp lba'
3723       write(ifhi,'(a)')       'xmod log ymod lin'
3724       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3725       write(ifhi,'(a,2e11.3)')'yrange auto auto'
3726        write(ifhi,'(a)')       'array 2'
3727
3728       do i=0,nptg
3729         x=xminr
3730         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3731 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3732         y=PhiExpo(.5,dsqrt(x),dsqrt(x),engy**2,biniDf)
3733      &       *dsqrt(x)**dble(-alplea(iclpro))
3734      &       *dsqrt(x)**dble(-alplea(icltar))
3735         write(ifhi,*) x,y
3736       enddo
3737
3738       write(ifhi,'(a)')    '  endarray'
3739       write(ifhi,'(a)')    'closehisto plot 0'
3740
3741       end
3742
3743
3744 cc----------------------------------------------------------------------
3745 c      subroutine xParPhi2
3746 cc----------------------------------------------------------------------
3747 c
3748 c      include 'epos.inc'
3749 c      include 'epos.incems'
3750 c      include 'epos.incsem'
3751 c      include 'epos.incpar'
3752 c      double precision x,xminr,xm,y,u(0:100),v(0:100),w(0:100)!,z
3753 c      double precision PhiExpo,omGam,PhiExpoK
3754 c      double precision PhiExact
3755 c
3756 c      nptg=30                  !number of point for the graphs
3757 c      biniDf=xpar2              !value of biniDf (impact parameter)
3758 c      xminr=1.d-3 !/dble(engy**2)  !value of xminr for plotting the function
3759 c      xm=dble(xpar4)
3760 c
3761 cc************************* yellow = PhiMExact ***************************
3762 c
3763 c      write(ifhi,'(a)')       'openhisto name Phi1Exact'
3764 c      write(ifhi,'(a)')       'htyp lru'
3765 c      write(ifhi,'(a)')      'xmod lin ymod lin'
3766 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3767 c      write(ifhi,'(a)')       'yrange auto auto'
3768 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3769 c      write(ifhi,'(a)') 'txt "yaxis  [F](x)/x^[a]!"'
3770 c      write(ifhi,'(a,i4,a,f4.1,a)')
3771 c     * 'txt  "title E=',nint(engy),' b=',biniDf,'"'
3772 c      write(ifhi,'(a)')       'array 2'
3773 c
3774 c      do i=0,nptg
3775 c       ! x=xminr
3776 c       ! if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3777 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3778 c        y=0.d0
3779 c        if(engy**2..lt.5.e06)
3780 c     &  y=PhiExact(1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
3781 c    ! &       *dsqrt(x)**dble(-alplea(iclpro))
3782 c    ! &       *dsqrt(x)**dble(-alplea(icltar))
3783 c        write(ifhi,*)x,y
3784 c      enddo
3785 c
3786 c      write(ifhi,'(a)')    '  endarray'
3787 c      write(ifhi,'(a)')    'closehisto plot 0-'
3788 c
3789 cc********************** blue = PhiMExpo ******************************
3790 c
3791 c      write(ifhi,'(a)')       'openhisto name Phi1Expo'
3792 c      write(ifhi,'(a)')       'htyp lbu'
3793 c      write(ifhi,'(a)')       'xmod lin ymod lin'
3794 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3795 c      write(ifhi,'(a,2e11.3)')'yrange auto auto'
3796 c       write(ifhi,'(a)')       'array 2'
3797 c
3798 c      do i=0,nptg
3799 c       ! x=xminr
3800 c       ! if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3801 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3802 c        y=PhiExpo(1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
3803 c  !   &       *dsqrt(x)**dble(-alplea(iclpro))
3804 c  !   &       *dsqrt(x)**dble(-alplea(icltar))
3805 c        write(ifhi,*) x,y
3806 c      enddo
3807 c
3808 c      write(ifhi,'(a)')    '  endarray'
3809 c      write(ifhi,'(a)')    'closehisto plot 0-'
3810 c
3811 cc**********************************************************************
3812 cc**********************************************************************
3813 c      do k=1,koll
3814 c        bk(k)=biniDf
3815 c      enddo
3816 c      call GfunPark(0)
3817 c      call integom1(0)
3818 c
3819 cc********************* points = PhiExpoK*********************************
3820 c
3821 c      if (biniDf.lt.1.) then
3822 c        k=int(10.*biniDf)
3823 c        write(ifhi,'(a,I1)')  'openhisto name PhiExpok-b0.',k
3824 c      else
3825 c        write(ifhi,'(a,f3.1)')  'openhisto name PhiExpok-b',biniDf
3826 c      endif
3827 c      write(ifhi,'(a)')       'htyp pfc'
3828 c      write(ifhi,'(a)')       'xmod log ymod lin'
3829 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3830 c      write(ifhi,'(a,2e11.3)')'yrange auto auto'
3831 c      write(ifhi,'(a)')    'text 0 0 "xaxis x"'
3832 c      write(ifhi,'(a)') 'text 0 0.1 "yaxis  [F](x+,x-)/x^[a]?remn!!"'
3833 c      if (xpar8.eq.1.) then
3834 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3835 c        write(ifhi,'(a,f5.2,a)')'text 0.5 0.9 "b=',biniDf,' fm"'
3836 c      endif
3837 c      write(ifhi,'(a)')       'array 2'
3838 c
3839 c      do i=0,nptg
3840 c        x=xminr
3841 c        if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3842 c        y=PhiExpoK(1,dsqrt(x),dsqrt(x))
3843 c        write(ifhi,*) x,y
3844 c      enddo
3845 c
3846 c      write(ifhi,'(a)')    '  endarray'
3847 c      write(ifhi,'(a)')    'closehisto plot 0'
3848 c
3849 cc************************* red = PhiMExact*omGam ************************
3850 c
3851 c      write(ifhi,'(a)')       'openhisto name GPhiExact'
3852 c      write(ifhi,'(a)')       'htyp lru'
3853 c      if(xpar5.eq.0.)then
3854 c        write(ifhi,'(a)')       'xmod lin ymod lin'
3855 c      else
3856 c        write(ifhi,'(a)')       'xmod log ymod log'
3857 c      endif
3858 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3859 c      write(ifhi,'(a)')       'yrange auto auto'
3860 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3861 c      write(ifhi,'(a,a)') 'text 0 0.1 "yaxis  G(x+,x-)*[F]'
3862 c     *,'(x+,x-)/x^[a]?remn!!"'
3863 c      if (xpar8.eq.1.) then
3864 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.2 "s=',engy**2,' GeV^2!"'
3865 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "b=',biniDf,' fm"'
3866 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "x-=',xm,'"'
3867 c      endif
3868 c      write(ifhi,'(a)')       'array 2'
3869 c
3870 c      do i=0,nptg
3871 c        x=xminr
3872 c        if(xpar5.ne.0.)then
3873 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3874 c        else
3875 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3876 c        endif
3877 cc        z=1.d0-dsqrt(x)
3878 c        v(i)=0.d0
3879 c        if(engy**2..lt.5.e06)
3880 c     *  v(i)=PhiExact(1.,1.d0-x,1.d0-xm,engy**2,biniDf)
3881 cc     *  v(i)=PhiExact(1.,z,z,engy**2,biniDf)
3882 c        u(i)=omGam(x,xm,biniDf)
3883 cc        u(i)=omGam(dsqrt(x),dsqrt(x),biniDf)
3884 c        y=u(i)*v(i)
3885 c        if(xpar5.ne.0.)y=dabs(y)
3886 c        write(ifhi,*)x,y
3887 c      enddo
3888 c
3889 c      write(ifhi,'(a)')    '  endarray'
3890 c      write(ifhi,'(a)')    'closehisto plot 0-'
3891 c
3892 cc************************* red = PhiMExpo*omGam ************************
3893 c
3894 c      write(ifhi,'(a)')       'openhisto name GPhiExpo'
3895 c      write(ifhi,'(a)')       'htyp lba'
3896 c      if(xpar5.eq.0.)then
3897 c        write(ifhi,'(a)')       'xmod lin ymod lin'
3898 c      else
3899 c        write(ifhi,'(a)')       'xmod log ymod log'
3900 c      endif
3901 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
3902 c      write(ifhi,'(a)')       'yrange auto auto'
3903 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
3904 c      write(ifhi,'(a,a)')
3905 c     * 'text 0 0.1 "yaxis  G(x+,x-)*[F]?'
3906 c     * ,'(1-x+,1-x-)/x^[a]?remn!!"'
3907 c      if (xpar8.eq.1.) then
3908 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
3909 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.8 "b=',biniDf,' fm"'
3910 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.7 "x-=',xm,'"'
3911 c      endif
3912 c      write(ifhi,'(a)')       'array 2'
3913 c
3914 c      do i=0,nptg
3915 c        x=xminr
3916 c        if(xpar5.ne.0.)then
3917 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3918 c        else
3919 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3920 c        endif
3921 cc        z=1.d0-dsqrt(x)
3922 cc        w(i)=PhiExpo(1.,z,z,engy**2,biniDf)
3923 c        w(i)=PhiExpo(1.,1.d0-x,1.d0-xm,engy**2,biniDf)
3924 c        y=u(i)*w(i)
3925 c        if(xpar5.ne.0.)y=dabs(y)
3926 c        write(ifhi,*)x,y
3927 c      enddo
3928 c
3929 c      write(ifhi,'(a)')    '  endarray'
3930 c      if(xpar5.ne.0.)then
3931 c      write(ifhi,'(a)')    'closehisto plot 0-'
3932 c
3933 cc************************* green = omGam ************************
3934 c
3935 c      write(ifhi,'(a)')       'openhisto name GM'
3936 c      write(ifhi,'(a)')       'htyp lgo'
3937 c      write(ifhi,'(a)')       'array 2'
3938 c
3939 c      do i=0,nptg
3940 c        x=xminr
3941 c        if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3942 c        write(ifhi,*)x,dabs(u(i))
3943 c      enddo
3944 c
3945 c      write(ifhi,'(a)')    '  endarray'
3946 c      write(ifhi,'(a)')    'closehisto plot 0-'
3947 c
3948 cc************************* circle = PhiMExact  ************************
3949 c
3950 c      write(ifhi,'(a)')       'openhisto name PhiExact'
3951 c      write(ifhi,'(a)')       'htyp poc'
3952 c      write(ifhi,'(a)')       'array 2'
3953 c
3954 c      do i=0,nptg
3955 c        x=xminr
3956 c        if(xpar5.ne.0.)then
3957 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3958 c        else
3959 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3960 c        endif
3961 c        y=v(i)
3962 c        if(xpar5.ne.0.)y=dabs(y)
3963 c        write(ifhi,*)x,y
3964 c      enddo
3965 c
3966 c      write(ifhi,'(a)')    '  endarray'
3967 c      write(ifhi,'(a)')    'closehisto plot 0-'
3968 c
3969 cc************************* triangle = PhiMExpo ************************
3970 c
3971 c      write(ifhi,'(a)')       'openhisto name PhiExpo'
3972 c      write(ifhi,'(a)')       'htyp pot'
3973 c      write(ifhi,'(a)')       'array 2'
3974 c
3975 c      do i=0,nptg
3976 c        x=xminr
3977 c        if(xpar5.ne.0.)then
3978 c          if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
3979 c        else
3980 c          x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
3981 c        endif
3982 c        y=w(i)
3983 c        if(xpar5.ne.0.)y=dabs(y)
3984 c        write(ifhi,*)x,y
3985 c      enddo
3986 c
3987 c      write(ifhi,'(a)')    '  endarray'
3988 c      endif
3989 c      write(ifhi,'(a)')    'closehisto plot 0'
3990 c
3991 c      end
3992 c
3993
3994 c----------------------------------------------------------------------
3995       subroutine xParPhi
3996 c----------------------------------------------------------------------
3997
3998       include 'epos.inc'
3999       include 'epos.incems'
4000       include 'epos.incsem'
4001       include 'epos.incpar'
4002       double precision x,xminr,y,z(0:200)!,Zn,Znorm
4003       double precision PhiExpo,PhiExact,PhiExpoK
4004
4005
4006       nptg=10                  !number of point for the graphs
4007       biniDf=xpar2              !value of biniDf (impact parameter)
4008       xminr=max(1.d-6,1.d0/dble(engy**2))  !value of xminr for plotting the function
4009
4010 c********************** full-red = PhiExact ***************************
4011
4012       if (biniDf.lt.1.) then
4013         k=int(10.*biniDf)
4014         write(ifhi,'(a,I1)')  'openhisto name PhiExact-b0.',k
4015       else
4016         write(ifhi,'(a,f3.1)')  'openhisto name PhiExact-b',biniDf
4017       endif
4018       write(ifhi,'(a)')       'htyp lru'
4019       write(ifhi,'(a)')       'xmod lin ymod lin'
4020       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4021       write(ifhi,'(a)')       'yrange auto auto'
4022       write(ifhi,'(a)')    'text 0 0 "xaxis x"'
4023       write(ifhi,'(a)') 'text 0 0 "yaxis  [F](x)/x^[a]"'
4024       write(ifhi,'(a,i4,a,f4.1,a)')
4025      * 'txt  "title E=',nint(engy),' b=',biniDf,'"'
4026       write(ifhi,'(a)')       'array 2'
4027
4028       do i=0,nptg
4029         !x=xminr
4030         !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4031         x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4032         y=PhiExact(1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
4033      &       *dsqrt(x)**dble(-alplea(iclpro))
4034      &       *dsqrt(x)**dble(-alplea(icltar))
4035         write(ifhi,*)x,y
4036       enddo
4037
4038       write(ifhi,'(a)')    '  endarray'
4039       write(ifhi,'(a)')    'closehisto plot 0-'
4040
4041 c******************** blue = PhiExpo ***************************
4042
4043       if (biniDf.lt.1.) then
4044         k=int(10.*biniDf)
4045         write(ifhi,'(a,I1)')  'openhisto name PhiExpo-b0.',k
4046       else
4047         write(ifhi,'(a,f3.1)')  'openhisto name PhiExpo-b',biniDf
4048       endif
4049       write(ifhi,'(a)')       'htyp lbu'
4050       write(ifhi,'(a)')       'xmod lin ymod lin'
4051       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4052       write(ifhi,'(a,2e11.3)')'yrange auto auto'
4053       write(ifhi,'(a)')       'array 2'
4054
4055       do i=0,nptg
4056         !x=xminr
4057         !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4058         x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4059         z(i)=PhiExpo(1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
4060      &       *dsqrt(x)**dble(-alplea(iclpro))
4061      &       *dsqrt(x)**dble(-alplea(icltar))
4062         write(ifhi,*) x,z(i)
4063       enddo
4064
4065       write(ifhi,'(a)')    '  endarray'
4066 c      write(ifhi,'(a)')    'closehisto plot 0-'
4067 c
4068 cc*********************yellow = PhiUnit*********************************
4069 c
4070 c      if (biniDf.lt.1.) then
4071 c        k=int(10.*biniDf)
4072 c        write(ifhi,'(a,I1)')  'openhisto name PhiUnit-b0.',k
4073 c      else
4074 c        write(ifhi,'(a,f3.1)')  'openhisto name PhiUnit-b',biniDf
4075 c      endif
4076 c      write(ifhi,'(a)')       'htyp lyu'
4077 c      write(ifhi,'(a)')       'xmod lin ymod lin'
4078 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4079 c      write(ifhi,'(a,2e11.3)')'yrange auto auto'
4080 c      write(ifhi,'(a)')       'array 2'
4081 c
4082 c      Zn=Znorm(engy**2,biniDf)
4083 c
4084 c      do i=0,nptg
4085 c        !x=xminr
4086 c        !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4087 c        x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4088 c        write(ifhi,*) x,z(i)/Zn
4089 c      enddo
4090 c
4091 c      write(ifhi,'(a)')    '  endarray'
4092
4093       if(koll.ge.1)then
4094
4095       write(ifhi,'(a)')    'closehisto plot 0-'
4096
4097 c**********************************************************************
4098 c**********************************************************************
4099       do k=1,koll
4100         bk(k)=biniDf
4101       enddo
4102       call GfunPark(0)
4103       call integom1(0)
4104
4105
4106 c*********************green = PhiExpoK*********************************
4107
4108       if (biniDf.lt.1.) then
4109         k=int(10.*biniDf)
4110         write(ifhi,'(a,I1)')  'openhisto name PhiExpok-b0.',k
4111       else
4112         write(ifhi,'(a,f3.1)')  'openhisto name PhiExpok-b',biniDf
4113       endif
4114       write(ifhi,'(a)')       'htyp lga'
4115       write(ifhi,'(a)')       'xmod lin ymod lin'
4116       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4117       write(ifhi,'(a,2e11.3)')'yrange auto auto'
4118       write(ifhi,'(a)')       'array 2'
4119
4120       do i=0,nptg
4121         !x=xminr
4122         !if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4123         x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4124         z(i)=PhiExpoK(1,dsqrt(x),dsqrt(x))
4125         write(ifhi,*) x,z(i)
4126       enddo
4127
4128       write(ifhi,'(a)')    '  endarray'
4129       endif
4130       write(ifhi,'(a)')    'closehisto plot 0'
4131
4132 c
4133       end
4134
4135 c----------------------------------------------------------------------
4136       subroutine xParH
4137 c----------------------------------------------------------------------
4138
4139       include 'epos.inc'
4140       include 'epos.incpar'
4141       parameter(idxD2=8)
4142       double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni
4143       common/DGamUni/GbetUni(  idxD0:idxD2),HbetUni(  idxD0:idxD2),
4144      &               GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2),
4145      &               HalpUni(idxD0:idxD2)
4146       double precision x,xminr,y,xm,utgam2
4147       double precision Hrst
4148
4149       nptg=20                  !number of point for the graphs
4150       biniDf=xpar2              !value of biniDf (impact parameter)
4151       xm=dble(xpar4)            !value of xminus
4152 c.....xminr=0.d0   !value of xminr for plotting the function
4153       xminr=1.d0/dble(engy**2)  !value of xminr for plotting the function
4154
4155       imax=idxD1
4156       if(iomega.eq.2)imax=1
4157       do i=idxDmin,imax
4158         call GfunPar(1,i,biniDf,smaxDf,alpx,betx,betpx,epsp,epst,epss
4159      &               ,gamv)
4160       enddo
4161       call GfomPar(biniDf,smaxDf)
4162       imax0=idxD1
4163       if(iomega.eq.2)imax0=1
4164       imax1=idxD2
4165       if(iomega.eq.2)imax1=imax1-1
4166       do i=idxDmin,imax0
4167         GbetUni(i)=utgam2(betUni(i,1)+1.d0)
4168         GbetpUni(i)=utgam2(betpUni(i,1)+1.d0)
4169         HbetUni(i)=utgam2(GbetUni(i))
4170         HbetpUni(i)=utgam2(GbetpUni(i))
4171         HalpUni(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
4172       enddo
4173       do i=0,1
4174         HbetUni(imax0+1+i)=betUni(i,1)+1.d0+betfom
4175         HbetUni(imax0+3+i)=betUni(i,1)+1.d0
4176         HbetUni(imax0+5+i)=betUni(i,1)+1.d0+betfom
4177         HbetpUni(imax0+1+i)=betpUni(i,1)+1.d0
4178         HbetpUni(imax0+3+i)=betpUni(i,1)+1.d0+betfom
4179         HbetpUni(imax0+5+i)=betpUni(i,1)+1.d0+betfom
4180         GbetUni(imax0+1+i)=utgam2(HbetUni(imax0+1+i))
4181         GbetUni(imax0+3+i)=utgam2(HbetUni(imax0+3+i))
4182         GbetUni(imax0+5+i)=utgam2(HbetUni(imax0+5+i))
4183         GbetpUni(imax0+1+i)=utgam2(HbetpUni(imax0+1+i))
4184         GbetpUni(imax0+3+i)=utgam2(HbetpUni(imax0+3+i))
4185         GbetpUni(imax0+5+i)=utgam2(HbetpUni(imax0+5+i))
4186         HalpUni(imax0+1+i)=ztUni*alpUni(i,1)
4187         HalpUni(imax0+3+i)=zpUni*alpUni(i,1)
4188         HalpUni(imax0+5+i)=zpUni*ztUni*alpUni(i,1)
4189       enddo
4190
4191 c***********************  red = Hrst  *********************************
4192
4193       write(ifhi,'(a)')       'openhisto name Hrst'
4194       write(ifhi,'(a)')       'htyp lru'
4195       write(ifhi,'(a)')       'xmod log ymod lin'
4196       write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4197       write(ifhi,'(a)')       'yrange auto auto'
4198       write(ifhi,'(a)')
4199      *     'text 0 0 "yaxis  H?2!(x+,x-)"'
4200       write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
4201       write(ifhi,'(a)')
4202      *     'text 0 0 "yaxis  H?2!(x+,x-)"'
4203       if (xpar8.eq.1.) then
4204         write(ifhi,'(a,e7.2,a)')'text 0.1 0.2 "s=',engy**2,' GeV^2!"'
4205         write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "b=',biniDf,' fm"'
4206         write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "x-=',xm,'"'
4207       endif
4208       write(ifhi,'(a)')       'array 2'
4209
4210       do i=0,nptg
4211         x=xminr
4212         if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4213 c.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4214         y=Hrst(smaxDf,biniDf,dsqrt(x),dsqrt(x))
4215         write(ifhi,*)x,y
4216       enddo
4217
4218       write(ifhi,'(a)')    '  endarray'
4219       write(ifhi,'(a)')    'closehisto plot 0'
4220
4221       end
4222
4223
4224
4225
4226 cc----------------------------------------------------------------------
4227 c      subroutine xParHPhiIntnew
4228 cc----------------------------------------------------------------------
4229 c
4230 c      include 'epos.inc'
4231 c      include 'epos.incsem'
4232 c      include 'epos.incems'
4233 c      include 'epos.incpar'
4234 c      double precision x,xminr,xm,y
4235 c      double precision PhiExact,omGam
4236 cc      double precision PhiExpo
4237 c
4238 c      nptg=30                  !number of point for the graphs
4239 c      biniDf=xpar2              !value of biniDf (impact parameter)
4240 c      xm=dble(xpar4)            !value of xminus
4241 cc.....xminr=0.d0   !value of xminr for plotting the function
4242 c      xminr=1.d-3 !/dble(engy**2)  !value of xminr for plotting the function
4243 cc************************* black = PhiExact ***************************
4244 c
4245 c      write(ifhi,'(a)')       'openhisto name Phi1Exact'
4246 c      write(ifhi,'(a)')       'htyp lru'
4247 c      if(xpar5.eq.0.)then
4248 c        write(ifhi,'(a)')       'xmod lin ymod lin'
4249 c      else
4250 c        write(ifhi,'(a)')       'xmod log ymod lin'
4251 c      endif
4252 c      write(ifhi,'(a,2e11.3)')'xrange',xminr,xmaxDf
4253 c      write(ifhi,'(a)')       'yrange auto auto'
4254 c      write(ifhi,'(a)')    'text 0 0 "xaxis x+"'
4255 c      write(ifhi,'(a)')
4256 c     * 'text 0 0.1 "yaxis  [F]?(x+,x-)/x^[a]?remn!!"'
4257 c      if (xpar8.eq.1.) then
4258 c        write(ifhi,'(a,e7.2,a)')'text 0.1 0.2 "s=',engy**2,' GeV^2!"'
4259 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.1 "b=',biniDf,' fm"'
4260 c        write(ifhi,'(a,f5.2,a)')'text 0.1 0.3 "x-=',xm,'"'
4261 c      endif
4262 c      write(ifhi,'(a)')       'array 2'
4263 c
4264 c      do i=0,nptg
4265 c        x=xminr
4266 c        if (i.ne.0) x=x*(xmaxDf/xminr)**(dble(i)/dble(nptg))
4267 cc.......x=xminr+(xmaxDf-xminr)*(dble(i)/dble(nptg))
4268 c      y=PhiExact(1.,dsqrt(x),dsqrt(x),engy**2,biniDf)
4269 c     &       *omGam(dsqrt(x),dsqrt(x),biniDf)
4270 c        write(ifhi,*)x,y
4271 c      enddo
4272 c
4273 c      write(ifhi,'(a)')    '  endarray'
4274 c      write(ifhi,'(a)')    'closehisto plot 0'
4275 c
4276 c      end
4277 c
4278 c----------------------------------------------------------------------
4279       subroutine xParHPhiInt
4280 c----------------------------------------------------------------------
4281
4282       include 'epos.inc'
4283       include 'epos.incpar'
4284       double precision y,HPhiInt
4285       common/geom/rmproj,rmtarg,bmax,bkmx
4286       parameter(nbib=32)
4287
4288
4289 c************************ dotted = gauss integration ******************
4290
4291       b1=0
4292       b2=max(abs(bkmx),3.)*1.2
4293       db=(b2-b1)/nbib
4294
4295       write(ifhi,'(a)')       'openhisto name HPhiExpoInt'
4296       write(ifhi,'(a)')       'htyp pfc xmod lin ymod lin'
4297       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4298       write(ifhi,'(a)')    'text 0 0 "yaxis  Int(H[F]?pp!)(s,b)"'
4299       if (xpar8.eq.1.) then
4300         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4301       endif
4302       write(ifhi,'(a)')       'array 2'
4303       imax=idxD1
4304       if(iomega.eq.2)imax=1
4305       do k=1,nbib
4306         b=b1+(k-0.5)*db
4307         do i=idxDmin,imax
4308           call GfunPar(1,i,b,smaxDf,alpx,betx,betpx,epsp,epst,epss,gamv)
4309           call GfunPar(2,i,b,smaxDf,alpx,betx,betpx,epsp,epst,epss,gamv)
4310         enddo
4311         y=HPhiInt(smaxDf,b)
4312         write(ifhi,*)b,y
4313       enddo
4314       write(ifhi,'(a)')    '  endarray'
4315       write(ifhi,'(a)')    'closehisto plot 0'
4316
4317       end
4318
4319 c----------------------------------------------------------------------
4320       subroutine xParZ
4321 c----------------------------------------------------------------------
4322       include 'epos.inc'
4323       include 'epos.incpar'
4324       double precision Znorm,y
4325       common/geom/rmproj,rmtarg,bmax,bkmx
4326       parameter(nbib=12)
4327
4328       b1=0
4329       b2=max(abs(bkmx),3.)*1.2
4330       db=(b2-b1)/nbib
4331
4332 c************************full-red = Znorm *****************************
4333
4334       write(ifhi,'(a)')       'openhisto name Znorm'
4335       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4336       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4337       write(ifhi,'(a)')    'text 0 0 "yaxis  Z(s,b) "'
4338       if (xpar8.eq.1.) then
4339         write(ifhi,'(a,e7.2,a)')'text 0.1 0.1 "s=',engy**2,' GeV^2!"'
4340       endif
4341       write(ifhi,'(a)')       'array 2'
4342 c        y=Znorm(engy**2,xpar2)
4343       do k=1,nbib
4344         b=b1+(k-0.5)*db
4345         y=Znorm(smaxDf,b)
4346         write(ifhi,*)b,y
4347       enddo
4348       write(ifhi,'(a)')    '  endarray'
4349       write(ifhi,'(a)')    'closehisto plot 0-'
4350
4351
4352 c**********************************************************************
4353
4354       write(ifhi,'(a)')       'openhisto name un'
4355       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4356       write(ifhi,'(a)')       'array 2'
4357       do k=1,nbib
4358         b=b1+(k-0.5)*db
4359         y=1
4360         write(ifhi,*)b,y
4361       enddo
4362       write(ifhi,'(a)')    '  endarray'
4363       write(ifhi,'(a)')    'closehisto plot 0'
4364
4365
4366       end
4367
4368 c----------------------------------------------------------------------
4369       subroutine xParPro
4370 c----------------------------------------------------------------------
4371       include 'epos.inc'
4372       include 'epos.incsem'
4373       include 'epos.incpar'
4374       double precision PhiExact,PhiExpo,y,om1intb,om1intbc,om1intgc
4375      &,om1intbi!,PhiUnit
4376       common/geom/rmproj,rmtarg,bmax,bkmx
4377       parameter(nbib=12)
4378
4379       b1=0
4380       b2=max(abs(bkmx),3.)*1.2
4381       db=(b2-b1)/nbib
4382
4383 c********************* full-red = 1-PhiExact **************************
4384
4385       write(ifhi,'(a)')       'openhisto name 1-PhiExact'
4386       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4387       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4388       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4389       if (xpar8.eq.1.) then
4390         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4391       endif
4392       write(ifhi,'(a)')       'array 2'
4393       do k=1,nbib
4394         b=b1+(k-0.5)*db
4395         y=1.d0-PhiExact(1.,1.d0,1.d0,engy**2,b)
4396         write(ifhi,*)b,y
4397       enddo
4398       write(ifhi,'(a)')    '  endarray'
4399       write(ifhi,'(a)')    'closehisto plot 0-'
4400
4401 c************************** blue-dashed = 1-PhiExpo *******************
4402
4403       write(ifhi,'(a)')       'openhisto name 1-PhiExpo'
4404       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4405       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4406       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4407       if (xpar8.eq.1.) then
4408         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4409       endif
4410       write(ifhi,'(a)')       'array 2'
4411       do k=1,nbib
4412         b=b1+(k-0.5)*db
4413         y=1.d0-PhiExpo(1.,1.d0,1.d0,engy**2,b)
4414         write(ifhi,*)b,y
4415       enddo
4416       write(ifhi,'(a)')    '  endarray'
4417       write(ifhi,'(a)')    'closehisto plot 0-'
4418
4419
4420 c**********************************************************************
4421
4422       write(ifhi,'(a)')       'openhisto'
4423       write(ifhi,'(a)')       'htyp lga xmod lin ymod lin'
4424       write(ifhi,'(a)')       'array 2'
4425       do k=1,nbib
4426         b=b1+(k-0.5)*db
4427         y=1
4428         write(ifhi,*)b,y
4429       enddo
4430       write(ifhi,'(a)')    '  endarray'
4431       write(ifhi,'(a)')    'closehisto plot 0'
4432
4433 c****************************** red = om1intbc ********************
4434
4435       write(ifhi,'(a)')       'openhisto name om1intbc'
4436       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4437       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4438       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1bc!(b) "'
4439       if (xpar8.eq.1.) then
4440         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4441       endif
4442       write(ifhi,'(a)')       'array 2'
4443
4444       do k=1,nbib
4445         b=b1+(k-0.5)*db
4446         y=om1intbc(b)
4447         write(ifhi,*)b,y
4448       enddo
4449       write(ifhi,'(a)')    '  endarray'
4450       write(ifhi,'(a)')    'closehisto plot 0-'
4451
4452
4453 c************************* blue dashed =  om1intb ********************
4454
4455       write(ifhi,'(a)')       'openhisto name om1intb'
4456       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4457       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4458       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1b!(b) "'
4459       if (xpar8.eq.1.) then
4460         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4461       endif
4462       write(ifhi,'(a)')       'array 2'
4463
4464       do k=1,nbib
4465         b=b1+(k-0.5)*db
4466         y=om1intb(b)
4467         write(ifhi,*)b,y
4468       enddo
4469       write(ifhi,'(a)')    '  endarray'
4470       write(ifhi,'(a)')    'closehisto plot 0-'
4471
4472 c****************************** green dot = om1intgc ********************
4473
4474       write(ifhi,'(a)')       'openhisto name om1intgc'
4475       write(ifhi,'(a)')       'htyp lgo xmod lin ymod lin'
4476       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4477       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1gc!(b) "'
4478       if (xpar8.eq.1.) then
4479         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4480       endif
4481       write(ifhi,'(a)')       'array 2'
4482
4483       do k=1,nbib
4484         b=b1+(k-0.5)*db
4485         y=om1intgc(b)
4486         write(ifhi,*)b,y
4487       enddo
4488       write(ifhi,'(a)')    '  endarray'
4489       write(ifhi,'(a)')    'closehisto plot 0'
4490
4491 c****************************** red = om1intbi(0) ********************
4492
4493       write(ifhi,'(a)')       'openhisto name om1intbc'
4494       write(ifhi,'(a)')       'htyp lru xmod lin ymod log'
4495       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4496       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1bc!(b) "'
4497       if (xpar8.eq.1.) then
4498         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4499       endif
4500       write(ifhi,'(a)')       'array 2'
4501
4502       do k=1,nbib
4503         b=b1+(k-0.5)*db
4504         y=om1intbi(b,0)
4505         write(ifhi,*)b,y
4506       enddo
4507       write(ifhi,'(a)')    '  endarray'
4508       write(ifhi,'(a)')    'closehisto plot 0-'
4509
4510
4511 c************************* blue dashed =  om1intbi(1) ********************
4512
4513       write(ifhi,'(a)')       'openhisto name om1intb'
4514       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4515       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4516       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1b!(b) "'
4517       if (xpar8.eq.1.) then
4518         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4519       endif
4520       write(ifhi,'(a)')       'array 2'
4521
4522       do k=1,nbib
4523         b=b1+(k-0.5)*db
4524         y=om1intbi(b,1)
4525         write(ifhi,*)b,y
4526       enddo
4527       write(ifhi,'(a)')    '  endarray'
4528       write(ifhi,'(a)')    'closehisto plot 0-'
4529
4530 c****************************** green dot = om1intbi(2) ********************
4531
4532       write(ifhi,'(a)')       'openhisto name om1intgc'
4533       write(ifhi,'(a)')       'htyp lgo xmod lin ymod lin'
4534       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4535       write(ifhi,'(a)')    'text 0 0 "yaxis  [w]?1gc!(b) "'
4536       if (xpar8.eq.1.) then
4537         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4538       endif
4539       write(ifhi,'(a)')       'array 2'
4540
4541       do k=1,nbib
4542         b=b1+(k-0.5)*db
4543         y=om1intbi(b,2)
4544         write(ifhi,*)b,y
4545       enddo
4546       write(ifhi,'(a)')    '  endarray'
4547       write(ifhi,'(a)')    'closehisto plot 0'
4548
4549
4550
4551       end
4552
4553 c----------------------------------------------------------------------
4554       subroutine xParPro1
4555 c----------------------------------------------------------------------
4556       include 'epos.inc'
4557       include 'epos.incsem'
4558       include 'epos.incpar'
4559       double precision PhiExact,PhiExpo,y
4560       common/geom/rmproj,rmtarg,bmax,bkmx
4561       parameter(nbib=12)
4562
4563       b1=0
4564       b2=max(abs(bkmx),3.)*1.2
4565       db=(b2-b1)/nbib
4566
4567 c********************* full-red = 1-PhiExact **************************
4568
4569       write(ifhi,'(a)')       'openhisto name 1-PhiExact'
4570       write(ifhi,'(a)')       'htyp lru xmod lin ymod lin'
4571       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4572       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
4573       if (xpar8.eq.1.) then
4574         write(ifhi,'(a,e7.2,a)')'text 0.5 0.9 "s=',engy**2,' GeV^2!"'
4575       endif
4576       write(ifhi,'(a)')       'array 2'
4577       do k=1,nbib
4578         b=b1+(k-0.5)*db
4579         y=1.d0-PhiExact(.5,1.d0,1.d0,engy**2,b)
4580         write(ifhi,*)b,y
4581       enddo
4582       write(ifhi,'(a)')    '  endarray'
4583       write(ifhi,'(a)')    'closehisto plot 0-'
4584
4585 c************************** blue-dashed = 1-PhiExpo *******************
4586
4587       write(ifhi,'(a)')       'openhisto name 1-PhiExpo'
4588       write(ifhi,'(a)')       'htyp lba xmod lin ymod lin'
4589       write(ifhi,'(a)')    'text 0 0 "xaxis  impact parameter b (fm)"'
4590       write(ifhi,'(a)')    'text 0 0 "yaxis  1-[F]?pp!(1,1) "'
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       do k=1,nbib
4596         b=b1+(k-0.5)*db
4597         y=1.d0-PhiExpo(.5,1.d0,1.d0,engy**2,b)
4598         write(ifhi,*)b,y
4599       enddo
4600       write(ifhi,'(a)')    '  endarray'
4601       write(ifhi,'(a)')    'closehisto plot 0'
4602
4603
4604
4605
4606       end
4607
4608 c----------------------------------------------------------------------
4609       subroutine xParGam
4610 c----------------------------------------------------------------------
4611
4612       include 'epos.inc'
4613       include 'epos.incsem'
4614       include 'epos.incpar'
4615       dimension bet(idxD0:idxD1)
4616       double precision utgam2,xgammag2!,xrem
4617       dimension ip(idxD0:idxD1),imax(idxD0:idxD1)
4618
4619       nptg=50                  !number of point for the graphs
4620       b=xpar2
4621       gamp=xpar6
4622       zmax=6.
4623 c      xrem=dble(xpar4)
4624       if(idxD0.ne.0.or.idxD1.ne.2) stop "Check xPargam"
4625
4626       do i=idxD0,idxD1
4627         imax(i)=4
4628         bet(i)=0
4629       enddo
4630       nmax=idxD1
4631       imax(idxD0)=int(zmax)
4632       imax(1)=imax(idxD0)
4633       imax(2)=imax(idxD0)
4634
4635       do i=idxD0,nmax
4636         gam=gamD(i,iclpro,icltar)*b**2
4637         bet(i)=gam+betDp(i,iclpro,icltar)-alppar+1.
4638       enddo
4639       write(ifhi,'(a)')       'openhisto name gExact'
4640       write(ifhi,'(a)')       'htyp pfs'
4641       write(ifhi,'(a)')       'xmod lin ymod log'
4642       write(ifhi,'(a,2e11.3)')'xrange',0.,zmax
4643       write(ifhi,'(a)')       'yrange auto auto'
4644 c      write(ifhi,'(a)')'yrange 1.e-10 2'
4645       write(ifhi,'(a)')    'text 0 0 "xaxis z"'
4646       write(ifhi,'(a)')    'text 0 0 "yaxis g(z) "'
4647       write(ifhi,'(a)')       'array 2'
4648
4649       do ip0=0,imax(0)
4650         ip(0)=ip0
4651         do ip1=0,imax(1)
4652           ip(1)=ip1
4653           do ip2=0,imax(2)
4654             ip(2)=ip2
4655           t=0.
4656           do i=idxD0,nmax
4657             t=t+real(ip(i))*bet(i)
4658           enddo
4659           write(ifhi,'(2e14.6)')t,utgam2(dble(alplea(2))+1.D0)
4660      &       /utgam2(dble(alplea(2))+1.D0+dble(t))
4661         enddo
4662       enddo
4663       enddo
4664
4665       write(ifhi,'(a)')    '  endarray'
4666       write(ifhi,'(a)')    'closehisto plot 0-'
4667
4668
4669 c**********************************************************************
4670
4671       write(ifhi,'(a)')       'openhisto name gExpo'
4672       write(ifhi,'(a)')       'htyp lbu'
4673       write(ifhi,'(a)')       'xmod lin ymod log'
4674       write(ifhi,'(a,2e11.3)')'xrange',0.,zmax
4675       write(ifhi,'(a)')'yrange auto auto'
4676       write(ifhi,'(a)')    'text 0 0 "xaxis z"'
4677       write(ifhi,'(a)')    'text 0 0 "yaxis g(z)"'
4678       write(ifhi,'(a)')       'array 2'
4679
4680       do i=0,nptg
4681         t=zmax*(real(i)/real(nptg))
4682         write(ifhi,'(2e14.6)') t,dexp(-dble(t))
4683       enddo
4684
4685       write(ifhi,'(a)')    '  endarray'
4686       write(ifhi,'(a)')    'closehisto plot 0-'
4687
4688
4689 c**********************************************************************
4690
4691
4692
4693       write(ifhi,'(a)')       'openhisto name gPower'
4694       write(ifhi,'(a)')       'htyp poc'
4695       write(ifhi,'(a)')       'xmod log ymod log'
4696       write(ifhi,'(a,2e11.3)')'xrange',0.,zmax
4697       write(ifhi,'(a)')       'yrange auto auto'
4698       write(ifhi,'(a)')    'text 0 0 "xaxis t"'
4699       write(ifhi,'(2a)')    'text 0 0 ',
4700      &     '"yaxis [P][G](1+[a]?L!)/[G](1+[a]?L!+[b])"'
4701       write(ifhi,'(a)')       'array 2'
4702
4703       do ip0=0,imax(0)
4704         ip(0)=ip0
4705         do ip1=0,imax(1)
4706           ip(1)=ip1
4707           do ip2=0,imax(2)
4708             ip(2)=ip2
4709 c$$$            do ip3=0,imax(3)
4710 c$$$              ip(3)=ip3
4711               t=0.
4712               do i=idxD0,nmax
4713                 t=t+real(ip(i))*bet(i)
4714               enddo
4715               write(ifhi,'(2e14.6)')t,xgammag2(iclpro,bet,ip,gamp)
4716 c$$$            enddo
4717           enddo
4718         enddo
4719       enddo
4720
4721       write(ifhi,'(a)')    '  endarray'
4722       write(ifhi,'(a)')    'closehisto plot 0'
4723
4724       end
4725
4726 c----------------------------------------------------------------------
4727       subroutine xParOmega1xy
4728 c----------------------------------------------------------------------
4729 c xpar2=b
4730 c xpar4=xh
4731 c xpar3=y
4732 c xpar7 : nucl coef
4733 c xpar8 : 2=xp/xm instead of xh/yp
4734 c----------------------------------------------------------------------
4735       include 'epos.incems'
4736       include 'epos.incsem'
4737       include 'epos.inc'
4738       include 'epos.incpar'
4739     !  double precision om1x,om1y
4740       double precision x,ranhis(0:51),y,ymax,xh
4741      &,om1xpk,om1xmk,t,om1xk,om1yk,xpr1,xmr1!,xp,xm,xmin,xmax
4742       common /psar7/ delx,alam3p,gam3p
4743
4744       nptg=50                  !number of point for the graphs
4745       biniDf=xpar2              !value of biniDf (impact paramter)
4746       xh=0.99d0
4747       xpr1=1.d0
4748       xmr1=1.d0
4749       if(xpar4.lt.1.)then
4750         xh=dble(xpar4**2)          !value of x
4751         xpr1=dble(xpar4)
4752         xmr1=dble(xpar4)
4753       endif
4754       do i=0,51
4755         ranhis(i)=0.d0
4756       enddo
4757 c$$$      xp=dsqrt(xh)*dble(exp(xpar3)) !y=xpar3
4758 c$$$      xm=1.d0
4759 c$$$      if(xp.ne.0.d0)xm=xh/xp
4760
4761
4762
4763       if(koll.eq.0)then
4764
4765       call xhistomran1(ranhis,biniDf)
4766
4767       stop'om1x not defined'
4768
4769
4770       write(ifhi,'(a)')       'openhisto name Om1x'
4771       write(ifhi,'(a)')       'htyp lin'
4772       write(ifhi,'(a)')       'xmod log ymod log'
4773       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4774       write(ifhi,'(a)')'yrange auto auto'
4775       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
4776       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x!"'
4777       if (xpar8.eq.1.) then
4778         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4779         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4780       endif
4781       write(ifhi,'(a)')       'array 2'
4782
4783       do i=0,nptg
4784         x=xminDf
4785         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
4786 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
4787         write(ifhi,*) x,   0   !om1x(x,biniDf)
4788       enddo
4789
4790       write(ifhi,'(a)')    '  endarray'
4791       write(ifhi,'(a)')    'closehisto plot 0-'
4792
4793       write(ifhi,'(a)')       'openhisto name Om1xRan'
4794       write(ifhi,'(a)')       'htyp his'
4795       write(ifhi,'(a)')       'xmod log ymod log'
4796       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4797       write(ifhi,'(a)')'yrange auto auto'
4798       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
4799       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x! random"'
4800       if (xpar8.eq.1.) then
4801         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4802         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4803       endif
4804       write(ifhi,'(a)')       'array 2'
4805
4806       do i=0,51
4807         x=xminDf
4808         if (i.ne.0) x=x*(xmaxDf/xminDf)**((dble(i)+.5d0)/51.d0)
4809 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
4810         write(ifhi,*) x,ranhis(i)
4811       enddo
4812
4813       write(ifhi,'(a)')    '  endarray'
4814       write(ifhi,'(a)')    'closehisto plot 0'
4815
4816       x=xh
4817       ymax=-.5D0*dlog(x)
4818       do i=0,51
4819         ranhis(i)=0.d0
4820       enddo
4821       call xhistomran2(ranhis,x,biniDf)
4822
4823       stop'om1y not defined'
4824
4825       write(ifhi,'(a)')       'openhisto name Om1y'
4826       write(ifhi,'(a)')       'htyp lin'
4827       write(ifhi,'(a)')       'xmod lin ymod lin'
4828       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
4829       write(ifhi,'(a)')'yrange auto auto'
4830       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
4831       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y!"'
4832       if (xpar8.eq.1.) then
4833         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4834         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4835       endif
4836       write(ifhi,'(a)')       'array 2'
4837
4838       do i=0,nptg
4839         y=-ymax+(2.d0*ymax)*(dble(i)/dble(nptg))
4840         write(ifhi,*) y, 0  !om1y(x,y,biniDf)
4841       enddo
4842
4843       write(ifhi,'(a)')    '  endarray'
4844       write(ifhi,'(a)')    'closehisto plot 0-'
4845
4846       write(ifhi,'(a)')       'openhisto name Om1yRan'
4847       write(ifhi,'(a)')       'htyp his'
4848       write(ifhi,'(a)')       'xmod lin ymod lin'
4849       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
4850       write(ifhi,'(a)')'yrange auto auto'
4851       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
4852       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y! random"'
4853       if (xpar8.eq.1.) then
4854         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4855         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4856       endif
4857       write(ifhi,'(a)')       'array 2'
4858
4859       do i=0,50
4860         y=-ymax+(2.d0*ymax)*(dble(i)/50.d0)
4861         write(ifhi,*) y,ranhis(i)
4862       enddo
4863
4864       write(ifhi,'(a)')    '  endarray'
4865       write(ifhi,'(a)')    'closehisto plot 0'
4866
4867
4868 c**********************************************************************
4869
4870       else
4871
4872       do k=1,koll
4873         bk(k)=biniDf
4874       enddo
4875       call GfunPark(0)
4876       call integom1(0)
4877
4878       if(xpar8.eq.2.)then
4879
4880         call xhistomran8(ranhis,xpr1,xmr1)
4881
4882
4883       write(ifhi,'(a)')       'openhisto name Om1xp'
4884       write(ifhi,'(a)')       'htyp lin'
4885       write(ifhi,'(a)')       'xmod log ymod log'
4886       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4887       write(ifhi,'(a)')'yrange auto auto'
4888       write(ifhi,'(a)')    'text 0 0 "xaxis X+"'
4889       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x+!"'
4890       if (xpar8.eq.1.) then
4891         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4892         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4893       endif
4894       write(ifhi,'(a)')       'array 2'
4895
4896       do i=0,nptg
4897         x=xminDf
4898         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
4899 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
4900         t=om1xpk(x,xpr1,xmr1,1)
4901         write(ifhi,*) x,t
4902       enddo
4903
4904       write(ifhi,'(a)')    '  endarray'
4905       write(ifhi,'(a)')    'closehisto plot 0-'
4906
4907       write(ifhi,'(a)')       'openhisto name Om1xpRan'
4908       write(ifhi,'(a)')       'htyp his'
4909       write(ifhi,'(a)')       'xmod log ymod log'
4910       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4911       write(ifhi,'(a)')'yrange auto auto'
4912       write(ifhi,'(a)')    'text 0 0 "xaxis X+"'
4913       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x+! random"'
4914       if (xpar8.eq.1.) then
4915         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4916         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4917       endif
4918       write(ifhi,'(a)')       'array 2'
4919
4920       do i=0,51
4921         x=xminDf
4922         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
4923 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
4924         write(ifhi,*) x,ranhis(i)
4925       enddo
4926
4927       write(ifhi,'(a)')    '  endarray'
4928       write(ifhi,'(a)')    'closehisto plot 0'
4929
4930       do i=0,51
4931         ranhis(i)=0.d0
4932       enddo
4933
4934       call xhistomran9(ranhis,xh,xpr1,xmr1)
4935
4936
4937       write(ifhi,'(a)')       'openhisto name Om1xm'
4938       write(ifhi,'(a)')       'htyp lin'
4939       write(ifhi,'(a)')       'xmod log ymod log'
4940       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4941       write(ifhi,'(a)')'yrange auto auto'
4942       write(ifhi,'(a)')    'text 0 0 "xaxis X-"'
4943       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x-!"'
4944       if (xpar8.eq.1.) then
4945         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4946         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4947         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x+=',xh,'"'
4948       endif
4949       write(ifhi,'(a)')       'array 2'
4950
4951       do i=0,nptg
4952         x=xminDf
4953         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
4954 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
4955         t=om1xmk(xh,x,xpr1,xmr1,1)
4956         write(ifhi,*) x,t
4957       enddo
4958
4959       write(ifhi,'(a)')    '  endarray'
4960       write(ifhi,'(a)')    'closehisto plot 0-'
4961
4962       write(ifhi,'(a)')       'openhisto name Om1xmRan'
4963       write(ifhi,'(a)')       'htyp his'
4964       write(ifhi,'(a)')       'xmod log ymod log'
4965       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4966       write(ifhi,'(a)')'yrange auto auto'
4967       write(ifhi,'(a)')    'text 0 0 "xaxis X-"'
4968       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x-! random"'
4969       if (xpar8.eq.1.) then
4970         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
4971         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
4972         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.8 "x+=',xh,'"'
4973       endif
4974       write(ifhi,'(a)')       'array 2'
4975
4976       do i=0,51
4977         x=xminDf
4978         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
4979 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
4980         write(ifhi,*) x,ranhis(i)
4981       enddo
4982
4983       write(ifhi,'(a)')    '  endarray'
4984       write(ifhi,'(a)')    'closehisto plot 0'
4985
4986 c**********************************************************************
4987       else
4988
4989       call xhistomran10(ranhis)
4990
4991
4992       write(ifhi,'(a)')       'openhisto name Om1xk'
4993       write(ifhi,'(a)')       'htyp lin'
4994       write(ifhi,'(a)')       'xmod log ymod log'
4995       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
4996       write(ifhi,'(a)')'yrange auto auto'
4997       write(ifhi,'(a)')    'text 0 0 "xaxis X"'
4998       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x!"'
4999       if (xpar8.eq.1.) then
5000         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5001         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5002       endif
5003       write(ifhi,'(a)')       'array 2'
5004
5005       do i=0,nptg
5006         x=xminDf
5007         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5008 c.......x=xminDf+(xmaxDf-xminDf)*(dble(i)/dble(nptg))
5009         t=om1xk(x,1)
5010         write(ifhi,*) x,t
5011       enddo
5012
5013       write(ifhi,'(a)')    '  endarray'
5014       write(ifhi,'(a)')    'closehisto plot 0-'
5015
5016       write(ifhi,'(a)')       'openhisto name Om1xpRan'
5017       write(ifhi,'(a)')       'htyp his'
5018       write(ifhi,'(a)')       'xmod log ymod log'
5019       write(ifhi,'(a,2e11.3)')'xrange',xminDf,xmaxDf
5020       write(ifhi,'(a)')'yrange auto auto'
5021       write(ifhi,'(a)')    'text 0 0 "xaxis X+"'
5022       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1x+! random"'
5023       if (xpar8.eq.1.) then
5024         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5025         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5026       endif
5027       write(ifhi,'(a)')       'array 2'
5028
5029       do i=0,51
5030         x=xminDf
5031         if (i.ne.0) x=x*(xmaxDf/xminDf)**(dble(i)/dble(nptg))
5032 c.......x=xminDf+(xmaxDf-xminDf)*((dble(i)+.5d0)/51.d0)
5033         write(ifhi,*) x,ranhis(i)
5034       enddo
5035
5036       write(ifhi,'(a)')    '  endarray'
5037       write(ifhi,'(a)')    'closehisto plot 0'
5038
5039       ymax=-.5D0*dlog(xh)
5040       do i=0,51
5041         ranhis(i)=0.d0
5042       enddo
5043
5044       call xhistomran11(ranhis,xh)
5045
5046
5047
5048       write(ifhi,'(a)')       'openhisto name Om1yk'
5049       write(ifhi,'(a)')       'htyp lin'
5050       write(ifhi,'(a)')       'xmod lin ymod lin'
5051       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
5052       write(ifhi,'(a)')'yrange auto auto'
5053       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
5054       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y!"'
5055       if (xpar8.eq.1.) then
5056         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5057         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5058       endif
5059       write(ifhi,'(a)')       'array 2'
5060
5061       do i=0,nptg
5062         y=-ymax+(2.d0*ymax)*(dble(i)/dble(nptg))
5063         t=om1yk(xh,y,1)
5064         write(ifhi,*) y,t
5065       enddo
5066
5067       write(ifhi,'(a)')    '  endarray'
5068       write(ifhi,'(a)')    'closehisto plot 0-'
5069
5070       write(ifhi,'(a)')       'openhisto name Om1yRan'
5071       write(ifhi,'(a)')       'htyp his'
5072       write(ifhi,'(a)')       'xmod lin ymod lin'
5073       write(ifhi,'(a,2e11.3)')'xrange',-ymax-1.d0,ymax+1.d0
5074       write(ifhi,'(a)')'yrange auto auto'
5075       write(ifhi,'(a)')    'text 0 0 "xaxis Y"'
5076       write(ifhi,'(a)')    'text 0 0 "yaxis [w]?1y! random"'
5077       if (xpar8.eq.1.) then
5078         write(ifhi,'(a,e7.2,a)')  'text 0.1 0.9 "s=',engy**2,' GeV^2!"'
5079         write(ifhi,'(a,f5.2,a)')  'text 0.5 0.9 "b=',biniDf,' fm"'
5080       endif
5081       write(ifhi,'(a)')       'array 2'
5082
5083       do i=0,50
5084         y=-ymax+(2.d0*ymax)*(dble(i)/50.d0)
5085         write(ifhi,*) y,ranhis(i)
5086       enddo
5087
5088       write(ifhi,'(a)')    '  endarray'
5089       write(ifhi,'(a)')    'closehisto plot 0'
5090
5091
5092       endif
5093
5094       endif
5095
5096       return
5097       end
5098
5099 c----------------------------------------------------------------------
5100       subroutine xRanPt
5101 c----------------------------------------------------------------------
5102 c xpar2=xcut
5103 c----------------------------------------------------------------------
5104       include 'epos.inc'
5105       double precision ranhis(0:101)
5106
5107       nptg=100                  !number of point for the graphs
5108       xcut=xpar2              !value of biniDf (impact paramter)
5109       if(xcut.le.0.)xcut=100.
5110       do i=0,101
5111         ranhis(i)=0.d0
5112       enddo
5113 c$$$      xp=dsqrt(xh)*dble(exp(xpar3)) !y=xpar3
5114 c$$$      xm=1.d0
5115 c$$$      if(xp.ne.0.d0)xm=xh/xp
5116
5117
5118
5119       call xranptg(ranhis,xcut)
5120
5121
5122
5123       write(ifhi,'(a)')       'openhisto name ranpt'
5124       write(ifhi,'(a)')       'htyp lin'
5125       write(ifhi,'(a)')       'xmod lin ymod log'
5126       write(ifhi,'(a,2e11.3)')'xrange',0.,min(10.,xcut+1.)
5127       write(ifhi,'(a)')'yrange auto auto'
5128       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5129       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5130       write(ifhi,'(a)')       'array 2'
5131
5132
5133       do i=0,nptg
5134         x=float(i)/10.
5135         write(ifhi,*) x,ranhis(i)
5136       enddo
5137
5138       write(ifhi,'(a)')    '  endarray'
5139       write(ifhi,'(a)')    'closehisto plot 0-'
5140
5141       call xranpte(ranhis,xcut)
5142
5143
5144
5145       write(ifhi,'(a)')       'openhisto name ranpt'
5146       write(ifhi,'(a)')       'htyp lin'
5147       write(ifhi,'(a)')       'xmod lin ymod log'
5148       write(ifhi,'(a,2e11.3)')'xrange',0.,min(10.,xcut+1.)
5149       write(ifhi,'(a)')'yrange auto auto'
5150       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5151       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5152       write(ifhi,'(a)')       'array 2'
5153
5154
5155       do i=0,nptg
5156         x=float(i)/10.
5157         write(ifhi,*) x,ranhis(i)
5158       enddo
5159
5160       write(ifhi,'(a)')    '  endarray'
5161       write(ifhi,'(a)')    'closehisto plot 0-'
5162
5163       call xranpts(ranhis,xcut)
5164
5165
5166
5167       write(ifhi,'(a)')       'openhisto name ranpt'
5168       write(ifhi,'(a)')       'htyp lin'
5169       write(ifhi,'(a)')       'xmod lin ymod log'
5170       write(ifhi,'(a,2e11.3)')'xrange',0.,min(10.,xcut+1.)
5171       write(ifhi,'(a)')'yrange auto auto'
5172       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5173       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5174       write(ifhi,'(a)')       'array 2'
5175
5176
5177       do i=0,nptg
5178         x=float(i)/10.
5179         write(ifhi,*) x,ranhis(i)
5180       enddo
5181
5182       write(ifhi,'(a)')    '  endarray'
5183       write(ifhi,'(a)')    'closehisto plot 0-'
5184
5185       call xranptc(ranhis,xcut)
5186
5187
5188
5189       write(ifhi,'(a)')       'openhisto name ranpt'
5190       write(ifhi,'(a)')       'htyp lin'
5191       write(ifhi,'(a)')       'xmod lin ymod log'
5192       write(ifhi,'(a,2e11.3)')'xrange',0.,min(10.,xcut+1.)
5193       write(ifhi,'(a)')'yrange auto auto'
5194       write(ifhi,'(a)')    'text 0 0 "xaxis pt"'
5195       write(ifhi,'(a)')    'text 0 0 "yaxis P"'
5196       write(ifhi,'(a)')       'array 2'
5197
5198
5199       do i=0,nptg
5200         x=float(i)/10.
5201         write(ifhi,*) x,ranhis(i)
5202       enddo
5203
5204       write(ifhi,'(a)')    '  endarray'
5205       write(ifhi,'(a)')    'closehisto plot 0'
5206
5207
5208       return
5209       end
5210
5211 c----------------------------------------------------------------------
5212
5213       double precision function xgammag2(iclrem,bet,ip,gamp)
5214
5215 c----------------------------------------------------------------------
5216       include 'epos.inc'
5217       include 'epos.incsem'
5218       double precision utgam2
5219       dimension bet(idxD0:idxD1),ip(idxD0:idxD1)
5220
5221       xgammag2=1.d0
5222
5223       imax=idxD1
5224
5225       do i=idxD0,imax
5226       if(ip(i).ne.0) xgammag2=xgammag2
5227      &   *(utgam2(dble(alplea(iclrem))+1.d0+dble(gamp))
5228      &   /(max(0.d0,dble(int(gamp+0.5))+1))
5229      &   /utgam2(dble(alplea(iclrem)+bet(i)+gamp)+1.D0))
5230      &                                          **dble(ip(i))
5231       enddo
5232
5233       return
5234       end
5235
5236 c----------------------------------------------------------------------
5237
5238       function xsigmafit(x)
5239
5240 c----------------------------------------------------------------------
5241
5242       include 'epos.incpar'
5243       double precision x,xDfit,sfsh,varifit,range,sig2
5244       double precision bf(maxdataDf),Db(maxdataDf)
5245       external varifit
5246
5247
5248       sig2=bmaxDf/2.
5249       range=sig2
5250       xp=real(dsqrt(x))
5251       xm=xp
5252
5253
5254       sfsh=xDfit(0,1,smaxDf,xp,xm,0.)
5255       if(dabs(sfsh).ge.1.d-5)then
5256       do i=0,nptf-1
5257         bf(i+1)=dble(-bmaxDf+real(i)*2.*bmaxDf/real(nptf-1))
5258         Db(i+1)=xDfit(0,1,smaxDf,xp,xm,real(bf(i+1)))/sfsh
5259       enddo
5260
5261 c.....Fit of D(X,b) between -bmaxDf and bmaxDf
5262       call minfit(varifit,bf,Db,nptf,sig2,range)
5263
5264       xsigmafit=real(sig2)
5265       else
5266       xsigmafit=0.
5267       endif
5268
5269
5270       return
5271       end
5272
5273
5274 c----------------------------------------------------------------------
5275
5276       subroutine xhistomran1(histo,b)
5277
5278 c----------------------------------------------------------------------
5279 c.....Make Histogram of om1xr
5280 c----------------------------------------------------------------------
5281
5282       include 'epos.incpar'
5283
5284       double precision histo(0:51),x,x1  !,om1xr
5285       integer*4 n
5286
5287
5288       n=100000
5289       do i=0,51
5290         histo(i)=0.d0
5291       enddo
5292       do j=1,n
5293         if(mod(j,10000).eq.0)write(*,*)"x1",j
5294         x=0  !om1xr(b)
5295         stop'om1xr(b) not defined'
5296         if(x.lt.xminDf)goto 111
5297 c.........Exponential
5298           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5299 c.........Linear
5300 c.........k=int(x*50.d0)
5301         histo(k)=histo(k)+1.d0
5302  111  enddo
5303       do i=0,51
5304
5305 c.......Exponential
5306
5307         x=xminDf
5308         x1=xminDf
5309         x=x**(1.d0-dble(i)/51.d0)
5310         x1=x1**(1.d0-dble(i+1)/51.d0)
5311         if(i.eq.51)then
5312           x1=1.d0
5313           x=0.d0
5314         endif
5315         histo(i)=histo(i)/dble(n)/(x1-x)
5316
5317 c.......Linear
5318 c        histo(i)=histo(i)/dble(n)*51.d0
5319       enddo
5320
5321       return
5322       end
5323
5324
5325 c----------------------------------------------------------------------
5326
5327       subroutine xhistomran2(histo,xh,b)
5328
5329 c----------------------------------------------------------------------
5330 c.....Make Histogram of om1yr
5331 c----------------------------------------------------------------------
5332
5333       double precision histo(0:51),x,xh,dx,ymax   !,om1yr
5334       integer*4 n
5335
5336       ymax=-.5D0*dlog(xh)
5337       dx=ymax/25.d0
5338
5339       n=100000
5340       do i=0,50
5341         histo(i)=0.d0
5342       enddo
5343       do j=1,n
5344         if(mod(j,10000).eq.0)write(*,*)"y1",j
5345         x= 0 !  om1yr(xh,b)
5346         stop'om1yr(xh,b) not defined'
5347         k=int((x/ymax+1.d0)*25.d0)
5348 c.......write(*,*)x,k
5349         histo(k)=histo(k)+1.d0
5350       enddo
5351       do i=0,50
5352         histo(i)=histo(i)/dble(n)/dx
5353       enddo
5354
5355       return
5356       end
5357
5358
5359 cc----------------------------------------------------------------------
5360 c
5361 c      subroutine xhistomran6(histo,bx,by,bmax,del)
5362 c
5363 cc----------------------------------------------------------------------
5364 cc.....Make Histogram of b1 (impact parameter of vertex in Y and X)
5365 cc----------------------------------------------------------------------
5366 c
5367 c      double precision histo(0:51),dx
5368 c      integer*4 n
5369 c
5370 c      n=100000
5371 c      dx=dble(bmax)/50.d0
5372 c      do i=0,50
5373 c        histo(i)=0.d0
5374 c      enddo
5375 c      do j=1,n
5376 c        if(mod(j,10000).eq.0)write(*,*)"b1",j
5377 c        z=rangen()
5378 c        zp=rangen()
5379 c        bb1x=(bx+sqrt(-del*log(z))*cos(2.*3.14*zp))/2.
5380 c        bb1y=(by+sqrt(-del*log(z))*sin(2.*3.14*zp))/2.
5381 c        x=sqrt((bx-bb1x)*(bx-bb1x)+(by-bb1y)*(by-bb1y))
5382 c        k=int(x/bmax*50.)
5383 c        if(k.le.50)then
5384 c          histo(k)=histo(k)+1.d0
5385 c        else
5386 c          histo(51)=histo(51)+1.d0
5387 c        endif
5388 c      enddo
5389 c      do i=0,50
5390 c        histo(i)=histo(i)/dble(n)/dx
5391 c      enddo
5392 c
5393 c      return
5394 c      end
5395 c
5396
5397 c----------------------------------------------------------------------
5398
5399       subroutine xhistomran8(histo,xpr,xmr)
5400
5401 c----------------------------------------------------------------------
5402 c.....Make Histogram of om1xprk
5403 c----------------------------------------------------------------------
5404
5405       include 'epos.incpar'
5406
5407       double precision histo(0:51),x,x1,om1xprk,xpr,xmr
5408       integer*4 n
5409
5410
5411       n=100000
5412       do i=0,51
5413         histo(i)=0.d0
5414       enddo
5415       do j=1,n
5416         if(mod(j,10000).eq.0)write(*,*)"x+",j
5417           x=om1xprk(1,xpr,xmr,1)
5418         if(x.lt.xminDf)goto 111
5419 c.........Exponential
5420           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5421 c.........Linear
5422 c.........k=int(x*50.d0)
5423         histo(k)=histo(k)+1.d0
5424  111  enddo
5425       do i=0,51
5426
5427 c.......Exponential
5428
5429         x=xminDf
5430         x1=xminDf
5431         x=x**(1.d0-dble(i)/51.d0)
5432         x1=x1**(1.d0-dble(i+1)/51.d0)
5433         if(i.eq.51)then
5434           x1=1.d0
5435           x=0.d0
5436         endif
5437         histo(i)=histo(i)/dble(n)/(x1-x)*xpr
5438
5439 c.......Linear
5440 c        histo(i)=histo(i)/dble(n)*51.d0
5441       enddo
5442
5443       return
5444       end
5445
5446 c----------------------------------------------------------------------
5447
5448       subroutine xhistomran9(histo,xp,xpr,xmr)
5449
5450 c----------------------------------------------------------------------
5451 c.....Make Histogram of om1xmrk
5452 c----------------------------------------------------------------------
5453
5454       include 'epos.incpar'
5455
5456       double precision histo(0:51),x,x1,om1xmrk,xp,xpr,xmr
5457       integer*4 n
5458
5459
5460       n=100000
5461       do i=0,51
5462         histo(i)=0.d0
5463       enddo
5464       do j=1,n
5465         if(mod(j,10000).eq.0)write(*,*)"x-",j
5466           x=om1xmrk(1,xp,xpr,xmr,1)
5467         if(x.lt.xminDf)goto 111
5468 c.........Exponential
5469           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5470 c.........Linear
5471 c.........k=int(x*50.d0)
5472         histo(k)=histo(k)+1.d0
5473  111  enddo
5474       do i=0,51
5475
5476 c.......Exponential
5477
5478         x=xminDf
5479         x1=xminDf
5480         x=x**(1.d0-dble(i)/51.d0)
5481         x1=x1**(1.d0-dble(i+1)/51.d0)
5482         if(i.eq.51)then
5483           x1=1.d0
5484           x=0.d0
5485         endif
5486         histo(i)=histo(i)/dble(n)/(x1-x)*xmr
5487
5488 c.......Linear
5489 c        histo(i)=histo(i)/dble(n)*51.d0
5490       enddo
5491
5492       return
5493       end
5494
5495 c----------------------------------------------------------------------
5496
5497       subroutine xhistomran10(histo)
5498
5499 c----------------------------------------------------------------------
5500 c.....Make Histogram of om1xrk
5501 c----------------------------------------------------------------------
5502
5503       include 'epos.incpar'
5504
5505       double precision histo(0:51),x,x1,om1xrk
5506       integer*4 n
5507
5508
5509       n=100000
5510       do i=0,51
5511         histo(i)=0.d0
5512       enddo
5513       do j=1,n
5514         if(mod(j,10000).eq.0)write(*,*)"xk",j
5515         x=om1xrk(1)
5516         if(x.lt.xminDf)goto 111
5517 c.........Exponential
5518           k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5519 c.........Linear
5520 c.........k=int(x*50.d0)
5521         histo(k)=histo(k)+1.d0
5522  111  enddo
5523       do i=0,51
5524
5525 c.......Exponential
5526
5527         x=xminDf
5528         x1=xminDf
5529         x=x**(1.d0-dble(i)/51.d0)
5530         x1=x1**(1.d0-dble(i+1)/51.d0)
5531         if(i.eq.51)then
5532           x1=1.d0
5533           x=0.d0
5534         endif
5535         histo(i)=histo(i)/dble(n)/(x1-x)
5536
5537 c.......Linear
5538 c        histo(i)=histo(i)/dble(n)*51.d0
5539       enddo
5540
5541       return
5542       end
5543
5544
5545 c----------------------------------------------------------------------
5546
5547       subroutine xhistomran11(histo,xh)
5548
5549 c----------------------------------------------------------------------
5550 c.....Make Histogram of om1yrk
5551 c----------------------------------------------------------------------
5552
5553       double precision histo(0:51),x,xh,dx,om1yrk,ymax
5554       integer*4 n
5555
5556       ymax=-.5D0*dlog(xh)
5557       dx=ymax/25.d0
5558
5559       n=100000
5560       do i=0,50
5561         histo(i)=0.d0
5562       enddo
5563       do j=1,n
5564         if(mod(j,10000).eq.0)write(*,*)"yk",j
5565         x=om1yrk(xh,1)
5566         k=int((x/ymax+1.d0)*25.d0)
5567 c.......write(*,*)x,k
5568         histo(k)=histo(k)+1.d0
5569       enddo
5570       do i=0,50
5571         histo(i)=histo(i)/dble(n)/dx
5572       enddo
5573
5574       return
5575       end
5576
5577 c----------------------------------------------------------------------
5578
5579       subroutine xranptg(histo,xcut)
5580
5581 c----------------------------------------------------------------------
5582 c.....Make Histogram of random distribution
5583 c----------------------------------------------------------------------
5584
5585       include 'epos.incpar'
5586
5587       double precision histo(0:101)
5588       integer*4 n
5589
5590
5591       n=100000
5592       do i=0,101
5593         histo(i)=0.d0
5594       enddo
5595       do j=1,n
5596         if(mod(j,10000).eq.0)write(*,*)"ptg",j
5597 c .........exp(-x**2)
5598  12   x=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
5599
5600       if(xcut.gt.0.)then
5601         if(rangen().lt.x/xcut)goto 12
5602       endif
5603 c.........Exponential
5604 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5605 c.........Linear
5606         k=int(x*10.)
5607         k=min(k,101)
5608         histo(k)=histo(k)+1.d0
5609  111  enddo
5610       do i=0,101
5611
5612 c.......Exponential
5613
5614 c        x=xminDf
5615 c        x1=xminDf
5616 c        x=x**(1.d0-dble(i)/51.d0)
5617 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5618 c        if(i.eq.51)then
5619 c          x1=1.d0
5620 c          x=0.d0
5621 c        endif
5622 c        histo(i)=histo(i)/dble(n)/(x1-x)
5623
5624 c.......Linear
5625         histo(i)=histo(i)/dble(n)*101.d0
5626       enddo
5627
5628       return
5629       end
5630
5631 c----------------------------------------------------------------------
5632
5633       subroutine xranpte(histo,xcut)
5634
5635 c----------------------------------------------------------------------
5636 c.....Make Histogram of random distribution
5637 c----------------------------------------------------------------------
5638
5639       include 'epos.incpar'
5640
5641       double precision histo(0:101)
5642       integer*4 n
5643
5644
5645       n=100000
5646       do i=0,101
5647         histo(i)=0.d0
5648       enddo
5649       do j=1,n
5650         if(mod(j,10000).eq.0)write(*,*)"pte",j
5651 c .........exp(-x)
5652   12  xmx=50
5653       r=2.
5654       do while (r.gt.1.)
5655   11    x=sqrt(exp(rangen()*log(1+xmx**2))-1)
5656         if(x.eq.0.)goto11
5657         r=rangen()  /  ( exp(-x)*(1+x**2) )
5658       enddo
5659       x=x/2.
5660
5661       if(xcut.gt.0.)then
5662         if(rangen().lt.x/xcut)goto 12
5663       endif
5664 c.........Exponential
5665 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5666 c.........Linear
5667         k=int(x*10.)
5668         k=min(k,101)
5669         histo(k)=histo(k)+1.d0
5670  111  enddo
5671       do i=0,101
5672
5673 c.......Exponential
5674
5675 c        x=xminDf
5676 c        x1=xminDf
5677 c        x=x**(1.d0-dble(i)/51.d0)
5678 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5679 c        if(i.eq.51)then
5680 c          x1=1.d0
5681 c          x=0.d0
5682 c        endif
5683 c        histo(i)=histo(i)/dble(n)/(x1-x)
5684
5685 c.......Linear
5686         histo(i)=histo(i)/dble(n)*101.d0
5687       enddo
5688
5689       return
5690       end
5691
5692 c----------------------------------------------------------------------
5693
5694       subroutine xranpts(histo,xcut)
5695
5696 c----------------------------------------------------------------------
5697 c.....Make Histogram of random distribution
5698 c----------------------------------------------------------------------
5699
5700       include 'epos.incpar'
5701
5702       double precision histo(0:101)
5703       integer*4 n
5704
5705
5706       n=100000
5707       do i=0,101
5708         histo(i)=0.d0
5709       enddo
5710       do j=1,n
5711         if(mod(j,10000).eq.0)write(*,*)"pts",j
5712 c .........exp(-sqrt(x))
5713  12   xmx=500
5714       r=2.
5715       do while (r.gt.1.)
5716         x=sqrt(exp(rangen()*log(1+xmx**2))-1)
5717         r=rangen()  /  ( exp(-sqrt(x))*(1+x**2)/5. )
5718       enddo
5719       x=x/20.
5720
5721       if(xcut.gt.0.)then
5722         if(rangen().lt.x/xcut)goto 12
5723       endif
5724 c.........Exponential
5725 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5726 c.........Linear
5727         k=int(x*10.)
5728         k=min(k,101)
5729         histo(k)=histo(k)+1.d0
5730  111  enddo
5731       do i=0,101
5732
5733 c.......Exponential
5734
5735 c        x=xminDf
5736 c        x1=xminDf
5737 c        x=x**(1.d0-dble(i)/51.d0)
5738 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5739 c        if(i.eq.51)then
5740 c          x1=1.d0
5741 c          x=0.d0
5742 c        endif
5743 c        histo(i)=histo(i)/dble(n)/(x1-x)
5744
5745 c.......Linear
5746         histo(i)=histo(i)/dble(n)*101.d0
5747       enddo
5748
5749       return
5750       end
5751
5752 c----------------------------------------------------------------------
5753
5754       subroutine xranptc(histo,xcut)
5755
5756 c----------------------------------------------------------------------
5757 c.....Make Histogram of random distribution
5758 c----------------------------------------------------------------------
5759
5760       include 'epos.incpar'
5761
5762       double precision histo(0:101)
5763       integer*4 n
5764
5765
5766       n=100000
5767       do i=0,101
5768         histo(i)=0.d0
5769       enddo
5770       do j=1,n
5771         if(mod(j,10000).eq.0)write(*,*)"ptc",j
5772
5773         x=ranptcut(xcut)
5774 c.........Exponential
5775 c          k=int((-dlog(x)/dlog(xminDf)+1.d0)*51.d0)
5776 c.........Linear
5777         k=int(x*10.)
5778         k=min(k,101)
5779         histo(k)=histo(k)+1.d0
5780  111  enddo
5781       do i=0,101
5782
5783 c.......Exponential
5784
5785 c        x=xminDf
5786 c        x1=xminDf
5787 c        x=x**(1.d0-dble(i)/51.d0)
5788 c        x1=x1**(1.d0-dble(i+1)/51.d0)
5789 c        if(i.eq.51)then
5790 c          x1=1.d0
5791 c          x=0.d0
5792 c        endif
5793 c        histo(i)=histo(i)/dble(n)/(x1-x)
5794
5795 c.......Linear
5796         histo(i)=histo(i)/dble(n)*101.d0
5797       enddo
5798
5799       return
5800       end
5801
5802
5803