Simple Integral

Integrate[((Sin[θ])^2),{θ,0,2π}]

Using 1000 points

points=RandomVariate[UniformDistribution[{0,2 π}],1000];
samples=(2 π -0 )(((Sin[#])^2))&/@points;

Ihat=Mean[samples]

V=Variance[samples]

δIhat=(1/(Sqrt[Length[samples]]))(Sqrt[V])

Abs[Ihat-π]

Using 100’000 points

points=RandomVariate[UniformDistribution[{0,2 π}],100000];
samples=(2 π -0 )(((Sin[#])^2))&/@points;

Ihat=Mean[samples]

V=Variance[samples]

δIhat=(1/(Sqrt[Length[samples]]))(Sqrt[V])

Abs[Ihat-π]

BES at LEP

Cross-section for μμ production
[see e.g. Cheng-Li ch 12.1]

Integrate[(1+((Cos[θ])^2)) Sin[θ],{θ,0,π}]

GeVtopb=(GeV -> (((10^(3 )) Mev)/(197  Mev fm)))/.fm-> (Sqrt[(((10^(12 ))pbarn)/(100 ))])

Block[{G=1.166 (10^(-5)) (GeV^(-2)),mz=91.2 GeV,Γz=2.495 GeV,α=1/128,sw=(Sqrt[0.22]),gV,gA,shat},
gV=-(1/2)+2 (sw^2);
gA=-(1/2);
shat=shatinGeV2 (GeV^2);
χ=((G (mz^2))/(2 (Sqrt[2])π α)) (shat/(shat-(mz^2)+I mz Γz));
σhat[shatinGeV2_]=Simplify[(8/3)((π (α^2))/(2 shat))(1+2 Re[χ](gV^2)+((Abs[χ])^2)((((gV^2)+(gA^2)))^2))]/.GeVtopb]

Plots of cross-sections before BES convolution

Plot[σhat[(Ecm^2)]/pbarn/1000,{Ecm,85,97},AxesLabel->{"(Subscript[E, cm])[GeV]","(Subscript[σ, μμ])[nb]"}]

LogPlot[σhat[(Ecm^2)]/pbarn/1000,{Ecm,50,200}]

BES Luminosity

We take a Gaussian BES luminosity with variable spread “δ”

dL[δ_]=NormalDistribution[1,δ]

PDF[dL[δ]][τ]

Cross-Section calculation around the Z-pole

Clear[samples,points]

We evaluate the cross-section on a list of “points” that we will sample from the BES luminosity.
Cross-section in nbarn

σ[Ecm_,points_]:=Block[{samples},
samples=(σhat[(Ecm^2) #]/pbarn/1000)&/@points;
Around[Mean[samples],(1/(Sqrt[Length[samples]]))(Sqrt[Variance[samples]])]
]

BES δ=1%

points=RandomVariate[dL[0.01],1000];
Show[Plot[σhat[(Ecm^2)]/pbarn/1000,{Ecm,85,97},AxesLabel->{"(Subscript[E, cm])[GeV]","(Subscript[σ, μμ])[nb]"},PlotStyle->Dotted],
ListPlot[{#,σ[#,points]}&/@{85,86,87,88,89,90.0,90.5,91,91.25,91.5,92,92.5,93,94},PlotRange->{All,{0,6 (10^(-6))}}],PlotLabel->"δ=1%"]

BES δ=2%

points=RandomVariate[dL[0.02],1000];
Show[Plot[σhat[(Ecm^2)]/pbarn/1000,{Ecm,85,97},AxesLabel->{"(Subscript[E, cm])[GeV]","(Subscript[σ, μμ])[nb]"},PlotStyle->Dotted],
ListPlot[{#,σ[#,points]}&/@{85,86,87,88,89,90.0,90.5,91,91.25,91.5,92,92.5,93,94},PlotRange->{All,{0,6 (10^(-6))}}],PlotLabel->"δ=2%"]

Distribution of (Subscript[M, μμ]) in collisions at (Subscript[E, cm])=196 GeV

We create a Monte Carlo dataset of events.
In each event we record (Sqrt[(OverHat[s])]) =(Sqrt[τ]) (Subscript[E, cm]) = (Subscript[M, μμ]) and the weight w equal to the cross-section at that value of  (Sqrt[(OverHat[s])])

With BES δ=2%

Block[{Ecm=196},
points=RandomVariate[dL[0.02],10000];
events=({Ecm (Sqrt[Abs[#]]),σhat[(Ecm^2) #]/pbarn})&/@points;
]

We plot a weighted histogram using the Monte Carlo dataset.
Notice that we might have computed the error on the prediction in each bin, but we did not for brevity. 
You can check qualitatively that the result is stable by running multiple times or raising the number of points.

wd=WeightedData@@Transpose[events];
Histogram[wd,{150,200,2},"PDF",AxesLabel->{"(Subscript[M, μμ])[GeV]","((dσ)/(Subscript[dM, μμ]))[pb/GeV]"},PlotLabel->"δ=2%"]

Histogram[wd,{50,200,2},"PDF",AxesLabel->{"(Subscript[M, μμ])[GeV]","((dσ)/(Subscript[dM, μμ]))[pb/GeV]"},PlotLabel->"δ=2%"]

With BES δ=10%

Block[{Ecm=196},
points=RandomVariate[dL[0.1],10000];
events=({Ecm (Sqrt[Abs[#]]),σhat[(Ecm^2) #]/pbarn})&/@points;
]

wd=WeightedData@@Transpose[events];
Histogram[wd,{150,200,2},"PDF",AxesLabel->{"(Subscript[M, μμ])[GeV]","((dσ)/(Subscript[dM, μμ]))[pb/GeV]"},PlotLabel->"δ=10%"]

Histogram[wd,{50,200,2},"PDF",AxesLabel->{"(Subscript[M, μμ])[GeV]","((dσ)/(Subscript[dM, μμ]))[pb/GeV]"},PlotLabel->"δ=10%"]

Electron PDF at LEP

Create 2d ((Subscript[x, 1]), (Subscript[x, 2])) samples drawn from the PDF for electron and positron, and from them a 1d sample with ((Sqrt[τ]))

This one is the contribution to the PDF that comes from the real radiation emission.
The radiation cutoff “E” has been set to 100 GeV

Log[((((100))^2)/(((.511 (10^(-3))))^2))]

feRealRad[x_]=(1/(128 2 π))Log[((100^2)/(((.511 (10^(-3))))^2))] ((1+(x^2))/(1-x))

We will use this fe[x] to model the pdf in the region x<1-ϵ .

The threshold ϵ must be large enough for the integral of fe not to exceed one. Otherwise, this simplified approach must be replaced by a more complete treatment. 
Fortunately, the integral would be below one even if ϵ was as small as (10^(-6)), as shown below.

RealInt[ϵ_]=Assuming[{ϵ>0},Integrate[feRealRad[x],{x,0,1-ϵ}]]

LogLinearPlot[RealInt[ϵ],{ϵ,(10^(-8)),.1}]

The complete PDF, including virtual effects, is:
(Subscript[f, e])(x)= ((Subscript[f, e])^(Real Radiation)) + (Subscript[p, Virtual])δ(1-x)
Since (Subscript[f, e])(x) must be normalized (this is ensured to be true, if we compute the virtual directly, by the KLN theorem up to small corrections),
(Subscript[p, Virtual])=1- ((∫_0)^(1-ϵ)) ((Subscript[f, e])^(Real Radiation))=1-RealInt[ϵ]

We thus generate the x points as a mixture: N *(1-(Subscript[p, Virtual])) points drawn from the pdf:
((((Subscript[f, e])^(Real Radiation)) )/(RealInt[ϵ]))
and N * (Subscript[p, Virtual]) points with x exactly equal to one to account for the delta function.

We consider ϵ=1/100 and ϵ=1/1000 to make sure that the result is stable

If you play with this, you will notice that the smaller ϵ, the more events you need to generate for a stable prediction of the Z-peak in (Subscript[M, μμ]). 
This is because less and less events are generated at small x when ϵ decreases.

pvirtual[(1/100)]=1-RealInt[(1/100)]

pvirtual[(1/1000)]=1-RealInt[(1/1000)]

RealPDF[(1/100)]=ProbabilityDistribution[(1/(RealInt[1/100])) feRealRad[x],{x,0,1-(1/100)}];
RealPDF[(1/1000)]=ProbabilityDistribution[(1/(RealInt[1/1000])) feRealRad[x],{x,0,1-(1/1000)}];

Clear[points2d,points1d]

points2d[(1/100),100000]=
Block[{Np=100000,ϵ=(1/100)},
Transpose[{RandomSample[Join[RandomVariate[RealPDF[ϵ],IntegerPart[Np(1-pvirtual[ϵ])]],Table[1,IntegerPart[Np*pvirtual[ϵ]]]]],
RandomSample[Join[RandomVariate[RealPDF[ϵ],IntegerPart[Np(1-pvirtual[ϵ])]],Table[1,IntegerPart[Np*pvirtual[ϵ]]]]]}]];
points1d[(1/100),100000]=(#[[1]]  #[[2]] )&/@points2d[(1/100),100000];

points2d[(1/1000),100000]=
Block[{Np=100000,ϵ=(1/1000)},
Transpose[{RandomSample[Join[RandomVariate[RealPDF[ϵ],IntegerPart[Np(1-pvirtual[ϵ])]],Table[1,IntegerPart[Np*pvirtual[ϵ]]]]],
RandomSample[Join[RandomVariate[RealPDF[ϵ],IntegerPart[Np(1-pvirtual[ϵ])]],Table[1,IntegerPart[Np*pvirtual[ϵ]]]]]}]];
points1d[(1/1000),100000]=(#[[1]]  #[[2]]) &/@points2d[(1/1000),100000];

This is the distribution of x we obtain:

Histogram[points2d[(1/100),100000][[All,1]],{0,1.01,.001},"PDF",ScalingFunctions->"Log",AxesLabel->{"x"},PlotLabel->"ϵ=(1/100)"]

Histogram[points2d[(1/1000),100000][[All,1]],{0,1.01,.001},"PDF",ScalingFunctions->"Log",AxesLabel->{"x"},PlotLabel->"ϵ=(1/1000)"]

Compute the cross-section and generate events.

The code is exactly the same as for BES, only that the points used are not from a gaussian but the points1d generated before:

σ[Ecm_,points_]:=Block[{samples},
samples=(σhat[(Ecm^2) #]/pbarn/1000)&/@points;
Around[Mean[samples],(1/(Sqrt[Length[samples]]))(Sqrt[Variance[samples]])]
]

Show[Plot[σhat[(Ecm^2)]/pbarn/1000,{Ecm,85,97},AxesLabel->{"(Subscript[E, cm])[GeV]","(Subscript[σ, μμ])[nb]"},PlotStyle->Dotted],
ListPlot[{#,σ[#,points1d[(1/100),100000]]}&/@{85,86,87,88,89,90.0,90.5,91,91.25,91.5,92,92.5,93,94,95,96},PlotRange->{All,{0,6 (10^(-6))}}],PlotLabel->"ISR, with ϵ=(1/100)"]

Show[Plot[σhat[(Ecm^2)]/pbarn/1000,{Ecm,85,97},AxesLabel->{"(Subscript[E, cm])[GeV]","(Subscript[σ, μμ])[nb]"},PlotStyle->Dotted],
ListPlot[{#,σ[#,points1d[(1/1000),100000]]}&/@{85,86,87,88,89,90.0,90.5,91,91.25,91.5,92,92.5,93,94,95,96},PlotRange->{All,{0,6 (10^(-6))}}],PlotLabel->"ISR, with ϵ=(1/1000)"]

Generate events at 196. Also in this case, the code is the same

Block[{Ecm=196},
events=({Ecm (Sqrt[Abs[#]]),σhat[(Ecm^2) #]/pbarn})&/@points1d[(1/100),100000];
];
wd=WeightedData@@Transpose[events];
Histogram[wd,{50,200,4},"PDF",AxesLabel->{"(Subscript[M, μμ])[GeV]","((dσ)/(Subscript[dM, μμ]))[pb/GeV]"},PlotLabel->"ISR, with ϵ=(1/100)",AspectRatio->1,Frame->True,PlotRange->{All,{0,0.1}}]

Block[{Ecm=196},
events=({Ecm (Sqrt[Abs[#]]),σhat[(Ecm^2) #]/pbarn})&/@points1d[(1/1000),100000];
];
wd=WeightedData@@Transpose[events];
Histogram[wd,{50,200,4},"PDF",AxesLabel->{"(Subscript[M, μμ])[GeV]","((dσ)/(Subscript[dM, μμ]))[pb/GeV]"},PlotLabel->"ISR, with ϵ=(1/1000)",AspectRatio->1,Frame->True,PlotRange->{All,{0,0.1}}]