python/physicsobjects/Jet.py

0001 from builtins import range
0002 from PhysicsTools.Heppy.physicsobjects.PhysicsObject import *
0003 from PhysicsTools.HeppyCore.utils.deltar import deltaPhi
0004 import math
0005
0006 loose_WP = [
0007     (0, 2.5, -0.8),
0008     (2.5, 2.75, -0.74),
0009     (2.75, 3.0, -0.68),
0010     (3.0, 5.0, -0.77),
0011     ]
0012
0013 # Working point 2 May 2013 (Phil via H2tau list)
0014 loose_53X_WP = [
0015     (0, 2.5, -0.63),
0016     (2.5, 2.75, -0.60),
0017     (2.75, 3.0, -0.55),
0018     (3.0, 5.2, -0.45),
0019     ]
0020
0021 _btagWPs = {
0022     "TCHEL": ("pfTrackCountingHighEffBJetTags", 1.7),
0023     "TCHEM": ("pfTrackCountingHighEffBJetTags", 3.3),
0024     "TCHPT": ("pfTrackCountingHighPurBJetTags", 3.41),
0025     "JPL": ("pfJetProbabilityBJetTags", 0.275),
0026     "JPM": ("pfJetProbabilityBJetTags", 0.545),
0027     "JPT": ("pfJetProbabilityBJetTags", 0.790),
0028     "CSVL": ("combinedSecondaryVertexBJetTags", 0.244),
0029     "CSVM": ("combinedSecondaryVertexBJetTags", 0.679),
0030     "CSVT": ("combinedSecondaryVertexBJetTags", 0.898),
0031 ###https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagRecommendation74X50ns
0032     "CMVAL": ("pfCombinedMVABJetTags", 0.630), # for same b-jet efficiency of CSVv2IVFL on ttbar MC, jet pt > 30
0033     "CMVAM": ("pfCombinedMVABJetTags", 0.732), # for same b-jet efficiency of CSVv2IVFM on ttbar MC, jet pt > 30
0034     "CMVAT": ("pfCombinedMVABJetTags", 0.813), # for same b-jet efficiency of CSVv2IVFT on ttbar MC, jet pt > 30
0035     "CMVAv2M": ("pfCombinedMVAV2BJetTags", 0.185), # for same b-jet efficiency of CSVv2IVFM on ttbar MC, jet pt > 30
0036 ###https://twiki.cern.ch/twiki/bin/viewauth/CMS/BtagRecommendation80X
0037     "CSVv2IVFL": ("pfCombinedInclusiveSecondaryVertexV2BJetTags", 0.460),
0038     "CSVv2IVFM": ("pfCombinedInclusiveSecondaryVertexV2BJetTags", 0.800),
0039     "CSVv2IVFT": ("pfCombinedInclusiveSecondaryVertexV2BJetTags", 0.935),
0040     "CMVAv2L": ("pfCombinedMVAV2BJetTags", -0.715), # for same b-jet efficiency of CSVv2IVFL on ttbar MC, jet pt > 30
0041     "CMVAv2M": ("pfCombinedMVAV2BJetTags", 0.185),  # for same b-jet efficiency of CSVv2IVFM on ttbar MC, jet pt > 30
0042     "CMVAv2T": ("pfCombinedMVAV2BJetTags", 0.875),  # for same b-jet efficiency of CSVv2IVFT on ttbar MC, jet pt > 30
0043
0044 }
0045
0046 class Jet(PhysicsObject):
0047     def __init__(self, *args, **kwargs):
0048         super(Jet, self).__init__(*args, **kwargs)
0049         self._physObjInit()
0050
0051     def _physObjInit(self):
0052         self._rawFactorMultiplier = 1.0
0053         self._recalibrated = False
0054         self._leadingTrack = None
0055         self._leadingTrackSearched = False
0056
0057     def rawEnergy(self):
0058         return self.energy() * self.rawFactor()
0059
0060     # these energy fraction methods need to be redefined here
0061     # because the pat::Jet's currentJECLevel data member cannot be update easily by the calibrator
0062     # and then the C++ methods would be broken when new jet corrections are applied in Heppy
0063     def chargedEmEnergyFraction(self):
0064         return self.chargedEmEnergy()/self.rawEnergy()
0065
0066     def chargedHadronEnergyFraction(self):
0067         return self.chargedHadronEnergy()/self.rawEnergy()
0068
0069     def chargedMuEnergyFraction(self):
0070         return self.chargedMuEnergy()/self.rawEnergy()
0071
0072     def electronEnergyFraction(self):
0073         return self.electronEnergy()/self.rawEnergy()
0074
0075     def muonEnergyFraction(self):
0076         return self.muonEnergy()/self.rawEnergy()
0077
0078     def neutralEmEnergyFraction(self):
0079         return self.neutralEmEnergy()/self.rawEnergy()
0080
0081     def neutralHadronEnergyFraction(self):
0082         return self.neutralHadronEnergy()/self.rawEnergy()
0083
0084     def photonEnergyFraction(self):
0085         return self.photonEnergy()/self.rawEnergy()
0086
0087
0088     def HFHadronEnergyFraction(self):
0089         return self.HFHadronEnergy()/self.rawEnergy()
0090
0091     def HFEMEnergyFraction(self):
0092         return self.HFEMEnergy()/self.rawEnergy()
0093
0094     def hoEnergyFraction(self):
0095         return self.hoEnergy()/self.rawEnergy()
0096
0097
0098
0099     def jetID(self,name=""):
0100         if not self.isPFJet():
0101             raise RuntimeError("jetID implemented only for PF Jets")
0102         eta = abs(self.eta());
0103         energy = self.rawEnergy();
0104         chf = self.chargedHadronEnergy()/energy;
0105         nhf = self.neutralHadronEnergy()/energy;
0106         phf = self.neutralEmEnergy()/energy;
0107         muf = self.muonEnergy()/energy;
0108         elf = self.chargedEmEnergy()/energy;
0109         chm = self.chargedHadronMultiplicity();
0110         npr = self.chargedMultiplicity() + self.neutralMultiplicity();
0111         npn = self.neutralMultiplicity();
0112         #if npr != self.nConstituents():
0113         #    import pdb; pdb.set_trace()
0114         if name == "POG_PFID":
0115
0116             if   self.jetID("PAG_monoID_Tight") and self.jetID("POG_PFID_Tight") : return 5;
0117             if   self.jetID("PAG_monoID_Loose") and self.jetID("POG_PFID_Tight") : return 4;
0118             if   self.jetID("POG_PFID_Tight")  : return 3;
0119             #elif self.jetID("POG_PFID_Medium") : return 2;  commented this line because this working point doesn't exist anymore (as 12/05/15)
0120             elif self.jetID("POG_PFID_Loose")  : return 1;
0121             else                               : return 0;
0122
0123         # jetID from here: https://twiki.cern.ch/twiki/bin/viewauth/CMS/JetID#Recommendations_for_13_TeV_data
0124         if name == "POG_PFID_Loose":    return ((eta<3.0 and ((npr>1 and phf<0.99 and nhf<0.99) and (eta>2.4 or (elf<0.99 and chf>0 and chm>0)))) or (eta>3.0 and (phf<0.90 and npn>10)));
0125         if name == "POG_PFID_Medium":   return (npr>1 and phf<0.95 and nhf<0.95 and muf < 0.8) and (eta>2.4 or (elf<0.99 and chf>0 and chm>0));
0126         if name == "POG_PFID_Tight":    return ((eta<3.0 and ((npr>1 and phf<0.90 and nhf<0.90) and (eta>2.4 or (elf<0.99 and chf>0 and chm>0)))) or (eta>3.0 and (phf<0.90 and npn>10)));
0127         if name == "VBFHBB_PFID_Loose":  return (npr>1 and phf<0.99 and nhf<0.99);
0128         if name == "VBFHBB_PFID_Medium": return (npr>1 and phf<0.99 and nhf<0.99) and ((eta<=2.4 and nhf<0.9 and phf<0.9 and elf<0.99 and muf<0.99 and chf>0 and chm>0) or eta>2.4);
0129         if name == "VBFHBB_PFID_Tight":  return (npr>1 and phf<0.99 and nhf<0.99) and ((eta<=2.4 and nhf<0.9 and phf<0.9 and elf<0.70 and muf<0.70 and chf>0 and chm>0) or eta>2.4);
0130         if name == "PAG_monoID_Loose":    return (eta<3.0 and chf>0.05 and nhf<0.7 and phf<0.8);
0131         if name == "PAG_monoID_Tight":    return (eta<3.0 and chf>0.2 and nhf<0.7 and phf<0.7);
0132
0133         raise RuntimeError("jetID '%s' not supported" % name)
0134
0135     def looseJetId(self):
0136         '''PF Jet ID (loose operation point) [method provided for convenience only]'''
0137         return self.jetID("POG_PFID_Loose")
0138
0139     def puMva(self, label="pileupJetId:fullDiscriminant"):
0140         if self.hasUserFloat(label):
0141             return self.userFloat(label)
0142         return -99
0143
0144     def puJetId(self, label="pileupJetId:fullDiscriminant"):
0145         '''Full mva PU jet id'''
0146
0147         puMva = self.puMva(label)
0148         wp = loose_53X_WP
0149         eta = abs(self.eta())
0150
0151         for etamin, etamax, cut in wp:
0152             if not(eta>=etamin and eta<etamax):
0153                 continue
0154             return puMva>cut
0155         return -99
0156
0157     def rawFactor(self):
0158         return self.jecFactor('Uncorrected') * self._rawFactorMultiplier
0159
0160     def setRawFactor(self, factor):
0161         self._rawFactorMultiplier = factor/self.jecFactor('Uncorrected')
0162
0163     def corrFactor(self):
0164         return 1.0/self.rawFactor()
0165
0166     def setCorrP4(self,newP4):
0167         self._recalibrated = True
0168         corr = newP4.Pt() / (self.pt() * self.rawFactor())
0169         self.setP4(newP4);
0170         self.setRawFactor(1.0/corr);
0171
0172     def l1corrFactor(self):
0173         if hasattr(self, 'CorrFactor_L1'):
0174             return self.CorrFactor_L1
0175         if self._recalibrated:
0176             raise RuntimeError("The jet was recalibrated, but without calculateSeparateCorrections. L1 is not available")
0177         jecLevels = self.physObj.availableJECLevels()
0178         for level in jecLevels:
0179             if "L1" in level:
0180                 return self.physObj.jecFactor(level)/self.physObj.jecFactor('Uncorrected')
0181         return 1.0 # Jet does not have any L1 correction
0182
0183     def btag(self,name):
0184         ret = self.bDiscriminator(name)
0185         if ret == -1000 and name.startswith("pf"):
0186             ret = self.bDiscriminator(name[2].lower()+name[3:])
0187         return ret
0188
0189     def btagWP(self,name):
0190         global _btagWPs
0191         (disc,val) = _btagWPs[name]
0192         return self.btag(disc) > val
0193
0194     def leadingTrack(self):
0195         if self._leadingTrackSearched :
0196             return self._leadingTrack
0197         self._leadingTrackSearched = True
0198         self._leadingTrack =  max( self.daughterPtrVector() , key = lambda x : x.pt() if  x.charge()!=0 else 0. )
0199         if self._leadingTrack.charge()==0: #in case of "all neutral"
0200             self._leadingTrack = None
0201         return self._leadingTrack
0202
0203     def leadTrackPt(self):
0204         lt=self.leadingTrack()
0205         if lt :
0206             return lt.pt()
0207         else :
0208             return 0.
0209     def qgl(self) :
0210         if not hasattr(self,"qgl_value") :
0211             if hasattr(self,"qgl_rho") : #check if qgl calculator is configured
0212                 self.computeQGvars()
0213                 self.qgl_value=self.qgl_calc(self,self.qgl_rho)
0214             else :
0215                 self.qgl_value=-1. #if no qgl calculator configured
0216
0217         return self.qgl_value
0218
0219     def computeQGvars(self):
0220         #return immediately if qgvars already computed or if qgl is disabled
0221         if not hasattr(self,"qgl_rho") or getattr(self,"hasQGVvars",False) :
0222             return self
0223         self.hasQGvars = True
0224
0225         jet = self
0226         jet.mult = 0
0227         sum_weight = 0.
0228         sum_pt = 0.
0229         sum_deta = 0.
0230         sum_dphi = 0.
0231         sum_deta2 = 0.
0232         sum_detadphi = 0.
0233         sum_dphi2 = 0.
0234
0235
0236
0237         for ii in range(0, jet.numberOfDaughters()) :
0238
0239             part = jet.daughter(ii)
0240
0241             if part.charge() == 0 : # neutral particles
0242
0243                 if part.pt() < 1.: continue
0244
0245             else : # charged particles
0246
0247                 if part.trackHighPurity()==False: continue
0248                 if part.fromPV()<=1: continue
0249
0250
0251             jet.mult += 1
0252
0253             deta = part.eta() - jet.eta()
0254             dphi = deltaPhi(part.phi(), jet.phi())
0255             partPt = part.pt()
0256             weight = partPt*partPt
0257             sum_weight += weight
0258             sum_pt += partPt
0259             sum_deta += deta*weight
0260             sum_dphi += dphi*weight
0261             sum_deta2 += deta*deta*weight
0262             sum_detadphi += deta*dphi*weight
0263             sum_dphi2 += dphi*dphi*weight
0264
0265
0266
0267
0268         a = 0.
0269         b = 0.
0270         c = 0.
0271
0272         if sum_weight > 0 :
0273             jet.ptd = math.sqrt(sum_weight)/sum_pt
0274             ave_deta = sum_deta/sum_weight
0275             ave_dphi = sum_dphi/sum_weight
0276             ave_deta2 = sum_deta2/sum_weight
0277             ave_dphi2 = sum_dphi2/sum_weight
0278             a = ave_deta2 - ave_deta*ave_deta
0279             b = ave_dphi2 - ave_dphi*ave_dphi
0280             c = -(sum_detadphi/sum_weight - ave_deta*ave_dphi)
0281         else: jet.ptd = 0.
0282
0283         delta = math.sqrt(math.fabs((a-b)*(a-b)+4.*c*c))
0284
0285         if a+b-delta > 0: jet.axis2 = -math.log(math.sqrt(0.5*(a+b-delta)))
0286         else: jet.axis2 = -1.
0287         return jet
0288
0289
0290
0291 class GenJet( PhysicsObject):
0292     pass
0293