Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-07-23 02:25:52 

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 )
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!