Back to home page

Project CMSSW displayed by LXR
 
 

    

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                     # if the two match, and we prefer the jet, then drop the lepton and be done
0037                     goodlep[il] = False
0038                     break 
0039                 elif choice == (j,l) or choice == (l,j):
0040                     # asked to keep both, so we don't consider this match
0041                     continue
0042             if d2i < d2m:
0043                 ibest, d2m = i, d2i
0044         # this lepton has been killed by a jet, then we clean the jet that best matches it
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             # Now take care of the optional arguments
0089             kwargs = { 'calculateSeparateCorrections':calculateSeparateCorrections,
0090                        'calculateType1METCorrection' :calculateType1METCorrection, }
0091             if kwargs['calculateType1METCorrection']: kwargs['type1METParams'] = cfg_ana.type1METParams
0092             # instantiate the jet re-calibrator
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         ## Read jets, if necessary recalibrate and shift MET
0131         if self.cfg_ana.copyJetsByValue: 
0132             import ROOT
0133             #from ROOT.heppy import JetUtils
0134             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
0135         else: 
0136             allJets = map(Jet, self.handles['jets'].product()) 
0137 
0138         #set dummy MC flavour for all jets in case we want to ntuplize discarded jets later
0139         for jet in allJets:
0140             jet.mcFlavour = 0
0141 
0142         self.deltaMetFromJEC = [0.,0.]
0143         self.type1METCorr    = [0.,0.,0.]
0144 #        print "before. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
0145         if self.doJEC:
0146             if not self.recalibrateJets:  # check point that things won't change
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 #        print "after. rho",self.rho,self.cfg_ana.collectionPostFix,'allJets len ',len(allJets),'pt', [j.pt() for j in allJets]
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 #                self.matchJets(event, allJets)
0170                 self.matchJets(event, [ j for j in allJets if j.pt()>self.cfg_ana.jetPt ]) # To match only jets above chosen threshold
0171             if getattr(self.cfg_ana, 'smearJets', False):
0172                 self.smearJets(event, allJets)
0173 
0174 
0175 
0176 
0177         ##Sort Jets by pT 
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         ## Apply jet selection
0189         self.jets = []
0190         self.jetsFailId = []
0191         self.jetsAllNoID = []
0192         self.jetsIdOnly = []
0193         for jet in allJets:
0194             #Check if lepton and jet have overlapping PF candidates 
0195             leps_with_overlaps = []
0196             if getattr(self.cfg_ana, 'checkLeptonPFOverlap', True):
0197                 for i in range(jet.numberOfSourceCandidatePtrs()):
0198                     p1 = jet.sourceCandidatePtr(i) #Ptr<Candidate> p1
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                     #print "overlap reco", lep.p4().pt(), lep.p4().eta(), lep.p4().phi(), lep.jetOverlap.p4().pt(), lep.jetOverlap.p4().eta(), lep.jetOverlap.p4().phi()
0234                     lep.jetOverlapIdx = self.cleanJetsAll.index(lep.jetOverlap)
0235                 elif lep.jetOverlap in self.discardedJets:
0236                     #print "overlap discarded", lep.p4().pt(), lep.p4().eta(), lep.p4().phi(), lep.jetOverlap.p4().pt(), lep.jetOverlap.p4().eta(), lep.jetOverlap.p4().phi()
0237                     lep.jetOverlapIdx = 1000 + self.discardedJets.index(lep.jetOverlap)
0238 
0239         ## First cleaning, then Jet Id
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         ## Clean Jets from photons (first cleaning, then Jet Id)
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         ## Jet Id, after jet/lepton cleaning
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         ## Jet Id, after jet/photon cleaning
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         ## Associate jets to leptons
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         ## Associate jets to taus 
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         #MC stuff
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         # 2 is loose pile-up jet id
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         # https://twiki.cern.ch/twiki/bin/viewauth/CMS/TWikiTopRefSyst#Jet_energy_resolution
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                 # from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
0454                 #8 TeV tables
0455                 factor = shiftJERfactor(self.shiftJER, aeta)
0456                 ptscale = max(0.0, (jetpt + (factor-1)*(jetpt-genpt))/jetpt)             
0457                 #print "get with pt %.1f (gen pt %.1f, ptscale = %.3f)" % (jetpt,genpt,ptscale)
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                     # leave the uncorrected unchanged for sync
0462                     jet.setRawFactor(jet.rawFactor()/ptscale)
0463             #else: print "jet with pt %.1d, eta %.2f is unmatched" % (jet.pt(), jet.eta())
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,      #Whether or not to copy the input jets or to work with references (should be 'True' if JetAnalyzer is run more than once)
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), # 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
0488     cleanSelectedLeptons = True, #Whether to clean 'selectedLeptons' after disambiguation. Treat with care (= 'False') if running Jetanalyzer more than once
0489     minLepPt = 10,
0490     lepSelCut = lambda lep : True,
0491     relaxJetId = False,  
0492     doPuId = False, # Not commissioned in 7.0.X
0493     doQG = False, 
0494     checkLeptonPFOverlap = True,
0495     recalibrateJets = False,
0496     applyL2L3Residual = 'Data', # if recalibrateJets, apply L2L3Residual to Data only
0497     recalibrationType = "AK4PFchs",
0498     shiftJEC = 0, # set to +1 or -1 to apply +/-1 sigma shift to the nominal jet energies
0499     addJECShifts = False, # if true, add  "corr", "corrJECUp", and "corrJECDown" for each jet (requires uncertainties to be available!)
0500     smearJets = True,
0501     shiftJER = 0, # set to +1 or -1 to get +/-1 sigma shifts    
0502     jecPath = "",
0503     calculateSeparateCorrections = False,
0504     calculateType1METCorrection  = False,
0505     type1METParams = { 'jetPtThreshold':15., 'skipEMfractionThreshold':0.9, 'skipMuons':True },
0506     addJERShifts = 0, # add +/-1 sigma shifts to jets, intended to be used with shiftJER=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, #FIXME: add here check for ispromptfinalstate
0518     collectionPostFix = ""
0519     )
0520 )
[ 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!