LHAPDF/lhapdf5.3.1/wrapcteq6.f

   1       subroutine CTEQ6evolve(x,Q,pdf)
   2       implicit real*8(a-h,o-z)
   3       include 'parmsetup.inc'
   4       character*16 name(nmxset)
   5       integer nmem(nmxset),ndef(nmxset),mmem
   6       common/NAME/name,nmem,ndef,mmem
   7       real*8 pdf(-6:6)
   8       integer nset
   9       Character Line*80
  10       PARAMETER (MXX = 96, MXQ = 20, MXF = 5, nhess = 40)
  11       PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX)
  12       Common
  13      > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX)
  14      > / CtqPar2 / Nx, Nt, NfMx
  15      > / XQrange / Qini, Qmax, Xmin
  16      > / QCDtable /  Alambda, Nfl, Iorder
  17      > / Masstbl / Amass(6)
  18       common/masses_LHA/cMass,bMass,tMass
  19       data pi / 3.141592653589793d0 /
  20       save
  21 c
  22       call getnset(iset)
  23       call getnmem(iset,imem)
  24 c
  25       U =         X * CtLhCtq6Pdf(imem,1,X,Q)
  26       D =         X * CtLhCtq6Pdf(imem,2,X,Q)
  27       USEA =      X * CtLhCtq6Pdf(imem,-1,X,Q)
  28       DSEA =      X * CtLhCtq6Pdf(imem,-2,X,Q)
  29       STR =       X * CtLhCtq6Pdf(imem,3,X,Q)
  30       CHM =       X * CtLhCtq6Pdf(imem,4,X,Q)
  31       BOT =       X * CtLhCtq6Pdf(imem,5,X,Q)
  32       GLU  =      X * CtLhCtq6Pdf(imem,0,X,Q)
  33 c
  34       pdf(0)  = glu
  35       pdf(1)  = d
  36       pdf(-1) = dsea
  37       pdf(2)  = u
  38       pdf(-2) = usea
  39       pdf(3)  = str
  40       pdf(-3) = str
  41       pdf(4)  = chm
  42       pdf(-4) = chm
  43       pdf(5)  = bot
  44       pdf(-5) = bot
  45       pdf(6)  = 0.0d0
  46       pdf(-6) = 0.0d0
  47
  48       return
  49 *
  50       entry CTEQ6read(nset)
  51
  52       call CtLhbldat1
  53       call CtLhbldat2
  54
  55       call LHct6set
  56
  57       read(1,*)nmem(nset),ndef(nset)    !*** nmem+1=number of members; ndef is not used for anything ***
  58       if(nmem(nset) .gt. nhess) then
  59          print *,'fatal error:  nmem=',nmem(nset),' > nhess=',nhess
  60          stop
  61       endif
  62       Read  (1, '(A)') Line
  63       Read  (1, '(A)') Line
  64       Read  (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
  65       Iorder = Nint(Dr)
  66       Nfl = Nint(Fl)
  67       Alambda = Al
  68
  69       cMass = Amass(4)
  70       bMass = Amass(5)
  71       tMass = Amass(6)
  72
  73       Read  (1, '(A)') Line
  74       Read  (1, *) NX,  NT, NfMx
  75
  76       Read  (1, '(A)') Line
  77       Read  (1, *) QINI, QMAX, (TV(I), I =0, NT)
  78
  79       Read  (1, '(A)') Line
  80       Read  (1, *) XMIN, (XV(I), I =0, NX)
  81
  82       Do 11 Iq = 0, NT
  83          TV(Iq) = Log(Log (TV(Iq) /Al))
  84    11 Continue
  85 C
  86 C                  Since quark = anti-quark for nfl>2 at this stage,
  87 C                  we Read  out only the non-redundent data points
  88 C     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
  89
  90       Nblk = (NX+1) * (NT+1)
  91       Npts =  Nblk  * (NfMx+3)
  92
  93       do ihess = 0,nmem(nset)           !*** new version: allows nmem < nhess ***
  94
  95         Read  (1, '(A)') Line
  96         Read  (1, '(A)') Line
  97         Read  (1, *, IOSTAT=IRET) (UPD(ihess,I), I=1,Npts)
  98
  99       enddo
 100       return
 101 *
 102       entry CTEQ6alfa(alfas,Qalfa)
 103       alfas = pi*CtLhALPI(Qalfa)
 104       return
 105 *
 106       entry CTEQ6init(Eorder,Q2fit)
 107
 108       return
 109 *
 110       entry CTEQ6pdf(mem)
 111 c        imem = mem
 112         call getnset(iset)
 113         call setnmem(iset,mem)
 114       return
 115 *
 116       end
 117
 118       subroutine LHct6set
 119       Implicit Double Precision (A-H,O-Z)
 120       common /ctq6co/ xlast, qlast, nxsave
 121       nxsave = -1000
 122       xlast = -2.
 123       qlast = -2.
 124       return
 125       end
 126
 127 c===========================================================================
 128       Function CtLhPartonX6 (iset,IPRTN, XX, QQ)
 129 c  Given the parton distribution function in the array U in
 130 c  COMMON / PEVLDT / , this routine interpolates to find
 131 c  the parton distribution at an arbitray point in x and q.
 132 c
 133       Implicit Double Precision (A-H,O-Z)
 134
 135       Parameter (MXX = 96, MXQ = 20, MXF = 5, nhess = 40)
 136       Parameter (MXQX= MXQ * MXX,   MXPQX = MXQX * (MXF+3))
 137
 138       Common
 139      > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX)
 140      > / CtqPar2 / Nx, Nt, NfMx
 141      > / XQrange / Qini, Qmax, Xmin
 142
 143       common /ctq6co/ xlast, qlast, nxsave
 144         parameter(nqvec = 4)
 145
 146       Dimension fvec(4), fij(4)
 147       Dimension xvpow(0:mxx)
 148       Data OneP / 1.00001 /
 149       Data xpow / 0.3d0 /       !**** choice of interpolation variable
 150       Save xvpow
 151
 152         save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
 153         save const1 , const2, const3, const4, const5, const6
 154         save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
 155
 156 c store the powers used for interpolation on first call...
 157       if(nx .ne. nxsave) then
 158          nxsave = nx
 159          xvpow(0) = 0.D0
 160          do i = 1, nx
 161             xvpow(i) = xv(i)**xpow
 162          enddo
 163       endif
 164
 165       X = XX
 166       Q = QQ
 167
 168       if((x.lt.xmin).or.(x.gt.1.d0)) print 98,x
 169   98  format(' WARNING:  X=',e12.5,' OUT OF RANGE')
 170       if((q.lt.qini).or.(q.gt.qmax)) print 99,q
 171   99  format(' WARNING:  Q=',e12.5,' OUT OF RANGE')
 172
 173 c skip the initialization in x if same as in the previous call.
 174         if(x .eq. xlast) goto 100
 175         xlast = x
 176
 177 c      -------------    find lower end of interval containing x, i.e.,
 178 c                       get jx such that xv(jx) .le. x .le. xv(jx+1)...
 179       JLx = -1
 180       JU = Nx+1
 181  11   If (JU-JLx .GT. 1) Then
 182          JM = (JU+JLx) / 2
 183          If (X .Ge. XV(JM)) Then
 184             JLx = JM
 185          Else
 186             JU = JM
 187          Endif
 188          Goto 11
 189       Endif
 190 C                     Ix    0   1   2      Jx  JLx         Nx-2     Nx
 191 C                           |---|---|---|...|---|-x-|---|...|---|---|
 192 C                     x     0  Xmin               x                 1
 193 C
 194       If     (JLx .LE. -1) Then
 195         Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX6 x=',x
 196         Stop
 197       ElseIf (JLx .Eq. 0) Then
 198          Jx = 0
 199       Elseif (JLx .LE. Nx-2) Then
 200
 201 C                For interior points, keep x in the middle, as shown above
 202          Jx = JLx - 1
 203       Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
 204
 205 C                  We tolerate a slight over-shoot of one (OneP=1.00001),
 206 C              perhaps due to roundoff or whatever, but not more than that.
 207 C                                      Keep at least 4 points >= Jx
 208          Jx = JLx - 2
 209       Else
 210         Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX6 x=',x
 211         Stop
 212       Endif
 213 C          ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
 214
 215 C                       This is the variable to be interpolated in
 216       ss = x**xpow
 217
 218       If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
 219
 220 c     initiation work for "interior bins": store the lattice points in s...
 221       svec1 = xvpow(jx)
 222       svec2 = xvpow(jx+1)
 223       svec3 = xvpow(jx+2)
 224       svec4 = xvpow(jx+3)
 225
 226       s12 = svec1 - svec2
 227       s13 = svec1 - svec3
 228       s23 = svec2 - svec3
 229       s24 = svec2 - svec4
 230       s34 = svec3 - svec4
 231
 232       sy2 = ss - svec2
 233       sy3 = ss - svec3
 234
 235 c constants needed for interpolating in s at fixed t lattice points...
 236       const1 = s13/s23
 237       const2 = s12/s23
 238       const3 = s34/s23
 239       const4 = s24/s23
 240       s1213 = s12 + s13
 241       s2434 = s24 + s34
 242       sdet = s12*s34 - s1213*s2434
 243       tmp = sy2*sy3/sdet
 244       const5 = (s34*sy2-s2434*sy3)*tmp/s12
 245       const6 = (s1213*sy2-s12*sy3)*tmp/s34
 246
 247       EndIf
 248
 249 100     continue
 250
 251 c skip the initialization in q if same as in the previous call.
 252         if(q .eq. qlast) goto 110
 253         qlast = q
 254
 255       tt = log(log(Q/Al))
 256
 257 c         --------------Now find lower end of interval containing Q, i.e.,
 258 c                          get jq such that qv(jq) .le. q .le. qv(jq+1)...
 259       JLq = -1
 260       JU = NT+1
 261  12   If (JU-JLq .GT. 1) Then
 262          JM = (JU+JLq) / 2
 263          If (tt .GE. TV(JM)) Then
 264             JLq = JM
 265          Else
 266             JU = JM
 267          Endif
 268          Goto 12
 269        Endif
 270
 271       If     (JLq .LE. 0) Then
 272          Jq = 0
 273       Elseif (JLq .LE. Nt-2) Then
 274 C                                  keep q in the middle, as shown above
 275          Jq = JLq - 1
 276       Else
 277 C                         JLq .GE. Nt-1 case:  Keep at least 4 points >= Jq.
 278         Jq = Nt - 3
 279
 280       Endif
 281 C                                   This is the interpolation variable in Q
 282
 283       If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
 284 c                                        store the lattice points in t...
 285       tvec1 = Tv(jq)
 286       tvec2 = Tv(jq+1)
 287       tvec3 = Tv(jq+2)
 288       tvec4 = Tv(jq+3)
 289
 290       t12 = tvec1 - tvec2
 291       t13 = tvec1 - tvec3
 292       t23 = tvec2 - tvec3
 293       t24 = tvec2 - tvec4
 294       t34 = tvec3 - tvec4
 295
 296       ty2 = tt - tvec2
 297       ty3 = tt - tvec3
 298
 299       tmp1 = t12 + t13
 300       tmp2 = t24 + t34
 301
 302       tdet = t12*t34 - tmp1*tmp2
 303
 304       EndIf
 305
 306 110     continue
 307
 308 c get the pdf function values at the lattice points...
 309 c In this code, we store 8 flavors: u,ubar,d,dbar,s=sbar,c=cbar,b=bbar,g
 310 c hence Iprtn=3,4,5 (s,c,b) are obtained from -3,-4,-5 (sbar,cbar,bbar)
 311
 312       If (Iprtn .GE. 3) Then
 313          Ip = - Iprtn
 314       Else
 315          Ip = Iprtn
 316       EndIf
 317       jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
 318
 319       Do it = 1, nqvec
 320
 321          J1  = jtmp + it*(NX+1)
 322
 323        If (Jx .Eq. 0) Then
 324 C                          For the first 4 x points, interpolate x^2*f(x,Q)
 325 C                           This applies to the two lowest bins JLx = 0, 1
 326 C            We cannot put the JLx.eq.1 bin into the "interior" section
 327 C                           (as we do for q), since Upd(J1) is undefined.
 328          fij(1) = 0
 329          fij(2) = Upd(iset,J1+1) * XV(1)**2
 330          fij(3) = Upd(iset,J1+2) * XV(2)**2
 331          fij(4) = Upd(iset,J1+3) * XV(3)**2
 332 C
 333 C                 Use CtLhPolint which allows x to be anywhere w.r.t. the grid
 334
 335          Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx)
 336
 337          If (x .GT. 0D0)  Fvec(it) =  Fx / x**2
 338 C                                              Pdf is undefined for x.eq.0
 339        ElseIf  (JLx .Eq. Nx-1) Then
 340 C                                                This is the highest x bin:
 341
 342 c** fix allow 4 consecutive elements with iset... mrw 19.9.2005
 343         fij(1) = Upd(iset,j1)
 344         fij(2) = Upd(iset,j1+1)
 345         fij(3) = Upd(iset,j1+2)
 346         fij(4) = Upd(iset,j1+3)
 347         Call CtLhPolint4 (XVpow(Nx-3), Fij(1), 4, ss, Fx, Dfx)
 348
 349         Fvec(it) = Fx
 350
 351        Else
 352
 353 C for all interior points, use Jon's in-line function
 354 C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
 355 c (This is cubic spline interpolation, as used by cteq; it was
 356 c changed to polint in previous Durham releases (jcp).)
 357          sf2 = Upd(iset,J1+1)
 358          sf3 = Upd(iset,J1+2)
 359
 360          Fvec(it) = (const5*(Upd(iset,J1)
 361      &                       - sf2*const1 + sf3*const2)
 362      &             + const6*(Upd(iset,J1+3)
 363      &                       + sf2*const3 - sf3*const4)
 364      &             + sf2*sy3 - sf3*sy2) / s23
 365
 366        Endif
 367
 368       enddo
 369 C                                   We now have the four values Fvec(1:4)
 370 c     interpolate in t...
 371
 372       If (JLq .LE. 0) Then
 373 C                         1st Q-bin, as well as extrapolation to lower Q
 374         Call CtLhPolint4(TV(0), Fvec(1), 4, tt, ff, Dfq)
 375
 376       ElseIf (JLq .GE. Nt-1) Then
 377 C                         Last Q-bin, as well as extrapolation to higher Q
 378         Call CtLhPolint4(TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
 379       Else
 380 C                         Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
 381 C       which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
 382 C                         the full range QV(0:Nt)  (in contrast to XV)
 383         tf2 = fvec(2)
 384         tf3 = fvec(3)
 385
 386         g1 = ( tf2*t13 - tf3*t12) / t23
 387         g4 = (-tf2*t34 + tf3*t24) / t23
 388
 389         h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
 390      &    +  (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
 391
 392         ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
 393       EndIf
 394
 395       CtLhPartonX6 = ff
 396
 397       Return
 398 C                                       ********************
 399       End
 400 c===========================================================================
 401       Function CtLhCtq6Pdf (iset,Iparton, X, Q)
 402       Implicit Double Precision (A-H,O-Z)
 403       Logical Warn
 404       Common
 405      > / CtqPar2 / Nx, Nt, NfMx
 406      > / QCDtable /  Alambda, Nfl, Iorder
 407
 408       Data Warn /.true./
 409       save Warn
 410
 411       If (X .lt. 0D0 .or. X .gt. 1D0) Then
 412         Print *, 'X out of range in CtLhCtq6Pdf: ', X
 413         Stop
 414       Endif
 415       If (Q .lt. Alambda) Then
 416         Print *, 'Q out of range in CtLhCtq6Pdf: ', Q
 417         Stop
 418       Endif
 419
 420 c added to force pdf = 0.0 at x=1.0 exactly - mrw
 421       if(x .eq. 1.0d0) then
 422           CtLhCtq6Pdf = 0.0d0
 423           return
 424       endif
 425 c
 426       If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
 427          If (Warn) Then
 428 C        put a warning for calling extra flavor.
 429              Warn = .false.
 430              Print *, 'Warning: Iparton out of range in CtLhCtq6Pdf: '
 431      >              , Iparton
 432          Endif
 433          CtLhCtq6Pdf = 0D0
 434          Return
 435       Endif
 436
 437       CtLhCtq6Pdf = CtLhPartonX6 (iset,Iparton, X, Q)
 438       if(CtLhCtq6Pdf.lt.0.D0)  CtLhCtq6Pdf = 0.D0
 439
 440       Return
 441
 442 C                             ********************
 443       End
 444       SUBROUTINE CtLhPOLINT4 (XA,YA,N,X,Y,DY)
 445 c fast version of polint, valid only for N=4
 446 c Have explicitly unrolled the loops.
 447
 448       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 449
 450       PARAMETER (NMAX=4)
 451       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 452
 453         if(n .ne. 4) then
 454            print *,'fatal CtLhPolint4 call',n
 455            stop
 456         endif
 457
 458       NS=1
 459       DIF=ABS(X-XA(1))
 460
 461         DIFT=ABS(X-XA(1))
 462         IF (DIFT.LT.DIF) THEN
 463           NS=1
 464           DIF=DIFT
 465         ENDIF
 466         C(1)=YA(1)
 467         D(1)=YA(1)
 468
 469         DIFT=ABS(X-XA(2))
 470         IF (DIFT.LT.DIF) THEN
 471           NS=2
 472           DIF=DIFT
 473         ENDIF
 474         C(2)=YA(2)
 475         D(2)=YA(2)
 476
 477         DIFT=ABS(X-XA(3))
 478         IF (DIFT.LT.DIF) THEN
 479           NS=3
 480           DIF=DIFT
 481         ENDIF
 482         C(3)=YA(3)
 483         D(3)=YA(3)
 484
 485         DIFT=ABS(X-XA(4))
 486         IF (DIFT.LT.DIF) THEN
 487           NS=4
 488           DIF=DIFT
 489         ENDIF
 490         C(4)=YA(4)
 491         D(4)=YA(4)
 492
 493
 494       Y=YA(NS)
 495       NS=NS-1
 496
 497
 498           HO=XA(1)-X
 499           HP=XA(2)-X
 500           W=C(2)-D(1)
 501           DEN=W/(HO-HP)
 502           D(1)=HP*DEN
 503           C(1)=HO*DEN
 504
 505
 506           HO=XA(2)-X
 507           HP=XA(3)-X
 508           W=C(3)-D(2)
 509           DEN=W/(HO-HP)
 510           D(2)=HP*DEN
 511           C(2)=HO*DEN
 512
 513
 514           HO=XA(3)-X
 515           HP=XA(4)-X
 516           W=C(4)-D(3)
 517           DEN=W/(HO-HP)
 518           D(3)=HP*DEN
 519           C(3)=HO*DEN
 520
 521         IF (2*NS.LT.3)THEN
 522           DY=C(NS+1)
 523         ELSE
 524           DY=D(NS)
 525           NS=NS-1
 526         ENDIF
 527         Y=Y+DY
 528
 529           HO=XA(1)-X
 530           HP=XA(3)-X
 531           W=C(2)-D(1)
 532           DEN=W/(HO-HP)
 533           D(1)=HP*DEN
 534           C(1)=HO*DEN
 535
 536           HO=XA(2)-X
 537           HP=XA(4)-X
 538           W=C(3)-D(2)
 539           DEN=W/(HO-HP)
 540           D(2)=HP*DEN
 541           C(2)=HO*DEN
 542
 543         IF (2*NS.LT.2)THEN
 544           DY=C(NS+1)
 545         ELSE
 546           DY=D(NS)
 547           NS=NS-1
 548         ENDIF
 549         Y=Y+DY
 550
 551           HO=XA(1)-X
 552           HP=XA(4)-X
 553           W=C(2)-D(1)
 554           DEN=W/(HO-HP)
 555           D(1)=HP*DEN
 556           C(1)=HO*DEN
 557
 558         IF (2*NS.LT.4-3)THEN
 559           DY=C(NS+1)
 560         ELSE
 561           DY=D(NS)
 562           NS=NS-1
 563         ENDIF
 564         Y=Y+DY
 565
 566       RETURN
 567       END
 568       SUBROUTINE CTLHPOLINT3 (XA,YA,N,X,Y,DY)
 569 c fast version of polint, valid only for N=3
 570 c Have explicitly unrolled the loops.
 571       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 572       PARAMETER (NMAX=3)
 573       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 574         if(n .ne. 3) then
 575            print *,'fatal CtLhPolint3 call',n
 576            stop
 577         endif
 578       NS=1
 579       DIF=ABS(X-XA(1))
 580         DIFT=ABS(X-XA(1))
 581         IF (DIFT.LT.DIF) THEN
 582           NS=1
 583           DIF=DIFT
 584         ENDIF
 585         C(1)=YA(1)
 586         D(1)=YA(1)
 587         DIFT=ABS(X-XA(2))
 588         IF (DIFT.LT.DIF) THEN
 589           NS=2
 590           DIF=DIFT
 591         ENDIF
 592         C(2)=YA(2)
 593         D(2)=YA(2)
 594         DIFT=ABS(X-XA(3))
 595         IF (DIFT.LT.DIF) THEN
 596           NS=3
 597           DIF=DIFT
 598         ENDIF
 599         C(3)=YA(3)
 600         D(3)=YA(3)
 601       Y=YA(NS)
 602       NS=NS-1
 603           HO=XA(1)-X
 604           HP=XA(2)-X
 605           W=C(2)-D(1)
 606           DEN=W/(HO-HP)
 607           D(1)=HP*DEN
 608           C(1)=HO*DEN
 609           HO=XA(2)-X
 610           HP=XA(3)-X
 611           W=C(3)-D(2)
 612           DEN=W/(HO-HP)
 613           D(2)=HP*DEN
 614           C(2)=HO*DEN
 615         IF (2*NS.LT.2)THEN
 616           DY=C(NS+1)
 617         ELSE
 618           DY=D(NS)
 619           NS=NS-1
 620         ENDIF
 621         Y=Y+DY
 622           HO=XA(1)-X
 623           HP=XA(3)-X
 624           W=C(2)-D(1)
 625           DEN=W/(HO-HP)
 626           D(1)=HP*DEN
 627           C(1)=HO*DEN
 628         IF (2*NS.LT.1)THEN
 629           DY=C(NS+1)
 630         ELSE
 631           DY=D(NS)
 632           NS=NS-1
 633         ENDIF
 634         Y=Y+DY
 635       RETURN
 636       END