diff --git a/test.py b/test.py
index 35797c7..82e69f2 100644
--- a/test.py
+++ b/test.py
@@ -8,6 +8,7 @@
 import pickle as pkl
 import sys
 import time
+from helperfunctions import display_time
 
 mmu = pdg['muon_M']
 mb = pdg["bquark_M"]
@@ -15,26 +16,6 @@
 mK = pdg["Ks_M"]
 mB = pdg["Bplus_M"]
 
-intervals = (
-    ('w', 604800),  # 60 * 60 * 24 * 7
-    ('d', 86400),    # 60 * 60 * 24
-    ('h', 3600),    # 60 * 60
-    ('min', 60),
-    ('s', 1),
-    )
-
-def display_time(seconds, granularity=2):
-    result = []
-
-    for name, count in intervals:
-        value = seconds // count
-        if value:
-            seconds -= value * count
-            if value == 1:
-                name = name.rstrip('s')
-            result.append("{} {}".format(value, name))
-    return ', '.join(result[:granularity])
-
 class model:
     
     def __init__(self):
@@ -233,7 +214,7 @@
         
         #Range of function in MeV
         
-        dq = np.linspace(x_min, x_max ,1000)
+        dq = np.linspace(x_min, x_max ,5000)
         
         #Translate to MeV**2
         
@@ -249,6 +230,8 @@
         
         #Generate random points
         
+        psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
+        
         _max = max(dgamma_tot)
 
         x = []
@@ -277,7 +260,7 @@
             x.append(ROOT.TRandom1().Uniform(x_min, x_max))
             y.append(ROOT.TRandom1().Uniform(0, _max*1.05))
             
-            if y[-1] < self.total_nonres(x[-1]**2):
+            if y[-1] < self.total_pdf(x[-1]**2):
                 x_part.append(x[-1])
                 y_part.append(y[-1])
                 
@@ -509,27 +492,54 @@
     
     def total_pdf(self, q2):
         
+        #Calculate the pdf with the added resonances
+        
         exec("_sum = " + self.total_pdf_string)
             
         return _sum
     
     
-    def resonance(self, q2, _gamma, _mass, amp, phase):
+    def resonance(self, q2, _mass, width, phase, scale): #returns [real, imaginary]
         
-        __gamma = np.sqrt(_mass**2 * (_mass**2 + _gamma**2))
+        #calculate the resonance ---------------------------------------------> Formula correct?
         
-        k = 2 * np.sqrt(2) * _mass * _gamma * __gamma / (np.pi * np.sqrt(_mass**2 + __gamma))
+        #if np.abs(np.sqrt(q2) - _mass) < 300:
+            #return 0., 0.
         
-        reso = amp * k / ((q2 + _mass**2)**2 + _mass**2 * _gamma**2)
+        np.sqrt(mB**2 + q2**2/mB**2 + mK**4/mB**2 - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
         
-        #reso = 1
+        #print(q2)
         
-        return reso
+        p = 0.5 * np.sqrt(q2 - 4*(mmu**2))
+        
+        p0 =  0.5 * np.sqrt(_mass**2 - 4*mmu**2)
+        
+        gamma_j = p / p0 * _mass /q2 * width
+        
+        _top_im = - _mass**2 * width * gamma_j
+        
+        _top_re = _mass * width * (_mass**2 - q2)
+        
+        _bottom = (_mass**2 - q2) + _mass**2 * gamma_j**2
+        
+        real = _top_re/_bottom
+        
+        imaginary = _top_im/_bottom
+        
+        length = np.sqrt(real**2 + imaginary**2)
+        
+        real = np.cos(phase)*length
+        
+        imaginary = np.sin(phase)*length
+        
+        return [scale * real, scale * imaginary]
     
         
-    def add_resonance(self, _gamma, _mass, amp, phase):        
+    def add_resonance(self, _mass, width, phase, scale):
         
-        self.total_pdf_string += "+ self.resonance(q2,{0},{1},{2},{3})".format(_gamma, _mass, amp, phase)
+        #Adds the resonace to the pdf in form of a string (to be executed later)
+        
+        self.total_pdf_string += "+ self.resonance(q2,{0},{1},{2},{3})[0]".format(_mass, width, phase, scale)
         
         
         
@@ -543,41 +553,51 @@
 draw = True
 
 set_size = 1e5
-#modl.add_resonance(0.09, 3096.0, 0.02, -1.5)
+
+jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg["jpsi"]
+#modl.add_resonance(jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale)
+
+psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
+#modl.add_resonance(psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale)
 
 
-#modl.add_resonance(0.09, 3096.0, 1.02, -1.5)
+##modl.add_resonance(0.09, 3096.0, 1.02, -1.5)
 
-if load_set:
+#if load_set:
 
-    with open(r"set.pkl", "rb") as input_file:
-        set_dic = pkl.load(input_file)
+    #with open(r"set.pkl", "rb") as input_file:
+        #set_dic = pkl.load(input_file)
     
-    part_set = (set_dic["x_part"], set_dic["y_part"], set_dic["x"], set_dic["y"])
+    #part_set = (set_dic["x_part"], set_dic["y_part"], set_dic["x"], set_dic["y"])
 
-else:
-    part_set = modl.generate_points(set_size)
+#else:
+    #part_set = modl.generate_points(set_size)
 
-if draw:
-    modl.draw_plots(part_set = part_set)
-    
+#if draw:
+    #modl.draw_plots(part_set = part_set)
+
 #print(modl.total_pdf_string)
 
-#w = np.linspace(modl.x_min, modl.x_max, 1000)
+test_x = np.linspace(modl.x_min, modl.x_max, 1000)
 
-#test_y = []
+#print(test_x[1]**2 - modl.x_min**2)
 
-#for i in range(len(w)):
-    #exec("_ = modl.resonance(w[i], 0.09, 3096.0, 1.02, -1.5)")
-    #test_y.append(_)
-    ##w[i] = np.sqrt(w[i])
-    ##print(test_y[i])
+
+
+test_y = []
+
+for i in range(len(test_x)):
+    #print(i**2 - 4*(mmu**2))
+    test_y.append(modl.resonance(test_x[i]**2, jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale))
+    #resonance(self, q2, _mass, width, phase, scale):
+    #w[i] = np.sqrt(w[i])
+    #print(test_y[i])
                   
-#plt.clf()
+plt.clf()
 
-#plt.plot(w, test_y)
+plt.plot(test_x, test_y)
 
-#plt.savefig("test.png")
+plt.savefig("test.png")
     
 
 print("Run finished")