Skip to content

Commit

Permalink
Merge pull request cms-sw#186 from gpetruc/CMGTools-from-CMSSW_7_2_3
Browse files Browse the repository at this point in the history
Upstream updates from cms-sw#185
  • Loading branch information
gpetruc committed Dec 19, 2014
2 parents 9e0742b + e55fbe1 commit 2440eb3
Show file tree
Hide file tree
Showing 19 changed files with 630 additions and 599 deletions.
Expand Up @@ -20,8 +20,6 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
if not getattr(self.cfg_ana, 'saveTLorentzVectors', False):
fourVectorType.removeVariable("p4")

## Declare how we store floats by default
self.tree.setDefaultFloatType("F"); # otherwise it's "D"

self.collections = {}
self.globalObjects = {}
Expand All @@ -33,6 +31,11 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
if hasattr(cfg_ana,"globalVariables"):
self.globalVariables=cfg_ana.globalVariables

def beginLoop(self, setup) :
super(AutoFillTreeProducer, self).beginLoop(setup)
## Declare how we store floats by default
self.tree.setDefaultFloatType("F"); # otherwise it's "D"

def declareHandles(self):
super(AutoFillTreeProducer, self).declareHandles()
# self.handles['TriggerResults'] = AutoHandle( ('TriggerResults','','HLT'), 'edm::TriggerResults' )
Expand Down
20 changes: 13 additions & 7 deletions PhysicsTools/Heppy/python/analyzers/core/TreeAnalyzerNumpy.py
Expand Up @@ -11,17 +11,22 @@ class TreeAnalyzerNumpy( Analyzer ):

def __init__(self, cfg_ana, cfg_comp, looperName):
super(TreeAnalyzerNumpy,self).__init__(cfg_ana, cfg_comp, looperName)
fileName = '/'.join([self.dirName,
'tree.root'])

isCompressed = self.cfg_ana.isCompressed if hasattr(cfg_ana,'isCompressed') else 1
print 'Compression', isCompressed

self.file = TFile( fileName, 'recreate', '', isCompressed )
self.tree = Tree('tree', self.name)

def beginLoop(self, setup) :
super(TreeAnalyzerNumpy, self).beginLoop(setup)
print setup.services
if "outputfile" in setup.services:
print "Using outputfile given in setup.outputfile"
self.file = setup.services["outputfile"].file
else :
fileName = '/'.join([self.dirName,
'tree.root'])
isCompressed = self.cfg_ana.isCompressed if hasattr(self.cfg_ana,'isCompressed') else 1
print 'Compression', isCompressed
self.file = TFile( fileName, 'recreate', '', isCompressed )
self.tree = Tree('tree', self.name)
self.declareVariables(setup)

def declareVariables(self,setup):
Expand All @@ -30,5 +35,6 @@ def declareVariables(self,setup):

def write(self, setup):
super(TreeAnalyzerNumpy, self).write(setup)
self.file.Write()
if "outputfile" not in setup.services:
self.file.Write()

55 changes: 54 additions & 1 deletion PhysicsTools/Heppy/python/analyzers/core/TriggerBitAnalyzer.py
Expand Up @@ -3,12 +3,15 @@
from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
from PhysicsTools.Heppy.analyzers.core.AutoFillTreeProducer import NTupleVariable
import PhysicsTools.HeppyCore.framework.config as cfg

class TriggerBitAnalyzer( Analyzer ):
def __init__(self, cfg_ana, cfg_comp, looperName ):
super(TriggerBitAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
self.processName = getattr(self.cfg_ana,"processName","HLT")
self.outprefix = getattr(self.cfg_ana,"outprefix", self.processName)
self.unrollbits = ( hasattr(self.cfg_ana,"unrollbits") and self.cfg_ana.unrollbits )


def declareHandles(self):
super(TriggerBitAnalyzer, self).declareHandles()
Expand All @@ -17,22 +20,72 @@ def declareHandles(self):
def beginLoop(self, setup):
super(TriggerBitAnalyzer,self).beginLoop(setup)
self.triggerBitCheckers = []
if self.unrollbits :
self.allPaths = set()
self.triggerBitCheckersSingleBits = []

for T, TL in self.cfg_ana.triggerBits.iteritems():
trigVec = ROOT.vector(ROOT.string)()
for TP in TL:
trigVec.push_back(TP)
outname="%s_%s"%(self.outprefix,T)
if self.unrollbits :
if TP not in self.allPaths :
self.allPaths.update(TP)
trigVecBit = ROOT.vector(ROOT.string)()
trigVecBit.push_back(TP)
outname="%s_BIT_%s"%(self.outprefix,TP)
if not hasattr(setup ,"globalVariables") :
setup.globalVariables = []
setup.globalVariables.append( NTupleVariable(outname, eval("lambda ev: ev.%s" % outname), help="Trigger bit %s"%TP) )
self.triggerBitCheckersSingleBits.append( (TP, ROOT.heppy.TriggerBitChecker(trigVecBit)) )

outname="%s_%s"%(self.outprefix,T)
if not hasattr(setup ,"globalVariables") :
setup.globalVariables = []
setup.globalVariables.append( NTupleVariable(outname, eval("lambda ev: ev.%s" % outname), help="OR of %s"%TL) )
self.triggerBitCheckers.append( (T, ROOT.heppy.TriggerBitChecker(trigVec)) )


def process(self, event):
self.readCollections( event.input )
triggerResults = self.handles['TriggerResults'].product()
for T,TC in self.triggerBitCheckers:
outname="%s_%s"%(self.outprefix,T)
setattr(event,outname, TC.check(event.input.object(), triggerResults))
if self.unrollbits :
for TP,TC in self.triggerBitCheckersSingleBits:
outname="%s_BIT_%s"%(self.outprefix,TP)
setattr(event,outname, TC.check(event.input.object(), triggerResults))


return True


setattr(TriggerBitAnalyzer,"defaultConfig",cfg.Analyzer(
TriggerBitAnalyzer, name="TriggerFlags",
processName = 'HLT',
triggerBits = {
# "<name>" : [ 'HLT_<Something>_v*', 'HLT_<SomethingElse>_v*' ]
}
)
)
setattr(TriggerBitAnalyzer,"defaultEventFlagsConfig",cfg.Analyzer(
TriggerBitAnalyzer, name="EventFlags",
processName = 'PAT',
outprefix = 'Flag',
triggerBits = {
"HBHENoiseFilter" : [ "Flag_HBHENoiseFilter" ],
"CSCTightHaloFilter" : [ "Flag_CSCTightHaloFilter" ],
"hcalLaserEventFilter" : [ "Flag_hcalLaserEventFilter" ],
"EcalDeadCellTriggerPrimitiveFilter" : [ "Flag_EcalDeadCellTriggerPrimitiveFilter" ],
"goodVertices" : [ "Flag_goodVertices" ],
"trackingFailureFilter" : [ "Flag_trackingFailureFilter" ],
"eeBadScFilter" : [ "Flag_eeBadScFilter" ],
"ecalLaserCorrFilter" : [ "Flag_ecalLaserCorrFilter" ],
"trkPOGFilters" : [ "Flag_trkPOGFilters" ],
"trkPOG_manystripclus53X" : [ "Flag_trkPOG_manystripclus53X" ],
"trkPOG_toomanystripclus53X" : [ "Flag_trkPOG_toomanystripclus53X" ],
"trkPOG_logErrorTooManyClusters" : [ "Flag_trkPOG_logErrorTooManyClusters" ],
}
)
)
2 changes: 1 addition & 1 deletion PhysicsTools/Heppy/python/analyzers/core/autovars.py
Expand Up @@ -33,7 +33,7 @@ def __init__(self,name,baseObjectTypes=[],mcOnly=[],variables=[]):
self.variables = variables
def ownVars(self,isMC):
"""Return only my vars, not including the ones from the bases"""
return [ v for v in self.variables if (isMC or not var.mcOnly) ]
return [ v for v in self.variables if (isMC or not v.mcOnly) ]
def allVars(self,isMC):
"""Return all vars, including the base ones. Duplicate bases are not added twice"""
ret = []; names = {}
Expand Down
1 change: 1 addition & 0 deletions PhysicsTools/Heppy/python/analyzers/objects/JetAnalyzer.py
Expand Up @@ -321,6 +321,7 @@ def smearJets(self, event, jets):
minLepPt = 10,
relaxJetId = False,
doPuId = False, # Not commissioned in 7.0.X
doQG = False,
recalibrateJets = False,
shiftJEC = 0, # set to +1 or -1 to get +/-1 sigma shifts
smearJets = True,
Expand Down
200 changes: 195 additions & 5 deletions PhysicsTools/Heppy/python/physicsobjects/Electron.py
@@ -1,6 +1,196 @@
from classes.Electron import Electron
#comment for old
from PhysicsTools.Heppy.physicsobjects.Lepton import Lepton
from PhysicsTools.Heppy.physicsutils.ElectronMVAID import ElectronMVAID_Trig, ElectronMVAID_NonTrig, ElectronMVAID_TrigNoIP
import ROOT
import FastObjects
FastObjects.decorate(ROOT.pat.Electron,Electron)
Electron=FastObjects.AddPhysObjAndCallInit

class Electron( Lepton ):

def __init__(self, *args, **kwargs):
'''Initializing tightIdResult to None. The user is responsible
for setting this attribute externally if he wants to use the tightId
function.'''
super(Electron, self).__init__(*args, **kwargs)
self._physObjInit()

def _physObjInit(self):
self.tightIdResult = None
self.associatedVertex = None
self.rho = None
self._mvaNonTrigV0 = {True:None, False:None}
self._mvaTrigV0 = {True:None, False:None}
self._mvaTrigNoIPV0 = {True:None, False:None}

def electronID( self, id, vertex=None, rho=None ):
if id is None or id == "": return True
if vertex == None and hasattr(self,'associatedVertex') and self.associatedVertex != None: vertex = self.associatedVertex
if rho == None and hasattr(self,'rho') and self.rho != None: rho = self.rho
if id == "POG_MVA_ID_NonTrig": return self.mvaIDLoose()
elif id == "POG_MVA_ID_Trig": return self.mvaIDTight()
elif id == "POG_MVA_ID_NonTrig_full5x5": return self.mvaIDLoose(full5x5=True)
elif id == "POG_MVA_ID_Trig_full5x5": return self.mvaIDTight(full5x5=True)
elif id.startswith("POG_Cuts_ID_"):
return self.cutBasedId(id.replace("POG_Cuts_ID_","POG_"))
for ID in self.electronIDs():
if ID.first == id:
return ID.second
raise RuntimeError, "Electron id '%s' not yet implemented in Electron.py" % id

def cutBasedId(self, wp, showerShapes="auto"):
if "_full5x5" in wp:
showerShapes = "full5x5"
wp = wp.replace("_full5x5","")
elif showerShapes == "auto":
if "POG_CSA14_25ns_v1" in wp or "POG_CSA14_50ns_v1" in wp:
showerShapes = "full5x5"
vars = {
'dEtaIn' : abs(self.physObj.deltaEtaSuperClusterTrackAtVtx()),
'dPhiIn' : abs(self.physObj.deltaPhiSuperClusterTrackAtVtx()),
'sigmaIEtaIEta' : self.physObj.full5x5_sigmaIetaIeta() if showerShapes == "full5x5" else self.physObj.sigmaIetaIeta(),
'H/E' : self.physObj.hadronicOverEm(),
#'1/E-1/p' : abs(1.0/self.physObj.ecalEnergy() - self.physObj.eSuperClusterOverP()/self.physObj.ecalEnergy()),
'1/E-1/p' : abs(1.0/self.physObj.ecalEnergy() - self.physObj.eSuperClusterOverP()/self.physObj.ecalEnergy()) if self.physObj.ecalEnergy()>0. else 9e9,
}
WP = {
## ------- https://twiki.cern.ch/twiki/bin/viewauth/CMS/EgammaCutBasedIdentification?rev=31
'POG_2012_Veto' : [('dEtaIn', [0.007, 0.01]), ('dPhiIn', [0.8, 0.7 ]), ('sigmaIEtaIEta', [0.01, 0.03]), ('H/E', [0.15, 9e9]), ('1/E-1/p', [9e9, 9e9])],
'POG_2012_Loose' : [('dEtaIn', [0.007, 0.009]), ('dPhiIn', [0.15, 0.1 ]), ('sigmaIEtaIEta', [0.01, 0.03]), ('H/E', [0.12, 0.1]), ('1/E-1/p', [0.05, 0.05])],
'POG_2012_Medium' : [('dEtaIn', [0.004, 0.007]), ('dPhiIn', [0.06, 0.03]), ('sigmaIEtaIEta', [0.01, 0.03]), ('H/E', [0.12, 0.1]), ('1/E-1/p', [0.05, 0.05])],
'POG_2012_Tight' : [('dEtaIn', [0.004, 0.005]), ('dPhiIn', [0.03, 0.02]), ('sigmaIEtaIEta', [0.01, 0.03]), ('H/E', [0.12, 0.1]), ('1/E-1/p', [0.05, 0.05])],
## ------- https://indico.cern.cH/Event/298242/contribution/1/material/slides/5.pdf (slide 5)
'POG_CSA14_25ns_v1_Veto' : [('dEtaIn', [0.012, 0.015]), ('dPhiIn', [0.8, 0.7 ]), ('sigmaIEtaIEta', [0.01 , 0.033]), ('H/E', [0.15, 9e9 ]), ('1/E-1/p', [9e9, 9e9])],
'POG_CSA14_25ns_v1_Loose' : [('dEtaIn', [0.012, 0.014]), ('dPhiIn', [0.15, 0.1 ]), ('sigmaIEtaIEta', [0.01 , 0.033]), ('H/E', [0.12, 0.12]), ('1/E-1/p', [0.05, 0.05])],
'POG_CSA14_25ns_v1_Medium' : [('dEtaIn', [0.009, 0.012]), ('dPhiIn', [0.06, 0.03]), ('sigmaIEtaIEta', [0.01 , 0.031]), ('H/E', [0.12, 0.12]), ('1/E-1/p', [0.05, 0.05])],
'POG_CSA14_25ns_v1_Tight' : [('dEtaIn', [0.009, 0.010]), ('dPhiIn', [0.03, 0.02]), ('sigmaIEtaIEta', [0.01 , 0.031]), ('H/E', [0.12, 0.12]), ('1/E-1/p', [0.05, 0.05])],
'POG_CSA14_50ns_v1_Veto' : [('dEtaIn', [0.012, 0.022]), ('dPhiIn', [0.8, 0.7 ]), ('sigmaIEtaIEta', [0.012, 0.033]), ('H/E', [0.15, 9e9 ]), ('1/E-1/p', [9e9, 9e9])],
'POG_CSA14_50ns_v1_Loose' : [('dEtaIn', [0.012, 0.021]), ('dPhiIn', [0.15, 0.1 ]), ('sigmaIEtaIEta', [0.012, 0.033]), ('H/E', [0.12, 0.12]), ('1/E-1/p', [0.05, 0.05])],
'POG_CSA14_50ns_v1_Medium' : [('dEtaIn', [0.009, 0.019]), ('dPhiIn', [0.06, 0.03]), ('sigmaIEtaIEta', [0.01 , 0.031]), ('H/E', [0.12, 0.12]), ('1/E-1/p', [0.05, 0.05])],
'POG_CSA14_50ns_v1_Tight' : [('dEtaIn', [0.009, 0.017]), ('dPhiIn', [0.03, 0.02]), ('sigmaIEtaIEta', [0.01 , 0.031]), ('H/E', [0.12, 0.12]), ('1/E-1/p', [0.05, 0.05])],
}
if wp not in WP:
raise RuntimeError, "Working point '%s' not yet implemented in Electron.py" % wp
for (cut_name,(cut_eb,cut_ee)) in WP[wp]:
if vars[cut_name] >= (cut_eb if self.physObj.isEB() else cut_ee):
return False
return True

def absEffAreaIso(self,rho,effectiveAreas):
'''MIKE, missing doc.
Should have the same name as the function in the mother class.
Can call the mother class function with super.
'''
return self.absIsoFromEA(rho,self.superCluster().eta(),effectiveAreas.eGamma)

def mvaId( self ):
return self.mvaNonTrigV0()

def tightId( self ):
return self.tightIdResult

def mvaNonTrigV0( self, full5x5=False, debug = False ):
if self._mvaNonTrigV0[full5x5] == None:
if self.associatedVertex == None: raise RuntimeError, "You need to set electron.associatedVertex before calling any MVA"
if self.rho == None: raise RuntimeError, "You need to set electron.rho before calling any MVA"
self._mvaNonTrigV0[full5x5] = ElectronMVAID_NonTrig(self.physObj, self.associatedVertex, self.rho, full5x5, debug)
return self._mvaNonTrigV0[full5x5]

def mvaTrigV0( self, full5x5=False, debug = False ):
if self._mvaTrigV0[full5x5] == None:
if self.associatedVertex == None: raise RuntimeError, "You need to set electron.associatedVertex before calling any MVA"
if self.rho == None: raise RuntimeError, "You need to set electron.rho before calling any MVA"
self._mvaTrigV0[full5x5] = ElectronMVAID_Trig(self.physObj, self.associatedVertex, self.rho, full5x5, debug)
return self._mvaTrigV0[full5x5]

def mvaTrigNoIPV0( self, full5x5=False, debug = False ):
if self._mvaTrigNoIPV0[full5x5] == None:
if self.associatedVertex == None: raise RuntimeError, "You need to set electron.associatedVertex before calling any MVA"
if self.rho == None: raise RuntimeError, "You need to set electron.rho before calling any MVA"
self._mvaTrigNoIPV0[full5x5] = ElectronMVAID_TrigNoIP(self.physObj, self.associatedVertex, self.rho, full5x5, debug)
return self._mvaTrigNoIPV0[full5x5]


def mvaIDTight(self, full5x5=False):
eta = abs(self.superCluster().eta())
if self.pt() < 20:
if (eta < 0.8) : return self.mvaTrigV0(full5x5) > +0.00;
elif (eta < 1.479): return self.mvaTrigV0(full5x5) > +0.10;
else : return self.mvaTrigV0(full5x5) > +0.62;
else:
if (eta < 0.8) : return self.mvaTrigV0(full5x5) > +0.94;
elif (eta < 1.479): return self.mvaTrigV0(full5x5) > +0.85;
else : return self.mvaTrigV0(full5x5) > +0.92;

def mvaIDLoose(self, full5x5=False):
eta = abs(self.superCluster().eta())
if self.pt() < 10:
if (eta < 0.8) : return self.mvaNonTrigV0(full5x5) > +0.47;
elif (eta < 1.479): return self.mvaNonTrigV0(full5x5) > +0.004;
else : return self.mvaNonTrigV0(full5x5) > +0.295;
else:
if (eta < 0.8) : return self.mvaNonTrigV0(full5x5) > -0.34;
elif (eta < 1.479): return self.mvaNonTrigV0(full5x5) > -0.65;
else : return self.mvaNonTrigV0(full5x5) > +0.60;

def mvaIDZZ(self):
return self.mvaIDLoose() and (self.gsfTrack().trackerExpectedHitsInner().numberOfLostHits()<=1)

def chargedHadronIsoR(self,R=0.4):
if R == 0.3: return self.physObj.pfIsolationVariables().sumChargedHadronPt
elif R == 0.4: return self.physObj.chargedHadronIso()
raise RuntimeError, "Electron chargedHadronIso missing for R=%s" % R

def neutralHadronIsoR(self,R=0.4):
if R == 0.3: return self.physObj.pfIsolationVariables().sumNeutralHadronEt
elif R == 0.4: return self.physObj.neutralHadronIso()
raise RuntimeError, "Electron neutralHadronIso missing for R=%s" % R

def photonIsoR(self,R=0.4):
if R == 0.3: return self.physObj.pfIsolationVariables().sumPhotonEt
elif R == 0.4: return self.physObj.photonIso()
raise RuntimeError, "Electron photonIso missing for R=%s" % R

def chargedAllIso(self,R=0.4):
if R == 0.3: return self.physObj.pfIsolationVariables().sumChargedParticlePt
raise RuntimeError, "Electron chargedAllIso missing for R=%s" % R

def puChargedHadronIsoR(self,R=0.4):
if R == 0.3: return self.physObj.pfIsolationVariables().sumPUPt
elif R == 0.4: return self.physObj.puChargedHadronIso()
raise RuntimeError, "Electron chargedHadronIso missing for R=%s" % R




def chargedAllIso(self):
'''This function is used in the isolation, see Lepton class.
Here, we replace the all charged isolation by the all charged isolation with cone veto'''
return self.chargedAllIsoWithConeVeto()


def dxy(self, vertex=None):
'''Returns dxy.
Computed using vertex (or self.associatedVertex if vertex not specified),
and the gsf track.
'''
if vertex is None:
vertex = self.associatedVertex
return self.gsfTrack().dxy( vertex.position() )
def p4(self):
return ROOT.reco.Candidate.p4(self.physObj)

# def p4(self):
# return self.physObj.p4(self.physObj.candidateP4Kind()) # if kind == None else kind)

def dz(self, vertex=None):
'''Returns dz.
Computed using vertex (or self.associatedVertex if vertex not specified),
and the gsf track.
'''
if vertex is None:
vertex = self.associatedVertex
return self.gsfTrack().dz( vertex.position() )

def lostInner(self) :
if hasattr(self.gsfTrack(),"trackerExpectedHitsInner") :
return self.gsfTrack().trackerExpectedHitsInner().numberOfLostHits()
else :
return self.gsfTrack().hitPattern().numberOfHits(ROOT.reco.HitPattern.MISSING_INNER_HITS)

0 comments on commit 2440eb3

Please sign in to comment.