LHAPDF/lhapdf5.3.1/wrapcteq65.f

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