diff --git a/__pycache__/pdg_const.cpython-37.pyc b/__pycache__/pdg_const.cpython-37.pyc
index 518ff0f..00f3d75 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 4309e9d..000160f 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 572a3db..01fc430 100644
--- a/pdg_const.py
+++ b/pdg_const.py
@@ -62,9 +62,9 @@
 "NR_BR": 4.37e-7,
 
 #Resonances format(mass, width, phase, scale)
-"jpsi": (3096, 0.09, -1.5, 2e-2),
+"jpsi": (3096.0, 0.09, -1.5, 2e-2),
 "jpsi_BR": 6.02e-5,
 
-"psi2s": (3686, 0.3, -1.5, 3.14e-3),
+"psi2s": (3686.0, 0.3, -1.5, 3.14e-3),
 "psi2s_BR": 4.97e-6
 }
diff --git a/raremodel-tf.py b/raremodel-tf.py
new file mode 100644
index 0000000..9d75d5d
--- /dev/null
+++ b/raremodel-tf.py
@@ -0,0 +1,341 @@
+import numpy as np
+from pdg_const import pdg
+import matplotlib
+import matplotlib.pyplot as plt
+import pickle as pkl
+import sys
+import time
+from helperfunctions import display_time, prepare_plot
+import cmath as c
+import scipy.integrate as integrate
+from scipy.optimize import fminbound
+from array import array as arr
+import collections
+from itertools import compress
+import tensorflow as tf
+import zfit
+from zfit import ztf
+
+#Create pdfs
+
+def formfactor( q2, subscript): #returns real value
+    #check if subscript is viable
+
+    if subscript != "0" and subscript != "+" and subscript != "T":
+        raise ValueError('Wrong subscript entered, choose either 0, + or T')
+
+    #get constants
+
+    mK = ztf.constant(pdg['Ks_M'])
+    mbstar0 = ztf.constant(pdg["mbstar0"])
+    mbstar = ztf.constant(pdg["mbstar"])
+    b0 = ztf.constant(pdg["b0"])
+    bplus = ztf.constant(pdg["bplus"])
+    bT = ztf.constant(pdg["bT"])
+
+    mmu = ztf.constant(pdg['muon_M'])
+    mb = ztf.constant(pdg['bquark_M'])
+    ms = ztf.constant(pdg['squark_M'])
+    mB = ztf.constant(pdg['Bplus_M'])
+
+    #N comes from derivation in paper
+
+    N = 3
+
+    #some helperfunctions
+
+    tpos = (mB - mK)**2
+    tzero = (mB + mK)*(ztf.sqrt(mB)-ztf.sqrt(mK))**2
+
+    z_oben = ztf.sqrt(tpos - q2) - ztf.sqrt(tpos - tzero)
+    z_unten = ztf.sqrt(tpos - q2) + ztf.sqrt(tpos - tzero)
+    z = tf.divide(z_oben, z_unten)
+
+    #calculate f0
+
+    if subscript == "0":
+        prefactor = 1/(1 - q2/(mbstar0**2))
+        _sum = 0
+
+        for i in range(N):
+            _sum += b0[i]*(tf.pow(z,i))
+
+        return prefactor * _sum
+
+    #calculate f+ or fT
+
+    else:
+        prefactor = 1/(1 - q2/(mbstar**2))
+        _sum = 0
+
+        if subscript == "T":
+            b = bT
+        else:
+            b = bplus
+
+        for i in range(N):
+            _sum += b[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))
+
+        return prefactor * _sum
+
+class resonance(zfit.pdf.ZPDF):
+    _N_OBS = 1  # dimension, can be omitted
+    _PARAMS = ['mass', 'width', 'phase']  # the name of the parameters
+
+    def _unnormalized_pdf(self, x):
+
+        x = x.unstack_x() # returns a list with the columns: do x, y, z = ztf.unstack_x(x) for 3D
+        q2 = tf.pow(x, 2)
+
+        _mass = self.params['mass']
+        width = self.params['width']
+        phase = self.params['phase']
+
+        mmu = ztf.constant(pdg['muon_M'])
+
+        p = 0.5 * ztf.sqrt(q2 - 4*(mmu**2))
+
+        p0 =  0.5 * ztf.sqrt(_mass**2 - 4*mmu**2)
+
+        gamma_j = tf.divide(p, q2) * _mass * width / p0
+
+        #Calculate the resonance
+
+        _top = tf.complex(_mass * width, ztf.constant(0.0))
+
+        _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)
+
+        com = _top/_bottom
+
+        #Rotate by the phase
+
+        r = tf.abs(com)
+
+        _phase = tf.angle(com)
+
+        _phase += phase
+
+        x = tf.cos(phase)*r
+        y = tf.sin(phase)*r
+
+        com = tf.complex(x, y)
+
+        return com
+
+class bifur_gauss(zfit.pdf.ZPDF):
+    _N_OBS = 1  # dimension, can be omitted
+    _PARAMS = ['mean', 'amp', 'sigma_L', 'sigma_R']  # the name of the parameters
+
+    def _unnormalized_pdf(self, x):
+        x = ztf.unstack_x(x)  # returns a list with the columns: do x, y, z = ztf.unstack_x(x) for 3D
+
+        mean = self.params['mean']
+        amp = self.params['amp']
+        sigma_L = self.params['sigma_L']
+        sigma_R = self.params['sigma_R']
+
+        x_left = tf.where(x < mean, x, False)
+
+        x_right = tf.where(q >= mean, x, False)
+
+        #Calculate the exponential part of the cusp
+
+        _exp_left = ztf.exp(- tf.pow((x_left-mean),2) / (2 * sigma_L**2))
+        _exp_right = ztf.exp(- tf.pow((x_right-mean),2) / (2 * sigma_R**2))
+
+
+        #Scale so the total area under curve is 1 and the top of the cusp is continuous
+
+        _exp = _exp_left + _exp_right
+
+        dgamma = scale*_exp/(ztf.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R)
+
+        com = ztf.complex(dgamma, 0)
+
+        return com
+
+def axiv_nonres(q2, scale_axiv):
+
+    GF = ztf.constant(pdg['GF'])
+    alpha_ew = ztf.constant(pdg['alpha_ew'])
+    Vtb = ztf.constant(pdg['Vtb'])
+    Vts = ztf.constant(pdg['Vts'])
+    C10eff = ztf.constant(pdg['C10eff'])
+
+    mmu = ztf.constant(pdg['muon_M'])
+    mb = ztf.constant(pdg['bquark_M'])
+    ms = ztf.constant(pdg['squark_M'])
+    mK = ztf.constant(pdg['Ks_M'])
+    mB = ztf.constant(pdg['Bplus_M'])
+
+    q2 = tf.pow(x, 2)
+
+    #Some helperfunctions
+
+    beta = ztf.sqrt(tf.abs(1. - 4. * mmu**2. / q2))
+
+    kabs = ztf.sqrt(mB**2. +tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)
+
+    #prefactor in front of whole bracket
+
+    prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)
+
+    #left term in bracket
+
+    bracket_left = 2./3. * kabs**2. * beta**2. *tf.abs(C10eff*formfactor(q2, "+"))**2.
+
+    #middle term in bracket
+
+    _top = 4. * mmu**2. * (mB**2. - mK**2.) * (mB**2. - mK**2.)
+
+    _under = q2 * mB**2.
+
+    bracket_middle = _top/_under *tf.pow(tf.abs(C10eff * formfactor(q2, "0")), 2)
+
+    #Note sqrt(q2) comes from derivation as we use q2 and plot q
+
+    return ztf.to_complex(scale_axiv * prefactor1 * (bracket_left + bracket_middle) * 2 *ztf.sqrt(q2))
+
+def vec(q2, scale_vec):
+
+    scale = self.params['scale']
+
+    q2 = tf.pow(x, 2)
+
+    GF = ztf.constant(pdg['GF'])
+    alpha_ew = ztf.constant(pdg['alpha_ew'])
+    Vtb = ztf.constant(pdg['Vtb'])
+    Vts = ztf.constant(pdg['Vts'])
+    C7eff = ztf.constant(pdg['C7eff'])
+    C9eff = ztf.constant(pdg['C9eff'])
+
+    mmu = ztf.constant(pdg['muon_M'])
+    mb = ztf.constant(pdg['bquark_M'])
+    ms = ztf.constant(pdg['squark_M'])
+    mK = ztf.constant(pdg['Ks_M'])
+    mB = ztf.constant(pdg['Bplus_M'])
+
+    #Some helperfunctions
+
+    beta = ztf.sqrt(tf.abs(1. - 4. * mmu**2. / q2))
+
+    kabs = ztf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
+
+    #prefactor in front of whole bracket
+
+    prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)
+
+    #right term in bracket
+
+    prefactor2 = kabs**2 * (1. - 1./3. * beta**2)
+
+    abs_bracket = tf.abs(C9eff * formfactor(q2, "+") + 2 * C7eff * (mb + ms)/(mB + mK) * formfactor(q2, "T"))**2
+
+    bracket_right = prefactor2 * abs_bracket
+
+    #Note sqrt(q2) comes from derivation as we use q2 and plot q
+
+    return ztf.to_complex(scale_vec * prefactor1 * bracket_right * 2 * ztf.sqrt(q2))
+
+class total_func(zfit.func.ZFunc):
+    _N_OBS = 1  # dimension, can be omitted
+    _PARAMS = ['jpsi_mass', 'jpsi_scale', 'jpsi_phase', 'jpsi_width',
+                'psi2s_mass', 'psi2s_scale', 'psi2s_phase', 'psi2s_width',
+                'scale_vec', 'scale_axiv',
+                'cusp_mass', 'sigma_L', 'sigma_R', 'cusp_scale'
+                ]  # the name of the parameters
+
+    def _func(self, x):
+        tot = tf.abs(vec_nonres + axiv_nonres)
+
+        return a
+
+#Load data
+
+x_min = 2*pdg['muon_M']
+x_max = (pdg["Bplus_M"]-pdg["Ks_M"]-0.1)
+
+obs = zfit.Space('q', limits = (x_min, x_max))
+
+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:
+    part_set = pkl.load(input_file)
+
+x_part = part_set['x_part']
+
+x_part = x_part.astype('float64')
+
+data = zfit.data.Data.from_numpy(array=x_part, obs=obs)
+
+#Build pdf
+
+#Old false nonres
+
+nr_of_pdfs = 2
+
+scale_test = zfit.Parameter("scale_test", ztf.constant(1.0))
+
+frac1 = zfit.Parameter("frac1", 1/nr_of_pdfs)
+frac2 = zfit.Parameter("frac2", 1/nr_of_pdfs)
+
+axiv_nr = axiv_nonres(obs = obs)
+
+vec_nr = vec_nonres(obs = obs, scale = scale_test)
+
+#jpsi
+
+jpsi_mass, jpsi_width, jpsi_phase, jpsi_scale = pdg["jpsi"]
+
+jpsi_m = zfit.Parameter("jpsi_m", ztf.constant(jpsi_mass))
+jpsi_w = zfit.Parameter("jpsi_w", ztf.constant(jpsi_width))
+jpsi_p = zfit.Parameter("jpsi_p", ztf.constant(jpsi_phase))
+
+# jpsi_res = resonance(obs = obs, mass = jpsi_m, width = jpsi_w, phase = jpsi_p)
+
+#psi2s
+
+psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
+
+
+total_pdf = total_func.as_pdf()
+print(total_pdf.obs)
+
+nll = zfit.loss.UnbinnedNLL(model=total_pdf, data=data)
+
+minimizer = zfit.minimize.MinuitMinimizer()
+result = minimizer.minimize(nll, )
+
+param_errors = result.error()
+
+for var, errors in param_errors.items():
+    print('{}: ^{{+{}}}_{{{}}}'.format(var.name, errors['upper'], errors['lower']))
+
+print("Function minimum:", result.fmin)
+#_____________________________________________________________________
+
+#
+# class c9pdf(BaseFunctor):
+#     __init__(c9pdf, param1, param2, param3, obs, name="c9pdf"):
+#         params =
+#
+#   super()__init__(pdfs=c9pdf,...)
+#
+# def _unnormalized_pdf(..):
+#     c9pdf = self.pdfs[0]
+#     blabla... = bla  c9pdf.pdf(x)  ...
+#
+# #___________
+#
+# c9 = SumPDF(....)
+# pdf = PDF(c9pdf=c9,...)
+#
+# #___________
+#
+# def _unnormalized_pdf(...):
+#     a, b, c = x.unstack_x()
+#     pdf1d = self.pdfs[0]
+#     pdf1d.pdf(x)
+#
+# #___________
+#
+# c9pdf = self.pdfs[0]
+# c9eff = c9pdf.pdf(x)
diff --git a/raremodel.py b/raremodel.py
index 3edcb68..5fff237 100644
--- a/raremodel.py
+++ b/raremodel.py
@@ -12,6 +12,8 @@
 from array import array as arr
 import collections
 from itertools import compress
+import tensorflow as tf
+import zfit as zf
 
 mmu = pdg['muon_M']
 mb = pdg["bquark_M"]
@@ -90,11 +92,11 @@
         #some helperfunctions
 
         tpos = (self.mB - self.mK)**2
-        tzero = (self.mB + self.mK)*(np.sqrt(self.mB)-np.sqrt(self.mK))**2
+        tzero = (self.mB + self.mK)*(tf.sqrt(self.mB)-tf.sqrt(self.mK))**2
 
-        z_oben = np.sqrt(tpos - q2) - np.sqrt(tpos - tzero)
-        z_unten = np.sqrt(tpos - q2) + np.sqrt(tpos - tzero)
-        z = np.divide(z_oben, z_unten)
+        z_oben = tf.sqrt(tpos - q2) - tf.sqrt(tpos - tzero)
+        z_unten = tf.sqrt(tpos - q2) + tf.sqrt(tpos - tzero)
+        z = tf.divide(z_oben, z_unten)
 
         #calculate f0
 
@@ -103,7 +105,7 @@
             _sum = 0
 
             for i in range(N):
-                _sum += b0[i]*(np.power(z,i))
+                _sum += b0[i]*(tf.pow(z,i))
 
             return prefactor * _sum
 
@@ -119,7 +121,7 @@
                 b = bplus
 
             for i in range(N):
-                _sum += b[i] * (np.power(z, i) - ((-1)**(i-N)) * (i/N) * np.power(z, N))
+                _sum += b[i] * (tf.pow(z, i) - ((-1)**(i-N)) * (i/N) * tf.pow(z, N))
 
             return prefactor * _sum
 
@@ -141,17 +143,17 @@
 
         #Some helperfunctions
 
-        beta = np.sqrt(np.abs(1. - 4. * self.mmu**2. / q2))
+        beta = tf.sqrt(tf.abs(1. - 4. * self.mmu**2. / q2))
 
-        kabs = np.sqrt(mB**2. + np.power(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)
+        kabs = tf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2. * (mB**2. * mK**2. + mK**2. * q2 + mB**2. * q2) / mB**2.)
 
         #prefactor in front of whole bracket
 
-        prefactor1 = GF**2. *alpha_ew**2. * (np.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)
+        prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2. * kabs * beta / (128. * np.pi**5.)
 
         #left term in bracket
 
-        bracket_left = 2./3. * kabs**2. * beta**2. * np.abs(C10eff*self.formfactor(q2, "+"))**2.
+        bracket_left = 2./3. * kabs**2. * beta**2. * tf.abs(C10eff*self.formfactor(q2, "+"))**2.
 
         #middle term in bracket
 
@@ -159,11 +161,11 @@
 
         _under = q2 * mB**2.
 
-        bracket_middle = _top/_under * np.power(np.abs(C10eff * self.formfactor(q2, "0")), 2)
+        bracket_middle = _top/_under * tf.pow(tf.abs(C10eff * self.formfactor(q2, "0")), 2)
 
         #Note sqrt(q2) comes from derivation as we use q2 and plot q
 
-        return prefactor1 * (bracket_left + bracket_middle) * 2 * np.sqrt(q2)
+        return prefactor1 * (bracket_left + bracket_middle) * 2 * tf.sqrt(q2)
 
     def vec_nonres(self, q2): #returns real value
 
@@ -184,25 +186,25 @@
 
         #Some helperfunctions
 
-        beta = np.sqrt(abs(1. - 4. * self.mmu**2. / q2))
+        beta = tf.sqrt(abs(1. - 4. * self.mmu**2. / q2))
 
-        kabs = np.sqrt(mB**2. + np.power(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
+        kabs = tf.sqrt(mB**2. + tf.pow(q2, 2)/mB**2. + mK**4./mB**2. - 2 * (mB**2 * mK**2 + mK**2 * q2 + mB**2 * q2) / mB**2)
 
         #prefactor in front of whole bracket
 
-        prefactor1 = GF**2. *alpha_ew**2. * (np.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)
+        prefactor1 = GF**2. *alpha_ew**2. * (tf.abs(Vtb*Vts))**2 * kabs * beta / (128. * np.pi**5.)
 
         #right term in bracket
 
         prefactor2 = kabs**2 * (1. - 1./3. * beta**2)
 
-        abs_bracket = np.abs(C9eff * self.formfactor(q2, "+") + 2 * C7eff * (mb + ms)/(mB + mK) * self.formfactor(q2, "T"))**2
+        abs_bracket = tf.abs(C9eff * self.formfactor(q2, "+") + 2 * C7eff * (mb + ms)/(mB + mK) * self.formfactor(q2, "T"))**2
 
         bracket_right = prefactor2 * abs_bracket
 
         #Note sqrt(q2) comes from derivation as we use q2 and plot q
 
-        return prefactor1 * bracket_right * 2 * np.sqrt(q2)
+        return prefactor1 * bracket_right * 2 * tf.sqrt(q2)
 
     def total_nonres(self, q2, parameters = 0, name = "", absolut = False): #returns complex value
 
@@ -245,7 +247,7 @@
             return factor*(vec_nonres_part + axiv_nonres_part)
 
         else:
-            return factor * np.vectorize(complex)(vec_nonres_part + axiv_nonres_part, 0.0)
+            return factor * tf.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, nr_of_toys = 1): #returns part_set (=dic)
 
@@ -304,7 +306,7 @@
 
                     # print(y_part_raw[:40])
 
-                    choose = np.where(y_part_raw < self.total_pdf(np.power(x_part_raw, 2)), True, False)
+                    choose = np.where(y_part_raw < self.total_pdf(tf.pow(x_part_raw, 2)), True, False)
 
                     # print(choose[:40])
 
@@ -366,7 +368,7 @@
 
                     y_part_raw = np.random.uniform(low = 0.0, high = maxi, size = 2000000)
 
-                    choose_pdf_tot = np.where(y_part_raw < self.total_pdf(np.power(x_part_raw, 2)), True, False)
+                    choose_pdf_tot = np.where(y_part_raw < self.total_pdf(tf.pow(x_part_raw, 2)), True, False)
 
                     choose_x_min = np.where(x_min < x_part_raw, True, False)
 
@@ -374,7 +376,7 @@
 
                     choose = choose_pdf_tot * choose_x_max * choose_x_min
 
-                    choose_nonres = np.where(y_part_raw < self.total_nonres(np.power(x_part_raw, 2)), True, False)
+                    choose_nonres = np.where(y_part_raw < self.total_nonres(tf.pow(x_part_raw, 2)), True, False)
 
 
                     # print(choose[:40])
@@ -407,90 +409,6 @@
 
             return part_set
 
-        elif mode == "fast_binned":
-
-            #Calculate number of bins
-
-            nbins = int((x_max-x_min)/resolution)
-
-            for toy_nr in range(nr_of_toys):
-
-                print("Generating toy set {1} of {2} with {0} bins...".format(nbins, toy_nr, nr_of_toys))
-
-                #Get the mean dq values for each bin and get value of pdf there
-
-                dq = np.linspace(x_min, x_max ,nbins+1)
-
-                bin_mean = []
-                bin_true_height = []
-
-                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)
-
-                #Vary the nonres_set_size
-
-                _nonres_set_size = int(ROOT.TRandom1().Gaus(nonres_set_size, np.sqrt(nonres_set_size)))
-
-                #Scale true bin size according to nonres_set_size
-
-                _area, _area_err = integrate.quad(lambda x: self.total_nonres(x ,absolut = True), self.x_min**2, self.x_max**2, limit = 3000)
-
-                _mean_scan = _area/(self.x_max - self.x_min)
-
-                _mean_histo_nonres = _nonres_set_size/nbins
-
-                _ =  _mean_histo_nonres/_mean_scan
-
-                # print(_)
-
-                for i in range(len(bin_true_height)):
-                    bin_true_height[i] = bin_true_height[i]*_
-
-                #Start time
-
-                start = time.time()
-
-                bin_height = []
-
-                #Draw random numbers
-
-                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]))))
-
-                #progress calls
-                if verbose != 0:
-                    end = time.time()
-                    print(" Time to generate set: {0}s".format(end-start))
-
-                print(" {0} bins simulated".format(nbins))
-
-                print(" Set generated!")
-
-                #Save the set
-
-                if save:
-
-                    _ = 0
-
-                    for i in bin_height:
-                        _ += i
-
-                    part_set = {"bin_mean" : bin_mean, "bin_height": bin_height, "nbins": nbins, "bin_true_height": bin_true_height, "_nonres_set_size": nonres_set_size, "x_min": x_min, "x_max": x_max}
-
-                    pkl.dump( part_set, open("./data/fast_binned/binned_toy_{0}.pkl".format(toy_nr) , "wb" ) )
-
-                    print(" Set saved!")
-
-                print
-
-            return part_set
-
     def draw_plots(self, part_set, mode, x_min = 2* mmu, x_max = (mB - mK) - 0.1, resolution = 7.0): #draws and saves plots, returns nothing
         #Resolution based on experiment chosen to be ~7MeV
         #check for invalid entries
@@ -524,7 +442,7 @@
 
         #Translate to MeV**2
 
-        dq2 = np.power(dq, 2)
+        dq2 = tf.pow(dq, 2)
 
         #calculate formfactors
 
@@ -574,8 +492,6 @@
 
         psi2s_mass, psi2s_width, psi2s_phase, psi2s_scale = pdg["psi2s"]
 
-        total_y = np.sum([pdf(dq2, self.param_val) for pdf in self.total_pdf_list])
-
         tot_y = self.total_pdf(dq2)
 
         plt.clf()
@@ -722,7 +638,7 @@
         for pdf in self.total_pdf_list:
             _sum += pdf(q2, param_val = self.param_val)
 
-        return np.abs(_sum)
+        return tf.abs(_sum)
 
     def resonance(self, q2, parameters, name = "", absolut = False): #returns complex value
 
@@ -730,36 +646,36 @@
 
         #Helper variables
 
-        p = 0.5 * np.sqrt(q2 - 4*(mmu**2))
+        p = 0.5 * tf.sqrt(q2 - 4*(mmu**2))
 
-        p0 =  0.5 * np.sqrt(_mass**2 - 4*mmu**2)
+        p0 =  0.5 * tf.sqrt(_mass**2 - 4*mmu**2)
 
-        gamma_j = np.divide(p, q2) * _mass * width / p0
+        gamma_j = tf.divide(p, q2) * _mass * width / p0
 
         #Calculate the resonance
 
-        _top = np.complex(_mass * width, 0.0)
+        _top = tf.complex(_mass * width, 0.0)
 
 
-        _bottom = np.vectorize(complex)(_mass**2 - q2, -_mass*gamma_j)
+        _bottom = tf.complex(_mass**2 - q2, -_mass*gamma_j)
 
         com = _top/_bottom * scale
 
         #Rotate by the phase
 
-        r = np.abs(com)
+        r = tf.abs(com)
 
-        _phase = np.angle(com)
+        _phase = tf.angle(com)
 
         _phase += phase
 
-        x = np.cos(phase)*r
-        y = np.sin(phase)*r
+        x = tf.cos(phase)*r
+        y = tf.sin(phase)*r
 
-        com = np.vectorize(complex)(x, y)
+        com = tf.complex(x, y)
 
         if absolut:
-            return np.abs(com)
+            return tf.abs(com)
         else:
             return com
 
@@ -788,27 +704,27 @@
 
         #Get q out of q2
 
-        q = np.sqrt(q2)
+        q = tf.sqrt(q2)
 
         #Choose the right sigma depending on side of the cusp
 
-        q_left = np.where(q < mean, q, False)
+        q_left = tf.where(q < mean, q, False)
 
-        q_right = np.where(q >= mean, q, False)
+        q_right = tf.where(q >= mean, q, False)
 
         #Calculate the exponential part of the cusp
 
-        _exp_left = np.exp(- np.power((q_left-mean),2) / (2 * sigma_L**2))
-        _exp_right = np.exp(- np.power((q_right-mean),2) / (2 * sigma_R**2))
+        _exp_left = tf.exp(- tf.pow((q_left-mean),2) / (2 * sigma_L**2))
+        _exp_right = tf.exp(- tf.pow((q_right-mean),2) / (2 * sigma_R**2))
 
 
         #Scale so the total area under curve is 1 and the top of the cusp is continuous
 
         _exp = _exp_left + _exp_right
 
-        dgamma = scale*_exp/(np.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R)
+        dgamma = scale*_exp/(tf.sqrt(2*np.pi))*2*(sigma_L*sigma_R)/(sigma_L+sigma_R)
 
-        com = np.vectorize(complex)(dgamma, 0)
+        com = tf.complex(dgamma, 0)
 
 
         # _sum = 0
@@ -822,7 +738,7 @@
         # com *= factor
 
         if absolut:
-            return np.abs(com)
+            return tf.abs(com)
         else:
             return com
 
@@ -905,6 +821,14 @@
 
         # self.normalize_pdf(verbose = 0)
 
-        ll = -1*np.sum(np.log(self.total_pdf(q2)))
+        ll = -1*tf.reduce_sum(tf.log(self.total_pdf(q2)))
 
         return ll
+
+    def C9eff(self, q2):
+        C9eff = pdg["C9eff"]
+
+        for pdf in self.total_pdf_list:
+            C9eff += pdf(q2, param_val = self.param_val)
+
+        return C9eff
diff --git a/test.py b/test.py
index 528618b..ef4f1bc 100644
--- a/test.py
+++ b/test.py
@@ -55,16 +55,16 @@
 # print(integrate.quad(modl.total_pdf, modl.x_min**2, modl.x_max**2, limit = 250))
 
 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)
+#
+# 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 = 3000)
+# 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
@@ -100,7 +100,7 @@
 
 # q = part_set["x_part"]
 #
-# q2 = np.power(q, 2)
+# q2 = tf.power(q, 2)
 #
 # jpsi_m = np.linspace(modl.param_val["jpsi phase"]-np.pi, modl.param_val["jpsi phase"]+np.pi, 100)
 #