ATLAS Offline Software
Loading...
Searching...
No Matches
dqt_zlumi_compute_lumi.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
4import ROOT
5import sys, os
6import logging
7logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.INFO)
8import argparse
9parser = argparse.ArgumentParser()
10parser.add_argument('infile', type=str, help='input HIST file')
11parser.add_argument('--grl', type=str, help='Specify an input GRL')
12parser.add_argument('--out', type=str, help='output ROOT file')
13parser.add_argument('--tag', type=str, help='Lumi tag',
14 default='OflLumiAcct-001')
15parser.add_argument('--useofficial', action='store_true', help='Use official lumi folder (otherwise, use OflLumiAcct')
16parser.add_argument('--lumifolder', type=str, help='Lumi folder', default='/TRIGGER/OFLLUMI/OflPrefLumi')
17parser.add_argument('--lumitag', type=str, help='Lumi tag', default='OflLumi-13TeV-009')
18parser.add_argument('--plotdir', type=str, help='Directory to dump plots',
19 default='plots')
20parser.add_argument('--mudep', type=int, help='Run mu-dependent efficiencies',
21 default=0)
22parser.add_argument('--dblivetime', action='store_true',
23 help='Look up livetime from DB')
24parser.add_argument('--mode', type=str, help='Zee or Zmumu')
25
26args = parser.parse_args()
27
28BINWIDTH=10
29
30ZPURITYFACTOR=0.9935
31ZXSEC=1.929
32#ZATIMESC=0.2578
33ZATIMESC=0.29632
34
35def mu_dep_eff(mu):
36 #make breakpoint at 8 match
37 if 0 <= mu < 8: return 0.3152 #0.3141
38 elif 8 <= mu < 27: return 0.3191 - 0.000493*mu
39 elif 27 <= mu: return 0.3443 - 0.00143*mu
40 else:
41 print('WTF??')
42 return ZATIMESC
43
44ROOT.gStyle.SetOptStat(0)
45
46fin = ROOT.TFile.Open(args.infile)
47runname = None
48for key in fin.GetListOfKeys():
49 if key.GetName().startswith('run_'):
50 runname = key.GetName()
51 break
52
53if args.grl:
54 import DQUtils
55 grl = DQUtils.grl.load_grl(args.grl)
56else:
57 grl = None
58
59if not runname:
60 logging.critical("Can't find run_* directory in input file %s", args.infile)
61 sys.exit(1)
62
63z_m = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_mu' % runname)
64if args.out:
65 outfname = args.out
66else:
67 outfname = '%s_data.root' % runname[4:]
68
69runmode = args.mode
70print('Running in', runmode, 'mode')
71if runmode == 'Zee':
72 z_m = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_el' % runname)
73 if not z_m:
74 logging.critical("Can't retrieve m_Z_Counter_el")
75 sys.exit(1)
76
77if runmode == 'Zmumu':
78 z_m = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_mu' % runname)
79 if not z_m:
80 logging.critical("Can't retrieve m_Z_Counter_mu")
81 sys.exit(1)
82
83
84fout = None
85t = None
86
87from array import array
88o_passgrl = array('i', [0])
89o_mu = array('f', [0.])
90
91o_lb_rb = array('I', [0,0])
92o_lbwhen_rb = array('d', [0., 0.])
93o_zlumi_rb = array('f', [0.])
94o_zlumistat_rb = array('f', [0.])
95o_offlumi_rb = array('f', [0.])
96o_mu_rb = array('f', [0.])
97o_lblive_rb = array('f', [0.])
98if True:
99 fout = ROOT.TFile.Open(outfname, 'RECREATE')
100 o_run = array('I', [int(runname[4:])])
101 o_lb = array('I', [0])
102 o_lbwhen = array('d', [0., 0.])
103 o_zraw = array('f', [0.])
104 o_zrawstat = array('f', [0.])
105 o_zlumi = array('f', [0.])
106 o_zlumistat = array('f', [0.])
107 o_offlumi = array('f', [0.])
108 o_lblive = array('f', [0.])
109 o_lhcfill = array('I', [0])
110 t = ROOT.TTree( 'lumitree', 'Luminosity tree' )
111 t.Branch('run', o_run, 'run/i')
112 t.Branch('lb', o_lb, 'lb/i')
113 t.Branch('lbwhen', o_lbwhen, 'lbwhen[2]/D')
114 t.Branch('zraw', o_zraw, 'zraw/F')
115 t.Branch('zrawstat', o_zrawstat, 'zrawstat/F')
116 t.Branch('zlumi', o_zlumi, 'zlumi/F')
117 t.Branch('zlumistat', o_zlumistat, 'zlumistat/F')
118 t.Branch('offlumi', o_offlumi, 'offlumi/F')
119 t.Branch('mu', o_mu, 'mu/F')
120 t.Branch('lblive', o_lblive, 'lblive/F')
121 t.Branch('pass_grl', o_passgrl, 'pass_grl/I')
122 t.Branch('lhcfill', o_lhcfill, 'lhcfill/i')
123
124 t_rb = ROOT.TTree( 'lumitree_rb', 'Luminosity tree, rebinned' )
125 t_rb.Branch('run', o_run, 'run/i')
126 t_rb.Branch('lb', o_lb_rb, 'lb[2]/i')
127 t_rb.Branch('lbwhen', o_lbwhen_rb, 'lbwhen[2]/D')
128 t_rb.Branch('zlumi', o_zlumi_rb, 'zlumi/F')
129 t_rb.Branch('zlumistat', o_zlumistat_rb, 'zlumistat/F')
130 t_rb.Branch('offlumi', o_offlumi_rb, 'offlumi/F')
131 t_rb.Branch('mu', o_mu_rb, 'mu/F')
132 t_rb.Branch('lblive', o_lblive_rb, 'lblive/F')
133
134 ROOT.gROOT.cd('/')
135
136
137lb_length = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_lblength_lb' % runname)
138lbmin, lbmax = lb_length.GetXaxis().GetXmin(), lb_length.GetXaxis().GetXmax()
139logging.info('low, high LBs: %s, %s', lbmin, lbmax)
140
141if args.dblivetime:
142 logging.info('Starting livetime lookup ... (remove when we have a proper in-file implementation ...)')
143 livetime = ROOT.TProfile('livetime', 'Livetime', int(lbmax-lbmin), lbmin, lbmax)
144else:
145 livetime = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_livetime_lb' % runname)
146
147official_lum = ROOT.TProfile('official_lum', 'official integrated luminosity', int(lbmax-lbmin), lbmin, lbmax)
148official_lum_zero = ROOT.TProfile('official_lum_zero', 'official inst luminosity', int(lbmax-lbmin), lbmin, lbmax)
149official_mu = ROOT.TProfile('official_mu', 'official mu', int(lbmax-lbmin), lbmin, lbmax)
150from DQUtils import fetch_iovs
151from DQUtils.iov_arrangement import inverse_lblb
152lblb = fetch_iovs("LBLB", runs=int(runname[4:]))
153lbtime = inverse_lblb(lblb)
154#print(list(lbtime))
155iovs_acct = fetch_iovs('COOLOFL_TRIGGER::/TRIGGER/OFLLUMI/LumiAccounting', lbtime.first.since, lbtime.last.until, tag=args.tag)
156if args.useofficial:
157 iovs_lum = fetch_iovs('COOLOFL_TRIGGER::%s' % args.lumifolder, lblb.first.since, lblb.last.until, tag=args.lumitag, channels=[0])
158 #print(list(iovs_lum))
159lb_start_end = {}
160lb_lhcfill = {}
161for iov in lblb:
162 lb_start_end[iov.since & 0xffffffff] = (iov.StartTime/1e9, iov.EndTime/1e9)
163
164for iov in iovs_acct:
165 if not lbmin < iov.LumiBlock < lbmax:
166 continue
167 lb_lhcfill[iov.LumiBlock] = iov.FillNumber
168 if args.dblivetime:
169 livetime.Fill(iov.LumiBlock, iov.LiveFraction)
170 #print(iov.InstLumi, iovs_lum[iov.LumiBlock-1].LBAvInstLumi)
171 if not args.useofficial:
172 official_lum_zero.Fill(iov.LumiBlock, iov.InstLumi/1e3)
173 official_lum.Fill(iov.LumiBlock, iov.InstLumi*iov.LBTime*iov.LiveFraction/1e3)
174 official_mu.Fill(iov.LumiBlock, iov.AvEvtsPerBX)
175 else:
176 offlumiov = [_ for _ in iovs_lum if _.since.lumi==iov.LumiBlock]
177 if len(offlumiov) != 1:
178 print('MAJOR PROBLEM, LUMI IOV MISMATCH')
179 print(len(offlumiov))
180 continue
181 offlumiov = offlumiov[0]
182 official_lum_zero.Fill(iov.LumiBlock, offlumiov.LBAvInstLumi/1e3)
183 official_lum.Fill(iov.LumiBlock, offlumiov.LBAvInstLumi*iov.LBTime*iov.LiveFraction/1e3)
184 official_mu.Fill(iov.LumiBlock, offlumiov.LBAvEvtsPerBX)
185
186divisor = lb_length.Clone('divisor').ProjectionX()
187px = livetime.ProjectionX()
188divisor.Multiply(px)
189
190nrebinned_bins = ((lbmax-lbmin) // BINWIDTH) + 1
191
192if runmode == 'Zee':
193 lumititle = 'Lumi, Z->ee (Run %s)' % runname[4:]
194 efftitle = 'eff #sigma, Z->ee'
195 lumirawtitle = 'Lumi, Z->ee per LB'
196if runmode == 'Zmumu':
197 lumititle = 'Lumi, Z->#mu#mu (Run %s)' % runname[4:]
198 efftitle = 'eff #sigma, Z->#mu#mu'
199 lumirawtitle = 'Lumi, Z->#mu#mu per LB'
200
201
202lumiplot_m = ROOT.TH1F('lumiplot_m', lumititle % runname[4:],
203 int(nrebinned_bins),
204 lbmin, lbmin+BINWIDTH*nrebinned_bins)
205lumiplot_m_ratio = ROOT.TH1F('lumiplot_m_ratio', 'Z/official lumi ratio (Run %s)' % runname[4:],
206 int(nrebinned_bins),
207 lbmin, lbmin+BINWIDTH*nrebinned_bins)
208lumiplot_m.SetXTitle('LB')
209lumiplot_m.SetYTitle('Luminosity (x 10^{33} cm^{-2} s^{-1})')
210xsec_m = ROOT.TH1F('xsec_m', efftitle, int(nrebinned_bins),
211 lbmin, lbmin+BINWIDTH*nrebinned_bins)
212lumiplot_raw_m = ROOT.TH1F('lumiplot_raw_m', lumirawtitle,
213 int(lbmax-lbmin),
214 lbmin, lbmax)
215
216num_m, lum, denom, weighted_mu = 0, 0, 0, 0
217tot_num_m, tot_denom, tot_lum = 0, 0, 0
218for ibin in xrange(1, int(lbmax-lbmin)+1):
219 profileflag=True
220 try:
221 z_m[ibin]
222 except IndexError as e:
223 logging.error('Something unfortunate has happened; LB %d missing from Z count' % (ibin + lbmin - 1))
224 profileflag=False
225 if args.mudep:
226 l_zatimesc = mu_dep_eff(official_mu[ibin])
227 else:
228 l_zatimesc = ZATIMESC
229 if grl and not DQUtils.grl.grl_contains_run_lb(grl, (int(runname[4:]), int(lumiplot_raw_m.GetBinCenter(ibin)))):
230 o_passgrl[0]=0
231 else:
232 o_passgrl[0]=1
233 if divisor[ibin] > 0 and profileflag:
234 lumiplot_raw_m.SetBinContent(ibin, z_m[ibin]/divisor[ibin]*ZPURITYFACTOR/l_zatimesc/ZXSEC)
235 lumiplot_raw_m.SetBinError(ibin, z_m[ibin]**.5/divisor[ibin]*ZPURITYFACTOR/l_zatimesc/ZXSEC)
236 o_mu[0] = official_mu[ibin]
237 if o_passgrl[0]:
238 if profileflag:
239 num_m += z_m[ibin]; tot_num_m += z_m[ibin]
240 denom += divisor[ibin]; tot_denom += divisor[ibin]
241 lum += official_lum[ibin]; tot_lum += official_lum[ibin]
242 weighted_mu += o_mu[0]*divisor[ibin]
243
244
245 # fill tree
246 if t:
247 #print(ibin, lumiplot_raw_m.GetBinCenter(ibin))
248 o_lb[0] = int(lumiplot_raw_m.GetBinCenter(ibin))
249 o_lbwhen[0] = lb_start_end[o_lb[0]][0]
250 o_lbwhen[1] = lb_start_end[o_lb[0]][1]
251 o_zraw[0] = z_m[ibin] if profileflag else 0
252 o_zrawstat[0] = z_m.GetBinError(ibin) if profileflag else 0
253 o_zlumi[0] = lumiplot_raw_m[ibin]
254 o_zlumistat[0] = lumiplot_raw_m.GetBinError(ibin)
255 o_offlumi[0] = official_lum_zero[ibin]
256 o_lblive[0] = divisor[ibin]
257 o_lhcfill[0] = lb_lhcfill[o_lb[0]]
258 t.Fill()
259
260 if (ibin % BINWIDTH) == 0:
261 ribin = int(ibin // BINWIDTH)
262 o_mu_rb[0] = weighted_mu/denom if denom > 0 else 0
263 if args.mudep:
264 l_zatimesc = mu_dep_eff(o_mu_rb[0])
265 else:
266 l_zatimesc = ZATIMESC
267 if denom > 0:
268 lumiplot_m.SetBinContent(ribin, num_m/denom*ZPURITYFACTOR/l_zatimesc/ZXSEC)
269 lumiplot_m.SetBinError(ribin, num_m**.5/denom*ZPURITYFACTOR/l_zatimesc/ZXSEC)
270 if lum > 0:
271 xsec_m.SetBinContent(ribin, num_m/lum*ZPURITYFACTOR/l_zatimesc)
272 xsec_m.SetBinError(ribin, num_m**.5/lum*ZPURITYFACTOR/l_zatimesc)
273 if denom > 0:
274 o_zlumi_rb[0] = num_m/denom*ZPURITYFACTOR/l_zatimesc/ZXSEC
275 o_zlumistat_rb[0] = num_m**.5/denom*ZPURITYFACTOR/l_zatimesc/ZXSEC
276 o_offlumi_rb[0] = lum/denom
277 if o_offlumi_rb[0] > 0:
278 lumiplot_m_ratio.SetBinContent(ribin, o_zlumi_rb[0]/o_offlumi_rb[0])
279 lumiplot_m_ratio.SetBinError(ribin, o_zlumistat_rb[0]/o_offlumi_rb[0])
280 else:
281 o_zlumi_rb[0] = 0.
282 o_zlumistat_rb[0] = 0.
283 o_offlumi_rb[0] = 0.
284 o_mu_rb[0] = 0.
285 o_lb_rb[1] = int(lumiplot_raw_m.GetBinCenter(ibin))
286 o_lb_rb[0] = int(lumiplot_raw_m.GetBinCenter(ibin)-BINWIDTH+1)
287 o_lbwhen_rb[0] = lb_start_end[o_lb_rb[0]][0]
288 o_lbwhen_rb[1] = lb_start_end[o_lb_rb[1]][1]
289 o_lblive_rb[0] = denom
290 if t:
291 t_rb.Fill()
292
293 num_m, lum, denom, weighted_mu = 0, 0, 0, 0
294
295if fout:
296 fout.cd()
297 t.Write()
298 t_rb.Write()
299
300 if tot_lum > 0:
301 tge = ROOT.TGraphErrors(1)
302 tge.SetPoint(0, int(runname[4:]), tot_num_m*ZPURITYFACTOR/ZATIMESC/ZXSEC/tot_lum)
303 tge.SetPointError(0, 0, tot_num_m**.5*ZPURITYFACTOR/ZATIMESC/ZXSEC/tot_lum)
304 tge.SetName('lum_ratio')
305 tge.Write()
306 fout.Close()
307
308c1 = ROOT.TCanvas()
309c1.SetTickx()
310c1.SetTicky()
311leg = ROOT.TLegend(0.6, 0.75, 0.89, 0.88)
312lumiplot_m.Draw()
313official_lum_zero.SetLineColor(ROOT.kRed)
314official_lum_zero.Draw('SAME,HIST')
315leg.AddEntry(lumiplot_m, 'Z luminosity')
316leg.AddEntry(official_lum_zero, 'ATLAS preferred lumi', 'L')
317leg.SetBorderSize(0)
318leg.Draw()
319c1.Print(os.path.join(args.plotdir, '%s_lumi.eps' % runname[4:]))
320c1.Print(os.path.join(args.plotdir, '%s_lumi.png' % runname[4:]))
321
322c1.Clear()
323lumiplot_m_ratio.Draw()
324c1.Print(os.path.join(args.plotdir, '%s_lumi_ratio.eps' % runname[4:]))
325c1.Print(os.path.join(args.plotdir, '%s_lumi_ratio.png' % runname[4:]))
print
void print(char *figname, TCanvas *c1)
Definition TRTCalib_StrawStatusPlots.cxx:26
array
STL class.
xrange
void xrange(TH1 *h, bool symmetric)
Definition computils.cxx:517
dqt_zlumi_compute_lumi.mu_dep_eff
mu_dep_eff(mu)
Definition dqt_zlumi_compute_lumi.py:35
Generated on for ATLAS Offline Software by doxygen 1.14.0