diff --git a/scripts/checkTrackSelection.py b/scripts/checkTrackSelection.py
index 4d5390b..4149ab9 100644
--- a/scripts/checkTrackSelection.py
+++ b/scripts/checkTrackSelection.py
@@ -41,6 +41,7 @@
 # r.gStyle.SetPalette(r.kBird)
 # r.gStyle.SetPalette(r.kCopper)
 r.gStyle.SetPalette(r.kViridis)
+r.gStyle.SetNumberContours(99)
 
 
 def sqrt_of_2d_histogram(h):
@@ -248,12 +249,18 @@
 
 
 class MyFitFunction():
-    def __init__(self, func):
-        # self.n_var = n_var
+    def __init__(self, func, n_var):
+        self.n_var = n_var
         self.func = func
 
     def wrap_func(self, x, par):
-        return self.func(*x, par=par)
+        args = []
+        for i in range(self.n_var):
+            args.append(x[i])
+        args = tuple(args)
+        # print(x[0], x[1], x[2], x[3], par)
+        # return self.func(*x, par=par)
+        return self.func(*args, par=par)
 
     def __call__(self, x, par):
         return self.wrap_func(x, par)
@@ -261,7 +268,7 @@
 
 def myTFunc(name, classname, func, func_range, par_names, par_inits=[], par_limits=[]):
     n_pars = len(par_names)
-    f = r.__getattr__(classname)(name, MyFitFunction(func), *(func_range + (n_pars,)))
+    f = r.__getattr__(classname)(name, MyFitFunction(func, int(classname[2])), *(func_range + (n_pars,)))
     f.SetParNames(*par_names)
     if par_inits:
         f.SetParameters(*par_inits)
@@ -299,12 +306,16 @@
     punzi = h_list[-1]
     conf = mydict({
         'names': ('norm', 'mu_x', 'sigma_x', 'mu_y', 'sigma_y', 'rho'),
-        'inits': (1., 1., 1., 0.01, 0.05, 0.3),
-        'limits': [],
+        'inits': (200, 1., 1., 0.01, 0.05, 0.3),
+        'limits': [(0, 100000), (0, 5), (0, 5), (-0.3, 1), (0, 0.5), (-1, 1)],
         'name': 'bivariate_gaussian',
         'range': (0., 5., 0., 1.),
     })
+    # r.gStyle.SetNumberContours(20)
     func = myTFunc(conf.name, 'TF2', bigauss, conf.range, conf.names, conf.inits, conf.limits)
     punzi.Fit(func)
     c = r.TCanvas()
     punzi.Draw('surf3')
+    print('Fitted parameters:')
+    print('x ~ N({}, {})'.format(func.GetParameter('mu_x'), func.GetParameter('sigma_x')))
+    print('y ~ N({}, {})'.format(func.GetParameter('mu_y'), func.GetParameter('sigma_y')))