File indexing completed on 2023-10-25 09:48:50
0001
0002
0003
0004
0005 subroutine paramini(imod)
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
0025
0026 call Class('paramini ')
0027
0028
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
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
0057
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
0078
0079
0080 call pompar(alpsf,betsf,0) !soft (taking into account hard)
0081 call pompar(alpsh,betsh,1) !sh
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
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
0104
0105
0106
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
0116
0117 parDf(1,3)=alpsf
0118 parDf(2,3)=betsf
0119 parDf(3,3)=alpsh
0120 parDf(4,3)=betsh
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
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
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
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
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
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
0226 betDp(2,iclpro,icltar)=-alpdifs+alppar
0227 betDpp(2,iclpro,icltar)=-alpdifs+alppar
0228
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
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 chisquare 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
