diff --git a/scripts/checkTrackSelection.py b/scripts/checkTrackSelection.py
index 4149ab9..1c23d4d 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-24 09:34:22
+# @Last Modified time: 2017-05-24 11:46:38
 import sys
 import os
 from math import sqrt
@@ -96,6 +96,27 @@
 csaver = []
 
 
+def putEllipse(c, det, fill):
+    if 'TT' in det:
+        if fill > 3478:
+            ellipse = r.TEllipse(1.26, 0.04, 0.9, 0.06, 0, 360, 1.2)
+        else:
+            ellipse = r.TEllipse(1.26, 0.04, 0.9, 0.06, 0, 360, 1.2)
+    elif 'IT' in det:
+        if fill > 4856:
+            ellipse = r.TEllipse(2.05, 0.04, 1.9, 0.04, 0, 360, 1.9)
+        elif fill > 3478:
+            ellipse = r.TEllipse(1.4, 0.025, 1.1, 0.04, 0, 360, 1.1)
+        else:
+            ellipse = r.TEllipse(1.2, 0.04, 0.9, 0.07, 0, 360, 4.7)
+    ellipse.SetFillStyle(0)
+    ellipse.SetLineWidth(4)
+    ellipse.SetLineColor(r.kRed - 4)
+    c.cd()
+    ellipse.Draw('same')
+    return c, ellipse
+
+
 def printSelection(f, t_sig, det, fill, stat, good_steps=True, return_all=False):
     r.gStyle.SetPadLeftMargin(0.20)
     r.gStyle.SetPadRightMargin(0.25)
@@ -129,6 +150,11 @@
         csaver[-1].Update()
         csaver[-1], line = putLine(csaver[-1], det)
         csaver[-1] = put_fill(csaver[-1], fill)
+
+        csaver[-1], ellipse = putEllipse(csaver[-1], det, int(fill))
+        h.line = line
+        h.ellipse = ellipse
+
         save_name = '%s/%s_%s_fill%s.pdf' % (out_location, det, h.name, fill)
         if good_steps:
             save_name = save_name.replace('.pdf', '_odinStep_lt_12.pdf')
@@ -237,6 +263,28 @@
     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 double_bigauss(x, y, par):
+    # Bivariate double gaussian
+    # 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
+    # 6 - frac
+    # 7 - sigma_x_2
+    # 8 - sigma_y_2
+    norm, mu_x, sigma_x = par[0], par[1], par[2]
+    mu_y, sigma_y, rho = par[3], par[4], par[5]
+    frac, sigma_x_2, sigma_y_2 = par[6], par[7], par[8]
+    return norm * exp(-0.5 / (1 - pow(rho, 2)) * ((frac * pow((x - mu_x) / sigma_x, 2) +
+           (1 - frac) * pow((x - mu_x) / sigma_x_2, 2)) + (frac * pow((y - mu_y) / sigma_y, 2) +
+           (1 - frac) * pow((y - mu_y) / sigma_y_2, 2)) - 2 * rho * frac * (x - mu_x) * (y - mu_y) / (sigma_x * sigma_y)) -
+           2 * rho * (1 - frac) * (x - mu_x) * (y - mu_y) / (sigma_x_2 * sigma_y_2))
+
+
 def wrap_bigauss(x, par):
     def get_TF2(xy, par):
         x, y = xy[0], xy[1]
@@ -298,24 +346,38 @@
 
 if __name__ == '__main__':
     # make_all_plots()
-    det = 'TT'
-    fill = '5448'
+    det = 'IT'
+    fill = '1616'
     f = r.TFile(location[det] + '%s.root' % fill, 'read')
     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]
-    conf = mydict({
-        'names': ('norm', 'mu_x', 'sigma_x', 'mu_y', 'sigma_y', 'rho'),
-        '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')))
+    # conf = mydict({
+    #     'names': ('norm', 'mu_x', 'sigma_x', 'mu_y', 'sigma_y', 'rho', 'frac', 'sigma_x_2', 'sigma_y_2'),
+    #     'inits': (200, 1., 0.6, 0.01, 0.1, 0.7, 0.5, 0.2, 0.2),
+    #     'limits': [(0, 100000), (-1.5, 5), (0, 5), (-0.5, 1), (0, 0.5), (-1, 1), (0, 1), (0, 5), (0, 1)],
+    #     'name': 'bivariate_gaussian',
+    #     'range': (0., 4., 0., 0.4),
+    # })
+    # # r.gStyle.SetNumberContours(20)
+    # func = myTFunc(conf.name, 'TF2', double_bigauss, conf.range, conf.names, conf.inits, conf.limits)
+    # res = punzi.Fit(func, 'RS')
+    # c = r.TCanvas()
+    # punzi.Draw('surf3')
+    # c2 = r.TCanvas()
+    # punzi.Draw('cont list')
+    # print('Fitted parameters:')
+    # print('x ~ N({mu}, {s1}) + N({mu}, {s2})'.format(mu=func.GetParameter('mu_x'), s1=func.GetParameter('sigma_x'), s2=func.GetParameter('sigma_x_2')))
+    # print('y ~ N({mu}, {s1}) + N({mu}, {s2})'.format(mu=func.GetParameter('mu_y'), s1=func.GetParameter('sigma_y'), s2=func.GetParameter('sigma_y_2')))
+    # print('The ratio between the Gaussians is ~ {f}'.format(f=func.GetParameter('frac')))
+    # print('Correlation rho is ~ {rho}'.format(rho=func.GetParameter('rho')))
+    # res.Print()
+    # conts = r.gROOT.GetListOfSpecials().FindObject('contours')
+    #
+    # c3 = r.TCanvas()
+    # punzi.Draw('colz')
+    # ellipse = r.TEllipse(1.46, 0.04, 0.99, 0.1, 0, 360, 3)
+    # ellipse.SetFillStyle(0)
+    # ellipse.SetLineWidth(3)
+    # ellipse.SetLineColor(r.kRed - 4)
+    # ellipse.Draw('same')