EPOS/epos167/epos-par-129.f

   1 c 15/02/2005 epos 1.03
   2
   3
   4 c----------------------------------------------------------------------
   5       subroutine paramini(imod)
   6 c----------------------------------------------------------------------
   7 c  Set  parameters of the parametrisation of the eikonals.
   8 c
   9 c xDfit=Sum(i=0,1)(alpD(i)*xp**betDp(i)*xm**betDpp(i)*s**betD(i)
  10 c                            *xs**(gamD(i)*b2)*exp(-b2/delD(i))
  11 c
  12 c Parameters stored in /Dparam/ (epos.inc)
  13 c if imod=0, do settings only for iclpro, if imod=1, do settings
  14 c for iclpro=2 and iclpro
  15 c----------------------------------------------------------------------
  16
  17       include 'epos.inc'
  18       include 'epos.incems'
  19       include 'epos.incsem'
  20       include 'epos.incpar'
  21       double precision PhiExact,y
  22       call utpri('parini',ish,ishini,3)
  23
  24 c Initialisation of the variables
  25
  26       call Class('paramini  ')
  27
  28 c Variables used for xparg (only graphical application)
  29
  30       spmin=4.*q2min         !Definition of spmin in psvin (epos-sha)
  31       sfshlim=5.*spmin          !transition energy for soft->hard
  32       emaxDf=engy            !energy for fit subroutines
  33       smaxDf=emaxDf**2       !energy squared for fit subroutines
  34 c      sfshlim=100.
  35       nptf=10                !number of point for the fit
  36       xmaxDf=1.d0            !maximum limit for the fit
  37       xminDf=1.d0/dble(smaxDf)   !minimum limit
  38       xshmin=3.d-2           !minimum limit for sh variance fit
  39       if(smaxDf.lt.sfshlim) xshmin=1.d0
  40       xfitmin=0.1d0   !sh fitted between f(1) and f(1)*xfitmin
  41       if(smaxDf.lt.sfshlim) xfitmin=1.d0
  42
  43       bmaxDf=2.              !maximum b for variance fit
  44
  45       idxDmin=idxD0          !minimum indice for parameters
  46       ntymin=ntymi           !minimum indice for diagram type
  47
  48       ucfpro=utgam1(1.+alplea(iclpro))
  49       ucftar=utgam1(1.+alplea(icltar))
  50
  51
  52
  53       iiiclegy=iclegy
  54
  55 c for pi or K - p crosse section calculation, we need alpD, ... for
  56 c iclpro=1 or 3 and iclpro=2
  57         iiiclpro1=iclpro
  58         iiidirec=1
  59         if(imod.eq.0)then
  60           iiiclpro2=iclpro
  61         else
  62           iiiclpro2=2
  63           if(iiiclpro1.lt.iiiclpro2)iiidirec=-1
  64         endif
  65         iclprosave=iclpro
  66
  67         do iiipro=iiiclpro2,iiiclpro1,iiidirec
  68
  69           iclpro=iiipro
  70
  71       if(ish.ge.4)write(ifch,*)'gamini'
  72      &          ,iscreen,iclpro,icltar,iclegy,engy,sfshlim
  73
  74       if(isetcs.le.1)then                      !if set mode, do fit
  75
  76 c First try for fit parameters
  77
  78 c linear fit of the 2 components of G as a function of x
  79         call pompar(alpsf,betsf,0) !soft (taking into account hard)
  80         call pompar(alpsh,betsh,1) !sh
  81 c Gaussian fit of the 2 components of G as a function of x and b
  82         call variance(delsf,gamsf,0)
  83         call variance(delsh,gamsh,1)
  84         gamsf=max(0.,gamsf)
  85         gamsh=max(0.,gamsh)
  86
  87
  88 c Fit GFF
  89
  90
  91 c       fit parameters
  92         numminDf=3             !minimum indice
  93         numparDf=4             !maximum indice
  94         betac=100.               !temperature for chi2
  95         betae=1.               !temperature for error
  96         fparDf=0.8             !variation amplitude for range
  97
  98         nmcxDf=20000          !number of try for x fit
  99
 100 c starting values
 101
 102         parDf(1,3)=alpsf
 103         parDf(2,3)=betsf
 104         parDf(3,3)=alpsh
 105         parDf(4,3)=betsh
 106
 107         if(smaxDf.ge.3.*sfshlim)then
 108           call paramx           !alpD and betD
 109
 110           call pompar(alpsf,betsf,-1) !soft (taking into account hard)
 111           parDf(1,3)=alpsf
 112           parDf(2,3)=max(-0.95+alppar,betsf)
 113
 114         endif
 115
 116         alpsf=parDf(1,3)
 117         betsf=parDf(2,3)
 118         alpsh=parDf(3,3)
 119         betsh=parDf(4,3)
 120
 121
 122       else                                     !else parameters from table (inirj)
 123         nbpsf=iDxD0
 124         if(iclegy2.gt.1)then
 125           al=1.+(log(engy)-log(egylow))/log(egyfac) !energy class
 126           i2=min(iiiclegy+1,iclegy2)
 127           i1=i2-1
 128         else
 129           i1=iclegy
 130           i2=iclegy
 131           al=float(iclegy)
 132         endif
 133         dl=al-i1
 134         dl1=max(0.,1.-dl)
 135                                 !linear interpolation
 136         alpsf=alpDs(nbpsf,i2,iclpro,icltar)*dl
 137      &       +alpDs(nbpsf,i1,iclpro,icltar)*dl1
 138         alpsh=alpDs(1,i2,iclpro,icltar)*dl
 139      &       +alpDs(1,i1,iclpro,icltar)*dl1
 140         betsf=betDs(nbpsf,i2,iclpro,icltar)*dl
 141      &       +betDs(nbpsf,i1,iclpro,icltar)*dl1
 142         betsh=betDs(1,i2,iclpro,icltar)*dl
 143      &       +betDs(1,i1,iclpro,icltar)*dl1
 144         gamsf=gamDs(nbpsf,i2,iclpro,icltar)*dl
 145      &       +gamDs(nbpsf,i1,iclpro,icltar)*dl1
 146         gamsh=gamDs(1,i2,iclpro,icltar)*dl
 147      &       +gamDs(1,i1,iclpro,icltar)*dl1
 148         delsf=delDs(nbpsf,i2,iclpro,icltar)*dl
 149      &       +delDs(nbpsf,i1,iclpro,icltar)*dl1
 150         delsh=delDs(1,i2,iclpro,icltar)*dl
 151      &       +delDs(1,i1,iclpro,icltar)*dl1
 152
 153
 154 c For the Plots
 155         parDf(1,3)=alpsf
 156         parDf(2,3)=betsf
 157         parDf(3,3)=alpsh
 158         parDf(4,3)=betsh
 159
 160       endif
 161
 162 c     if energy too small to have semi-hard interaction -> only soft diagram
 163
 164       if(smaxDf.lt.sfshlim.and.idxD0.eq.0)then !no hard: soft+hard=soft
 165         alpsf=alpsf/2.
 166         alpsh=alpsf
 167         betsh=betsf
 168         gamsh=gamsf
 169         delsh=delsf
 170       endif
 171
 172
 173
 174 c Print results
 175
 176       if(ish.ge.4)then
 177         write(ifch,*)"parameters for iclpro:",iclpro
 178         write(ifch,*)"alp,bet,gam,del sf:",alpsf,betsf,gamsf,delsf
 179         write(ifch,*)"alp,bet,gam,del sh:",alpsh,betsh,gamsh,delsh
 180       endif
 181
 182
 183 c Record parameters
 184
 185
 186       alpD(idxD0,iclpro,icltar)=alpsf
 187       alpDp(idxD0,iclpro,icltar)=0.
 188       alpDpp(idxD0,iclpro,icltar)=0.
 189       betD(idxD0,iclpro,icltar)=betsf
 190       betDp(idxD0,iclpro,icltar)=betsf
 191       betDpp(idxD0,iclpro,icltar)=betsf
 192       gamD(idxD0,iclpro,icltar)=gamsf
 193       delD(idxD0,iclpro,icltar)=delsf
 194
 195       alpD(1,iclpro,icltar)=alpsh
 196       alpDp(1,iclpro,icltar)=0.
 197       alpDpp(1,iclpro,icltar)=0.
 198       betD(1,iclpro,icltar)=betsh
 199       betDp(1,iclpro,icltar)=betsh
 200       betDpp(1,iclpro,icltar)=betsh
 201       gamD(1,iclpro,icltar)=gamsh
 202       delD(1,iclpro,icltar)=delsh
 203
 204       if(iomega.lt.2.and.alpdif.ne.1.)then
 205         alpDp(2,iclpro,icltar)=0.
 206         alpDpp(2,iclpro,icltar)=0.
 207         betD(2,iclpro,icltar)=0. !max(0.,betsf)
 208 c        alpdifs=alpdif
 209 c        alpdif=0.99
 210         betDp(2,iclpro,icltar)=-alpdif+alppar
 211         betDpp(2,iclpro,icltar)=-alpdif+alppar
 212 c        alpD(2,iclpro,icltar)=(alpsf+alpsh)*wdiff(iclpro)*wdiff(icltar)
 213         alpD(2,iclpro,icltar)=wdiff(iclpro)*wdiff(icltar)
 214      &            /utgam1(1.-alpdif)**2
 215      &            *utgam1(2.-alpdif+alplea(iclpro))
 216      &            *utgam1(2.-alpdif+alplea(icltar))
 217      &            /chad(iclpro)/chad(icltar)
 218 c        alpdif=alpdifs
 219         gamD(2,iclpro,icltar)=0.
 220         delD(2,iclpro,icltar)=4.*.0389*(gwidth*(r2had(iclpro)
 221      &                    +r2had(icltar))+slopoms*log(smaxDf))
 222       else
 223         alpD(2,iclpro,icltar)=0.
 224         alpDp(2,iclpro,icltar)=0.
 225         alpDpp(2,iclpro,icltar)=0.
 226         betD(2,iclpro,icltar)=0.
 227         betDp(2,iclpro,icltar)=0.
 228         betDpp(2,iclpro,icltar)=0.
 229         gamD(2,iclpro,icltar)=0.
 230         delD(2,iclpro,icltar)=1.
 231       endif
 232       if(ish.ge.4)write(ifch,*)"alp,bet,betp,del dif:"
 233      &   ,alpD(2,iclpro,icltar),betD(2,iclpro,icltar)
 234      &   ,betDp(2,iclpro,icltar),delD(2,iclpro,icltar)
 235
 236
 237       bmxdif(iclpro,icltar)=conbmxdif()     !important to do it before kfit,  because it's used in.
 238 c          call Kfit(-1)    !xkappafit not used : if arg=-1, set xkappafit to 1
 239       if(isetcs.le.1)then
 240         if(isetcs.eq.0)then
 241           call Kfit(-1)         !xkappafit not used : if arg=-1, set xkappafit to 1)
 242         else
 243 c          call Kfit(-1)    !xkappafit not used : if arg=-1, set xkappafit to 1)
 244           call Kfit(iclegy)
 245         endif
 246 c for plots record alpDs, betDs, etc ...
 247         alpDs(idxD0,iclegy,iclpro,icltar)=alpsf
 248         alpDps(idxD0,iclegy,iclpro,icltar)=0.
 249         alpDpps(idxD0,iclegy,iclpro,icltar)=0.
 250         betDs(idxD0,iclegy,iclpro,icltar)=betsf
 251         betDps(idxD0,iclegy,iclpro,icltar)=betsf
 252         betDpps(idxD0,iclegy,iclpro,icltar)=betsf
 253         gamDs(idxD0,iclegy,iclpro,icltar)=gamsf
 254         delDs(idxD0,iclegy,iclpro,icltar)=delsf
 255
 256         alpDs(1,iclegy,iclpro,icltar)=alpsh
 257         alpDps(1,iclegy,iclpro,icltar)=0.
 258         alpDpps(1,iclegy,iclpro,icltar)=0.
 259         betDs(1,iclegy,iclpro,icltar)=betsh
 260         betDps(1,iclegy,iclpro,icltar)=betsh
 261         betDpps(1,iclegy,iclpro,icltar)=betsh
 262         gamDs(1,iclegy,iclpro,icltar)=gamsh
 263         delDs(1,iclegy,iclpro,icltar)=delsh
 264       endif
 265
 266       enddo
 267
 268       if(iclpro.ne.iclprosave)stop'problem in parini with iclpro'
 269
 270       if(ish.ge.4)then             !check PhiExact value for x=1
 271         y=PhiExact(1.,1.d0,1.d0,smaxDf,0.)
 272         write(ifch,*)'PhiExact=',y
 273       endif
 274
 275       call utprix('parini',ish,ishini,3)
 276
 277       return
 278       end
 279
 280 c----------------------------------------------------------------------
 281       subroutine Class(text)
 282 c----------------------------------------------------------------------
 283       include 'epos.inc'
 284       include 'epos.incpar'
 285       parameter (eps=1.e-5)    !to correct for precision problem)
 286       character*10 text
 287       if(iclegy1.eq.iclegy2)then
 288         iclegy=iclegy1
 289       else
 290       iclegy=1+int( (log(engy)-log(egylow))/log(egyfac) + eps ) !energy class
 291       if(iclegy.gt.iclegy2)then
 292          write(ifch,*)'***********************************************'
 293          write(ifch,*)'Warning in ',text
 294          write(ifch,*)'Energy above the range used for the fit of D:'
 295          write(ifch,*)egylow*egyfac**(iclegy1-1),egylow*egyfac**iclegy2
 296          write(ifch,*)'***********************************************'
 297          iclegy=iclegy2
 298       endif
 299       if(iclegy.lt.iclegy1)then
 300          write(ifch,*)'***********************************************'
 301          write(ifch,*)'Warning in ',text
 302          write(ifch,*)'Energy below the range used for the fit of D:'
 303          write(ifch,*)egylow*egyfac**(iclegy1-1),egylow*egyfac**iclegy2
 304          write(ifch,*)'***********************************************'
 305          iclegy=iclegy1
 306       endif
 307       endif
 308       end
 309
 310 c----------------------------------------------------------------------
 311       subroutine param
 312 c----------------------------------------------------------------------
 313 c  Set the parameter of the parametrisation of the eikonale.
 314 c  We group the parameters into 4 array with a dimension of idxD1(=1)
 315 c  to define xDfit (see below).
 316 c
 317 c xDfit=Sum(i,0,1)(alpD(i)*xp**betDp(i)*xm**betDpp(i)*s**betD(i)
 318 c                            *xs**(gamD(i)*b2)*exp(-b2/delD(i))
 319 c
 320 c subroutine used for tabulation.
 321 c----------------------------------------------------------------------
 322       include 'epos.inc'
 323       include 'epos.incsem'
 324       include 'epos.incpar'
 325
 326 c Initialisation of the variables
 327
 328       emaxDf=egyfac**(iclegy-1)*egylow
 329       smaxDf=emaxDf**2
 330
 331       spmin=4.*q2min         !Definition of spmin in psvin (epos-sha)
 332       sfshlim=5.*spmin
 333       nptf=10
 334       xmaxDf=1.d0
 335       xminDf=1d0/dble(smaxDf)
 336       xshmin=3.d-2           !minimum limit for sh variance fit
 337       if(smaxDf.lt.sfshlim) xshmin=1.d0
 338       xfitmin=0.1d0   !sh fitted between f(1) and f(1)*xfitmin
 339       if(smaxDf.lt.sfshlim) xfitmin=1.d0
 340       bmaxDf=2.
 341
 342
 343       if(idxD0.ne.0.and.idxD1.ne.1) stop "* idxD0/1 are not good! *"
 344
 345       engytmp=engy
 346       engy=emaxDf
 347
 348 c Initialisation of the parameters
 349
 350       do i=1,nbpf
 351         do j=1,4
 352           parDf(i,j)=1.
 353         enddo
 354       enddo
 355
 356 c.......Calculations of the parameters
 357
 358 c First try for fit parameters
 359
 360 c linear fit of the 2 components of G as a function of x
 361       call pompar(alpsf,betsf,0) !soft
 362       call pompar(alpsh,betsh,1) !sh
 363 c Gaussian fit of the 2 components of G as a function of x and b
 364       call variance(delsf,gamsf,0)
 365       call variance(delsh,gamsh,1)
 366       gamsf=max(0.,gamsf)
 367       gamsh=max(0.,gamsh)
 368
 369 c Fit GFF
 370
 371 c      fit parameters
 372       numminDf=3                !minimum indice
 373       numparDf=4                !maximum indice
 374       betac=100.                 !temperature for chi2
 375       betae=1.                 !temperature for error
 376       fparDf=0.8                 !variation amplitude for range
 377
 378       nmcxDf=20000              !number of try for x fit
 379
 380 c starting values
 381
 382       parDf(1,3)=alpsf
 383       parDf(2,3)=betsf
 384       parDf(3,3)=alpsh
 385       parDf(4,3)=betsh
 386
 387       if(smaxDf.ge.3.*sfshlim)then
 388         call paramx             !alpD and betD
 389
 390         call pompar(alpsf,betsf,-1) !soft (taking into account hard)
 391         parDf(1,3)=alpsf
 392         parDf(2,3)=max(-0.95+alppar,betsf)
 393
 394       endif
 395
 396       alpsf=parDf(1,3)
 397       betsf=parDf(2,3)
 398       alpsh=parDf(3,3)
 399       betsh=parDf(4,3)
 400
 401       if(ish.ge.4)then
 402         write(ifch,*)"param: fit parameters :",iscreen,iclpro,icltar
 403      *                                        ,iclegy,engy
 404         write(ifch,*)"alp,bet,gam,del sf:",alpsf,betsf,gamsf,delsf
 405         write(ifch,*)"alp,bet,gam,del sh:",alpsh,betsh,gamsh,delsh
 406       endif
 407
 408       alpDs(idxD0,iclegy,iclpro,icltar)=alpsf
 409       alpDps(idxD0,iclegy,iclpro,icltar)=betsf
 410       alpDpps(idxD0,iclegy,iclpro,icltar)=0.
 411       betDs(idxD0,iclegy,iclpro,icltar)=betsf
 412       betDps(idxD0,iclegy,iclpro,icltar)=betsf
 413       betDpps(idxD0,iclegy,iclpro,icltar)=betsf
 414       gamDs(idxD0,iclegy,iclpro,icltar)=gamsf
 415       delDs(idxD0,iclegy,iclpro,icltar)=delsf
 416
 417       alpDs(1,iclegy,iclpro,icltar)=alpsh
 418       alpDps(1,iclegy,iclpro,icltar)=betsh
 419       alpDpps(1,iclegy,iclpro,icltar)=0.
 420       betDs(1,iclegy,iclpro,icltar)=betsh
 421       betDps(1,iclegy,iclpro,icltar)=betsh
 422       betDpps(1,iclegy,iclpro,icltar)=betsh
 423       gamDs(1,iclegy,iclpro,icltar)=gamsh
 424       delDs(1,iclegy,iclpro,icltar)=delsh
 425
 426       engy=engytmp
 427
 428       return
 429       end
 430
 431
 432 c----------------------------------------------------------------------
 433
 434         subroutine pompar(alpha,beta,iqq)
 435
 436 c----------------------------------------------------------------------
 437 c  Return the power beta and the factor alpha of the fit of the eikonal
 438 c of a pomeron of type iqq : D(X)=alpha*(X)**beta
 439 c----------------------------------------------------------------------
 440
 441       include 'epos.inc'
 442       include 'epos.incsem'
 443       include 'epos.incpar'
 444       double precision X,D1,D0,X0,D,droot
 445       double precision Dsoftshval,xmax
 446       dimension xlnXs(maxdataDf),xlnD(maxdataDf),sigma(maxdataDf)
 447
 448       do i=1,nptf
 449         sigma(i)=1.e-2
 450       enddo
 451
 452       if(iqq.le.0) then
 453
 454         iscr=iqq
 455         xmax=min(0.1d0,10.d0*xminDf)
 456         X0=xminDf
 457         if(ish.ge.4)write (ifch,*) 'pompar (0) x0,xmax=',X0,xmax
 458
 459         do i=0,nptf-1
 460           X=X0
 461           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
 462           D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
 463           if(D.eq.1.d-10)then
 464             write(ifch,*)
 465      &    "Warning in pompar ! Dsoftshval(0) could be negative"
 466             sigma(i+1)=1.e5
 467           endif
 468           xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
 469           xlnD(i+1)=real(dlog(D))
 470         enddo
 471
 472
 473 c Fit of D(X) between X0 and xmax
 474
 475         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
 476         if(beta.gt.10.)beta=10.
 477
 478         alpha=real(Dsoftshval(real(X0)*smaxDf,X0,0.d0,0.,iscr))
 479      &       *(real(X0)*smaxDf)**(-beta)
 480
 481
 482       elseif(iqq.eq.1.and.xfitmin.ne.1.d0) then
 483
 484         iscr=2
 485         xmax=1.d0
 486
 487 c Definition of D0=D(X0)
 488
 489         D1=Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr)
 490
 491         D0=xfitmin*D1
 492
 493 c Calculation of X0 and D(X)
 494
 495         X0=droot(D0,D1,xmax,iscr)
 496         if(ish.ge.4)write (ifch,*) 'pompar (1) x0,xmax=',X0,xmax
 497
 498         do i=0,nptf-1
 499           X=X0
 500           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
 501           D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
 502           if(D.eq.1.d-10)then
 503             write(ifch,*)
 504      &    "Warning in pompar ! Dsoftshval(1) could be negative"
 505             sigma(i+1)=1.e5
 506           endif
 507           xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
 508           xlnD(i+1)=real(dlog(D))
 509         enddo
 510
 511
 512 c Fit of D(X) between X0 and xmax
 513
 514         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
 515         if(beta.gt.10.)beta=10.
 516
 517         alpha=real(Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr))
 518      &       *(real(xmax)*smaxDf)**(-beta)
 519
 520
 521       elseif(iqq.eq.10.and.xfitmin.ne.1.d0) then    ! iqq=10
 522
 523         iscr=0
 524         xmax=1.d0 !2.d0/max(2.d0,dlog(dble(smaxDf)/1.d3))
 525
 526 c Definition of D0=D(X0)
 527
 528         D1=Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr)
 529
 530         D0=xfitmin*D1
 531
 532 c Calculation of X0 and D(X)
 533
 534         X0=droot(D0,D1,xmax,iscr)
 535         if(ish.ge.4)write (ifch,*) 'pompar (1) x0,xmax=',X0,xmax
 536
 537         do i=0,nptf-1
 538           X=X0
 539           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
 540           D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
 541           if(D.eq.1.d-10)then
 542             write(ifch,*)
 543      &    "Warning in pompar ! Dsoftshval(10) could be negative"
 544             sigma(i+1)=1.e5
 545           endif
 546           xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
 547           xlnD(i+1)=real(dlog(D))
 548         enddo
 549
 550
 551 c Fit of D(X) between X0 and xmax
 552
 553         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
 554         if(beta.gt.10.)beta=10.
 555
 556         alpha=real(Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr))
 557      &       *(real(xmax)*smaxDf)**(-beta)
 558
 559
 560
 561       else                      !iqq=-1 or iqq=1 and xfitmin=1
 562
 563 c Calculation of X0 and D(X)
 564         iscr=0
 565
 566         X0=1.d0/dble(smaxDf)
 567         xmax=max(2.d0/dble(smaxDf),
 568      &       min(max(0.03d0,dble(smaxDf)/2.d5),0.1d0))
 569
 570         if(ish.ge.4)write (ifch,*) 'pompar (-1) x0,xmax=',X0,xmax
 571
 572         do i=0,nptf-1
 573           X=X0
 574           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
 575           D=max(1.d-10,Dsoftshval(real(X)*smaxDf,X,0.d0,0.,iscr))
 576           if(D.eq.1.d-10)then
 577             write(ifch,*)
 578      &    "Warning in pompar ! Dsoftshval(-1) could be negative"
 579             sigma(i+1)=1.e5
 580           endif
 581           xlnXs(i+1)=real(dlog(X*dble(smaxDf)))
 582           xlnD(i+1)=real(dlog(D))
 583         enddo
 584
 585 c Fit of D(X) between X0 and xmax
 586
 587         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
 588         if(beta.gt.10.)beta=10.
 589
 590
 591         alpha=real(Dsoftshval(real(xmax)*smaxDf,xmax,0.d0,0.,iscr))
 592      &           *(real(xmax)*smaxDf)**(-beta)
 593       endif
 594
 595         if(ish.ge.4)write(ifch,*) '%%%%%%%%%%%%% pompar %%%%%%%%%%%%%'
 596         if(ish.ge.4)write(ifch,*) 'alpD ini =',alpha,' betD ini=',beta
 597
 598       return
 599       end
 600
 601 c----------------------------------------------------------------------
 602
 603       double precision function droot(d0,d1,xmax,iscr)
 604
 605 c----------------------------------------------------------------------
 606 c Find x0 which gives f(x0)=D(x0*S)-d0=0
 607 c iqq=0 soft pomeron
 608 c iqq=1 semi-hard pomeron
 609 c----------------------------------------------------------------------
 610       include 'epos.inc'
 611       include 'epos.incsem'
 612       include 'epos.incpar'
 613       double precision Dsoftshval,d0,d1,d2,x0,x1,x2,f0,f1,f2,xmax
 614       parameter (kmax=1000)
 615
 616
 617       k=0
 618       x0=min(xfitmin,100.d0*xminDf)
 619       x1=xmax
 620  5    d2=dabs(Dsoftshval(real(x0)*smaxDf,x0,0.d0,0.,iscr))
 621       if(d2.lt.1.d-10.and.x0.lt.x1)then
 622         x0=dsqrt(x0*x1)
 623 c        write(ifch,*)"droot",x0,x1,d0,d1,d2
 624         goto 5
 625         elseif(d2.gt.d0)then
 626         droot=max(x0,xfitmin)
 627 c        write(ifch,*)"droot",x0,x1,d0,d1,d2
 628         return
 629       endif
 630       f0=d2-d0
 631       f1=d1-d0
 632       if(f0*f1.lt.0.d0)then
 633
 634
 635  10   x2=dsqrt(x0*x1)
 636       d2=dabs(Dsoftshval(real(x2)*smaxDf,x2,0.d0,0.,iscr))
 637       f2=d2-d0
 638       k=k+1
 639 c        write (ifch,*) '******************* droot **************'
 640 c        write (ifch,*) x0,x1,x2,f0,f1,f2
 641
 642       if (f0*f2.lt.0.D0) then
 643         x1=x2
 644         f1=f2
 645       else
 646         x0=x2
 647         f0=f2
 648       endif
 649
 650       if (dabs((x1-x0)/x1).gt.(1.D-5).and.k.le.kmax.and.x1.ne.x0) then
 651         goto 10
 652       else
 653         if (k.gt.kmax) then
 654           write(ifch,*)'??? Warning in Droot: Delta=',dabs((x1-x0)/x1)
 655 c.........stop 'Error in Droot, too many steps'
 656         endif
 657         droot=dsqrt(x1*x0)
 658       endif
 659
 660       else
 661         droot=dsqrt(x1*x0)
 662       endif
 663
 664       return
 665       end
 666
 667 c----------------------------------------------------------------------
 668
 669       double precision function drootom(d0,dmax,bmax,eps)
 670
 671 c----------------------------------------------------------------------
 672 c Find b0 which gives f(b0)=(1-exp(-om(b0,iqq)))/dmax-d0=0
 673       include 'epos.inc'
 674       include 'epos.incsem'
 675       include 'epos.incpar'
 676       double precision om1intbc,d0,d1,d2,f0,f1,f2,dmax
 677       parameter (kmax=1000)
 678
 679
 680       k=0
 681       b0=0.
 682       b1=bmax
 683       d2=(1.d0-exp(-om1intbc(b1)))/dmax
 684       if(d2.gt.d0)then
 685         drootom=b1
 686 c        write(*,*)"drootom exit (1)",b0,b1,d0,d1,d2
 687         return
 688       endif
 689       d1=(1.d0-exp(-om1intbc(b0)))/dmax
 690       f0=d1-d0
 691       f1=d2-d0
 692       if(f0*f1.lt.0.d0)then
 693
 694
 695  10   b2=0.5*(b0+b1)
 696       d2=(1.d0-dexp(-om1intbc(b2)))/dmax
 697       f2=d2-d0
 698       k=k+1
 699 c      write (*,*) '******************* drootom **************'
 700 c      write (*,*) b0,b1,b2,f0,f1,f2
 701
 702       if (f1*f2.lt.0.D0) then
 703         b0=b2
 704         f0=f2
 705       else
 706         b1=b2
 707         f1=f2
 708       endif
 709
 710       if (abs(f2).gt.eps.and.k.le.kmax.and.b1.ne.b0) then
 711         goto 10
 712       else
 713         if (k.gt.kmax) then
 714           write(ifch,*)'??? Warning in Drootom: Delta=',abs((b1-b0)/b1)
 715 c.........stop 'Error in Droot, too many steps'
 716         endif
 717         drootom=0.5*(b1+b0)
 718       endif
 719
 720       else
 721 c        write(*,*)"drootom exit (2)",b0,b1,d0,d1,d2
 722         drootom=0.5*(b1+b0)
 723       endif
 724
 725       return
 726       end
 727
 728 c----------------------------------------------------------------------
 729       subroutine variance(r2,alp,iqq)
 730 c----------------------------------------------------------------------
 731 c fit sigma2 into : 1/sigma2(x)=1/r2-alp*log(x*s)
 732 c  iqq=0 -> soft pomeron
 733 c  iqq=1 -> semi-hard pomeron
 734 c  iqq=2 -> sum
 735 c----------------------------------------------------------------------
 736       include 'epos.inc'
 737       include 'epos.incsem'
 738       include 'epos.incpar'
 739       dimension Xs(maxdataDf),vari(maxdataDf),sigma(maxdataDf)
 740       double precision X,X0,xmax
 741
 742       do i=1,nptf
 743         sigma(i)=1.e-2
 744       enddo
 745
 746       if(iqq.eq.0.or.xshmin.gt.0.95d0)then
 747         X0=xminDf
 748         xmax=xmaxDf
 749       elseif(iqq.eq.2)then
 750         X0=xshmin
 751         xmax=xmaxDf
 752       else
 753         X0=.1d0/dlog(max(dexp(2.d0),dble(smaxDf)/1.d3))
 754         if(smaxDf.lt.100.*q2min)X0=.95d0
 755         xmax=xmaxDf
 756       endif
 757       if(iqq.ne.3.and.iqq.ne.4)then
 758
 759         do i=0,nptf-1
 760           X=X0
 761           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
 762           Xs(i+1)=log(real(X)*smaxDf)
 763           sig2=sigma2(X,iqq)
 764           if(sig2.le.0.) call utstop
 765      &   ('In variance, initial(1) sigma2 not def!&')
 766           vari(i+1)=1./sig2
 767         enddo
 768
 769 c Fit of the variance of D(X,b) between X0 and xmaxDf
 770
 771         call fit(Xs,vari,nptf,sigma,0,tr2,talp)
 772         r2=1./tr2
 773         alp=-talp
 774 c in principle, the formula to convert 1/(del+eps*log(sy)) into
 775 c  1/del*(1-eps/del*log(sy)) is valid only if eps/del*log(sy)=alp*r2*log(sy)
 776 c is small. In practice, since the fit of G(x) being an approximation, each
 777 c component of the fit should not be taken separatly but we should consider
 778 c G as one function. Then it works even with large alp (gamD).
 779 c        ttt=alp*r2*log(smaxDf)
 780 c        if(ttt.gt.0.5)
 781 c     &    write(ifmt,*)'Warning, G(b) parametrization not optimal : ',
 782 c     &          'gamD too large compared to delD !',ttt
 783       else
 784         if(iqq.eq.3)r2=sigma2(xmaxDf,3)
 785         if(iqq.eq.4)r2=sigma2(xshmin,3)
 786         if(r2.le.0.) call utstop
 787      &('In variance, initial(2) sigma2 not def!&')
 788         alp=0.
 789       endif
 790
 791       if(ish.ge.4)then
 792         write(ifch,*) '%%%%%%%%%% variance ini %%%%%%%%%%%%'
 793         write(ifch,*) 'X0=',X0
 794         write(ifch,*) 'delD ini=',r2
 795         write(ifch,*) 'gamD ini=',alp
 796       endif
 797
 798       return
 799       end
 800
 801
 802
 803 c----------------------------------------------------------------------
 804
 805         function sigma2(x,iqq)
 806
 807 c----------------------------------------------------------------------
 808 c Return the variance for a given x of :
 809 c For G :
 810 c iqq=0 the soft pomeron
 811 c iqq=1 the semi-hard and valence quark pomeron
 812 c iqq=2 the sum
 813 c----------------------------------------------------------------------
 814
 815       include 'epos.inc'
 816       include 'epos.incpar'
 817       double precision x,Dsoftshval,sfsh,om51p,eps,range,sig2!,omNpuncut
 818       external varifit
 819       double precision varifit,Db(maxdataDf),bf(maxdataDf)
 820
 821       bmax=bmaxDf
 822       sig2=bmax*0.5
 823       bmin=-bmax
 824       eps=1.d-10
 825       ierro=0
 826
 827       if(iqq.eq.0)then
 828         range=sig2
 829         sfsh=om51p(real(x)*smaxDf,x,0.d0,0.,0)
 830         if (dabs(sfsh).gt.eps) then
 831         do i=0,nptf-1
 832           bf(i+1)=dble(bmin+real(i)*(bmax-bmin)/real(nptf-1))
 833           Db(i+1)=om51p(real(x)*smaxDf,x,0.d0,real(bf(i+1)),0)/sfsh
 834         enddo
 835         else
 836           ierro=1
 837         endif
 838       elseif(iqq.eq.1.and.xshmin.lt..95d0)then
 839         range=sig2
 840         sfsh=0.d0
 841         do j=1,4
 842           sfsh=sfsh+om51p(real(x)*smaxDf,x,0.d0,0.,j)
 843         enddo
 844         if (dabs(sfsh).gt.eps) then
 845         do i=0,nptf-1
 846           bf(i+1)=dble(bmin+real(i)*(bmax-bmin)/real(nptf-1))
 847           Db(i+1)=0.d0
 848           do j=1,4
 849            Db(i+1)=Db(i+1)+om51p(real(x)*smaxDf,x,0.d0,real(bf(i+1)),j)
 850           enddo
 851           Db(i+1)=Db(i+1)/sfsh
 852         enddo
 853         else
 854           ierro=1
 855         endif
 856       else
 857         sig2=2.d0*sig2
 858         range=sig2
 859         iscr=0
 860         sfsh=Dsoftshval(real(x)*smaxDf,x,0.d0,0.,iscr)
 861         if (dabs(sfsh).gt.eps) then
 862         do i=0,nptf-1
 863           bf(i+1)=dble(bmin+real(i)*(bmax-bmin)/real(nptf-1))
 864           Db(i+1)=Dsoftshval(real(x)*smaxDf,x,0.d0,real(bf(i+1)),iscr)
 865      &              /sfsh
 866         enddo
 867         else
 868           ierro=1
 869         endif
 870       endif
 871
 872 c Fit of D(X,b) between -bmaxDf and bmaxDf
 873
 874       if(ierro.ne.1)then
 875         nptft=nptf
 876         call minfit(varifit,bf,Db,nptft,sig2,range)
 877         sigma2=real(sig2)
 878       else
 879         sigma2=0.
 880       endif
 881
 882       return
 883       end
 884
 885 c----------------------------------------------------------------------
 886
 887         subroutine paramx
 888
 889 c----------------------------------------------------------------------
 890 c updates the 4 parameters alpsf,betsf,alpsh,betsh by fitting GFF
 891 c  parDf(1,3) parDf(2,3) ... alp, bet soft
 892 c  parDf(3,3) parDf(4,3) ... alp, bet semihard
 893 c----------------------------------------------------------------------
 894
 895       include 'epos.inc'
 896       include 'epos.incpar'
 897
 898       double precision Dsoftshpar
 899
 900       external Dsoftshpar
 901
 902       dimension range(nbpf)
 903
 904       call givedatax
 905
 906       !determine parameter range
 907       do i=numminDf,numparDf
 908         range(i)=fparDf*parDf(i,3)
 909         parDf(i,1)=parDf(i,3)-range(i)
 910         parDf(i,2)=parDf(i,3)+range(i)
 911       enddo
 912
 913
 914    !   write(ifch,*) '%%%%%%%%%%%%%%%%%%% fitx %%%%%%%%%%%%%%%%%%%%%%%'
 915
 916       call fitx(Dsoftshpar,nmcxDf,chi2,err)
 917
 918    !   write(ifch,*) 'chi2=',chi2
 919    !   write(ifch,*) 'err=',err
 920    !   write(ifch,*) 'alpD(1)=',parDf(1,3),' betD(1)=',parDf(2,3)
 921    !   write(ifch,*) 'alpD(2)=',parDf(3,3),' betD(2)=',parDf(4,3)
 922
 923       return
 924       end
 925
 926
 927 c----------------------------------------------------------------------
 928       subroutine givedatax
 929 c----------------------------------------------------------------------
 930
 931       include 'epos.inc'
 932       include 'epos.incsem'
 933       include 'epos.incpar'
 934       double precision X,X0,X1,Dsoftshval,Xseuil
 935
 936       numdataDf=nptf
 937
 938       X0=max(xminDf,dble(sfshlim/smaxDf))
 939       X1=xmaxDf
 940       Xseuil=.0001d0
 941
 942 c Fit of G(X) between X0 and X1
 943       do i=0,nptf-1
 944         X=X0
 945         if (i.ne.0) X=X*(X1/X0)**(dble(i)/dble(nptf-1))
 946         datafitD(i+1,2)=max(1.e-10,
 947      &                  real(Dsoftshval(real(X)*smaxDf,X,0.d0,0.,1)))
 948         datafitD(i+1,1)=real(X)
 949         datafitD(i+1,3)=1.
 950         if (X.gt.Xseuil) datafitD(i+1,3)=exp((Xseuil/X)-1.)
 951       enddo
 952
 953       return
 954       end
 955
 956
 957
 958
 959
 960 c----------------------------------------------------------------------
 961
 962       function sigma1i(x)
 963
 964 c----------------------------------------------------------------------
 965 c Return the variance of the sum of the soft pomeron and the semi-hard
 966 c pomeron for a given x.
 967 c----------------------------------------------------------------------
 968
 969       include 'epos.inc'
 970       include 'epos.incpar'
 971       double precision x,Dsoftshval,Dint
 972
 973
 974       iscr=0
 975       Dint=Dsoftshval(real(x)*smaxDf,x,0.d0,0.,iscr)
 976
 977       sigma1i=0.
 978       if(Dint.ne.0.)
 979      &sigma1i=real(-1.d0/dlog(Dsoftshval(real(x)*smaxDf,x,0.d0,1.,iscr)
 980      &   /Dint))
 981
 982
 983       return
 984       end
 985
 986
 987 c----------------------------------------------------------------------
 988
 989       SUBROUTINE minfit(func,x,y,ndata,a,range)
 990
 991 c----------------------------------------------------------------------
 992 c Given a set of data points x(1:ndata),y(1:ndata), and the range of
 993 c the parameter a, fit it on function func by minimizing chi2.
 994 c In input a define the expected value of a, and on output they
 995 c correspond to  the fited value.
 996 c ---------------------------------------------------------------------
 997       include 'epos.inc'
 998       double precision x(ndata),y(ndata),func,a,range,Smin,Som,a1,a2,eps
 999      *,amin,rr,yp
1000       parameter (eps=1.d-5)
1001       external func
1002
1003
1004       Smin=1.d20
1005       amin=a
1006
1007
1008
1009       a1=a-range
1010       a2=a+range
1011
1012       do j=1,2000
1013         rr=dble(rangen())
1014         a=a1+(a2-a1)*rr
1015         k=0
1016
1017  10     if(a.lt.0.d0.and.k.lt.100) then
1018           rr=dble(rangen())
1019           a=a1+(a2-a1)*rr
1020           k=k+1
1021           goto 10
1022         endif
1023         if(k.ge.100) call utstop
1024      &('Always negative variance in minfit ...&')
1025
1026         Som=0.d0
1027         do k=1,ndata
1028              yp=min(1.d10,func(x(k),a))  !proposal function
1029               Som=Som+(yp-y(k))**2.d0
1030         enddo
1031         if(Som.lt.Smin)then
1032
1033           if(Smin.lt.1.)then
1034             if(a.gt.amin)then
1035               a1=amin
1036             else
1037               a2=amin
1038             endif
1039           endif
1040           amin=a
1041           Smin=Som
1042         endif
1043         if(Smin.lt.eps)goto 20
1044       enddo
1045
1046  20   continue
1047       a=amin
1048
1049       return
1050       end
1051
1052
1053
1054 c----------------------------------------------------------------------
1055       subroutine fitx(func,nmc,chi2,err)
1056 c----------------------------------------------------------------------
1057 c  Determines parameters of the funcion func
1058 c  representing the best fit of the data.
1059 c  At the end of the run, the "best" parameters are stored in parDf(n,3).
1060 c  The function func has to be defined via "function" using the parameters
1061 c  parDf(n,3), n=1,numparDf .
1062 c  Parameters as well as data are stored on /fitpar/:
1063 c    numparDf: number of parameters  (input)
1064 c    parDf: array containing parameters:
1065 c         parDf(n,1): lower limit    (input)
1066 c         parDf(n,2): upper limit    (input)
1067 c         parDf(n,3): current parameter (internal and output = final result)
1068 c         parDf(n,4): previous parameter (internal)
1069 c    numdataDf: number of data points  (input)
1070 c    datafitD: array containing data:
1071 c         datafitD(i,1): x value       (input)
1072 c         datafitD(i,2): y value       (input)
1073 c         datafitD(i,3): error         (input)
1074 c----------------------------------------------------------------------
1075
1076       include 'epos.inc'
1077       include 'epos.incpar'
1078       double precision func,x
1079       external func
1080
1081  !     write (ifch,*) 'numparDf,numminDf',numparDf,numminDf
1082
1083
1084 c initial configuration (better if one start directly with initial one)
1085
1086 c      do n=numminDf,numparDf
1087 c              parDf(n,3)=parDf(n,1)+rangen()*(parDf(n,2)-parDf(n,1))
1088 c      enddo
1089
1090       chi2=0.
1091       err=0.
1092       do i=1,numdataDf
1093         x=dble(datafitD(i,1))
1094         fx=real(func(x))
1095         chi2=chi2+(log(fx)-log(datafitD(i,2)))**2/datafitD(i,3)**2
1096         err=err+(log(fx)-log(datafitD(i,2)))/datafitD(i,3)**2
1097       enddo
1098       err=abs(err)/real(numdataDf)
1099
1100 c metropolis iteration
1101
1102       do i=1,nmc
1103 c        if(mod(i,int(real(nmc)/1000.)).eq.0)then
1104           betac=betac*(1.+1./real(nmc))!1.05
1105           betae=betae*(1.+1./real(nmc))!1.05
1106 c        endif
1107 c        if(mod(i,int(real(nmc)/20.)).eq.0)write(ifch,*)i,chi2,err
1108
1109         do n=numminDf,numparDf
1110           parDf(n,4)=parDf(n,3)
1111         enddo
1112         chi2x=chi2
1113         errx=err
1114
1115         n=numminDf+int(rangen()*(numparDf-numminDf+1))
1116         n=max(n,numminDf)
1117         n=min(n,numparDf)
1118 c              if(mod(i,int(real(nmc)/20.)).eq.0)write(ifch,*)n
1119
1120  10     parDf(n,3)=parDf(n,1)+rangen()*(parDf(n,2)-parDf(n,1))
1121         chi2=0
1122         err=0
1123         do j=1,numdataDf
1124           x=dble(datafitD(j,1))
1125           fx=real(func(x))
1126           chi2=chi2+(log(fx)-log(datafitD(j,2)))**2/datafitD(j,3)**2
1127           err=err+(log(fx)-log(datafitD(j,2)))/datafitD(j,3)**2
1128         enddo
1129         err=abs(err)/real(numdataDf)
1130
1131         if(chi2.gt.chi2x.and.rangen()
1132      $             .gt.exp(-min(50.,max(-50.,betac*(chi2-chi2x))))
1133      &             .or.err.gt.errx.and.rangen()
1134      $             .gt.exp(-min(50.,max(-50.,betae*(err-errx))))
1135      &                                                     ) then
1136           do n=numminDf,numparDf
1137             parDf(n,3)=parDf(n,4)
1138           enddo
1139           chi2=chi2x
1140           err=errx
1141         endif
1142
1143       enddo
1144
1145       return
1146       end
1147
1148
1149 c----------------------------------------------------------------------
1150
1151         SUBROUTINE fit(x,y,ndata,sig,mwt,a,b)
1152
1153 c----------------------------------------------------------------------
1154 c Given a set of data points x(1:ndata),y(1:ndata) with individual standard
1155 c deviations sig(1:ndata), fit them to a straight line y = a + bx by
1156 c minimizing chi2 .
1157 c Returned are a,b and their respective probable uncertainties siga and sigb,
1158 c the chisquare chi2, and the goodness-of-fit probability q (that the fit
1159 c would have chi2 this large or larger). If mwt=0 on input, then the standard
1160 c deviations are assumed to be unavailable: q is returned as 1.0 and the
1161 c normalization of chi2 is to unit standard deviation on all points.
1162 c ---------------------------------------------------------------------
1163
1164         implicit none
1165         INTEGER mwt,ndata
1166         REAL sig(ndata),x(ndata),y(ndata)
1167         REAL a,b,siga,sigb,chi2 !,q
1168         INTEGER i
1169         REAL sigdat,ss,st2,sx,sxoss,sy,t,wt
1170
1171
1172         sx=0.                 !Initialize sums to zero.
1173         sy=0.
1174         st2=0.
1175         b=0.
1176         if(mwt.ne.0) then ! Accumulate sums ...
1177           ss=0.
1178           do i=1,ndata          !...with weights
1179             wt=1./(sig(i)**2)
1180             ss=ss+wt
1181             sx=sx+x(i)*wt
1182             sy=sy+y(i)*wt
1183           enddo
1184         else
1185           do i=1,ndata          !...or without weights.
1186             sx=sx+x(i)
1187             sy=sy+y(i)
1188           enddo
1189           ss=float(ndata)
1190         endif
1191         sxoss=sx/ss
1192         if(mwt.ne.0) then
1193           do i=1,ndata
1194             t=(x(i)-sxoss)/sig(i)
1195             st2=st2+t*t
1196             b=b+t*y(i)/sig(i)
1197           enddo
1198         else
1199           do i=1,ndata
1200             t=x(i)-sxoss
1201             st2=st2+t*t
1202             b=b+t*y(i)
1203           enddo
1204         endif
1205         b=b/st2                 !Solve for a, b, oe a , and oe b .
1206         a=(sy-sx*b)/ss
1207         siga=sqrt((1.+sx*sx/(ss*st2))/ss)
1208         sigb=sqrt(1./st2)
1209         chi2=0.                 !Calculate chi2 .
1210 c        q=1.
1211         if(mwt.eq.0) then
1212           do i=1,ndata
1213             chi2=chi2+(y(i)-a-b*x(i))**2
1214           enddo
1215
1216 c For unweighted data evaluate typical sig using chi2, and adjust
1217 c the standard deviations.
1218
1219           sigdat=sqrt(chi2/(ndata-2))
1220           siga=siga*sigdat
1221           sigb=sigb*sigdat
1222         else
1223           do i=1,ndata
1224             chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
1225           enddo
1226         endif
1227
1228         if(chi2.ge.0.2)then
1229           b=(y(ndata)-y(1))/(x(ndata)-x(1))
1230           a=y(ndata)-b*x(ndata)
1231         endif
1232
1233
1234 c        write(ifch,*) '$$$$$$$$$$ fit : a,b,siga,sigb,chi2,q $$$$$$$$$$$'
1235 c        write(ifch,*) a,b,siga,sigb,chi2!???????????????
1236
1237
1238         return
1239         END
1240
1241
1242
1243 c----------------------------------------------------------------------
1244
1245       double precision function varifit(x,var)
1246
1247 c----------------------------------------------------------------------
1248       double precision x,var
1249
1250       varifit=dexp(-min(50.d0,x**2.d0/var))
1251
1252       return
1253       end
1254
1255
1256
1257 c----------------------------------------------------------------------
1258
1259       double precision function Dsoftshval(sy,x,y,b,iscr)
1260
1261 c----------------------------------------------------------------------
1262 c iscr=-1 sum of om5p (i), i from 0 to 4 - fit of hard
1263 c iscr=0 sum of om5p (i), i from 0 to 4
1264 c iscr=1 sum of om5p (i), i from 0 to 4 * F * F
1265 c iscr=2 sum of om5p (i), i from 1 to 4 (semihard + valence quark)
1266 c----------------------------------------------------------------------
1267       double precision x,om51p,y,xp,xm
1268       include 'epos.inc'
1269       include 'epos.incsem'
1270       include 'epos.incpar'
1271
1272       Dsoftshval=0.d0
1273
1274       if(iscr.le.0)then
1275         do i=0,4
1276           Dsoftshval=Dsoftshval+om51p(sy,x,y,b,i)
1277         enddo
1278
1279
1280       elseif(iscr.eq.1)then
1281         xp=dsqrt(x)*dexp(y)
1282         if(dabs(xp).ge.1.d-15)then
1283           xm=x/xp
1284         else
1285           xm=1.d0
1286           write(ifch,*)'Warning in Dsoftshval in epos-par'
1287         endif
1288         do i=0,4
1289           Dsoftshval=Dsoftshval+om51p(sy,x,y,b,i)
1290         enddo
1291         Dsoftshval=Dsoftshval*(1.d0-xm)**dble(alplea(icltar))
1292      &                       *(1.d0-xp)**dble(alplea(iclpro))
1293       elseif(iscr.eq.2)then
1294         do i=1,4
1295           Dsoftshval=Dsoftshval+om51p(sy,x,y,b,i)
1296         enddo
1297       endif
1298
1299
1300       Dsoftshval=2.d0*Dsoftshval
1301      &     /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
1302
1303       if(iscr.eq.-1.and.parDf(3,3).lt.parDf(1,3))Dsoftshval=Dsoftshval
1304      &     -dble(parDf(3,3)*sy**parDf(4,3))
1305
1306       return
1307       end
1308
1309 c----------------------------------------------------------------------
1310
1311       double precision function Dsoftshpar(x)
1312
1313 c----------------------------------------------------------------------
1314       double precision x,xp,xm
1315       include 'epos.inc'
1316       include 'epos.incpar'
1317
1318       Dsoftshpar=dble(
1319      &        parDf(1,3)*(real(x)*smaxDf)**parDf(2,3)
1320      &       +parDf(3,3)*(real(x)*smaxDf)**parDf(4,3) )
1321       xp=dsqrt(x)
1322       xm=xp
1323       Dsoftshpar=Dsoftshpar*(1.d0-xm)**dble(alplea(icltar))
1324      &                       *(1.d0-xp)**dble(alplea(iclpro))
1325       Dsoftshpar=min(max(1.d-15,Dsoftshpar),1.d15)
1326
1327       return
1328       end