ATLAS Offline Software
Loading...
Searching...
No Matches
dqt_zlumi_pandas.py
Go to the documentation of this file.
1#!/usr/bin/env python
2# Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3
4import numpy as np
5import csv
6from array import array
7import ROOT
8import sys, os
9import logging
10import argparse
11import math
12from DQUtils import fetch_iovs
13from DQUtils.iov_arrangement import inverse_lblb
14import ZLumiScripts.tools.zlumi_alleff as dq_eff
15import ZLumiScripts.tools.zlumi_mc_cf as dq_cf
16
17
18from DataQualityUtils import doZLumi
19
20# testing toy sampling
21import ZLumiScripts.tools.toys as toys
22
23ROOT.gROOT.SetBatch(ROOT.kTRUE)
24ROOT.gStyle.SetOptStat(0)
25
26logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.INFO)
27
28parser = argparse.ArgumentParser()
29parser.add_argument('--infile', type=str, help='input HIST file')
30parser.add_argument('--grl', type=str, help='Specify an input GRL')
31parser.add_argument('--tag', type=str, help='Lumi tag', default='OflLumiAcct-Run3-003')
32parser.add_argument('--useofficial', action='store_true', help='Use official lumi folder (otherwise, use OflLumiAcct')
33parser.add_argument('--lumifolder', type=str, help='Lumi folder', default='/TRIGGER/OFLLUMI/OflPrefLumi')
34parser.add_argument('--lumitag', type=str, help='Lumi tag', default='OflLumi-Run3-003')
35parser.add_argument('--outdir', type=str, help='Directory to dump plots', default='plots')
36parser.add_argument('--dblivetime', action='store_true', help='Look up livetime from DB')
37parser.add_argument('--campaign', type=str, help='mc16a/d/e, mc21, mc23a/d/e')
38
39args = parser.parse_args()
40campaign = args.campaign
41run3 = "mc2" in campaign
42
43ZPURITYFACTOR = 0.9935
44if run3:
45 ZXSEC = 2.0675
46else:
47 ZXSEC = 1.929
48
49zee_missing_lbs = []
50zmumu_missing_lbs = []
51
52ntoys = 10000000
53do_toys = False
54
55fin = ROOT.TFile.Open(args.infile)
56runname = None
57for key in fin.GetListOfKeys():
58 if key.GetName().startswith('run_'):
59 runname = key.GetName()
60 runnumber = int(runname.replace('run_',''))
61 break
62
63if not runname:
64 logging.critical("Can't find run_* directory in input file %s", args.infile)
65 sys.exit(1)
66
67print("Starting HIST to CSV conversion for Run ", runnumber)
68
69if args.outdir:
70 out_dir = args.outdir
71 os.system("mkdir -p " + out_dir)
72 out_dir += "/" + runname + ".csv"
73else:
74 out_dir = runname + ".csv"
75
76if run3:
77 lb_length_name = '%s/GLOBAL/DQTGlobalWZFinder/duration_vs_LB' % runname
78 livetime_name = '%s/GLOBAL/DQTGlobalWZFinder/avgLiveFrac_vs_LB' % runname
79else:
80 lb_length_name = '%s/GLOBAL/DQTGlobalWZFinder/m_lblength_lb' % runname
81 livetime_name = '%s/GLOBAL/DQTGlobalWZFinder/m_livetime_lb' % runname
82
83
84if args.grl:
85 grlReader = ROOT.Root.TGoodRunsListReader(args.grl)
86 grlReader.Interpret()
87 grl = grlReader.GetMergedGRLCollection()
88else:
89 grlname = 'grl_'+str(runnumber)+'.xml'
90 grl_file = doZLumi.makeGRL(runnumber, 'PHYS_StandardGRL_All_Good', grlname)
91
92lb_length_old = fin.Get(lb_length_name)
93lbmin, lbmax = lb_length_old.FindFirstBinAbove(0, 1, 0, -1), lb_length_old.FindLastBinAbove(0, 1, 0, -1)
94lb_length = ROOT.TProfile('lb_length', 'LB length', int(lbmax-lbmin), lbmin, lbmax)
95
96for i in range(lbmin, lbmax):
97 lb_length.Fill(i, lb_length_old.GetBinContent(i))
98
99logging.info('low, high LBs: %s, %s', lbmin, lbmax)
100
101if args.dblivetime:
102 logging.info('Starting livetime lookup ... (remove when we have a proper in-file implementation ...)')
103 livetime = ROOT.TProfile('livetime', 'Livetime', int(lbmax-lbmin), lbmin, lbmax)
104else:
105 livetime = fin.Get(livetime_name)
106
107official_lum_zero = ROOT.TProfile('official_lum_zero', 'official inst luminosity', int(lbmax-lbmin), lbmin, lbmax)
108official_mu = ROOT.TProfile('official_mu', 'official mu', int(lbmax-lbmin), lbmin, lbmax)
109
110lblb = fetch_iovs("LBLB", runs=runnumber)
111lbtime = inverse_lblb(lblb)
112iovs_acct = fetch_iovs('COOLOFL_TRIGGER::/TRIGGER/OFLLUMI/LumiAccounting', lbtime.first.since, lbtime.last.until, tag=args.tag)
113
114if args.useofficial:
115 print("Using official lumitag", args.lumitag)
116 iovs_lum = fetch_iovs('COOLOFL_TRIGGER::%s' % args.lumifolder, lblb.first.since, lblb.last.until, tag=args.lumitag, channels=[0])
117
118lb_start_end = {}
119lb_lhcfill = {}
120for iov in lblb:
121 lb_start_end[iov.since & 0xffffffff] = (iov.StartTime/1e9, iov.EndTime/1e9)
122
123for iov in iovs_acct:
124 if not lbmin-1 < iov.LumiBlock < lbmax:
125 continue
126 lb_lhcfill[iov.LumiBlock] = iov.FillNumber
127 if args.dblivetime:
128 livetime.Fill(iov.LumiBlock, iov.LiveFraction)
129
130 if not args.useofficial:
131 official_lum_zero.Fill(iov.LumiBlock, iov.InstLumi/1e3)
132 official_mu.Fill(iov.LumiBlock, iov.AvEvtsPerBX)
133 else:
134 offlumiov = [_ for _ in iovs_lum if _.since.lumi==iov.LumiBlock]
135 if len(offlumiov) != 1:
136 print('MAJOR PROBLEM, LUMI IOV MISMATCH', iov.LumiBlock)
137 continue
138
139 offlumiov = offlumiov[0]
140 official_lum_zero.Fill(iov.LumiBlock, offlumiov.LBAvInstLumi/1e3)
141 official_mu.Fill(iov.LumiBlock, offlumiov.LBAvEvtsPerBX)
142
143lb_full = lb_length.Clone('lb_full').ProjectionX()
144divisor = lb_length.Clone('divisor').ProjectionX()
145px = livetime.ProjectionX()
146divisor.Multiply(px)
147
148# Get run-wise electron channel histos outside of loop
149# Also do the container efficiency bkg fit here as it is done
150# once per run, then normalised to the ratio of LB/run in the function
151# container_efficiency
152hto = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_ele_template_os' % (runname))
153hts = fin.Get('%s/GLOBAL/DQTGlobalWZFinder/m_ele_template_ss' % (runname))
154
155# ==== Set signal region == 0, then fit ====
156# hphotontotal = fin.Get("%s/GLOBAL/DQTGlobalWZFinder/m_elContainertp_nomatch" % (runname))
157# hphotontotal.GetXaxis().SetRangeUser(66000, 250000)
158# h_fit = hphotontotal.Clone()
159# for xbin in range(1, h_fit.GetNbinsX()+1):
160# mass = h_fit.GetBinLowEdge(xbin)
161# if mass > 75000 and mass < 100000:
162# h_fit.SetBinContent(xbin, 0)
163# h_fit.SetBinError(xbin, 0)
164# h_fit.Fit("pol2", "q")
165
166
167o_recoeff_fit = {}
168o_recoerr_fit = {}
169# Do electron channel reco eff calculation once,
170# then apply iterative correction to bkg subtraction in the main loop
171lb_minus_one_reco_eff = [1.0, 1.0, 1.0]
172for ibin in range(1, int(lbmax-lbmin)+1):
173 this_lb = int(lb_full.GetBinCenter(ibin))
174 lb = "lb_" + str(this_lb)
175 pileup = round(official_mu[ibin])
176
177 hmo = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
178 hms = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
179 hno = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
180 hns = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
181
182 try:
183 eff, err = dq_eff.template_method(hmo, hms, hno, hns, hto, hts, False, lb_minus_one_reco_eff) # do_scale set to false here as we are calculating it, when applying set True (next loop)
184 lb_minus_one_reco_eff = [eff, err, this_lb]
185 except AttributeError:
186 eff = 0
187 err = 0
188
189 if err != 0 and eff != 0:
190 weight = 1/pow(err, 2)
191 if pileup not in o_recoeff_fit:
192 o_recoeff_fit[pileup] = weight * eff
193 o_recoerr_fit[pileup] = weight
194 else:
195 o_recoeff_fit[pileup] += weight * eff
196 o_recoerr_fit[pileup] += weight
197
198arr_rec_eff = array('d')
199arr_rec_err = array('d')
200arr_mu = array('d')
201
202for pileup in o_recoeff_fit:
203 arr_mu.append(float(pileup))
204 arr_rec_eff.append(o_recoeff_fit[pileup]/o_recoerr_fit[pileup])
205 arr_rec_err.append(1/pow(o_recoerr_fit[pileup], 0.5))
206
207# If no pileup data is available, add dummy entry with all parameters set to 1
208
209if len(arr_mu) == 0:
210 arr_mu.append(1)
211 arr_rec_eff.append(1)
212 arr_rec_err.append(1)
213
214tg_fit = ROOT.TGraphErrors(len(arr_mu), arr_mu, arr_rec_eff, ROOT.nullptr, arr_rec_err)
215if len(o_recoeff_fit) == 0 or len(o_recoeff_fit) == 1:
216 fit_type = "pol0"
217elif len(o_recoeff_fit) == 2:
218 fit_type = "pol1"
219elif len(o_recoeff_fit) > 2:
220 fit_type = "pol2"
221
222tg_fit.Fit(fit_type, "q")
223
224csvfile = open(out_dir, 'w')
225csvwriter = csv.writer(csvfile, delimiter=',')
226csvwriter.writerow(['FillNum','RunNum','LBNum','LBStart','LBEnd','LBLive','LBFull','OffLumi','OffMu', 'PassGRL',
227 'ZeeRaw','ZeeRawErr','ZeeN1','ZeeN2','ZeeEffTrig','ZeeErrTrig','ZeeEffReco','ZeeErrReco','ZeeEffComb','ZeeErrComb','ZeeEffAComb','ZeeErrAComb','ZeeDefTrig','ZeeDefReco','ZeeLumi','ZeeLumiErr','ZeeRate',
228 'ZmumuRaw','ZmumuRawErr','ZmumuN1','ZmumuN2','ZmumuEffTrig','ZmumuErrTrig','ZmumuEffReco','ZmumuErrReco','ZmumuEffComb','ZmumuErrComb','ZmumuEffAComb','ZmumuErrAComb','ZmumuDefTrig','ZmumuDefReco','ZmumuLumi','ZmumuLumiErr','ZmumuRate',
229 'ZllLumi','ZllLumiErr'])
230
231lb_minus_one_reco_eff = {}
232lb_minus_one_reco_eff["Zee"] = [1.0, 1.0, 1]
233lb_minus_one_reco_eff["Zmumu"] = [1.0, 1.0, 1]
234
235lb_minus_one_trig_eff = {}
236lb_minus_one_trig_eff["Zee"] = [1.0, 1.0, 1]
237lb_minus_one_trig_eff["Zmumu"] = [1.0, 1.0, 1]
238
239bad_database = False
240
241for ibin in range(1, int(lbmax-lbmin)+1):
242 out_dict = {}
243 out_dict["Zee"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
244 out_dict["Zmumu"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
245
246 this_lb = int(lb_full.GetBinCenter(ibin))
247 loclivetime = divisor[ibin]
248 # some ad-hoc Run 2 corrections
249 if runnumber == 302831 and this_lb < 11:
250 loclivetime = 0
251 elif runnumber == 329835 and this_lb < 554:
252 loclivetime = 0
253 elif runnumber == 310247 and (this_lb == 442 or this_lb == 462):
254 loclivetime = 0
255 elif runnumber == 281385 and this_lb < 197:
256 loclivetime *= 4.0/6.0
257 elif runnumber == 281385 and this_lb < 375:
258 loclivetime *= 5.0/6.0
259 elif runnumber == 286367:
260 loclivetime *= 5.0/6.0
261
262 try:
263 this_fill = lb_lhcfill[this_lb]
264 except KeyError:
265 bad_database = True
266 this_fill = "NA"
267
268 passgrl = 1
269 for channel in ["Zee", "Zmumu"]:
270 if grl and not grl.HasRunLumiBlock(runnumber, this_lb):
271 passgrl = 0
272 continue
273 # at less than 9s LB livetime: all Z counting zeroed
274 if loclivetime < 9:
275 continue
276
277 lb = "lb_" + str(this_lb)
278 if channel == 'Zee':
279 # Nominal signal histogram
280 hname = runname + '/lb_'+str(int(ibin+lbmin-0.5))+'/GLOBAL/DQTGlobalWZFinder/m_Z_mass_opsele'
281
282 # Tag-and-probe histos
283 h = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_eltrigtp_matches_os' % (runname, lb))
284 hmo = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
285 hms = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
286 hno = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
287 hns = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
288 #hphoton = fin.Get("%s/%s/GLOBAL/DQTGlobalWZFinder/m_elContainertp_nomatch" % (runname, lb))
289 #hpass = fin.Get("%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_passkine" % (runname, lb))
290 #hphoton.GetXaxis().SetRangeUser(66000, 250000)
291 #hpass.GetXaxis().SetRangeUser(66000, 250000)
292 if run3:
293 ACCEPTANCE = 0.2971
294 else:
295 ACCEPTANCE = 0.2996
296 elif channel == 'Zmumu':
297 # Nominal signal histogram
298 hname = runname + '/lb_'+str(int(ibin+lbmin-0.5))+'/GLOBAL/DQTGlobalWZFinder/m_Z_mass_opsmu'
299
300 # Tag-and-probe histos
301 h = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_mutrigtp_matches' % (runname, lb))
302 hmo = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_os' % (runname, lb))
303 hms = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_ss' % (runname, lb))
304 hno = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_os' % (runname, lb))
305 hns = fin.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_ss' % (runname, lb))
306 if run3:
307 ACCEPTANCE = 0.3292
308 else:
309 ACCEPTANCE = 0.3323224
310
311 try:
312 z_m = fin.Get(hname).Integral()
313 z_merr = math.sqrt(z_m)
314 except AttributeError:
315 if channel == "Zee":
316 zee_missing_lbs.append(this_lb)
317 elif channel == "Zmumu":
318 zmumu_missing_lbs.append(this_lb)
319
320 continue
321
322
323 N1 = N2 = 0
324 try:
325 N1 = h.GetBinContent(2)
326 N2 = h.GetBinContent(3)
327 if do_toys:
328 eff_trig, err_trig, arr_trig, arr_NZ = toys.toy_trigeff(N1, N2, ntoys)
329 else:
330 eff_trig, err_trig = dq_eff.trig_tag_and_probe(h, lb_minus_one_trig_eff[channel])
331 except TypeError:
332 eff_trig, err_trig = 0.0, 0.0
333 except AttributeError:
334 eff_trig, err_trig = 0.0, 0.0
335 if channel == "Zmumu":
336 try:
337 bin1 = hmo.GetXaxis().FindBin(86000)
338 bin2 = hmo.GetXaxis().FindBin(95000)
339 if do_toys:
340 eff_reco, err_reco, arr_reco = toys.muon_toy_recoeff(hmo.Integral(bin1, bin2), hms.Integral(bin1, bin2), hno.Integral(bin1, bin2), hns.Integral(bin1, bin2), ntoys)
341 else:
342 eff_reco, err_reco = dq_eff.reco_tag_and_probe(hmo, hms, hno, hns, lb_minus_one_reco_eff[channel])
343 except TypeError:
344 eff_reco, err_reco = 0.0, 0.0
345 except AttributeError:
346 eff_reco, err_reco = 0.0, 0.0
347 elif channel == "Zee":
348 try:
349 if do_toys:
350 bin1 = hmo.GetXaxis().FindBin(75000)
351 bin2 = hmo.GetXaxis().FindBin(104000)
352 bin3 = hmo.GetXaxis().FindBin(120000)
353 bin4 = hmo.GetXaxis().FindBin(250000)
354 matchos_peak = hmo.Integral(bin1, bin2)
355 matchos_tail = hmo.Integral(bin3, bin4)
356 matchss_tail = hms.Integral(bin3, bin4)
357 nomatchos_peak = hno.Integral(bin1, bin2)
358 nomatchos_tail = hno.Integral(bin3, bin4)
359 templateos_peak = hto.Integral(bin1, bin2)
360 templateos_tail = hto.Integral(bin3, bin4)
361 templatess_tail = hts.Integral(bin3, bin4)
362 eff_reco, err_reco, arr_reco = toys.electron_toy_recoeff(matchos_peak, matchos_tail, matchss_tail, nomatchos_peak, nomatchos_tail, templateos_peak, templateos_tail, templatess_tail, ntoys)
363 else:
364 eff_reco, err_reco = dq_eff.template_method(hmo, hms, hno, hns, hto, hts, True, lb_minus_one_reco_eff[channel], tg_fit.GetFunction(fit_type).Eval(official_mu[ibin]))
365
366 except TypeError:
367 eff_reco, err_reco = 0.0, 0.0
368 except AttributeError:
369 eff_reco, err_reco = 0.0, 0.0
370
371 #dq_eff.container_efficiency(hphoton, hphotontotal, h_fit, hpass, hto, hts)
372
373
374 defaulted_reco_eff = 0
375 defaulted_trig_eff = 0
376 if eff_reco == lb_minus_one_reco_eff[channel][0]:
377 defaulted_reco_eff = 1
378 if eff_trig == lb_minus_one_trig_eff[channel][0]:
379 defaulted_trig_eff = 1
380
381 if eff_reco != 0.0:
382 lb_minus_one_reco_eff[channel] = [eff_reco, err_reco, this_lb]
383 if eff_trig != 0.0:
384 lb_minus_one_trig_eff[channel] = [eff_trig, err_trig, this_lb]
385
386 if do_toys:
387 arr_comb = (1.0 - (1.0 - arr_trig)**2) * (arr_reco)**2
388 nonan_arr_comb = arr_comb[~np.isnan(arr_comb)]
389 eff_comb = np.median(nonan_arr_comb)
390 err_comb = nonan_arr_comb.std()
391 else:
392 eff_comb = (1-(1-eff_trig)**2)*(eff_reco)**2
393 err_comb = ((eff_reco**2*2*(1-eff_trig)*err_trig)**2+(2*eff_reco*(1-(1-eff_trig)**2)*err_reco)**2)**.5
394
395 eff_Acomb = ACCEPTANCE * eff_comb
396 err_Acomb = ACCEPTANCE * err_comb
397
398 if do_toys:
399 effcy = arr_comb * dq_cf.correction(official_mu[ibin], channel, campaign, runnumber)
400 else:
401 effcy = eff_comb * dq_cf.correction(official_mu[ibin], channel, campaign, runnumber)
402 effcyerr = err_comb * dq_cf.correction(official_mu[ibin], channel, campaign, runnumber)
403
404 zlumi = zlumistat = zrate = 0.0
405 CORRECTIONS = ZPURITYFACTOR/ACCEPTANCE/ZXSEC
406
407 if do_toys and loclivetime != 0.0:
408 arr_zlumi = np.divide(arr_NZ, effcy) * (CORRECTIONS)/loclivetime
409 arr_zlumi = arr_zlumi[~np.isnan(arr_zlumi)]
410 zlumi = np.median(arr_zlumi)
411 zlumistat = arr_zlumi.std()
412 zrate = zlumi / CORRECTIONS
413 elif (loclivetime != 0.0 and effcy != 0.0):
414 zlumi = (z_m/effcy)*(CORRECTIONS)/loclivetime
415 zlumistat = math.sqrt(pow(z_merr/effcy, 2) + pow(z_m/effcy**2*effcyerr, 2))*CORRECTIONS/loclivetime
416 zrate = zlumi / CORRECTIONS
417
418
419 out_dict[channel] = [z_m, z_merr, N1, N2, eff_trig, err_trig, eff_reco, err_reco, eff_comb, err_comb, eff_Acomb, err_Acomb, defaulted_trig_eff, defaulted_reco_eff, zlumi, zlumistat, zrate]
420
421 lumi_index = len(out_dict['Zee'])-3
422 error_index = len(out_dict['Zee'])-2
423 zll_lumi = (out_dict['Zee'][lumi_index] + out_dict['Zmumu'][lumi_index])/2
424 zll_lumi_err = 0.5 * math.sqrt( pow(out_dict['Zee'][error_index], 2) + pow(out_dict['Zmumu'][error_index], 2) )
425 try:
426 lb_start, lb_end = lb_start_end[this_lb][0], lb_start_end[this_lb][1]
427 except KeyError:
428 lb_start, lb_end = 0, 0
429 print("WARNING: database not filled completely for LB", this_lb)
430 out_write = [this_fill, runnumber, this_lb, lb_start, lb_end, loclivetime, lb_full[ibin], official_lum_zero[ibin], official_mu[ibin], passgrl] + out_dict["Zee"] + out_dict["Zmumu"] + [zll_lumi, zll_lumi_err]
431 csvwriter.writerow(out_write)
432
433if (len(zee_missing_lbs) > 0 or len(zmumu_missing_lbs) > 0):
434 print("Missing LBs in Zee channel: ", zee_missing_lbs)
435 print("Missing LBs in Zmumu channel: ", zmumu_missing_lbs)
436
437if bad_database:
438 print("WARNING: There was an error retrieving information from the lumi database, likely need to update the tags.")
439
440csvfile.close()
print
void print(char *figname, TCanvas *c1)
Definition TRTCalib_StrawStatusPlots.cxx:26
pow
constexpr int pow(int base, int exp) noexcept
Definition ap_fixedTest.cxx:15
array
STL class.
Generated on for ATLAS Offline Software by doxygen 1.14.0