Newer
Older
STAging / scripts / checkTrackSelection.py
@Elena Graverini Elena Graverini on 29 May 2017 19 KB [scripts/trackSel] Add latest plots
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Author: Elena Graverini
# @Date:   2015-10-27 18:26:49
# @Last Modified by:   Elena Graverini
# @Last Modified time: 2017-05-29 11:24:56
import sys
import os
from math import sqrt, radians
import ROOT as r

location = {'TT': os.path.expandvars("$DISK/data/ST/Aging_Tuples/TT/TTaU/"),
            'IT': os.path.expandvars("$DISK/data/ST/Aging_Tuples/IT/T3X2/")}
out_location = os.getcwd() + '/trackSel'
os.system('mkdir -p %s' % out_location)

# Load fill numbers
macros = os.path.expandvars('$CCEHOME/macros/CCEScan/')
with open(macros + 'Fills.dat', 'rb') as f:
    fills = f.read().splitlines()

# fills = [3478, 4518, 4643, 4856, 5162, 5448]
dets = ['TT', 'IT']
val = {'TT': 5, 'IT': 7}
layer = {'TT': 'TTaU', 'IT': 'T3X2'}
cut = {'TT': 12, 'IT': 20}
stat_2d = 200000000
stat = stat_2d  # 2000000

# r.gStyle.SetPadRightMargin(0.16)
# r.gStyle.SetPadLeftMargin(r.gStyle.GetPadLeftMargin() * 0.7)

r.gROOT.SetBatch(True)
r.gROOT.ProcessLine('.X %s/../include/lhcbStyle.C' % os.getcwd())
r.gStyle.SetPadLeftMargin(0.20)
r.gStyle.SetPadRightMargin(0.25)
r.gStyle.SetPadBottomMargin(0.15)
r.gStyle.SetPadTopMargin(0.05)

# 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):
    new = h.Clone()
    nbins = h.GetNbinsX() * h.GetNbinsY()
    for i in range(nbins):
        new.SetBinContent(i, sqrt(h.GetBinContent(i)))
    return new


def s_over_b(s, b):
    h = s.Clone()
    h.Divide(b)
    h.name = 's_over_b'
    return h


def s_over_sqrtb(s, b):
    h = s.Clone()
    sqrtb = sqrt_of_2d_histogram(b)
    h.Divide(sqrtb)
    h.name = 's_over_sqrtb'
    return h


def punzi_significance(s, b):
    h = s.Clone()
    sb = s.Clone()
    sb.Add(b)
    sqrt_sb = sqrt_of_2d_histogram(sb)
    h.Divide(sqrt_sb)
    h.name = 'punzi_significance'
    return h


def put_fill(c, fill):
    label = r.TPaveText( 0.74 - r.gStyle.GetPadRightMargin(), 0.87 - r.gStyle.GetPadTopMargin(),
                         0.95 - r.gStyle.GetPadRightMargin(), 0.95 - r.gStyle.GetPadTopMargin(), "BRNDC")
    label.SetFillColor(0)
    label.SetTextAlign(12)
    label.SetBorderSize(0)
    label.SetTextFont(132)
    label.SetTextSize(0.06)
    label.SetTextAlign(22)
    label.AddText('Fill %s' % fill)
    label.Draw('same')
    c.Modified()
    c.Update()
    return c


def sqrt_histogram(h):
    sh = h.Clone()
    for i in range(h.GetNbinsX()):
        sh.SetBinContent(i, h.GetBinError(i))
    sh.Sumw2()
    return sh


csaver = []


def putEllipse(c, det, fill):
    if 'TT' in det:
        if fill > 3478:
            ellipse = r.TEllipse(1.35, 0.04, 1.1, 0.08, 0, 360, 1.5)
        elif fill > 2252:
            ellipse = r.TEllipse(1.26, 0.04, 1.0, 0.07, 0, 360, 2.7)
        elif fill > 2083:
            ellipse = r.TEllipse(1.15, 0.05, 1.0, 0.1, 0, 360, 5.)
        else:
            ellipse = r.TEllipse(1.3, 0.05, 1.4, 0.11, 0, 360, 4.1)
    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)
        elif fill > 1616:
            ellipse = r.TEllipse(1.2, 0.04, 0.9, 0.07, 0, 360, 4.7)
        else:
            ellipse = r.TEllipse(1.3, 0.04, 1.8, 0.11, 0, 360, 6.2)
    ellipse.SetFillStyle(0)
    ellipse.SetLineWidth(4)
    ellipse.SetLineColor(r.kRed - 4)
    c.cd()
    ellipse.Draw('same')
    return c, ellipse


def elliptical_cut(cx, cy, rx, ry, theta):
    theta = radians(theta)
    cus = "(((xprime-{cx})^2/{rx}^2 + (yprime-{cy})^2/{ry}^2) < 1.)".format(
          cx=cx, rx=rx, cy=cy, ry=ry)
    cus = cus.replace('xprime', 'x*TMath::Cos({t})-y*TMath::Sin({t})'.format(t=theta))
    cus = cus.replace('yprime', 'y*TMath::Cos({t})-x*TMath::Sin({t})'.format(t=theta))
    cus = cus.replace('x', '(TrChi2/TrNDoF)')
    cus = cus.replace('y', 'GhostP')
    return cus


ellptical_cuts_dict = {
    'TT': {
        (3478 + 1, 9999): elliptical_cut(1.35, 0.04, 1.1, 0.08, 1.5),
        (2252 + 1, 3478): elliptical_cut(1.26, 0.04, 1.0, 0.07, 2.7),
        (2083 + 1, 2252): elliptical_cut(1.15, 0.05, 1.0, 0.1, 5.),
        (0000 + 1, 2083): elliptical_cut(1.3, 0.05, 1.4, 0.11, 4.1)},
    'IT': {
        (4856 + 1, 9999): elliptical_cut(2.05, 0.04, 1.9, 0.04, 1.9),
        (3478 + 1, 4856): elliptical_cut(1.4, 0.025, 1.1, 0.04, 1.1),
        (1616 + 1, 3478): elliptical_cut(1.2, 0.04, 0.9, 0.07, 4.7),
        (0000 + 1, 1616): elliptical_cut(1.3, 0.04, 1.8, 0.11, 6.2)}
}

house_cuts = {
    'TT': "(TrChi2/TrNDoF<3.0) && (GhostP<0.01+0.1/1.6*TrChi2/TrNDoF) && (GhostP<0.01-0.1/3.4*(TrChi2/TrNDoF-5.0))",
    'IT': "(GhostP<0.1*TrChi2/TrNDoF) && (GhostP<0.2-0.1*(TrChi2/TrNDoF-5.0))",
}


def test_selection(t, fill, det, stat):
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_tot_s(40,0.0,100.0)", "(odinStep==0) && (val{v}>{sc})".format(v=val[det], sc=cut[det]), 'goff', stat, 0)
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_tot_n(40,0.0,100.0)", "(odinStep==0) && (val{v}<{sc})".format(v=val[det], sc=cut[det]), 'goff', stat, 0)
    h_tot_s = r.gDirectory.Get('h_tot_s')
    h_tot_n = r.gDirectory.Get('h_tot_n')
    fill_ranges = ellptical_cuts_dict[det].keys()
    fill_ranges_good = [ra for ra in fill_ranges if ra[0] <= fill <= ra[1]]
    fill_range = fill_ranges_good[0]
    assert(fill_range[0] <= fill <= fill_range[1])
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_s(40,0.0,100.0)", "(odinStep==0) && (val{v}>{sc}) && ({cus})".format(v=val[det], sc=cut[det], cus=ellptical_cuts_dict[det][fill_range]), 'goff', stat, 0)
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_n(40,0.0,100.0)", "(odinStep==0) && (val{v}<{sc}) && !({cus})".format(v=val[det], sc=cut[det], cus=ellptical_cuts_dict[det][fill_range]), 'goff', stat, 0)
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_tot(40,0.0,100.0)", "(odinStep==0) && ({cus})".format(cus=ellptical_cuts_dict[det][fill_range]), 'goff', stat, 0)
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_s_christian(40,0.0,100.0)", "(odinStep==0) && (val{v}>{sc}) && ({cus})".format(v=val[det], sc=cut[det], cus=house_cuts[det]), 'goff', stat, 0)
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_n_christian(40,0.0,100.0)", "(odinStep==0) && (val{v}<{sc}) && !({cus})".format(v=val[det], sc=cut[det], cus=house_cuts[det]), 'goff', stat, 0)
    t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_tot_christian(40,0.0,100.0)", "(odinStep==0) && ({cus})".format(cus=house_cuts[det]), 'goff', stat, 0)
    h_sel_s = r.gDirectory.Get('h_sel_s')
    h_sel_n = r.gDirectory.Get('h_sel_n')
    h_sel_tot = r.gDirectory.Get('h_sel_tot')
    h_sel_s_christian = r.gDirectory.Get('h_sel_s_christian')
    h_sel_n_christian = r.gDirectory.Get('h_sel_n_christian')
    h_sel_tot_christian = r.gDirectory.Get('h_sel_tot_christian')
    h_eff_s = r.TH1F('h_eff_s', 'h_eff_s', 40, 0.0, 100.0)
    h_pur_n = r.TH1F('h_pur_n', 'h_pur_n', 40, 0.0, 100.0)
    h_snr = r.TH1F('h_snr', 'h_snr', 40, 0.0, 100.0)
    h_sign = r.TH1F('h_sign', 'h_sign', 40, 0.0, 100.0)
    h_punzi = r.TH1F('h_punzi', 'h_punzi', 40, 0.0, 100.0)
    h_eff_s.Divide(h_sel_s, h_tot_s, 1.0, 1.0, "b")
    h_pur_n.Divide(h_sel_n, h_tot_n, 1.0, 1.0, "b")
    h_snr.Divide(h_sel_s, h_sel_n, 1.0, 1.0, "b")
    h_sign.Divide(h_sel_s, sqrt_histogram(h_sel_n), 1.0, 1.0, "b")
    h_punzi.Divide(h_sel_s, sqrt_histogram(h_sel_tot), 1.0, 1.0, "b")
    h_eff_s_christian = r.TH1F('h_eff_s_christian', 'h_eff_s_christian', 40, 0.0, 100.0)
    h_pur_n_christian = r.TH1F('h_pur_n_christian', 'h_pur_n_christian', 40, 0.0, 100.0)
    h_snr_christian = r.TH1F('h_snr_christian', 'h_snr_christian', 40, 0.0, 100.0)
    h_sign_christian = r.TH1F('h_sign_christian', 'h_sign_christian', 40, 0.0, 100.0)
    h_punzi_christian = r.TH1F('h_punzi_christian', 'h_punzi_christian', 40, 0.0, 100.0)
    h_snr_christian.Divide(h_sel_s_christian, h_sel_n_christian, 1.0, 1.0, "b")
    h_sign_christian.Divide(h_sel_s_christian, sqrt_histogram(h_sel_n_christian), 1.0, 1.0, "b")
    h_punzi_christian.Divide(h_sel_s_christian, sqrt_histogram(h_sel_tot_christian), 1.0, 1.0, "b")
    h_eff_s_christian.Divide(h_sel_s_christian, h_tot_s, 1.0, 1.0, "b")
    h_pur_n_christian.Divide(h_sel_n_christian, h_tot_n, 1.0, 1.0, "b")
    h_pur_n.SetLineColor(r.kRed)
    h_pur_n_christian.SetLineColor(r.kRed)
    h_eff_s_christian.SetLineStyle(r.kDashed)
    h_pur_n_christian.SetLineStyle(r.kDashed)
    c_eff = r.TCanvas("c_eff", "Eff./Pur. Canvas", 800, 600)
    h_eff_s.SetMaximum(1.1)
    h_eff_s.SetMinimum(-0.1)
    h_pur_n.SetMaximum(1.1)
    h_pur_n.SetMinimum(-0.1)
    h_eff_s_christian.SetXTitle("#it{r} [cm]")
    h_eff_s_christian.SetYTitle("(Rejection) Efficiency")
    h_eff_s_christian.SetMaximum(1.0)
    c_eff.cd()
    h_eff_s_christian.Draw()
    h_pur_n_christian.Draw("same")
    h_eff_s.Draw("same")
    h_pur_n.Draw("same")
    c_eff.SaveAs('trackSel/c_eff_{}_{}_elena.pdf'.format(fill, det))
    c_snr = r.TCanvas()
    c_snr.cd()
    h_snr_christian.SetXTitle("#it{r} [cm]")
    h_snr_christian.SetYTitle("S/N ratio")
    h_snr_christian.SetLineStyle(r.kDashed)
    h_snr_christian.Draw('hist')
    h_snr.Draw('hist same')
    c_snr.SaveAs('trackSel/c_snr_{}_{}_elena.pdf'.format(fill, det))
    c_sign = r.TCanvas()
    c_sign.cd()
    h_sign_christian.SetXTitle("#it{r} [cm]")
    h_sign_christian.SetYTitle("S/#sqrt{N}")
    h_sign_christian.SetLineStyle(r.kDashed)
    h_sign_christian.Draw('hist')
    h_sign.Draw('hist same')
    c_sign.SaveAs('trackSel/c_sign_{}_{}_elena.pdf'.format(fill, det))
    c_punzi = r.TCanvas()
    c_punzi.cd()
    h_punzi_christian.SetXTitle("#it{r} [cm]")
    h_punzi_christian.SetYTitle("S/#sqrt{S+N}")
    h_punzi_christian.SetLineStyle(r.kDashed)
    h_punzi_christian.Draw('hist')
    h_punzi.Draw('hist same')
    c_punzi.SaveAs('trackSel/c_punzi_{}_{}_elena.pdf'.format(fill, det))


def printSelection(f, t_sig, det, fill, stat, good_steps=True, return_all=False):
    r.gStyle.SetPadLeftMargin(0.20)
    r.gStyle.SetPadRightMargin(0.25)
    r.gStyle.SetPadBottomMargin(0.15)
    r.gStyle.SetPadTopMargin(0.05)
    cut_sig = "(val%s>%s)" % (val[det], cut[det])
    cut_noise = "(val%s<%s)" % (val[det], cut[det])
    if good_steps:
        cut_sig += ' && (odinStep==0)'
        cut_noise += ' && (odinStep==0)'
    print(t_sig.Draw("GhostP : TrChi2/TrNDoF>>hsig(100,0., 5,100, 0.,1.)", cut_sig, "goff", stat, 0))
    h_sig = r.gDirectory.Get('hsig')
    print(t_sig.Draw("GhostP : TrChi2/TrNDoF>>hnoise(100,0., 5,100, 0.,1.)", cut_noise, "goff", stat, 0))
    h_noise = r.gDirectory.Get('hnoise')
    # c = r.TCanvas()
    h_s_over_b = s_over_b(h_sig, h_noise)
    h_s_over_b.GetZaxis().SetRangeUser(0, 10)
    h_s_over_b.GetZaxis().SetTitle('S/B ratio')
    h_s_over_sqrtb = s_over_sqrtb(h_sig, h_noise)
    h_s_over_sqrtb.GetZaxis().SetRangeUser(0, 100)
    h_s_over_sqrtb.GetZaxis().SetTitle('S/#sqrt{B}')
    h_punzi_significance = punzi_significance(h_sig, h_noise)
    h_punzi_significance.GetZaxis().SetRangeUser(0, 100)
    h_punzi_significance.GetZaxis().SetTitle('S/#sqrt{S+B}')
    for h in [h_s_over_b, h_s_over_sqrtb, h_punzi_significance]:
        csaver.append(r.TCanvas())
        csaver[-1].cd()
        h.GetXaxis().SetTitle(r'#chi^{2}_{track}/ndf')
        h.GetYaxis().SetTitle(r'Ghost Prob.')
        h.Draw('colz')
        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')
        csaver[-1].SaveAs(save_name)
        if not return_all:
            h.Delete()
    # h_sig.Divide(h_noise)
    # h_sig.GetZaxis().SetRangeUser(0, 10)
    # h_sig.SetTitle('TTaU S/B ratio for fill %s' % fill)
    # h_sig.GetXaxis().SetTitle(r'#chi^{2}_{track}/ndf')
    # h_sig.GetYaxis().SetTitle(r'Ghost Prob.')
    # h_sig.GetZaxis().SetTitle(r'S/B ratio')
    # h_sig.Draw("colz")
    # c.Update()
    # c, line = putLine(c, det)
    # c = put_fill(c, fill)
    # c.SaveAs('%s/%s_fill%s.pdf' % (out_location, det, fill))
    if not return_all:
        h_sig.Delete()
        h_noise.Delete()
        return f, t_sig
    else:
        return f, t_sig, [h_s_over_b, h_s_over_sqrtb, h_punzi_significance]


def printADC(f, t_sig, det, fill, stat, step):
    if step > 10:
        print(t_sig.Draw("val%s>>hsig%s(210,-35, 175)" % (val[det], step), "odinStep==%s" % step, "goff"))  # , stat, 0))
    else:
        print(t_sig.Draw("val%s>>hsig%s(210,-35, 175)" % (val[det], step), "odinStep==%s" % step, "goff", stat, 0))
    h_sig = r.gDirectory.Get('hsig%s' % step)
    r.gROOT.ProcessLine('.X %s/../include/lhcbstyle.C' % os.getcwd())
    r.gStyle.SetPadLeftMargin(0.15)
    c = r.TCanvas()
    h_sig.GetXaxis().SetTitle('Signal height [ADC value]')
    h_sig.GetYaxis().SetTitle('Fraction of hits / ADC value')
    h_sig.DrawNormalized()
    c.Update()
    label = r.TPaveText( 0.74 - r.gStyle.GetPadRightMargin(), 0.80 - r.gStyle.GetPadTopMargin(),
                         0.95 - r.gStyle.GetPadRightMargin(), 0.95 - r.gStyle.GetPadTopMargin(), "BRNDC")
    label.SetFillColor(0)
    label.SetTextAlign(12)
    label.SetBorderSize(0)
    label.SetTextFont(132)
    label.SetTextSize(0.06)
    label.SetTextAlign(22)
    label.AddText('Fill %s\n' % fill)
    from checkPulseData import get_vmap
    v_bias = get_vmap(det)[int(step) / 6]
    label.AddText('V_{bias} = %s V' % v_bias)
    label.Draw('same')
    line = r.TLine()
    line.SetLineWidth(3)
    line.SetLineColor(r.kRed)
    line.SetLineStyle(7)
    line.DrawLine(float(cut[det]), 0.0, float(cut[det]), r.gPad.GetUymax())
    c.Modified()
    c.Update()
    c.SaveAs('%s/%s_odinStep%s_fill%s.pdf' % (out_location, det, step, fill))
    h_sig.Delete()
    return f, t_sig


def putLine(c, det):
    line = r.TLine()
    line.SetLineWidth(4)
    line.SetLineColor(r.kWhite)
    c.cd()
    if 'TT' in det:
        line.DrawLine(0.0, 0.01, 1.6, 0.16)
        line.DrawLine(3.0, 0.07, 1.6, 0.16)
        line.DrawLine(3.0, 0.07, 3.0, 0.0)
    else:
        line.DrawLine(0.0, 0.0, 3.5, 0.35)
        line.DrawLine(3.5, 0.35, 5.0, 0.20)
    return c, line


# def printByStep(fill, stat, det):
#     f = r.TFile(location[det] + '%s.root' % fill, 'read')
#     t_sig = f.Get('STADCTrackMonitor/HitInfo/TTaU')
#
#     # Look at Landau distribution
#     for step in range(67):
#         c = r.TCanvas()
#         t_sig.Draw('val5>>hnew(70,-30,90)', 'odinStep==%s' % step)
#         c.SaveAs()
#
#     draw_sig = "GhostP : TrChi2/TrNDoF>>hsig(100, 0., 10., 100, 0., 1.)"
#     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 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]
        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, n_var):
        self.n_var = n_var
        self.func = func

    def wrap_func(self, x, 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)


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, int(classname[2])), *(func_range + (n_pars,)))
    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:
        for fill in fills:
            f = r.TFile(location[det] + '%s.root' % fill, 'read')
            t_sig = f.Get('STADCTrackMonitor/HitInfo/%s' % layer[det])

            f, t_sig = printSelection(f, t_sig, det, fill, stat, good_steps=True)
            if '--steps' in sys.argv:
                for step in [0, 12, 48, 54, 60]:
                    f, t_sig = printADC(f, t_sig, det, fill, stat, step)

            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'
    # fill = '2252'
    for det in dets:
        for fill in fills:
            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_2d, good_steps=True, return_all=True)
            test_selection(t_sig, int(fill), det, stat)
            f.Close()
    # punzi = h_list[-1]
    # 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')