Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:32:27

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)