Newer
Older
Master_thesis / test.py
  1. import numpy as np
  2. from pdg_const import pdg
  3. import matplotlib
  4. import matplotlib.pyplot as plt
  5. import pickle as pkl
  6. import sys
  7. import time
  8. from helperfunctions import display_time
  9. import cmath as c
  10. import raremodel as rm
  11. import scipy.integrate as integrate
  12.  
  13. modl = rm.model()
  14.  
  15. load_set = True
  16. draw = True
  17. fit = False
  18.  
  19. mode = "slim_points"
  20.  
  21. modl.mode = mode
  22.  
  23. set_size = 4e5
  24. nonres_set_size = 44000
  25. # nonres_set_size = 1000
  26. nr_of_toys = 1
  27.  
  28. x_min = 3150.0
  29. x_max= 3650.0
  30.  
  31. x_max = modl.x_max
  32. x_min = modl.x_min
  33.  
  34. NR_BR = pdg["NR_BR"]
  35. modl.add_nonres(NR_BR)
  36.  
  37. jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg["jpsi"]
  38. jpsi_BR = pdg["jpsi_BR"]
  39. modl.add_resonance(jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale, "jpsi", jpsi_BR)
  40.  
  41. psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
  42. psi2s_BR = pdg["psi2s_BR"]
  43. modl.add_resonance(psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale, "psi2s", psi2s_BR)
  44.  
  45. # modl.add_cusp(3550, 3e-7, 200, 7, 3e-7)
  46.  
  47. # print(integrate.quad(modl.total_pdf, modl.x_min**2, modl.x_max**2, limit = 3000))
  48.  
  49. # modl.param_val[0] = 1e7
  50.  
  51. # modl.normalize_pdf()
  52.  
  53. # print(modl.param_val[0])
  54.  
  55. # print(integrate.quad(modl.total_pdf, modl.x_min**2, modl.x_max**2, limit = 250))
  56.  
  57. modl.param_list()
  58. #
  59. # area_NR = modl.integrate_pdf(pdf_name = "nonres")
  60. # area_jpsi = modl.integrate_pdf(pdf_name = "jpsi")
  61. # area_jpsi2 = modl.integrate_pdf(x_min = 2906, x_max = 3196, pdf_name = "jpsi", limit = 3000)
  62. # area_psi2s = modl.integrate_pdf(pdf_name = "psi2s")
  63. #
  64. # print(area_jpsi, area_jpsi2)
  65. #
  66. # print(area_jpsi/area_NR)
  67. # print(area_psi2s/area_NR)
  68.  
  69. # modl.param_val["jpsi scale"] *= jpsi_BR/NR_BR * area_NR/area_jpsi# /100
  70. # modl.param_val["psi2s scale"] *= psi2s_BR/NR_BR * area_NR/area_psi2s# /10
  71.  
  72. modl.mode = "no_data"
  73.  
  74. modl.draw_plots(part_set = 1, x_min = x_min, x_max = x_max, mode = "no_data")
  75.  
  76. modl.mode = mode
  77. # print(modl.total_scale_amp)
  78.  
  79. if load_set:
  80.  
  81. if mode == "true_data":
  82.  
  83. with open(r"./data/true_data/true_data_toy_{0}_range({1}-{2}).pkl".format(0, int(x_min), int(x_max)), "rb") as input_file:
  84. part_set = pkl.load(input_file)
  85.  
  86. elif mode == "slim_points":
  87.  
  88. with open(r"./data/slim_points/slim_points_toy_0_range({0}-{1}).pkl".format(int(x_min), int(x_max)), "rb") as input_file:
  89. part_set = pkl.load(input_file)
  90.  
  91. else:
  92. start = time.time()
  93. part_set = modl.generate_points(set_size, x_min = x_min, x_max = x_max, mode = mode, verbose = 1, nr_of_toys = nr_of_toys, nonres_set_size = nonres_set_size)
  94. print(time.time() - start)
  95. if draw:
  96. modl.draw_plots(part_set = part_set, x_min = x_min, x_max = x_max, mode = mode)
  97.  
  98. if fit:
  99. modl.fit_pdf_to_data(part_set)
  100.  
  101. # q = part_set["x_part"]
  102. #
  103. # q2 = tf.power(q, 2)
  104. #
  105. # jpsi_m = np.linspace(modl.param_val["jpsi phase"]-np.pi, modl.param_val["jpsi phase"]+np.pi, 100)
  106. #
  107. # neg_ll = []
  108. #
  109. # for mass in jpsi_m:
  110. # modl.param_val["jpsi phase"] = mass
  111. # neg_ll.append(modl.neg_log_likelihood(q2))
  112. # print(mass)
  113. #
  114. # plt.clf()
  115. #
  116. # plt.plot(jpsi_m, neg_ll, label = "jpsi phase scan")
  117. #
  118. # plt.savefig("test_ll.png")
  119.  
  120. print("Run finished")