mbstar=5415.4;
mbstar0=5711.0;
B0M=5279.5;
BsM=5366.7;
BplusM=5279.3;
LbM=5619.4;
D0M=1864.8;
DstM=2010;
piM=139.6;
JpsiM=3096.9;
Psi2sM=3685.6;
kaonM=493.7;
KsM=497.6;
phiM=1019.5;
rhoM=775.26;
rhowidth=149.1;
omegaM=782.65;
omegawidth= 8.49;
muonM = 105.7;

squarkM = 95.0;
bquarkM = 4180.0;
cquarkM = 1275.0;

Bplustau = 1.638*10^(-12);

C7eff = -0.306;

C9eff = 4.211;
C10eff = -4.103;

GF = 1.1663787*10^(-5);
alphaew = 1.0/137.0;
Vts = 0.0394;
Vtb = 1.019;
numberofdecays= 5404696;

b0= {0.292 ,0.281, 0.150};
bplus= {0.466,-0.885,-0.213};
bT= {0.460 ,-1.089 ,-1.114};

NRBR= 4.37*10^(-7);
NRauc= 0.00133;

jpsim= 3096.0;
jpsiw=0.09 ;
jpsip=-1.5;
jpsia=184.39;
jpsiBR= 6.02*10^(-5);
jpsiauc= 0.2126825758464027;

psi2sm= 3686.0;
psi2sw=-0.3;
psi2sp = -1.5;
psi2sa = 23.07;
psi2sBR= 4.97*10^(-6);
psi2sauc= 2.802257483178487*10^(-10);

xmin = 2*muonM;
xmax = mB-mK-0.1;

mB = BplusM;
mK = kaonM;
tpos=(mB-mK)^2;
tzero=(mB+mK)*(Sqrt[mB]-Sqrt[mK])^2;
z[x_]:=(Sqrt[tpos-x^2]-Sqrt[tpos-tzero])/(Sqrt[tpos-x^2]+Sqrt[tpos-tzero])
Pplus[x_, y_]:=1-x^2/y^2
fT[x_]:= 1/Pplus[x,mbstar]*(∑_(j=0)^2)(bT[[j+1]](z[x]^j-(-1)^(j-3)*(j/3)*z[x]^3))
fplus[x_]:= 1/Pplus[x,mbstar]*(∑_(j=0)^2)(bplus[[j+1]](z[x]^j-(-1)^(j-3)*(j/3)*z[x]^3))
f0[x_]:= 1/Pplus[x,mbstar0]*(∑_(j=0)^2)(b0[[j+1]]z[x]^j)

Plot[{f0[x],fplus[x],fT[x]},{x,xmin,xmax},PlotRange->{0,5}]

p[x_]:=1/2*Sqrt[x^2-4*muonM^2]
delta[x_,m_,w_]:=p[x]/p[m]*m/x*w
resonance[x_,m_,w_,p_,a_]:=(a*m*w)/((m^2-x^2)-I*m*delta[x,m,w])*E^(I*p)

------------------------>

Plot[Abs[resonance[x, jpsim, jpsiw, jpsip, jpsia]+resonance[x, psi2sm,psi2sw, psi2sp,psi2sa]],{x,3000,3800}]

beta[x_]:= Sqrt[Abs[1-4*muonM^2/x^2]];
kabs[x_]:=Abs[Sqrt[mB^2+x^4/mB^2+mK^4/mB^2-2*(mB^2*mK^2+mK^2*x^2+mB^2*x^2)/mB^2]];
prefactor[x_] = GF^2*alphaew^2*Abs[Vtb*Vts]^2/(128*Pi^5)*kabs[x]*beta[x]

bracketleft[x_]:=2/3*kabs[x]^2*beta[x]^2*Abs[C10eff*fplus[x]]^2

bracketmiddle[x_]:=4*muonM^2*(mB^2-mK^2)^2/(x^2*mB^2)*Abs[C10eff*f0[x]]^2

bracketright[x_]:=kabs[x]^2*(1-1/3*beta[x]^2)*Abs[(C9eff+resonance[x, jpsim, jpsiw, jpsip, jpsia]+resonance[x, psi2sm,psi2sw, psi2sp,psi2sa])*fplus[x]+2*C7eff*(bquarkM+squarkM)/(mB+mK)*fT[x]]^2

dgamma[x_]:=prefactor[x]*(bracketleft[x]+bracketmiddle[x]+bracketright[x])*2*x

dgamma[500.0]

Plot[dgamma[x],{x,3000,4000},MaxRecursion->50000, Mesh -> All]

Integrate[dgamma[x],x,Assumptions->xmin<x<xmax]

FortranForm[%]