Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pfvalidation #11

Merged
merged 2 commits into from
Jan 24, 2019
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 62 additions & 0 deletions Validation/RecoParticleFlow/test/IO.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
# This file contains classes and methods for basic input and ouput operations
# for the TeraJet jet particology package.

import ROOT as r

# Not used yet, let's see if we need one later
class Event:

def __init__(self, evt):
#self.__filename = filename
self._njet = evt.nJet
self._eventno = evt.event
self._run = evt.run
self._lumi = evt.luminosityBlock

# Getters
@property
def njet(self):
return self._njet
@property
def eventno(self):
return self._eventno
@property
def run(self):
return self._run
@property
def lumi(self):
return self._lumi

#def load_file(self, filename):
# self.__filename = filename
# f = r.TFile(filename,"read")

#def load_vars(self, evt):
# self.njet = evt.nJet
# self.



def load_file(dtype, filename):

f = r.TFile(filename, "read")
if dtype == "SMPJ":
treepath = "ak5/ProcessedTree"
elif dtype == "NANO":
treepath = "Events"
else:
print "Fatal error: no known tree path for data type '%s'" % dtype
exit()
tree = r.TTree()
tree = f.Get(treepath)

# Thanks ROOT/C++, do this to return an object with it's original type
tree.SetDirectory(0)

#tree.AddDirectory(kFalse)
#r.gInterpreter.ProcessLine('TTree::AddDirectory(kFalse)')
entries = tree.GetEntriesFast()

return tree


193 changes: 193 additions & 0 deletions Validation/RecoParticleFlow/test/Jet.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
import IO
import sys

class Jet(IO.Event):

def __init__(self, evt, idx, isMC, dtier):

# Jet NanoAOD
if (dtier == "naod"):
self._idx = idx
self._pt = evt.Jet_pt[idx]
self._rawpt = evt.Jet_pt[idx]*(1-evt.Jet_rawFactor[idx])
self._isMC = isMC
self._eta = evt.Jet_eta[idx]
self._phi = evt.Jet_phi[idx]
self._mass = evt.Jet_mass[idx]
self._abseta = abs(self._eta)
self._area = evt.Jet_area[idx]
self._cemf = evt.Jet_chEmEF[idx]
self._nemf = evt.Jet_neEmEF[idx]
self._chf = evt.Jet_chHEF[idx]
self._nhf = evt.Jet_neHEF[idx]

self._jetid = evt.Jet_jetId[idx]
self._rawfac = evt.Jet_rawFactor[idx]
self._nconst = evt.Jet_nConstituents[idx]
self._nelec = evt.Jet_nElectrons[idx]
self._nmuon = evt.Jet_nMuons[idx]
self._bgtags = [evt.Jet_btagCMVA[idx],evt.Jet_btagCSVV2[idx],
evt.Jet_btagDeepB[idx],
evt.Jet_btagDeepC[idx]]#evt.Jet_btagDeepFlavB[idx]]
self._qgl = evt.Jet_qgl[idx]

# Vars only in MC
if self._isMC:
self._pflav = evt.Jet_partonFlavour[idx]
self._hflav = evt.Jet_hadronFlavour[idx]

self._genpt = evt.GenJet_pt[idx]
self._geneta = evt.GenJet_eta[idx]
self._genphi = evt.GenJet_phi[idx]
self._genmass = evt.GenJet_mass[idx]
self._genabseta = abs(self._geneta)




# True muon fraction only in Sal's data for now
# (Except for now it is in official data too, thanks to meeee!!)
#if not self._isMC:
self._muf = evt.Jet_muEF[idx]

# Calculate tight jetID boolean
foo = [4,5,6,7] # Hack for getting ID from binary
# (someone please tell me how to do it less ugly!)


# This is now TightLeptonVetoID
if (evt.Jet_jetId[idx] not in foo):
self._jetidt = False
else:
self._jetidt = True

elif (dtier == "djred"):

self._idx = idx
self._pt = evt.pTWJ_j1 if idx == 0 else evt.pTWJ_j2
self._isMC = isMC
self._eta = evt.Jet_eta[idx]
self._abseta = abs(self._eta)
self._area = evt.Jet_area[idx]
self._cemf = evt.Jet_chEmEF[idx]
self._nemf = evt.Jet_neEmEF[idx]
self._chf = evt.Jet_chHEF[idx]
self._nhf = evt.Jet_neHEF[idx]
self._muf = evt.Jet_muEF[idx]


else:
print "Fatal error: unknown data tier '%s'. Cya kook" % dtier
sys.exit()

# Getters
@property
def idx(self):
return self._idx
@property
def isMC(self):
return self._isMC
@property
def pt(self):
return self._pt
@property
def rawpt(self):
return self._rawpt
@property
def phi(self):
return self._phi
@property
def area(self):
return self._area
@property
def cemf(self):
return self._cemf
@property
def nemf(self):
return self._nemf
@property
def chf(self):
return self._chf
@property
def nhf(self):
return self._nhf
@property
def muf(self):
return self._muf
@property
def jetidt(self):
return self._jetidt
@property
def rawfac(self):
return self._rawfac
@property
def jetidt(self):
return self._jetidt
@property
def eta(self):
return self._eta
@property
def abseta(self):
return self._abseta
@property
def nconst(self):
return self._nconst
@property
def nelec(self):
return self._nelec
@property
def nmuon(self):
return self._nmuon
@property
def btags(self):
return self._btags
@property
def qgl(self):
return self._qgl
@property
def mass(self):
return self._mass

# MC-only vars need special care
@property
def genpt(self):
if self._isMC:
return self._genpt
else:
return None
@property
def geneta(self):
if self._isMC:
return self._geneta
else:
return None
@property
def genphi(self):
if self._isMC:
return self._genphi
else:
return None
@property
def genmass(self):
if self._isMC:
return self._genmass
else:
return None
@property
def genabseta(self):
if self._isMC:
return self._genabseta
else:
return None
@property
def pflav(self):
if self._isMC:
return self._pflav
else:
return None


# Setters
@pt.setter
def pt(self, val):
self._pt = val
133 changes: 133 additions & 0 deletions Validation/RecoParticleFlow/test/pfreso.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@

import ROOT as r
import sys
import array as ar


shf=-7; shf2=1
colors = [r.kRed-4,r.kRed+shf,r.kBlue+shf,r.kGreen+shf+1,r.kCyan+shf+1,
r.kMagenta+shf,r.kGray+1,r.kOrange,r.kCyan]
darkcolors = [r.kRed-2,r.kRed+shf2,r.kBlue+shf2,r.kGreen+shf2+1,r.kCyan+shf2+1,
r.kMagenta+shf2,r.kGray+2,r.kOrange+1,r.kCyan+1]

fullmarkers = [r.kFullCircle,r.kFullSquare,r.kFullDiamond,r.kFullTriangleUp,
r.kFullTriangleDown,r.kFullStar,r.kFullCross,31]
openmarkers = [r.kOpenCircle,r.kOpenSquare,r.kOpenDiamond,
r.kOpenTriangleUp,r.kOpenTriangleDown,r.kOpenStar,r.kOpenCross]

def pf_resolution(forig, fdev, path, apdx):

# Eta bin from apdx
if (apdx == "eta05"):
info = "|#eta| < 0.5"
if (apdx == "eta13"):
info = "|#eta| < 1.3"
if (apdx == "eta21"):
info = "1.3 < |#eta| < 2.1"
if (apdx == "eta25"):
info = "2.1 < |#eta| < 2.5"
if (apdx == "eta30"):
info = "2.5 < |#eta| < 3.0"

creso = r.TCanvas("resolution","Jet pT resolution", 600, 600)

# Loop responses i.e. dist widths to an array. Or TProfile.
ptbins = ar.array('d',[10,24,32,43,56,74,97,133,174,245,300,362,430,
507,592,686,846,1032,1248,1588,2000,2500,3000,4000,6000])
nptbins = len(ptbins)-1

#etabins = [[0,0.5],[0,1.3],[1.3,2.1],[2.1,2.5],[2.5,3.0]]

preso = r.TProfile("preso","Jet pT resolution",nptbins,ptbins)
presodev = r.TProfile("presodev","Jet pT resolution",nptbins,ptbins)

leg = r.TLegend(0.5,0.75,0.9,0.85)
leg.SetFillStyle(0);
leg.SetBorderSize(0);
leg.SetEntrySeparation(0);
leg.AddEntry(preso,"Original PF","pl")
leg.AddEntry(presodev,"Modified PF","pl")



for idx in range(nptbins):
elow = ptbins[idx]
ehigh = ptbins[idx+1]
h = forig.Get("%sreso_dist_%i_%i_%s" % (path,elow,ehigh,apdx))
hdev = fdev.Get("%sreso_dist_%i_%i_%s" % (path,elow,ehigh,apdx))
std = h.GetStdDev()
err = h.GetStdDevError()
stddev = hdev.GetStdDev()
errdev = hdev.GetStdDevError()
preso.Fill(elow,std)
presodev.Fill(elow,stddev)
# Set errors
preso.SetBinError(idx,err)
presodev.SetBinError(idx,errdev)

preso.SetXTitle("Jet pT [GeV]")
preso.SetYTitle("#sigma(p_{T}^{reco}/p_{T}^{gen})")
preso.GetYaxis().SetTitleOffset(1.3)

preso.SetMarkerStyle(fullmarkers[0])
preso.SetMarkerColor(colors[1])
preso.SetLineColor(colors[1])

presodev.SetMarkerStyle(openmarkers[1])
presodev.SetMarkerColor(colors[2])
presodev.SetLineColor(colors[2])


r.gStyle.SetOptStat(0)
preso.SetAxisRange(10,3000,"X")
preso.SetAxisRange(0,0.5,"Y")
preso.Draw("p")
presodev.Draw("psame")
leg.Draw()
creso.SetLogx()
preso.GetXaxis().SetMoreLogLabels()
preso.GetXaxis().SetNoExponent()

t = r.TLatex()
t.SetTextColor(1)
t.SetTextSize(0.05)
t.SetNDC()

t.DrawLatex(0.15,0.15,info)

creso.SaveAs("res/pfreso_%s.pdf" % apdx)



def main():

forig = r.TFile("/home/juska/Dropbox/postdoc/pfdev/files/"
"DQM_V0001_R000000001__Global__CMSSW_X_Y_Z__RECO.root","read")

path = "DQMData/Run 1/Physics/Run summary/JetResponse/"

#dresoorig = f.Get(path)
#print dresoorig

fdev = r.TFile("/home/juska/Dropbox/postdoc/pfdev/files/"
"DQM_V0001_R000000001__Global__CMSSW_X_Y_Z__RECO.root","read")

#print fdev

dresodev = fdev.Get(path)


#print dresodev
# h0 = fdev.Get("DQMData/Run 1/Physics/Run summary/JetResponse/reso_dist_1032_1248_eta05")




pf_resolution(forig, fdev, path, "eta05")
pf_resolution(forig, fdev, path, "eta13")
pf_resolution(forig, fdev, path, "eta21")
pf_resolution(forig, fdev, path, "eta25")
pf_resolution(forig, fdev, path, "eta30")

#input("Press enter to quit")
main()