Comparison of Perturbative and Extracted b dependence
Used for Aidala, Field, Gamberg, Rogers:  arXiv:1401.2654

Clear["Global *"]
Needs["PlotLegends`"]

Clear[bmax]
C1 = 1.123;
BB1 = .1669;    (** Gaussian slope from COMPASS for the initial Q.**)
BB2 = .2325; (** Gaussian slope from COMPASS for the final Q. **)
AA = (9/(4 Pi)); (** Parameter for alphaS.  Normalization from 1/beta0 with flavors = 3. Corresponds to 1 / A in paper. ***)
LamQCD = .2123; (** Parameter for alphaS.  Lambda_QCD.***)
Cevol = .0234; (** Maximum gaussian evolution parameter extracted from Gaussian model.***)
Q0 = 1.049; (*** Initial Scale of CS evolution. ***)
Qf = 2.114; (** Final Scale of CS evolution.***)
BB1b = .1717; (** Our Gaussian slope for the initial Q.**)
BB2b = .2477;   (** Our Gaussian slope for the final Q.**)
Kap = 1.14; (** Width parameter for our Kaplan fit - see FFinakap[bt_] below**)
nn = 2.593;  (**Exponent for our Kaplan fit - see FFinakap[bt_] below**)

bstar[bt_,bmax_] := bt/ Sqrt[1 + bt^2 / bmax^2]  (** Prescription for matching to b_max. **)
mub[bt_,bmax_] := C1 / bstar[bt,bmax];
AlphaS[mu_]:= 1 / (2 AA Log[mu/LamQCD]);

gK1[bt_,g2_] :=(g2/2)bt^2

gK2[bt_,g2_,aa_] :=(aa^2) (g2/2)Log[1 + bt^2 / aa^2]

Initial and Final Gaussian Fits:

Finit[bt_] := Exp[-((BB1 bt^2)/4)] 
Ffina[bt_] := Exp[-((BB2 bt^2)/4)]
Finit2[bt_] := Exp[-((BB1b bt^2)/4)] 
Ffina2[bt_] := Exp[-((BB2b bt^2)/4)]
FFinakap[bt_] :=((bt)^(nn^2 - 1)) BesselK[1 - nn^2,Kap bt](Gamma[nn^2] (1 / 2 / Kap)(2 / Kap)^(nn^2))^(-1)

Evolution Exponents

evolfactPDF[Q_,bt_,bmax_,g2_] :=(2 /AA/ Pi) (Log[Log[Q/LamQCD]/Log[mub[bt,bmax]/LamQCD]] - (4/3) Log[Q/LamQCD](Log[Log[Q/LamQCD]/Log[mub[bt,bmax]/LamQCD]])+ (4/3)Log[Q/mub[bt,bmax]]) - 2 gK1[bt,g2] Log[Q/Q0] - (2 / AA /Pi) (Log[Log[Q0/LamQCD]/Log[mub[bt,bmax]/LamQCD]] - (4/3) Log[Q0/LamQCD](Log[Log[Q0/LamQCD]/Log[mub[bt,bmax]/LamQCD]])+ (4/3)Log[Q0/mub[bt,bmax]]) 
(** Evolution factor with Gaussian g_K.**)

evolfactPDF2[Q_,bt_,bmax_,g2_,aa_] :=(2 / Pi/ AA) (Log[Log[Q/LamQCD]/Log[mub[bt,bmax]/LamQCD]] - (4/3) Log[Q/LamQCD](Log[Log[Q/LamQCD]/Log[mub[bt,bmax]/LamQCD]])+ (4/3)Log[Q/mub[bt,bmax]]) - 2 gK2[bt,g2,aa] Log[Q/Q0] - (2 / Pi / AA) (Log[Log[Q0/LamQCD]/Log[mub[bt,bmax]/LamQCD]] - (4/3) Log[Q0/LamQCD](Log[Log[Q0/LamQCD]/Log[mub[bt,bmax]/LamQCD]])+ (4/3)Log[Q0/mub[bt,bmax]]) 
(** Evolution factor version with new g_K.**)

Kapp[bt_] := Cevol bt^2  (** function for evolving in the Gaussian approximation.**)
evolfactapprox[Q_,bt_] := Exp[-Kapp[bt] Log[Q/ Q0]] (**Evolution Factor in the Gaussian approximation.**)

Multiplying Gaussian by the exponential of evolution terms.

fullPDF1[bt_,Q_,bmax_,g2_] := Finit[bt] Exp[evolfactPDF[Q,bt,bmax,g2]]  (**Evolved cross section.**)
fullPDF2[bt_,Q_,bmax_,g2_,aa_] := Finit[bt] Exp[evolfactPDF2[Q,bt,bmax,g2,aa]]  (**Evolved cross section.**)
fullPDFapprox[bt_,Q_] := Finit[bt] evolfactapprox[Q,bt]                        (** Using Gaussian approximation.**)

Checking the evolution parameter.

Plot[{bt fullPDFapprox[bt,Q0],bt Finit[bt]},{bt,.01,12}]
Plot[{bt fullPDFapprox[bt,Qf],bt Ffina[bt]},{bt,.01,12}]

Plot 3: Plotting the normalized initial and final Gaussians :
	- Including both the original COMPASS fits (Solid) and, 
	- Our refits (Dotted)
Note that everything is normalized to unity in the integration over bT.

This plot should go in the paper. (Fig. 5 in arXiv:1401.2654)

Plot[{bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}],bt Finit2[bt]/NIntegrate[bbt Finit2[bbt],{bbt,0,15}],bt Ffina2[bt]/NIntegrate[bbt Ffina2[bbt],{bbt,0,15}]},{bt,.01,12},PlotStyle-> {{Red,Thick},{Blue,Thick},{Dotted,Thick},{Dotted,Thick},{Dashed,Purple,Thick}},LabelStyle-> {Large},PlotLegend-> {"(Q^2) = 1.1 (GeV^2) COMPASS","(Q^2) = 4.47 (GeV^2) COMPASS","(Q^2) = 1.1 (GeV^2) Our Fit","(Q^2) = 4.47 (GeV^2) Our Fit"},LegendPosition->{0.8,-0.3},LegendShadow-> None,PlotLabel-> "Initial and Final Gaussian Fits",AxesLabel->{"b",""}]

Plotting the initial and final COMPASS Guassian fits and our Kaplan fit (Dashed Curve) for the final Q.

This plot should go in the paper.  (Fig. 5 in arXiv:1401.2654)

Plot[{bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}],bt FFinakap[bt]},{bt,.01,12},PlotStyle-> {{Red,Thick},{Blue,Thick},{Dashed,Purple,Thick}},LabelStyle-> {Large},PlotLegend-> {"(Q^2) = 1.1 (GeV^2)","(Q^2) = 4.47 (GeV^2)","(Q^2) = 4.47 (GeV^2), Kaplan Fit"},LegendPosition->{0.8,-0.3},LegendShadow-> None,AxesLabel->{"b",""}]

Plot 4: Putting all the last few plots together.

Plot[{bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}],bt FFinakap[bt],bt Finit2[bt]/NIntegrate[bbt Finit2[bbt],{bbt,0,15}],bt Ffina2[bt]/NIntegrate[bbt Ffina2[bbt],{bbt,0,15}]},{bt,.01,12},PlotStyle-> {{Red,Thick},{Blue,Thick},{Dashed,Purple,Thick},{Dotted,Cyan,Thick},{Dotted,Green,Thick}},LabelStyle-> {Large},AxesLabel->{"b",""}]

Plot 5: Checking that the perturbation theory evolution kernals function with g2 = 0.  Note that for Q = Q0, it matches the Q0 Gaussian fit.

Plot[{bt fullPDF1[bt,Q0,.5,0]/NIntegrate[bbt fullPDF1[bbt,Q0,.5,0],{bbt,0,15}],bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt fullPDF1[bt,Qf,.5,0]/NIntegrate[bbt fullPDF1[bbt,Qf,.5,0],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}]},{bt,0.01,15},PlotStyle-> {{Purple,Thick,Dashed},{Red,Thick},{Green,Thick,Dashed},{Blue,Thick}},LabelStyle-> {Large},PlotLegend-> {"(Q^2) = 1.1 (GeV^2), Pert. Evolution","(Q^2) = 1.1 (GeV^2), Guassian Fit","(Q^2) = 4.47 (GeV^2), Pert. Evolution","(Q^2) = 4.47 (GeV^2), Gaussian Fit"},LegendPosition->{0.8,-0.3},LegendShadow-> None,AxesLabel->{"b",""}]

Plot 6: Plotting the initial and final COMPASS fits and our Kaplan fit along with the evolved distribution for a selection of values for (Subscript[g, 2]).
Here we use a conservative (Subscript[b, max  ])= 0.5 (GeV^(-2))  .

This plot should go in the paper. (Left plots in Fig. 8 in arXiv:1401.2654)

Plot[{bt fullPDF1[bt,Qf,.5,0]/NIntegrate[bbt fullPDF1[bbt,Qf,.5,0],{bbt,0,15}],bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt fullPDF1[bt,Qf,.5,Cevol]/NIntegrate[bbt fullPDF1[bbt,Qf,.5,Cevol],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}],bt FFinakap[bt],bt fullPDF1[bt,Qf,.5,0.1]/NIntegrate[bbt fullPDF1[bbt,Qf,.5,0.1],{bbt,0,15}],bt fullPDF1[bt,Qf,.5,0.7]/NIntegrate[bbt fullPDF1[bbt,Qf,.5,0.7],{bbt,0,15}]},{bt,0.01,15},PlotRange->{0,.25},PlotStyle-> {{Purple,Dotted},{Red,Thick},{Green,Dotted},{Blue,Thick},{Brown,Dashed,Thick},{Cyan,Dotted},{Yellow,Dotted}},LabelStyle-> {Large},PlotLegend-> {"(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0","(Q^2) = 1.1 (GeV^2), Guassian Fit","(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = (Subscript[C, evol])","(Q^2) = 4.47 (GeV^2), Gaussian Fit","(Q^2) = 4.47 (GeV^2), Kaplan Fit","(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0.1 (GeV^2)","(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0.7 (GeV^2)"},LegendPosition->{0.5,-0.5},LegendShadow-> None,AxesLabel->{"b",""}]

Plot 7: Plotting the initial and final COMPASS fits and our Kaplan fit along with the evolved distribution for a selection of values for (Subscript[g, 2]).
Here we use a liberal (Subscript[b, max  ])= 2.0 (GeV^(-2))  .

This plot should go in the paper. (Right plots in Fig. 8 in arXiv:1401.2654)

Plot[{bt fullPDF1[bt,Qf,2,0]/NIntegrate[bbt fullPDF1[bbt,Qf,2,0],{bbt,0,15}],bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt fullPDF1[bt,Qf,2,Cevol]/NIntegrate[bbt fullPDF1[bbt,Qf,2,Cevol],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}],bt FFinakap[bt],bt fullPDF1[bt,Qf,2,0.1]/NIntegrate[bbt fullPDF1[bbt,Qf,2,0.1],{bbt,0,15}],bt fullPDF1[bt,Qf,2,0.7]/NIntegrate[bbt fullPDF1[bbt,Qf,2,0.7],{bbt,0,15}]},{bt,0.01,15},PlotRange->{0,.25},PlotStyle-> {{Purple,Dotted},{Red,Thick},{Green,Dotted},{Blue,Thick},{Brown,Dashed,Thick},{Cyan,Dotted},{Yellow,Dotted}},LabelStyle-> {Large},PlotLegend-> {"(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0","(Q^2) = 1.1 (GeV^2), Guassian Fit","(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = (Subscript[C, evol])","(Q^2) = 4.47 (GeV^2), Gaussian Fit","(Q^2) = 4.47 (GeV^2), Kaplan Fit","(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0.1 (GeV^2)","(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0.7 (GeV^2)"},LegendPosition->{0.5,-0.5},LegendShadow-> None,AxesLabel->{"b",""}]

Plot 8: Plotting the initial and final COMPASS fits and our 
Kaplan fit along with the evolved distribution using for the logarithmic model of (Subscript[g, K]).
We use a  (Subscript[b, max  ])= 0.5 (GeV^(-2))  , and the parameters are (Subscript[g, 2 ])=0.1 (GeV^(-2  ))and aa = (Subscript[b, NP]) =2 (GeV^(-1))

This plot should go in the paper. (Fig. 5 in arXiv:1401.2654)

Plot[{bt fullPDF2[bt,Qf,.5,.1,2]/NIntegrate[bbt fullPDF2[bbt,Qf,.5,.1,2],{bbt,0,15}],bt Finit[bt]/NIntegrate[bbt Finit[bbt],{bbt,0,15}],bt Ffina[bt]/NIntegrate[bbt Ffina[bbt],{bbt,0,15}],bt FFinakap[bt]},{bt,0.01,15},PlotRange->{0,.25},PlotStyle-> {{Purple,Dotted,Thick},{Red,Thick},{Blue,Thick},{Brown,Dashed,Thick}},LabelStyle-> {Large},PlotLegend-> {"(Q^2) = 4.47 (GeV^2), (Subscript[g, 2]) = 0","(Q^2) = 1.1 (GeV^2), Guassian Fit","(Q^2) = 4.47 (GeV^2), Gaussian Fit","(Q^2) = 4.47 (GeV^2), Kaplan Fit"},LegendPosition->{0.5,-0.5},LegendShadow-> None,AxesLabel->{"b",""}]