Quark TMD FF calculation at Q = (Sqrt[2.4]), 5 and 91.19 GeV
with bmax=0.5 GeV^-1.

Clear["Global *"]

Setting up the collinear parts.

collpath="/Users/mert/Work/ted/TMD";
SetDirectory[collpath];
<<tmd.m;
?FFCollinear

Calling MRST2008 NLO PDF Set. We need this for AlphaS
If you do not have MRST2008 PDF sets installed skip this and use our simple implementation of 
AlphaS below.

pdfpath="/Users/mert/Work/C_files/MRST_2008";
SetDirectory[pdfpath];
<<mstwpdf.m
prefix=pdfpath<>"/Grids/mstw2008nlo";
Timing[ReadPDFGrid[prefix,0]]
Install["alphaS.exe"]
InitAlphaS[1,1.,1.,alphaSQ0,mCharm,mBottom,1.*^10]
?AlphaS

If you do not have MRST2008 PDF sets installed uncomment and use the following implementation of 
AlphaS. However, note that the results below and in arXiv:1101.5057 are obtained using MSTW AlphaS.

(*
AlphaS[mu_]:=Module[{LamQCD,nf,beta0,beta1},
nf=5;
beta0=11.-nf 2./3.;
beta1=102.-nf 38./3.;
LamQCD=0.226;

4 Pi/beta0/Log[mu^2/LamQCD^2] (1-beta1/beta0^2 Log[Log[mu^2/LamQCD^2]]/Log[mu^2/LamQCD^2])
]
*)

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];

Setting up CS evolution 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.31 in arXiv:1101.5057 that can be calculated purely within perturbation theory.

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

Combine this with the part calculated using integrated FFs, so everything that can be calculated without extra non - perturbative input is calcuted.

CompletePert[za_,bt_,mu_,xi_] := PertEvol[mu,bt,xi] FFCollinear[up,za,bt];

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.

NonPert[xa_,bt_,xi_,M_,cpar1_,cpar2_] := Exp[-cpar1 bt^2 - cpar2 bt^2 Log[Sqrt[xi]/M ]];

Choose the value of z :

zval = .6;

Setting up the parameters for the Brock - Landry - Nadolsky - Yuan fit (BLNY), modified to fit the Schweitzer, Teckentrup, Metz (STM) Gaussian at largre z and small Q.

Q0=3.2;
(*Note that Q0 here is 2*Q_0 of Eq.32*)
g2=0.68;
g1=0.08-g2 Log[Sqrt[2.4]/Q0] zval^2;(*Matching the inverse fourier transform of Eq.39 with Eq. 33 at Q=Sqrt[2.4] GeV*)

C2 = g2/2;
C3 = g1/2 / zval^2;

Set up the full form of Eq.31 from arXiv:1101.5057 in coordinate b - space.

NLOcoord[xa_,bt_,xi_,mu_,M_,cpar1_,cpar2_] := NonPert[xa,bt,xi,M,cpar1,cpar2] CompletePert[xa,bt,mu,xi];

For doing the numerical integration, we need to make tables for NLOcoord :

Clear[NLOmuLow];
NLOmuLow = {};
Do[NLOmuLow = Append[NLOmuLow,{.01 i,NLOcoord[zval,.01 i,2.4,Sqrt[2.4],Q0,C3,C2]}],{i,1000}];
Clear[NLOmuMed];
NLOmuMed = {};
Do[NLOmuMed = Append[NLOmuMed,{.01 i,NLOcoord[zval,.01 i,5^2,5,Q0,C3,C2]}],{i,1000}];
Clear[NLOmuHigh];
NLOmuHigh = {};
Do[NLOmuHigh = Append[NLOmuHigh,{.01 i,NLOcoord[zval,.01 i,91.19^2,91.19,Q0,C3,C2]}],{i,1000}];
Clear[NLOmuLowplot];
NLOmuLowplot := Interpolation[NLOmuLow]
Clear[NLOmuMedplot];
NLOmuMedplot := Interpolation[NLOmuMed]
Clear[NLOmuHighplot];
NLOmuHighplot := Interpolation[NLOmuHigh]

Calculate the momentum space TMD FFs for Q = (Sqrt[2.4]) , 5, and 91.19 GeV:

Clear[TMDNLOmuLow];
TMDNLOmuLow[kt_] := (1 / 2 / Pi) NIntegrate[bt NLOmuLowplot[bt]BesselJ[0,bt kt],{bt,0.01,10},AccuracyGoal-> 5]
Clear[TMDNLOmuMed];
TMDNLOmuMed[kt_] := (1 / 2 / Pi) NIntegrate[bt NLOmuMedplot[bt]BesselJ[0,bt kt],{bt,0.01,10},AccuracyGoal-> 5]
Clear[TMDNLOmuHigh];
TMDNLOmuHigh[kt_] := (1 / 2 / Pi) NIntegrate[bt NLOmuHighplot[bt]BesselJ[0,bt kt],{bt,0.01,10},AccuracyGoal-> 7]

Set up tables for plotting and exporting:

Clear[HighFit]
HighFit = {};
HighFit = Flatten[Table[{{.01j,TMDNLOmuHigh[.01 j]}},{j,1,2000}],1];
Clear[MedFit]
MedFit = {};
MedFit = Flatten[Table[{{.01j,TMDNLOmuMed[.01 j]}},{j,1,600}],1];
Clear[LowFit]
LowFit = {};
LowFit = Flatten[Table[{{.01j,TMDNLOmuLow[.01 j]}},{j,1,600}],1];

listforplot = ListLogPlot[{LowFit,MedFit,HighFit},PlotRange->{0.0001,10}]

We will now fit TMD FF to a Gaussian function  to demonstrate
quantitatively the effects of evolution.

Setting up functional form to fit

form1 = a Exp[-b kt^2 ];

Making the fit to the TMD FF at Q = 91.19 GeV

gfithigh =FindFit[HighFit,form1,{a,b},kt]
gHigh[kt_]:=form1/.gfithigh

Show[{LogPlot[{gHigh[kt]},{kt,0,20},PlotRange->{0.0001,1.5}], ListLogPlot[{HighFit},PlotRange->{0.0001,1.5}]}]

Calculate  (Sqrt[<kT^2>]) in both cases to see the difference

ktCutoff=11;
avktG=Sqrt[2 Pi NIntegrate[gHigh[kt]*kt^3,{kt,0,ktCutoff}]]
highinterpolate=Interpolation[HighFit];
avktK=Sqrt[2 Pi NIntegrate[highinterpolate[kt]*kt^3,{kt,0,ktCutoff}]]
Abs[(avktG-avktK)]/avktG*100 (*This is the percentage difference *)