JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
mstwpdf.f
Go to the documentation of this file.
1 C----------------------------------------------------------------------
2 C-- Fortran interpolation code for MSTW PDFs, building on existing
3 C-- MRST Fortran code and Jeppe Andersen's C++ code.
4 C-- Three user interfaces:
5 C-- call GetAllPDFs(prefix,ih,x,q,upv,dnv,usea,dsea,
6 C-- str,sbar,chm,cbar,bot,bbar,glu,phot)
7 C-- call GetAllPDFsAlt(prefix,ih,x,q,xpdf,xphoton)
8 C-- xf = GetOnePDF(prefix,ih,x,q,f)
9 C-- See enclosed example.f for usage.
10 C-- Comments to Graeme Watt <Graeme.Watt(at)cern.ch>.
11 C-- Updated 25/06/2010: Enlarge allowed range for m_c and m_b.
12 C-- Updated 25/01/2011: Fix "NaN" bug for q <= m_c when m_c^2 < 1.25 GeV^2.
13 C----------------------------------------------------------------------
14 
15 C----------------------------------------------------------------------
16 
17 C-- Traditional MRST-like interface: return all flavours.
18 C-- (Note the additional "sbar", "cbar", "bbar" and "phot"
19 C-- compared to previous MRST releases.)
20  subroutine getallpdfs(prefix,ih,x,q,
21  & upv,dnv,usea,dsea,str,sbar,chm,cbar,bot,bbar,glu,phot)
22  implicit none
23  integer ih
24  double precision x,q,upv,dnv,usea,dsea,str,sbar,chm,cbar,
25  & bot,bbar,glu,phot,getonepdf,up,dn,sv,cv,bv
26  character*(*) prefix
27 
28 C-- Quarks.
29  dn = getonepdf(prefix,ih,x,q,1)
30  up = getonepdf(prefix,ih,x,q,2)
31  str = getonepdf(prefix,ih,x,q,3)
32  chm = getonepdf(prefix,ih,x,q,4)
33  bot = getonepdf(prefix,ih,x,q,5)
34 
35 C-- Valence quarks.
36  dnv = getonepdf(prefix,ih,x,q,7)
37  upv = getonepdf(prefix,ih,x,q,8)
38  sv = getonepdf(prefix,ih,x,q,9)
39  cv = getonepdf(prefix,ih,x,q,10)
40  bv = getonepdf(prefix,ih,x,q,11)
41 
42 C-- Antiquarks = quarks - valence quarks.
43  dsea = dn - dnv
44  usea = up - upv
45  sbar = str - sv
46  cbar = chm - cv
47  bbar = bot - bv
48 
49 C-- Gluon.
50  glu = getonepdf(prefix,ih,x,q,0)
51 
52 C-- Photon (= zero unless considering QED contributions).
53  phot = getonepdf(prefix,ih,x,q,13)
54 
55  return
56  end
57 
58 C----------------------------------------------------------------------
59 
60 C-- Alternative LHAPDF-like interface: return PDFs in an array.
61  subroutine getallpdfsalt(prefix,ih,x,q,xpdf,xphoton)
62  implicit none
63  integer ih,f
64  double precision x,q,xpdf(-6:6),xphoton,xvalence,GetOnePDF
65  character*(*) prefix
66 
67  do f = 1, 6
68  xpdf(f) = getonepdf(prefix,ih,x,q,f) ! quarks
69  xvalence = getonepdf(prefix,ih,x,q,f+6) ! valence quarks
70  xpdf(-f) = xpdf(f) - xvalence ! antiquarks
71  end do
72  xpdf(0) = getonepdf(prefix,ih,x,q,0) ! gluon
73  xphoton = getonepdf(prefix,ih,x,q,13) ! photon
74 
75  return
76  end
77 
78 C----------------------------------------------------------------------
79 
80 C-- Get only one parton flavour 'f', using PDG notation
81 C-- (apart from gluon has f=0, not 21):
82 C-- f = -6, -5, -4, -3, -2, -1,0,1,2,3,4,5,6
83 C-- = tbar,bbar,cbar,sbar,ubar,dbar,g,d,u,s,c,b,t.
84 C-- Can also get valence quarks directly:
85 C-- f = 7, 8, 9,10,11,12.
86 C-- = dv,uv,sv,cv,bv,tv.
87 C-- Photon: f = 13.
88  double precision function getonepdf(prefix,ih,x,q,f)
89  implicit none
90  logical warn,fatal
91  parameter(warn=.false.,fatal=.true.)
92 C-- Set warn=.true. to turn on warnings when extrapolating.
93 C-- Set fatal=.false. to return zero instead of terminating when
94 C-- invalid input values of x and q are used.
95  integer ih,f,nhess,nx,nq,np,nqc0,nqb0,n,m,ip,io,
96  & alphasorder,alphasnfmax,nextraflavours,lentrim
97  double precision x,q,xmin,xmax,qsqmin,qsqmax,mc2,mb2,eps,
98  & dummy,qsq,xlog,qsqlog,res,res1,anom,extrapolatepdf,
99  & interpolatepdf,distance,tolerance,
100  & mcharm,mbottom,alphasq0,alphasmz
101  parameter(nx=64,nq=48,np=12)
102  parameter(xmin=1d-6,xmax=1d0,qsqmin=1d0,qsqmax=1d9,eps=1d-6)
103  parameter(nhess=2*20)
104  character set*2,prefix*(*),filename*60,oldprefix(0:nhess)*50
105  character dummychar,dummyword*50
106  double precision ff(np,nx,nq)
107  double precision qqorig(nq),qq(nq),xx(nx),cc(np,0:nhess,nx,nq,4,4)
108  double precision xxl(nx),qql(nq)
109 C-- Store distance along each eigenvector, tolerance,
110 C-- heavy quark masses and alphaS parameters in COMMON block.
111  common/mstwcommon/distance,tolerance,
112  & mcharm,mbottom,alphasq0,alphasmz,alphasorder,alphasnfmax
113  save
114  data xx/1d-6,2d-6,4d-6,6d-6,8d-6,
115  & 1d-5,2d-5,4d-5,6d-5,8d-5,
116  & 1d-4,2d-4,4d-4,6d-4,8d-4,
117  & 1d-3,2d-3,4d-3,6d-3,8d-3,
118  & 1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
119  & .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
120  & .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
121  & .5d0,.525d0,.55d0,.575d0,.6d0,.625d0,.65d0,.675d0,
122  & .7d0,.725d0,.75d0,.775d0,.8d0,.825d0,.85d0,.875d0,
123  & .9d0,.925d0,.95d0,.975d0,1d0/
124  data qqorig/1.d0,
125  & 1.25d0,1.5d0,0.d0,0.d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,
126  & 1d1,1.2d1,0.d0,0.d0,2.6d1,4d1,6.4d1,1d2,
127  & 1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
128  & 1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
129  & 1.8d6,3.2d6,5.6d6,1d7,1.8d7,3.2d7,5.6d7,1d8,
130  & 1.8d8,3.2d8,5.6d8,1d9/
131 
132  if (f.lt.-6.or.f.gt.13) then
133  print *,"Error: invalid parton flavour = ",f
134  stop
135  end if
136 
137  if (ih.lt.0.or.ih.gt.nhess) then
138  print *,"Error: invalid eigenvector number = ",ih
139  stop
140  end if
141 
142 C-- Check if the requested parton set is already in memory.
143  if (oldprefix(ih).ne.prefix) then
144 
145 C-- Start of initialisation for eigenvector set "i" ...
146 C-- Do this only the first time the set "i" is called,
147 C-- OR if the prefix has changed from the last time.
148 
149 C-- Check that the character arrays "oldprefix" and "filename"
150 C-- are large enough.
151  if (lentrim(prefix).gt.len(oldprefix(ih))) then
152  print *,"Error in GetOnePDF: increase size of oldprefix"
153  stop
154  else if (lentrim(prefix)+7.gt.len(filename)) then
155  print *,"Error in GetOnePDF: increase size of filename"
156  stop
157  end if
158 
159  write(set,'(I2.2)') ih ! convert integer to string
160 C-- Remove trailing blanks from prefix before assigning filename.
161  filename = prefix(1:lentrim(prefix))//'.'//set//'.dat'
162 C-- Line below can be commented out if you don't want this message.
163  print *,"Reading PDF grid from ",filename(1:lentrim(filename))
164  open(unit=33,file=filename,iostat=io,status='old')
165  if (io.ne.0) then
166  print *,"Error in GetOnePDF: can't open ",
167  & filename(1:lentrim(filename))
168  stop
169  end if
170 
171 C-- Read header containing heavy quark masses and alphaS values.
172  read(33,*)
173  read(33,*)
174  read(33,*) dummychar,dummyword,dummyword,dummychar,
175  & distance,tolerance
176  read(33,*) dummychar,dummyword,dummychar,mcharm
177  read(33,*) dummychar,dummyword,dummychar,mbottom
178  read(33,*) dummychar,dummyword,dummychar,alphasq0
179  read(33,*) dummychar,dummyword,dummychar,alphasmz
180  read(33,*) dummychar,dummyword,dummyword,dummychar,
181  & alphasorder,alphasnfmax
182  read(33,*) dummychar,dummyword,dummychar,nextraflavours
183  read(33,*)
184  read(33,*)
185  mc2=mcharm**2
186  mb2=mbottom**2
187 
188 C-- Check that the heavy quark masses are sensible.
189 C-- Redistribute grid points if not in usual range.
190  do m=1,nq
191  qq(m) = qqorig(m)
192  end do
193  if (mc2.le.qq(1).or.mc2+eps.ge.qq(8)) then
194  print *,"Error in GetOnePDF: invalid mCharm = ",mcharm
195  stop
196  else if (mc2.lt.qq(2)) then
197  nqc0=2
198  qq(4)=qq(2)
199  qq(5)=qq(3)
200  else if (mc2.lt.qq(3)) then
201  nqc0=3
202  qq(5)=qq(3)
203  else if (mc2.lt.qq(6)) then
204  nqc0=4
205  else if (mc2.lt.qq(7)) then
206  nqc0=5
207  qq(4)=qq(6)
208  else
209  nqc0=6
210  qq(4)=qq(6)
211  qq(5)=qq(7)
212  end if
213  if (mb2.le.qq(12).or.mb2+eps.ge.qq(17)) then
214  print *,"Error in GetOnePDF: invalid mBottom = ",mbottom
215  stop
216  else if (mb2.lt.qq(13)) then
217  nqb0=13
218  qq(15)=qq(13)
219  else if (mb2.lt.qq(16)) then
220  nqb0=14
221  else
222  nqb0=15
223  qq(14)=qq(16)
224  end if
225  qq(nqc0)=mc2
226  qq(nqc0+1)=mc2+eps
227  qq(nqb0)=mb2
228  qq(nqb0+1)=mb2+eps
229 
230 C-- The nExtraFlavours variable is provided to aid compatibility
231 C-- with future grids where, for example, a photon distribution
232 C-- might be provided (cf. the MRST2004QED PDFs).
233  if (nextraflavours.lt.0.or.nextraflavours.gt.1) then
234  print *,"Error in GetOnePDF: invalid nExtraFlavours = ",
235  & nextraflavours
236  stop
237  end if
238 
239 C-- Now read in the grids from the grid file.
240  do n=1,nx-1
241  do m=1,nq
242  if (nextraflavours.gt.0) then
243  if (alphasorder.eq.2) then ! NNLO
244  read(33,'(12(1pe12.4))',iostat=io)
245  & (ff(ip,n,m),ip=1,12)
246  else ! LO or NLO
247  ff(10,n,m) = 0.d0 ! = chm-cbar
248  ff(11,n,m) = 0.d0 ! = bot-bbar
249  read(33,'(10(1pe12.4))',iostat=io)
250  & (ff(ip,n,m),ip=1,9),ff(12,n,m)
251  end if
252  else ! nExtraFlavours = 0
253  if (alphasorder.eq.2) then ! NNLO
254  ff(12,n,m) = 0.d0 ! = photon
255  read(33,'(11(1pe12.4))',iostat=io)
256  & (ff(ip,n,m),ip=1,11)
257  else ! LO or NLO
258  ff(10,n,m) = 0.d0 ! = chm-cbar
259  ff(11,n,m) = 0.d0 ! = bot-bbar
260  ff(12,n,m) = 0.d0 ! = photon
261  read(33,'(9(1pe12.4))',iostat=io)
262  & (ff(ip,n,m),ip=1,9)
263  end if
264  end if
265  if (io.ne.0) then
266  print *,"Error in GetOnePDF reading ",filename
267  stop
268  end if
269  enddo
270  enddo
271 
272 C-- Check that ALL the file contents have been read in.
273  read(33,*,iostat=io) dummy
274  if (io.eq.0) then
275  print *,"Error in GetOnePDF: not at end of ",filename
276  stop
277  end if
278  close(unit=33)
279 
280 C-- PDFs are identically zero at x = 1.
281  do m=1,nq
282  do ip=1,np
283  ff(ip,nx,m)=0d0
284  enddo
285  enddo
286 
287  do n=1,nx
288  xxl(n)=log10(xx(n))
289  enddo
290  do m=1,nq
291  qql(m)=log10(qq(m))
292  enddo
293 
294 C-- Initialise all parton flavours.
295  do ip=1,np
296  call initialisepdf(ip,np,ih,nhess,nx,nq,nqc0,nqb0,
297  & xxl,qql,ff,cc)
298  enddo
299 
300  oldprefix(ih) = prefix
301 
302 C-- ... End of initialisation for eigenvector set "ih".
303 
304  end if ! oldprefix(ih).ne.prefix
305 
306 C----------------------------------------------------------------------
307 
308  qsq=q*q
309 C-- If mc2 < qsq < mc2+eps, then qsq = mc2+eps.
310  if (qsq.gt.qq(nqc0).and.qsq.lt.qq(nqc0+1)) qsq = qq(nqc0+1)
311 C-- If mb2 < qsq < mb2+eps, then qsq = mb2+eps.
312  if (qsq.gt.qq(nqb0).and.qsq.lt.qq(nqb0+1)) qsq = qq(nqb0+1)
313 
314  xlog=log10(x)
315  qsqlog=log10(qsq)
316 
317  res = 0.d0
318 
319  if (f.eq.0) then ! gluon
320  ip = 1
321  else if (f.ge.1.and.f.le.5) then ! quarks
322  ip = f+1
323  else if (f.le.-1.and.f.ge.-5) then ! antiquarks
324  ip = -f+1
325  else if (f.ge.7.and.f.le.11) then ! valence quarks
326  ip = f
327  else if (f.eq.13) then ! photon
328  ip = 12
329  else if (abs(f).ne.6.and.f.ne.12) then
330  if (warn.or.fatal) print *,"Error in GetOnePDF: f = ",f
331  if (fatal) stop
332  end if
333 
334  if (x.le.0.d0.or.x.gt.xmax.or.q.le.0.d0) then
335 
336  if (warn.or.fatal) print *,"Error in GetOnePDF: x,qsq = ",
337  & x,qsq
338  if (fatal) stop
339 
340  else if (abs(f).eq.6.or.f.eq.12) then ! set top quarks to zero
341 
342  res = 0.d0
343 
344  else if (qsq.lt.qsqmin) then ! extrapolate to low Q^2
345 
346  if (warn) then
347  print *, "Warning in GetOnePDF, extrapolating: f = ",f,
348  & ", x = ",x,", q = ",q
349  end if
350 
351  if (x.lt.xmin) then ! extrapolate to low x
352 
353  res = extrapolatepdf(ip,np,ih,nhess,xlog,
354  & log10(qsqmin),nx,nq,xxl,qql,cc)
355  res1 = extrapolatepdf(ip,np,ih,nhess,xlog,
356  & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
357  if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
358  res = res - extrapolatepdf(ip+5,np,ih,nhess,xlog,
359  & log10(qsqmin),nx,nq,xxl,qql,cc)
360  res1 = res1 - extrapolatepdf(ip+5,np,ih,nhess,xlog,
361  & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
362  end if
363 
364  else ! do usual interpolation
365 
366  res = interpolatepdf(ip,np,ih,nhess,xlog,
367  & log10(qsqmin),nx,nq,xxl,qql,cc)
368  res1 = interpolatepdf(ip,np,ih,nhess,xlog,
369  & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
370  if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
371  res = res - interpolatepdf(ip+5,np,ih,nhess,xlog,
372  & log10(qsqmin),nx,nq,xxl,qql,cc)
373  res1 = res1 - interpolatepdf(ip+5,np,ih,nhess,xlog,
374  & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
375  end if
376 
377  end if
378 
379 C-- Calculate the anomalous dimension, dlog(xf)/dlog(qsq),
380 C-- evaluated at qsqmin. Then extrapolate the PDFs to low
381 C-- qsq < qsqmin by interpolating the anomalous dimenion between
382 C-- the value at qsqmin and a value of 1 for qsq << qsqmin.
383 C-- If value of PDF at qsqmin is very small, just set
384 C-- anomalous dimension to 1 to prevent rounding errors.
385 C-- Impose minimum anomalous dimension of -2.5.
386  if (abs(res).ge.1.d-5) then
387  anom = max( -2.5d0, (res1-res)/res/0.01d0 )
388  else
389  anom = 1.d0
390  end if
391  res = res*(qsq/qsqmin)**(anom*qsq/qsqmin+1.d0-qsq/qsqmin)
392 
393  else if (x.lt.xmin.or.qsq.gt.qsqmax) then ! extrapolate
394 
395  if (warn) then
396  print *, "Warning in GetOnePDF, extrapolating: f = ",f,
397  & ", x = ",x,", q = ",q
398  end if
399 
400  res = extrapolatepdf(ip,np,ih,nhess,xlog,
401  & qsqlog,nx,nq,xxl,qql,cc)
402 
403  if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
404  res = res - extrapolatepdf(ip+5,np,ih,nhess,xlog,
405  & qsqlog,nx,nq,xxl,qql,cc)
406  end if
407 
408  else ! do usual interpolation
409 
410  res = interpolatepdf(ip,np,ih,nhess,xlog,
411  & qsqlog,nx,nq,xxl,qql,cc)
412 
413  if (f.le.-1.and.f.ge.-5) then ! antiquark = quark - valence
414  res = res - interpolatepdf(ip+5,np,ih,nhess,xlog,
415  & qsqlog,nx,nq,xxl,qql,cc)
416  end if
417 
418  end if
419 
420  getonepdf = res
421 
422  return
423  end
424 
425 C----------------------------------------------------------------------
426 
427  subroutine initialisepdf(ip,np,ih,nhess,nx,my,myc0,myb0,
428  & xx,yy,ff,cc)
429  implicit none
430  integer nhess,ih,nx,my,myc0,myb0,j,k,l,m,n,ip,np
431  double precision xx(nx),yy(my),ff(np,nx,my),
432  & ff1(nx,my),ff2(nx,my),ff12(nx,my),ff21(nx,my),
433  & yy0(4),yy1(4),yy2(4),yy12(4),z(16),
434  & cl(16),cc(np,0:nhess,nx,my,4,4),iwt(16,16),
435  & polderiv1,polderiv2,polderiv3,d1,d2,d1d2,xxd
436 
437  data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
438  & 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
439  & -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
440  & 2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
441  & 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
442  & 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
443  & 0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
444  & 0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
445  & -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
446  & 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
447  & 9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
448  & -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
449  & 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
450  & 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
451  & -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
452  & 4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
453 
454  do m=1,my
455  ff1(1,m)=polderiv1(xx(1),xx(2),xx(3),
456  & ff(ip,1,m),ff(ip,2,m),ff(ip,3,m))
457  ff1(nx,m)=polderiv3(xx(nx-2),xx(nx-1),xx(nx),
458  & ff(ip,nx-2,m),ff(ip,nx-1,m),ff(ip,nx,m))
459  do n=2,nx-1
460  ff1(n,m)=polderiv2(xx(n-1),xx(n),xx(n+1),
461  & ff(ip,n-1,m),ff(ip,n,m),ff(ip,n+1,m))
462  enddo
463  enddo
464 
465 C-- Calculate the derivatives at qsq=mc2,mc2+eps,mb2,mb2+eps
466 C-- in a similar way as at the endpoints qsqmin and qsqmax.
467  do n=1,nx
468  do m=1,my
469  if (myc0.eq.2.and.m.eq.1) then
470  ff2(n,m)=(ff(ip,n,m+1)-ff(ip,n,m))/(yy(m+1)-yy(m))
471  else if (myc0.eq.2.and.m.eq.2) then
472  ff2(n,m)=(ff(ip,n,m)-ff(ip,n,m-1))/(yy(m)-yy(m-1))
473  else if (m.eq.1.or.m.eq.myc0+1.or.m.eq.myb0+1) then
474  ff2(n,m)=polderiv1(yy(m),yy(m+1),yy(m+2),
475  & ff(ip,n,m),ff(ip,n,m+1),ff(ip,n,m+2))
476  else if (m.eq.my.or.m.eq.myc0.or.m.eq.myb0) then
477  ff2(n,m)=polderiv3(yy(m-2),yy(m-1),yy(m),
478  & ff(ip,n,m-2),ff(ip,n,m-1),ff(ip,n,m))
479  else
480  ff2(n,m)=polderiv2(yy(m-1),yy(m),yy(m+1),
481  & ff(ip,n,m-1),ff(ip,n,m),ff(ip,n,m+1))
482  end if
483  end do
484  end do
485 
486 C-- Calculate the cross derivatives (d/dx)(d/dy).
487  do m=1,my
488  ff12(1,m)=polderiv1(xx(1),xx(2),xx(3),
489  & ff2(1,m),ff2(2,m),ff2(3,m))
490  ff12(nx,m)=polderiv3(xx(nx-2),xx(nx-1),xx(nx),
491  & ff2(nx-2,m),ff2(nx-1,m),ff2(nx,m))
492  do n=2,nx-1
493  ff12(n,m)=polderiv2(xx(n-1),xx(n),xx(n+1),
494  & ff2(n-1,m),ff2(n,m),ff2(n+1,m))
495  enddo
496  enddo
497 
498 C-- Calculate the cross derivatives (d/dy)(d/dx).
499  do n=1,nx
500  do m = 1, my
501  if (myc0.eq.2.and.m.eq.1) then
502  ff21(n,m)=(ff1(n,m+1)-ff1(n,m))/(yy(m+1)-yy(m))
503  else if (myc0.eq.2.and.m.eq.2) then
504  ff21(n,m)=(ff1(n,m)-ff1(n,m-1))/(yy(m)-yy(m-1))
505  else if (m.eq.1.or.m.eq.myc0+1.or.m.eq.myb0+1) then
506  ff21(n,m)=polderiv1(yy(m),yy(m+1),yy(m+2),
507  & ff1(n,m),ff1(n,m+1),ff1(n,m+2))
508  else if (m.eq.my.or.m.eq.myc0.or.m.eq.myb0) then
509  ff21(n,m)=polderiv3(yy(m-2),yy(m-1),yy(m),
510  & ff1(n,m-2),ff1(n,m-1),ff1(n,m))
511  else
512  ff21(n,m)=polderiv2(yy(m-1),yy(m),yy(m+1),
513  & ff1(n,m-1),ff1(n,m),ff1(n,m+1))
514  end if
515  end do
516  end do
517 
518 C-- Take the average of (d/dx)(d/dy) and (d/dy)(d/dx).
519  do n=1,nx
520  do m = 1, my
521  ff12(n,m)=0.5*(ff12(n,m)+ff21(n,m))
522  end do
523  end do
524 
525  do n=1,nx-1
526  do m=1,my-1
527  d1=xx(n+1)-xx(n)
528  d2=yy(m+1)-yy(m)
529  d1d2=d1*d2
530 
531  yy0(1)=ff(ip,n,m)
532  yy0(2)=ff(ip,n+1,m)
533  yy0(3)=ff(ip,n+1,m+1)
534  yy0(4)=ff(ip,n,m+1)
535 
536  yy1(1)=ff1(n,m)
537  yy1(2)=ff1(n+1,m)
538  yy1(3)=ff1(n+1,m+1)
539  yy1(4)=ff1(n,m+1)
540 
541  yy2(1)=ff2(n,m)
542  yy2(2)=ff2(n+1,m)
543  yy2(3)=ff2(n+1,m+1)
544  yy2(4)=ff2(n,m+1)
545 
546  yy12(1)=ff12(n,m)
547  yy12(2)=ff12(n+1,m)
548  yy12(3)=ff12(n+1,m+1)
549  yy12(4)=ff12(n,m+1)
550 
551  do k=1,4
552  z(k)=yy0(k)
553  z(k+4)=yy1(k)*d1
554  z(k+8)=yy2(k)*d2
555  z(k+12)=yy12(k)*d1d2
556  enddo
557 
558  do l=1,16
559  xxd=0.d0
560  do k=1,16
561  xxd=xxd+iwt(k,l)*z(k)
562  enddo
563  cl(l)=xxd
564  enddo
565  l=0
566  do k=1,4
567  do j=1,4
568  l=l+1
569  cc(ip,ih,n,m,k,j)=cl(l)
570  enddo
571  enddo
572  enddo
573  enddo
574  return
575  end
576 
577 C----------------------------------------------------------------------
578 
579  double precision function interpolatepdf(ip,np,ih,nhess,x,y,
580  & nx,my,xx,yy,cc)
581  implicit none
582  integer ih,nx,my,nhess,locx,l,m,n,ip,np
583  double precision xx(nx),yy(my),cc(np,0:nhess,nx,my,4,4),
584  & x,y,z,t,u
585 
586  n=locx(xx,nx,x)
587  m=locx(yy,my,y)
588 
589  t=(x-xx(n))/(xx(n+1)-xx(n))
590  u=(y-yy(m))/(yy(m+1)-yy(m))
591 
592  z=0.d0
593  do l=4,1,-1
594  z=t*z+((cc(ip,ih,n,m,l,4)*u+cc(ip,ih,n,m,l,3))*u
595  & +cc(ip,ih,n,m,l,2))*u+cc(ip,ih,n,m,l,1)
596  enddo
597 
598  interpolatepdf = z
599 
600  return
601  end
602 
603 C----------------------------------------------------------------------
604 
605  double precision function extrapolatepdf(ip,np,ih,nhess,x,y,
606  & nx,my,xx,yy,cc)
607  implicit none
608  integer ih,nx,my,nhess,locx,n,m,ip,np
609  double precision xx(nx),yy(my),cc(np,0:nhess,nx,my,4,4),
610  & x,y,z,f0,f1,z0,z1,interpolatepdf
611 
612  n=locx(xx,nx,x) ! 0: below xmin, nx: above xmax
613  m=locx(yy,my,y) ! 0: below qsqmin, my: above qsqmax
614 
615 C-- If extrapolation in small x only:
616  if (n.eq.0.and.m.gt.0.and.m.lt.my) then
617  f0 = interpolatepdf(ip,np,ih,nhess,xx(1),y,nx,my,xx,yy,cc)
618  f1 = interpolatepdf(ip,np,ih,nhess,xx(2),y,nx,my,xx,yy,cc)
619  if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
620  z = exp(log(f0)+(log(f1)-log(f0))/(xx(2)-xx(1))*(x-xx(1)))
621  else
622  z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
623  end if
624 C-- If extrapolation into large q only:
625  else if (n.gt.0.and.m.eq.my) then
626  f0 = interpolatepdf(ip,np,ih,nhess,x,yy(my),nx,my,xx,yy,cc)
627  f1 = interpolatepdf(ip,np,ih,nhess,x,yy(my-1),nx,my,xx,yy,cc)
628  if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
629  z = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
630  & (y-yy(my)))
631  else
632  z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
633  end if
634 C-- If extrapolation into large q AND small x:
635  else if (n.eq.0.and.m.eq.my) then
636  f0 = interpolatepdf(ip,np,ih,nhess,xx(1),yy(my),nx,my,xx,yy,cc)
637  f1 = interpolatepdf(ip,np,ih,nhess,xx(1),yy(my-1),nx,my,xx,yy,
638  & cc)
639  if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
640  z0 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
641  & (y-yy(my)))
642  else
643  z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
644  end if
645  f0 = interpolatepdf(ip,np,ih,nhess,xx(2),yy(my),nx,my,xx,yy,cc)
646  f1 = interpolatepdf(ip,np,ih,nhess,xx(2),yy(my-1),nx,my,xx,yy,
647  & cc)
648  if (f0.gt.1.d-3.and.f1.gt.1.d-3) then
649  z1 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
650  & (y-yy(my)))
651  else
652  z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
653  end if
654  if (z0.gt.1.d-3.and.z1.gt.1.d-3) then
655  z = exp(log(z0)+(log(z1)-log(z0))/(xx(2)-xx(1))*(x-xx(1)))
656  else
657  z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
658  end if
659  else
660  print *,"Error in ExtrapolatePDF"
661  stop
662  end if
663 
664  extrapolatepdf = z
665 
666  return
667  end
668 
669 C----------------------------------------------------------------------
670 
671  integer function locx(xx,nx,x)
672 C-- returns an integer j such that x lies inbetween xx(j) and xx(j+1).
673 C-- nx is the length of the array with xx(nx) the highest element.
674  implicit none
675  integer nx,jl,ju,jm
676  double precision x,xx(nx)
677  if(x.eq.xx(1)) then
678  locx=1
679  return
680  endif
681  if(x.eq.xx(nx)) then
682  locx=nx-1
683  return
684  endif
685  ju=nx+1
686  jl=0
687  1 if((ju-jl).le.1) go to 2
688  jm=(ju+jl)/2
689  if(x.ge.xx(jm)) then
690  jl=jm
691  else
692  ju=jm
693  endif
694  go to 1
695  2 locx=jl
696  return
697  end
698 
699 C----------------------------------------------------------------------
700 
701  double precision function polderiv1(x1,x2,x3,y1,y2,y3)
702 C-- returns the estimate of the derivative at x1 obtained by a
703 C-- polynomial interpolation using the three points (x_i,y_i).
704  implicit none
705  double precision x1,x2,x3,y1,y2,y3
706  polderiv1=(x3*x3*(y1-y2)+2.d0*x1*(x3*(-y1+y2)+x2*(y1-y3))
707  & +x2*x2*(-y1+y3)+x1*x1*(-y2+y3))/((x1-x2)*(x1-x3)*(x2-x3))
708  return
709  end
710 
711  double precision function polderiv2(x1,x2,x3,y1,y2,y3)
712 C-- returns the estimate of the derivative at x2 obtained by a
713 C-- polynomial interpolation using the three points (x_i,y_i).
714  implicit none
715  double precision x1,x2,x3,y1,y2,y3
716  polderiv2=(x3*x3*(y1-y2)-2.d0*x2*(x3*(y1-y2)+x1*(y2-y3))
717  & +x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
718  return
719  end
720 
721  double precision function polderiv3(x1,x2,x3,y1,y2,y3)
722 C-- returns the estimate of the derivative at x3 obtained by a
723 C-- polynomial interpolation using the three points (x_i,y_i).
724  implicit none
725  double precision x1,x2,x3,y1,y2,y3
726  polderiv3=(x3*x3*(-y1+y2)+2.d0*x2*x3*(y1-y3)+x1*x1*(y2-y3)
727  & +x2*x2*(-y1+y3)+2.d0*x1*x3*(-y2+y3))/
728  & ((x1-x2)*(x1-x3)*(x2-x3))
729  return
730  end
731 
732 C----------------------------------------------------------------------
733 
734 C-- G.W. 05/07/2010 Avoid using Fortran 90 intrinsic function
735 C-- "len_trim" since not supported by all Fortran 77 compilers.
736  integer function lentrim(s)
737  implicit none
738  character*(*) s
739  do lentrim = len(s), 1, -1
740  if (s(lentrim:lentrim) .ne. ' ') return
741  end do
742  return
743  end
744 
745 C----------------------------------------------------------------------
getallpdfs
subroutine getallpdfs(prefix, ih, x, q, upv, dnv, usea, dsea, str, sbar, chm, cbar, bot, bbar, glu, phot)
Definition: mstwpdf.f:22
polderiv1
double precision function polderiv1(x1, x2, x3, y1, y2, y3)
Definition: mstwpdf.f:702
polderiv3
double precision function polderiv3(x1, x2, x3, y1, y2, y3)
Definition: mstwpdf.f:722
polderiv2
double precision function polderiv2(x1, x2, x3, y1, y2, y3)
Definition: mstwpdf.f:712
set
set(CMAKE_MACOSX_RPATH 1) set(CMAKE_CACHEFILE_DIR CMakeFiles/) cmake_minimum_required(VERSION 2.8) project(collier) enable_language(Fortran) get_filename_component(Fortran_COMPILER_NAME $
Definition: CMakeLists.txt:65
initialisepdf
subroutine initialisepdf(ip, np, ih, nhess, nx, my, myc0, myb0, xx, yy, ff, cc)
Definition: mstwpdf.f:429
locx
integer function locx(xx, nx, x)
Definition: mstwpdf.f:672
extrapolatepdf
double precision function extrapolatepdf(ip, np, ih, nhess, x, y, nx, my, xx, yy, cc)
Definition: mstwpdf.f:607
getonepdf
double precision function getonepdf(prefix, ih, x, q, f)
Definition: mstwpdf.f:89
interpolatepdf
double precision function interpolatepdf(ip, np, ih, nhess, x, y, nx, my, xx, yy, cc)
Definition: mstwpdf.f:581
getallpdfsalt
subroutine getallpdfsalt(prefix, ih, x, q, xpdf, xphoton)
Definition: mstwpdf.f:62
lentrim
integer function lentrim(s)
Definition: mstwpdf.f:737
if
CBT if(CBT MATCHES debug) message("Build type
Definition: CMakeLists.txt:146