Extracting the BLOCK DATA in a separate file. Changes to make it working on macosx
[u/mrichter/AliRoot.git] / 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 ***************************************************************