ATLAS Offline Software
dqt_zlumi_alleff_HIST.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3 import sys, os, glob
4 import ROOT
5 
6 ROOT.gStyle.SetOptStat(0)
7 
8 #ACCEPTANCE = 3.173927e-01
9 ACCEPTANCE = 3.323224e-01
10 
11 import argparse
12 parser = argparse.ArgumentParser()
13 parser.add_argument('infile', type=str, help='input HIST file')
14 parser.add_argument('--out', type=str, help='output ROOT file')
15 parser.add_argument('--plotdir', type=str, help='Directory to dump plots',
16  default='plots')
17 parser.add_argument('--debug', action='store_true', help='Be verbose in output')
18 parser.add_argument('--mode', type=str, help='Zee or Zmumu')
19 
20 
21 args = parser.parse_args()
22 
23 infilename = args.infile
24 infile = ROOT.TFile.Open(infilename, 'READ')
25 
26 runmode = args.mode
27 print('Running in', runmode, 'mode')
28 
29 
30 runname = None
31 for k in infile.GetListOfKeys():
32  if k.GetName().startswith('run_'):
33  runname = k.GetName()
34  break
35 if not runname:
36  print('Cannot find run directory in input file')
37  sys.exit(1)
38 else:
39  print('Found runname', runname)
40 
41 lbdirs = []
42 for k in infile.Get(runname).GetListOfKeys():
43  if k.GetName().startswith('lb_'):
44  lbdirs.append(k.GetName())
45 
46 print('Now to dump')
47 lbnums = sorted([int(_[3:]) for _ in lbdirs])
48 
49 effcyt = ROOT.TH1F('effcyt', 'Trigger efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
50  lbnums[-1]+0.5)
51 effcyr = ROOT.TH1F('effcyr', 'Loose muon reco efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
52  lbnums[-1]+0.5)
53 effcya = ROOT.TH1F('effcya', 'Combined acc x efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
54  lbnums[-1]+0.5)
55 effcydir = ROOT.TH1F('effcydir', 'Direct acc x efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
56  lbnums[-1]+0.5)
57 
58 from array import array
59 fout = ROOT.TFile(args.out if args.out else '%s_all.root' % runname[4:], 'RECREATE')
60 o_run = array('I', [0])
61 o_lb = array('I', [0])
62 o_lbwhen = array('d', [0., 0.])
63 o_z_one = array('f', [0.])
64 o_z_two = array('f', [0.])
65 o_trigeff = array('f', [0.])
66 o_trigeffstat = array('f', [0.])
67 o_recoeff = array('f', [0.])
68 o_recoeffstat = array('f', [0.])
69 o_alleff = array('f', [0.])
70 o_alleffstat = array('f', [0.])
71 o_ae = array('f', [0.])
72 o_aestat = array('f', [0.])
73 tl = ROOT.TTree( 'lumitree', 'Luminosity tree' )
74 tl.Branch('run', o_run, 'run/i')
75 tl.Branch('lb', o_lb, 'lb/i')
76 tl.Branch('lbwhen', o_lbwhen, 'lbwhen[2]/D')
77 tl.Branch('z_one', o_z_one, 'z_one/F')
78 tl.Branch('z_two', o_z_two, 'z_two/F')
79 tl.Branch('trigeff', o_trigeff, 'trigeff/F')
80 tl.Branch('trigeffstat', o_trigeffstat, 'trigeffstat/F')
81 tl.Branch('recoeff', o_recoeff, 'recoeff/F')
82 tl.Branch('recoeffstat', o_recoeffstat, 'recoeffstat/F')
83 tl.Branch('alleff', o_alleff, 'alleff/F')
84 tl.Branch('alleffstat', o_alleffstat, 'alleffstat/F')
85 tl.Branch('ae', o_ae, 'ae/F')
86 tl.Branch('aestat', o_aestat, 'aestat/F')
87 
88 
89 from DQUtils import fetch_iovs
90 #rset=set(_[0] for _ in rlb)
91 #print(rset)
92 lblb = fetch_iovs("LBLB", runs=int(runname[4:])).by_run
93 for lb in sorted(lbdirs):
94  if runmode == "Zee":
95  h = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_eltrigtp_matches' % (runname, lb))
96  hmo = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
97  hms = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
98  hno = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
99  hns = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
100  if runmode == "Zmumu":
101  h = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_mutrigtp_matches' % (runname, lb))
102  hmo = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_os' % (runname, lb))
103  hms = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_ss' % (runname, lb))
104  hno = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_os' % (runname, lb))
105  hns = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_ss' % (runname, lb))
106  lbnum = int(lb[3:])
107  yld = (h[2], h[3])
108  ylderr = (h.GetBinError(2), h.GetBinError(3))
109  #print(yld, ylderr)
110  A, B = yld
111  o_z_one[0], o_z_two[0] = yld
112  if B == 0: continue
113  eff = 1./(float(A)/B/2.+1.)
114  inverrsq = ((1/2./B)*ylderr[0])**2+((A/2./B**2)*ylderr[1])**2
115  o_trigeff[0] = eff
116  o_trigeffstat[0] = (inverrsq**.5)*(eff**2)
117  o_run[0], o_lb[0] = int(runname[4:]), lbnum
118  try:
119  iov = lblb[int(runname[4:])][lbnum-1]
120  o_lbwhen[0], o_lbwhen[1] = iov.StartTime/1e9, iov.EndTime/1e9
121  except Exception as e:
122  o_lbwhen[0], o_lbwhen[1] = 0, 0
123  effcyt.SetBinContent(lbnum-lbnums[0]+1, eff)
124  effcyt.SetBinError(lbnum-lbnums[0]+1, o_trigeffstat[0])
125 
126  def extract(histogram):
127  dbl = ROOT.Double()
128  rv1 = histogram.IntegralAndError(21, 30, dbl)
129  return (rv1, float(dbl))
130  matchos, matchoserr = extract(hmo)
131  matchss, matchsserr = extract(hms)
132  nomatchos, nomatchoserr = extract(hno)
133  nomatchss, nomatchsserr = extract(hns)
134  if args.debug:
135  print(lb)
136  print(' ->', matchos, matchoserr)
137  print(' ->', matchss, matchsserr)
138  print(' ->', nomatchos, nomatchoserr)
139  print(' ->', nomatchss, nomatchsserr)
140  A = float(matchos-matchss)
141  Aerr = (matchoserr**2+matchsserr**2)**.5
142  B = float(nomatchos-nomatchss)
143  Berr = (nomatchoserr**2+nomatchsserr**2)**.5
144  if Berr == 0: Berr = 1.
145  if A == 0 or B/A == -1:
146  eff = 1.
147  inverrsq = 1.
148  else:
149  eff = 1./(1.+B/A)
150  inverrsq = ((-B/A**2)*Aerr)**2+((1./A)*Berr)**2
151  o_recoeff[0] = eff
152  o_recoeffstat[0] = (inverrsq**.5)*(eff**2)
153  effcyr.SetBinContent(lbnum-lbnums[0]+1, eff)
154  effcyr.SetBinError(lbnum-lbnums[0]+1, o_recoeffstat[0])
155 
156  o_ae[0] = ACCEPTANCE*(1-(1-o_trigeff[0])**2)*(o_recoeff[0])**2
157  o_aestat[0] = ACCEPTANCE*((o_recoeff[0]**2*2*(1-o_trigeff[0])*o_trigeffstat[0])**2+(2*o_recoeff[0]*(1-(1-o_trigeff[0])**2)*o_recoeffstat[0])**2)**.5
158  o_alleff[0] = (1-(1-o_trigeff[0])**2)*(o_recoeff[0])**2
159  o_alleffstat[0] = ((o_recoeff[0]**2*2*(1-o_trigeff[0])*o_trigeffstat[0])**2+(2*o_recoeff[0]*(1-(1-o_trigeff[0])**2)*o_recoeffstat[0])**2)**.5
160  effcya.SetBinContent(lbnum-lbnums[0]+1, o_ae[0])
161  effcya.SetBinError(lbnum-lbnums[0]+1, o_aestat[0])
162 
163 
164  tl.Fill()
165 tl.Write()
166 print('Done')
167 
168 c1 = ROOT.TCanvas()
169 effcya.SetMarkerStyle(21)
170 effcya.SetMarkerColor(ROOT.kBlue)
171 effcya.GetYaxis().SetRangeUser(0.25,0.31)
172 effcya.Draw('PE')
173 c1.Print(os.path.join(args.plotdir, '%s_combined_efficiency.eps' % runname[4:]))
174 fout.WriteTObject(effcya)
175 c1.Clear()
176 effcyt.SetMarkerStyle(21)
177 effcyt.SetMarkerColor(ROOT.kBlue)
178 effcyt.GetYaxis().SetRangeUser(0.66,0.86)
179 effcyt.Draw('PE')
180 c1.Print(os.path.join(args.plotdir, '%s_trigger_efficiency.eps' % runname[4:]))
181 fout.WriteTObject(effcyt)
182 c1.Clear()
183 effcyr.SetMarkerStyle(21)
184 effcyr.SetMarkerColor(ROOT.kBlue)
185 effcyr.GetYaxis().SetRangeUser(0.9,1.0)
186 effcyr.Draw('PE')
187 c1.Print(os.path.join(args.plotdir, '%s_reco_efficiency.eps' % runname[4:]))
188 fout.WriteTObject(effcyr)
189 fout.Close()
190 
191 if runmode == "Zee":
192  sumweights = infile.Get('%s/GLOBAL/DQTDataFlow/m_sumweights' % runname)
193  ctr = infile.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_el' % runname)
194 if runmode == "Zmumu":
195  sumweights = infile.Get('%s/GLOBAL/DQTDataFlow/m_sumweights' % runname)
196  ctr = infile.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_mu' % runname)
197 if sumweights:
198  for ibin in xrange(1,sumweights.GetNbinsX()+1):
199  o_lb[0] = int(sumweights.GetBinCenter(ibin))
200  ctrbin = ctr.FindBin(o_lb[0])
201  print(ibin, o_lb[0], sumweights[ibin], ctr[ctrbin])
202  if sumweights[ibin] == 0: continue
203  p = ctr[ctrbin]/sumweights[ibin]
204  o_alleff[0]=p
205  try:
206  o_alleffstat[0]=(p*(1-p))**.5*(sumweights.GetBinError(ibin)/sumweights[ibin])
207  except ValueError:
208  o_alleffstat[0]=(sumweights.GetBinError(ibin)/sumweights[ibin])
209  effcydir.SetBinContent(effcydir.FindBin(o_lb[0]), p)
210  effcydir.SetBinError(effcydir.FindBin(o_lb[0]), o_alleffstat[0])
211 
212  effcya.GetYaxis().SetRangeUser(0.27,0.31)
213  effcya.Draw('PE')
214  effcydir.SetMarkerStyle(20)
215  effcydir.SetMarkerColor(ROOT.kRed)
216  effcydir.Draw('SAME,PE')
217  leg=ROOT.TLegend(0.65, 0.7, 0.89, 0.89)
218  leg.AddEntry(effcya, 'Predicted A#epsilon', 'PE')
219  leg.AddEntry(effcydir, 'Actual A#epsilon', 'PE')
220  leg.Draw()
221  c1.Print(os.path.join(args.plotdir, '%s_tp_comparison.eps' % runname[4:]))
222 
223  effcyrat=effcydir.Clone()
224  effcyrat.Divide(effcya)
225  effcyrat.SetTitle('MC Correction Factor')
226  effcyrat.SetXTitle('<#mu>')
227  effcyrat.Draw('PE')
228  effcyrat.Fit('pol1')
229  c1.Print(os.path.join(args.plotdir, '%s_tp_correction.eps' % runname[4:]))
xrange
void xrange(TH1 *h, bool symmetric)
Definition: computils.cxx:515
python.db.fetch_iovs
def fetch_iovs(folder_name, since=None, until=None, channels=None, tag="", what="all", max_records=-1, with_channel=True, loud=False, database=None, convert_time=False, named_channels=False, selection=None, runs=None, with_time=False, unicode_strings=False)
Definition: DQUtils/python/db.py:67
DerivationFramework::TriggerMatchingUtils::sorted
std::vector< typename T::value_type > sorted(T begin, T end)
Helper function to create a sorted vector from an unsorted one.
array
dqt_zlumi_alleff_HIST.extract
def extract(histogram)
Definition: dqt_zlumi_alleff_HIST.py:126
readCCLHist.int
int
Definition: readCCLHist.py:84
Muon::print
std::string print(const MuPatSegment &)
Definition: MuonTrackSteering.cxx:28
readCCLHist.float
float
Definition: readCCLHist.py:83
Generated on Fri May 17 2024 21:14:08 for ATLAS Offline Software by   doxygen 1.8.18