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