MSTW PDF example code

Clear all definitions.

Clear["Global`*"]

Load Mathematica package containing global variables and functions to access grids.
Change the pdfpath to the location where the mstwpdf.m file is stored.

pdfpath="/home/watt/MSTW/Public/";
SetDirectory[pdfpath];

<<mstwpdf.m

?ReadPDFGrid

?xf

Use of the central PDF set

The prefix should be set to the location of the grid files.

prefix=pdfpath<>"Grids/mstw2008nnlo";
Timing[ReadPDFGrid[prefix,0]]

Write out the central PDFs for a fixed momentum scale x and scale q.

ih=0;x=1.*^-3;q=10.;
upv=xf[ih,x,q,8]
dnv=xf[ih,x,q,7]
usea=xf[ih,x,q,-2]
dsea=xf[ih,x,q,-1]
str=xf[ih,x,q,3]
sbar=xf[ih,x,q,-3]
chm=xf[ih,x,q,4]
cbar=xf[ih,x,q,-4]
bot=xf[ih,x,q,5]
bbar=xf[ih,x,q,-5]
glu=xf[ih,x,q,0]
phot=xf[ih,x,q,13]

Clear the above definitions.

Clear[ih,x,q,upv,dnv,usea,dsea,str,sbar,chm,cbar,bot,bbar,glu,phot];

Check that the number and momentum sum rules are satisfied to reasonable accuracy.  (Note that cv = c - cbar and bv = b - bbar are only non-zero at NNLO.)

NIntegrate[xf[0,x,10.,7]/x,{x,1.*^-9,1.},MaxRecursion->100](* dv integrated over x: should be 1 *)

NIntegrate[xf[0,x,10.,8]/x,{x,1.*^-9,1.},MaxRecursion->100](* uv integrated over x: should be 2 *)

NIntegrate[xf[0,x,10.,9]/x,{x,1.*^-9,1.},MaxRecursion->100](* sv integrated over x: should be zero *)

NIntegrate[xf[0,x,10.,10]/x,{x,1.*^-9,1.},MaxRecursion->100](* cv integrated over x: should be zero *)

NIntegrate[xf[0,x,10.,11]/x,{x,1.*^-9,1.},MaxRecursion->100](* bv integrated over x: should be zero *)

NIntegrate[xf[0,x,10.,9],{x,1.*^-6,1.},MaxRecursion->100](* xsv integrated over x: should be non-zero *)

NIntegrate[xf[0,x,10.,10],{x,1.*^-6,1.},MaxRecursion->100](* xcv integrated over x: should be non-zero only at NNLO *)

NIntegrate[xf[0,x,10.,11],{x,1.*^-6,1.},MaxRecursion->100](* xbv integrated over x: should be non-zero only at NNLO *)

NIntegrate[Sum[xf[0,x,q,f]/.{q->10.},{f,-6,6}],{x,1.*^-6,1.},MaxRecursion->100] (* momentum sum rule: should be 1 *)

Plot the valence quark distributions of the central set.
(Note that cv = c - cbar and bv = b - bbar are only non-zero at NNLO.)

Plot[xf[0,10.^ix,q,7]/.{q->10.},{ix,-6,0},AxesLabel->{"Log[10,x]","xdv = xd - xdbar"}]

Plot[xf[0,10.^ix,q,8]/.{q->10.},{ix,-6,0},PlotRange->All,AxesLabel->{"Log[10,x]","xuv = xu - xubar"}]

Plot[xf[0,10.^ix,q,9]/.{q->10.},{ix,-6,0},AxesLabel->{"Log[10,x]","xsv = xs - xsbar"}]

Plot[xf[0,10.^ix,q,10]/.{q->10.},{ix,-6,0},AxesLabel->{"Log[10,x]","xcv = xc - xcbar"}]

Plot[xf[0,10.^ix,q,11]/.{q->10.},{ix,-6,0},AxesLabel->{"Log[10,x]","xbv = xb - xbbar"}]

Use of the eigenvector PDF sets

?distance

?tolerance

Get the central value of the gluon distribution.

ih=0;x=1.*^-3;q=10.;
glu=xf[ih,x,q,0]
distance
tolerance

Now get the value of the gluon distribution at a distance t from the global minimum along the first eigenvector direction in the '+' direction.  The distance t is chosen to give a tolerance T ensuring all data sets are described within their 90% confidence level (C.L.) limits.

ReadPDFGrid[prefix<>".90cl",1]
glu=xf[1,x,q,0]
distance
tolerance

Now the same thing, but with the tolerance chosen to ensure that all data sets are described only within their one-sigma (68%) C.L. limits.

ReadPDFGrid[prefix<>".68cl",1]
glu=xf[1,x,q,0]
distance
tolerance

Clear[ih,x,q,glu];

?nhess
nhess

Load the eigenvector sets into memory.  Choose either the 90% C.L. or the 68% C.L. sets.  Warning : the next line will take a long time to evaluate!  (For an estimate look at the time taken to read the central fit above and multiply by nhess.)

Timing[Do[ReadPDFGrid[prefix<>".90cl",ih],{ih,1,nhess}]]

Define functions to give the uncertainty on the parton distributions using both the formula for asymmetric errors [eqs.(51,52) of arXiv:0901.0002] and the formula for symmetric errors [eq.(50) of same paper].

delp[x_Real,q_Real,f_Integer]:=Sqrt[Sum[Max[xf[2i-1,x,q,f]-xf[0,x,q,f],xf[2i,x,q,f]-xf[0,x,q,f],0.]^2,{i,1,nhess/2}]]
delm[x_Real,q_Real,f_Integer]:=Sqrt[Sum[Max[xf[0,x,q,f]-xf[2i-1,x,q,f],xf[0,x,q,f]-xf[2i,x,q,f],0.]^2,{i,1,nhess/2}]]
del[x_Real,q_Real,f_Integer]:=0.5*Sqrt[Sum[(xf[2i-1,x,q,f]-xf[2i,x,q,f])^2,{i,1,nhess/2}]]

Plot gluon distribution with uncertainty band at 10 GeV.
(The Plot option Filling requires Mathematica 6.)

Plot[Evaluate[{xf[0,x,q,f],xf[0,x,q,f]+delp[x,q,f],xf[0,x,q,f]-delm[x,q,f]}/.{q->10.,x->10.^ix,f->0}],{ix,-6,0},AxesLabel->{"Log[10,x]","xg at 10 GeV"},Filling->{2->{3}}]

Plot fractional uncertainty in gluon distribution at 10 GeV.
(The Plot option "Filling" requires Mathematica 6.)

Plot[Evaluate[{0.,delp[x,q,f]/xf[0,x,q,f],-delm[x,q,f]/xf[0,x,q,f]}/.{q->10.,x->10.^ix,f->0}],{ix,-6,0},AxesLabel->{"Log[10,x]","Fractional uncertainty in xg at 10 GeV"},Filling->{2->{3}}]

Define a function to convert a real number to equivalent of "1PE12.4" Fortran format or"%12.4E" C format.

format[number_Real]:=ToString[ScientificForm[number,10,NumberFormat->(ToString[PaddedForm[ToExpression[#1],{6,4}]]<>"E"<>ToString[PaddedForm[If[Head[ToExpression[#3]]==Integer,ToExpression[#3],0,0],2,NumberPadding->{"0",""},SignPadding->True,NumberSigns->{"-","+"}]]&)]]

Define a function to pad a string to a certain length.

PadString[str_String,len_Integer]:=Module[{temp=str},While[StringLength[temp]<len,temp=" "<>temp]; Return[temp];]

Define a function to write the PDFs as a function of x to a file at a fixed value of q.
Extrapolation will be used for x < 10^-6.

WritePDFvsX[q_Real,f_Integer]:=Module[{flavours,ix,nx=100,xmin=1.*^-7,xmax=0.99,x,line},
flavours={"bbar","cbar","sbar","ubar","dbar","glu", "dn","up","str","chm","bot"};
filename="x"<>flavours[[f+6]]<>"_vs_x_math.dat";
file=OpenWrite[filename];
WriteString[file," # q = "<>format[q]<>"\n"];
WriteString[file," #"<>PadString["x",10]<>PadString["x"<>flavours[[f+6]],12]<>PadString["error+",12]<>PadString["error-",12]<>PadString["error",12]<>"\n"];
Do[{x=10.^(Log[10,xmin]+(ix-1.)/(nx-1.)*(Log[10,xmax]-Log[10,xmin]));
line=format[x]<>format[xf[0,x,q,f]]<>format[delp[x,q,f]]<>format[delm[x,q,f]]<>format[del[x,q,f]]<>"\n";
WriteString[file,line];
},{ix,1,nx}];
Close[file];
Print["x"<>flavours[[f+6]]," vs. x for q2 = ",q^2," written to ",filename];
]

Do[WritePDFvsX[q,f]/.{q->10.},{f,-5,5}]

Define a function to write the PDFs as a function of q^2 to a file at a fixed value of x.
Extrapolation will be used for q^2 < 1 GeV^2 and for q^2 > 10^9 GeV^2.

WritePDFvsQ2[x_Real,f_Integer]:=Module[{flavours,iq,nq=100,q2min=0.5,q2max=1.*^10,q,line},
flavours={"bbar","cbar","sbar","ubar","dbar","glu", "dn","up","str","chm","bot"};
filename="x"<>flavours[[f+6]]<>"_vs_q2_math.dat";
file=OpenWrite[filename];
WriteString[file," # x = "<>format[x]<>"\n"];
WriteString[file," #"<>PadString["q2",10]<>PadString["x"<>flavours[[f+6]],12]<>PadString["error+",12]<>PadString["error-",12]<>PadString["error",12]<>"\n"];
Do[{q=Sqrt[10.^(Log[10,q2min]+(iq-1.)/(nq-1.)*(Log[10,q2max]-Log[10,q2min]))];
line=format[q^2]<>format[xf[0,x,q,f]]<>format[delp[x,q,f]]<>format[delm[x,q,f]]<>format[del[x,q,f]]<>"\n";
WriteString[file,line];
},{iq,1,nq}];
Close[file];
Print["x"<>flavours[[f+6]]," vs. q2 for x = ",x," written to ",filename];
]

Do[WritePDFvsQ2[x,f]/.{x->1.*^-3},{f,-5,5}]

Mathlink interface to Fortran code for (Subscript[α, S])

Values of mCharm, mBottom, alphaSQ0, alphaSMZ, alphaSorder,alphaSnfmax are filled when the ReadPDFGrid function is called.

?mCharm
mCharm

?mBottom
mBottom

?alphaSQ0
alphaSQ0

?alphaSMZ
alphaSMZ

?alphaSorder
alphaSorder

?alphaSnfmax
alphaSnfmax

Install Mathematica definitions to call Fortran code for alphaS via Mathlink interface.
Requires the binary alphaS.exe to be compiled from alphaS.c, alphaS.tm and alphaS.f.  Under Linux, do this with "mcc -o alphaS.exe alphaS.c alphaS.tm alphaS.f -lg2c".

Install["alphaS.exe"]

?InitAlphaS

InitAlphaS[alphaSorder,1.,1.,alphaSQ0,mCharm,mBottom,1.*^10]

?AlphaS

AlphaS[1.]

AlphaS[91.1876]

Alternatively, initialise alphaS at the Z mass.

InitAlphaS[alphaSorder,1.,91.1876,alphaSMZ,mCharm,mBottom,1.*^10]

AlphaS[91.1876]

AlphaS[1.]

Plot[AlphaS[mur],{mur,1,100},AxesLabel->{"mur (GeV)","alphaS(mur)"}]