import ROOT
#from ROOT import TTree, TFile, Double
import numpy as np
from pdg_const import pdg
import matplotlib
matplotlib.use("Qt5Agg")
import matplotlib.pyplot as plt
import pickle as pkl
class model:
def __init__(self):
self.mmu = pdg['muon_M']
self.mb = pdg["bquark_M"]
self.ms = pdg["squark_M"]
self.mK = pdg["Ks_M"]
self.mB = pdg["Bplus_M"]
self.C7eff = pdg["C7eff"]
self.C9eff = pdg["C9eff"]
self.C10eff = pdg["C10eff"]
#self.C1 = pdg["C1"]
#self.C2 = pdg["C2"]
#self.C3 = pdg["C3"]
#self.C4 = pdg["C4"]
self.GF = pdg["GF"] #Fermi coupling const.
self.alpha_ew = pdg["alpha_ew"]
self.Vts = pdg["Vts"]
self.Vtb = pdg["Vtb"]
def formfactor(self, q2, subscript):
#check if subscript is viable
if subscript != "0" and subscript != "+" and subscript != "T":
raise ValueError('Wrong subscript entered, choose either 0, + or T')
#get constants
mh = self.mK
mbstar0 = pdg["mbstar0"]
mbstar = pdg["mbstar"]
b0 = pdg["b0"]
bplus = pdg["bplus"]
bT = pdg["bT"]
N = 3
#some helperfunctions
tpos = (self.mB - self.mK)**2
tzero = (self.mB + self.mK)*(np.sqrt(self.mB)-np.sqrt(self.mK))**2
z_oben = np.sqrt(tpos - q2) - np.sqrt(tpos - tzero)
z_unten = np.sqrt(tpos - q2) + np.sqrt(tpos - tzero)
z = z_oben/z_unten
#calculate f0
if subscript == "0":
prefactor = 1/(1 - q2/(mbstar0**2))
_sum = 0
for i in range(N):
_sum += b0[i]*(z**i)
return prefactor * _sum
#calculate f+ or fT
else:
prefactor = 1/(1 - q2/(mbstar**2))
_sum = 0
if subscript == "T":
b = bT
else:
b = bplus
for i in range(N):
_sum += b[i] * (z**i - ((-1)**(i-N)) * (i/N) * z**N)
return prefactor * _sum
def axiv_nonres(self, q2):
GF = self.GF
alpha_ew = self.alpha_ew
Vtb = self.Vtb
Vts = self.Vts
C10eff = self.C10eff
mmu = self.mmu
mb = self.mb
ms = self.ms
mK = self.mK
mB = self.mB
#Some helperfunctions
beta = np.sqrt(np.abs(1. - 4. * self.mmu**2. / q2))
kabs = np.sqrt(mB**2 + q2**2/mB**2 + mK**4/mB**2 - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
#prefactor in front of whole bracket
prefactor1 = GF**2. *alpha_ew**2. * (np.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)
#left term in bracket
bracket_left = 2./3. * kabs**2 * beta**2 * np.abs(C10eff*self.formfactor(q2, "+"))**2
#middle term in bracket
_top = 4. * mmu**2 * (mB**2 - mK**2) * (mB**2 - mK**2)
_under = q2 * mB**2
bracket_middle = _top/_under * np.abs(C10eff * self.formfactor(q2, "0"))**2
return prefactor1 * (bracket_left + bracket_middle) * 2 * np.sqrt(q2)
def vec_nonres(self, q2):
GF = self.GF
alpha_ew = self.alpha_ew
Vtb = self.Vtb
Vts = self.Vts
C7eff = self.C7eff
C9eff = self.C9eff
mmu = self.mmu
mb = self.mb
ms = self.ms
mK = self.mK
mB = self.mB
#Some helperfunctions
beta = np.sqrt(np.abs(1. - 4. * self.mmu**2. / q2))
kabs = np.sqrt(mB**2 + q2**2/mB**2 + mK**4/mB**2 - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
#prefactor in front of whole bracket
prefactor1 = GF**2. *alpha_ew**2. * (np.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)
#right term in bracket
prefactor2 = kabs**2 * (1. - 1./3. * beta**2)
abs_bracket = np.abs(C9eff * self.formfactor(q2, "+") + 2 * C7eff * (mb + ms)/(mB + mK) * self.formfactor(q2, "T"))**2
bracket_right = prefactor2 * abs_bracket
return prefactor1 * bracket_right * 2 * np.sqrt(q2)
def total_nonres(self, q2):
#Get constants
GF = self.GF
alpha_ew = self.alpha_ew
Vtb = self.Vtb
Vts = self.Vts
C10eff = self.C10eff
C9eff = self.C9eff
C7eff = self.C7eff
mmu = self.mmu
mb = self.mb
ms = self.ms
mK = self.mK
mB = self.mB
#vector nonresonant part
vec_nonres_part = self.vec_nonres(q2)
#axial verctor nonresonant part including C7
axiv_nonres_part = self.axiv_nonres(q2)
#Complete term
return vec_nonres_part + axiv_nonres_part
modl = model()
#f0 = modl.formfactor(10000, "+")
#print(f0)
_dgammadq2 = modl.total_nonres(10)
mB = modl.mB
mK = modl.mK
mmu = modl.mmu
#print((mB-mK))
dq = np.linspace(2*mmu, (mB-mK) - 0.1 ,1000)
dq2 = []
for i in dq:
dq2.append(i**2)
ff_plus = []
ff_T = []
ff_0 = []
dgamma_axiv_nonres = []
dgamma_vec_nonres = []
dgamma_tot = []
for i in dq:
ff_0.append(modl.formfactor(i**2, "0"))
ff_T.append(modl.formfactor(i**2, "T"))
ff_plus.append(modl.formfactor(i**2, "+"))
for i in dq:
dgamma_axiv_nonres.append(modl.axiv_nonres(i**2))
dgamma_vec_nonres.append(modl.vec_nonres(i**2))
dgamma_tot.append(modl.total_nonres(i**2))
plt.plot(dq2, ff_0, label = "0")
plt.plot(dq2, ff_T, label = "T")
plt.plot(dq2, ff_plus, label = "+")
plt.grid()
plt.legend()
plt.savefig("ff.png")
plt.clf()
plt.plot(dq, dgamma_axiv_nonres, label = "axiv")
plt.plot(dq, dgamma_vec_nonres, label = "vec")
plt.grid()
plt.legend()
plt.savefig("vec_axiv.png")
plt.clf()
plt.plot(dq, dgamma_tot, label = "total")
#print(max(dgamma_tot))
plt.grid()
plt.legend()
plt.savefig("tot.png")
#Generate a random set
set_size = 10000
x = []
y = []
x_part = []
y_part = []
print("Generating set of size {}...".format(set_size))
#print(len(y))
#ROOT.TRandom1().SetSeed(0)
while len(x_part) < set_size:
x.append(ROOT.TRandom1().Uniform(2*mmu, (mB-mK)))
y.append(ROOT.TRandom1().Uniform(0, max(dgamma_tot)*1.05))
if y[-1] < modl.total_nonres(x[-1]**2):
x_part.append(x[-1])
y_part.append(y[-1])
#progress calls
if len(x_part)%int(set_size/10) == 0:
print("{0}/{1} completed".format(len(x_part), set_size))
print("{0} out of {1} particles chosen!".format(len(x_part), len(x)))
print("Set generated!")
part_set = {"q" : x_part, "y": y_part}
pkl.dump( part_set, open("set.pkl" , "wb" ) )
print("Set saved!")
plt.clf()
plt.plot(x, y, label = "total", marker = ".", linewidth = 0)
plt.plot(dq, dgamma_tot, label = "pdf")
plt.plot(x_part, y_part, label = "chosen", marker = ".", linewidth = 0)
plt.grid()
plt.legend()
plt.savefig("set.png")
print("test finished")
