21 & upv,dnv,usea,dsea,str,sbar,chm,cbar,bot,bbar,glu,phot)
24 double precision x,q,upv,dnv,usea,dsea,str,sbar,chm,cbar,
25 & bot,bbar,glu,phot,
getonepdf,up,dn,sv,cv,bv
64 double precision x,q,xpdf(-6:6),xphoton,xvalence,GetOnePDF
68 xpdf(f) = getonepdf(prefix,ih,x,q,f)
69 xvalence = getonepdf(prefix,ih,x,q,f+6)
70 xpdf(-f) = xpdf(f) - xvalence
72 xpdf(0) = getonepdf(prefix,ih,x,q,0)
73 xphoton = getonepdf(prefix,ih,x,q,13)
88 double precision function getonepdf(prefix,ih,x,q,f)
91 parameter(warn=.false.,fatal=.true.)
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,
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)
111 common/mstwcommon/distance,tolerance,
112 & mcharm,mbottom,alphasq0,alphasmz,alphasorder,alphasnfmax
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/
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/
132 if (f.lt.-6.or.f.gt.13)
then
133 print *,
"Error: invalid parton flavour = ",f
137 if (ih.lt.0.or.ih.gt.nhess)
then
138 print *,
"Error: invalid eigenvector number = ",ih
143 if (oldprefix(ih).ne.prefix)
then
151 if (
lentrim(prefix).gt.len(oldprefix(ih)))
then
152 print *,
"Error in GetOnePDF: increase size of oldprefix"
154 else if (
lentrim(prefix)+7.gt.len(filename))
then
155 print *,
"Error in GetOnePDF: increase size of filename"
159 write(
set,
'(I2.2)') ih
161 filename = prefix(1:
lentrim(prefix))//
'.'//
set//
'.dat'
163 print *,
"Reading PDF grid from ",filename(1:
lentrim(filename))
164 open(unit=33,file=filename,iostat=io,status=
'old')
166 print *,
"Error in GetOnePDF: can't open ",
167 & filename(1:
lentrim(filename))
174 read(33,*) dummychar,dummyword,dummyword,dummychar,
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
193 if (mc2.le.qq(1).or.mc2+eps.ge.qq(8))
then
194 print *,
"Error in GetOnePDF: invalid mCharm = ",mcharm
196 else if (mc2.lt.qq(2))
then
200 else if (mc2.lt.qq(3))
then
203 else if (mc2.lt.qq(6))
then
205 else if (mc2.lt.qq(7))
then
213 if (mb2.le.qq(12).or.mb2+eps.ge.qq(17))
then
214 print *,
"Error in GetOnePDF: invalid mBottom = ",mbottom
216 else if (mb2.lt.qq(13))
then
219 else if (mb2.lt.qq(16))
then
233 if (nextraflavours.lt.0.or.nextraflavours.gt.1)
then
234 print *,
"Error in GetOnePDF: invalid nExtraFlavours = ",
242 if (nextraflavours.gt.0)
then
243 if (alphasorder.eq.2)
then
244 read(33,
'(12(1pe12.4))',iostat=io)
245 & (ff(ip,n,m),ip=1,12)
249 read(33,
'(10(1pe12.4))',iostat=io)
250 & (ff(ip,n,m),ip=1,9),ff(12,n,m)
253 if (alphasorder.eq.2)
then
255 read(33,
'(11(1pe12.4))',iostat=io)
256 & (ff(ip,n,m),ip=1,11)
261 read(33,
'(9(1pe12.4))',iostat=io)
262 & (ff(ip,n,m),ip=1,9)
266 print *,
"Error in GetOnePDF reading ",filename
273 read(33,*,iostat=io) dummy
275 print *,
"Error in GetOnePDF: not at end of ",filename
300 oldprefix(ih) = prefix
310 if (qsq.gt.qq(nqc0).and.qsq.lt.qq(nqc0+1)) qsq = qq(nqc0+1)
312 if (qsq.gt.qq(nqb0).and.qsq.lt.qq(nqb0+1)) qsq = qq(nqb0+1)
321 else if (f.ge.1.and.f.le.5)
then
323 else if (f.le.-1.and.f.ge.-5)
then
325 else if (f.ge.7.and.f.le.11)
then
327 else if (f.eq.13)
then
329 else if (abs(f).ne.6.and.f.ne.12)
then
330 if (warn.or.fatal) print *,
"Error in GetOnePDF: f = ",f
334 if (x.le.0.d0.or.x.gt.xmax.or.q.le.0.d0)
then
336 if (warn.or.fatal) print *,
"Error in GetOnePDF: x,qsq = ",
340 else if (abs(f).eq.6.or.f.eq.12)
then
344 else if (qsq.lt.qsqmin)
then
347 print *,
"Warning in GetOnePDF, extrapolating: f = ",f,
348 &
", x = ",x,
", q = ",q
354 & log10(qsqmin),nx,nq,xxl,qql,cc)
356 & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
357 if (f.le.-1.and.f.ge.-5)
then
359 & log10(qsqmin),nx,nq,xxl,qql,cc)
361 & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
367 & log10(qsqmin),nx,nq,xxl,qql,cc)
369 & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
370 if (f.le.-1.and.f.ge.-5)
then
372 & log10(qsqmin),nx,nq,xxl,qql,cc)
374 & log10(1.01d0*qsqmin),nx,nq,xxl,qql,cc)
386 if (abs(res).ge.1.d-5)
then
387 anom = max( -2.5d0, (res1-res)/res/0.01d0 )
391 res = res*(qsq/qsqmin)**(anom*qsq/qsqmin+1.d0-qsq/qsqmin)
393 else if (x.lt.xmin.or.qsq.gt.qsqmax)
then
396 print *,
"Warning in GetOnePDF, extrapolating: f = ",f,
397 &
", x = ",x,
", q = ",q
401 & qsqlog,nx,nq,xxl,qql,cc)
403 if (f.le.-1.and.f.ge.-5)
then
405 & qsqlog,nx,nq,xxl,qql,cc)
411 & qsqlog,nx,nq,xxl,qql,cc)
413 if (f.le.-1.and.f.ge.-5)
then
415 & qsqlog,nx,nq,xxl,qql,cc)
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),
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/
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))
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))
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))
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))
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))
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))
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))
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))
521 ff12(n,m)=0.5*(ff12(n,m)+ff21(n,m))
533 yy0(3)=ff(ip,n+1,m+1)
548 yy12(3)=ff12(n+1,m+1)
561 xxd=xxd+iwt(k,l)*z(k)
569 cc(ip,ih,n,m,k,j)=cl(l)
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),
589 t=(x-xx(n))/(xx(n+1)-xx(n))
590 u=(y-yy(m))/(yy(m+1)-yy(m))
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)
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),
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)))
622 z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
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))*
632 z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
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,
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))*
643 z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
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,
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))*
652 z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
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)))
657 z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
660 print *,
"Error in ExtrapolatePDF"
671 integer function locx(xx,nx,x)
676 double precision x,xx(nx)
687 1
if((ju-jl).le.1)
go to 2
701 double precision function polderiv1(x1,x2,x3,y1,y2,y3)
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))
711 double precision function polderiv2(x1,x2,x3,y1,y2,y3)
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))
721 double precision function polderiv3(x1,x2,x3,y1,y2,y3)
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))