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