Comparison of A(b) - function using different treatments of  gK.
	From arXiv:1412.3820 (Updated Dec 14, 2014) 
References:
	 	[1]   Strong Coupling S. Bethke, Nucl.Phys.Proc.Suppl. 234 (2013) 229-234 
	 	[2]  “BLNY” Phys. Rev. D67, 073016 (2003), hep-ph/0212159
		[3]  “KN”     Phys. Lett. B633, 710 (2006), arXiv:hep-ph/0506225
		[4]  “BSY”  (B) Nucl. Phys. B806, 23 (2009), arXiv:0804.2408, (SY) Phys. Rev. D88, 034016 (2013), arXiv:1304.5037; Phys. Rev. D88, 114012 (2013), arXiv:1308.5003 
		[5]  “QZ”     Phys.Rev. D63 (2001) 114011, arXiv:hep-ph/0012348; arXiv:hep- ph/0205115

Clear["Global *"]

(* ============ Utilities ============ *)

changeRule[list_,rule:(sym_->_)]:=Append[FilterRules[list,Except[sym]],rule]

changeRules[list_, newRules :(_->_)..] := Sort[Join[FilterRules[list,Except@Alternatives[newRules][[All,1]]],{newRules}]]

(* ============ Parameterizations/fits ============*)

Some non-perturbative parametrizations from the literature.
	- PT: No non-perturbative gK, LO K calculation extended to bT =  ∞ .
          - BLNY: Brock-Landry-Nadolsky-Yuan fits, bmax =  0.5 GeV-1. [2]
	- KN: Konychev-Nadolsky fits, bmax =  1.5 GeV-1. [3]
	- BSY: Boer, Sun-Yuan, modified form of K with no bmax . [4]
	- QZ: Qiu-Zhang fits, bmax =  0.5 GeV-1. [5]

PT = {bmax->Infinity,
         gAfn-> Function[b,0],
         Afn-> Function[b, ACSS[b,bmax,gAfn]],
         PlotStyle->{Red,Thick}};

PTFixed1 = changeRules[{PT}, mu1->2,Lfn-> Function[b,L1[b,mu1] ], PlotStyle->Darker[Darker[Darker[Green]]] ];

KN = {bmax->1.5, a1-> 0.201, a2->0.184, a3-> -0.026,
        gAfn-> Function[b, 0.5*a2*b^2],
        Afn-> Function[b, ACSS[b,bmax,gAfn]],
        PlotStyle->Black};

BLNY = changeRules[{KN},bmax->0.5, a1-> 0.21, a2->0.68, a3-> -0.6*0.21, PlotStyle->Blue];

QZ = {bmax-> 0.5, alpha->0.61, g1->0.4, g2->10, PlotStyle->{Brown}};

SY = {Afn->Function[b,Anormal1[Q]],  PlotStyle->Brown};

(* ============ QCD coupling etc ============*)

Lambda = 0.339; (* === [1] === *)
nf = 3;
CF = 4/3;
beta0 := 11-2*nf/3;
beta1 := 102-38*nf/3;
stdloops = 2;
alpiN[mu_,loops_] := 2/( beta0 * Log[mu/Lambda]) +
                                        If[ loops > 1, (-beta1*Log[2*Log[mu/Lambda]] / (beta0^3 * Log[mu/Lambda]^2 )), 0];    
alpi[mu_] := alpiN[mu,stdloops]
mub[b_] := 2*Exp[-EulerGamma]/b
muQ[Q_] := Q

(* ============ CSS functions, etc ============*)

bstar[b_,bmax_] := If[ bmax == Infinity, b, b/Sqrt[1+b^2/bmax^2] ]

Expressions for A(b) - different depending on non-perturbative treatment

(* ==== Purely NLO perturbative expression from Eq.(37). ==== *)

Anormal1[b_] := CF * alpi[mub[b]] ;
Anormal2[b_] := 0.5*CF * alpi[mub[b]]^2 * ( 3 * (67/18-Pi^2/6) - nf*10/18);

(* ==== RG improved L. ==== *)

Anormal[b_] := Anormal1[b] + Anormal2[b]

(* ==== The standard CSS version, with an arbitrary form for (Subscript[g, K]). ==== *)

ACSS[b_, bmax_, gLfn_] := (Anormal[bstar[b,bmax]]*bstar[b,bmax]^2/b^2 +If[bmax==Infinity,0, gLfn[b]])

Aperturb[b_,bmax_] := Anormal[bstar[b,bmax]]*bstar[b,bmax]^2/b^2

(* ==== BSY form, but no (Subscript[g, K]) and no RG improved (Subscript[α, s]).[4] ==== *)

ABSY[b_,Q_] := CF * alpi[Q]

(* ==== Qiu-Zhang form, discrete jump from small to large (Subscript[b, T]). [5] ==== *)

AQZ[b_,bmax_,alpha_,g1_,g2_] := If[b<bmax,
                                                                                           Anormal[b],
                                                                                           Anormal[bmax]  + alpha* g1*b^(2*alpha)+g2*(b^2)]

A[b_,params_] := (Afn[b])//.params

(* ============ Plot definitions ============*)

OnePlot[A_, opts_] := Plot[ A,  {b,0.01,5}, PlotRange->{0,1}, 
                                                        BaseStyle->{FontSize->10},AxesLabel->{"(Subscript[b, T]) ((GeV^(-1)))","A((Subscript[b, T]))"}, opts ]

OnePlotLowb[A_, opts_] := LogLinearPlot[ A,  {b,0.01,1.0}, PlotRange->{0,0.4}, 
                                                        BaseStyle->{FontSize->10},AxesLabel->{"(Subscript[b, T]) ((GeV^(-1)))","A((Subscript[b, T]))"}, opts ]

OneLogPlot[A_, opts_] := LogLinearPlot[ A,  {b,0.01,4},AxesLabel->{"(Subscript[b, T]) ((GeV^(-1)))","A((Subscript[b, T]))"}, PlotRange->{0,1}, 
                                                        BaseStyle->{FontSize->10}, opts ]

OnePlotLong[A_, opts_] := Plot[ A,  {b,0.01,10}, PlotRange->{0,1},AxesLabel->{"(Subscript[b, T]) ((GeV^(-1)))","A((Subscript[b, T]))"}, opts ]

plotNormalLong = OnePlotLong[ Anormal[b], PlotStyle->{Red,Thick} ]

OnePlotA[params_] :=                                                    
   Plot[ A[b,params],  {b,0.01,5},PlotRange->{0,1}, PlotStyle-> {Red,Thick},BaseStyle->{FontSize->10},AxesLabel->{"(Subscript[b, T]) ((GeV^(-1)))","A((Subscript[b, T]))"}, Evaluate[ FilterRules[{params},Options[Plot]] ]  ]

OnePlotLowbA[params_] :=                                                    
   LogLinearPlot[ A[b,params],  {b,0.01,2.0},PlotRange->{0,0.4},AxesLabel->{"(Subscript[b, T]) ((GeV^(-1)))","A((Subscript[b, T]))"},  PlotStyle-> {Red,Thick},BaseStyle->{FontSize->10}, Evaluate[ FilterRules[{params},Options[Plot]] ]  ]

OnePlotA[PT]

OnePlotLowbA[PT]

Some values of 1- and 2- loop A:

Ashow[b_] := Print[ "b = ", b, " => A1, A2 = ", N[Anormal1[b]], ", ", N[Anormal2[b]] ]

Ashow[0.1]; Ashow[0.5] ; Ashow[1]; Ashow[2]

Plots of A(b):

(* ============= purely RG improved perturbative expression ========= *)
plotNormal = OnePlot[ Anormal[b], PlotStyle->{Red,Thickness[.005]}];

(* ============= purely RG improved perturbative expression at low (Subscript[b, T]) ========= *)
plotNormalLowb = OnePlotLowb[ Anormal[b], PlotStyle->{Red,Thickness[.005]} ];

(* ============= purely perturbative expression at fixed μ = (Subscript[C, 1]) / 2 ========= *)
plot1 = OnePlot[ Aperturb[b,2] , PlotStyle->{Darker[Darker[Darker[Green]]],Thick}];

(* ============= CSS expression with (Subscript[b, max]) = 1.5 (GeV^(-1)) and KN parameters [3] ========= *)
plotKN = OnePlot[ A[b, KN] , PlotStyle->{Black,Thick,DotDashed}];

(* ============= CSS expression with (Subscript[b, max]) = 1.5 (GeV^(-1)) and KN parameters [3], low (Subscript[b, T]) range ========= *)
plotKNLowb = OnePlotLowb[ A[b,KN] , PlotStyle->{Black,Thick,DotDashed}];

(* ============= CSS expression with (Subscript[b, max]) = 0.5 (GeV^(-1)) and BLNY parameters [2] ========= *)
plotBLNY = OnePlot[ A[b, BLNY] , PlotStyle->{Blue,DotDashed}];

(* ============= CSS expression with (Subscript[b, max]) = 0.5 (GeV^(-1)) and BLNY parameters [2], low (Subscript[b, T]) range ========= *)
plotBLNYLowb = OnePlotLowb[ A[b,BLNY] , PlotStyle->{Blue,DotDashed}];

(* ============= CSS expression with (Subscript[b, max]) = 0.5 (GeV^(-1)) and no non-perturbative part. ========= *)
plotPBLNY = OnePlot[ Aperturb[b,0.5] , PlotStyle->{Blue,Dashed}];

(* ============= CSS expression with (Subscript[b, max]) = 0.5 (GeV^(-1)) and no non-perturbative part, low (Subscript[b, T]) range ========= *)
plotPBLNYLowb = OnePlotLowb[ Aperturb[b,0.5] , PlotStyle->{Blue, Dashed}];

(* ============= CSS expression with (Subscript[b, max]) = 1.5 (GeV^(-1)) and no non-perturbative part. ========= *)
plotPKN = OnePlot[ Aperturb[b,1.5], PlotStyle->{Black, Dashed,Thick} ];

(* ============= CSS expression with (Subscript[b, max]) = 1.5 (GeV^(-1)) and no non-perturbative part, low (Subscript[b, T]) range ========= *)
plotPKNLowb = OnePlotLowb[Aperturb[ b,1.5], PlotStyle->{Black, Dashed, Thick} ];

(* ============= BSY expression with Q = 1 GeV ========= *)

plotBSY = OnePlot[ABSY[b,1.0], PlotStyle->{Thick,Brown}];

(* ============= BSY expression [4] with Q = 1 GeV, low (Subscript[b, T]) range ========= *)

plotBSYLowb = OnePlotLowb[ABSY[b,1.0], PlotStyle->{Thick,Brown}];

(* ============= BSY expression [4] with Q = 10 GeV ========= *)

plotBSY10 = OnePlot[ABSY[b,10], PlotStyle->{Blue, Dashed,Thick}];

(* ============= BSY expression [4] with Q = 10 GeV, low (Subscript[b, T]) range ========= *)

plotBSY10Lowb = OnePlotLowb[ABSY[b,10], PlotStyle->{Blue, Dashed,Thick}];

Print[ "Thick Red = perterbatuve CSS, no taming of large (Subscript[b, T]).\n" ,"Black dot-dashed = KN treatment, (Subscript[b, max])=1.5 (GeV^(-1)) .\n","Blue dot-dashed = BLNY treatment, (Subscript[b, max])= 0.5 (GeV^(-1)).\n","Black dashed = perturbative treatment but with (SubStar[b]) prescription and (Subscript[b, max])=1.5 (GeV^(-1)), no (Subscript[g, K]).\n",
             "Blue dashed = perturbative treatment but with (SubStar[b]) prescription and (Subscript[b, max])=0.5 (GeV^(-1)), no (Subscript[g, K])."]

Parametrizations with a quadratic (Subscript[g, K]) [1,2] compared with expressions with no (Subscript[g, K]) .

QuadgKplots = Show[plotKN, plotNormal, plotBLNY,plotPKN,plotPBLNY]
Export["QuadgKplots.pdf",QuadgKplots]

Same but with small (Subscript[b, T]) magnified.

QuadgKplotslowb =Show[plotKNLowb, plotNormalLowb, plotBLNYLowb,plotPKNLowb,plotPBLNYLowb]
Export["QuadgKplotslowb.pdf",QuadgKplotslowb]

BSY [4] treatment to be compared with quadratic [1,2] treatments.

Print[ "Thick Red = perterbatuve CSS, no taming of large (Subscript[b, T]).\n" ,"Black dot-dashed = KN treatment, (Subscript[b, max])=1.5 (GeV^(-1)) .\n","Solid Brown = BSY treatment, Q = 1.0 GeV.\n","Blue dashed = BSY treatment, Q = 10.0 GeV.\n"]

KNvsBSYplots = Show[plotKN, plotBLNY,plotNormal, plotBSY, plotBSY10]

Same but with small (Subscript[b, T]) magnified

KNvsBSYplotsLowbt =Show[plotKNLowb, plotBLNYLowb,plotNormalLowb, plotBSYLowb, plotBSY10Lowb]

Export["KNvsSYplots.pdf",KNvsBSYplots]
Export["KNvsSYplotsLowbt.pdf",KNvsBSYplotsLowbt]

Qiu - Zhang treatment [5], to be compared with quadratic [2,3] treatments.

plotQZ = OnePlot[AQZ[b,0.5,0.65,0.4,0], PlotStyle->{Brown,Thick}];

plotQZtogether =Show[plotQZ, plotKN,plotBLNY,plotNormal]
Export["plotQZ.pdf",plotQZtogether]