6 from array
import array
12 from DQUtils
import fetch_iovs
13 from DQUtils.iov_arrangement
import inverse_lblb
14 import ZLumiScripts.tools.zlumi_alleff
as dq_eff
15 import ZLumiScripts.tools.zlumi_mc_cf
as dq_cf
18 from DataQualityUtils
import doZLumi
21 import ZLumiScripts.tools.toys
as toys
23 ROOT.gROOT.SetBatch(ROOT.kTRUE)
24 ROOT.gStyle.SetOptStat(0)
26 logging.basicConfig(format=
'%(levelname)s:%(message)s', level=logging.INFO)
28 parser = argparse.ArgumentParser()
29 parser.add_argument(
'--infile', type=str, help=
'input HIST file')
30 parser.add_argument(
'--grl', type=str, help=
'Specify an input GRL')
31 parser.add_argument(
'--tag', type=str, help=
'Lumi tag', default=
'OflLumiAcct-Run3-003')
32 parser.add_argument(
'--useofficial', action=
'store_true', help=
'Use official lumi folder (otherwise, use OflLumiAcct')
33 parser.add_argument(
'--lumifolder', type=str, help=
'Lumi folder', default=
'/TRIGGER/OFLLUMI/OflPrefLumi')
34 parser.add_argument(
'--lumitag', type=str, help=
'Lumi tag', default=
'OflLumi-Run3-003')
35 parser.add_argument(
'--outdir', type=str, help=
'Directory to dump plots', default=
'plots')
36 parser.add_argument(
'--dblivetime', action=
'store_true', help=
'Look up livetime from DB')
37 parser.add_argument(
'--campaign', type=str, help=
'mc16a/d/e, mc21, mc23a')
39 args = parser.parse_args()
40 campaign = args.campaign
42 ZPURITYFACTOR = 0.9935
43 if args.campaign
in [
"mc21",
"mc23a"]:
49 zmumu_missing_lbs = []
54 fin = ROOT.TFile.Open(args.infile)
56 for key
in fin.GetListOfKeys():
57 if key.GetName().startswith(
'run_'):
58 runname = key.GetName()
59 runnumber =
int(runname.replace(
'run_',
''))
63 logging.critical(
"Can't find run_* directory in input file %s", args.infile)
66 print(
"Starting HIST to CSV conversion for Run ", runnumber)
70 os.system(
"mkdir -p " + out_dir)
71 out_dir +=
"/" + runname +
".csv"
73 out_dir = runname +
".csv"
75 if args.campaign
in [
"mc21",
"mc23a"]:
76 lb_length_name =
'%s/GLOBAL/DQTGlobalWZFinder/duration_vs_LB' % runname
77 livetime_name =
'%s/GLOBAL/DQTGlobalWZFinder/avgLiveFrac_vs_LB' % runname
79 lb_length_name =
'%s/GLOBAL/DQTGlobalWZFinder/m_lblength_lb' % runname
80 livetime_name =
'%s/GLOBAL/DQTGlobalWZFinder/m_livetime_lb' % runname
84 grlReader = ROOT.Root.TGoodRunsListReader(args.grl)
86 grl = grlReader.GetMergedGRLCollection()
88 grlname =
'grl_'+
str(runnumber)+
'.xml'
89 grl_file = doZLumi.makeGRL(runnumber,
'PHYS_StandardGRL_All_Good', grlname)
91 lb_length_old = fin.Get(lb_length_name)
92 lbmin, lbmax = lb_length_old.FindFirstBinAbove(0, 1, 0, -1), lb_length_old.FindLastBinAbove(0, 1, 0, -1)
93 lb_length = ROOT.TProfile(
'lb_length',
'LB length',
int(lbmax-lbmin), lbmin, lbmax)
95 for i
in range(lbmin, lbmax):
96 lb_length.Fill(i, lb_length_old.GetBinContent(i))
98 logging.info(
'low, high LBs: %s, %s', lbmin, lbmax)
101 logging.info(
'Starting livetime lookup ... (remove when we have a proper in-file implementation ...)')
102 livetime = ROOT.TProfile(
'livetime',
'Livetime',
int(lbmax-lbmin), lbmin, lbmax)
104 livetime = fin.Get(livetime_name)
106 official_lum_zero = ROOT.TProfile(
'official_lum_zero',
'official inst luminosity',
int(lbmax-lbmin), lbmin, lbmax)
107 official_mu = ROOT.TProfile(
'official_mu',
'official mu',
int(lbmax-lbmin), lbmin, lbmax)
111 iovs_acct =
fetch_iovs(
'COOLOFL_TRIGGER::/TRIGGER/OFLLUMI/LumiAccounting', lbtime.first.since, lbtime.last.until, tag=args.tag)
114 print(
"Using official lumitag", args.lumitag)
115 iovs_lum =
fetch_iovs(
'COOLOFL_TRIGGER::%s' % args.lumifolder, lblb.first.since, lblb.last.until, tag=args.lumitag, channels=[0])
120 lb_start_end[iov.since & 0xffffffff] = (iov.StartTime/1e9, iov.EndTime/1e9)
122 for iov
in iovs_acct:
123 if not lbmin-1 < iov.LumiBlock < lbmax:
125 lb_lhcfill[iov.LumiBlock] = iov.FillNumber
127 livetime.Fill(iov.LumiBlock, iov.LiveFraction)
129 if not args.useofficial:
130 official_lum_zero.Fill(iov.LumiBlock, iov.InstLumi/1e3)
131 official_mu.Fill(iov.LumiBlock, iov.AvEvtsPerBX)
133 offlumiov = [_
for _
in iovs_lum
if _.since.lumi==iov.LumiBlock]
134 if len(offlumiov) != 1:
135 print(
'MAJOR PROBLEM, LUMI IOV MISMATCH', iov.LumiBlock)
138 offlumiov = offlumiov[0]
139 official_lum_zero.Fill(iov.LumiBlock, offlumiov.LBAvInstLumi/1e3)
140 official_mu.Fill(iov.LumiBlock, offlumiov.LBAvEvtsPerBX)
142 lb_full = lb_length.Clone(
'lb_full').ProjectionX()
143 divisor = lb_length.Clone(
'divisor').ProjectionX()
144 px = livetime.ProjectionX()
151 hto = fin.Get(
'%s/GLOBAL/DQTGlobalWZFinder/m_ele_template_os' % (runname))
152 hts = fin.Get(
'%s/GLOBAL/DQTGlobalWZFinder/m_ele_template_ss' % (runname))
170 lb_minus_one_reco_eff = [1.0, 1.0, 1.0]
171 for ibin
in range(1,
int(lbmax-lbmin)+1):
172 this_lb =
int(lb_full.GetBinCenter(ibin))
176 hmo = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
177 hms = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
178 hno = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
179 hns = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
182 eff, err = dq_eff.template_method(hmo, hms, hno, hns, hto, hts,
False, lb_minus_one_reco_eff)
183 lb_minus_one_reco_eff = [eff, err, this_lb]
184 except AttributeError:
188 if err != 0
and eff != 0:
190 if pileup
not in o_recoeff_fit:
191 o_recoeff_fit[pileup] = weight * eff
192 o_recoerr_fit[pileup] = weight
194 o_recoeff_fit[pileup] += weight * eff
195 o_recoerr_fit[pileup] += weight
201 for pileup
in o_recoeff_fit:
202 arr_mu.append(
float(pileup))
203 arr_rec_eff.append(o_recoeff_fit[pileup]/o_recoerr_fit[pileup])
204 arr_rec_err.append(1/
pow(o_recoerr_fit[pileup], 0.5))
210 arr_rec_eff.append(1)
211 arr_rec_err.append(1)
213 tg_fit = ROOT.TGraphErrors(len(arr_mu), arr_mu, arr_rec_eff, ROOT.nullptr, arr_rec_err)
214 if len(o_recoeff_fit) == 0
or len(o_recoeff_fit) == 1:
216 elif len(o_recoeff_fit) == 2:
218 elif len(o_recoeff_fit) > 2:
221 tg_fit.Fit(fit_type,
"q")
224 csvwriter = csv.writer(csvfile, delimiter=
',')
225 csvwriter.writerow([
'FillNum',
'RunNum',
'LBNum',
'LBStart',
'LBEnd',
'LBLive',
'LBFull',
'OffLumi',
'OffMu',
'PassGRL',
226 'ZeeRaw',
'ZeeRawErr',
'ZeeN1',
'ZeeN2',
'ZeeEffTrig',
'ZeeErrTrig',
'ZeeEffReco',
'ZeeErrReco',
'ZeeEffComb',
'ZeeErrComb',
'ZeeEffAComb',
'ZeeErrAComb',
'ZeeDefTrig',
'ZeeDefReco',
'ZeeLumi',
'ZeeLumiErr',
'ZeeRate',
227 'ZmumuRaw',
'ZmumuRawErr',
'ZmumuN1',
'ZmumuN2',
'ZmumuEffTrig',
'ZmumuErrTrig',
'ZmumuEffReco',
'ZmumuErrReco',
'ZmumuEffComb',
'ZmumuErrComb',
'ZmumuEffAComb',
'ZmumuErrAComb',
'ZmumuDefTrig',
'ZmumuDefReco',
'ZmumuLumi',
'ZmumuLumiErr',
'ZmumuRate',
228 'ZllLumi',
'ZllLumiErr'])
230 lb_minus_one_reco_eff = {}
231 lb_minus_one_reco_eff[
"Zee"] = [1.0, 1.0, 1]
232 lb_minus_one_reco_eff[
"Zmumu"] = [1.0, 1.0, 1]
234 lb_minus_one_trig_eff = {}
235 lb_minus_one_trig_eff[
"Zee"] = [1.0, 1.0, 1]
236 lb_minus_one_trig_eff[
"Zmumu"] = [1.0, 1.0, 1]
240 for ibin
in range(1,
int(lbmax-lbmin)+1):
242 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]
243 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 this_lb =
int(lb_full.GetBinCenter(ibin))
246 loclivetime = divisor[ibin]
248 if runnumber == 302831
and this_lb < 11:
250 elif runnumber == 329835
and this_lb < 554:
252 elif runnumber == 310247
and (this_lb == 442
or this_lb == 462):
254 elif runnumber == 281385
and this_lb < 197:
255 loclivetime *= 4.0/6.0
256 elif runnumber == 281385
and this_lb < 375:
257 loclivetime *= 5.0/6.0
258 elif runnumber == 286367:
259 loclivetime *= 5.0/6.0
262 this_fill = lb_lhcfill[this_lb]
268 for channel
in [
"Zee",
"Zmumu"]:
269 if grl
and not grl.HasRunLumiBlock(runnumber, this_lb):
276 lb =
"lb_" +
str(this_lb)
279 hname = runname +
'/lb_'+
str(
int(ibin+lbmin-0.5))+
'/GLOBAL/DQTGlobalWZFinder/m_Z_mass_opsele'
282 h = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_eltrigtp_matches_os' % (runname, lb))
283 hmo = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_os' % (runname, lb))
284 hms = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_good_ss' % (runname, lb))
285 hno = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_os' % (runname, lb))
286 hns = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_ele_tight_bad_ss' % (runname, lb))
291 if args.campaign
in [
"mc21",
"mc23a"]:
295 elif channel ==
'Zmumu':
297 hname = runname +
'/lb_'+
str(
int(ibin+lbmin-0.5))+
'/GLOBAL/DQTGlobalWZFinder/m_Z_mass_opsmu'
300 h = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_mutrigtp_matches' % (runname, lb))
301 hmo = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_os' % (runname, lb))
302 hms = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_ss' % (runname, lb))
303 hno = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_os' % (runname, lb))
304 hns = fin.Get(
'%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_ss' % (runname, lb))
305 if args.campaign
in [
"mc21",
"mc23a"]:
308 ACCEPTANCE = 0.3323224
311 z_m = fin.Get(hname).Integral()
312 z_merr = math.sqrt(z_m)
313 except AttributeError:
315 zee_missing_lbs.append(this_lb)
316 elif channel ==
"Zmumu":
317 zmumu_missing_lbs.append(this_lb)
324 N1 = h.GetBinContent(2)
325 N2 = h.GetBinContent(3)
327 eff_trig, err_trig, arr_trig, arr_NZ = toys.toy_trigeff(N1, N2, ntoys)
329 eff_trig, err_trig = dq_eff.trig_tag_and_probe(h, lb_minus_one_trig_eff[channel])
331 eff_trig, err_trig = 0.0, 0.0
332 except AttributeError:
333 eff_trig, err_trig = 0.0, 0.0
334 if channel ==
"Zmumu":
339 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 eff_reco, err_reco = dq_eff.reco_tag_and_probe(hmo, hms, hno, hns, lb_minus_one_reco_eff[channel])
343 eff_reco, err_reco = 0.0, 0.0
344 except AttributeError:
345 eff_reco, err_reco = 0.0, 0.0
346 elif channel ==
"Zee":
349 bin1 = hmo.GetXaxis().
FindBin(75000)
350 bin2 = hmo.GetXaxis().
FindBin(104000)
353 matchos_peak = hmo.Integral(bin1, bin2)
354 matchos_tail = hmo.Integral(bin3, bin4)
355 matchss_tail = hms.Integral(bin3, bin4)
356 nomatchos_peak = hno.Integral(bin1, bin2)
357 nomatchos_tail = hno.Integral(bin3, bin4)
358 templateos_peak = hto.Integral(bin1, bin2)
359 templateos_tail = hto.Integral(bin3, bin4)
360 templatess_tail = hts.Integral(bin3, bin4)
361 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 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]))
366 eff_reco, err_reco = 0.0, 0.0
367 except AttributeError:
368 eff_reco, err_reco = 0.0, 0.0
373 defaulted_reco_eff = 0
374 defaulted_trig_eff = 0
375 if eff_reco == lb_minus_one_reco_eff[channel][0]:
376 defaulted_reco_eff = 1
377 if eff_trig == lb_minus_one_trig_eff[channel][0]:
378 defaulted_trig_eff = 1
381 lb_minus_one_reco_eff[channel] = [eff_reco, err_reco, this_lb]
383 lb_minus_one_trig_eff[channel] = [eff_trig, err_trig, this_lb]
386 arr_comb = (1.0 - (1.0 - arr_trig)**2) * (arr_reco)**2
387 nonan_arr_comb = arr_comb[~np.isnan(arr_comb)]
388 eff_comb = np.median(nonan_arr_comb)
389 err_comb = nonan_arr_comb.std()
391 eff_comb = (1-(1-eff_trig)**2)*(eff_reco)**2
392 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 eff_Acomb = ACCEPTANCE * eff_comb
395 err_Acomb = ACCEPTANCE * err_comb
398 effcy = arr_comb * dq_cf.correction(official_mu[ibin], channel, campaign, runnumber)
400 effcy = eff_comb * dq_cf.correction(official_mu[ibin], channel, campaign, runnumber)
401 effcyerr = err_comb * dq_cf.correction(official_mu[ibin], channel, campaign, runnumber)
403 zlumi = zlumistat = zrate = 0.0
404 CORRECTIONS = ZPURITYFACTOR/ACCEPTANCE/ZXSEC
406 if do_toys
and loclivetime != 0.0:
407 arr_zlumi = np.divide(arr_NZ, effcy) * (CORRECTIONS)/loclivetime
408 arr_zlumi = arr_zlumi[~np.isnan(arr_zlumi)]
409 zlumi = np.median(arr_zlumi)
410 zlumistat = arr_zlumi.std()
411 zrate = zlumi / CORRECTIONS
412 elif (loclivetime != 0.0
and effcy != 0.0):
413 zlumi = (z_m/effcy)*(CORRECTIONS)/loclivetime
414 zlumistat = math.sqrt(
pow(z_merr/effcy, 2) +
pow(z_m/effcy**2*effcyerr, 2))*CORRECTIONS/loclivetime
415 zrate = zlumi / CORRECTIONS
418 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 lumi_index = len(out_dict[
'Zee'])-3
421 error_index = len(out_dict[
'Zee'])-2
422 zll_lumi = (out_dict[
'Zee'][lumi_index] + out_dict[
'Zmumu'][lumi_index])/2
423 zll_lumi_err = 0.5 * math.sqrt(
pow(out_dict[
'Zee'][error_index], 2) +
pow(out_dict[
'Zmumu'][error_index], 2) )
425 lb_start, lb_end = lb_start_end[this_lb][0], lb_start_end[this_lb][1]
427 lb_start, lb_end = 0, 0
428 print(
"WARNING: database not filled completely for LB", this_lb)
429 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]
430 csvwriter.writerow(out_write)
432 if (len(zee_missing_lbs) > 0
or len(zmumu_missing_lbs) > 0):
433 print(
"Missing LBs in Zee channel: ", zee_missing_lbs)
434 print(
"Missing LBs in Zmumu channel: ", zmumu_missing_lbs)
437 print(
"WARNING: There was an error retrieving information from the lumi database, likely need to update the tags.")