LHAPDF/lhapdf5.5.1/src/lhaglue.f

   1 ! -*- F90 -*-
   2
   3
   4 !   LHAGLUE Interface to LHAPDF library of modern parton
   5 !    density functions (PDF) with uncertainties
   6 !
   7 ! Authors for v4: Dimitri Bourilkov, Craig Group, Mike Whalley
   8 !
   9 ! Authors for v3: Dimitri Bourilkov, Craig Group, Mike Whalley
  10 !
  11 ! Author for v1 and v2: Dimitri Bourilkov  bourilkov@mailaps.org
  12 !                       University of Florida
  13 !
  14 ! HERWIG interface by Dimitri Bourilkov and Craig Group
  15 !
  16 ! New numbering scheme and upgrade for LHAPDF v2.1
  17 ! by Dimitri Bourilkov and Mike Whalley
  18 !
  19 ! For more information, or when you cite this interface, currently
  20 ! the official reference is:
  21 ! D.Bourilkov, "Study of Parton Density Function Uncertainties with
  22 ! LHAPDF and PYTHIA at LHC", hep-ph/0305126.
  23 !
  24 ! The official LHAPDF page is:
  25 !
  26 !    http://durpdg.dur.ac.uk/lhapdf/index.html
  27 !
  28 ! The interface contains four subroutines (similar to PDFLIB).
  29 ! It can be used seamlessly by Monte Carlo generators
  30 ! interfaced to PDFLIB or in stand-alone mode.
  31 !
  32 !     For initialization (called once)
  33 !
  34 !     PDFSET(PARM,VALUE)
  35 !
  36 !     For the proton/pion structure functions
  37 !
  38 !     STRUCTM(X,Q,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
  39 !
  40 !     For the photon structure functions
  41 !
  42 !     STRUCTP(X,Q2,P2,IP2,UPV,DNV,USEA,DSEA,STR,CHM,BOT,TOP,GLU)
  43 !
  44 !     For statistics ON structure functions (under/over-flows)
  45 !
  46 !     PDFSTA
  47 !
  48 ! This interface can be invoked in 3 ways depending
  49 ! on the value of PARM(1) provided by the user when
  50 ! calling PDFSET(PARM,VALUE):
  51 !
  52 !     For PYTHIA:         PARM(1).EQ.'NPTYPE'
  53 !       (this is set automatically by PYTHIA)
  54 !
  55 !     For HERWIG:         PARM(1).EQ.'HWLHAPDF'
  56 !       (set by the USER e.g. in the main program like this:
  57 !           AUTPDF(1) = 'HWLHAPDF'
  58 !           AUTPDF(2) = 'HWLHAPDF'                         )
  59 !
  60 !     For Stand-alone:    PARM(1).EQ.'DEFAULT'
  61 !       (can be used for PDF studies or when interfacing
  62 !        new generators)
  63 !
  64 ! The LHAPDF set/member is selected depending on the value of:
  65 !
  66 !         PYTHIA:   ABS(MSTP(51)) - proton
  67 !                   ABS(MSTP(53)) - pion
  68 !                   ABS(MSTP(55)) - photon
  69 !
  70 !         HERWIG:   ABS(INT(VALUE(1)))
  71 !
  72 !    STAND-ALONE:   ABS(INT(VALUE(1)))
  73 !
  74 !
  75 !         CONTROL switches:
  76 !        ==================
  77 !
  78 !      THE LOCATION OF THE LHAPDF LIBRARY HAS TO BE SPECIFIED
  79 !      AS DESCRIBED BELOW (the rest is optional)
  80 !
  81 !      if the user does nothing, sensible defaults
  82 !      are active; to change the behaviour, the corresponding
  83 !      values of LHAPARM() should be set to the values given below
  84 !
  85 !   Location of the LHAPDF library of PDFs (pathname):
  86 !      uses common block /LHAPDFC/
  87 !
  88 !    If the user does nothing => default = subdir PDFsets of the
  89 !                               current directory (can be real subdir
  90 !                               OR a soft link to the real location)
  91 !    If the user sets LHAPATH => supplied by the USER who defines the
  92 !                         path in common block COMMON/LHAPDFC/LHAPATH
  93 !                         BEFORE calling PDFSET
  94 !
  95 !    Other controls:
  96 !    ===============
  97 !      use common block /LHACONTROL/
  98 !
  99 !   Collect statistics on under/over-flow requests for PDFs
 100 !   outside their validity ranges in X and Q**2
 101 !   (call PDFSTA at end of run to print it out)
 102 !
 103 !       LHAPARM(16).EQ.'NOSTAT' => No statistics (faster)
 104 !       LHAPARM(16).NE.'NOSTAT' => Default: collect statistics
 105 !
 106 !   Option to use the values for the strong coupling alpha_s
 107 !   as computed in LHAPDF in the MC generator
 108 !   (to ensure uniformity between the MC generator and the PDF set)
 109 !   WARNING: implemented ONLY for PYTHIA in LHAPDFv4
 110 !
 111 !       LHAPARM(17).EQ.'LHAPDF' => Use alpha_s from LHAPDF
 112 !       LHAPARM(17).NE.'LHAPDF' => Default (same as LHAPDF v1/v3)
 113 !
 114 !   Extrapolation of PDFs outside LHAPDF validity range given by
 115 !   [Xmin,Xmax] and [Q2min,Q2max]; DEFAULT => PDFs "frozen" at the
 116 !   boundaries
 117 !
 118 !       LHAPARM(18).EQ.'EXTRAPOLATE' => Extrapolate PDFs on OWN RISK
 119 !                           WARNING: Crazy values can be returned
 120 !
 121 !   Printout of initialization information in PDFSET (by default)
 122 !
 123 !       LHAPARM(19).EQ.'SILENT' => No printout (silent mode)
 124 !       LHAPARM(19).EQ.'LOWKEY' => Print 5 times (almost silent mode)
 125 !
 126 !
 127 !*********************************************************************
 128 !
 129 ! $Id: lhaglue.f 356 2008-08-28 15:58:02Z buckley $
 130 !
 131 ! $Log$
 132 ! Revision 1.7  2005/12/02 14:50:54  whalley
 133 ! Changes for new CTEQ code/AB sets
 134 !
 135 ! Revision 1.6  2005/10/18 15:35:48  whalley
 136 ! fix to allow LHAPATH to be user defined as well as lhapdf-config
 137 !
 138 ! Revision 1.5  2005/10/18 11:47:48  whalley
 139 ! Change to only set LHAPATH once per run
 140 !
 141 ! Revision 1.1.1.2  1996/10/30 08:29:06  cernlib
 142 ! Version 7.04
 143 !
 144 ! Revision 1.1.1.1  1996/04/12 15:29:26  plothow
 145 ! Version 7.01
 146 !
 147 ! v5.0  06-Oct-2005  Major change to allow multiset-initializations
 148 ! v4.0  28-Apr-2005  PDFSTA routine; option to use Alfa_s from LHAPDF
 149 ! v4.0  21-Mar-2005  Photon/pion/new p PDFs, updated for LHAPDF v4
 150 ! v3.1  26-Apr-2004  New numbering scheme, updated for LHAPDF v2/v3
 151 ! v3.0  23-Jan-2004  HERWIG interface added
 152 ! v2.0  20-Sep-2003  PDFLIB style adopted
 153 ! v1.0  05-Mar-2003  First working version from PYTHIA to LHAPDF v1
 154 !
 155 ! interface to LHAPDF library
 156 !*********************************************************************
 157
 158
 159 ! PDFSET
 160 ! Initialization for use of parton distributions
 161 !  according to the LHAPDF interface.
 162 !
 163 ! v4.0  28-Apr-2005  Option to use Alfa_s from LHAPDF
 164 ! v4.0  21-Mar-2005  Photon/pion/new p PDFs, updated for LHAPDF v4
 165 ! v3.1  26-Apr-2004  New numbering scheme
 166 ! v3.0  23-Jan-2004  HERWIG interface added
 167 !
 168 ! Interface to LHAPDF library
 169  subroutine pdfset(parm,value, mstu11, mstp51, mstp53, mstp55, qcdl4, qcdl5, axmin, axmax, aq2min, aq2max)
 170  entry pdfset_herwig(parm, value)
 171   ! Double precision and integer declarations.
 172   implicit double precision(a-h, o-z)
 173   implicit integer(i-n)
 174   ! Common blocks
 175   include 'commonlhapdf.inc'
 176   include 'commonlhasets.inc'
 177   include 'commonlhapdfc.inc'
 178   include 'commonlhacontrol.inc'
 179   include 'commonlhaglsta.inc'
 180   ! additions for multiset use
 181   double precision xxmin(nmxset),xxmax(nmxset),qq2min(nmxset),qq2max(nmxset)
 182   save xxmin,xxmax,qq2min,qq2max
 183   ! Interface to LHAPDFLIB.
 184   double precision qcdlha4, qcdlha5
 185   integer nfllha
 186   common/lhapdfr/qcdlha4, qcdlha5, nfllha
 187   save /lhapdfr/
 188   integer lhaextrp
 189   common/lhapdfe/lhaextrp
 190   save /lhapdfe/
 191   integer lhasilent
 192   common/lhasilent/lhasilent
 193   save /lhasilent/
 194   ! Interface to PDFLIB.
 195   common/w50511/ nptypepdfl,ngrouppdfl,nsetpdfl,modepdfl,nflpdfl,lopdfl,tmaspdfl
 196   save /w50511/
 197   double precision tmaspdfl
 198   double precision qcdl4,qcdl5
 199   double precision axmin, axmax, aq2min, aq2max
 200   ! Interface to PDFLIB.
 201   common/w50513/xmin,xmax,q2min,q2max
 202   save /w50513/
 203   double precision xmin,xmax,q2min,q2max
 204   ! Local arrays and character variables (NOT USED here DB)
 205   character*20 parm(20)
 206   double precision value(20)
 207   integer lhapathlen
 208   integer :: lhainput = 1
 209   !integer lhaselect
 210   integer lhaprint
 211   integer lhaonce
 212   integer lhafive
 213   save lhaonce
 214   save lhafive
 215   data lhaonce/0/
 216   data lhafive/0/
 217   logical first
 218   data first/.true./
 219   character*512 dirpath
 220   save first
 221   character*1000 chroot
 222   chroot=' '
 223   ! Initialise common blocks
 224   call commoninit()
 225
 226 !  if (first) then
 227 !     call getdirpath(dirpath)
 228 !     first = .FALSE.
 229 !  endif
 230   if(first .AND. (LHAPARM(20).NE.'LHAPATH')) then
 231 !...overide the default PDFsets path
 232 ! ... check first if the environmental variable LHAPATH is set
 233          call getenv('LHAPATH',lhapath)
 234          if(lhapath.eq.'') then
 235 !     The environment variable LHAPATH is not set.
 236 !     Take the data from $ALICE_ROOT/LHAPDF/PDFsets
 237             CALL GETENV('ALICE_ROOT',CHROOT)
 238             LNROOT = LNBLNK(CHROOT)
 239             IF(LNROOT.LE.0) THEN
 240                LHAPATH='PDFsets' ! Default value
 241             ELSE
 242                LHAPATH=CHROOT(1:LNROOT)//'/LHAPDF/PDFsets'
 243             ENDIF
 244          endif
 245       first=.FALSE.
 246       endif
 247
 248
 249   ! Init
 250   lhaextrp = 0
 251   if(lhaparm(18).EQ.'EXTRAPOLATE') then  ! Extrapolate PDFs on own risk
 252      lhaextrp = 1
 253   endif
 254   lhasilent = 0
 255   if (lhaparm(19).EQ.'SILENT') then    !  No printout (silent MODE)
 256      lhasilent = 1
 257   elseif (lhaparm(19).EQ.'LOWKEY') then ! print 5 times (lowkey mode)
 258      if (lhafive .lt. 6) then
 259         lhafive = lhafive + 1
 260      else
 261         lhasilent = 1
 262      endif
 263   endif
 264   if (parm(1).EQ.'NPTYPE') then        !  pythia
 265      lhaprint = mstu11
 266      if(value(1) .eq. 1) then         !   nucleon
 267         lhainput = abs(mstp51)
 268      elseif(value(1) .eq. 2) then     !   pion
 269         lhainput = abs(mstp53)
 270      elseif(value(1) .eq. 3) then     !   photon
 271         lhainput = abs(mstp55)
 272      endif
 273      if (lhasilent .ne. 1) print *,'==== PYTHIA WILL USE LHAPDF ===='
 274   elseif(parm(1).EQ.'HWLHAPDF') then  !  herwig
 275      lhainput = abs(int(value(1)))
 276      if(lhaonce.eq.lhainput) return
 277      if(lhasilent .ne. 1) print *,'==== HERWIG WILL USE LHAPDF ===='
 278      lhaprint = 6
 279      lhaonce = lhainput
 280   elseif(parm(1).EQ.'DEFAULT') then  !  stand-alone
 281      lhainput = abs(int(value(1)))
 282      if(lhaonce.eq.lhainput) return
 283      if(lhasilent .ne. 1) print *,'==== STAND-ALONE LHAGLUE MODE TO USE LHAPDF ===='
 284      lhaprint = 6
 285      lhaonce = lhainput
 286   else
 287      print *,'== UNKNOWN LHAPDF INTERFACE CALL! STOP EXECUTION! =='
 288      stop
 289   endif
 290   ! Initialize parton distributions: LHAPDFLIB.
 291   lhapathlen=index(lhapath,' ') - 1
 292   lhaset = lhainput
 293   xmin = 1.0D-6      ! X_min for current PDF set
 294   xmax = 1.0D0       ! X_max for current PDF set
 295   q2min = 1.0D0**2   ! Q**2_min scale for current PDF set [GeV]
 296   q2max = 1.0D5**2   ! Q**2_max scale for current PDF set [GeV]
 297   !
 298   ! Protons
 299   !
 300   ! CTEQ Family
 301   if ((lhainput .ge. 10000) .and. (lhainput .le. 19999)) then
 302      q2max = 1.0d08
 303      if ((lhainput .ge. 10000) .and. (lhainput .le. 10040)) then
 304         lhaset = 10000
 305         lhaname=lhapath(1:lhapathlen)//'/cteq6.LHpdf'
 306         q2min = 1.69d0
 307      elseif((lhainput .ge. 10041) .and. (lhainput .le. 10041)) then
 308         lhaset = 10041
 309         lhaname=lhapath(1:lhapathlen)//'/cteq6l.LHpdf'
 310         q2min = 1.69d0
 311      elseif((lhainput .ge. 10042) .and. (lhainput .le. 10042)) then
 312         lhaset = 10042
 313         lhaname=lhapath(1:lhapathlen)//'/cteq6ll.LHpdf'
 314         q2min = 1.69d0
 315      elseif((lhainput .ge. 10050) .and. (lhainput .le. 10090)) then
 316         lhaset = 10050
 317         lhaname=lhapath(1:lhapathlen)//'/cteq6mE.LHgrid'
 318         q2min = 1.69d0
 319      elseif((lhainput .ge. 10100) .and. (lhainput .le. 10140)) then
 320         lhaset = 10100
 321         lhaname=lhapath(1:lhapathlen)//'/cteq61.LHpdf'
 322         q2min = 1.69d0
 323      elseif((lhainput .ge. 10150) .and. (lhainput .le. 10190)) then
 324         lhaset = 10150
 325         lhaname=lhapath(1:lhapathlen)//'/cteq61.LHgrid'
 326         q2min = 1.69d0
 327      elseif((lhainput .ge. 10250) .and. (lhainput .le. 10269)) then
 328         lhaset = 10250
 329         lhaname=lhapath(1:lhapathlen)//'/cteq6AB.LHgrid'
 330         q2min = 1.69d0
 331      elseif((lhainput .ge. 10350) .and. (lhainput .le. 10390)) then
 332         lhaset = 10350
 333         lhaname=lhapath(1:lhapathlen)//'/cteq65.LHgrid'
 334         q2min = 1.69d0
 335         q2max = 1.0d10
 336         xmin =  1.0d-7
 337      elseif((lhainput .ge. 10450) .and. (lhainput .le. 10456)) then
 338         lhaset = 10450
 339         lhaname=lhapath(1:lhapathlen)//'/cteq65c.LHgrid'
 340         q2min = 1.69d0
 341         q2max = 1.0d10
 342         xmin =  1.0d-7
 343      elseif((lhainput .ge. 10460) .and. (lhainput .le. 10467)) then
 344         lhaset = 10460
 345         lhaname=lhapath(1:lhapathlen)//'/cteq65s.LHgrid'
 346         q2min = 1.69d0
 347         q2max = 1.0d10
 348         xmin =  1.0d-7
 349      elseif((lhainput .ge. 10550) .and. (lhainput .le. 10594)) then
 350         lhaset = 10550
 351         lhaname=lhapath(1:lhapathlen)//'/cteq66.LHgrid'
 352         q2min = 1.69d0
 353         q2max = 1.0d10
 354         xmin =  1.0d-8
 355      elseif((lhainput .ge. 10650) .and. (lhainput .le. 10653)) then
 356         lhaset = 10650
 357         lhaname=lhapath(1:lhapathlen)//'/cteq66c.LHgrid'
 358         q2min = 1.69d0
 359         q2max = 1.0d10
 360         xmin =  1.0d-8
 361      elseif((lhainput .ge. 10660) .and. (lhainput .le. 10663)) then
 362         lhaset = 10660
 363         lhaname=lhapath(1:lhapathlen)//'/cteq66a.LHgrid'
 364         q2min = 1.69d0
 365         q2max = 1.0d10
 366         xmin =  1.0d-8
 367      elseif((lhainput .ge. 10670) .and. (lhainput .le. 10677)) then
 368         lhaset = 10670
 369         lhaname=lhapath(1:lhapathlen)//'/cteq6lg.LHgrid'
 370         q2min = 1.69d0
 371      elseif((lhainput .ge. 19050) .and. (lhainput .le. 19050)) then
 372         lhaset = 19050
 373         lhaname=lhapath(1:lhapathlen)//'/cteq5m.LHgrid'
 374         xmin=1.0d-5
 375      elseif((lhainput .ge. 19051) .and. (lhainput .le. 19051)) then
 376         lhaset = 19051
 377         lhaname=lhapath(1:lhapathlen)//'/cteq5m1.LHgrid'
 378         xmin=1.0d-5
 379      elseif((lhainput .ge. 19053) .and. (lhainput .le. 19053)) then
 380         lhaset = 19053
 381         lhaname=lhapath(1:lhapathlen)//'/cteq5f3.LHgrid'
 382         xmin=1.0d-5
 383      elseif((lhainput .ge. 19054) .and. (lhainput .le. 19054)) then
 384         lhaset = 19054
 385         lhaname=lhapath(1:lhapathlen)//'/cteq5f4.LHgrid'
 386         xmin=1.0d-5
 387      elseif((lhainput .ge. 19060) .and. (lhainput .le. 19060)) then
 388         lhaset = 19060
 389         lhaname=lhapath(1:lhapathlen)//'/cteq5d.LHgrid'
 390         xmin=1.0d-5
 391      elseif((lhainput .ge. 19070) .and. (lhainput .le. 19070)) then
 392         lhaset = 19070
 393         lhaname=lhapath(1:lhapathlen)//'/cteq5l.LHgrid'
 394         xmin=1.0d-5
 395      elseif((lhainput .ge. 19150) .and. (lhainput .le. 19150)) then
 396         lhaset = 19150
 397         lhaname=lhapath(1:lhapathlen)//'/cteq4m.LHgrid'
 398         q2min = 2.56d0
 399         xmin=1.0d-5
 400      elseif((lhainput .ge. 19160) .and. (lhainput .le. 19160)) then
 401         lhaset = 19160
 402         lhaname=lhapath(1:lhapathlen)//'/cteq4d.LHgrid'
 403         q2min = 2.56d0
 404         xmin=1.0d-5
 405      elseif((lhainput .ge. 19170) .and. (lhainput .le. 19170)) then
 406         lhaset = 19170
 407         lhaname=lhapath(1:lhapathlen)//'/cteq4l.LHgrid'
 408         q2min = 2.56d0
 409         xmin=1.0d-5
 410      else
 411         write(lhaprint,5150)  lhaset
 412         stop
 413      endif
 414      ! MRST Family
 415   elseif((lhainput .ge. 20000) .and. (lhainput .le. 29999)) then
 416      q2min = 1.25d0
 417      q2max = 1.0d07
 418      xmin = 1.0d-5
 419      if((lhainput .ge. 20000) .and. (lhainput .le. 20004)) then
 420         lhaset = 20000
 421         lhaname=lhapath(1:lhapathlen)//'/MRST2001nlo.LHpdf'
 422      elseif((lhainput .ge. 20050) .and. (lhainput .le. 20054)) then
 423         lhaset = 20050
 424         lhaname=lhapath(1:lhapathlen)//'/MRST2001nlo.LHgrid'
 425      elseif((lhainput .ge. 20060) .and. (lhainput .le. 20061)) then
 426         lhaset = 20060
 427         lhaname=lhapath(1:lhapathlen)//'/MRST2001lo.LHgrid'
 428      elseif((lhainput .ge. 20070) .and. (lhainput .le. 20074)) then
 429         lhaset = 20070
 430         lhaname=lhapath(1:lhapathlen)//'/MRST2001nnlo.LHgrid'
 431      elseif((lhainput .ge. 20100) .and. (lhainput .le. 20130)) then
 432         lhaset = 20100
 433         lhaname=lhapath(1:lhapathlen)//'/MRST2001E.LHpdf'
 434      elseif((lhainput .ge. 20150) .and. (lhainput .le. 20180)) then
 435         lhaset = 20150
 436         lhaname=lhapath(1:lhapathlen)//'/MRST2001E.LHgrid'
 437      elseif((lhainput .ge. 20200) .and. (lhainput .le. 20201)) then
 438         lhaset = 20200
 439         lhaname=lhapath(1:lhapathlen)//'/MRST2002nlo.LHpdf'
 440      elseif((lhainput .ge. 20250) .and. (lhainput .le. 20251)) then
 441         lhaset = 20250
 442         lhaname=lhapath(1:lhapathlen)//'/MRST2002nlo.LHgrid'
 443      elseif((lhainput .ge. 20270) .and. (lhainput .le. 20271)) then
 444         lhaset = 20270
 445         lhaname=lhapath(1:lhapathlen)//'/MRST2002nnlo.LHgrid'
 446      elseif((lhainput .ge. 20300) .and. (lhainput .le. 20301)) then
 447         lhaset = 20300
 448         lhaname=lhapath(1:lhapathlen)//'/MRST2003cnlo.LHpdf'
 449         q2min = 10.0d0
 450         xmin = 1.0d-3
 451      elseif((lhainput .ge. 20350) .and. (lhainput .le. 20351)) then
 452         lhaset = 20350
 453         lhaname=lhapath(1:lhapathlen)//'/MRST2003cnlo.LHgrid'
 454         q2min = 10.0d0
 455         xmin = 1.0d-3
 456      elseif((lhainput .ge. 20370) .and. (lhainput .le. 20371)) then
 457         lhaset = 20370
 458         lhaname=lhapath(1:lhapathlen)//'/MRST2003cnnlo.LHgrid'
 459         q2min = 7.0d0
 460         xmin = 1.0d-3
 461      elseif((lhainput .ge. 20400) .and. (lhainput .le. 20401)) then
 462         lhaset = 20400
 463         lhaname=lhapath(1:lhapathlen)//'/MRST2004nlo.LHpdf'
 464      elseif((lhainput .ge. 20406) .and. (lhainput .le. 20407)) then
 465         lhaset = 20406
 466         lhaname=lhapath(1:lhapathlen)//'/MRST2004FF3nlo.LHpdf'
 467      elseif((lhainput .ge. 20408) .and. (lhainput .le. 20409)) then
 468         lhaset = 20408
 469         lhaname=lhapath(1:lhapathlen)//'/MRST2004FF4nlo.LHpdf'
 470      elseif((lhainput .ge. 20450) .and. (lhainput .le. 20451)) then
 471         lhaset = 20450
 472         lhaname=lhapath(1:lhapathlen)//'/MRST2004nlo.LHgrid'
 473      elseif((lhainput .ge. 20452) .and. (lhainput .le. 20453)) then
 474         lhaset = 20452
 475         lhaname=lhapath(1:lhapathlen)//'/MRST2004FF3lo.LHgrid'
 476      elseif((lhainput .ge. 20454) .and. (lhainput .le. 20455)) then
 477         lhaset = 20454
 478         lhaname=lhapath(1:lhapathlen)//'/MRST2004FF4lo.LHgrid'
 479      elseif((lhainput .ge. 20456) .and. (lhainput .le. 20457)) then
 480         lhaset = 20456
 481         lhaname=lhapath(1:lhapathlen)//'/MRST2004FF3nlo.LHgrid'
 482      elseif((lhainput .ge. 20458) .and. (lhainput .le. 20459)) then
 483         lhaset = 20458
 484         lhaname=lhapath(1:lhapathlen)//'/MRST2004FF4nlo.LHgrid'
 485      elseif((lhainput .ge. 20460) .and. (lhainput .le. 20462)) then
 486         lhaset = 20460
 487         lhaname=lhapath(1:lhapathlen)//'/MRST2004qed.LHgrid'
 488      elseif((lhainput .ge. 20470) .and. (lhainput .le. 20471)) then
 489         lhaset = 20470
 490         lhaname=lhapath(1:lhapathlen)//'/MRST2004nnlo.LHgrid'
 491      elseif((lhainput .ge. 20550) .and. (lhainput .le. 20580)) then
 492         lhaset = 20550
 493         lhaname=lhapath(1:lhapathlen)//'/MRST2006nnlo.LHgrid'
 494         q2min = 1.0d0
 495         q2max = 1.0d09
 496         xmin = 1.0d-6
 497      elseif((lhainput .ge. 20650) .and. (lhainput .le. 20650)) then
 498         lhaset = 20650
 499         lhaname=lhapath(1:lhapathlen)//'/MRST2007lomod.LHgrid'
 500      elseif((lhainput .ge. 20651) .and. (lhainput .le. 20651)) then
 501         lhaset = 20651
 502         lhaname=lhapath(1:lhapathlen)//'/MRSTMCal.LHgrid'
 503      elseif((lhainput .ge. 29000) .and. (lhainput .le. 29003)) then
 504         lhaset = 29000
 505         lhaname=lhapath(1:lhapathlen)//'/MRST98.LHpdf'
 506      elseif((lhainput .ge. 29040) .and. (lhainput .le. 29045)) then
 507         lhaset = 29040
 508         lhaname=lhapath(1:lhapathlen)//'/MRST98lo.LHgrid'
 509      elseif((lhainput .ge. 29050) .and. (lhainput .le. 29055)) then
 510         lhaset = 29050
 511         lhaname=lhapath(1:lhapathlen)//'/MRST98nlo.LHgrid'
 512      elseif((lhainput .ge. 29060) .and. (lhainput .le. 29065)) then
 513         lhaset = 29060
 514         lhaname=lhapath(1:lhapathlen)//'/MRST98dis.LHgrid'
 515      elseif((lhainput .ge. 29070) .and. (lhainput .le. 29071)) then
 516         lhaset = 29070
 517         lhaname=lhapath(1:lhapathlen)//'/MRST98ht.LHgrid'
 518      else
 519         write(lhaprint,5150)  lhaset
 520         stop
 521      endif
 522      ! Fermi Family
 523   elseif((lhainput .ge. 30000) .and. (lhainput .le. 39999)) then
 524      if((lhainput .ge. 30100) .and. (lhainput .le. 30200)) then
 525         lhaset = 30100
 526         lhaname=lhapath(1:lhapathlen)//'/Fermi2002_100.LHpdf'
 527      elseif((lhainput .ge. 31000) .and. (lhainput .le. 32000)) then
 528         lhaset = 31000
 529         lhaname=lhapath(1:lhapathlen)//'/Fermi2002_1000.LHpdf'
 530      else
 531         write(lhaprint,5150)  lhaset
 532         stop
 533      endif
 534      ! Alekhin Family
 535   elseif((lhainput .ge. 40000) .and. (lhainput .le. 49999)) then
 536      if((lhainput .ge. 40100) .and. (lhainput .le. 40200)) then
 537         lhaset = 40100
 538         lhaname=lhapath(1:lhapathlen)//'/Alekhin_100.LHpdf'
 539      elseif((lhainput .ge. 41000) .and. (lhainput .le. 41999)) then
 540         lhaset = 41000
 541         lhaname=lhapath(1:lhapathlen)//'/Alekhin_1000.LHpdf'
 542      elseif((lhainput .ge. 40350) .and. (lhainput .le. 40367)) then
 543         lhaset = 40350
 544         lhaname=lhapath(1:lhapathlen)//'/a02m_lo.LHgrid'
 545         xmin = 1.0d-7
 546         q2min = 0.8d0
 547         q2max = 2.0d08
 548      elseif((lhainput .ge. 40450) .and. (lhainput .le. 40467)) then
 549         lhaset = 40450
 550         lhaname=lhapath(1:lhapathlen)//'/a02m_nlo.LHgrid'
 551         xmin = 1.0d-7
 552         q2min = 0.8d0
 553         q2max = 2.0d08
 554      elseif((lhainput .ge. 40550) .and. (lhainput .le. 40567)) then
 555         lhaset = 40550
 556         lhaname=lhapath(1:lhapathlen)//'/a02m_nnlo.LHgrid'
 557         xmin = 1.0d-7
 558         q2min = 0.8d0
 559         q2max = 2.0d08
 560      else
 561         write(lhaprint,5150)  lhaset
 562         stop
 563      endif
 564      ! Botje Family
 565   elseif((lhainput .ge. 50000) .and. (lhainput .le. 59999)) then
 566      if((lhainput .ge. 50100) .and. (lhainput .le. 50200)) then
 567         lhaset = 50100
 568         lhaname=lhapath(1:lhapathlen)//'/Botje_100.LHpdf'
 569      elseif((lhainput .ge. 51000) .and. (lhainput .le. 51999)) then
 570         lhaset = 51000
 571         lhaname=lhapath(1:lhapathlen)//'/Botje_1000.LHpdf'
 572      else
 573         write(lhaprint,5150)  lhaset
 574         stop
 575      endif
 576      ! ZEUS Family
 577   elseif((lhainput .ge. 60000) .and. (lhainput .le. 69999)) then
 578      q2min = 0.3d0
 579      q2max = 2.0d05
 580      if((lhainput .ge. 60000) .and. (lhainput .le. 60022)) then
 581         lhaset = 60000
 582         lhaname=lhapath(1:lhapathlen)//'/ZEUS2002_TR.LHpdf'
 583      elseif((lhainput .ge. 60100) .and. (lhainput .le. 60122)) then
 584         lhaset = 60100
 585         lhaname=lhapath(1:lhapathlen)//'/ZEUS2002_ZM.LHpdf'
 586      elseif((lhainput .ge. 60200) .and. (lhainput .le. 60222)) then
 587         lhaset = 60200
 588         lhaname=lhapath(1:lhapathlen)//'/ZEUS2002_FF.LHpdf'
 589      elseif((lhainput .ge. 60300) .and. (lhainput .le. 60322)) then
 590         lhaset = 60300
 591         lhaname=lhapath(1:lhapathlen)//'/ZEUS2005_ZJ.LHpdf'
 592      else
 593         write(lhaprint,5150)  lhaset
 594         stop
 595      endif
 596      ! H1 Family
 597   elseif((lhainput .ge. 70000) .and. (lhainput .le. 79999)) then
 598      q2min = 1.5d0
 599      q2max = 3.5d04
 600      xmin = 5.7d-5
 601      if((lhainput .ge. 70050) .and. (lhainput .le. 70050)) then
 602         lhaset = 70050
 603         lhaname=lhapath(1:lhapathlen)//'/H12000ms.LHgrid'
 604      elseif((lhainput .ge. 70051) .and. (lhainput .le. 70070)) then
 605         lhaset = 70050
 606         lhaname=lhapath(1:lhapathlen)//'/H12000msE.LHgrid'
 607      elseif((lhainput .ge. 70150) .and. (lhainput .le. 70150)) then
 608         lhaset = 70150
 609         lhaname=lhapath(1:lhapathlen)//'/H12000dis.LHgrid'
 610      elseif((lhainput .ge. 70151) .and. (lhainput .le. 70170)) then
 611         lhaset = 70150
 612         lhaname=lhapath(1:lhapathlen)//'/H12000disE.LHgrid'
 613      elseif((lhainput .ge. 70250) .and. (lhainput .le. 70250)) then
 614         lhaset = 70250
 615         lhaname=lhapath(1:lhapathlen)//'/H12000lo.LHgrid'
 616      elseif((lhainput .ge. 70251) .and. (lhainput .le. 70270)) then
 617         lhaset = 70250
 618         lhaname=lhapath(1:lhapathlen)//'/H12000loE.LHgrid'
 619         ! Temporarily removed on returning to original H12000 files
 620         ! elseif((lhainput .ge. 70350) .and. (lhainput .le. 70350)) then
 621         ! lhaset = 70350
 622         ! lhaname=lhapath(1:lhapathlen)//'/H12000lo2.LHgrid'
 623         ! elseif((lhainput .ge. 70351) .and. (lhainput .le. 70370)) then
 624         ! lhaset = 70350
 625         ! lhaname=lhapath(1:lhapathlen)//'/H12000lo2E.LHgrid'
 626      else
 627         write(lhaprint,5150)  lhaset
 628         stop
 629      endif
 630      ! GRV Family
 631   elseif((lhainput .ge. 80000) .and. (lhainput .le. 89999)) then
 632      q2min = 0.8d0
 633      q2max = 2.0d06
 634      xmin = 1.0d-9
 635      if((lhainput .ge. 80050) .and. (lhainput .le. 80051)) then
 636         lhaset = 80050
 637         lhaname=lhapath(1:lhapathlen)//'/GRV98nlo.LHgrid'
 638      elseif((lhainput .ge. 80060) .and. (lhainput .le. 80060)) then
 639         lhaset = 80060
 640         lhaname=lhapath(1:lhapathlen)//'/GRV98lo.LHgrid'
 641      else
 642         write(lhaprint,5150)  lhaset
 643         stop
 644      endif
 645      !
 646      ! Pions
 647      !
 648      ! OW-PI Family
 649   elseif((lhainput .ge. 210) .and. (lhainput .le. 212)) then
 650      q2min = 4.0d0
 651      q2max = 2.0d03
 652      xmin = 5.0d-03
 653      xmax = 0.9998d0
 654      if((lhainput .ge. 210) .and. (lhainput .le. 212)) then
 655         lhaset = 210
 656         lhaname=lhapath(1:lhapathlen)//'/OWPI.LHgrid'
 657      else
 658         write(lhaprint,5150)  lhaset
 659         stop
 660      endif
 661      ! SMRS-PI Family
 662   elseif((lhainput .ge. 230) .and. (lhainput .le. 233)) then
 663      q2min = 5.0d0
 664      q2max = 1.31d06
 665      xmin = 1.0d-05
 666      xmax = 0.9998d0
 667      if((lhainput .ge. 230) .and. (lhainput .le. 233)) then
 668         lhaset = 230
 669         lhaname=lhapath(1:lhapathlen)//'/SMRSPI.LHgrid'
 670      else
 671         write(lhaprint,5150)  lhaset
 672         stop
 673      endif
 674      ! GRV-PI Family
 675   elseif((lhainput .ge. 250) .and. (lhainput .le. 252)) then
 676      q2max = 1.00d06
 677      xmin = 1.0d-05
 678      xmax = 0.9998d0
 679      if((lhainput .ge. 250) .and. (lhainput .le. 251)) then
 680         q2min = 3.0d-1
 681         lhaset = 250
 682         lhaname=lhapath(1:lhapathlen)//'/GRVPI1.LHgrid'
 683      elseif((lhainput .ge. 252) .and. (lhainput .le. 252)) then
 684         q2min = 2.5d-1
 685         lhaset = 252
 686         lhaname=lhapath(1:lhapathlen)//'/GRVPI0.LHgrid'
 687      else
 688         write(lhaprint,5150)  lhaset
 689         stop
 690      endif
 691      ! ABFKW-PI Family
 692   elseif((lhainput .ge. 260) .and. (lhainput .le. 263)) then
 693      q2min = 2.0d0
 694      q2max = 1.00d08
 695      xmin = 1.0d-03
 696      xmax = 0.9998d0
 697      if((lhainput .ge. 260) .and. (lhainput .le. 263)) then
 698         lhaset = 260
 699         lhaname=lhapath(1:lhapathlen)//'/ABFKWPI.LHgrid'
 700      else
 701         write(lhaprint,5150)  lhaset
 702         stop
 703      endif
 704      !
 705      ! photons
 706      !
 707      ! DO-G Family
 708   elseif((lhainput .ge. 310) .and. (lhainput .le. 312)) then
 709      q2min = 1.0d01
 710      q2max = 1.00d04
 711      xmin = 1.0d-05
 712      xmax = 0.9d0
 713      if((lhainput .ge. 310) .and. (lhainput .le. 311)) then
 714         lhaset = 310
 715         lhaname=lhapath(1:lhapathlen)//'/DOG0.LHgrid'
 716      elseif((lhainput .ge. 312) .and. (lhainput .le. 312)) then
 717         lhaset = 312
 718         lhaname=lhapath(1:lhapathlen)//'/DOG1.LHgrid'
 719      else
 720         write(lhaprint,5150)  lhaset
 721         stop
 722      endif
 723      ! DG-G Family
 724   elseif((lhainput .ge. 320) .and. (lhainput .le. 324)) then
 725      xmin = 1.0d-05
 726      xmax = 0.9998d0
 727      lhaset = 320
 728      if((lhainput .ge. 320) .and. (lhainput .le. 321)) then
 729         q2min = 1.0d0
 730         q2max = 1.0d04
 731         ! lhaset = 320
 732         lhaname=lhapath(1:lhapathlen)//'/DGG.LHgrid'
 733      elseif((lhainput .ge. 322) .and. (lhainput .le. 322)) then
 734         q2min = 1.0d0
 735         q2max = 5.0d01
 736         ! lhaset = 322
 737         lhaname=lhapath(1:lhapathlen)//'/DGG.LHgrid'
 738      elseif((lhainput .ge. 323) .and. (lhainput .le. 323)) then
 739         q2min = 2.0d1
 740         q2max = 5.0d02
 741         ! lhaset = 323
 742         lhaname=lhapath(1:lhapathlen)//'/DGG.LHgrid'
 743      elseif((lhainput .ge. 324) .and. (lhainput .le. 324)) then
 744         q2min = 2.0d2
 745         q2max = 1.0d04
 746         ! lhaset = 324
 747         lhaname=lhapath(1:lhapathlen)//'/DGG.LHgrid'
 748      else
 749         write(lhaprint,5150)  lhaset
 750         stop
 751      endif
 752      ! LAC/GAL-G Family
 753   elseif((lhainput .ge. 330) .and. (lhainput .le. 334)) then
 754      q2min = 4.0d00
 755      q2max = 1.0d05
 756      xmin = 1.0d-04
 757      xmax = 0.9998d0
 758      lhaset = 330
 759      if((lhainput .ge. 330) .and. (lhainput .le. 332)) then
 760         ! lhaset = 330
 761         lhaname=lhapath(1:lhapathlen)//'/LACG.LHgrid'
 762      elseif((lhainput .ge. 333) .and. (lhainput .le. 333)) then
 763         q2min = 1.0d00
 764         ! lhaset = 333
 765         lhaname=lhapath(1:lhapathlen)//'/LACG.LHgrid'
 766      elseif((lhainput .ge. 334) .and. (lhainput .le. 334)) then
 767         q2min = 4.0d00
 768         ! lhaset = 334
 769         lhaname=lhapath(1:lhapathlen)//'/LACG.LHgrid'
 770      else
 771         write(lhaprint,5150)  lhaset
 772         stop
 773      endif
 774      ! GSG/GSG96-G Family
 775   elseif((lhainput .ge. 340) .and. (lhainput .le. 345)) then
 776      q2min = 5.3d00
 777      q2max = 1.0d08
 778      xmin = 5.0d-04
 779      xmax = 0.9998d0
 780      if((lhainput .ge. 340) .and. (lhainput .le. 341)) then
 781         lhaset = 340
 782         lhaname=lhapath(1:lhapathlen)//'/GSG1.LHgrid'
 783      elseif((lhainput .ge. 342) .and. (lhainput .le. 343)) then
 784         lhaset = 342
 785         lhaname=lhapath(1:lhapathlen)//'/GSG0.LHgrid'
 786      elseif((lhainput .ge. 344) .and. (lhainput .le. 344)) then
 787         lhaset = 344
 788         lhaname=lhapath(1:lhapathlen)//'/GSG961.LHgrid'
 789      elseif((lhainput .ge. 345) .and. (lhainput .le. 345)) then
 790         lhaset = 345
 791         lhaname=lhapath(1:lhapathlen)//'/GSG960.LHgrid'
 792      else
 793         write(lhaprint,5150)  lhaset
 794         stop
 795      endif
 796      ! GRV-G Family
 797   elseif((lhainput .ge. 350) .and. (lhainput .le. 354)) then
 798      q2min = 3.0d-1
 799      q2max = 1.0d06
 800      xmin = 1.0d-05
 801      xmax = 0.9998d0
 802      if((lhainput .ge. 350) .and. (lhainput .le. 352)) then
 803         lhaset = 350
 804         lhaname=lhapath(1:lhapathlen)//'/GRVG1.LHgrid'
 805      elseif((lhainput .ge. 353) .and. (lhainput .le. 353)) then
 806         q2min = 2.5d-1
 807         lhaset = 353
 808         lhaname=lhapath(1:lhapathlen)//'/GRVG0.LHgrid'
 809      elseif((lhainput .ge. 354) .and. (lhainput .le. 354)) then
 810         q2min = 6.0d-1
 811         q2max = 5.0d04
 812         lhaset = 354
 813         lhaname=lhapath(1:lhapathlen)//'/GRVG0.LHgrid'
 814      else
 815         write(lhaprint,5150)  lhaset
 816         stop
 817      endif
 818      ! ACFGP-G Family
 819   elseif((lhainput .ge. 360) .and. (lhainput .le. 363)) then
 820      q2min = 2.0d00
 821      q2max = 5.5d05
 822      xmin = 1.37d-03
 823      xmax = 0.9998d0
 824      if((lhainput .ge. 360) .and. (lhainput .le. 363)) then
 825         lhaset = 360
 826         lhaname=lhapath(1:lhapathlen)//'/ACFGPG.LHgrid'
 827      else
 828         write(lhaprint,5150)  lhaset
 829         stop
 830      endif
 831      ! WHIT-G Family
 832   elseif((lhainput .ge. 380) .and. (lhainput .le. 386)) then
 833      q2min = 4.0d00
 834      q2max = 2.5d03
 835      xmin = 1.0d-03
 836      xmax = 0.9998d0
 837      if((lhainput .ge. 380) .and. (lhainput .le. 386)) then
 838         lhaset = 380
 839         lhaname=lhapath(1:lhapathlen)//'/WHITG.LHgrid'
 840      else
 841         write(lhaprint,5150)  lhaset
 842         stop
 843      endif
 844      ! SAS-G Family
 845   elseif ((lhainput .ge. 390) .and. (lhainput .le. 398)) then
 846      q2max = 5.0d04
 847      xmin = 1.0d-05
 848      xmax = 0.9998d0
 849      lhaset = 390
 850      if ((lhainput .ge. 390) .and. (lhainput .le. 392)) then
 851         q2min = 3.6d-1
 852         ! lhaset = 390
 853         lhaname=lhapath(1:lhapathlen)//'/SASG.LHgrid'
 854      elseif((lhainput .ge. 393) .and. (lhainput .le. 394)) then
 855         q2min = 4.0d00
 856         !           lhaset = 393
 857         lhaname=lhapath(1:lhapathlen)//'/SASG.LHgrid'
 858      elseif((lhainput .ge. 395) .and. (lhainput .le. 396)) then
 859         q2min = 3.6d-1
 860         !           lhaset = 395
 861         lhaname=lhapath(1:lhapathlen)//'/SASG.LHgrid'
 862      elseif((lhainput .ge. 397) .and. (lhainput .le. 398)) then
 863         q2min = 4.0d00
 864         !           lhaset = 397
 865         lhaname=lhapath(1:lhapathlen)//'/SASG.LHgrid'
 866      else
 867         write(lhaprint,5150)  lhaset
 868         stop
 869      endif
 870      ! Unknown Family ?! Giving up
 871   else
 872      write(lhaprint,5150)  lhaset
 873      stop
 874   endif
 875
 876   lhamemb=lhainput-lhaset
 877   ! Now work out if we have already called this set/member
 878   iset = 0
 879   do j=1,nsets
 880      if (lhaname.eq.lhanames(j).and. &
 881           lhamemb.eq.lhamembers(j)) then
 882         iset = j
 883      endif
 884   enddo
 885   if (iset.eq.0) then
 886      nsets=nsets+1
 887      if (nsets.gt.nmxset) then
 888         if (LHASILENT.ne.1) then
 889            print *, "WARNING: too many sets initialised"
 890            print *,"overwriting from set 1 again"
 891         endif
 892         nsets = 1
 893         ! stop
 894      endif
 895      iset=nsets
 896      lhanames(iset)=lhaname
 897      lhanumbers(iset)=lhainput
 898      lhamembers(iset)=lhamemb
 899      xxmin(iset)=xmin
 900      xxmax(iset)=xmax
 901      qq2min(iset)=q2min
 902      qq2max(iset)=q2max
 903      call initpdfsetm(iset,lhaname)
 904      call numberpdfm(iset,lhaallmem)
 905      if(lhasilent .ne. 1) then
 906         write(lhaprint,5151)
 907         write(lhaprint,5152) lhaname
 908         write(lhaprint,5153) lhaallmem
 909         write(lhaprint,5154)
 910      endif
 911      if ((lhamemb.lt.0) .or. (lhamemb.gt.lhaallmem)) then
 912         write(lhaprint,5155)  lhamemb
 913         write(lhaprint,5156)  lhaallmem
 914         stop
 915      endif
 916
 917      ! print *,'calling initpdf',lhamemb
 918      ! print *,'calling initpdfm ',iset,lhaname,lhamemb
 919      ! print *,'LHAGLUE .... initializing set,member ',iset,lhamemb
 920      call initpdfm(iset,lhamemb)
 921   endif
 922   !  the rest is done every time pdfset is called
 923   !print *,'setting nset to:',iset
 924   call setnset(iset)
 925   call setnmem(iset,lhamemb)
 926   xmin = xxmin(iset)
 927   xmax = xxmax(iset)
 928   q2min=qq2min(iset)
 929   q2max=qq2max(iset)
 930   call GetLam4M(iset,LHAMEMB,qcdlha4)
 931   call GetLam5M(iset,LHAMEMB,qcdlha5)
 932
 933   QMZ = 91.1876D0
 934   alphasLHA = alphasPDFM(iset,QMZ)
 935   if(lhasilent .ne. 1) write(lhaprint,5158) alphasLHA
 936
 937   if(lhaparm(17).EQ.'LHAPDF') then
 938      nptypepdfl = 1      ! Proton PDFs
 939      nflpdfl = 4
 940      qcdl4 = qcdlha4
 941      qcdl5 = qcdlha5
 942      if (LHASILENT .NE. 1) write(lhaprint,5159) qcdl4, qcdl5
 943   else
 944      nptypepdfl = 1      ! Proton PDFs
 945      nflpdfl = 4
 946      alambda = 0.192d0
 947      qcdlha4 = alambda
 948      qcdlha5 = alambda
 949      if (parm(1).EQ.'NPTYPE') then        !  PYTHIA
 950         qcdl4 = alambda
 951         qcdl5 = alambda
 952      endif
 953   endif
 954
 955   if(parm(1).EQ.'NPTYPE') then  !  herwig
 956      axmin  = xmin
 957      axmax  = xmax
 958      aq2min = q2min
 959      aq2max = q2max
 960   endif
 961   ! Formats for initialization information.
 962 5150 format(1X,'WRONG LHAPDF set number =',I12,' given! STOP EXE!')
 963 5151 format(1X,'==============================================')
 964 5152 format(1X,'PDFset name ',A80)
 965 5153 format(1X,'with ',I10,' members')
 966 5154 format(1X,'====  initialized. ===========================')
 967 5155 format(1X,'LHAPDF problem => YOU asked for member = ',I10)
 968 5156 format(1X,'Valid range is: 0 - ',I10,' Execution stopped.')
 969   !5157 format(1X,'Number of flavors for PDF is:',I4)
 970 5158 format(1X,'Strong coupling at Mz for PDF is:',F9.5)
 971 5159 format(1X,'Will use for PYTHIA QCDL4, QCDL5:',2F9.5)
 972
 973   return
 974 end subroutine pdfset
 975
 976
 977 !********************************************************************
 978 ! -- STRUCTA
 979 ! -- copy of PDFLIB to use the eks98 nuclear correction factors
 980
 981 subroutine structa(x,q,a,upv,dnv,usea,dsea,str,chm,bot,top,glu)
 982   implicit double precision (a-h,o-z)
 983   character*20 lparm
 984   call getlhaparm(15,lparm)
 985
 986   if(lparm.eq.'EPS08') then
 987      call eps08(x,q,a,ruv,rdv,ru,rd,rs,rc,rb,rt,rg)
 988   else if (lparm.eq.'EPS09LO') then
 989      call eps09(1, 1, a, x, q, ruv, rdv, ru, rd, rs, rc, rb, rt, rg)
 990   else if (lparm.eq.'EPS09NLO') then
 991      call eps09(2, 1, a, x, q, ruv, rdv, ru, rd, rs, rc, rb, rt, rg)
 992   else
 993      call eks98(x,q,a,ruv,rdv,ru,rd,rs,rc,rb,rt,rg)
 994   endif
 995   call structm(x,q,upv,dnv,usea,dsea,str,chm,bot,top,glu)
 996   upv = ruv*upv
 997   dnv = rdv*dnv
 998   usea = ru*usea
 999   dsea = rd*dsea
1000   str = rs*str
1001   chm = rc*chm
1002   bot = rb*bot
1003   top = rt*top
1004   glu = rg*glu
1005   return
1006 end subroutine structa
1007
1008
1009 !*********************************************************************
1010 ! STRUCTM
1011 ! Gives parton distributions according to the LHAPDF interface.
1012 ! Two evolution codes used:
1013 !   EVLCTEQ for CTEQ PDF sets
1014 !   QCDNUM  for Other PDF sets
1015 !
1016 ! Author: Dimitri Bourilkov  bourilkov@mailaps.org
1017 !
1018 ! v4.0  21-Mar-2005  Photon/pion/new p PDFs, updated for LHAPDF v4
1019 ! v3.0  23-Jan-2004
1020 !
1021 ! interface to LHAPDF library
1022 subroutine structm(dx,dq,upv,dnv,usea,dsea,str,chm,bot,top,glu)
1023   ! double precision and integer declarations.
1024   implicit double precision(a-h, o-z)
1025   implicit integer(i-n)
1026   ! Common blocks
1027   include 'commonlhapdf.inc'
1028   include 'commonlhasets.inc'
1029   include 'commonlhacontrol.inc'
1030   include 'commonlhaglsta.inc'
1031   ! interface to lhapdflib.
1032   double precision qcdlha4, qcdlha5
1033   integer nfllha
1034   common/lhapdfr/qcdlha4, qcdlha5, nfllha
1035   save /lhapdfr/
1036   integer lhaextrp
1037   common/lhapdfe/lhaextrp
1038   save /lhapdfe/
1039   ! interface to pdflib.
1040   common/w50513/xmin,xmax,q2min,q2max
1041   save /w50513/
1042   double precision xmin,xmax,q2min,q2max
1043   ! local variables
1044   double precision upv,dnv,usea,dsea,str,chm,bot,top,glu
1045   double precision dx,dq,x,q,f(-6:6),photon,gluino
1046
1047   x = dx
1048   q = dq
1049   q2 = q**2
1050   ! statistics
1051   if(lhaparm(16).ne.'NOSTAT') then
1052      totnum = totnum+1.d0
1053      if(x .lt. xmin) xminnum = xminnum+1.d0
1054      if(x .gt. xmax) xmaxnum = xmaxnum+1.d0
1055      if(q2 .lt. q2min) q2minnum = q2minnum+1.d0
1056      if(q2 .gt. q2max) q2maxnum = q2maxnum+1.d0
1057   endif
1058
1059   ! range of validity e.g. 10^-6 < x < 1, q2min < q^2 extended by
1060   ! freezing x*f(x,q2) at borders.
1061   if(lhaextrp .ne. 1) then    ! safe mode == "freeze"
1062      xin=max(xmin,min(xmax,x))
1063      q=sqrt(max(0d0,q2min,min(q2max,q2)))
1064   else                        ! adventurous mode == own risk !
1065      xin=x
1066   endif
1067
1068   call getnset(iset)
1069   !print *,'calling evolvepdfm:',iset
1070
1071   ! fix to allow STRUCTM to work for photon PDFs (Herwig does this)
1072   ! set P2 = 0.0d0 and IP2 = 0
1073   if(lhanumbers(iset).ge.300.and.lhanumbers(iset).le.399) then
1074      p2 = 0.0d0
1075      ip2 = 0
1076      call evolvepdfpm(iset,xin,q,p2,ip2,f)
1077   else if (lhanumbers(iset).ge.20460.and.lhanumbers(iset).le.20462) then
1078      call evolvepdfphotonm(iset,xin,q,f,photon)
1079   else if (lhanumbers(iset).ge.20670.and.lhanumbers(iset).le.20677) then
1080      call evolvepdfgluinom(iset,xin,q,f,gluino)
1081   else
1082      call evolvepdfm(iset,xin,q,f)
1083   endif
1084   glu = f(0)
1085   dsea = f(-1)
1086   dnv = f(1) - dsea
1087   usea = f(-2)
1088   upv = f(2) - usea
1089   str = f(3)
1090   chm = f(4)
1091   bot = f(5)
1092   top = f(6)
1093
1094   return
1095 end subroutine structm
1096
1097
1098 !*********************************************************************
1099 ! STRUCTP
1100 ! Gives parton distributions according to the LHAPDF interface.
1101 ! Used for photons.
1102 !
1103 ! v4.0  21-Mar-2005  Photon/pion/new p PDFs, updated for LHAPDF v4
1104 !
1105 ! Interface to LHAPDF library
1106 subroutine structp(dx,dq2,p2,ip2,upv,dnv,usea,dsea,str,chm,bot,top,glu)
1107   ! Double precision and integer declarations.
1108   implicit double precision(a-h, o-z)
1109   implicit integer(i-n)
1110   ! Common blocks
1111   include 'parmsetup.inc'
1112   include 'commonlhapdf.inc'
1113   include 'commonlhacontrol.inc'
1114   include 'commonlhaglsta.inc'
1115   ! Interface to LHAPDFLIB.
1116   double precision qcdlha4, qcdlha5
1117   integer nfllha
1118   common/lhapdfr/qcdlha4, qcdlha5, nfllha
1119   save /lhapdfr/
1120   integer lhaextrp
1121   common/lhapdfe/lhaextrp
1122   save /lhapdfe/
1123   ! Interface to PDFLIB.
1124   common/w50513/xmin,xmax,q2min,q2max
1125   save /w50513/
1126   double precision xmin,xmax,q2min,q2max
1127   ! Local variables
1128   double precision upv,dnv,usea,dsea,str,chm,bot,top,glu
1129   double precision dx,dq2,q2,x,q,f(-6:6)
1130
1131   x = dx
1132   q2 = dq2
1133   ! Statistics
1134   if(lhaparm(16).ne.'NOSTAT') then
1135      totnup = totnup+1.d0
1136      if(x .lt. xmin) xminnup = xminnup+1.d0
1137      if(x .gt. xmax) xmaxnup = xmaxnup+1.d0
1138      if(q2 .lt. q2min) q2minnup = q2minnup+1.d0
1139      if(q2 .gt. q2max) q2maxnup = q2maxnup+1.d0
1140   endif
1141
1142   ! Range of validity e.g. 10^-6 < x < 1, Q2MIN < Q^2 extended by
1143   ! freezing x*f(x,Q2) at borders.
1144   q = dsqrt(q2)
1145   if(lhaextrp .ne. 1) then    ! safe mode == "freeze"
1146      xin=max(xmin,min(xmax,x))
1147      q=sqrt(max(0d0,q2min,min(q2max,q2)))
1148   else                        ! adventurous mode == OWN RISK !
1149      xin=x
1150   endif
1151   call getnset(iset)
1152   call evolvepdfpm(iset,xin,q,p2,ip2,f)
1153   glu = f(0)
1154   dsea = f(-1)
1155   dnv = f(1) - dsea
1156   usea = f(-2)
1157   upv = f(2) - usea
1158   str = f(3)
1159   chm = f(4)
1160   bot = f(5)
1161   top = f(6)
1162   return
1163 end subroutine structp
1164
1165
1166 !*********************************************************************
1167 ! PDFSTA
1168 ! For statistics ON structure functions (under/over-flows)
1169 !
1170 ! Author: Dimitri Bourilkov  bourilkov@mailaps.org
1171 !
1172 !
1173 ! first introduced in v4.0  28-Apr-2005
1174 !
1175 subroutine pdfsta
1176   ! Double precision and integer declarations.
1177   implicit double precision(a-h, o-z)
1178   implicit integer(i-n)
1179   ! Common blocks
1180   include 'commonlhaglsta.inc'
1181   ! Interface to LHAPDFLIB.
1182
1183   print *
1184   print *,'===== PDFSTA statistics for PDF under/over-flows ===='
1185   print *
1186   print *,'====== STRUCTM statistics for nucleon/pion PDFs ====='
1187   print *
1188   print *,'  total # of calls ',TOTNUM
1189   if(totnum .gt. 0.d0) then
1190      percbelow = 100.d0*xminnum/totnum
1191      percabove = 100.d0*xmaxnum/totnum
1192      print *,'  X  below PDF min ',xminnum,' or ',percbelow, ' %'
1193      print *,'  X  above PDF max ',xmaxnum,' or ',percabove, ' %'
1194      percbelow = 100.d0*q2minnum/totnum
1195      percabove = 100.d0*q2maxnum/totnum
1196      print *,'  Q2 below PDF min ',q2minnum,' or ',percbelow, ' %'
1197      print *,'  Q2 above PDF max ',q2maxnum,' or ',percabove, ' %'
1198   endif
1199   print *
1200   print *,'========= STRUCTP statistics for photon PDFs ========'
1201   print *
1202   print *,'  total # of calls ',totnup
1203   if(totnup .gt. 0.d0) then
1204      percbelow = 100.d0*xminnup/totnup
1205      percabove = 100.d0*xmaxnup/totnup
1206      print *,'  X  below PDF min ',xminnup,' or ',percbelow, ' %'
1207      print *,'  X  above PDF max ',xmaxnup,' or ',percabove, ' %'
1208      percbelow = 100.d0*q2minnup/totnup
1209      percabove = 100.d0*q2maxnup/totnup
1210      print *,'  Q2 below PDF min ',q2minnup,' or ',percbelow, ' %'
1211      print *,'  Q2 above PDF max ',q2maxnup,' or ',percabove, ' %'
1212   endif
1213   print *
1214   return
1215 end subroutine pdfsta
1216
1217
1218 subroutine pftopdg(dx,dscale,dxpdf)
1219   !include "pdf/expdp.inc"
1220   double precision dx,dscale,dupv,ddnv,dusea,ddsea,dstr,dchm,dbot,dtop,dgl,dxpdf(-6:6)
1221   ! Call STRUCTM in PDFLIB to get flavour content
1222   call structm(dx,dscale,dupv,ddnv,dusea,ddsea,dstr,dchm,dbot,dtop,dgl)
1223   ! Convert flavour convention of PDFLIB to PDG convention
1224   dxpdf(0) = dgl
1225   dxpdf(1) = ddnv + ddsea
1226   dxpdf(2) = dupv + dusea
1227   dxpdf(3) = dstr
1228   dxpdf(4) = dchm
1229   dxpdf(5) = dbot
1230   dxpdf(6) = dtop
1231   dxpdf(-1) = ddsea
1232   dxpdf(-2) = dusea
1233   dxpdf(-3) = dstr
1234   dxpdf(-4) = dchm
1235   dxpdf(-5) = dbot
1236   dxpdf(-6) = dtop
1237   return
1238 end subroutine pftopdg