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