Skip to content
Snippets Groups Projects
Commit 04ed4101 authored by Markus Kuhlmann's avatar Markus Kuhlmann
Browse files

BBUnstable phasespace factor numerically stable

parent 7a38283c
No related branches found
No related tags found
No related merge requests found
......@@ -25,8 +25,16 @@ function cm(real, imag, m, µ)
function bbcm(sreal,simag,m11,m21,m12,m22,mR1,fR1,mR2,fR2,epsilon)
a1 = (m11+m21)^2
a2 = (m12+m22)^2
function scaler(u)
if u == 1.0
return (1.0/((1.0000000000001-u)^2))
end
return (1.0/((1.0-u)^2))
end
function sprime(a, t)
if t == 1.0
return (a+(t/(1.00000000001-t)))
end
return (a+(t/(1.0-t)))
end
function ds(t, m1, m2, fR, mR, a)
......@@ -37,19 +45,19 @@ function bbcm(sreal,simag,m11,m21,m12,m22,mR1,fR1,mR2,fR2,epsilon)
end
function integrand(x,f)
if fR1 == 0.0
tmp = (((1.0/pi)*fR2^2*imag(Sigma(a2,x[2], m12, m22)))/abs2(ds(x[2],m12, m22, fR2, mR2, a2)))*cm(sreal, simag, sqrt(sprime(a2, x[2])), mR1)*(1.0/((1.0-x[2])^2))
tmp = (((1.0/pi)*fR2^2*imag(Sigma(a2,x[2], m12, m22)))/abs2(ds(x[2],m12, m22, fR2, mR2, a2)))*cm(sreal, simag, sqrt(sprime(a2, x[2])), mR1)*scaler(x[1])
elseif fR2 == 0.0
tmp = (((1.0/pi)*fR1^2*imag(Sigma(a1,x[1], m11, m21)))/abs2(ds(x[1],m11, m21, fR1, mR1, a1)))*cm(sreal, simag, sqrt(sprime(a1, x[1])), mR2)*(1.0/((1.0-x[1])^2))
tmp = (((1.0/pi)*fR1^2*imag(Sigma(a1,x[1], m11, m21)))/abs2(ds(x[1],m11, m21, fR1, mR1, a1)))*cm(sreal, simag, sqrt(sprime(a1, x[1])), mR2)*scaler(x[2])
else
tmp = ((1.0/(pi^2))*fR1^2*imag(Sigma(a1,x[1], m11, m21)))/abs2(ds(x[1],m11, m21, fR1, mR1, a1))*cm(sreal, simag, sqrt(sprime(a1, x[1])), sqrt(sprime(a2, x[2])))*(1.0/((1.0-x[1])^2))*(1.0/((1.0-x[2])^2))*fR2^2*imag(Sigma(a2,x[2], m12, m22))/abs2(ds(x[2],m12, m22, fR2, mR2, a2))
tmp = ((1.0/(pi^2))*fR1^2*imag(Sigma(a1,x[1], m11, m21)))/abs2(ds(x[1],m11, m21, fR1, mR1, a1))*cm(sreal, simag, sqrt(sprime(a1, x[1])), sqrt(sprime(a2, x[2])))*scaler(x[1])*scaler(x[2])*fR2^2*imag(Sigma(a2,x[2], m12, m22))/abs2(ds(x[2],m12, m22, fR2, mR2, a2))
end
f[1], f[2] = reim(tmp)
end
if fR1 == 0.0 && fR2!=0.0
if fR1 == 0.0 && fR2 != 0.0
result = cuhre(integrand,1,2)
elseif fR2 == 0.0 && fR1!=0.0
elseif fR2 == 0.0 && fR1!= 0.0
result = cuhre(integrand,1,2)
elseif fR2 == 0.0 && fR1==0.0
elseif fR2 == 0.0 && fR1 == 0.0
result = cm(sreal, simag, mR1, mR2)
else
result = cuhre(integrand,2,2)
......@@ -60,17 +68,3 @@ function bbcm(sreal,simag,m11,m21,m12,m22,mR1,fR1,mR2,fR2,epsilon)
return resultCmplx
end
end
using PhaseSpace
for i = 0:400
for j = 0:0
iscale = i/100.0
jscale = j/100.0
#result = bbcm(iscale,jscale,0.14,0.14,0.14,0.14,0.758287544,0.442718872,0.14,0.0,0.0)
result = bbcm(iscale,jscale,0.14,0.14,0.14,0.14,0.77,0.0,0.14,0.0,0.0)
betrag=abs(complex(result[1],result[2]))
println(sqrt(iscale), " ",sqrt(jscale), " ", betrag)
end
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment