JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
Functions/Subroutines
mstwpdf.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine getallpdfs (prefix, ih, x, q, upv, dnv, usea, dsea, str, sbar, chm, cbar, bot, bbar, glu, phot)
 
subroutine getallpdfsalt (prefix, ih, x, q, xpdf, xphoton)
 
double precision function getonepdf (prefix, ih, x, q, f)
 
subroutine initialisepdf (ip, np, ih, nhess, nx, my, myc0, myb0, xx, yy, ff, cc)
 
double precision function interpolatepdf (ip, np, ih, nhess, x, y, nx, my, xx, yy, cc)
 
double precision function extrapolatepdf (ip, np, ih, nhess, x, y, nx, my, xx, yy, cc)
 
integer function locx (xx, nx, x)
 
double precision function polderiv1 (x1, x2, x3, y1, y2, y3)
 
double precision function polderiv2 (x1, x2, x3, y1, y2, y3)
 
double precision function polderiv3 (x1, x2, x3, y1, y2, y3)
 
integer function lentrim (s)
 

Function/Subroutine Documentation

◆ extrapolatepdf()

double precision function extrapolatepdf ( integer  ip,
integer  np,
integer  ih,
integer  nhess,
double precision  x,
double precision  y,
integer  nx,
integer  my,
double precision, dimension(nx)  xx,
double precision, dimension(my)  yy,
double precision, dimension(np,0:nhess,nx,my,4,4)  cc 
)

Definition at line 607 of file mstwpdf.f.

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

◆ getallpdfs()

subroutine getallpdfs ( character*(*)  prefix,
integer  ih,
double precision  x,
double precision  q,
double precision  upv,
double precision  dnv,
double precision  usea,
double precision  dsea,
double precision  str,
double precision  sbar,
double precision  chm,
double precision  cbar,
double precision  bot,
double precision  bbar,
double precision  glu,
double precision  phot 
)

Definition at line 22 of file mstwpdf.f.

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

◆ getallpdfsalt()

subroutine getallpdfsalt ( character*(*)  prefix,
integer  ih,
double precision  x,
double precision  q,
double precision, dimension(-6:6)  xpdf,
double precision  xphoton 
)

Definition at line 62 of file mstwpdf.f.

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

◆ getonepdf()

double precision function getonepdf ( character, dimension(*)  prefix,
integer  ih,
double precision  x,
double precision  q,
integer  f 
)

Definition at line 89 of file mstwpdf.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

◆ initialisepdf()

subroutine initialisepdf ( integer  ip,
integer  np,
integer  ih,
integer  nhess,
integer  nx,
integer  my,
integer  myc0,
integer  myb0,
double precision, dimension(nx)  xx,
double precision, dimension(my)  yy,
double precision, dimension(np,nx,my)  ff,
double precision, dimension(np,0:nhess,nx,my,4,4)  cc 
)

Definition at line 429 of file mstwpdf.f.

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

◆ interpolatepdf()

double precision function interpolatepdf ( integer  ip,
integer  np,
integer  ih,
integer  nhess,
double precision  x,
double precision  y,
integer  nx,
integer  my,
double precision, dimension(nx)  xx,
double precision, dimension(my)  yy,
double precision, dimension(np,0:nhess,nx,my,4,4)  cc 
)

Definition at line 581 of file mstwpdf.f.

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

◆ lentrim()

integer function lentrim ( character*(*)  s)

Definition at line 737 of file mstwpdf.f.

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

◆ locx()

integer function locx ( double precision, dimension(nx)  xx,
integer  nx,
double precision  x 
)

Definition at line 672 of file mstwpdf.f.

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

◆ polderiv1()

double precision function polderiv1 ( double precision  x1,
double precision  x2,
double precision  x3,
double precision  y1,
double precision  y2,
double precision  y3 
)

Definition at line 702 of file mstwpdf.f.

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

◆ polderiv2()

double precision function polderiv2 ( double precision  x1,
double precision  x2,
double precision  x3,
double precision  y1,
double precision  y2,
double precision  y3 
)

Definition at line 712 of file mstwpdf.f.

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

◆ polderiv3()

double precision function polderiv3 ( double precision  x1,
double precision  x2,
double precision  x3,
double precision  y1,
double precision  y2,
double precision  y3 
)

Definition at line 722 of file mstwpdf.f.

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
initialisepdf
subroutine initialisepdf(ip, np, ih, nhess, nx, my, myc0, myb0, xx, yy, ff, cc)
Definition: mstwpdf.f:429
extrapolatepdf
double precision function extrapolatepdf(ip, np, ih, nhess, x, y, nx, my, xx, yy, cc)
Definition: mstwpdf.f:607