Package for calculating the evolved Sivers functions based on a range of fits and models (Updated Oct 18, 2011).
	Bochum fits are from: (Phys.Rev.D73:014021(2006)).
	Torino fits are from:  (Eur.Phys J. A 39:89 (2009)).
	Soft evolution data are from:  (Phys.Lett.B 562, 182 (2003))

BeginPackage["Sivers`"]

BochumUpSivers::usage ="BochumupSivers[xx,kk,qq] is the Sivers function for the up quark at Bjorken-x xx, transverse momentum kk (in GeV), and hard scale qq (in GeV) evolved from the Bochum fits."
BochumDownSivers::usage ="BochumDownSivers[xx,kk,qq] is the Sivers function for the down quark at Bjorken-x xx, transverse momentum kk (in GeV), and hard scale qq (in GeV) evolved from the Bochum fits."
TorinoUpSivers::usage ="TorinoUpSivers[xx,kk,qq] is the Sivers function for the up quark at Bjorken-x xx, transverse momentum kk (in GeV), and hard scale qq (in GeV) evolved from the Torino fits."
TorinoDownSivers::usage ="TorinoDownSivers[xx,kk,qq] is the Sivers function for the down quark at Bjorken-x xx, transverse momentum kk (in GeV), and hard scale qq (in GeV) evolved from the Torino fits."

Need MSTW PDFs for collinear factors in Torino fits.

pdfpath="/local12/rogers/mstw";
SetDirectory[pdfpath];
<<mstwpdf.m
prefix=pdfpath<>"/Grids/mstw2008lo";
Timing[ReadPDFGrid[prefix,0]]

Basic pQCD parameters and integrated collinear PDFs.

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];
uppdf[x_,q_] := (1 / x) xf[0,x,q,2];
pdfqup[x_,bt_] := (1 / x) uppdf[x,mub[bt]];
zerothorderup[xa_] := xa pdfqup[xa,1] ;
downpdf[x_,q_] := (1 / x) xf[0,x,q,1];
pdfqdown[x_,bt_] := (1 / x) downpdf[x,mub[bt]];
zerothorderdown[xa_] := xa pdfqdown[xa,1] ;

Basic  input for Torino fits.

Nu = .35;
alphau = .73;
beta = 3.46;
Nup[xa_] := Nu xa^alphau (1 - xa)^beta (alphau + beta)^(alphau + beta) / alphau^alphau / beta^beta
fup[xa_] := zerothorderup[xa]
Nd = .90; (** Note this is actually negative.  Remember to include minus sign in the end.**)
alphad = 1.08;
Ndown[xa_] := Nd xa^alphad (1 - xa)^beta (alphad + beta)^(alphad + beta) / alphad^alphad / beta^beta
fdown[xa_] := zerothorderdown[xa]
M1 = Sqrt[.34];
kt2 = .25;
Mp = .938;

Begin["`Private`"]

xval = .1;
g2val = .68;

kt2zerovalbochum =.2 
kt2zerovaltorino = M1^2 kt2 / (M1^2 + kt2)

fonepertbochumup[xa_,kt2zero_,g2_] := -(2 Mp^2 / kt2zero)(.17 xa^(.66 - 1) (1 - xa)^5)
fonepertturinup[xa_,kt2zero_,g2_] :=  -(4 Pi Mp Sqrt[2 Exp[1]]  / M1 / Pi / kt2) Nup[xa] fup[xa] (kt2zero / 2 )^2 / kt2zero

fonepertbochumdown[xa_,kt2zero_,g2_] := -(2 Mp^2 / kt2zero)(.17 xa^(.66 - 1) (1 - xa)^5)
fonepertturindown[xa_,kt2zero_,g2_] :=  -(4 Pi Mp Sqrt[2 Exp[1]]  / M1 / Pi / kt2) Ndown[xa] fdown[xa] (kt2zero / 2 )^2 / kt2zero

Fits of normalization and slope as a function of hard scale q:

gaussnormfitlowbochumup[q_]:=1.76248/(-28.0094 q+20.0822 q^2-0.374729 q^4+0.00296727 q^6)
gaussnormfithighbochumup[q_]:= 0.000141945 + 1145.69 / (-44440.5 Sqrt[q] + 43335.0 q)
gaussparfitbochumup[q_] := 1.45226/Log[0.557029 q^2]

gaussnormfitlowtorinoup[q_] := 0.856642/(-7.92053 q+5.3952 q^2- 0.0307745 q^4-0.00134312 q^6)
gaussnormfithightorinoup[q_] := 0.00020356 + 1635.84/(-47186.2 Sqrt[q] + 43453.6 q)
gaussparfittorinoup[q_] := 1.4556/Log[0.513906 q^2]

gaussnormfitlowtorinodown[q_] := 4.67378/(-32.0522 q+21.8329 q^2-0.124536 q^4-0.00543525 q^6)
gaussnormfithightorinodown[q_] := 0.000274447 + 746.554 / (-15972.4 Sqrt[q] + 14708.9 q) 
gaussparfittorinodown[q_] := 1.4556/Log[0.513906 q^2]

Extrapolations.  Note that we have a (Mp / 2 / kt) to convert to F^{\perp} following Eq.(3) of Ref. (Eur.Phys J. A 39:89 (2009)).
These extrapolations come from SiversCalcBochu.nb, SiversCalcTorinoUp.nb, SiversCalcTorinoDown.nb.

BochumUpSivers[xa_,kt_,q_] :=  (-Mp / 2 / kt / xval)(fonepertbochumup[xa,kt2zerovalbochum,g2val] / fonepertbochumup[xval,kt2zerovalbochum,g2val])Piecewise[{{gaussnormfitlowbochumup[q] kt Exp[-gaussparfitbochumup[q] kt^2],q≤ 5},{gaussnormfithighbochumup[q] kt Exp[-gaussparfitbochumup[q] kt^2],q>5}}];
TorinoUpSivers[xa_,kt_,q_] := (-Mp / 2 / kt / xval) (fonepertturinup[xa,kt2zerovaltorino,g2val] / fonepertturinup[xval,kt2zerovaltorino,g2val])Piecewise[{{gaussnormfitlowtorinoup[q] kt Exp[-gaussparfittorinoup[q] kt^2],q≤ 5},{gaussnormfithightorinoup[q] kt Exp[-gaussparfittorinoup[q] kt^2],q>5}}];

(*** Note that the Bochum down Sivers is just the negative of the up Sivers.***)
BochumDownSivers[xa_,kt_,q_] :=  (Mp / 2 / kt / xval) (fonepertbochumup[xa,kt2zerovalbochum,g2val] / fonepertbochumup[xval,kt2zerovalbochum,g2val])Piecewise[{{gaussnormfitlowbochumup[q] kt Exp[-gaussparfitbochumup[q] kt^2],q≤ 5},{gaussnormfithighbochumup[q] kt Exp[-gaussparfitbochumup[q] kt^2],q>5}}];
(*** Note we change the sign for the down functions ***)
TorinoDownSivers[xa_,kt_,q_] :=  (Mp / 2 / kt / xval) (fonepertturindown[xa,kt2zerovaltorino,g2val] / fonepertturindown[xval,kt2zerovaltorino,g2val])Piecewise[{{gaussnormfitlowtorinodown[q] kt Exp[-gaussparfittorinodown[q] kt^2],q≤ 5},{gaussnormfithightorinodown[q] kt Exp[-gaussparfittorinodown[q] kt^2],q>5}}];

End[] 
EndPackage[]