- 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
-
-
- 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, subscript):
-
- # x = q2
-
- ffplus = ""
-
- 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"]
- mB = self.mB
- mK = self.mK
-
- N = 3
-
- #some helperfunctions
-
- tpos = "(({0} - {1})*({0} - {1}))".format(mB, mK)
- tzero = "(({0} + {1})*(sqrt({0})-sqrt({1}))^2)".format(mB, mK)
-
- z_oben = "(sqrt({0} - pow(x,2)) - sqrt({0} - {1}))".format(tpos, tzero)
- z_unten = "(sqrt({0} - pow(x,2)) + sqrt({0} - {1}))".format(tpos, tzero)
- z = "({0})/({1})".format(z_oben, z_unten)
-
- #calculate f0
-
- if subscript == "0":
-
- prefactor_ff0 = "1/(1 - pow(x,2)/({0}^2))".format(mbstar0)
- ff0 = ""
-
- for i in range(N):
- ff0 += "{0}*({1}^{2})".format(b0[i], z, i)
- if i < (N -1):
- ff0 += "+"
-
- ff0 = prefactor_ff0 + "*" + ff0
-
- return "(" + ff0 + ")"
-
-
- #calculate f+
-
- elif subscript == "+":
-
- ffplus = ""
- prefactor_ffplus = "1/(1 - pow(x,2)/(pow({0},2)))".format(mbstar)
-
- for i in range(N):
- ffplus += "{0} * (pow({1},{2}) - (pow((-1),({2}-3)) * ({2}/3) * pow({1},3)))".format(bplus[i], z, i)
- if i < (N -1):
- ffplus += "+"
-
- ffplus = prefactor_ffplus + "*" + ffplus
-
- return "(" + ffplus + ")"
-
- #calculate fT
-
- else:
-
- ffT = ""
- prefactor_ffT = "1/(1 - pow(x,2)/(pow({0},2)))".format(mbstar)
-
- for i in range(N):
- ffplus += "{0} * (pow({1},{2}) - (pow((-1),({2}-3)) * ({2}/3) * pow({1},3)))".format(bT[i], z, i)
- if i < (N -1):
- ffT += "+"
-
- ffT = prefactor_ffT + "*" + ffT
- print(ffT)
-
- return "(" + ffT + ")"
-
-
-
- def axiv_nonres(self):
-
- # [0] = amp
-
- 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 = "sqrt(1. - 4. * pow({0},2) / x))".format(mmu)
-
- kabs = "(sqrt(pow({0},2) + pow(x,2)/pow({0},2) + pow({1},4)/pow({0},2) - 2 * (pow({0},2) * pow({1},2) + pow({1},2) * x + pow({0},2) * x) / pow({0},2)))".format(mB, mK)
-
- #prefactor in front of whole bracket
-
- prefactor1 = "(pow({0},2) * pow({1},2) * pow({2}*{3},2) * {4} * {5} / (128. * pow({6},5)))".format(GF, alpha_ew, Vtb, Vts, kabs, beta, np.pi)
-
- #left term in bracket
-
- bracket_left = "(2./3. * pow({0},2) * pow({1},2) * pow({2} * ({3}),2))".format(kabs, beta, C10eff, self.formfactor("+"))
-
- #middle term in bracket
-
- _top = "(4. * pow({0},2) * pow((pow({1},2) - pow({2},2)),2))".format(mmu, mB, mK)
-
- _under = "(x * pow({0},2))".format(mB)
-
- bracket_middle = "({0}/{1} * (sqrt((({2} * {3}) * ({2} * {3}))) *sqrt((({2} * {3}) * ({2} * {3})))))".format(_top, _under, C10eff, self.formfactor("0"))
-
- axiv = prefactor1 + "*" + "(" + bracket_left + "+" + bracket_middle + ") * 2 * sqrt(x) * [0]"
-
- return "(" + axiv + ")"
-
-
- def vec_nonres(self):
-
- # [1] = amp
-
- 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 = "sqrt(1. - 4. * pow({0},2) / x))".format(mmu)
-
- kabs = "(sqrt(pow({0},2) + pow(x,2)/pow({0},2) + pow({1},4)/pow({0},2) - 2 * (pow({0},2) * pow({1},2) + pow({1},2) * x + pow({0},2) * x) / pow({0},2)))".format(mB, mK)
-
- #prefactor in front of whole bracket
-
- prefactor1 = "(pow({0},2) * pow({1},2) * pow({2}*{3},2) * {4} * {5} / (128. * pow({6},5)))".format(GF, alpha_ew, Vtb, Vts, kabs, beta, np.pi)
-
- #right term in bracket
-
- prefactor2 = "(pow({0},2) * (1. - 1./3. * pow({1},2)))".format(kabs, beta)
-
- abs_bracket = "pow(abs({0} * {1} + 2 * {2} * ({3} + {4})/({5} + {6}) * {7}),2)".format(C9eff, self.formfactor("+"), C7eff, mb, ms, mB, mK, self.formfactor("T"))
-
- bracket_right = "{0} * {1}".format(prefactor2, abs_bracket)
-
- vec = "[1] * {0} * {1} * 2 * sqrt(x)".format(prefactor1, bracket_right)
-
- return "(" + vec + ")"
-
- def total_nonres(self):
-
- #[0] = amp axiv
- #[1] = vec amp
- #[2] = tot nonres amp
-
- #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()
-
- #axial verctor nonresonant part including C7
-
- axiv_nonres_part = self.axiv_nonres()
-
- #Complete term
-
- return "[2] * ({0} + {1})".format(vec_nonres_part, axiv_nonres_part)
-
-
- modl = model()
-
- mB = modl.mB
- mK = modl.mK
- mmu = modl.mmu
-
- tpos = "(({0} - {1})*({0} - {1}))".format(mB, mK)
- tzero = "(({0} + {1})*(sqrt({0})-sqrt({1}))^2)".format(mB, mK)
-
- z_oben = "(sqrt({0} - x) - sqrt({0} - {1}))".format(tpos, tzero)
- z_unten = "(sqrt({0} - x) + sqrt({0} - {1}))".format(tpos, tzero)
- z = "({0})/({1})".format(z_oben, z_unten)
-
- #print(z_oben)
-
- #z_oben = "sqrt(" + tpos + "-x) - sqrt(" + tpos + "-" + tzero + ")"
- #print(z_oben)
-
-
-
-
- dq = np.linspace(2*mmu, (mB-mK) - 0.1 ,1000)
-
- dq2 = []
-
- for i in dq:
- dq2.append(i**2)
-
- vec_fit = ROOT.TF1("vec_fit", modl.formfactor("T"), 0, (mB-mK) - 0.1)
-
- vec_fit.SetParameter(0,1)
- vec_fit.SetParameter(1,1)
- vec_fit.SetParameter(2,1)
-
- canvas = ROOT.TCanvas("canvas")
-
- canvas.cd()
-
- vec_fit.Draw()
-
- canvas.Print("root_test.pdf")
-
- canvas.Close()
-
-