Evaluation of Sudakov radiator

idycoupling = {as->a[μR]/csi*(1-a[μR]*β1/β0*Log[csi]/csi+(a[μR]*β1/β0*Log[csi]/csi)^2-a[μR]^2*(β1^2/β0^2*(Log[csi]-a[μR]*β0*t)/csi^2+β2*a[μR]*t/csi^2))/.csi->1-ρ/.t-> -ρ/a[μR]/β0(*a[μR]*β0*t/.t->2 Log[μ/μR]*)}

asCMW = A1 as (1+as/(2 Pi) A2/A1 + as^2/(2 Pi)^2 A3/A1)/.idycoupling/.a[μR]->α

(* integrate the "collinear" term: *)

ker1 = (asCMW/Pi*(1/2 ρ/α/β0)+as/Pi A1( B1/A1/2+ as/Pi B2/A1/4))/α/β0/.idycoupling/.a[μR]->α;
int1=Integrate[ker1,ρ];
int1=(int1/.ρ->0)-Limit[int1,ρ->  λ,Assumptions->{λ>0,α>0,β0<0}];

R1= Series[int1,{α,0,1}]//FullSimplify;

(* integrate the "soft" term: *)

ker2 = asCMW/Pi (-1/2  ρ/α/β0 +λ/α/β0)/α/β0;
int2 = Integrate[ker2,ρ];
int2 = Limit[int2,ρ-> λ,Assumptions->{λ>0}]-Limit[int2,ρ-> 2λ,Assumptions->{λ>0}];

R2 = Series[int2,{α,0,1}]//FullSimplify;

(* obtain the Sudakov radiator R *)

R = R1+R2;

(* restore resummation scale dependence *)

R = Series[R,{α,0,1}];

h1[λ_] = Collect[SeriesCoefficient[R,{α,0,-1}]//Expand,A1,FullSimplify]/λ β0

Anomalous dimensions and constants

constants={
(*beta function*)
β0->(1/(4π))(-4/3*nf*Tf+11/3*Ca),
β1->(1/(((π))^2))(((17*Ca^2-5*Ca*nf-3*Cf*nf))/24),
β2->(1/(((4π))^3))(325/54*nf^2-5033/18*nf+2857/2),
A1->2Cf,
A2-> Cf(((67/9)-((π^2)/3))Ca-(20/9)Tf nf),
B1-> -3 Cf,
(*Tf*)
Tf->(1/2),
Ca->3,
Cf->4/3,
nf->5,
as->0.118 (* strong coupling *)
};

exp[-R(τ)]

h1[λ_]=((A1 ((-1+2 λ) Log[1-2 λ]-2 (-1+λ) Log[1-λ]))/(2 π β0 λ));

Rp[L_]=-(D[X h1[as β0 X],X])/.X->L//.constants;

(Subscript[Σ, max])[L_] = Exp[L h1[as β0 L]]//.constants;

Analytic Result for the transfer function

Fanalytic[L_] = Exp[-EulerGamma Rp[L]]/Gamma[1+Rp[L]];

Computation of the transfer function (efficient way)

number of events to generate, and IR cutoff

Nevents = 1000;  
ln1Cutoff=10;

Array containing observable’s values at which the cross section is computed

v=Array[# 1/100&,30]

start with ln1zeta = Log[1/(Subscript[ζ, i])] = Log[1/(Subscript[ζ, 1])]=0, only the first emission is present

F=Array[0&,Length[v]];
errF=Array[0&,Length[v]];
AbsoluteTiming[For[ibin=1,ibin≤ Length[v],ibin++,
For[iev=1,iev≤ Nevents, iev++,
ln1zeta= 0;
zetaSum = 0;
nems=0;
While[ln1zeta < ln1Cutoff,
zetaSum = zetaSum + Exp[-ln1zeta];
nems = nems + 1;
r=0;
While[r == 0,r=RandomReal[]];
u = -Log[r]/Rp[Log[1/v[[ibin]]]];
ln1zeta = ln1zeta + u;
];
F[[ibin]]= F[[ibin]]+ Exp[-Rp[Log[1/v[[ibin]]]] Log[zetaSum]]/Nevents;
errF[[ibin]]= errF[[ibin]]+ (Exp[-Rp[Log[1/v[[ibin]]]] Log[zetaSum]])^2/Nevents;
];
errF[[ibin]] = Sqrt[Abs[errF[[ibin]]-F[[ibin]]^2]]/Nevents;]]

Compare to the analytic solution

Histo=Table[{v[[i]],F[[i]]},{i,1,Length[v]}];
Show[{ListPlot[Histo],Plot[Fanalytic[Log[1/τ]],{τ,0,v[[Length[v]]]}]}]

Computation of the transfer function (inefficient way)

number of events to generate, and IR cutoff

Nevents = 1000;  
ln1Cutoff=10;

Array containing observable’s values at which the cross section is computed

v=Array[# 1/100&,30]

v=Array[# 1/1000&,100]//N

start with ln1zeta = Log[1/(Subscript[ζ, i])] = Log[1/(Subscript[ζ, 1])]=0, only the first emission is present

F=Array[0&,Length[v]];
AbsoluteTiming[For[ibin=1,ibin≤ Length[v],ibin++,
For[iev=1,iev≤ Nevents, iev++,
ln1zeta= 0;
zetaSum = 0;
nems=0;
(* generate first emission *)
r=0;
While[r == 0,r=RandomReal[]];
ln1zeta= -Log[r]/Rp[Log[1/v[[ibin]]]];
While[ln1zeta < ln1Cutoff,
zetaSum = zetaSum + Exp[-ln1zeta];
nems = nems + 1;
r=0;
While[r == 0,r=RandomReal[]];
u = -Log[r]/Rp[Log[1/v[[ibin]]]];
ln1zeta = ln1zeta + u;
];
If[zetaSum ≤ 1, F[[ibin]]= F[[ibin]]+ 1/Nevents];
]]]

Compare to the analytic solution

Histo=Table[{v[[i]],F[[i]]},{i,1,Length[v]}];
Show[{ListLogPlot[Histo,PlotRange->All],LogPlot[Fanalytic[Log[1/τ]],{τ,0,v[[Length[v]]]},PlotRange->All]}]

Thrust resummed distribution and comparison to fixed order

(Subscript[Σ, res])[τ_] := (Subscript[Σ, max])[Log[1/τ]] Fanalytic[Log[1/τ]];

compute the differential distribution and compare to the fixed-order at LO (from Ellis-Stirling-Webber Eq. 3.44)

(Subscript[dΣ, res])[τ_]:=D[(Subscript[Σ, res])[x],x]/.x->τ

(Subscript[dΣ, LO])[τ_] := Cf as/(2 Pi) (2 (3 T^2-3 T+2)/(T(1-T))Log[(2T-1)/(1-T)]-3 (3T -2)(2-T)/(1-T)) /.T->1-τ//.constants

Plot[{(Subscript[dΣ, LO])[τ],(Subscript[dΣ, res])[τ]},{τ,0.001,1/3},PlotRange->{{0,40}},AxesLabel->{"τ","(1/σ)\!\(\*FractionBox[\(dσ\), \(d\\ τ\)]\)"},PlotLegends->Placed[{"LO distribution","Resummed distribution"},{Center}]]