LHAPDF/lhapdf5.2.2/lhaglue.f

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