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