analyzers/objects/JetAnalyzer.py

0001 from __future__ import print_function
0002 from builtins import range
0003 import math, os
0004 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
0005 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
0006 from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import Jet
0007 from PhysicsTools.HeppyCore.utils.deltar import deltaR2, deltaPhi, matchObjectCollection, matchObjectCollection2, bestMatch,matchObjectCollection3
0008 from PhysicsTools.Heppy.physicsutils.JetReCalibrator import JetReCalibrator
0009 import PhysicsTools.HeppyCore.framework.config as cfg
0010
0011 from PhysicsTools.Heppy.physicsutils.QGLikelihoodCalculator import QGLikelihoodCalculator
0012
0013 import copy
0014 def cleanNearestJetOnly(jets,leptons,deltaR):
0015     dr2 = deltaR**2
0016     good = [ True for j in jets ]
0017     for l in leptons:
0018         ibest, d2m = -1, dr2
0019         for i,j in enumerate(jets):
0020             d2i = deltaR2(l.eta(),l.phi(), j.eta(),j.phi())
0021             if d2i < d2m:
0022                 ibest, d2m = i, d2i
0023         if ibest != -1: good[ibest] = False
0024     return [ j for (i,j) in enumerate(jets) if good[i] == True ]
0025
0026 def cleanJetsAndLeptons(jets,leptons,deltaR,arbitration):
0027     dr2 = deltaR**2
0028     goodjet = [ True for j in jets ]
0029     goodlep = [ True for l in leptons ]
0030     for il, l in enumerate(leptons):
0031         ibest, d2m = -1, dr2
0032         for i,j in enumerate(jets):
0033             d2i = deltaR2(l.eta(),l.phi(), j.eta(),j.phi())
0034             if d2i < dr2:
0035                 choice = arbitration(j,l)
0036                 if choice == j:
0037                     # if the two match, and we prefer the jet, then drop the lepton and be done
0038                     goodlep[il] = False
0039                     break
0040                 elif choice == (j,l) or choice == (l,j):
0041                     # asked to keep both, so we don't consider this match
0042                     continue
0043             if d2i < d2m:
0044                 ibest, d2m = i, d2i
0045         # this lepton has been killed by a jet, then we clean the jet that best matches it
0046         if not goodlep[il]: continue
0047         if ibest != -1: goodjet[ibest] = False
0048     return ( [ j for (i ,j) in enumerate(jets)    if goodjet[i ] == True ],
0049              [ l for (il,l) in enumerate(leptons) if goodlep[il] == True ] )
0050
0051
0052
0053 def shiftJERfactor(JERShift, aeta):
0054     factor = 1.079 + JERShift*0.026
0055     if   aeta > 3.2: factor = 1.056 + JERShift * 0.191
0056     elif aeta > 2.8: factor = 1.395 + JERShift * 0.063
0057     elif aeta > 2.3: factor = 1.254 + JERShift * 0.062
0058     elif aeta > 1.7: factor = 1.208 + JERShift * 0.046
0059     elif aeta > 1.1: factor = 1.121 + JERShift * 0.029
0060     elif aeta > 0.5: factor = 1.099 + JERShift * 0.028
0061     return factor
0062
0063
0064
0065
0066
0067 class JetAnalyzer( Analyzer ):
0068     """Taken from RootTools.JetAnalyzer, simplified, modified, added corrections    """
0069     def __init__(self, cfg_ana, cfg_comp, looperName):
0070         super(JetAnalyzer,self).__init__(cfg_ana, cfg_comp, looperName)
0071         mcGT   = cfg_ana.mcGT   if hasattr(cfg_ana,'mcGT')   else "PHYS14_25_V2"
0072         dataGT = cfg_ana.dataGT if hasattr(cfg_ana,'dataGT') else "GR_70_V2_AN1"
0073         self.shiftJEC = self.cfg_ana.shiftJEC if hasattr(self.cfg_ana, 'shiftJEC') else 0
0074         self.recalibrateJets = self.cfg_ana.recalibrateJets
0075         self.addJECShifts = self.cfg_ana.addJECShifts if hasattr(self.cfg_ana, 'addJECShifts') else 0
0076         if   self.recalibrateJets == "MC"  : self.recalibrateJets =     self.cfg_comp.isMC
0077         elif self.recalibrateJets == "Data": self.recalibrateJets = not self.cfg_comp.isMC
0078         elif self.recalibrateJets not in [True,False]: raise RuntimeError("recalibrateJets must be any of { True, False, 'MC', 'Data' }, while it is %r " % self.recalibrateJets)
0079
0080         calculateSeparateCorrections = getattr(cfg_ana,"calculateSeparateCorrections", False);
0081         calculateType1METCorrection  = getattr(cfg_ana,"calculateType1METCorrection",  False);
0082         self.doJEC = self.recalibrateJets or (self.shiftJEC != 0) or self.addJECShifts or calculateSeparateCorrections or calculateType1METCorrection
0083         if self.doJEC:
0084             doResidual = getattr(cfg_ana, 'applyL2L3Residual', 'Data')
0085             if   doResidual == "MC":   doResidual = self.cfg_comp.isMC
0086             elif doResidual == "Data": doResidual = not self.cfg_comp.isMC
0087             elif doResidual not in [True,False]: raise RuntimeError("If specified, applyL2L3Residual must be any of { True, False, 'MC', 'Data'(default)}")
0088             GT = getattr(cfg_comp, 'jecGT', mcGT if self.cfg_comp.isMC else dataGT)
0089             # Now take care of the optional arguments
0090             kwargs = { 'calculateSeparateCorrections':calculateSeparateCorrections,
0091                        'calculateType1METCorrection' :calculateType1METCorrection, }
0092             if kwargs['calculateType1METCorrection']: kwargs['type1METParams'] = cfg_ana.type1METParams
0093             # instantiate the jet re-calibrator
0094             self.jetReCalibrator = JetReCalibrator(GT, cfg_ana.recalibrationType, doResidual, cfg_ana.jecPath, **kwargs)
0095         self.doPuId = getattr(self.cfg_ana, 'doPuId', True)
0096         self.jetLepDR = getattr(self.cfg_ana, 'jetLepDR', 0.4)
0097         self.jetLepArbitration = getattr(self.cfg_ana, 'jetLepArbitration', lambda jet,lepton: lepton)
0098         self.lepPtMin = getattr(self.cfg_ana, 'minLepPt', -1)
0099         self.lepSelCut = getattr(self.cfg_ana, 'lepSelCut', lambda lep : True)
0100         self.jetGammaDR =  getattr(self.cfg_ana, 'jetGammaDR', 0.4)
0101         self.jetGammaLepDR =  getattr(self.cfg_ana, 'jetGammaLepDR', 0.4)
0102         self.cleanFromLepAndGammaSimultaneously = getattr(self.cfg_ana, 'cleanFromLepAndGammaSimultaneously', False)
0103         if self.cleanFromLepAndGammaSimultaneously:
0104             if hasattr(self.cfg_ana, 'jetGammaLepDR'):
0105                 self.jetGammaLepDR =  self.jetGammaLepDR
0106             elif (self.jetGammaDR == self.jetLepDR):
0107                 self.jetGammaLepDR = self.jetGammaDR
0108             else:
0109                 raise RuntimeError("DR for simultaneous cleaning of jets from leptons and photons is not defined, and dR(gamma, jet)!=dR(lep, jet)")
0110         if(self.cfg_ana.doQG):
0111             qgdefname="{CMSSW_BASE}/src/PhysicsTools/Heppy/data/pdfQG_AK4chs_13TeV_v2b.root"
0112             self.qglcalc = QGLikelihoodCalculator(getattr(self.cfg_ana,"QGpath",qgdefname).format(CMSSW_BASE= os.environ['CMSSW_BASE']))
0113         if not hasattr(self.cfg_ana ,"collectionPostFix"):self.cfg_ana.collectionPostFix=""
0114
0115     def declareHandles(self):
0116         super(JetAnalyzer, self).declareHandles()
0117         self.handles['jets']   = AutoHandle( self.cfg_ana.jetCol, 'std::vector<pat::Jet>' )
0118         self.handles['genJet'] = AutoHandle( self.cfg_ana.genJetCol, 'vector<reco::GenJet>' )
0119         self.shiftJER = self.cfg_ana.shiftJER if hasattr(self.cfg_ana, 'shiftJER') else 0
0120         self.addJERShifts = self.cfg_ana.addJERShifts if hasattr(self.cfg_ana, 'addJERShifts') else 0
0121         self.handles['rho'] = AutoHandle( self.cfg_ana.rho, 'double' )
0122
0123     def beginLoop(self, setup):
0124         super(JetAnalyzer,self).beginLoop(setup)
0125
0126     def process(self, event):
0127         self.readCollections( event.input )
0128         rho  = float(self.handles['rho'].product()[0])
0129         self.rho = rho
0130
0131         ## Read jets, if necessary recalibrate and shift MET
0132         if self.cfg_ana.copyJetsByValue:
0133             import ROOT
0134             #from ROOT.heppy import JetUtils
0135             allJets = map(lambda j:Jet(ROOT.heppy.JetUtils.copyJet(j)), self.handles['jets'].product())  #copy-by-value is safe if JetAnalyzer is ran more than once
0136         else:
0137             allJets = map(Jet, self.handles['jets'].product())
0138
0139         #set dummy MC flavour for all jets in case we want to ntuplize discarded jets later
0140         for jet in allJets:
0141             jet.mcFlavour = 0
0142
0143         self.deltaMetFromJEC = [0.,0.]
0144         self.type1METCorr    = [0.,0.,0.]
0145 #        print "before. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
0146         if self.doJEC:
0147             if not self.recalibrateJets:  # check point that things won't change
0148                 jetsBefore = [ (j.pt(),j.eta(),j.phi(),j.rawFactor()) for j in allJets ]
0149             self.jetReCalibrator.correctAll(allJets, rho, delta=self.shiftJEC,
0150                                                 addCorr=True, addShifts=self.addJECShifts,
0151                                                 metShift=self.deltaMetFromJEC, type1METCorr=self.type1METCorr )
0152             if not self.recalibrateJets:
0153                 jetsAfter = [ (j.pt(),j.eta(),j.phi(),j.rawFactor()) for j in allJets ]
0154                 if len(jetsBefore) != len(jetsAfter):
0155                     print("ERROR: I had to recompute jet corrections, and they rejected some of the jets:\nold = %s\n new = %s\n" % (jetsBefore,jetsAfter))
0156                 else:
0157                     for (told, tnew) in zip(jetsBefore,jetsAfter):
0158                         if (deltaR2(told[1],told[2],tnew[1],tnew[2])) > 0.0001:
0159                             print("ERROR: I had to recompute jet corrections, and one jet direction moved: old = %s, new = %s\n" % (told, tnew))
0160                         elif abs(told[0]-tnew[0])/(told[0]+tnew[0]) > 0.5e-3 or abs(told[3]-tnew[3]) > 0.5e-3:
0161                             print("ERROR: I had to recompute jet corrections, and one jet pt or corr changed: old = %s, new = %s\n" % (told, tnew))
0162         self.allJetsUsedForMET = allJets
0163 #        print "after. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
0164
0165         if self.cfg_comp.isMC:
0166             self.genJets = [ x for x in self.handles['genJet'].product() ]
0167             if self.cfg_ana.do_mc_match:
0168                 for igj, gj in enumerate(self.genJets):
0169                     gj.index = igj
0170 #                self.matchJets(event, allJets)
0171                 self.matchJets(event, [ j for j in allJets if j.pt()>self.cfg_ana.jetPt ]) # To match only jets above chosen threshold
0172             if getattr(self.cfg_ana, 'smearJets', False):
0173                 self.smearJets(event, allJets)
0174
0175
0176
0177
0178         ##Sort Jets by pT
0179         allJets.sort(key = lambda j : j.pt(), reverse = True)
0180
0181         leptons = []
0182         if hasattr(event, 'selectedLeptons'):
0183             leptons = [ l for l in event.selectedLeptons if l.pt() > self.lepPtMin and self.lepSelCut(l) ]
0184         if self.cfg_ana.cleanJetsFromTaus and hasattr(event, 'selectedTaus'):
0185             leptons = leptons[:] + event.selectedTaus
0186         if self.cfg_ana.cleanJetsFromIsoTracks and hasattr(event, 'selectedIsoCleanTrack'):
0187             leptons = leptons[:] + event.selectedIsoCleanTrack
0188
0189         ## Apply jet selection
0190         self.jets = []
0191         self.jetsFailId = []
0192         self.jetsAllNoID = []
0193         self.jetsIdOnly = []
0194         for jet in allJets:
0195             #Check if lepton and jet have overlapping PF candidates
0196             leps_with_overlaps = []
0197             if getattr(self.cfg_ana, 'checkLeptonPFOverlap', True):
0198                 for i in range(jet.numberOfSourceCandidatePtrs()):
0199                     p1 = jet.sourceCandidatePtr(i) #Ptr<Candidate> p1
0200                     for lep in leptons:
0201                         for j in range(lep.numberOfSourceCandidatePtrs()):
0202                             p2 = lep.sourceCandidatePtr(j)
0203                             has_overlaps = p1.key() == p2.key() and p1.refCore().id().productIndex() == p2.refCore().id().productIndex() and p1.refCore().id().processIndex() == p2.refCore().id().processIndex()
0204                             if has_overlaps:
0205                                 leps_with_overlaps += [lep]
0206             if len(leps_with_overlaps)>0:
0207                 for lep in leps_with_overlaps:
0208                     lep.jetOverlap = jet
0209             if self.testJetNoID( jet ):
0210                 self.jetsAllNoID.append(jet)
0211                 if(self.cfg_ana.doQG):
0212                     jet.qgl_calc =  self.qglcalc.computeQGLikelihood
0213                     jet.qgl_rho =  rho
0214                 if self.testJetID( jet ):
0215                     self.jets.append(jet)
0216                     self.jetsIdOnly.append(jet)
0217                 else:
0218                     self.jetsFailId.append(jet)
0219             elif self.testJetID (jet ):
0220                 self.jetsIdOnly.append(jet)
0221
0222         jetsEtaCut = [j for j in self.jets if abs(j.eta()) <  self.cfg_ana.jetEta ]
0223         self.cleanJetsAll, cleanLeptons = cleanJetsAndLeptons(jetsEtaCut, leptons, self.jetLepDR, self.jetLepArbitration)
0224
0225         self.cleanJets    = [j for j in self.cleanJetsAll if abs(j.eta()) <  self.cfg_ana.jetEtaCentral ]
0226         self.cleanJetsFwd = [j for j in self.cleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
0227         self.discardedJets = [j for j in self.jets if j not in self.cleanJetsAll]
0228         if hasattr(event, 'selectedLeptons') and self.cfg_ana.cleanSelectedLeptons:
0229             event.discardedLeptons = [ l for l in leptons if l not in cleanLeptons ]
0230             event.selectedLeptons  = [ l for l in event.selectedLeptons if l not in event.discardedLeptons ]
0231         for lep in leptons:
0232             if hasattr(lep, "jetOverlap"):
0233                 if lep.jetOverlap in self.cleanJetsAll:
0234                     #print "overlap reco", lep.p4().pt(), lep.p4().eta(), lep.p4().phi(), lep.jetOverlap.p4().pt(), lep.jetOverlap.p4().eta(), lep.jetOverlap.p4().phi()
0235                     lep.jetOverlapIdx = self.cleanJetsAll.index(lep.jetOverlap)
0236                 elif lep.jetOverlap in self.discardedJets:
0237                     #print "overlap discarded", lep.p4().pt(), lep.p4().eta(), lep.p4().phi(), lep.jetOverlap.p4().pt(), lep.jetOverlap.p4().eta(), lep.jetOverlap.p4().phi()
0238                     lep.jetOverlapIdx = 1000 + self.discardedJets.index(lep.jetOverlap)
0239
0240         ## First cleaning, then Jet Id
0241         self.noIdCleanJetsAll, cleanLeptons = cleanJetsAndLeptons(self.jetsAllNoID, leptons, self.jetLepDR, self.jetLepArbitration)
0242         self.noIdCleanJets = [j for j in self.noIdCleanJetsAll if abs(j.eta()) <  self.cfg_ana.jetEtaCentral ]
0243         self.noIdCleanJetsFwd = [j for j in self.noIdCleanJetsAll if abs(j.eta()) >=  self.cfg_ana.jetEtaCentral ]
0244         self.noIdDiscardedJets = [j for j in self.jetsAllNoID if j not in self.noIdCleanJetsAll]
0245
0246         ## Clean Jets from photons (first cleaning, then Jet Id)
0247         photons = []
0248         if hasattr(event, 'selectedPhotons'):
0249             if self.cfg_ana.cleanJetsFromFirstPhoton:
0250                 photons = event.selectedPhotons[:1]
0251             else:
0252                 photons = [ g for g in event.selectedPhotons ]
0253
0254         self.gamma_cleanJetaAll = []
0255         self.gamma_noIdCleanJetsAll = []
0256
0257         if self.cleanFromLepAndGammaSimultaneously:
0258             self.gamma_cleanJetsAll = cleanNearestJetOnly(jetsEtaCut, photons+leptons, self.jetGammaLepDR)
0259             self.gamma_noIdCleanJetsAll = cleanNearestJetOnly(self.jetsAllNoID, photons+leptons, self.jetGammaLepDR)
0260         else:
0261             self.gamma_cleanJetsAll = cleanNearestJetOnly(self.cleanJetsAll, photons, self.jetGammaDR)
0262             self.gamma_noIdCleanJetsAll = cleanNearestJetOnly(self.noIdCleanJetsAll, photons, self.jetGammaDR)
0263
0264         self.gamma_cleanJets    = [j for j in self.gamma_cleanJetsAll if abs(j.eta()) <  self.cfg_ana.jetEtaCentral ]
0265         self.gamma_cleanJetsFwd = [j for j in self.gamma_cleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
0266
0267         self.gamma_noIdCleanJets    = [j for j in self.gamma_noIdCleanJetsAll if abs(j.eta()) <  self.cfg_ana.jetEtaCentral ]
0268         self.gamma_noIdCleanJetsFwd = [j for j in self.gamma_noIdCleanJetsAll if abs(j.eta()) >= self.cfg_ana.jetEtaCentral ]
0269         ###
0270
0271         if self.cfg_ana.alwaysCleanPhotons:
0272             self.cleanJets = self.gamma_cleanJets
0273             self.cleanJetsAll = self.gamma_cleanJetsAll
0274             self.cleanJetsFwd = self.gamma_cleanJetsFwd
0275             #
0276             self.noIdCleanJets = self.gamma_noIdCleanJets
0277             self.noIdCleanJetsAll = self.gamma_noIdCleanJetsAll
0278             self.noIdCleanJetsFwd = self.gamma_noIdCleanJetsFwd
0279
0280         ## Jet Id, after jet/lepton cleaning
0281         self.cleanJetsFailIdAll = []
0282         for jet in self.noIdCleanJetsAll:
0283             if not self.testJetID( jet ):
0284                 self.cleanJetsFailIdAll.append(jet)
0285
0286         self.cleanJetsFailId = [j for j in self.cleanJetsFailIdAll if abs(j.eta()) <  self.cfg_ana.jetEtaCentral ]
0287
0288         ## Jet Id, after jet/photon cleaning
0289         self.gamma_cleanJetsFailIdAll = []
0290         for jet in self.gamma_noIdCleanJetsAll:
0291             if not self.testJetID( jet ):
0292                 self.gamma_cleanJetsFailIdAll.append(jet)
0293
0294         self.gamma_cleanJetsFailId = [j for j in self.gamma_cleanJetsFailIdAll if abs(j.eta()) <  self.cfg_ana.jetEtaCentral ]
0295
0296         ## Associate jets to leptons
0297         incleptons = event.inclusiveLeptons if hasattr(event, 'inclusiveLeptons') else event.selectedLeptons
0298         jlpairs = matchObjectCollection(incleptons, allJets, self.jetLepDR**2)
0299
0300         for jet in allJets:
0301             jet.leptons = [l for l in jlpairs if jlpairs[l] == jet ]
0302         for lep in incleptons:
0303             jet = jlpairs[lep]
0304             if jet is None:
0305                 setattr(lep,"jet"+self.cfg_ana.collectionPostFix,lep)
0306             else:
0307                 setattr(lep,"jet"+self.cfg_ana.collectionPostFix,jet)
0308         ## Associate jets to taus
0309         taus = getattr(event,'selectedTaus',[])
0310         jtaupairs = matchObjectCollection( taus, allJets, self.jetLepDR**2)
0311
0312         for jet in allJets:
0313             jet.taus = [l for l in jtaupairs if jtaupairs[l] == jet ]
0314         for tau in taus:
0315             setattr(tau,"jet"+self.cfg_ana.collectionPostFix,jtaupairs[tau])
0316
0317         #MC stuff
0318         if self.cfg_comp.isMC:
0319             self.deltaMetFromJetSmearing = [0, 0]
0320             for j in self.cleanJetsAll:
0321                 if hasattr(j, 'deltaMetFromJetSmearing'):
0322                     self.deltaMetFromJetSmearing[0] += j.deltaMetFromJetSmearing[0]
0323                     self.deltaMetFromJetSmearing[1] += j.deltaMetFromJetSmearing[1]
0324
0325             self.cleanGenJets = cleanNearestJetOnly(self.genJets, leptons, self.jetLepDR)
0326
0327             if self.cfg_ana.cleanGenJetsFromPhoton:
0328                 self.cleanGenJets = cleanNearestJetOnly(self.cleanGenJets, photons, self.jetLepDR)
0329
0330             if getattr(self.cfg_ana, 'attachNeutrinos', True) and hasattr(self.cfg_ana,"genNuSelection") :
0331                 jetNus=[x for x in event.genParticles if abs(x.pdgId()) in [12,14,16] and self.cfg_ana.genNuSelection(x) ]
0332                 pairs= matchObjectCollection (jetNus, self.genJets, 0.4**2)
0333
0334                 for (nu,genJet) in pairs.items() :
0335                     if genJet is not None :
0336                         if not hasattr(genJet,"nu") :
0337                             genJet.nu=nu.p4()
0338                         else :
0339                             genJet.nu+=nu.p4()
0340
0341
0342             if self.cfg_ana.do_mc_match:
0343                 self.jetFlavour(event)
0344
0345         if hasattr(event,"jets"+self.cfg_ana.collectionPostFix): raise RuntimeError("Event already contains a jet collection with the following postfix: "+self.cfg_ana.collectionPostFix)
0346         setattr(event,"rho"                    +self.cfg_ana.collectionPostFix, self.rho                    )
0347         setattr(event,"deltaMetFromJEC"        +self.cfg_ana.collectionPostFix, self.deltaMetFromJEC        )
0348         setattr(event,"type1METCorr"           +self.cfg_ana.collectionPostFix, self.type1METCorr           )
0349         setattr(event,"allJetsUsedForMET"      +self.cfg_ana.collectionPostFix, self.allJetsUsedForMET      )
0350         setattr(event,"jets"                   +self.cfg_ana.collectionPostFix, self.jets                   )
0351         setattr(event,"jetsFailId"             +self.cfg_ana.collectionPostFix, self.jetsFailId             )
0352         setattr(event,"jetsAllNoID"            +self.cfg_ana.collectionPostFix, self.jetsAllNoID            )
0353         setattr(event,"jetsIdOnly"             +self.cfg_ana.collectionPostFix, self.jetsIdOnly             )
0354         setattr(event,"cleanJetsAll"           +self.cfg_ana.collectionPostFix, self.cleanJetsAll           )
0355         setattr(event,"cleanJets"              +self.cfg_ana.collectionPostFix, self.cleanJets              )
0356         setattr(event,"cleanJetsFwd"           +self.cfg_ana.collectionPostFix, self.cleanJetsFwd           )
0357         setattr(event,"cleanJetsFailIdAll"           +self.cfg_ana.collectionPostFix, self.cleanJetsFailIdAll           )
0358         setattr(event,"cleanJetsFailId"              +self.cfg_ana.collectionPostFix, self.cleanJetsFailId              )
0359         setattr(event,"discardedJets"          +self.cfg_ana.collectionPostFix, self.discardedJets          )
0360         setattr(event,"gamma_cleanJetsAll"     +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsAll     )
0361         setattr(event,"gamma_cleanJets"        +self.cfg_ana.collectionPostFix, self.gamma_cleanJets        )
0362         setattr(event,"gamma_cleanJetsFwd"     +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsFwd     )
0363         setattr(event,"gamma_cleanJetsFailIdAll"     +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsFailIdAll     )
0364         setattr(event,"gamma_cleanJetsFailId"        +self.cfg_ana.collectionPostFix, self.gamma_cleanJetsFailId        )
0365
0366
0367         if self.cfg_comp.isMC:
0368             setattr(event,"deltaMetFromJetSmearing"+self.cfg_ana.collectionPostFix, self.deltaMetFromJetSmearing)
0369             setattr(event,"cleanGenJets"           +self.cfg_ana.collectionPostFix, self.cleanGenJets           )
0370             setattr(event,"genJets"                +self.cfg_ana.collectionPostFix, self.genJets                )
0371             if self.cfg_ana.do_mc_match:
0372                 setattr(event,"bqObjects"              +self.cfg_ana.collectionPostFix, self.bqObjects              )
0373                 setattr(event,"cqObjects"              +self.cfg_ana.collectionPostFix, self.cqObjects              )
0374                 setattr(event,"partons"                +self.cfg_ana.collectionPostFix, self.partons                )
0375                 setattr(event,"heaviestQCDFlavour"     +self.cfg_ana.collectionPostFix, self.heaviestQCDFlavour     )
0376
0377
0378         return True
0379
0380
0381
0382     def testJetID(self, jet):
0383         jet.puJetIdPassed = jet.puJetId()
0384         jet.pfJetIdPassed = jet.jetID('POG_PFID_Loose')
0385         if self.cfg_ana.relaxJetId:
0386             return True
0387         else:
0388             return jet.pfJetIdPassed and (jet.puJetIdPassed or not(self.doPuId))
0389
0390     def testJetNoID( self, jet ):
0391         # 2 is loose pile-up jet id
0392         return jet.pt() > self.cfg_ana.jetPt and \
0393                abs( jet.eta() ) < self.cfg_ana.jetEta;
0394
0395     def jetFlavour(self,event):
0396         def isFlavour(x,f):
0397             id = abs(x.pdgId())
0398             if id > 999: return (id/1000)%10 == f
0399             if id >  99: return  (id/100)%10 == f
0400             return id % 100 == f
0401
0402
0403
0404         self.bqObjects = [ p for p in event.genParticles if (p.status() == 2 and isFlavour(p,5)) ]
0405         self.cqObjects = [ p for p in event.genParticles if (p.status() == 2 and isFlavour(p,4)) ]
0406
0407         self.partons   = [ p for p in event.genParticles if ((p.status() == 23 or p.status() == 3) and abs(p.pdgId())>0 and (abs(p.pdgId()) in [1,2,3,4,5,21]) ) ]
0408         match = matchObjectCollection2(self.cleanJetsAll,
0409                                        self.partons,
0410                                        deltaRMax = 0.3)
0411
0412         for jet in self.cleanJetsAll:
0413             parton = match[jet]
0414             jet.partonId = (parton.pdgId() if parton != None else 0)
0415             jet.partonMotherId = (parton.mother(0).pdgId() if parton != None and parton.numberOfMothers()>0 else 0)
0416
0417         for jet in self.jets:
0418             (bmatch, dr) = bestMatch(jet, self.bqObjects)
0419             if dr < 0.4:
0420                 jet.mcFlavour = 5
0421             else:
0422                 (cmatch, dr) = bestMatch(jet, self.cqObjects)
0423                 if dr < 0.4:
0424                     jet.mcFlavour = 4
0425                 else:
0426                     jet.mcFlavour = 0
0427
0428         self.heaviestQCDFlavour = 5 if len(self.bqObjects) else (4 if len(self.cqObjects) else 1);
0429
0430     def matchJets(self, event, jets):
0431         match = matchObjectCollection2(jets,
0432                                        event.genbquarks + event.genwzquarks,
0433                                        deltaRMax = 0.3)
0434         for jet in jets:
0435             gen = match[jet]
0436             jet.mcParton    = gen
0437             jet.mcMatchId   = (gen.sourceId     if gen != None else 0)
0438             jet.mcMatchFlav = (abs(gen.pdgId()) if gen != None else 0)
0439
0440         match = matchObjectCollection2(jets,
0441                                        self.genJets,
0442                                        deltaRMax = 0.3)
0443         for jet in jets:
0444             jet.mcJet = match[jet]
0445
0446
0447
0448     def smearJets(self, event, jets):
0449         # https://twiki.cern.ch/twiki/bin/viewauth/CMS/TWikiTopRefSyst#Jet_energy_resolution
0450         for jet in jets:
0451             gen = jet.mcJet
0452             if gen != None:
0453                 genpt, jetpt, aeta = gen.pt(), jet.pt(), abs(jet.eta())
0454                 # from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
0455                 #8 TeV tables
0456                 factor = shiftJERfactor(self.shiftJER, aeta)
0457                 ptscale = max(0.0, (jetpt + (factor-1)*(jetpt-genpt))/jetpt)
0458                 #print "get with pt %.1f (gen pt %.1f, ptscale = %.3f)" % (jetpt,genpt,ptscale)
0459                 jet.deltaMetFromJetSmearing = [ -(ptscale-1)*jet.rawFactor()*jet.px(), -(ptscale-1)*jet.rawFactor()*jet.py() ]
0460                 if ptscale != 0:
0461                     jet.setP4(jet.p4()*ptscale)
0462                     # leave the uncorrected unchanged for sync
0463                     jet.setRawFactor(jet.rawFactor()/ptscale)
0464             #else: print "jet with pt %.1d, eta %.2f is unmatched" % (jet.pt(), jet.eta())
0465                 if (self.shiftJER==0) and (self.addJERShifts):
0466                     setattr(jet, "corrJER", ptscale )
0467                     factorJERUp= shiftJERfactor(1, aeta)
0468                     ptscaleJERUp = max(0.0, (jetpt + (factorJERUp-1)*(jetpt-genpt))/jetpt)
0469                     setattr(jet, "corrJERUp", ptscaleJERUp)
0470                     factorJERDown= shiftJERfactor(-1, aeta)
0471                     ptscaleJERDown = max(0.0, (jetpt + (factorJERDown-1)*(jetpt-genpt))/jetpt)
0472                     setattr(jet, "corrJERDown", ptscaleJERDown)
0473
0474
0475
0476
0477
0478 setattr(JetAnalyzer,"defaultConfig", cfg.Analyzer(
0479     class_object = JetAnalyzer,
0480     jetCol = 'slimmedJets',
0481     copyJetsByValue = False,      #Whether or not to copy the input jets or to work with references (should be 'True' if JetAnalyzer is run more than once)
0482     genJetCol = 'slimmedGenJets',
0483     rho = ('fixedGridRhoFastjetAll','',''),
0484     jetPt = 25.,
0485     jetEta = 4.7,
0486     jetEtaCentral = 2.4,
0487     jetLepDR = 0.4,
0488     jetLepArbitration = (lambda jet,lepton : lepton), # you can decide which to keep in case of overlaps; e.g. if the jet is b-tagged you might want to keep the jet
0489     cleanSelectedLeptons = True, #Whether to clean 'selectedLeptons' after disambiguation. Treat with care (= 'False') if running Jetanalyzer more than once
0490     minLepPt = 10,
0491     lepSelCut = lambda lep : True,
0492     relaxJetId = False,
0493     doPuId = False, # Not commissioned in 7.0.X
0494     doQG = False,
0495     checkLeptonPFOverlap = True,
0496     recalibrateJets = False,
0497     applyL2L3Residual = 'Data', # if recalibrateJets, apply L2L3Residual to Data only
0498     recalibrationType = "AK4PFchs",
0499     shiftJEC = 0, # set to +1 or -1 to apply +/-1 sigma shift to the nominal jet energies
0500     addJECShifts = False, # if true, add  "corr", "corrJECUp", and "corrJECDown" for each jet (requires uncertainties to be available!)
0501     smearJets = True,
0502     shiftJER = 0, # set to +1 or -1 to get +/-1 sigma shifts
0503     jecPath = "",
0504     calculateSeparateCorrections = False,
0505     calculateType1METCorrection  = False,
0506     type1METParams = { 'jetPtThreshold':15., 'skipEMfractionThreshold':0.9, 'skipMuons':True },
0507     addJERShifts = 0, # add +/-1 sigma shifts to jets, intended to be used with shiftJER=0
0508     cleanJetsFromFirstPhoton = False,
0509     cleanJetsFromTaus = False,
0510     cleanJetsFromIsoTracks = False,
0511     alwaysCleanPhotons = False,
0512     do_mc_match=True,
0513     cleanGenJetsFromPhoton = False,
0514     jetGammaDR=0.4,
0515     cleanFromLepAndGammaSimultaneously = False,
0516     jetGammaLepDR=0.4,
0517     attachNeutrinos = True,
0518     genNuSelection = lambda nu : True, #FIXME: add here check for ispromptfinalstate
0519     collectionPostFix = ""
0520     )
0521 )