diff --git a/__pycache__/pdg_const.cpython-37.pyc b/__pycache__/pdg_const.cpython-37.pyc index 648d836..518ff0f 100644 --- a/__pycache__/pdg_const.cpython-37.pyc +++ b/__pycache__/pdg_const.cpython-37.pyc Binary files differ diff --git a/__pycache__/raremodel.cpython-37.pyc b/__pycache__/raremodel.cpython-37.pyc index c33f983..4309e9d 100644 --- a/__pycache__/raremodel.cpython-37.pyc +++ b/__pycache__/raremodel.cpython-37.pyc Binary files differ diff --git a/pdg_const.py b/pdg_const.py index a553ebc..572a3db 100644 --- a/pdg_const.py +++ b/pdg_const.py @@ -59,9 +59,12 @@ "bplus" : [0.466, -0.885, -0.213], "bT" : [0.460, -1.089, -1.114], +"NR_BR": 4.37e-7, + #Resonances format(mass, width, phase, scale) -# "jpsi": (3096, 0.09, -1.5, 4.63977949519e-09), "jpsi": (3096, 0.09, -1.5, 2e-2), -# "psi2s": (3686, 0.3, -1.5, 7.28445380745e-10) -"psi2s": (3686, 0.3, -1.5, 3.14e-3) +"jpsi_BR": 6.02e-5, + +"psi2s": (3686, 0.3, -1.5, 3.14e-3), +"psi2s_BR": 4.97e-6 } diff --git a/plots/no_data/ff.png b/plots/no_data/ff.png index 4dc00f8..64a7f57 100644 --- a/plots/no_data/ff.png +++ b/plots/no_data/ff.png Binary files differ diff --git a/plots/no_data/pdf_and_parts.png b/plots/no_data/pdf_and_parts.png index c2acc0e..ad0bc7d 100644 --- a/plots/no_data/pdf_and_parts.png +++ b/plots/no_data/pdf_and_parts.png Binary files differ diff --git a/plots/no_data/vec_axiv.png b/plots/no_data/vec_axiv.png index c8daa9a..f95ed31 100644 --- a/plots/no_data/vec_axiv.png +++ b/plots/no_data/vec_axiv.png Binary files differ diff --git a/plots/points/ff.png b/plots/points/ff.png index 4dc00f8..64a7f57 100644 --- a/plots/points/ff.png +++ b/plots/points/ff.png Binary files differ diff --git a/plots/points/histo.png b/plots/points/histo.png index 15e218e..72a440f 100644 --- a/plots/points/histo.png +++ b/plots/points/histo.png Binary files differ diff --git a/plots/points/histo_raw.png b/plots/points/histo_raw.png index 1136109..aec7bbb 100644 --- a/plots/points/histo_raw.png +++ b/plots/points/histo_raw.png Binary files differ diff --git a/plots/points/pdf_and_parts.png b/plots/points/pdf_and_parts.png index c2acc0e..ad0bc7d 100644 --- a/plots/points/pdf_and_parts.png +++ b/plots/points/pdf_and_parts.png Binary files differ diff --git a/plots/points/vec_axiv.png b/plots/points/vec_axiv.png index c8daa9a..f95ed31 100644 --- a/plots/points/vec_axiv.png +++ b/plots/points/vec_axiv.png Binary files differ diff --git a/raremodel.py b/raremodel.py index 51b21d9..3edcb68 100644 --- a/raremodel.py +++ b/raremodel.py @@ -37,7 +37,7 @@ self.number_of_decays = pdg["number_of_decays"] - self.nr_of_part_list = {} + self.list_of_BR = {} self.resolution = 7.0 @@ -506,6 +506,8 @@ folder = mode elif mode == "fast_binned": folder = mode + elif mode == "no_data": + folder = mode else: folder = "points" @@ -585,6 +587,9 @@ plt.plot(dq, _pdf, label = name) _ += 1 + if x_min < 3096 < x_max: + plt.ylim(0, 4e-2) + plt.plot(dq, tot_y, label = "total pdf") prepare_plot("All pdfs") @@ -758,7 +763,7 @@ else: return com - def add_resonance(self, _mass, width, phase, scale, name, nr_of_part): #adds string to total_pdf and adds constant places for fit + def add_resonance(self, _mass, width, phase, scale, name, BR): #adds string to total_pdf and adds constant places for fit #Adds the resonace to the pdf in form of a string (to be executed later) @@ -766,7 +771,7 @@ self.param_val[name + " width"] = width self.param_val[name + " phase"] = phase self.param_val[name + " scale"] = scale - self.nr_of_part_list[name] = nr_of_part + self.list_of_BR[name] = BR self.param_val_orig = self.param_val @@ -821,7 +826,7 @@ else: return com - def add_cusp(self, mean, amp, sigma_L, sigma_R, nr_of_part): #adds string to total_pdf and adds constant places for fit + def add_cusp(self, mean, amp, sigma_L, sigma_R, BR): #adds string to total_pdf and adds constant places for fit name = "cusp" @@ -838,28 +843,33 @@ self.param_val[name + " sigma_L"] = sigma_L self.param_val[name + " sigma_R"] = sigma_R self.param_val[name + " scale"] = amp - self.nr_of_part_list[name] = nr_of_part + self.list_of_BR[name] = BR self.param_val_orig = self.param_val - def normalize_pdf(self, value = 1.0, verbose = 1): #Normalizes pdf with a global factor in front, saves mean and global factor + def integrate_pdf(self, pdf_name, x_min = None, x_max = None, limit = 3000): #Normalizes pdf with a global factor in front, saves mean and global factor - if verbose == 1: - print("Normalizing pdf...") + if not x_min: + x_min = self.x_min + if not x_max: + x_max = self.x_max - area, err = integrate.quad(self.total_pdf, self.x_min**2, self.x_max**2, limit = 3000) + if pdf_name == "total": + area, err = integrate.quad(self.total_pdf, x_min**2, x_max**2, limit = limit) - _mean = area/(self.x_max-self.x_min) + else: + _ = 0 + for pdf in self.total_pdf_list: + name = self.total_pdf_names[_] + if name == pdf_name: + area, err = integrate.quad(lambda q2: pdf(q2, param_val = self.param_val, absolut = True), x_min**2, x_max**2, limit = limit) + _ += 1 - self.param_val[0] = self.param_val[0]/area*value - self.param_val[1] = _mean/area*value + if err/area>0.05: + print("Big error on integration of {0}, increasing the limit may help!") - self.param_val_orig = self.param_val - - if verbose == 1: - print(" pdf normalized!") - print + return area def param_list(self, orig = True): #prints a list of all parameters of the form: number, name, value @@ -880,16 +890,15 @@ print("") - def add_nonres(self): + def add_nonres(self, BR): - nr_of_part = 44000 - name = "Nonres" + name = "nonres" def nonres_func(q2, param_val, absolut = False): return self.total_nonres(q2, parameters = param_val, name = name, absolut = absolut) self.total_pdf_list += [nonres_func] - # self.nr_of_part_listnr_of_part) + self.list_of_BR[name] = BR self.total_pdf_names += [name] def neg_log_likelihood(self, q2): @@ -899,84 +908,3 @@ ll = -1*np.sum(np.log(self.total_pdf(q2))) return ll - - # def fit_pdf_to_data(self, part_set, fitting_param_index = [2,3,4,5,9,13], unbinned_data = True, resolution = 0): - - # # x_part = part_set["x_part"] - # - # if resolution == 0: - # resolution = self.resolution - # - # if unbinned_data: - # - # #Load data - # - # x_part = part_set["x_part"] - # x_min = part_set["x_min"] - # x_max = part_set["x_max"] - # - # #Fill data into bins - # - # nbins = int((x_max-x_min)/resolution) - # - # bin_height , _x, _ = plt.hist(x_part, bins=nbins, range=(x_min, x_max)) - # plt.clf() - # - # bin_mean = [] - # - # for i in range(len(_x)-1): - # bin_mean.append((_x[i]+_x[i+1])/2) - # - # _part_set = {"bin_height": bin_height, "bin_mean": bin_mean, "_x": _x, "nbins": nbins, "x_max": x_max, "x_min": x_min} - # - # x_part = arr("f", x_part) - # - # nPoints = len(x_part) - # - # name = self.param_str - # - # npar = len(fitting_param_index) - # - # vstart = arr("d", npar*[0]) - # - # # print(len(vstart)) - # - # amp_err = self.param_val[0]/20 - # - # cusp_m_err = 10 - # - # step = arr( 'd', ( amp_err, cusp_m_err, 1.0, 0.1, 1e-7, self.param_val[9]/50, self.param_val[13] )) - # npar = len(vstart) - # - # def fcn(npar, deriv, f, apar, iflag): - # f[0] = self.log_likelihood(npar, apar, _part_set, resolution = 7.0, fitting_param_index = fitting_param_index) - # # return f[0] - # - # gMinuit = ROOT.TMinuit(npar) - # gMinuit.SetFCN( fcn ) - # - # gMinuit.SetErrorDef(0.5) - # - # arglist = arr( 'd', npar*[0.] ) - # ierflg = ROOT.Long(1998) - # - # arglist[0] = 1 - # gMinuit.mnexcm( "SET ERR", arglist, 1, ierflg ) - # - # _ = 0 - # - # for j in fitting_param_index: - # # print(_, fitting_param_index[_]) - # vstart[_] = self.param_val[fitting_param_index[_]] - # if _ == 1: - # vstart[_] -= 50 - # gMinuit.mnparm( _, self.param_str[fitting_param_index[_]], vstart[_], step[_], 0, 0, ierflg ) - # _ += 1 - # - # arglist[0] = 6000 - # arglist[1] = 1. - # - # gMinuit.mnexcm( "MIGRAD", arglist, 2, ierflg ) - # - # amin, edm, errdef = 0.18, 0.19, 0.20 - # nvpar, nparx, icstat = 1983, 1984, 1985 diff --git a/test.py b/test.py index 63871f8..528618b 100644 --- a/test.py +++ b/test.py @@ -12,11 +12,11 @@ modl = rm.model() -load_set = False +load_set = True draw = True fit = False -mode = "true_data" +mode = "slim_points" modl.mode = mode @@ -28,20 +28,21 @@ x_min = 3150.0 x_max= 3650.0 -#x_max = modl.x_max -#x_min = modl.x_min +x_max = modl.x_max +x_min = modl.x_min -modl.add_nonres() +NR_BR = pdg["NR_BR"] +modl.add_nonres(NR_BR) jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg["jpsi"] -jpsi_scale = jpsi_scale -modl.add_resonance(jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale, "jpsi", 5e6) +jpsi_BR = pdg["jpsi_BR"] +modl.add_resonance(jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale, "jpsi", jpsi_BR) psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"] -psi2s_scale = psi2s_scale -modl.add_resonance(psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale, "psi2s", 2e6) +psi2s_BR = pdg["psi2s_BR"] +modl.add_resonance(psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale, "psi2s", psi2s_BR) -modl.add_cusp(3550, 3e-7, 200, 7, 20000) +# modl.add_cusp(3550, 3e-7, 200, 7, 3e-7) # print(integrate.quad(modl.total_pdf, modl.x_min**2, modl.x_max**2, limit = 3000)) @@ -55,6 +56,19 @@ modl.param_list() +area_NR = modl.integrate_pdf(pdf_name = "nonres") +area_jpsi = modl.integrate_pdf(pdf_name = "jpsi") +area_jpsi2 = modl.integrate_pdf(x_min = 2906, x_max = 3196, pdf_name = "jpsi", limit = 300000) +area_psi2s = modl.integrate_pdf(pdf_name = "psi2s") + +print(area_jpsi, area_jpsi2) + +print(area_jpsi/area_NR) +print(area_psi2s/area_NR) + +# modl.param_val["jpsi scale"] *= jpsi_BR/NR_BR * area_NR/area_jpsi# /100 +# modl.param_val["psi2s scale"] *= psi2s_BR/NR_BR * area_NR/area_psi2s# /10 + modl.mode = "no_data" modl.draw_plots(part_set = 1, x_min = x_min, x_max = x_max, mode = "no_data") @@ -84,23 +98,23 @@ if fit: modl.fit_pdf_to_data(part_set) -q = part_set["x_part"] - -q2 = np.power(q, 2) - -jpsi_m = np.linspace(modl.param_val["jpsi phase"]-np.pi, modl.param_val["jpsi phase"]+np.pi, 100) - -neg_ll = [] - -for mass in jpsi_m: - modl.param_val["jpsi phase"] = mass - neg_ll.append(modl.neg_log_likelihood(q2)) - print(mass) - -plt.clf() - -plt.plot(jpsi_m, neg_ll, label = "jpsi phase scan") - -plt.savefig("test_ll.png") +# q = part_set["x_part"] +# +# q2 = np.power(q, 2) +# +# jpsi_m = np.linspace(modl.param_val["jpsi phase"]-np.pi, modl.param_val["jpsi phase"]+np.pi, 100) +# +# neg_ll = [] +# +# for mass in jpsi_m: +# modl.param_val["jpsi phase"] = mass +# neg_ll.append(modl.neg_log_likelihood(q2)) +# print(mass) +# +# plt.clf() +# +# plt.plot(jpsi_m, neg_ll, label = "jpsi phase scan") +# +# plt.savefig("test_ll.png") print("Run finished")