Up quark Sivers function from small to large Q  using Bochum Gaussian fits (Oct 18, 2011).

Clear["Global *"]

Setting up the MSTW PDF Path

pdfpath="/local12/rogers/mstw";
SetDirectory[pdfpath]; 
<<mstwpdf.m
prefix=pdfpath<>"/Grids/mstw2008lo";
Timing[ReadPDFGrid[prefix,0]]
Install["alphaS.exe"]
InitAlphaS[1,1.,1.,alphaSQ0,mCharm,mBottom,1.*^10]
?AlphaS

Setting up constant parameters, color factors, etc. needed in the calculation of perturbative parts.

CF = 4 /3 ;
TF = 1 / 2;
C1 = 1.123;   
bmax = 0.5;
bstar[bt_] := bt / Sqrt[1 + bt^2 / bmax^2]; 
mub[bt_] := C1 / bstar[bt];
Mp = .938;

Setting up CS kernel, anomalous dimensions, etc.

K[mu_,bt_] := -(AlphaS[mu] / Pi) CF (Log[mu^2 bt^2] - Log[4] + 2 EulerGamma);
Clear[gammaK]; Clear[Dgam];
gammaK[mu_] := 2 AlphaS[mu] CF / Pi;
Dgam[mu_,xi_] := (AlphaS[mu] CF / Pi) ( 3 / 2 - Log[xi/mu^2]);

Combine all of the parts of Eq... in Ref... that can be calculated purely within perturbation theory.

PertEvol[mu_,bt_,xi_,mu0_,Q0_] := N[Exp[Log[Sqrt[xi] / Q0] K[mub[bt],bstar[bt]]
+ NIntegrate[Dgam[muprime,muprime^2] / muprime - Log[Sqrt[xi]/muprime] gammaK[muprime] / muprime,{muprime,mu0,mu},AccuracyGoal-> 5] + NIntegrate[Log[Sqrt[xi]/Q0] gammaK[muprime] / muprime,{muprime,mu0,mub[bt]},AccuracyGoal-> 5] ],10]

Give the basic form of the parametrization of the non - perturbative part.  Assume a basic Gaussian form.  The cpar1 and cpar2 are the parameters of the scale dependent and non scale dependent parts.

NonPertEvol[g2_,bt_,xi_,Q0_] := Exp[-g2 Log[Sqrt[xi]/Q0] bt^2 / 2]

The Bochum fit from (Phys.Ref.D73:014021(2006)).

fonepert[xa_,kt2zero_] := -(2 Mp^2 / kt2zero)(.17 xa^(.66 - 1) (1 - xa)^5)
InitFunc[xa_,bt_,kt2zero_] := - (kt2zero fonepert[xa,kt2zero]/ 2) bt Exp[-kt2zero bt^2 / 4];
Evolved[xa_,bt_,kt2zero_,g2_,Q0_,mu_,xi_] :=InitFunc[xa,bt,kt2zero] PertEvol[mu,bt,xi,Q0,Q0]NonPertEvol[g2,bt,xi,Q0]

Choose the value of x for all subsequent calculations:
Input the parameters relavant to the evolution (see Phys.Rev.D83:114042(2011)).

xval = .1;
g2val = .68;
Q0val = Sqrt[2.4];
kt2zeroval =.2 ;
Q = Sqrt[2.4];

Make a table to set up an interpolation function for later manipulations.

Clear[Evolvedrt2p4];
Evolvedrt2p4 = {}
Do[Evolvedrt2p4 = Append[Evolvedrt2p4,{.01 i,Evolved[.1,.01 i,kt2zeroval,g2val,Q0val,Q,Q^2] }],{i,1400}]

Clear[Evolvedrt2p4forplot];
Evolvedrt2p4forplot = Interpolation[Evolvedrt2p4]

Momentum space Sivers effect evolved to HERMES scale sqrt{2.4}.  
Note that the Sivers function is represented in the xΔF convention of Eq.(3) and Fig.6 in Ref. (Eur.Phys J. A 39 89 (2009))!

Momspacert2p4[kt_] := -(2 xval kt / Mp) (-1 / 2 / Pi/ kt) NIntegrate[bt BesselJ[1,bt kt] Evolvedrt2p4forplot[bt],{bt,.01,14}];

plot2p4 = LogPlot[Momspacert2p4[kt],{kt,.01,1},PlotRange-> {.000001,2}]

Setup the fit to a Gaussian.  form1 is the Gaussian parametrization.

form1 = a Exp[-b kt^2]

Plots for paper : Repeat the above steps for different values of Q and output to data files.
Note that the Sivers function HERE is expressed in the F^\perp_1T convention of Eq.() in Ref...!

Q=91.19;
Clear[Evolvedq];
Evolvedq={};
Do[Evolvedq=Append[Evolvedq,{.01 i,Evolved[xval,.01 i,kt2zeroval,g2val,Q0val,Q,Q^2]}],{i,1400}];
Clear[Evolvedqforplot];
Evolvedqforplot=Interpolation[Evolvedq];
Momspaceq91p19[kt_]:=(-1/2/Pi/kt) NIntegrate[bt BesselJ[1,bt kt] Evolvedqforplot[bt],{bt,.01,14}];

plotq91p19=LogPlot[-Momspaceq91p19[kt],{kt,.01,10},PlotRange->{.000001,2}];
Clear[Fitq];
Fitq={};
Fitq=Flatten[Table[{{.01 j,Momspaceq91p19[.01 j]}},{j,1,300}],1];
Clear[gaussfitq]
gaussfitq=FindFit[Fitq,form1,{a,b},kt]
gaussfitqplot91p19[kt_]:=form1/.a->gaussfitq[[1,2]]/.b->gaussfitq[[2,2]];
fitted91p19=LogPlot[-gaussfitqplot91p19[kt],{kt,.01,10}];
Show[plotq91p19,fitted91p19]

Plot[gaussfitqplot91p19[kt]/Momspaceq91p19[kt],{kt,0,10}]

Q=5.0;
Clear[Evolvedq];
Evolvedq={};
Do[Evolvedq=Append[Evolvedq,{.01 i,Evolved[xval,.01 i,kt2zeroval,g2val,Q0val,Q,Q^2]}],{i,1400}];
Clear[Evolvedqforplot];
Evolvedqforplot=Interpolation[Evolvedq];
Momspaceq5[kt_]:=(-1/2/Pi/kt) NIntegrate[bt BesselJ[1,bt kt] Evolvedqforplot[bt],{bt,.01,14}];

plotq5=LogPlot[-Momspaceq5[kt],{kt,.01,10},PlotRange->{.000001,2}];
Clear[Fitq];
Fitq={};
Fitq=Flatten[Table[{{.01 j,Momspaceq5[.01 j]}},{j,1,300}],1];
Clear[gaussfitq];
gaussfitq=FindFit[Fitq,form1,{a,b},kt]
gaussfitqplot5[kt_]:=form1/.a->gaussfitq[[1,2]]/.b->gaussfitq[[2,2]];
fitted5=LogPlot[-gaussfitqplot5[kt],{kt,.01,10}];
Show[plotq5,fitted5]
Plot[gaussfitqplot5[kt]/Momspaceq5[kt],{kt,0,5}]

Q=Sqrt[2.4];
Clear[Evolvedq];
Evolvedq={};
Do[Evolvedq=Append[Evolvedq,{.01 i,Evolved[xval,.01 i,kt2zeroval,g2val,Q0val,Q,Q^2]}],{i,1400}];
Clear[Evolvedqforplot];
Evolvedqforplot=Interpolation[Evolvedq];
Momspaceqrt2p4[kt_]:=(-1/2/Pi/kt) NIntegrate[bt BesselJ[1,bt kt] Evolvedqforplot[bt],{bt,.01,14}];

plotqrt2p4=LogPlot[-Momspaceqrt2p4[kt],{kt,.01,1},PlotRange->{.0001,60}];
Clear[Fitq];
Fitq={};
Fitq=Flatten[Table[{{.01 j,Momspaceqrt2p4[.01 j]}},{j,1,70}],1];
Clear[gaussfitq];
gaussfitq=FindFit[Fitq,form1,{a,b},kt]
gaussfitqplotrt2p4[kt_]:=form1/.a->gaussfitq[[1,2]]/.b->gaussfitq[[2,2]];
fittedrt2p4=LogPlot[-gaussfitqplotrt2p4[kt],{kt,.01,1}];
Show[plotqrt2p4,fittedrt2p4]
Plot[gaussfitqplotrt2p4[kt]/Momspaceqrt2p4[kt],{kt,0,1}]