Newer
Older
Master_thesis / jobs / 3583786 / test.py
@saslie saslie on 2 Apr 2019 1 KB ...
  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. import pickle as pkl
  9. import sys
  10. import time
  11. from helperfunctions import display_time
  12. import cmath as c
  13. import raremodel as rm
  14.  
  15. modl = rm.model()
  16.  
  17. load_set = True
  18.  
  19. draw = False
  20.  
  21. mode = "slim_points"
  22.  
  23. modl.mode = mode
  24.  
  25. min_bin_scaling = 100
  26.  
  27. set_size = 1e3
  28.  
  29. x_min = 211.0
  30.  
  31. x_max= 4781.0
  32.  
  33. jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg["jpsi"]
  34. modl.add_resonance(jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale)
  35.  
  36. psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
  37. modl.add_resonance(psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale)
  38.  
  39. modl.add_cusp(3525, 3e-7, 200, 7)
  40.  
  41. modl.normalize_pdf()
  42.  
  43. if load_set:
  44.  
  45. with open(r"./data/set_{0}_range({1}-{2}).pkl".format(int(set_size), int(x_min), int(x_max)), "rb") as input_file:
  46. set_dic = pkl.load(input_file)
  47. part_set = (set_dic["x_part"], set_dic["y_part"])
  48. counter_tot = set_dic["counter_tot"]
  49.  
  50. else:
  51. x_part, y_part, counter = modl.generate_points(set_size, mode = mode, min_bin_scaling = min_bin_scaling, verbose = 1)
  52. part_set = (x_part, y_part)
  53.  
  54. if draw:
  55. modl.draw_plots(part_set = part_set, counter = counter_tot, mode = mode, min_bin_scaling = min_bin_scaling)
  56. cusp_amp_scan_y = []
  57. cusp_amp_scan_x = np.linspace(0, 2*modl.cusp_amp, 100)
  58.  
  59. for i in cusp_amp_scan_x:
  60. cusp_amp_scan_y.append(modl.log_likelihood(x_part = set_dic["x_part"], cusp_amp = i))
  61. plt.clf()
  62.  
  63. plt.title("Cusp amp log_likelihood scan")
  64.  
  65. #plt.yscale("log")
  66.  
  67. #plt.ylim(0, 2*self._mean)
  68.  
  69. plt.grid()
  70.  
  71. plt.plot(cusp_amp_scan_x, cusp_amp_scan_y, label = "log_likelihood")
  72.  
  73. plt.legend()
  74.  
  75. plt.savefig("./cusp_amp_scan.png")
  76.  
  77. print(" pdf_and_parts.png created")
  78.  
  79. print(modl.log_likelihood(x_part = set_dic["x_part"], cusp_amp = modl.cusp_amp))
  80.  
  81. print("Run finished")