Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:48:50

0001 c 15/02/2005 epos 1.03
0002 
0003 
0004 c----------------------------------------------------------------------
0005       subroutine paramini(imod)
0006 c----------------------------------------------------------------------
0007 c  Set  parameters of the parametrisation of the eikonals.
0008 c
0009 c xDfit=Sum(i=0,1)(alpD(i)*xp**betDp(i)*xm**betDpp(i)*s**betD(i)
0010 c                            *xs**(gamD(i)*b2)*exp(-b2/delD(i))
0011 c
0012 c Parameters stored in /Dparam/ (epos.inc)
0013 c if imod=0, do settings only for iclpro, if imod=1, do settings
0014 c for iclpro=2 and iclpro
0015 c----------------------------------------------------------------------
0016 
0017       include 'epos.inc'
0018       include 'epos.incems'
0019       include 'epos.incsem'
0020       include 'epos.incpar'
0021       double precision PhiExact,y!,Dsoftshval
0022       call utpri('parini',ish,ishini,3)
0023 
0024 c Initialisation of the variables
0025 
0026       call Class('paramini  ')
0027 
0028 c Variables used for xparg (only graphical application)
0029 
0030       spmin=4.*q2min         !Definition of spmin in psvin (epos-sha)
0031       sfshlim=3.*spmin          !transition energy for soft->hard
0032       emaxDf=engy            !energy for fit subroutines
0033       smaxDf=emaxDf**2       !energy squared for fit subroutines
0034 c      sfshlim=100.
0035       nptf=10                !number of point for the fit
0036       xmaxDf=1.d0            !maximum limit for the fit
0037       xminDf=1.d0/dble(smaxDf)   !minimum limit
0038       xggfit=xminDf*sfshlim
0039       xshmin=3.d-2           !minimum limit for sh variance fit
0040       if(smaxDf.lt.sfshlim) xshmin=1.d0
0041       xfitmin=0.1d0   !sh fitted between f(1) and f(1)*xfitmin
0042       if(smaxDf.lt.sfshlim) xfitmin=1.d0
0043 
0044       bmaxDf=2.              !maximum b for variance fit
0045 
0046       idxDmin=idxD0          !minimum indice for parameters
0047       ntymin=ntymi           !minimum indice for diagram type
0048 
0049       ucfpro=utgam1(1.+alplea(iclpro))
0050       ucftar=utgam1(1.+alplea(icltar))
0051 
0052 
0053 
0054       iiiclegy=iclegy
0055 
0056 c for pi or K - p crosse section calculation, we need alpD, ... for
0057 c iclpro=1 or 3 and iclpro=2
0058         iiiclpro1=iclpro
0059         iiidirec=1
0060         if(imod.eq.0)then
0061           iiiclpro2=iclpro
0062         else
0063           iiiclpro2=2
0064           if(iiiclpro1.lt.iiiclpro2)iiidirec=-1
0065         endif
0066         iclprosave=iclpro
0067 
0068         do iiipro=iiiclpro2,iiiclpro1,iiidirec
0069 
0070           iclpro=iiipro
0071 
0072       if(ish.ge.4)write(ifch,*)'gamini'
0073      &          ,iscreen,iclpro,icltar,iclegy,smaxDf,sfshlim,spmin
0074 
0075       if(isetcs.le.1)then                      !if set mode, do fit
0076 
0077 c First try for fit parameters
0078 
0079 c linear fit of the 2 components of G as a function of x
0080         call pompar(alpsf,betsf,0) !soft (taking into account hard)
0081         call pompar(alpsh,betsh,1) !sh
0082 c        betsh=0.31
0083 c        alpsh0=sngl(Dsoftshval(smaxDf,1d0,0.d0,0.,2))*smaxDf**(-betsh)
0084 c        alpsh1=sngl(Dsoftshval(smaxDf,1d0,0.d0,0.,0))*smaxDf**(-betsh)
0085 c        if(alpsh0.lt.alpsf)alpsh1=alpsh0
0086 c        if(alpsh0*smaxDf**betsh.lt.alpsf*smaxDf**betsf)alpsh1=alpsh0
0087 c        if(smaxDf.gt.100.*sfshlim)then
0088 c          xfmin=1.e-2
0089 c          alpsh2=sngl(Dsoftshval(xfmin*smaxDf,dble(xfmin),0.d0,0.,0))
0090 c     &             *(xfmin*smaxDf)**(-betsh)
0091 c        else
0092 c          alpsh2=-1e10
0093 c        endif
0094 c        alpsh=max(alpsh1,alpsh2)
0095 c        alpsh=max(alpsh0,alpsh2)
0096 c Gaussian fit of the 2 components of G as a function of x and b
0097         call variance(delsf,gamsf,0)
0098         call variance(delsh,gamsh,1)
0099         gamsf=max(0.,gamsf)
0100         gamsh=max(0.,gamsh)
0101 
0102 
0103 c Fit GFF
0104 
0105 
0106 c       fit parameters
0107         numminDf=3             !minimum indice
0108         numparDf=4             !maximum indice
0109         betac=100.               !temperature for chi2
0110         betae=1.               !temperature for error
0111         fparDf=0.8             !variation amplitude for range
0112 
0113         nmcxDf=20000          !number of try for x fit
0114 
0115 c starting values
0116 
0117         parDf(1,3)=alpsf
0118         parDf(2,3)=betsf
0119         parDf(3,3)=alpsh
0120         parDf(4,3)=betsh
0121 
0122 c        if(smaxDf.ge.sfshlim)then
0123 c          call paramx           !alpD and betD
0124 
0125 c          call pompar(alpsf,betsf,-1) !soft (taking into account hard)
0126 c          parDf(1,3)=alpsf
0127 c          parDf(2,3)=max(-0.95+alppar,betsf)
0128 c          parDf(2,3)=max(0.,betsf)
0129 
0130 c        endif
0131 
0132         alpsf=parDf(1,3)
0133         betsf=parDf(2,3)
0134         alpsh=parDf(3,3)
0135         betsh=parDf(4,3)
0136 
0137 
0138       else                                     !else parameters from table (inirj)
0139         nbpsf=iDxD0
0140         if(iclegy2.gt.1)then
0141           al=1.+(log(engy)-log(egylow))/log(egyfac) !energy class
0142           i2=min(iiiclegy+1,iclegy2)
0143           i1=i2-1
0144         else
0145           i1=iclegy
0146           i2=iclegy
0147           al=float(iclegy)
0148         endif
0149         dl=al-i1
0150         dl1=max(0.,1.-dl)
0151                                 !linear interpolation
0152         alpsf=alpDs(nbpsf,i2,iclpro,icltar)*dl
0153      &       +alpDs(nbpsf,i1,iclpro,icltar)*dl1
0154         alpsh=alpDs(1,i2,iclpro,icltar)*dl
0155      &       +alpDs(1,i1,iclpro,icltar)*dl1
0156         betsf=betDs(nbpsf,i2,iclpro,icltar)*dl
0157      &       +betDs(nbpsf,i1,iclpro,icltar)*dl1
0158         betsh=betDs(1,i2,iclpro,icltar)*dl
0159      &       +betDs(1,i1,iclpro,icltar)*dl1
0160         gamsf=gamDs(nbpsf,i2,iclpro,icltar)*dl
0161      &       +gamDs(nbpsf,i1,iclpro,icltar)*dl1
0162         gamsh=gamDs(1,i2,iclpro,icltar)*dl
0163      &       +gamDs(1,i1,iclpro,icltar)*dl1
0164         delsf=delDs(nbpsf,i2,iclpro,icltar)*dl
0165      &       +delDs(nbpsf,i1,iclpro,icltar)*dl1
0166         delsh=delDs(1,i2,iclpro,icltar)*dl
0167      &       +delDs(1,i1,iclpro,icltar)*dl1
0168 
0169 
0170 c For the Plots
0171         parDf(1,3)=alpsf
0172         parDf(2,3)=betsf
0173         parDf(3,3)=alpsh
0174         parDf(4,3)=betsh
0175 
0176       endif
0177 
0178 c     if energy too small to have semi-hard interaction -> only soft diagram
0179 
0180       if(smaxDf.lt.sfshlim.and.idxD0.eq.0)then !no hard: soft+hard=soft
0181         alpsf=alpsf/2.
0182         alpsh=alpsf
0183         betsh=betsf
0184         gamsh=gamsf
0185         delsh=delsf
0186       endif
0187 
0188 
0189 
0190 c Print results
0191 
0192       if(ish.ge.4)then
0193         write(ifch,*)"parameters for iclpro:",iclpro
0194         write(ifch,*)"alp,bet,gam,del sf:",alpsf,betsf,gamsf,delsf
0195         write(ifch,*)"alp,bet,gam,del sh:",alpsh,betsh,gamsh,delsh
0196       endif
0197 
0198 
0199 c Record parameters
0200 
0201 
0202       alpD(idxD0,iclpro,icltar)=alpsf
0203       alpDp(idxD0,iclpro,icltar)=0.
0204       alpDpp(idxD0,iclpro,icltar)=0.
0205       betD(idxD0,iclpro,icltar)=betsf
0206       betDp(idxD0,iclpro,icltar)=betsf
0207       betDpp(idxD0,iclpro,icltar)=betsf
0208       gamD(idxD0,iclpro,icltar)=gamsf
0209       delD(idxD0,iclpro,icltar)=delsf
0210 
0211       alpD(1,iclpro,icltar)=alpsh
0212       alpDp(1,iclpro,icltar)=0.
0213       alpDpp(1,iclpro,icltar)=0.
0214       betD(1,iclpro,icltar)=betsh
0215       betDp(1,iclpro,icltar)=betsh
0216       betDpp(1,iclpro,icltar)=betsh
0217       gamD(1,iclpro,icltar)=gamsh
0218       delD(1,iclpro,icltar)=delsh
0219 
0220       if(iomega.lt.2.and.alpdif.ne.1.)then
0221         alpDp(2,iclpro,icltar)=0.
0222         alpDpp(2,iclpro,icltar)=0.
0223         betD(2,iclpro,icltar)=0. !max(0.,betsf)
0224         alpdifs=alpdif
0225 c        alpdif=0.99
0226         betDp(2,iclpro,icltar)=-alpdifs+alppar
0227         betDpp(2,iclpro,icltar)=-alpdifs+alppar
0228 c        alpD(2,iclpro,icltar)=(alpsf+alpsh)*wdiff(iclpro)*wdiff(icltar)
0229         alpD(2,iclpro,icltar)=wdiff(iclpro)*wdiff(icltar)
0230      &            /utgam1(1.-alpdifs)**2
0231      &            *utgam1(2.-alpdifs+alplea(iclpro))
0232      &            *utgam1(2.-alpdifs+alplea(icltar))
0233      &            /chad(iclpro)/chad(icltar)
0234 c        alpdif=alpdifs
0235         gamD(2,iclpro,icltar)=0.
0236         delD(2,iclpro,icltar)=4.*.0389*(gwidth*(r2had(iclpro)
0237      &                    +r2had(icltar))+slopoms*log(smaxDf))
0238       else
0239         alpD(2,iclpro,icltar)=0.
0240         alpDp(2,iclpro,icltar)=0.
0241         alpDpp(2,iclpro,icltar)=0.
0242         betD(2,iclpro,icltar)=0.
0243         betDp(2,iclpro,icltar)=0.
0244         betDpp(2,iclpro,icltar)=0.
0245         gamD(2,iclpro,icltar)=0.
0246         delD(2,iclpro,icltar)=1.
0247       endif
0248       if(ish.ge.4)write(ifch,*)"alp,bet,betp,del dif:"
0249      &   ,alpD(2,iclpro,icltar),betD(2,iclpro,icltar)
0250      &   ,betDp(2,iclpro,icltar),delD(2,iclpro,icltar)
0251 
0252 
0253       bmxdif(iclpro,icltar)=conbmxdif()     !important to do it before kfit,  because it's used in.
0254 c          call Kfit(-1)    !xkappafit not used : if arg=-1, set xkappafit to 1
0255       if(isetcs.le.1)then
0256         if(isetcs.eq.0)then
0257           call Kfit(-1)         !xkappafit not used : if arg=-1, set xkappafit to 1)
0258         else
0259 c          call Kfit(-1)    !xkappafit not used : if arg=-1, set xkappafit to 1)
0260           call Kfit(iclegy)
0261         endif
0262 c for plots record alpDs, betDs, etc ...
0263         alpDs(idxD0,iclegy,iclpro,icltar)=alpsf
0264         alpDps(idxD0,iclegy,iclpro,icltar)=0.
0265         alpDpps(idxD0,iclegy,iclpro,icltar)=0.
0266         betDs(idxD0,iclegy,iclpro,icltar)=betsf
0267         betDps(idxD0,iclegy,iclpro,icltar)=betsf
0268         betDpps(idxD0,iclegy,iclpro,icltar)=betsf
0269         gamDs(idxD0,iclegy,iclpro,icltar)=gamsf
0270         delDs(idxD0,iclegy,iclpro,icltar)=delsf
0271 
0272         alpDs(1,iclegy,iclpro,icltar)=alpsh
0273         alpDps(1,iclegy,iclpro,icltar)=0.
0274         alpDpps(1,iclegy,iclpro,icltar)=0.
0275         betDs(1,iclegy,iclpro,icltar)=betsh
0276         betDps(1,iclegy,iclpro,icltar)=betsh
0277         betDpps(1,iclegy,iclpro,icltar)=betsh
0278         gamDs(1,iclegy,iclpro,icltar)=gamsh
0279         delDs(1,iclegy,iclpro,icltar)=delsh
0280       endif
0281 
0282       enddo
0283 
0284       if(iclpro.ne.iclprosave)stop'problem in parini with iclpro'
0285 
0286 c initialize some variable for screening
0287       if(iscreen.ne.0)then
0288         fegypp=epscrw*fscra(engy/egyscr)
0289         b2xscr=2.*epscrp*4.*.0389
0290      &       *(r2had(iclpro)+r2had(icltar)+slopom*log(engy**2))
0291 c caculate the radius where Z is saturated at epscrx to define the bases
0292 c of nuclear shadowing
0293         satrad=0.
0294         if(fegypp.gt.0.)satrad=-b2xscr*log(epscrx/fegypp)
0295         bglaubx=zbrads*sqrt(max(0.1,satrad))
0296         zbcutx=zbcut*bglaubx
0297 c        print *,'---->',bglaubx,zbcutx
0298       endif
0299 
0300       if(ish.ge.4)then             !check PhiExact value for x=1
0301         y=Phiexact(0.,0.,1.,1.d0,1.d0,smaxDf,0.)
0302         write(ifch,*)'PhiExact=',y
0303       endif
0304 
0305       call utprix('parini',ish,ishini,3)
0306 
0307       return
0308       end
0309 
0310 c----------------------------------------------------------------------
0311       subroutine Class(text)
0312 c----------------------------------------------------------------------
0313       include 'epos.inc'
0314       include 'epos.incpar'
0315       parameter (eps=1.e-5)    !to correct for precision problem)
0316       character*10 text
0317       if(iclegy1.eq.iclegy2)then
0318         iclegy=iclegy1
0319       else
0320       iclegy=1+int( (log(engy)-log(egylow))/log(egyfac) + eps ) !energy class
0321       if(iclegy.gt.iclegy2)then
0322          write(ifch,*)'***********************************************'
0323          write(ifch,*)'Warning in ',text
0324          write(ifch,*)'Energy above the range used for the fit of D:'
0325          write(ifch,*)egylow*egyfac**(iclegy1-1),egylow*egyfac**iclegy2
0326          write(ifch,*)'***********************************************'
0327          iclegy=iclegy2
0328       endif
0329       if(iclegy.lt.iclegy1)then
0330          write(ifch,*)'***********************************************'
0331          write(ifch,*)'Warning in ',text
0332          write(ifch,*)'Energy below the range used for the fit of D:'
0333          write(ifch,*)egylow*egyfac**(iclegy1-1),egylow*egyfac**iclegy2
0334          write(ifch,*)'***********************************************'
0335          iclegy=iclegy1
0336       endif
0337       endif
0338       end
0339 
0340 c----------------------------------------------------------------------
0341       subroutine param
0342 c----------------------------------------------------------------------
0343 c  Set the parameter of the parametrisation of the eikonale.
0344 c  We group the parameters into 4 array with a dimension of idxD1(=1)
0345 c  to define xDfit (see below).
0346 c
0347 c xDfit=Sum(i,0,1)(alpD(i)*xp**betDp(i)*xm**betDpp(i)*s**betD(i)
0348 c                            *xs**(gamD(i)*b2)*exp(-b2/delD(i))
0349 c
0350 c subroutine used for tabulation.
0351 c----------------------------------------------------------------------
0352       include 'epos.inc'
0353       include 'epos.incsem'
0354       include 'epos.incpar'
0355 
0356 c Initialisation of the variables
0357 
0358       emaxDf=egyfac**(iclegy-1)*egylow
0359       smaxDf=emaxDf**2
0360 
0361       spmin=4.*q2min         !Definition of spmin in psvin (epos-sha)
0362       sfshlim=3.*spmin
0363       nptf=10
0364       xmaxDf=1.d0
0365       xminDf=1d0/dble(smaxDf)
0366       xshmin=3.d-2           !minimum limit for sh variance fit
0367       if(smaxDf.lt.sfshlim) xshmin=1.d0
0368       xfitmin=0.1d0   !sh fitted between f(1) and f(1)*xfitmin
0369       if(smaxDf.lt.sfshlim) xfitmin=1.d0
0370       bmaxDf=2.
0371 
0372 
0373       if(idxD0.ne.0.and.idxD1.ne.1) stop "* idxD0/1 are not good! *"
0374 
0375       engytmp=engy
0376       engy=emaxDf
0377 
0378 c Initialisation of the parameters
0379 
0380       do i=1,nbpf
0381         do j=1,4
0382           parDf(i,j)=1.
0383         enddo
0384       enddo
0385 
0386 c.......Calculations of the parameters
0387 
0388 c First try for fit parameters
0389 
0390 c linear fit of the 2 components of G as a function of x
0391       call pompar(alpsf,betsf,0) !soft
0392       call pompar(alpsh,betsh,1) !sh
0393 c Gaussian fit of the 2 components of G as a function of x and b
0394       call variance(delsf,gamsf,0)
0395       call variance(delsh,gamsh,1)
0396       gamsf=max(0.,gamsf)
0397       gamsh=max(0.,gamsh)
0398 
0399 c Fit GFF
0400 
0401 c      fit parameters
0402       numminDf=3                !minimum indice
0403       numparDf=4                !maximum indice
0404       betac=100.                 !temperature for chi2
0405       betae=1.                 !temperature for error
0406       fparDf=0.8                 !variation amplitude for range
0407 
0408       nmcxDf=20000              !number of try for x fit
0409 
0410 c starting values
0411 
0412       parDf(1,3)=alpsf
0413       parDf(2,3)=betsf
0414       parDf(3,3)=alpsh
0415       parDf(4,3)=betsh
0416 
0417 c      if(smaxDf.ge.3.*sfshlim)then
0418 c        call paramx             !alpD and betD
0419 c
0420 c        call pompar(alpsf,betsf,-1) !soft (taking into account hard)
0421 c        parDf(1,3)=alpsf
0422 c        parDf(2,3)=max(-0.95+alppar,betsf)
0423 c
0424 c      endif
0425 
0426       alpsf=parDf(1,3)
0427       betsf=parDf(2,3)
0428       alpsh=parDf(3,3)
0429       betsh=parDf(4,3)
0430 
0431       if(ish.ge.4)then
0432         write(ifch,*)"param: fit parameters :",iscreen,iclpro,icltar
0433      *                                        ,iclegy,engy
0434         write(ifch,*)"alp,bet,gam,del sf:",alpsf,betsf,gamsf,delsf
0435         write(ifch,*)"alp,bet,gam,del sh:",alpsh,betsh,gamsh,delsh
0436       endif
0437 
0438       alpDs(idxD0,iclegy,iclpro,icltar)=alpsf
0439       alpDps(idxD0,iclegy,iclpro,icltar)=betsf
0440       alpDpps(idxD0,iclegy,iclpro,icltar)=0.
0441       betDs(idxD0,iclegy,iclpro,icltar)=betsf
0442       betDps(idxD0,iclegy,iclpro,icltar)=betsf
0443       betDpps(idxD0,iclegy,iclpro,icltar)=betsf
0444       gamDs(idxD0,iclegy,iclpro,icltar)=gamsf
0445       delDs(idxD0,iclegy,iclpro,icltar)=delsf
0446 
0447       alpDs(1,iclegy,iclpro,icltar)=alpsh
0448       alpDps(1,iclegy,iclpro,icltar)=betsh
0449       alpDpps(1,iclegy,iclpro,icltar)=0.
0450       betDs(1,iclegy,iclpro,icltar)=betsh
0451       betDps(1,iclegy,iclpro,icltar)=betsh
0452       betDpps(1,iclegy,iclpro,icltar)=betsh
0453       gamDs(1,iclegy,iclpro,icltar)=gamsh
0454       delDs(1,iclegy,iclpro,icltar)=delsh
0455 
0456       engy=engytmp
0457 
0458       return
0459       end
0460 
0461 
0462 c----------------------------------------------------------------------
0463 
0464         subroutine pompar(alpha,beta,iqq)
0465 
0466 c----------------------------------------------------------------------
0467 c  Return the power beta and the factor alpha of the fit of the eikonal
0468 c of a pomeron of type iqq : D(X)=alpha*(X)**beta
0469 c----------------------------------------------------------------------
0470 
0471       include 'epos.inc'
0472       include 'epos.incsem'
0473       include 'epos.incpar'
0474       double precision X,D1,D0,X0,D,droot
0475       double precision Dsoftshval,xmax
0476       dimension xlnXs(maxdataDf),xlnD(maxdataDf),sigma(maxdataDf)
0477 
0478       do i=1,nptf
0479         sigma(i)=1.e-2
0480       enddo
0481 
0482       if(iqq.le.0) then
0483 
0484         iscr=iqq
0485         xmax=min(0.1d0,10.d0*xminDf)
0486         X0=xminDf
0487         if(ish.ge.4)write (ifch,*) 'pompar (0) x0,xmax=',X0,xmax
0488 
0489         do i=0,nptf-1
0490           X=X0
0491           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
0492           D=max(1.d-10,Dsoftshval(sngl(X)*smaxDf,X,0.d0,0.,iscr))
0493           if(D.eq.1.d-10)then
0494             write(ifch,*)
0495      &    "Warning in pompar ! Dsoftshval(0) could be negative"
0496             sigma(i+1)=1.e5
0497           endif
0498           xlnXs(i+1)=sngl(dlog(X*dble(smaxDf)))
0499           xlnD(i+1)=sngl(dlog(D))
0500         enddo
0501 
0502 
0503 c Fit of D(X) between X0 and xmax
0504 
0505         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
0506         if(beta.gt.10.)beta=10.
0507 
0508         alpha=exp(a)
0509 c        alpha=sngl(Dsoftshval(sngl(X0)*smaxDf,X0,0.d0,0.,iscr))
0510 c     &       *(sngl(X0)*smaxDf)**(-beta)
0511 
0512 
0513       elseif(iqq.eq.1.and.xfitmin.ne.1.d0) then
0514 
0515         iscr=2
0516 c        xmax=max(0.01d0,min(1d0,dble(sfshlim*100./smaxDf)))
0517         xmax=1d0!min(1d0,dble(sfshlim*1000./smaxDf))
0518 
0519 c Definition of D0=D(X0)
0520 
0521         D1=Dsoftshval(sngl(xmax)*smaxDf,xmax,0.d0,0.,iscr)
0522 
0523         D0=xfitmin*D1
0524 
0525 c Calculation of X0 and D(X)
0526 
0527         X0=droot(D0,D1,xmax,iscr)
0528         if(ish.ge.4)write (ifch,*) 'pompar (1) x0,xmax=',X0,xmax
0529 
0530         do i=0,nptf-1
0531           X=X0
0532           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
0533           D=max(1.d-10,Dsoftshval(sngl(X)*smaxDf,X,0.d0,0.,iscr))
0534           if(D.eq.1.d-10)then
0535             write(ifch,*)
0536      &    "Warning in pompar ! Dsoftshval(1) could be negative"
0537             sigma(i+1)=1.e5
0538           endif
0539           xlnXs(i+1)=sngl(dlog(X*dble(smaxDf)))
0540           xlnD(i+1)=sngl(dlog(D))
0541         enddo
0542 
0543 
0544 c Fit of D(X) between X0 and xmax
0545 
0546         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
0547         if(beta.gt.10.)beta=10.
0548 
0549         alpha=exp(a)
0550 c        alpha=sngl(Dsoftshval(sngl(xmax)*smaxDf,xmax,0.d0,0.,iscr))
0551 c     &       *(sngl(xmax)*smaxDf)**(-beta)
0552 
0553 
0554       elseif(iqq.eq.10.and.xfitmin.ne.1.d0) then    ! iqq=10
0555 
0556         iscr=0
0557         xmax=1.d0 !2.d0/max(2.d0,dlog(dble(smaxDf)/1.d3))
0558 
0559 c Definition of D0=D(X0)
0560 
0561         D1=Dsoftshval(sngl(xmax)*smaxDf,xmax,0.d0,0.,iscr)
0562 
0563         D0=xfitmin*D1
0564 
0565 c Calculation of X0 and D(X)
0566 
0567         X0=droot(D0,D1,xmax,iscr)
0568         if(ish.ge.4)write (ifch,*) 'pompar (1) x0,xmax=',X0,xmax
0569 
0570         do i=0,nptf-1
0571           X=X0
0572           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
0573           D=max(1.d-10,Dsoftshval(sngl(X)*smaxDf,X,0.d0,0.,iscr))
0574           if(D.eq.1.d-10)then
0575             write(ifch,*)
0576      &    "Warning in pompar ! Dsoftshval(10) could be negative"
0577             sigma(i+1)=1.e5
0578           endif
0579           xlnXs(i+1)=sngl(dlog(X*dble(smaxDf)))
0580           xlnD(i+1)=sngl(dlog(D))
0581         enddo
0582 
0583 
0584 c Fit of D(X) between X0 and xmax
0585 
0586         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
0587         if(beta.gt.10.)beta=10.
0588 
0589         alpha=sngl(Dsoftshval(sngl(xmax)*smaxDf,xmax,0.d0,0.,iscr))
0590      &       *(sngl(xmax)*smaxDf)**(-beta)
0591 
0592 
0593 
0594       else                      !iqq=-1 or iqq=1 and xfitmin=1
0595 
0596 c Calculation of X0 and D(X)
0597         iscr=0
0598 
0599         X0=1.d0/dble(smaxDf)
0600         xmax=max(2.d0/dble(smaxDf),
0601      &       min(max(0.03d0,dble(smaxDf)/2.d5),0.1d0))
0602 
0603         if(ish.ge.4)write (ifch,*) 'pompar (-1) x0,xmax=',X0,xmax
0604 
0605         do i=0,nptf-1
0606           X=X0
0607           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
0608           D=max(1.d-10,Dsoftshval(sngl(X)*smaxDf,X,0.d0,0.,iscr))
0609           if(D.eq.1.d-10)then
0610             write(ifch,*)
0611      &    "Warning in pompar ! Dsoftshval(-1) could be negative"
0612             sigma(i+1)=1.e5
0613           endif
0614           xlnXs(i+1)=sngl(dlog(X*dble(smaxDf)))
0615           xlnD(i+1)=sngl(dlog(D))
0616         enddo
0617 
0618 c Fit of D(X) between X0 and xmax
0619 
0620         call fit(xlnXs,xlnD,nptf,sigma,0,a,beta)
0621         if(beta.gt.10.)beta=10.
0622 
0623 
0624         alpha=exp(a)
0625 c        alpha=sngl(Dsoftshval(sngl(xmax)*smaxDf,xmax,0.d0,0.,iscr))
0626 c     &           *(sngl(xmax)*smaxDf)**(-beta)
0627       endif
0628 
0629         if(ish.ge.4)write(ifch,*) '%%%%%%%%%%%%% pompar %%%%%%%%%%%%%'
0630         if(ish.ge.4)write(ifch,*) 'alpD ini =',alpha,' betD ini=',beta
0631 
0632       return
0633       end
0634 
0635 c----------------------------------------------------------------------
0636 
0637       double precision function droot(d0,d1,xmax,iscr)
0638 
0639 c----------------------------------------------------------------------
0640 c Find x0 which gives f(x0)=D(x0*S)-d0=0
0641 c iqq=0 soft pomeron
0642 c iqq=1 semi-hard pomeron
0643 c----------------------------------------------------------------------
0644       include 'epos.inc'
0645       include 'epos.incsem'
0646       include 'epos.incpar'
0647       double precision Dsoftshval,d0,d1,d2,x0,x1,x2,f0,f1,f2,xmax
0648       parameter (kmax=1000)
0649 
0650 
0651       k=0
0652       x0=max(1.d0/dble(sfshlim),1d-5)
0653       x1=xmax
0654  5    d2=dabs(Dsoftshval(sngl(x0)*smaxDf,x0,0.d0,0.,iscr))
0655       if(d2.lt.1.d-10.and.x0.lt.x1)then
0656         x0=dsqrt(x0*x1)
0657 c        write(ifch,*)"droot",x0,x1,d0,d1,d2
0658         goto 5
0659         elseif(d2.gt.d0)then
0660         droot=x0
0661 c        write(ifch,*)"droot",x0,x1,d0,d1,d2
0662         return
0663       endif
0664       f0=d2-d0
0665       f1=d1-d0
0666       if(f0*f1.lt.0.d0)then
0667 
0668 
0669  10   x2=dsqrt(x0*x1)
0670       d2=dabs(Dsoftshval(sngl(x2)*smaxDf,x2,0.d0,0.,iscr))
0671       f2=d2-d0
0672       k=k+1
0673 c        write (ifch,*) '******************* droot **************'
0674 c        write (ifch,*) x0,x1,x2,f0,f1,f2
0675 
0676       if (f0*f2.lt.0.D0) then
0677         x1=x2
0678         f1=f2
0679       else
0680         x0=x2
0681         f0=f2
0682       endif
0683 
0684       if (dabs((x1-x0)/x1).gt.(1.D-5).and.k.le.kmax.and.x1.ne.x0) then
0685         goto 10
0686       else
0687         if (k.gt.kmax) then
0688           write(ifch,*)'??? Warning in Droot: Delta=',dabs((x1-x0)/x1)
0689 c.........stop 'Error in Droot, too many steps'
0690         endif
0691         droot=dsqrt(x1*x0)
0692       endif
0693 
0694       else
0695         droot=dsqrt(x1*x0)
0696       endif
0697 
0698       return
0699       end
0700 
0701 c----------------------------------------------------------------------
0702 
0703       double precision function drootom(d0,dmax,bmax,eps)
0704 
0705 c----------------------------------------------------------------------
0706 c Find b0 which gives f(b0)=(1-exp(-om(b0,iqq)))/dmax-d0=0
0707       include 'epos.inc'
0708       include 'epos.incsem'
0709       include 'epos.incpar'
0710       double precision om1intbc,d0,d1,d2,f0,f1,f2,dmax
0711       parameter (kmax=1000)
0712 
0713 
0714       k=0
0715       b0=0.
0716       b1=bmax
0717       d2=(1.d0-exp(-om1intbc(b1)))/dmax
0718       if(d2.gt.d0)then
0719         drootom=b1
0720 c        write(*,*)"drootom exit (1)",b0,b1,d0,d1,d2
0721         return
0722       endif
0723       d1=(1.d0-exp(-om1intbc(b0)))/dmax
0724       f0=d1-d0
0725       f1=d2-d0
0726       if(f0*f1.lt.0.d0)then
0727 
0728 
0729  10   b2=0.5*(b0+b1)
0730       d2=(1.d0-dexp(-om1intbc(b2)))/dmax
0731       f2=d2-d0
0732       k=k+1
0733 c      write (*,*) '******************* drootom **************'
0734 c      write (*,*) b0,b1,b2,f0,f1,f2
0735 
0736       if (f1*f2.lt.0.D0) then
0737         b0=b2
0738         f0=f2
0739       else
0740         b1=b2
0741         f1=f2
0742       endif
0743 
0744       if (abs(f2).gt.eps.and.k.le.kmax.and.b1.ne.b0) then
0745         goto 10
0746       else
0747         if (k.gt.kmax) then
0748           write(ifch,*)'??? Warning in Drootom: Delta=',abs((b1-b0)/b1)
0749 c.........stop 'Error in Droot, too many steps'
0750         endif
0751         drootom=0.5*(b1+b0)
0752       endif
0753 
0754       else
0755 c        write(*,*)"drootom exit (2)",b0,b1,d0,d1,d2
0756         drootom=0.5*(b1+b0)
0757       endif
0758 
0759       return
0760       end
0761 
0762 c----------------------------------------------------------------------
0763       subroutine variance(r2,alp,iqq)
0764 c----------------------------------------------------------------------
0765 c fit sigma2 into : 1/sigma2(x)=1/r2-alp*log(x*s)
0766 c  iqq=0 -> soft pomeron
0767 c  iqq=1 -> semi-hard pomeron
0768 c  iqq=2 -> sum
0769 c----------------------------------------------------------------------
0770       include 'epos.inc'
0771       include 'epos.incsem'
0772       include 'epos.incpar'
0773       dimension Xs(maxdataDf),vari(maxdataDf),sigma(maxdataDf)
0774       double precision X,X0,xmax
0775 
0776       do i=1,nptf
0777         sigma(i)=1.e-2
0778       enddo
0779 
0780       if(iqq.eq.0.or.xshmin.gt.0.95d0)then
0781         X0=xminDf
0782         xmax=xmaxDf
0783       elseif(iqq.eq.2)then
0784         X0=xshmin
0785         xmax=xmaxDf
0786       else
0787         X0=1d0/dlog(max(exp(2.d0),dble(smaxDf)/1.d3))
0788 c        if(smaxDf.lt.100.*q2min)X0=.95d0
0789         xmax=xmaxDf
0790       endif
0791       if(iqq.ne.3.and.iqq.ne.4)then
0792 
0793         do i=0,nptf-1
0794           X=X0
0795           if (i.ne.0) X=X*(xmax/X0)**(dble(i)/dble(nptf-1))
0796           Xs(i+1)=log(sngl(X)*smaxDf)
0797           sig2=sigma2(X,iqq)
0798           if(sig2.le.0.) call utstop
0799      &   ('In variance, initial(1) sigma2 not def!&')
0800           vari(i+1)=1./sig2
0801         enddo
0802 
0803 c Fit of the variance of D(X,b) between X0 and xmaxDf
0804 
0805         call fit(Xs,vari,nptf,sigma,0,tr2,talp)
0806         r2=1./tr2
0807         alp=-talp
0808 c in principle, the formula to convert 1/(del+eps*log(sy)) into
0809 c  1/del*(1-eps/del*log(sy)) is valid only if eps/del*log(sy)=alp*r2*log(sy)
0810 c is small. In practice, since the fit of G(x) being an approximation, each
0811 c component of the fit should not be taken separatly but we should consider
0812 c G as one function. Then it works even with large alp (gamD).
0813 c        ttt=alp*r2*log(smaxDf)
0814 c        if(ttt.gt.0.5)
0815 c     &    write(ifmt,*)'Warning, G(b) parametrization not optimal : ',
0816 c     &          'gamD too large compared to delD !',ttt
0817       else
0818         if(iqq.eq.3)r2=sigma2(xmaxDf,3)
0819         if(iqq.eq.4)r2=sigma2(xshmin,3)
0820         if(r2.le.0.) call utstop
0821      &('In variance, initial(2) sigma2 not def!&')
0822         alp=0.
0823       endif
0824 
0825       if(ish.ge.4)then
0826         write(ifch,*) '%%%%%%%%%% variance ini %%%%%%%%%%%%'
0827         write(ifch,*) 'X0=',X0
0828         write(ifch,*) 'delD ini=',r2
0829         write(ifch,*) 'gamD ini=',alp
0830       endif
0831 
0832       return
0833       end
0834 
0835 
0836 
0837 c----------------------------------------------------------------------
0838 
0839         function sigma2(x,iqq)
0840 
0841 c----------------------------------------------------------------------
0842 c Return the variance for a given x of :
0843 c For G :
0844 c iqq=0 the soft pomeron
0845 c iqq=1 the semi-hard and valence quark pomeron
0846 c iqq=2 the sum
0847 c----------------------------------------------------------------------
0848 
0849       include 'epos.inc'
0850       include 'epos.incpar'
0851       double precision x,Dsoftshval,sfsh,om51p,eps,range,sig2!,omNpuncut
0852       external varifit
0853       double precision varifit,Db(maxdataDf),bf(maxdataDf)
0854 
0855       bmax=bmaxDf
0856       sig2=bmax*0.5
0857       bmin=-bmax
0858       eps=1.d-10
0859       ierro=0
0860 
0861       if(iqq.eq.0)then
0862         range=sig2
0863         sfsh=om51p(sngl(x)*smaxDf,x,0.d0,0.,0)
0864         if (dabs(sfsh).gt.eps) then
0865         do i=0,nptf-1
0866           bf(i+1)=dble(bmin+float(i)*(bmax-bmin)/float(nptf-1))
0867           Db(i+1)=om51p(sngl(x)*smaxDf,x,0.d0,sngl(bf(i+1)),0)/sfsh
0868         enddo
0869         else
0870           ierro=1
0871         endif
0872       elseif(iqq.eq.1.and.xshmin.lt..95d0)then
0873         range=sig2
0874         sfsh=0.d0
0875         do j=1,4
0876           sfsh=sfsh+om51p(sngl(x)*smaxDf,x,0.d0,0.,j)
0877         enddo
0878         if (dabs(sfsh).gt.eps) then
0879         do i=0,nptf-1
0880           bf(i+1)=dble(bmin+float(i)*(bmax-bmin)/float(nptf-1))
0881           Db(i+1)=0.d0
0882           do j=1,4
0883            Db(i+1)=Db(i+1)+om51p(sngl(x)*smaxDf,x,0.d0,sngl(bf(i+1)),j)
0884           enddo
0885           Db(i+1)=Db(i+1)/sfsh
0886         enddo
0887         else
0888           ierro=1
0889         endif
0890       else
0891         sig2=2.d0*sig2
0892         range=sig2
0893         iscr=0
0894         sfsh=Dsoftshval(sngl(x)*smaxDf,x,0.d0,0.,iscr)
0895         if (dabs(sfsh).gt.eps) then
0896         do i=0,nptf-1
0897           bf(i+1)=dble(bmin+float(i)*(bmax-bmin)/float(nptf-1))
0898           Db(i+1)=Dsoftshval(sngl(x)*smaxDf,x,0.d0,sngl(bf(i+1)),iscr)
0899      &              /sfsh
0900         enddo
0901         else
0902           ierro=1
0903         endif
0904       endif
0905 
0906 c Fit of D(X,b) between -bmaxDf and bmaxDf
0907 
0908       if(ierro.ne.1)then
0909         nptft=nptf
0910         call minfit(varifit,bf,Db,nptft,sig2,range)
0911         sigma2=sngl(sig2)
0912       else
0913         sigma2=0.
0914       endif
0915 
0916       return
0917       end
0918 
0919 c----------------------------------------------------------------------
0920 
0921         subroutine paramx
0922 
0923 c----------------------------------------------------------------------
0924 c updates the 4 parameters alpsf,betsf,alpsh,betsh by fitting GFF
0925 c  parDf(1,3) parDf(2,3) ... alp, bet soft
0926 c  parDf(3,3) parDf(4,3) ... alp, bet semihard
0927 c----------------------------------------------------------------------
0928 
0929       include 'epos.inc'
0930       include 'epos.incpar'
0931 
0932       double precision Dsoftshpar
0933 
0934       external Dsoftshpar
0935 
0936       dimension range(nbpf)
0937 
0938       call givedatax
0939 
0940       !determine parameter range
0941       do i=numminDf,numparDf
0942         range(i)=fparDf*parDf(i,3)
0943         parDf(i,1)=parDf(i,3)-range(i)
0944         parDf(i,2)=parDf(i,3)+range(i)
0945       enddo
0946 
0947 
0948    !   write(ifch,*) '%%%%%%%%%%%%%%%%%%% fitx %%%%%%%%%%%%%%%%%%%%%%%'
0949 
0950       call fitx(Dsoftshpar,nmcxDf,chi2,err)
0951 
0952    !   write(ifch,*) 'chi2=',chi2
0953    !   write(ifch,*) 'err=',err
0954    !   write(ifch,*) 'alpD(1)=',parDf(1,3),' betD(1)=',parDf(2,3)
0955    !   write(ifch,*) 'alpD(2)=',parDf(3,3),' betD(2)=',parDf(4,3)
0956 
0957       return
0958       end
0959 
0960 
0961 c----------------------------------------------------------------------
0962       subroutine givedatax
0963 c----------------------------------------------------------------------
0964 
0965       include 'epos.inc'
0966       include 'epos.incsem'
0967       include 'epos.incpar'
0968       double precision X,X0,X1,Dsoftshval,Xseuil
0969 
0970       numdataDf=nptf
0971 
0972       X0=xminDf
0973       X1=xmaxDf
0974       Xseuil=1d0 !min(1.d0,dble(sfshlim*1e4)*XminDf)
0975 c      print *,'--------------->',Xseuil
0976 
0977 c Fit of G(X) between X0 and X1
0978       do i=0,nptf-1
0979         X=X0
0980         if (i.ne.0) X=X*(X1/X0)**(dble(i)/dble(nptf-1))
0981         datafitD(i+1,2)=max(1.e-10,
0982      &                  sngl(Dsoftshval(sngl(X)*smaxDf,X,0.d0,0.,1)))
0983         datafitD(i+1,1)=sngl(X)
0984         datafitD(i+1,3)=1.
0985         if (X.gt.Xseuil)
0986      &  datafitD(i+1,3)=exp(-min((sngl(Xseuil/X)-1.),50.))
0987       enddo
0988 
0989       return
0990       end
0991 
0992 
0993 
0994 
0995 
0996 c----------------------------------------------------------------------
0997 
0998       function sigma1i(x)
0999 
1000 c----------------------------------------------------------------------
1001 c Return the variance of the sum of the soft pomeron and the semi-hard
1002 c pomeron for a given x.
1003 c----------------------------------------------------------------------
1004 
1005       include 'epos.inc'
1006       include 'epos.incpar'
1007       double precision x,Dsoftshval,Dint
1008 
1009 
1010       iscr=0
1011       Dint=Dsoftshval(sngl(x)*smaxDf,x,0.d0,0.,iscr)
1012 
1013       sigma1i=0.
1014       if(Dint.ne.0.)
1015      &sigma1i=sngl(-1.d0/dlog(Dsoftshval(sngl(x)*smaxDf,x,0.d0,1.,iscr)
1016      &   /Dint))
1017 
1018 
1019       return
1020       end
1021 
1022 
1023 c----------------------------------------------------------------------
1024 
1025       SUBROUTINE minfit(func,x,y,ndata,a,range)
1026 
1027 c----------------------------------------------------------------------
1028 c Given a set of data points x(1:ndata),y(1:ndata), and the range of
1029 c the parameter a, fit it on function func by minimizing chi2.
1030 c In input a define the expected value of a, and on output they
1031 c correspond to  the fited value.
1032 c ---------------------------------------------------------------------
1033       include 'epos.inc'
1034       double precision x(ndata),y(ndata),func,a,range,Smin,Som,a1,a2,eps
1035      *,amin,rr,yp
1036       parameter (eps=1.d-5)
1037       external func
1038 
1039 
1040       Smin=1.d20
1041       amin=a
1042 
1043 
1044 
1045       a1=a-range
1046       a2=a+range
1047 
1048       do j=1,2000
1049         rr=dble(rangen())
1050         a=a1+(a2-a1)*rr
1051         k=0
1052 
1053  10     if(a.lt.0.d0.and.k.lt.100) then
1054           rr=dble(rangen())
1055           a=a1+(a2-a1)*rr
1056           k=k+1
1057           goto 10
1058         endif
1059         if(k.ge.100) call utstop
1060      &('Always negative variance in minfit ...&')
1061 
1062         Som=0.d0
1063         do k=1,ndata
1064              yp=min(1.d10,func(x(k),a))  !proposal function
1065               Som=Som+(yp-y(k))**2.d0
1066         enddo
1067         if(Som.lt.Smin)then
1068 
1069           if(Smin.lt.1.)then
1070             if(a.gt.amin)then
1071               a1=amin
1072             else
1073               a2=amin
1074             endif
1075           endif
1076           amin=a
1077           Smin=Som
1078         endif
1079         if(Smin.lt.eps)goto 20
1080       enddo
1081 
1082  20   continue
1083       a=amin
1084 
1085       return
1086       end
1087 
1088 
1089 
1090 c----------------------------------------------------------------------
1091       subroutine fitx(func,nmc,chi2,err)
1092 c----------------------------------------------------------------------
1093 c  Determines parameters of the funcion func
1094 c  representing the best fit of the data.
1095 c  At the end of the run, the "best" parameters are stored in parDf(n,3).
1096 c  The function func has to be defined via "function" using the parameters
1097 c  parDf(n,3), n=1,numparDf .
1098 c  Parameters as well as data are stored on /fitpar/:
1099 c    numparDf: number of parameters  (input)
1100 c    parDf: array containing parameters:
1101 c         parDf(n,1): lower limit    (input)
1102 c         parDf(n,2): upper limit    (input)
1103 c         parDf(n,3): current parameter (internal and output = final result)
1104 c         parDf(n,4): previous parameter (internal)
1105 c    numdataDf: number of data points  (input)
1106 c    datafitD: array containing data:
1107 c         datafitD(i,1): x value       (input)
1108 c         datafitD(i,2): y value       (input)
1109 c         datafitD(i,3): error         (input)
1110 c----------------------------------------------------------------------
1111 
1112       include 'epos.inc'
1113       include 'epos.incpar'
1114       double precision func,x
1115       external func
1116 
1117  !     write (ifch,*) 'numparDf,numminDf',numparDf,numminDf
1118 
1119 
1120 c initial configuration (better if one start directly with initial one)
1121 
1122 c      do n=numminDf,numparDf
1123 c              parDf(n,3)=parDf(n,1)+rangen()*(parDf(n,2)-parDf(n,1))
1124 c      enddo
1125 
1126       chi2=0.
1127       err=0.
1128       do i=1,numdataDf
1129         x=dble(datafitD(i,1))
1130         fx=sngl(func(x))
1131         chi2=chi2+(log(fx)-log(datafitD(i,2)))**2/datafitD(i,3)**2
1132         err=err+(log(fx)-log(datafitD(i,2)))/datafitD(i,3)**2
1133       enddo
1134       err=abs(err)/float(numdataDf)
1135 
1136 c metropolis iteration
1137 
1138       do i=1,nmc
1139 c        if(mod(i,int(float(nmc)/1000.)).eq.0)then
1140           betac=betac*(1.+1./float(nmc))!1.05
1141           betae=betae*(1.+1./float(nmc))!1.05
1142 c        endif
1143 c        if(mod(i,int(float(nmc)/20.)).eq.0)write(ifch,*)i,chi2,err
1144 
1145         do n=numminDf,numparDf
1146           parDf(n,4)=parDf(n,3)
1147         enddo
1148         chi2x=chi2
1149         errx=err
1150 
1151         n=numminDf+int(rangen()*(numparDf-numminDf+1))
1152         n=max(n,numminDf)
1153         n=min(n,numparDf)
1154 c              if(mod(i,int(float(nmc)/20.)).eq.0)write(ifch,*)n
1155 
1156         parDf(n,3)=parDf(n,1)+rangen()*(parDf(n,2)-parDf(n,1))
1157         chi2=0
1158         err=0
1159         do j=1,numdataDf
1160           x=dble(datafitD(j,1))
1161           fx=sngl(func(x))
1162           chi2=chi2+(log(fx)-log(datafitD(j,2)))**2/datafitD(j,3)**2
1163           err=err+(log(fx)-log(datafitD(j,2)))/datafitD(j,3)**2
1164         enddo
1165         err=abs(err)/float(numdataDf)
1166 
1167         if(chi2.gt.chi2x.and.rangen()
1168      $             .gt.exp(-min(50.,max(-50.,betac*(chi2-chi2x))))
1169      &             .or.err.gt.errx.and.rangen()
1170      $             .gt.exp(-min(50.,max(-50.,betae*(err-errx))))
1171      &                                                     ) then
1172           do n=numminDf,numparDf
1173             parDf(n,3)=parDf(n,4)
1174           enddo
1175           chi2=chi2x
1176           err=errx
1177         endif
1178 
1179       enddo
1180 
1181       return
1182       end
1183 
1184 
1185 c----------------------------------------------------------------------
1186 
1187         SUBROUTINE fit(x,y,ndata,sig,mwt,a,b)
1188 
1189 c----------------------------------------------------------------------
1190 c Given a set of data points x(1:ndata),y(1:ndata) with individual standard
1191 c deviations sig(1:ndata), fit them to a straight line y = a + bx by
1192 c minimizing chi2 .
1193 c Returned are a,b and their respective probable uncertainties siga and sigb,
1194 c the chi­square chi2, and the goodness-of-fit probability q (that the fit
1195 c would have chi2 this large or larger). If mwt=0 on input, then the standard
1196 c deviations are assumed to be unavailable: q is returned as 1.0 and the
1197 c normalization of chi2 is to unit standard deviation on all points.
1198 c ---------------------------------------------------------------------
1199 
1200         implicit none
1201         INTEGER mwt,ndata
1202         REAL sig(ndata),x(ndata),y(ndata)
1203         REAL a,b,siga,sigb,chi2 !,q
1204         INTEGER i
1205         REAL sigdat,ss,st2,sx,sxoss,sy,t,wt
1206 
1207 
1208         sx=0.                 !Initialize sums to zero.
1209         sy=0.
1210         st2=0.
1211         b=0.
1212         if(mwt.ne.0) then ! Accumulate sums ...
1213           ss=0.
1214           do i=1,ndata          !...with weights
1215             wt=1./(sig(i)**2)
1216             ss=ss+wt
1217             sx=sx+x(i)*wt
1218             sy=sy+y(i)*wt
1219           enddo
1220         else
1221           do i=1,ndata          !...or without weights.
1222             sx=sx+x(i)
1223             sy=sy+y(i)
1224           enddo
1225           ss=float(ndata)
1226         endif
1227         sxoss=sx/ss
1228         if(mwt.ne.0) then
1229           do i=1,ndata
1230             t=(x(i)-sxoss)/sig(i)
1231             st2=st2+t*t
1232             b=b+t*y(i)/sig(i)
1233           enddo
1234         else
1235           do i=1,ndata
1236             t=x(i)-sxoss
1237             st2=st2+t*t
1238             b=b+t*y(i)
1239           enddo
1240         endif
1241         b=b/st2                 !Solve for a, b, oe a , and oe b .
1242         a=(sy-sx*b)/ss
1243         siga=sqrt((1.+sx*sx/(ss*st2))/ss)
1244         sigb=sqrt(1./st2)
1245         chi2=0.                 !Calculate chi2 .
1246 c        q=1.
1247         if(mwt.eq.0) then
1248           do i=1,ndata
1249             chi2=chi2+(y(i)-a-b*x(i))**2
1250           enddo
1251 
1252 c For unweighted data evaluate typical sig using chi2, and adjust
1253 c the standard deviations.
1254 
1255           sigdat=sqrt(chi2/(ndata-2))
1256           siga=siga*sigdat
1257           sigb=sigb*sigdat
1258         else
1259           do i=1,ndata
1260             chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2
1261           enddo
1262         endif
1263 
1264         if(chi2.ge.0.2)then
1265           b=(y(ndata)-y(1))/(x(ndata)-x(1))
1266           a=y(ndata)-b*x(ndata)
1267         endif
1268 
1269 
1270 c        write(*,*) '$$$$$$$$ fit : a,b,siga,sigb,chi2,q $$$$$$$$$'
1271 c        write(*,*) a,b,siga,sigb,chi2!???????????????
1272 
1273 
1274         return
1275         END
1276 
1277 
1278 
1279 c----------------------------------------------------------------------
1280 
1281       double precision function varifit(x,var)
1282 
1283 c----------------------------------------------------------------------
1284       double precision x,var
1285 
1286       varifit=dexp(-min(50.d0,x**2.d0/var))
1287 
1288       return
1289       end
1290 
1291 
1292 
1293 c----------------------------------------------------------------------
1294 
1295       double precision function Dsoftshval(sy,x,y,b,iscr)
1296 
1297 c----------------------------------------------------------------------
1298 c iscr=-1 sum of om5p (i), i from 0 to 4 - fit of hard
1299 c iscr=0 sum of om5p (i), i from 0 to 4
1300 c iscr=1 sum of om5p (i), i from 0 to 4 * F * F
1301 c iscr=2 sum of om5p (i), i from 1 to 4 (semihard + valence quark)
1302 c----------------------------------------------------------------------
1303       double precision x,om51p,y,xp,xm
1304       include 'epos.inc'
1305       include 'epos.incsem'
1306       include 'epos.incpar'
1307 
1308       Dsoftshval=0.d0
1309 
1310       if(iscr.le.0)then
1311         do i=0,4
1312           Dsoftshval=Dsoftshval+om51p(sy,x,y,b,i)
1313         enddo
1314 
1315 
1316       elseif(iscr.eq.1)then
1317         xp=dsqrt(x)*dexp(y)
1318         if(dabs(xp).ge.1.d-15)then
1319           xm=x/xp
1320         else
1321           xm=1.d0
1322           write(ifch,*)'Warning in Dsoftshval in epos-par'
1323         endif
1324         do i=0,4
1325           Dsoftshval=Dsoftshval+om51p(sy,x,y,b,i)
1326         enddo
1327         Dsoftshval=Dsoftshval*(1.d0-xm)**dble(alplea(icltar))
1328      &                       *(1.d0-xp)**dble(alplea(iclpro))
1329       elseif(iscr.eq.2)then
1330         do i=1,4
1331           Dsoftshval=Dsoftshval+om51p(sy,x,y,b,i)
1332         enddo
1333       endif
1334 
1335 
1336       Dsoftshval=2.d0*Dsoftshval
1337      &     /(x**dble(-alppar)*dble(chad(iclpro)*chad(icltar)))
1338 
1339       if(iscr.eq.-1.and.parDf(3,3).lt.parDf(1,3))Dsoftshval=Dsoftshval
1340      &     -dble(parDf(3,3)*sy**parDf(4,3))
1341 
1342       return
1343       end
1344 
1345 c----------------------------------------------------------------------
1346 
1347       double precision function Dsoftshpar(x)
1348 
1349 c----------------------------------------------------------------------
1350       double precision x,xp,xm
1351       include 'epos.inc'
1352       include 'epos.incpar'
1353 
1354       Dsoftshpar=dble(
1355      &        parDf(1,3)*(sngl(x)*smaxDf)**parDf(2,3)
1356      &       +parDf(3,3)*(sngl(x)*smaxDf)**parDf(4,3) )
1357       xp=dsqrt(x)
1358       xm=xp
1359       Dsoftshpar=Dsoftshpar*(1.d0-xm)**dble(alplea(icltar))
1360      &                       *(1.d0-xp)**dble(alplea(iclpro))
1361       Dsoftshpar=min(max(1.d-15,Dsoftshpar),1.d15)
1362 
1363       return
1364       end