diff --git a/scripts/checkTrackSelection.py b/scripts/checkTrackSelection.py
index 6980fce..9fe079e 100644
--- a/scripts/checkTrackSelection.py
+++ b/scripts/checkTrackSelection.py
@@ -3,7 +3,7 @@
 # @Author: Elena Graverini
 # @Date:   2015-10-27 18:26:49
 # @Last Modified by:   Elena Graverini
-# @Last Modified time: 2017-05-22 15:33:52
+# @Last Modified time: 2017-05-22 16:54:29
 import sys
 import os
 from math import sqrt
@@ -216,6 +216,54 @@
 #     draw_bg = "GhostP : TrChi2/TrNDoF>>hnoise(100, 0., 10., 100, 0., 1.)"
 #     f.Close()
 
+from math import exp, pow
+
+
+def bigauss(x, y, par):
+    # Definition from:
+    # https://en.wikipedia.org/wiki/Multivariate_normal_distribution
+    # Params:
+    # 0 - norm
+    # 1 - mu_x
+    # 2 - sigma_x
+    # 3 - mu_y
+    # 4 - sigma_y
+    # 5 - rho
+    return par[0] * exp(-0.5 / (1 - pow(par[5], 2)) * (pow((x - par[1]) / par[2], 2) + pow((y - par[3]) / par[4], 2) - 2 * par[5] * (x - par[1]) * (y - par[3]) / par[2] / par[4]))
+
+
+def wrap_bigauss(x, par):
+    def get_TF2(xy, par):
+        x, y = xy[0], xy[1]
+        return bigauss(x, y, par)
+    func = get_TF2(x, par)
+    func.SetParNames('norm', 'mu_x', 'sigma_x',
+                     'mu_y', 'sigma_y', 'rho')
+    func.SetParameters(1., 1., 1., 0.01, 0.05, 0.3)
+    return r.TF2('bivariate_gauss', func, 0., 5., 0., 1.)
+
+
+class MyFitFunction():
+    def __init__(self, func):
+        # self.n_var = n_var
+        self.func = func
+
+    def wrap_func(self, x, par):
+        return self.func(*x, par=par)
+
+    def __call__(self, x, par):
+        return self.wrap_func(x, par)
+
+
+def myTFunc(name, classname, func, func_range, par_names, par_inits=[], par_limits=[]):
+    f = r.__getattr__(classname)(name, MyFitFunction(func), *func_range)
+    f.SetParNames(*par_names)
+    if par_inits:
+        f.SetParameters(*par_inits)
+    for i in range(len(par_limits)):
+        f.SetParLimits(i, *par_limits[i])
+    return f
+
 
 def make_all_plots():
     for det in dets:
@@ -231,6 +279,11 @@
             f.Close()
 
 
+# Dictionary which items are accessible with d.key
+class mydict(dict):
+    __getattr__ = dict.__getitem__
+    __setattr__ = dict.__setitem__
+
 if __name__ == '__main__':
     # make_all_plots()
     det = 'TT'
@@ -239,4 +292,14 @@
     t_sig = f.Get('STADCTrackMonitor/HitInfo/%s' % layer[det])
     f, t_sig, h_list = printSelection(f, t_sig, det, fill, stat, good_steps=True, return_all=True)
     punzi = h_list[-1]
-    punzi.Fit('bigaus')
+    conf = mydict({
+        'names': ('norm', 'mu_x', 'sigma_x', 'mu_y', 'sigma_y', 'rho'),
+        'inits': (1., 1., 1., 0.01, 0.05, 0.3),
+        'limits': [],
+        'name': 'bivariate_gaussian',
+        'range': (0., 5., 0., 1.),
+    })
+    func = myTFunc(conf.name, 'TF2', bigauss, conf.range, conf.names, conf.inits, conf.limits)
+    punzi.Fit(func)
+    c = r.TCanvas()
+    punzi.Draw('surf3')