Newer
Older
Master_thesis / root_test.py
  1. import ROOT
  2. #from ROOT import TTree, TFile, Double
  3. import numpy as np
  4. from pdg_const import pdg
  5. import matplotlib
  6. matplotlib.use("Qt5Agg")
  7. import matplotlib.pyplot as plt
  8.  
  9. class model:
  10. def __init__(self):
  11.  
  12. self.mmu = pdg['muon_M']
  13. self.mb = pdg["bquark_M"]
  14. self.ms = pdg["squark_M"]
  15. self.mK = pdg["Ks_M"]
  16. self.mB = pdg["Bplus_M"]
  17.  
  18. self.C7eff = pdg["C7eff"]
  19. self.C9eff = pdg["C9eff"]
  20. self.C10eff = pdg["C10eff"]
  21. #self.C1 = pdg["C1"]
  22. #self.C2 = pdg["C2"]
  23. #self.C3 = pdg["C3"]
  24. #self.C4 = pdg["C4"]
  25. self.GF = pdg["GF"] #Fermi coupling const.
  26. self.alpha_ew = pdg["alpha_ew"]
  27. self.Vts = pdg["Vts"]
  28. self.Vtb = pdg["Vtb"]
  29. def formfactor(self, subscript):
  30. # x = q2
  31. ffplus = ""
  32. if subscript != "0" and subscript != "+" and subscript != "T":
  33. raise ValueError('Wrong subscript entered, choose either 0, + or T')
  34. #get constants
  35. mh = self.mK
  36. mbstar0 = pdg["mbstar0"]
  37. mbstar = pdg["mbstar"]
  38. b0 = pdg["b0"]
  39. bplus = pdg["bplus"]
  40. bT = pdg["bT"]
  41. mB = self.mB
  42. mK = self.mK
  43. N = 3
  44. #some helperfunctions
  45. tpos = "(({0} - {1})*({0} - {1}))".format(mB, mK)
  46. tzero = "(({0} + {1})*(sqrt({0})-sqrt({1}))^2)".format(mB, mK)
  47. z_oben = "(sqrt({0} - pow(x,2)) - sqrt({0} - {1}))".format(tpos, tzero)
  48. z_unten = "(sqrt({0} - pow(x,2)) + sqrt({0} - {1}))".format(tpos, tzero)
  49. z = "({0})/({1})".format(z_oben, z_unten)
  50. #calculate f0
  51. if subscript == "0":
  52. prefactor_ff0 = "1/(1 - pow(x,2)/({0}^2))".format(mbstar0)
  53. ff0 = ""
  54. for i in range(N):
  55. ff0 += "{0}*({1}^{2})".format(b0[i], z, i)
  56. if i < (N -1):
  57. ff0 += "+"
  58. ff0 = prefactor_ff0 + "*" + ff0
  59. return "(" + ff0 + ")"
  60. #calculate f+
  61. elif subscript == "+":
  62. ffplus = ""
  63. prefactor_ffplus = "1/(1 - pow(x,2)/(pow({0},2)))".format(mbstar)
  64. for i in range(N):
  65. ffplus += "{0} * (pow({1},{2}) - (pow((-1),({2}-3)) * ({2}/3) * pow({1},3)))".format(bplus[i], z, i)
  66. if i < (N -1):
  67. ffplus += "+"
  68. ffplus = prefactor_ffplus + "*" + ffplus
  69. return "(" + ffplus + ")"
  70. #calculate fT
  71. else:
  72. ffT = ""
  73. prefactor_ffT = "1/(1 - pow(x,2)/(pow({0},2)))".format(mbstar)
  74. for i in range(N):
  75. ffplus += "{0} * (pow({1},{2}) - (pow((-1),({2}-3)) * ({2}/3) * pow({1},3)))".format(bT[i], z, i)
  76. if i < (N -1):
  77. ffT += "+"
  78. ffT = prefactor_ffT + "*" + ffT
  79. print(ffT)
  80. return "(" + ffT + ")"
  81. def axiv_nonres(self):
  82. # [0] = amp
  83. GF = self.GF
  84. alpha_ew = self.alpha_ew
  85. Vtb = self.Vtb
  86. Vts = self.Vts
  87. C10eff = self.C10eff
  88. mmu = self.mmu
  89. mb = self.mb
  90. ms = self.ms
  91. mK = self.mK
  92. mB = self.mB
  93. #Some helperfunctions
  94. beta = "sqrt(1. - 4. * pow({0},2) / x))".format(mmu)
  95. 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)
  96. #prefactor in front of whole bracket
  97. 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)
  98. #left term in bracket
  99. bracket_left = "(2./3. * pow({0},2) * pow({1},2) * pow({2} * ({3}),2))".format(kabs, beta, C10eff, self.formfactor("+"))
  100. #middle term in bracket
  101. _top = "(4. * pow({0},2) * pow((pow({1},2) - pow({2},2)),2))".format(mmu, mB, mK)
  102. _under = "(x * pow({0},2))".format(mB)
  103. bracket_middle = "({0}/{1} * (sqrt((({2} * {3}) * ({2} * {3}))) *sqrt((({2} * {3}) * ({2} * {3})))))".format(_top, _under, C10eff, self.formfactor("0"))
  104. axiv = prefactor1 + "*" + "(" + bracket_left + "+" + bracket_middle + ") * 2 * sqrt(x) * [0]"
  105. return "(" + axiv + ")"
  106.  
  107. def vec_nonres(self):
  108. # [1] = amp
  109. GF = self.GF
  110. alpha_ew = self.alpha_ew
  111. Vtb = self.Vtb
  112. Vts = self.Vts
  113. C7eff = self.C7eff
  114. C9eff = self.C9eff
  115. mmu = self.mmu
  116. mb = self.mb
  117. ms = self.ms
  118. mK = self.mK
  119. mB = self.mB
  120. #Some helperfunctions
  121. beta = "sqrt(1. - 4. * pow({0},2) / x))".format(mmu)
  122. 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)
  123. #prefactor in front of whole bracket
  124. 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)
  125. #right term in bracket
  126. prefactor2 = "(pow({0},2) * (1. - 1./3. * pow({1},2)))".format(kabs, beta)
  127. 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"))
  128. bracket_right = "{0} * {1}".format(prefactor2, abs_bracket)
  129. vec = "[1] * {0} * {1} * 2 * sqrt(x)".format(prefactor1, bracket_right)
  130. return "(" + vec + ")"
  131. def total_nonres(self):
  132. #[0] = amp axiv
  133. #[1] = vec amp
  134. #[2] = tot nonres amp
  135. #Get constants
  136. GF = self.GF
  137. alpha_ew = self.alpha_ew
  138. Vtb = self.Vtb
  139. Vts = self.Vts
  140. C10eff = self.C10eff
  141. C9eff = self.C9eff
  142. C7eff = self.C7eff
  143. mmu = self.mmu
  144. mb = self.mb
  145. ms = self.ms
  146. mK = self.mK
  147. mB = self.mB
  148. #vector nonresonant part
  149. vec_nonres_part = self.vec_nonres()
  150. #axial verctor nonresonant part including C7
  151. axiv_nonres_part = self.axiv_nonres()
  152. #Complete term
  153. return "[2] * ({0} + {1})".format(vec_nonres_part, axiv_nonres_part)
  154.  
  155. modl = model()
  156.  
  157. mB = modl.mB
  158. mK = modl.mK
  159. mmu = modl.mmu
  160.  
  161. tpos = "(({0} - {1})*({0} - {1}))".format(mB, mK)
  162. tzero = "(({0} + {1})*(sqrt({0})-sqrt({1}))^2)".format(mB, mK)
  163.  
  164. z_oben = "(sqrt({0} - x) - sqrt({0} - {1}))".format(tpos, tzero)
  165. z_unten = "(sqrt({0} - x) + sqrt({0} - {1}))".format(tpos, tzero)
  166. z = "({0})/({1})".format(z_oben, z_unten)
  167.  
  168. #print(z_oben)
  169.  
  170. #z_oben = "sqrt(" + tpos + "-x) - sqrt(" + tpos + "-" + tzero + ")"
  171. #print(z_oben)
  172.  
  173.  
  174.  
  175.  
  176. dq = np.linspace(2*mmu, (mB-mK) - 0.1 ,1000)
  177.  
  178. dq2 = []
  179.  
  180. for i in dq:
  181. dq2.append(i**2)
  182. vec_fit = ROOT.TF1("vec_fit", modl.formfactor("T"), 0, (mB-mK) - 0.1)
  183.  
  184. vec_fit.SetParameter(0,1)
  185. vec_fit.SetParameter(1,1)
  186. vec_fit.SetParameter(2,1)
  187.  
  188. canvas = ROOT.TCanvas("canvas")
  189.  
  190. canvas.cd()
  191.  
  192. vec_fit.Draw()
  193.  
  194. canvas.Print("root_test.pdf")
  195.  
  196. canvas.Close()
  197.  
  198.