diff --git a/__pycache__/raremodel.cpython-37.pyc b/__pycache__/raremodel.cpython-37.pyc index a7263b4..c33f983 100644 --- a/__pycache__/raremodel.cpython-37.pyc +++ b/__pycache__/raremodel.cpython-37.pyc Binary files differ diff --git "a/data/slim_points/slim_points_toy_0_range\050211-4781\051.pkl" "b/data/slim_points/slim_points_toy_0_range\050211-4781\051.pkl" index e0d7e05..c7b303e 100644 --- "a/data/slim_points/slim_points_toy_0_range\050211-4781\051.pkl" +++ "b/data/slim_points/slim_points_toy_0_range\050211-4781\051.pkl" Binary files differ diff --git "a/data/slim_points/slim_points_toy_0_range\0503150-3650\051.pkl" "b/data/slim_points/slim_points_toy_0_range\0503150-3650\051.pkl" index e347ff3..2b1cb5c 100644 --- "a/data/slim_points/slim_points_toy_0_range\0503150-3650\051.pkl" +++ "b/data/slim_points/slim_points_toy_0_range\0503150-3650\051.pkl" Binary files differ diff --git a/plots/points/ff.png b/plots/points/ff.png index 64a7f57..4dc00f8 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 8830c35..15e218e 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 68a3044..1136109 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 da61548..c2acc0e 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 f95ed31..c8daa9a 100644 --- a/plots/points/vec_axiv.png +++ b/plots/points/vec_axiv.png Binary files differ diff --git a/plots/true_data/ff.png b/plots/true_data/ff.png new file mode 100644 index 0000000..4dc00f8 --- /dev/null +++ b/plots/true_data/ff.png Binary files differ diff --git a/plots/true_data/pdf_and_parts.png b/plots/true_data/pdf_and_parts.png new file mode 100644 index 0000000..c2acc0e --- /dev/null +++ b/plots/true_data/pdf_and_parts.png Binary files differ diff --git a/plots/true_data/vec_axiv.png b/plots/true_data/vec_axiv.png new file mode 100644 index 0000000..c8daa9a --- /dev/null +++ b/plots/true_data/vec_axiv.png Binary files differ diff --git a/raremodel.py b/raremodel.py index b459a91..51b21d9 100644 --- a/raremodel.py +++ b/raremodel.py @@ -281,7 +281,7 @@ maxi = self.total_pdf(np.array(min_x)) - if np.abs(min_x-3096) < 10 or np.abs(min_x-3686): + if np.abs(min_x-3096) < 10 or np.abs(min_x-3686) < 10: maxi /= 1000 for toy in range(nr_of_toys): @@ -338,131 +338,75 @@ elif mode == "true_data": - #Change Seed here if necessary + set_sizes = np.random.poisson(nonres_set_size, nr_of_toys) - #ROOT.TRandom1().SetSeed(0) + def total_pdf_neg(q2): + return -1*self.total_pdf(q2) - #Range of the whole nonres spectrum (used to count total hits in the nonres part) + min_x = fminbound(total_pdf_neg, x_min**2, x_max**2) - non_res_sizes = np.random.normal(nonres_set_size, np.sqrt(nonres_set_size), nr_of_toys) + maxi = self.total_pdf(np.array(min_x)) - if x_min < 3096.0 and x_max > 3096.0: - _max = self.total_pdf(3096.0**2) + if np.abs(min_x-3096) < 10 or np.abs(min_x-3686) < 10: + maxi /= 1000 - elif self.total_pdf(x_min**2) < self.total_pdf(x_max**2): - _max = self.total_pdf(x_max**2) - else: - _max = self.total_pdf(x_min**2) + for toy in range(nr_of_toys): - x_min_nr = self.x_min - x_max_nr = self.x_max + set_size_intermed = 0 - _start = time.time() - _end = _start + print("Generating true sized toy {0} of nonres size {1}...".format(toy+1, int(set_sizes[toy]))) - tot = [] - for toy_nr in range(nr_of_toys): - tot.append(0) - - for toy_nr in range(nr_of_toys): - - #Prepare variables - - o = 0 x_part = [] - _nonres_set_size = int(ROOT.TRandom1().Gaus(nonres_set_size, np.sqrt(nonres_set_size))) - print("Generating toy set {1} of {2} with {0} total rare nonresonant particles...".format(int(_nonres_set_size), toy_nr, nr_of_toys - 1)) + _ = 0 - #Take starting time and other variables for projected time + while set_size_intermed < set_sizes[toy]: - start = time.time() + x_part_raw = np.random.uniform(low = self.x_min, high = self.x_max, size = 2000000) - counter = 0 - counter_x = 0 - counter_x_nr = 0 + y_part_raw = np.random.uniform(low = 0.0, high = maxi, size = 2000000) - #Prepare verbose calls + choose_pdf_tot = np.where(y_part_raw < self.total_pdf(np.power(x_part_raw, 2)), True, False) - if verbose != 0: - verbose_calls = [] - o = 0 - for j in range(100): - verbose_calls.append(int(_nonres_set_size*(j+1)/100)) + choose_x_min = np.where(x_min < x_part_raw, True, False) - #Loop until enough particles in nonresonant part + choose_x_max = np.where(x_max > x_part_raw, True, False) - while counter_x_nr < _nonres_set_size: + choose = choose_pdf_tot * choose_x_max * choose_x_min - #Generate random points + choose_nonres = np.where(y_part_raw < self.total_nonres(np.power(x_part_raw, 2)), True, False) - x = ROOT.TRandom1().Uniform(x_min_nr, x_max_nr) - y = ROOT.TRandom1().Uniform(0, _max) - #Add the part to the set if within range(x_min,x_max) and under pdf + # print(choose[:40]) - if x > x_min and x < x_max and y < self.total_pdf(x**2): - x_part.append(x) - counter_x += 1 + x_part_intermed = list(compress(x_part_raw, choose)) - #Count the hits in the nonres part -> Range of the whole pdf! + x_part += x_part_intermed - if y < self.total_nonres(x**2, absolut = True): - counter_x_nr += 1 + set_size_intermed += np.sum(choose_nonres) - #Calculate total estimated time to produce all toys + x_part = x_part[:int(set_sizes[toy])] - end = time.time() - tot[toy_nr] = (end-start)*_nonres_set_size/(counter_x_nr+1) - total_estim = nr_of_toys*np.mean(tot[:toy_nr+1]) + x_part = np.array(x_part) - #progress calls - if verbose != 0: - if o < 100 and counter%10000 == 0: - print(" {0}{1} completed of this toy".format(o, "%")) - print(" {0} rare nonresonant particles".format(counter_x_nr)) - print(" {0} total particles".format(counter_x)) - print(" Projected time left for this toy: {0}".format(display_time(int((end-start)*_nonres_set_size/(counter_x_nr+1)-(end-start))))) - print(" Projected time left for all toys: {0}".format(display_time(int( total_estim - (end-_start) )))) - sys.stdout.write("\033[F") - sys.stdout.write("\x1b[2K") - sys.stdout.write("\033[F") - sys.stdout.write("\x1b[2K") - sys.stdout.write("\033[F") - sys.stdout.write("\x1b[2K") - sys.stdout.write("\033[F") - sys.stdout.write("\x1b[2K") - sys.stdout.write("\033[F") - sys.stdout.write("\x1b[2K") - if o >=100: - print(" Time to generate this toy set: {0}".format(display_time(int(end-start)))) - - if counter_x_nr == verbose_calls[o]: - o += 1 - - counter += 1 - - _end = time.time() - - print(" {0} out of {1} particles chosen!".format(int(counter_x), counter)) - - print(" Toy set {0} generated!".format(toy_nr)) + print(" Toy {0} of {1} generated!".format(toy+1, nr_of_toys)) #Save the set if save: - part_set = {"x_part" : x_part, "counter_x_nr": counter_x_nr, "x_min": x_min, "x_max": x_max, "mode": mode} + part_set = {"x_part" : x_part} - pkl.dump( part_set, open("./data/true_data/true_data_toy_{0}_range({1}-{2}).pkl".format(toy_nr ,int(x_min), int(x_max)) , "wb" ) ) + pkl.dump( part_set, open("./data/slim_points/slim_points_toy_{0}_range({1}-{2}).pkl".format(toy, int(x_min), int(x_max)) , "wb" ) ) - print(" Toy set {0} saved!".format(toy_nr)) + print(" Set saved!\n") print #returns all the chosen points (x_part, y_part) and all the random points generated (x, y) return part_set + elif mode == "fast_binned": #Calculate number of bins @@ -948,7 +892,7 @@ # self.nr_of_part_listnr_of_part) self.total_pdf_names += [name] - def log_likelihood(self, q2): + def neg_log_likelihood(self, q2): # self.normalize_pdf(verbose = 0) diff --git a/test.py b/test.py index b1a4139..63871f8 100644 --- a/test.py +++ b/test.py @@ -16,7 +16,7 @@ draw = True fit = False -mode = "slim_points" +mode = "true_data" modl.mode = mode @@ -25,11 +25,11 @@ # nonres_set_size = 1000 nr_of_toys = 1 -# x_min = 3150.0 -# x_max= 3650.0 +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() @@ -84,6 +84,23 @@ if fit: modl.fit_pdf_to_data(part_set) -print() +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") diff --git a/test_ll.png b/test_ll.png new file mode 100644 index 0000000..da5f544 --- /dev/null +++ b/test_ll.png Binary files differ