File indexing completed on 2024-04-06 12:23:26
0001 import random
0002
0003 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0004 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0005 from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import Jet, GenJet
0006
0007 from PhysicsTools.HeppyCore.utils.deltar import cleanObjectCollection, matchObjectCollection
0008 from PhysicsTools.HeppyCore.statistics.counter import Counter, Counters
0009 from PhysicsTools.HeppyCore.utils.deltar import deltaR2
0010
0011
0012 from PhysicsTools.Heppy.physicsutils.BTagSF import BTagSF
0013 from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import GenParticle
0014 from PhysicsTools.Heppy.utils.cmsswRelease import isNewerThan
0015
0016 class JetAnalyzer( Analyzer ):
0017 """Analyze jets ;-)
0018
0019 This analyzer filters the jets that do not correspond to the leptons
0020 stored in event.selectedLeptons, and puts in the event:
0021 - jets: all jets passing the pt and eta cuts
0022 - cleanJets: the collection of jets away from the leptons
0023 - cleanBJets: the jets passing testBJet, and away from the leptons
0024
0025 Example configuration:
0026
0027 jetAna = cfg.Analyzer(
0028 'JetAnalyzer',
0029 jetCol = 'slimmedJets'
0030 # cmg jet input collection
0031 # pt threshold
0032 jetPt = 30,
0033 # eta range definition
0034 jetEta = 5.0,
0035 # seed for the btag scale factor
0036 btagSFseed = 0xdeadbeef,
0037 # if True, the PF and PU jet ID are not applied, and the jets get flagged
0038 relaxJetId = False,
0039 )
0040 """
0041
0042 def __init__(self, cfg_ana, cfg_comp, looperName):
0043 super(JetAnalyzer,self).__init__(cfg_ana, cfg_comp, looperName)
0044 self.btagSF = BTagSF (cfg_ana.btagSFseed)
0045 self.is2012 = isNewerThan('CMSSW_5_2_0')
0046
0047 def declareHandles(self):
0048 super(JetAnalyzer, self).declareHandles()
0049
0050 self.handles['jets'] = AutoHandle( self.cfg_ana.jetCol,
0051 'std::vector<pat::Jet>' )
0052 if self.cfg_comp.isMC:
0053
0054 self.mchandles['genParticles'] = AutoHandle( 'packedGenParticles',
0055 'std::vector<pat::PackedGenParticle>' )
0056 self.mchandles['genJets'] = AutoHandle('slimmedGenJets',
0057 'std::vector<reco::GenJet>')
0058
0059 def beginLoop(self, setup):
0060 super(JetAnalyzer,self).beginLoop(setup)
0061 self.counters.addCounter('jets')
0062 count = self.counters.counter('jets')
0063 count.register('all events')
0064 count.register('at least 2 good jets')
0065 count.register('at least 2 clean jets')
0066 count.register('at least 1 b jet')
0067 count.register('at least 2 b jets')
0068
0069 def process(self, event):
0070
0071 self.readCollections( event.input )
0072 miniaodjets = self.handles['jets'].product()
0073
0074 allJets = []
0075 event.jets = []
0076 event.bJets = []
0077 event.cleanJets = []
0078 event.cleanBJets = []
0079
0080 leptons = []
0081 if hasattr(event, 'selectedLeptons'):
0082 leptons = event.selectedLeptons
0083
0084 genJets = None
0085 if self.cfg_comp.isMC:
0086 genJets = map( GenJet, self.mchandles['genJets'].product() )
0087
0088 for maodjet in miniaodjets:
0089 jet = Jet( maodjet )
0090 allJets.append( jet )
0091 if self.cfg_comp.isMC and hasattr( self.cfg_comp, 'jetScale'):
0092 scale = random.gauss( self.cfg_comp.jetScale,
0093 self.cfg_comp.jetSmear )
0094 jet.scaleEnergy( scale )
0095 if genJets:
0096
0097 pairs = matchObjectCollection( [jet], genJets, 0.25*0.25)
0098 if pairs[jet] is None:
0099 pass
0100 else:
0101 jet.matchedGenJet = pairs[jet]
0102
0103 if self.cfg_comp.isMC and hasattr(self.cfg_ana, 'jerCorr') and self.cfg_ana.jerCorr:
0104 self.jerCorrection(jet)
0105
0106 if self.cfg_comp.isMC and hasattr(self.cfg_ana, 'jesCorr'):
0107 self.jesCorrection(jet, self.cfg_ana.jesCorr)
0108 if self.testJet( jet ):
0109 event.jets.append(jet)
0110 if self.testBJet(jet):
0111 event.bJets.append(jet)
0112
0113 self.counters.counter('jets').inc('all events')
0114
0115 event.cleanJets, dummy = cleanObjectCollection( event.jets,
0116 masks = leptons,
0117 deltaRMin = 0.5 )
0118 event.cleanBJets, dummy = cleanObjectCollection( event.bJets,
0119 masks = leptons,
0120 deltaRMin = 0.5 )
0121
0122 pairs = matchObjectCollection( leptons, allJets, 0.5*0.5)
0123
0124 for lepton in leptons:
0125 jet = pairs[lepton]
0126 if jet is None:
0127 lepton.jet = lepton
0128 else:
0129 lepton.jet = jet
0130
0131
0132 invpairs = matchObjectCollection( event.cleanJets, leptons, 99999. )
0133 for jet in event.cleanJets:
0134 leg = invpairs[jet]
0135 jet.leg = leg
0136
0137 for jet in event.cleanJets:
0138 jet.matchGenParton=999.0
0139
0140 if self.cfg_comp.isMC and "BB" in self.cfg_comp.name:
0141 genParticles = self.mchandles['genParticles'].product()
0142 event.genParticles = map( GenParticle, genParticles)
0143 for gen in genParticles:
0144 if abs(gen.pdgId())==5 and gen.mother() and abs(gen.mother().pdgId())==21:
0145 for jet in event.cleanJets:
0146 dR=deltaR2(jet.eta(), jet.phi(), gen.eta(), gen.phi() )
0147 if dR<jet.matchGenParton:
0148 jet.matchGenParton=dR
0149
0150 event.jets30 = [jet for jet in event.jets if jet.pt()>30]
0151 event.cleanJets30 = [jet for jet in event.cleanJets if jet.pt()>30]
0152 if len( event.jets30 )>=2:
0153 self.counters.counter('jets').inc('at least 2 good jets')
0154 if len( event.cleanJets30 )>=2:
0155 self.counters.counter('jets').inc('at least 2 clean jets')
0156 if len(event.cleanBJets)>0:
0157 self.counters.counter('jets').inc('at least 1 b jet')
0158 if len(event.cleanBJets)>1:
0159 self.counters.counter('jets').inc('at least 2 b jets')
0160 return True
0161
0162 def jerCorrection(self, jet):
0163 ''' Adds JER correction according to first method at
0164 https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
0165
0166 Requires some attention when genJet matching fails.
0167 '''
0168 if not hasattr(jet, 'matchedGenJet'):
0169 return
0170
0171 corrections = [0.052, 0.057, 0.096, 0.134, 0.288]
0172 maxEtas = [0.5, 1.1, 1.7, 2.3, 5.0]
0173 eta = abs(jet.eta())
0174 for i, maxEta in enumerate(maxEtas):
0175 if eta < maxEta:
0176 pt = jet.pt()
0177 deltaPt = (pt - jet.matchedGenJet.pt()) * corrections[i]
0178 totalScale = (pt + deltaPt) / pt
0179
0180 if totalScale < 0.:
0181 totalScale = 0.
0182 jet.scaleEnergy(totalScale)
0183 break
0184
0185 def jesCorrection(self, jet, scale=0.):
0186 ''' Adds JES correction in number of sigmas (scale)
0187 '''
0188
0189 if scale == 0.:
0190 return
0191 unc = jet.uncOnFourVectorScale()
0192 totalScale = 1. + scale * unc
0193 if totalScale < 0.:
0194 totalScale = 0.
0195 jet.scaleEnergy(totalScale)
0196
0197 def testJetID(self, jet):
0198 jet.puJetIdPassed = jet.puJetId()
0199 jet.pfJetIdPassed = jet.jetID("POG_PFID_Loose")
0200 if self.cfg_ana.relaxJetId:
0201 return True
0202 else:
0203 return jet.puJetIdPassed and jet.pfJetIdPassed
0204
0205
0206 def testJet( self, jet ):
0207 return jet.pt() > self.cfg_ana.jetPt and \
0208 abs( jet.eta() ) < self.cfg_ana.jetEta and \
0209 self.testJetID(jet)
0210
0211 def testBJet(self, jet):
0212
0213
0214 jet.btagMVA = jet.btag("combinedSecondaryVertexBJetTags")
0215 jet.btagFlag = self.btagSF.BTagSFcalc.isbtagged(
0216 jet.pt(),
0217 jet.eta(),
0218 jet.btag("combinedSecondaryVertexBJetTags"),
0219 abs(jet.partonFlavour()),
0220 not self.cfg_comp.isMC,
0221 0,0,
0222 self.is2012
0223 )
0224 return jet.pt()>20 and \
0225 abs( jet.eta() ) < 2.4 and \
0226 self.testJetID(jet)