Skip to content
Snippets Groups Projects
Commit f50e82cb authored by Peter Onyisi's avatar Peter Onyisi
Browse files

Refactor packages for Z lumi scripts 

Former-commit-id: 7546c8f203172c5a86b8027e3861504bac98a5f3
parent 64e1c3f3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
DataQuality/DataQualityTools/scripts/dqt_zlumi_alleff_HIST.py deleted 100755 → 0
+ 0
− 193
View file @ 64e1c3f3
#!/usr/bin/env python
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
import sys, glob
import ROOT
ROOT.gStyle.SetOptStat(0)
ACCEPTANCE = 3.173927e-01
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('infile', type=str, help='input HIST file')
parser.add_argument('--out', type=str, help='output ROOT file')
args = parser.parse_args()
infilename = args.infile
infile = ROOT.TFile.Open(infilename, 'READ')
rundir = None
for k in infile.GetListOfKeys():
if k.GetName().startswith('run_'):
rundir = k.GetName()
break
if not rundir:
print 'Cannot find run directory in input file'
sys.exit(1)
else:
print 'Found rundir', rundir
lbdirs = []
for k in infile.Get(rundir).GetListOfKeys():
if k.GetName().startswith('lb_'):
lbdirs.append(k.GetName())
print 'Now to dump'
lbnums = sorted([int(_[3:]) for _ in lbdirs])
effcyt = ROOT.TH1F('effcyt', 'Trigger efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
lbnums[-1]+0.5)
effcyr = ROOT.TH1F('effcyr', 'Loose muon reco efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
lbnums[-1]+0.5)
effcya = ROOT.TH1F('effcya', 'Combined acc x efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
lbnums[-1]+0.5)
effcydir = ROOT.TH1F('effcydir', 'Direct acc x efficiency', lbnums[-1]-lbnums[0]+1, lbnums[0]-0.5,
lbnums[-1]+0.5)
from array import array
fout = ROOT.TFile(args.out if args.out else '%s_all.root' % rundir[4:], 'RECREATE')
o_run = array('I', [0])
o_lb = array('I', [0])
o_lbwhen = array('d', [0., 0.])
o_z_one = array('f', [0.])
o_z_two = array('f', [0.])
o_trigeff = array('f', [0.])
o_trigeffstat = array('f', [0.])
o_recoeff = array('f', [0.])
o_recoeffstat = array('f', [0.])
o_alleff = array('f', [0.])
o_alleffstat = array('f', [0.])
o_ae = array('f', [0.])
o_aestat = array('f', [0.])
tl = ROOT.TTree( 'lumitree', 'Luminosity tree' )
tl.Branch('run', o_run, 'run/i')
tl.Branch('lb', o_lb, 'lb/i')
tl.Branch('lbwhen', o_lbwhen, 'lbwhen[2]/D')
tl.Branch('z_one', o_z_one, 'z_one/F')
tl.Branch('z_two', o_z_two, 'z_two/F')
tl.Branch('trigeff', o_trigeff, 'trigeff/F')
tl.Branch('trigeffstat', o_trigeffstat, 'trigeffstat/F')
tl.Branch('recoeff', o_recoeff, 'recoeff/F')
tl.Branch('recoeffstat', o_recoeffstat, 'recoeffstat/F')
tl.Branch('alleff', o_alleff, 'alleff/F')
tl.Branch('alleffstat', o_alleffstat, 'alleffstat/F')
tl.Branch('ae', o_ae, 'ae/F')
tl.Branch('aestat', o_aestat, 'aestat/F')
from DQUtils import fetch_iovs
#rset=set(_[0] for _ in rlb)
#print rset
lblb = fetch_iovs("LBLB", runs=int(rundir[4:])).by_run
for lb in sorted(lbdirs):
h = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_mutrigtp_matches' % (rundir, lb))
hmo = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_os' % (rundir, lb))
hms = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_match_ss' % (rundir, lb))
hno = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_os' % (rundir, lb))
hns = infile.Get('%s/%s/GLOBAL/DQTGlobalWZFinder/m_muloosetp_nomatch_ss' % (rundir, lb))
lbnum = int(lb[3:])
yld = (h[2], h[3])
ylderr = (h.GetBinError(2), h.GetBinError(3))
#print yld, ylderr
A, B = yld
o_z_one[0], o_z_two[0] = yld
if B == 0: continue
eff = 1./(float(A)/B/2.+1.)
inverrsq = ((1/2./B)*ylderr[0])**2+((A/2./B**2)*ylderr[1])**2
o_trigeff[0] = eff
o_trigeffstat[0] = (inverrsq**.5)*(eff**2)
o_run[0], o_lb[0] = int(rundir[4:]), lbnum
try:
iov = lblb[int(rundir[4:])][lbnum-1]
o_lbwhen[0], o_lbwhen[1] = iov.StartTime/1e9, iov.EndTime/1e9
except Exception, e:
o_lbwhen[0], o_lbwhen[1] = 0, 0
effcyt.SetBinContent(lbnum-lbnums[0]+1, eff)
effcyt.SetBinError(lbnum-lbnums[0]+1, o_trigeffstat[0])
def extract(histogram):
dbl = ROOT.Double()
rv1 = histogram.IntegralAndError(21, 30, dbl)
return (rv1, float(dbl))
matchos, matchoserr = extract(hmo)
matchss, matchsserr = extract(hms)
nomatchos, nomatchoserr = extract(hno)
nomatchss, nomatchsserr = extract(hns)
print lb
print ' ->', matchos, matchoserr
print ' ->', matchss, matchsserr
print ' ->', nomatchos, nomatchoserr
print ' ->', nomatchss, nomatchsserr
A = float(matchos-matchss)
Aerr = (matchoserr**2+matchsserr**2)**.5
B = float(nomatchos-nomatchss)
Berr = (nomatchoserr**2+nomatchsserr**2)**.5
if B == 0: Berr = 1.
if A == 0:
eff = 1.
inverrsq = 1.
else:
eff = 1./(1.+B/A)
inverrsq = ((-B/A**2)*Aerr)**2+((1./A)*Berr)**2
o_recoeff[0] = eff
o_recoeffstat[0] = (inverrsq**.5)*(eff**2)
effcyr.SetBinContent(lbnum-lbnums[0]+1, eff)
effcyr.SetBinError(lbnum-lbnums[0]+1, o_recoeffstat[0])
o_ae[0] = ACCEPTANCE*(1-(1-o_trigeff[0])**2)*(o_recoeff[0])**2
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
o_alleff[0] = (1-(1-o_trigeff[0])**2)*(o_recoeff[0])**2
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
effcya.SetBinContent(lbnum-lbnums[0]+1, o_ae[0])
effcya.SetBinError(lbnum-lbnums[0]+1, o_aestat[0])
tl.Fill()
tl.Write()
print 'Done'
c1 = ROOT.TCanvas()
effcya.SetMarkerStyle(21)
effcya.SetMarkerColor(ROOT.kBlue)
effcya.GetYaxis().SetRangeUser(0.27,0.29)
effcya.Draw('PE')
c1.Print('%s_combined_efficiency.eps' % rundir)
fout.WriteTObject(effcya)
fout.Close()
sumweights = infile.Get('%s/GLOBAL/DQTDataFlow/m_sumweights' % rundir)
ctr = infile.Get('%s/GLOBAL/DQTGlobalWZFinder/m_Z_Counter_mu' % rundir)
if sumweights:
for ibin in xrange(1,sumweights.GetNbinsX()+1):
o_lb[0] = int(sumweights.GetBinCenter(ibin))
ctrbin = ctr.FindBin(o_lb[0])
print ibin, o_lb[0], sumweights[ibin], ctr[ctrbin]
if sumweights[ibin] == 0: continue
p = ctr[ctrbin]/sumweights[ibin]
o_alleff[0]=p
try:
o_alleffstat[0]=(p*(1-p))**.5*(sumweights.GetBinError(ibin)/sumweights[ibin])
except ValueError:
o_alleffstat[0]=(sumweights.GetBinError(ibin)/sumweights[ibin])
effcydir.SetBinContent(effcydir.FindBin(o_lb[0]), p)
effcydir.SetBinError(effcydir.FindBin(o_lb[0]), o_alleffstat[0])
effcya.GetYaxis().SetRangeUser(0.27,0.31)
effcya.Draw('PE')
effcydir.SetMarkerStyle(20)
effcydir.SetMarkerColor(ROOT.kRed)
effcydir.Draw('SAME,PE')
leg=ROOT.TLegend(0.65, 0.7, 0.89, 0.89)
leg.AddEntry(effcya, 'Predicted A#epsilon', 'PE')
leg.AddEntry(effcydir, 'Actual A#epsilon', 'PE')
leg.Draw()
c1.Print('%s_tp_comparison.eps' % rundir)
effcyrat=effcydir.Clone()
effcyrat.Divide(effcya)
effcyrat.SetTitle('MC Correction Factor')
effcyrat.SetXTitle('<#mu>')
effcyrat.Draw('PE')
effcyrat.Fit('pol1')
c1.Print('%s_tp_correction.eps' % rundir)
DataQuality/DataQualityTools/scripts/dqt_zlumi_combine_lumi.py deleted 100755 → 0
+ 0
− 126
View file @ 64e1c3f3
#!/usr/bin/env python
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
import ROOT
import sys
import argparse
parser = argparse.ArgumentParser()
parser.add_argument('recofile', type=str, help='File with per-LB yields')
parser.add_argument('efffile', type=str, help='File with efficiencies')
parser.add_argument('outfile', type=str, help='Output file')
parser.add_argument('--nlb', type=int, help='# of LBs to combine',
default=20)
args = parser.parse_args()
recozfname = args.recofile
effzfname = args.efffile
outfname = args.outfile
LUMIBLOCKS = args.nlb
ACCEPTANCE = 3.173927e-01
ZXSEC=1.929
ZPURITYFACTOR=0.9935
def correction(mu):
#return 1.04589-0.000158187*mu
return 1.04524-0.000108956*mu
#return 1.04589-0.000158187*20
recozfile = ROOT.TFile.Open(recozfname)
effzfile = ROOT.TFile.Open(effzfname)
recoztree = recozfile.lumitree
effztree = effzfile.lumitree
entrydict = {}
for i in xrange(recoztree.GetEntries()):
recoztree.GetEntry(i)
if not recoztree.pass_grl: continue
# If livetime less than 10 sec, ignore
if recoztree.lblive < 10 : continue
effztree.Draw('alleff:alleffstat', 'run==%s&&lb==%s' % (recoztree.run, recoztree.lb), 'goff')
if effztree.GetSelectedRows() == 0:
print 'Broken for run, lb %s %s' % (recoztree.run, recoztree.lb)
print 'We THINK there are %d events here ...' % (recoztree.zraw)
continue
lbzero = (recoztree.lb // LUMIBLOCKS)*LUMIBLOCKS
run = recoztree.run
if (run, lbzero) not in entrydict: entrydict[(run, lbzero)] = {'zcount': 0., 'zcounterrsq': 0., 'livetime': 0., 'lbwhen': [-1, -1], 'mu': 0., 'offlumi': 0., 'rolleff': 0., 'rollefferrsq': 0., 'lhcfill': recoztree.lhcfill}
thisdict = entrydict[(run, lbzero)]
effcy = (effztree.GetV1()[0]*correction(recoztree.mu))
#thisdict['zcount'] += recoztree.zraw/effcy
#thisdict['zcounterrsq'] += (1/effcy*recoztree.zrawstat)**2+(recoztree.zraw/effcy**2*effztree.GetV2()[0])**2
thisdict['zcount'] += recoztree.zraw
thisdict['zcounterrsq'] += recoztree.zrawstat**2
effcywght = (effztree.GetV2()[0]*correction(recoztree.mu))**2
thisdict['rolleff'] += effcy/effcywght
thisdict['rollefferrsq'] += 1/effcywght
loclivetime = recoztree.lblive
#loclivetime = (recoztree.lbwhen[1]-recoztree.lbwhen[0])
thisdict['livetime'] += loclivetime
thisdict['mu'] += recoztree.mu*loclivetime
thisdict['offlumi'] += recoztree.offlumi*loclivetime
if thisdict['lbwhen'][0] > recoztree.lbwhen[0] or thisdict['lbwhen'][0] == -1:
thisdict['lbwhen'][0] = recoztree.lbwhen[0]
if thisdict['lbwhen'][1] < recoztree.lbwhen[1] or thisdict['lbwhen'][1] == -1:
thisdict['lbwhen'][1] = recoztree.lbwhen[1]
from array import array
fout = ROOT.TFile.Open(outfname, 'RECREATE')
o_run = array('I', [0])
o_lb = array('I', [0])
o_lbwhen = array('d', [0., 0.])
o_zrate = array('f', [0.])
o_zratestat = array('f', [0.])
o_zlumi = array('f', [0.])
o_zlumistat = array('f', [0.])
o_mu = array('f', [0.])
o_alleffcorr = array('f', [0.])
o_alleffcorrstat = array('f', [0.])
o_offlumi = array('f', [0.])
o_lblive = array('f', [0.])
o_lhcfill = array('I', [0])
t = ROOT.TTree( 'lumitree', 'Luminosity tree' )
t.Branch('run', o_run, 'run/i')
t.Branch('lb', o_lb, 'lb/i')
t.Branch('lbwhen', o_lbwhen, 'lbwhen[2]/D')
t.Branch('zrate', o_zrate, 'zrate/F')
t.Branch('zratestat', o_zratestat, 'zratestat/F')
t.Branch('zlumi', o_zlumi, 'zlumi/F')
t.Branch('zlumistat', o_zlumistat, 'zlumistat/F')
t.Branch('offlumi', o_offlumi, 'offlumi/F')
t.Branch('mu', o_mu, 'mu/F')
t.Branch('alleffcorr', o_alleffcorr, 'alleffcorr/F')
t.Branch('alleffcorrstat', o_alleffcorrstat, 'alleffcorrstat/F')
t.Branch('lblive', o_lblive, 'lblive/F')
t.Branch('lhcfill', o_lhcfill, 'lhcfill/i')
for entry, entryval in sorted(entrydict.items()):
if entryval['livetime'] > 0:
entryval['mu'] /= entryval['livetime']
entryval['offlumi'] /= entryval['livetime']
eff = entryval['rolleff']/entryval['rollefferrsq']
efferr = 1/entryval['rollefferrsq']**.5
#print 'LIVETIME2', entryval['livetime']
entryval['zrate'] = entryval['zcount']/eff/entryval['livetime']
entryval['zratestat'] = (entryval['zcounterrsq']/eff + (entryval['zcount']/eff**2*efferr)**2)**.5/entryval['livetime']
o_run[0], o_lb[0] = entry
o_lbwhen[0], o_lbwhen[1] = entryval['lbwhen']
o_zrate[0] = entryval['zrate']
o_zratestat[0] = entryval['zratestat']
o_zlumi[0] = o_zrate[0]*ZPURITYFACTOR/ACCEPTANCE/ZXSEC
o_zlumistat[0] = o_zratestat[0]*ZPURITYFACTOR/ACCEPTANCE/ZXSEC
o_mu[0] = entryval['mu']
o_alleffcorr[0] = eff
o_alleffcorrstat[0] = efferr
o_offlumi[0] = entryval['offlumi']
o_lblive[0] = entryval['livetime']
o_lhcfill[0] = entryval['lhcfill']
if o_zlumi[0] < 4 or o_zlumi[0] > 15:
print o_lb[0], o_zlumi[0], entryval['zcount'], eff, entryval['livetime']
t.Fill()
t.Write()
fout.Close()
DataQuality/DataQualityTools/scripts/dqt_zlumi_display_z_rate.py deleted 100755 → 0
+ 0
− 67
View file @ 64e1c3f3
#!/usr/bin/env python
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
import ROOT
import sys
runnum = sys.argv[1].split('_')[0]
f = ROOT.TFile.Open(sys.argv[1], 'UPDATE')
c1 = ROOT.TCanvas()
lumitree = f.lumitree
lumitree.Draw("zrate:lb+10:zratestat", "", "goff")
print 'Selected rows', lumitree.GetSelectedRows()
if lumitree.GetSelectedRows() > 0:
gr = ROOT.TGraphErrors(lumitree.GetSelectedRows(), lumitree.GetV2(), lumitree.GetV1(), ROOT.nullptr, lumitree.GetV3())
gr.Draw("ap")
gr.GetHistogram().SetXTitle('LB')
gr.GetHistogram().SetYTitle('Fiducial Z yield/second')
gr.SetMarkerStyle(20)
gr.SetTitle('')
f.WriteTObject(gr, 'fid_z_rate')
c1.Update()
c1.Print('%s_fidyield.eps' % runnum)
# dump CSV
import time
csvout = open('%s_zrate.csv' % runnum, 'w')
lumitree.Draw("zrate:lbwhen[0]:lbwhen[1]:lhcfill:lblive:offlumi", "", "goff,para")
timeformat = '%y/%m/%d %H:%M:%S'
#timeformat = '%m/%d %H:%M:%S'
for i in xrange(lumitree.GetSelectedRows()):
zrate = lumitree.GetV1()[i]
instlumi = lumitree.GetVal(5)[i]
livetime = lumitree.GetVal(4)[i]
print >>csvout, '%d, %s, %s, %6f, %6f, %4f, %6f' % (lumitree.GetV4()[i],
time.strftime(timeformat, time.gmtime(lumitree.GetV2()[i])),
time.strftime(timeformat, time.gmtime(lumitree.GetV3()[i])),
lumitree.GetV1()[i],
instlumi/1e3,
instlumi*livetime/1e3,
lumitree.GetV1()[i]*livetime
)
csvout.close()
lumitree.Draw("zlumi:lb+10:zlumistat", "", "goff")
if lumitree.GetSelectedRows() > 0:
gr = ROOT.TGraphErrors(lumitree.GetSelectedRows(), lumitree.GetV2(), lumitree.GetV1(), ROOT.nullptr, lumitree.GetV3())
gr.Draw("ap")
gr.GetHistogram().SetXTitle('LB')
gr.GetHistogram().SetYTitle('Luminosity (x10^{33})')
gr.SetMarkerStyle(20)
gr.SetTitle('')
f.WriteTObject(gr, 'z_lumi')
lumitree.Draw("offlumi:lb+10", "", "goff")
gr2 = ROOT.TGraphErrors(lumitree.GetSelectedRows(), lumitree.GetV2(), lumitree.GetV1())
gr2.SetMarkerStyle(21)
gr2.SetMarkerSize(0.5)
gr2.SetMarkerColor(ROOT.kRed)
gr2.Draw('same,pl')
f.WriteTObject(gr, 'official_lumi')
leg = ROOT.TLegend(0.65, 0.7, 0.89, 0.89)
leg.SetBorderSize(0)
leg.AddEntry(gr, 'Z luminosity', 'pl')
leg.AddEntry(gr2, '-005 calibration', 'pl')
leg.Draw()
c1.Update()
c1.Print('%s_lumicomp.eps' % runnum)
f.WriteTObject(c1, 'lumicomp_canvas')
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment