diff --git a/scripts/checkTrackSelection.py b/scripts/checkTrackSelection.py index 1c23d4d..b283eb2 100644 --- a/scripts/checkTrackSelection.py +++ b/scripts/checkTrackSelection.py @@ -3,10 +3,10 @@ # @Author: Elena Graverini # @Date: 2015-10-27 18:26:49 # @Last Modified by: Elena Graverini -# @Last Modified time: 2017-05-24 11:46:38 +# @Last Modified time: 2017-05-24 17:26:51 import sys import os -from math import sqrt +from math import sqrt, radians import ROOT as r location = {'TT': os.path.expandvars("$DISK/data/ST/Aging_Tuples/TT/TTaU/"), @@ -29,14 +29,12 @@ # r.gStyle.SetPadRightMargin(0.16) # r.gStyle.SetPadLeftMargin(r.gStyle.GetPadLeftMargin() * 0.7) -''' This causes segmentation violation in ROOT 6.06 -# r.gROOT.SetBatch(True) -r.gROOT.ProcessLine('.X %s/../include/lhcbstyle.C' % os.getcwd()) +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) @@ -99,16 +97,22 @@ 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) + 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.26, 0.04, 0.9, 0.06, 0, 360, 1.2) + 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) - else: + 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) @@ -117,6 +121,62 @@ 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)} +} + + +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)", "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)", "val{v}<{sc}".format(v=val[det], sc=cut[det]), 'goff', stat, 0) + fill_ranges = ellptical_cuts_dict[det].keys() + print(fill_ranges) + fill_ranges_good = [ra for ra in fill_ranges if ra[0] <= fill <= ra[1]] + print(fill_ranges_good, fill) + fill_range = fill_ranges_good[0] + print(fill_range) + t.Draw("TMath::Sqrt(xHit*xHit+yHit*yHit)/10.0>>h_sel_s(40,0.0,100.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)", "(val{v}<{sc}) && !({cus})".format(v=val[det], sc=cut[det], cus=ellptical_cuts_dict[det][fill_range]), 'goff', stat, 0) + h_tot_s = r.gDirectory.Get('h_tot_s') + h_tot_n = r.gDirectory.Get('h_tot_n') + h_sel_s = r.gDirectory.Get('h_sel_s') + h_sel_n = r.gDirectory.Get('h_sel_n') + 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_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") + 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.SetXTitle("#it{r} [cm]") + h_eff_s.SetYTitle("(Rejection) Efficiency") + h_eff_s.Draw() + h_pur_n.Draw("same") + c_eff.SaveAs('trackSel/c_eff_{}_{}_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) @@ -346,12 +406,16 @@ if __name__ == '__main__': # make_all_plots() - 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] + # 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, 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),