analyzers/examples/JetAnalyzer.py

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             # and ("BB" in self.cfg_comp.name):
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                 # Use DeltaR = 0.25 matching like JetMET
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             #Add JER correction for MC jets. Requires gen-jet matching.
0103             if self.cfg_comp.isMC and hasattr(self.cfg_ana, 'jerCorr') and self.cfg_ana.jerCorr:
0104                 self.jerCorrection(jet)
0105             #Add JES correction for MC jets.
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         # associating a jet to each lepton
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         # associating a leg to each clean jet
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         #import pdb; pdb.set_trace()
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         # Do nothing if nothing to change
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         # medium csv working point
0213         # https://twiki.cern.ch/twiki/bin/viewauth/CMS/BTagPerformanceOP#B_tagging_Operating_Points_for_3
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)