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