diff --git a/raremodel.py b/raremodel.py
index b74851d..2621771 100644
--- a/raremodel.py
+++ b/raremodel.py
@@ -47,6 +47,7 @@
 
         self.param_str = []
         self.param_val = []
+        self.param_val_orig = []
 
         #helperconstants
 
@@ -72,6 +73,8 @@
         self.param_str.append("cusp amp")
         self.param_val.append(0.0)
 
+        self.param_val_orig = self.param_val
+
     def formfactor(self, q2, subscript): #returns real value
 
         #check if subscript is viable
@@ -243,7 +246,7 @@
         else:
             return complex(vec_nonres_part + axiv_nonres_part, 0.0)
 
-    def generate_points(self, set_size = 10000, x_min = 2* mmu, x_max = (mB - mK) - 0.1, save = True, verbose = 1, mode = "true_data", resolution = 7.0, nonres_set_size = 44000): #returns part_set (=dic)
+    def generate_points(self, set_size = 10000, x_min = 2* mmu, x_max = (mB - mK) - 0.1, save = True, verbose = 1, mode = "true_data", resolution = 7.0, nonres_set_size = 44000, nr_of_toys = 1): #returns part_set (=dic)
 
         #Check if entered mode and verbose are valid
 
@@ -263,7 +266,7 @@
 
         #Prepare some variables used in all modes
 
-        _max = self.total_pdf(3096**2)
+        _max = self.total_pdf(3096.0**2)
 
         #slim_mode: Generates data of a set size in range(x_min,x_max) according to total_pdf -> Only saves the valid q values (Hit and miss method)
 
@@ -372,172 +375,210 @@
 
         elif mode == "true_data":
 
-            #Prepare variables
-
-            x_part = []
-
-            print("Generating set with {0} total rare nonresonant particles...".format(int(nonres_set_size)))
-
             #Change Seed here if necessary
 
             #ROOT.TRandom1().SetSeed(0)
 
             #Range of the whole nonres spectrum (used to count total hits in the nonres part)
 
+            if x_min < 3096.0 and x_max > 3096.0:
+                _max = self.total_pdf(3096.0**2)
+
+            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)
+
             x_min_nr = self.x_min
             x_max_nr = self.x_max
 
-            #Prepare verbose calls
+            _start = time.time()
+            _end = _start
 
-            if verbose != 0:
-                verbose_calls = []
-                j = 0
+            tot = []
+            for toy_nr in range(nr_of_toys):
+                tot.append(0)
+
+            for toy_nr in range(nr_of_toys):
+
+                #Prepare variables
+
                 o = 0
-                while j < 100:
-                    j += verbose
-                    verbose_calls.append(int(nonres_set_size*j/100))
+                x_part = []
+                _nonres_set_size = int(ROOT.TRandom1().Gaus(nonres_set_size, np.sqrt(nonres_set_size)))
 
-            #Take starting time and other variables for projected time
+                print("Generating toy set {1} of {2} with {0} total rare nonresonant particles...".format(int(_nonres_set_size), toy_nr, nr_of_toys - 1))
 
-            start = time.time()
+                #Take starting time and other variables for projected time
 
-            counter = 0
-            counter_x = 0
-            counter_x_nr = 0
+                start = time.time()
 
-            #Loop until enough particles in nonresonant part
+                counter = 0
+                counter_x = 0
+                counter_x_nr = 0
 
-            while counter_x_nr < nonres_set_size:
+                #Prepare verbose calls
 
-                #Generate random points
-
-                counter += 1
-                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
-
-                if x > x_min and x < x_max and y < self.total_pdf(x**2):
-                    x_part.append(x)
-                    counter_x += 1
-
-                #Count the hits in the nonres part -> Range of the whole pdf!
-
-                if y < abs(self.total_nonres(x**2)):
-                    counter_x_nr += 1
-
-                #progress calls
                 if verbose != 0:
+                    verbose_calls = []
+                    o = 0
+                    for j in range(100):
+                        verbose_calls.append(int(_nonres_set_size*(j+1)/100))
+
+                #Loop until enough particles in nonresonant part
+
+                while counter_x_nr < _nonres_set_size:
+
+                    #Generate random points
+
+                    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
+
+                    if x > x_min and x < x_max and y < self.total_pdf(x**2):
+                        x_part.append(x)
+                        counter_x += 1
+
+                    #Count the hits in the nonres part -> Range of the whole pdf!
+
+                    if x > x_min and x < x_max and y < self.total_nonres(x**2, absolut = True):
+                        counter_x_nr += 1
+
+                    #Calculate total estimated time to produce all toys
+
                     end = time.time()
-                    if o < 100 and counter%5000 == 0:
-                        print(" {0}{1} completed".format(o, "%"))
-                        print(" Projected time left: {0}".format(display_time(int((end-start)*nonres_set_size/(counter_x_nr+1)-(end-start)))))
-                        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 set: {0}".format(display_time(int(end-start))))
+                    tot[toy_nr] = (end-start)*_nonres_set_size/(counter_x_nr+1)
+                    total_estim = nr_of_toys*np.mean(tot[:toy_nr+1])
 
-                    if counter_x_nr == verbose_calls[o]:
-                        o += 1
+                    #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
 
-            print(" {0} out of {1} particles chosen!".format(int(counter_x), counter))
+                        counter += 1
 
-            print(" Set generated!")
+                _end = time.time()
 
-            #Save the set
+                print(" {0} out of {1} particles chosen!".format(int(counter_x), counter))
 
-            if save:
+                print(" Toy set {0} generated!".format(toy_nr))
 
-                part_set = {"x_part" : x_part, "counter_x_nr": counter_x_nr}
+                #Save the set
 
-                pkl.dump( part_set, open("./data/set_true_tot_nonres_{0}_range({1}-{2}).pkl".format(int(nonres_set_size),int(x_min), int(x_max)) , "wb" ) )
+                if save:
 
-                print(" Set saved!")
+                    part_set = {"x_part" : x_part, "counter_x_nr": counter_x_nr}
 
-            print
+                    pkl.dump( part_set, open("./data/true_data_toy_{0}_range({1}-{2}).pkl".format(_nonres_set_size ,int(x_min), int(x_max)) , "wb" ) )
+
+                    print(" Toy set {0} saved!".format(toy_nr))
+
+                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
 
             nbins = int((x_max-x_min)/resolution)
 
-            print("Generating set with {} bins...".format(nbins))
+            for toy_nr in range(nr_of_toys):
 
-            #Get the mean dq values for each bin and get value of pdf there
+                print("Generating toy set {1} of {2} with {0} bins...".format(nbins, toy_nr, nr_of_toys))
 
-            dq = np.linspace(x_min, x_max ,nbins+1)
+                #Get the mean dq values for each bin and get value of pdf there
 
-            bin_mean = []
-            bin_true_height = []
+                dq = np.linspace(x_min, x_max ,nbins+1)
 
-            for i in range(len(dq)-1):
-                a = dq[i]
-                b = dq[i+1]
-                bin_mean.append((a+b)/2)
-                _a, _e = integrate.quad(self.total_pdf, a**2, b**2, limit = 200)
-                height = _a/(b-a)
-                bin_true_height.append(height)
+                bin_mean = []
+                bin_true_height = []
 
-            #Scale true bin size according to nonres_set_size
+                for i in range(len(dq)-1):
+                    a = dq[i]
+                    b = dq[i+1]
+                    bin_mean.append((a+b)/2)
+                    _a, _e = integrate.quad(self.total_pdf, a**2, b**2, limit = 200)
+                    height = _a/(b-a)
+                    bin_true_height.append(height)
 
-            _area, _area_err = integrate.quad(lambda x: self.total_nonres(x ,absolut = True), self.x_min**2, self.x_max**2, limit = 3000)
+                #Vary the nonres_set_size
 
-            _mean_scan = _area/(self.x_max - self.x_min)
+                _nonres_set_size = int(ROOT.TRandom1().Gaus(nonres_set_size, np.sqrt(nonres_set_size)))
 
-            _mean_histo_nonres = nonres_set_size/nbins
+                #Scale true bin size according to nonres_set_size
 
-            _ =  _mean_histo_nonres/_mean_scan
+                _area, _area_err = integrate.quad(lambda x: self.total_nonres(x ,absolut = True), self.x_min**2, self.x_max**2, limit = 3000)
 
-            # print(_)
+                _mean_scan = _area/(self.x_max - self.x_min)
 
-            for i in range(len(bin_true_height)):
-                bin_true_height[i] = bin_true_height[i]*_
+                _mean_histo_nonres = _nonres_set_size/nbins
 
-            #Start time
+                _ =  _mean_histo_nonres/_mean_scan
 
-            start = time.time()
+                # print(_)
 
-            bin_height = []
+                for i in range(len(bin_true_height)):
+                    bin_true_height[i] = bin_true_height[i]*_
 
-            #Draw random numbers
+                #Start time
 
-            for i in range(len(bin_mean)):
-                bin_height.append(int(ROOT.TRandom1().Gaus(bin_true_height[i], np.sqrt(bin_true_height[i]))))
-                #print(int(ROOT.TRandom1().Gaus(bin_true_height[i], np.sqrt(bin_true_height[i]))))
+                start = time.time()
 
-            #progress calls
-            if verbose != 0:
-                end = time.time()
-                print(" Time to generate set: {0}s".format(end-start))
+                bin_height = []
 
-            print(" {0} bins simulated".format(nbins))
+                #Draw random numbers
 
-            print(" Set generated!")
+                for i in range(len(bin_mean)):
+                    bin_height.append(int(ROOT.TRandom1().Gaus(bin_true_height[i], np.sqrt(bin_true_height[i]))))
+                    #print(int(ROOT.TRandom1().Gaus(bin_true_height[i], np.sqrt(bin_true_height[i]))))
 
-            #Save the set
+                #progress calls
+                if verbose != 0:
+                    end = time.time()
+                    print(" Time to generate set: {0}s".format(end-start))
 
-            if save:
+                print(" {0} bins simulated".format(nbins))
 
-                _ = 0
+                print(" Set generated!")
 
-                for i in bin_height:
-                    _ += i
+                #Save the set
 
-                part_set = {"bin_mean" : bin_mean, "bin_height": bin_height, "nbins": nbins, "bin_true_height": bin_true_height}
+                if save:
 
-                pkl.dump( part_set, open("./data/binned_set_{0}bins_{1}part.pkl".format(nbins, _) , "wb" ) )
+                    _ = 0
 
-                print(" Set saved!")
+                    for i in bin_height:
+                        _ += i
 
-            print
+                    part_set = {"bin_mean" : bin_mean, "bin_height": bin_height, "nbins": nbins, "bin_true_height": bin_true_height, "nonres_set_size": nonres_set_size}
+
+                    pkl.dump( part_set, open("./data/fast_binned/binned_toy_{0}.pkl".format(toy_nr) , "wb" ) )
+
+                    print(" Set saved!")
+
+                print
 
             return part_set
 
@@ -681,7 +722,7 @@
             bin_mean = part_set["bin_mean"]
             bin_height = part_set["bin_height"]
 
-            #Scale the bins to the pdf
+            #Scale the pdf to the histo
 
             _sum = np.sum(bin_height)
 
@@ -689,15 +730,15 @@
             #
             nbins = part_set["nbins"]
             #
-            # _bin_mean = np.mean(bin_height)
-            # pdf_mean = self.param_val[1]
-            #
-            # for i in range(len(bin_height)):
-            #     bin_height[i] = bin_height[i] * pdf_mean/_bin_mean
+            _ = self.param_val[0]
+
+            area, err = integrate.quad(self.total_pdf, x_min**2, x_max**2, limit = 300)
+            _m = area/(x_max-x_min)
+            self.param_val[0] = self.param_val[0]*_mean/_m
 
             plt.clf()
 
-            plt.ylim(0.0,_mean*2)
+            plt.ylim(0.0, 2*_mean)
 
             #Choose range in between the 2 resonances
             # plt.xlim(3150, 3650)
@@ -712,6 +753,8 @@
 
             plt.savefig("./plots/fast_binned/histo.png")
 
+            self.param_val[0] = _
+
             print(" histo.png created")
 
             print(" All plots drawn \n")
@@ -801,6 +844,8 @@
 
             _mean_histo = float(np.mean(_y))
 
+            plt.ylim(0.0, 2*_mean_histo)
+
             prepare_plot("Binned toy data")
 
             plt.savefig("./plots/points/histo_raw.png")
@@ -814,19 +859,28 @@
 
             #Plot histo
 
-            #Scale the histo to the pdf
+            #Scale the pdf to the histo
 
             _ = []
 
+            __ = self.param_val[0]
+
             for i in range(len(_y)):
-                _y[i] = _y[i]*np.mean(dgamma_tot)/_mean_histo
                 _.append(x_min+(x_max-x_min)/bins*i)
 
 
+            area, err = integrate.quad(self.total_pdf, x_min**2, x_max**2, limit = 300)
+            _mean = area/(x_max-x_min)
+            self.param_val[0] = self.param_val[0]*_mean_histo/_mean
+
+            plt.ylim(0.0, 2*_mean_histo)
+
             plt.hist(_, bins=bins, range=(x_min, x_max), weights = _y, label = "toy data binned  (scaled down to pdf)")
 
             plt.plot(dq, dgamma_tot, label = "pdf")
 
+            self.param_val[0] = __
+
             #Only show from 0 to twice the mean of pdf
             # plt.ylim(0, 2*self.param_val[1])
 
@@ -886,7 +940,7 @@
 
         return _sum
 
-    def resonance(self, q2, _mass, width, phase, scale): #returns complex value
+    def resonance(self, q2, _mass, width, phase, scale, absolut = False): #returns complex value
 
         #Helper variables
 
@@ -917,7 +971,10 @@
 
         com = complex(x,y)
 
-        return self.param_val[0]*com
+        if absolut:
+            return abs(self.param_val[0]*com)
+        else:
+            return self.param_val[0]*com
 
     def add_resonance(self, _mass, width, phase, scale, name): #adds string to total_pdf and adds constant places for fit
 
@@ -936,11 +993,13 @@
         self.param_str.append(name + " scale")
         self.param_val.append(scale)
 
+        self.param_val_orig = self.param_val
+
         self.total_pdf_string += "+ self.resonance(q2,{0},{1},{2},{3})".format(self.param_val[6+_*4], self.param_val[7+_*4],self.param_val[8+_*4], self.param_val[9+_*4])
 
         self.res_counter += 1
 
-    def bifur_gauss(self, q2, mean, amp, sigma_L, sigma_R): #returns complex value
+    def bifur_gauss(self, q2, mean, amp, sigma_L, sigma_R, absolut = False): #returns complex value
 
         #Get q out of q2
 
@@ -963,7 +1022,10 @@
 
         com = complex(dgamma, 0)
 
-        return self.param_val[0]*com
+        if absolut:
+            return abs(self.param_val[0]*com)
+        else:
+            return self.param_val[0]*com
 
     def add_cusp(self, mean, amp, sigma_L, sigma_R): #adds string to total_pdf and adds constant places for fit
 
@@ -978,7 +1040,9 @@
         self.param_val[4] = sigma_R
         self.param_val[5] = amp
 
-    def normalize_pdf(self): #Normalizes pdf with a global factor in front, saves mean and global factor
+        self.param_val_orig = self.param_val
+
+    def normalize_pdf(self, value = 1.0): #Normalizes pdf with a global factor in front, saves mean and global factor
 
         print("Normalizing pdf...")
 
@@ -986,21 +1050,29 @@
 
         _mean = area/(self.x_max-self.x_min)
 
-        self.param_val[0] = self.param_val[0]/area
+        self.param_val[0] = self.param_val[0]/area*value
 
-        self.param_val[1] = _mean
+        self.param_val[1] = _mean/area*value
+
+        self.param_val_orig = self.param_val
 
         print(" pdf normalized!")
 
         print
 
-    def param_list(self): #prints a list of all parameters of the form: number, name, value
+    def param_list(self, orig = True): #prints a list of all parameters of the form: number, name, value
 
-        print("List of parameters")
+        if orig:
+            print("List of parameters (true values)")
+        else:
+            print("List of parameters (fitted values)")
         print(" Nr. Name Value")
 
-        for i in range(len(self.param_val)):
-            print(" {0}. {1}: {2}".format(i, self.param_str[i], self.param_val[i]))
+        for i in range(len(self.param_val_orig)):
+            if orig:
+                print(" {0}. {1}: {2}".format(i, self.param_str[i], self.param_val_orig[i]))
+            else:
+                print(" {0}. {1}: {2}".format(i, self.param_str[i], self.param_val[i]))
 
         print
 
@@ -1008,13 +1080,145 @@
 
         self.total_pdf_string += " + self.total_nonres(q2)"
 
+    def log_likelihood(self, npar, par, part_set, unbinned_data = True, resolution = 7.0): #Replaced soon with TMinuit
 
-    def log_likelihood(self, x_part, cusp_amp): #Replaced soon with TMinuit
+        self.param_val = par[1,:]
 
-        _sum = 0.0
+        if unbinned_data:
 
-        for i in x_part:
+            #Load data
 
-            _sum += np.log(self.total_pdf(i**2, cusp_amp = cusp_amp, scan = True))
+            x_part = part_set["x_part"]
+            x_min = part_set["x_min"]
+            x_max = part_set["x_max"]
 
-        return _sum
+            #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)
+
+        ll = 0.0
+
+        for i in bin_height:
+            area, _ = integrate.quad(self.total_pdf, self.x_min**2, self.x_max**2, limit = 300)
+            y_true = area/(_x[i+1]-_x[i])
+
+            ll += np.log(bin_height[i]*y_true)
+
+        ll = -ll
+
+        return ll
+
+    def fit_pdf_to_data(self, part_set, par, err):
+
+        x_part = part_set["x_part"]
+
+        x_part = arr("f", x_part)
+
+        nPoints = len(x_part)
+
+        name = self.param_str
+
+        vstart = arr("d", self.param_val)
+
+        step   = arr( 'd', ( amp_err, 0, cusp_m_err, 1.0, 0.1, 1e-7, 0, 0, 0, self.param_val[9]/50, 0, 0, 0, self.param_val[13] ))
+        npar = len(vstart)
+
+        def fcn(npar, f, x_part, iflag):
+            self.log_likelihood(npar, par, err, part_set, unbinned_data = True, resolution = 7.0)
+
+        #__________________________________________________________________________________
+
+        gMinuit = ROOT.TMinuit(13)
+        gMinuit.SetFCN( self.log_likelihood(part_set = part_set, ) )
+
+        arglist = arr( 'd', 15*[0.] )
+        gMinuit.mnexcm( "SET ERR", arglist, 1, ierflg )
+
+        ierflg = 0
+
+        arglist = arr( 'd', (6000, 0.3) )
+
+        amp_err = np.sqrt(self.param_val[0])
+        cusp_m_err = np.sqrt(self.param_val[2])
+
+        vstart = arr( 'd', self.param_val )
+        step   = arr( 'd', ( amp_err, 0, cusp_m_err, 1.0, 0.1, 1e-7, 0, 0, 0, self.param_val[9]/50, 0, 0, 0, self.param_val[13] ))
+
+        for i in range(len(self.param_val)):
+            gMinuit.mnparm( i, "a{0}".format(i), vstart[i], step[i], 0, 0, ierflg )
+
+        gMinuit.mnexcm( "MIGRAD", arglist, 2, ierflg )
+
+        #__________________________________________________________________________________
+
+        Error = 0;
+        z = array( 'f', ( 1., 0.96, 0.89, 0.85, 0.78 ) )
+        errorz = array( 'f', 5*[0.01] )
+
+        x = array( 'f', ( 1.5751, 1.5825,  1.6069,  1.6339,   1.6706  ) )
+        y = array( 'f', ( 1.0642, 0.97685, 1.13168, 1.128654, 1.44016 ) )
+
+        ncount = 0
+
+        ##______________________________________________________________________________
+        def testfit():
+           #
+           # gMinuit = TMinuit(5)
+           # gMinuit.SetFCN( fcn )
+           #
+           # arglist = array( 'd', 10*[0.] )
+           # ierflg = 1982
+           #
+           # arglist[0] = 1
+           # gMinuit.mnexcm( "SET ERR", arglist, 1, ierflg )
+
+         # Set starting values and step sizes for parameters
+           vstart = array( 'd', ( 3,  1,  0.1,  0.01  ) )
+           step   = array( 'd', ( 0.1, 0.1, 0.01, 0.001 ) )
+           gMinuit.mnparm( 0, "a1", vstart[0], step[0], 0, 0, ierflg )
+           gMinuit.mnparm( 1, "a2", vstart[1], step[1], 0, 0, ierflg )
+           gMinuit.mnparm( 2, "a3", vstart[2], step[2], 0, 0, ierflg )
+           gMinuit.mnparm( 3, "a4", vstart[3], step[3], 0, 0, ierflg )
+
+         # Now ready for minimization step
+           arglist[0] = 500
+           arglist[1] = 1.
+           gMinuit.mnexcm( "MIGRAD", arglist, 2, ierflg )
+
+         # Print results
+           amin, edm, errdef = 0.18, 0.19, 0.20
+           nvpar, nparx, icstat = 1983, 1984, 1985
+           gMinuit.mnstat( amin, edm, errdef, nvpar, nparx, icstat )
+           gMinuit.mnprin( 3, amin )
+
+
+        ##______________________________________________________________________________
+        # def fcn( npar, gin, f, par, iflag ):
+           global ncount
+         #   nbins = 5
+         #
+         # # calculate chisquare
+         #   chisq, delta = 0., 0.
+         #   for i in range(nbins):
+         #      delta  = (z[i]-func(x[i],y[i],par))/errorz[i]
+         #      chisq += delta*delta
+         #
+         #   f[0] = chisq
+           ncount += 1
+        #
+        # def func( x, y, par ):
+        #    value = ( (par[0]*par[0])/(x*x)-1)/ ( par[1]+par[2]*y-par[3]*y*y)
+        #    return value
+
+
+        ##______________________________________________________________________________
+        if __name__ == '__main__':
+           testfit()
diff --git a/raremodel.pyc b/raremodel.pyc
index 0571b9a..bad569f 100644
--- a/raremodel.pyc
+++ b/raremodel.pyc
Binary files differ
diff --git a/test.py b/test.py
index d7b80be..31bdf6f 100644
--- a/test.py
+++ b/test.py
@@ -21,19 +21,21 @@
 draw = True
 step_scan = False
 
-mode = "fast_binned"
+mode = "true_data"
 
 modl.mode = mode
 
 set_size = 1e5
 nonres_set_size = 44000
-modl.reco_steps = 100000
+# nonres_set_size = 1000
+modl.reco_steps = 10000
+nr_of_toys = 10
 
-# x_min = 3150.0
-# x_max= 3650.0
-
-x_min = modl.x_min
-x_max = modl.x_max
+x_min = 3150.0
+x_max= 3650.0
+#
+# x_min = modl.x_min
+# x_max = modl.x_max
 
 modl.add_nonres()
 
@@ -73,7 +75,7 @@
         counter = set_dic["counter_tot"]
 
 else:
-    part_set = modl.generate_points(set_size, x_min = x_min, x_max = x_max, mode = mode, verbose = 1)
+    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)
 
 if draw:
     modl.draw_plots(part_set = part_set, x_min = x_min, x_max = x_max, mode = mode)