ATLAS Offline Software
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 
4 import ROOT
5 import sys, os
6 import logging
7 logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.INFO)
8 import argparse
9 parser = argparse.ArgumentParser()
10 parser.add_argument('infile', type=str, help='input HIST file')
11 parser.add_argument('--grl', type=str, help='Specify an input GRL')
12 parser.add_argument('--out', type=str, help='output ROOT file')
13 parser.add_argument('--tag', type=str, help='Lumi tag',
14  default='OflLumiAcct-001')
15 parser.add_argument('--useofficial', action='store_true', help='Use official lumi folder (otherwise, use OflLumiAcct')
16 parser.add_argument('--lumifolder', type=str, help='Lumi folder', default='/TRIGGER/OFLLUMI/OflPrefLumi')
17 parser.add_argument('--lumitag', type=str, help='Lumi tag', default='OflLumi-13TeV-009')
18 parser.add_argument('--plotdir', type=str, help='Directory to dump plots',
19  default='plots')
20 parser.add_argument('--mudep', type=int, help='Run mu-dependent efficiencies',
21  default=0)
22 parser.add_argument('--dblivetime', action='store_true',
23  help='Look up livetime from DB')
24 parser.add_argument('--mode', type=str, help='Zee or Zmumu')
25 
26 args = parser.parse_args()
27 
28 BINWIDTH=10
29 
30 ZPURITYFACTOR=0.9935
31 ZXSEC=1.929
32 #ZATIMESC=0.2578
33 ZATIMESC=0.29632
34 
35 def 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 
44 ROOT.gStyle.SetOptStat(0)
45 
46 fin = ROOT.TFile.Open(args.infile)
47 runname = None
48 for key in fin.GetListOfKeys():
49  if key.GetName().startswith('run_'):
50  runname = key.GetName()
51  break
52 
53 if args.grl:
54  import DQUtils
55  grl = DQUtils.grl.load_grl(args.grl)
56 else:
57  grl = None
58 
59 if not runname:
60  logging.critical("Can't find run_* directory in input file %s", args.infile)
61  sys.exit(1)
62 
63 z_m = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_mu' % runname)
64 if args.out:
65  outfname = args.out
66 else:
67  outfname = '%s_data.root' % runname[4:]
68 
69 runmode = args.mode
70 print('Running in', runmode, 'mode')
71 if 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 
77 if 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 
84 fout = None
85 t = None
86 
87 from array import array
88 o_passgrl = array('i', [0])
89 o_mu = array('f', [0.])
90 
91 o_lb_rb = array('I', [0,0])
92 o_lbwhen_rb = array('d', [0., 0.])
93 o_zlumi_rb = array('f', [0.])
94 o_zlumistat_rb = array('f', [0.])
95 o_offlumi_rb = array('f', [0.])
96 o_mu_rb = array('f', [0.])
97 o_lblive_rb = array('f', [0.])
98 if 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 
137 lb_length = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_lblength_lb' % runname)
138 lbmin, lbmax = lb_length.GetXaxis().GetXmin(), lb_length.GetXaxis().GetXmax()
139 logging.info('low, high LBs: %s, %s', lbmin, lbmax)
140 
141 if 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)
144 else:
145  livetime = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_livetime_lb' % runname)
146 
147 official_lum = ROOT.TProfile('official_lum', 'official integrated luminosity', int(lbmax-lbmin), lbmin, lbmax)
148 official_lum_zero = ROOT.TProfile('official_lum_zero', 'official inst luminosity', int(lbmax-lbmin), lbmin, lbmax)
149 official_mu = ROOT.TProfile('official_mu', 'official mu', int(lbmax-lbmin), lbmin, lbmax)
150 from DQUtils import fetch_iovs
151 from DQUtils.iov_arrangement import inverse_lblb
152 lblb = fetch_iovs("LBLB", runs=int(runname[4:]))
153 lbtime = inverse_lblb(lblb)
154 #print(list(lbtime))
155 iovs_acct = fetch_iovs('COOLOFL_TRIGGER::/TRIGGER/OFLLUMI/LumiAccounting', lbtime.first.since, lbtime.last.until, tag=args.tag)
156 if 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))
159 lb_start_end = {}
160 lb_lhcfill = {}
161 for iov in lblb:
162  lb_start_end[iov.since & 0xffffffff] = (iov.StartTime/1e9, iov.EndTime/1e9)
163 
164 for 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 
186 divisor = lb_length.Clone('divisor').ProjectionX()
187 px = livetime.ProjectionX()
188 divisor.Multiply(px)
189 
190 nrebinned_bins = ((lbmax-lbmin) // BINWIDTH) + 1
191 
192 if runmode == 'Zee':
193  lumititle = 'Lumi, Z->ee (Run %s)' % runname[4:]
194  efftitle = 'eff #sigma, Z->ee'
195  lumirawtitle = 'Lumi, Z->ee per LB'
196 if 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 
202 lumiplot_m = ROOT.TH1F('lumiplot_m', lumititle % runname[4:],
203  int(nrebinned_bins),
204  lbmin, lbmin+BINWIDTH*nrebinned_bins)
205 lumiplot_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)
208 lumiplot_m.SetXTitle('LB')
209 lumiplot_m.SetYTitle('Luminosity (x 10^{33} cm^{-2} s^{-1})')
210 xsec_m = ROOT.TH1F('xsec_m', efftitle, int(nrebinned_bins),
211  lbmin, lbmin+BINWIDTH*nrebinned_bins)
212 lumiplot_raw_m = ROOT.TH1F('lumiplot_raw_m', lumirawtitle,
213  int(lbmax-lbmin),
214  lbmin, lbmax)
215 
216 num_m, lum, denom, weighted_mu = 0, 0, 0, 0
217 tot_num_m, tot_denom, tot_lum = 0, 0, 0
218 for 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 
295 if 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 
308 c1 = ROOT.TCanvas()
309 c1.SetTickx()
310 c1.SetTicky()
311 leg = ROOT.TLegend(0.6, 0.75, 0.89, 0.88)
312 lumiplot_m.Draw()
313 official_lum_zero.SetLineColor(ROOT.kRed)
314 official_lum_zero.Draw('SAME,HIST')
315 leg.AddEntry(lumiplot_m, 'Z luminosity')
316 leg.AddEntry(official_lum_zero, 'ATLAS preferred lumi', 'L')
317 leg.SetBorderSize(0)
318 leg.Draw()
319 c1.Print(os.path.join(args.plotdir, '%s_lumi.eps' % runname[4:]))
320 c1.Print(os.path.join(args.plotdir, '%s_lumi.png' % runname[4:]))
321 
322 c1.Clear()
323 lumiplot_m_ratio.Draw()
324 c1.Print(os.path.join(args.plotdir, '%s_lumi_ratio.eps' % runname[4:]))
325 c1.Print(os.path.join(args.plotdir, '%s_lumi_ratio.png' % runname[4:]))
xrange
void xrange(TH1 *h, bool symmetric)
Definition: computils.cxx:515
python.iov_arrangement.inverse_lblb
def inverse_lblb(iovs)
Definition: iov_arrangement.py:34
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
dqt_zlumi_compute_lumi.mu_dep_eff
def mu_dep_eff(mu)
Definition: dqt_zlumi_compute_lumi.py:35
array
dqt_zlumi_compute_lumi.int
int
Definition: dqt_zlumi_compute_lumi.py:20
dbg::print
void print(std::FILE *stream, std::format_string< Args... > fmt, Args &&... args)
Definition: SGImplSvc.cxx:70
Generated on Thu Nov 7 2024 21:13:49 for ATLAS Offline Software by   doxygen 1.8.18