Merge pull request cms-sw#50 from arizzi/tauvhbb

Tauvhbb

PhysicsTools/Heppy/python/analyzers/objects/TauAnalyzer.py

@@ -20,54 +20,53 @@ def declareHandles(self):
         super(TauAnalyzer, self).declareHandles()
         self.handles['taus'] = AutoHandle( ('slimmedTaus',''),'std::vector<pat::Tau>')
 
-
     def beginLoop(self, setup):
         super(TauAnalyzer,self).beginLoop(setup)
         self.counters.addCounter('events')
         count = self.counters.counter('events')
         count.register('all events')
         count.register('has >=1 tau at preselection')
         count.register('has >=1 selected taus')
-        count.register('has >=1 loose taus')
-        count.register('has >=1 inclusive taus')
+        count.register('has >=1 other taus')
 
     #------------------
     # MAKE LEPTON LISTS
     #------------------
     def makeTaus(self, event):
-        event.selectedTaus = []
-        event.looseTaus = []
         event.inclusiveTaus = []
+        event.selectedTaus = []
+        event.otherTaus = []
 
         #get all
         alltaus = map( Tau, self.handles['taus'].product() )
 
-        foundTau = False
+        #make inclusive taus
         for tau in alltaus:
             tau.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
             tau.lepVeto = False
             tau.idDecayMode = tau.tauID("decayModeFinding")
             tau.idDecayModeNewDMs = tau.tauID("decayModeFindingNewDMs")
-            if hasattr(self.cfg_ana, 'decayModeID') and self.cfg_ana.decayModeID and not tau.tauID(self.cfg_ana.decayModeID):
+
+            if hasattr(self.cfg_ana, 'inclusive_decayModeID') and self.cfg_ana.inclusive_decayModeID and not tau.tauID(self.cfg_ana.inclusive_decayModeID):
                 continue
 
-            if self.cfg_ana.vetoLeptons:
+            tau.inclusive_lepVeto = False
+            if self.cfg_ana.inclusive_vetoLeptons:
                 for lep in event.selectedLeptons:
-                    if deltaR(lep.eta(), lep.phi(), tau.eta(), tau.phi()) < self.cfg_ana.leptonVetoDR:
-                        tau.lepVeto = True
-                if tau.lepVeto: continue
-            if self.cfg_ana.vetoLeptonsPOG:
-                if not tau.tauID(self.cfg_ana.tauAntiMuonID):
-                        tau.lepVeto = True
-                if not tau.tauID(self.cfg_ana.tauAntiElectronID):
-                        tau.lepVeto = True
-                if tau.lepVeto: continue
-
-            if tau.pt() < self.cfg_ana.ptMin: continue
-            if abs(tau.eta()) > self.cfg_ana.etaMax: continue
-            if abs(tau.dxy()) > self.cfg_ana.dxyMax or abs(tau.dz()) > self.cfg_ana.dzMax: continue
-
-            foundTau = True
+                    if deltaR(lep.eta(), lep.phi(), tau.eta(), tau.phi()) < self.cfg_ana.inclusive_leptonVetoDR:
+                        tau.inclusive_lepVeto = True
+                if tau.inclusive_lepVeto: continue
+            if self.cfg_ana.inclusive_vetoLeptonsPOG:
+                if not tau.tauID(self.cfg_ana.inclusive_tauAntiMuonID):
+                        tau.inclusive_lepVeto = True
+                if not tau.tauID(self.cfg_ana.inclusive_tauAntiElectronID):
+                        tau.inclusive_lepVeto = True
+                if tau.inclusive_lepVeto: continue
+
+            if tau.pt() < self.cfg_ana.inclusive_ptMin: continue
+            if abs(tau.eta()) > self.cfg_ana.inclusive_etaMax: continue
+            if abs(tau.dxy()) > self.cfg_ana.inclusive_dxyMax or abs(tau.dz()) > self.cfg_ana.inclusive_dzMax: continue
+
             def id3(tau,X):
                 """Create an integer equal to 1-2-3 for (loose,medium,tight)"""
                 return tau.tauID(X%"Loose") + tau.tauID(X%"Medium") + tau.tauID(X%"Tight")
@@ -86,21 +85,38 @@ def id6(tau,X):
             tau.idAntiMu = tau.tauID("againstMuonLoose") + tau.tauID("againstMuonTight")
             tau.idAntiE = id5(tau, "againstElectron%sMVA5")
             #print "Tau pt %5.1f: idMVA2 %d, idCI3hit %d, %s, %s" % (tau.pt(), tau.idMVA2, tau.idCI3hit, tau.tauID(self.cfg_ana.tauID), tau.tauID(self.cfg_ana.tauLooseID))
-            if tau.tauID(self.cfg_ana.tauID):
-                event.selectedTaus.append(tau)
-                event.inclusiveTaus.append(tau)
-            elif tau.tauID(self.cfg_ana.tauLooseID):
-                event.looseTaus.append(tau)
+
+            if tau.tauID(self.cfg_ana.inclusive_tauID):
                 event.inclusiveTaus.append(tau)
+
+        for tau in event.inclusiveTaus:
 
+            tau.loose_lepVeto = False
+            if self.cfg_ana.loose_vetoLeptons:
+                for lep in event.selectedLeptons:
+                    if deltaR(lep.eta(), lep.phi(), tau.eta(), tau.phi()) < self.cfg_ana.loose_leptonVetoDR:
+                        tau.loose_lepVeto = True
+            if self.cfg_ana.loose_vetoLeptonsPOG:
+                if not tau.tauID(self.cfg_ana.loose_tauAntiMuonID):
+                        tau.loose_lepVeto = True
+                if not tau.tauID(self.cfg_ana.loose_tauAntiElectronID):
+                        tau.loose_lepVeto = True
+
+            if tau.tauID(self.cfg_ana.loose_decayModeID) and \
+                      tau.pt() > self.cfg_ana.loose_ptMin and abs(tau.eta()) < self.cfg_ana.loose_etaMax and \
+                      abs(tau.dxy()) < self.cfg_ana.loose_dxyMax and abs(tau.dz()) < self.cfg_ana.loose_dzMax and \
+                      tau.tauID(self.cfg_ana.loose_tauID) and not tau.loose_lepVeto:
+                event.selectedTaus.append(tau)
+            else:
+                event.otherTaus.append(tau)
+
+        event.inclusiveTaus.sort(key = lambda l : l.pt(), reverse = True)        
         event.selectedTaus.sort(key = lambda l : l.pt(), reverse = True)
-        event.looseTaus.sort(key = lambda l : l.pt(), reverse = True)
+        event.otherTaus.sort(key = lambda l : l.pt(), reverse = True)
         self.counters.counter('events').inc('all events')
-        if foundTau: self.counters.counter('events').inc('has >=1 tau at preselection')
+        if len(event.inclusiveTaus): self.counters.counter('events').inc('has >=1 tau at preselection')
         if len(event.selectedTaus): self.counters.counter('events').inc('has >=1 selected taus')
-        if len(event.looseTaus): self.counters.counter('events').inc('has >=1 loose taus')
-        if len(event.inclusiveTaus): self.counters.counter('events').inc('has >=1 inclusive taus')
-
+        if len(event.otherTaus): self.counters.counter('events').inc('has >=1 other taus')
 
     def matchTaus(self, event):
         match = matchObjectCollection3(event.inclusiveTaus, event.gentaus, deltaRMax = 0.5)
@@ -125,18 +141,31 @@ def process(self, event):
 # http://cmslxr.fnal.gov/lxr/source/PhysicsTools/PatAlgos/python/producersLayer1/tauProducer_cfi.py
 
 setattr(TauAnalyzer,"defaultConfig",cfg.Analyzer(
-    class_object=TauAnalyzer,
-    ptMin = 20,
-    etaMax = 9999,
-    dxyMax = 1000.,
-    dzMax = 0.2,
-    vetoLeptons = True,
-    leptonVetoDR = 0.4,
-    decayModeID = "decayModeFindingNewDMs", # ignored if not set or ""
-    tauID = "byLooseCombinedIsolationDeltaBetaCorr3Hits",
-    vetoLeptonsPOG = False, # If True, the following two IDs are required
-    tauAntiMuonID = "againstMuonLoose3",
-    tauAntiElectronID = "againstElectronLooseMVA5",
-    tauLooseID = "decayModeFinding",
+    class_object = TauAnalyzer,
+    # inclusive very loose hadronic tau selection
+    inclusive_ptMin = 18,
+    inclusive_etaMax = 9999,
+    inclusive_dxyMax = 1000.,
+    inclusive_dzMax = 0.4,
+    inclusive_vetoLeptons = False,
+    inclusive_leptonVetoDR = 0.4,
+    inclusive_decayModeID = "decayModeFindingNewDMs", # ignored if not set or ""
+    inclusive_tauID = "decayModeFindingNewDMs",
+    inclusive_vetoLeptonsPOG = False, # If True, the following two IDs are required
+    inclusive_tauAntiMuonID = "",
+    inclusive_tauAntiElectronID = "",
+    # loose hadronic tau selection
+    loose_ptMin = 18,
+    loose_etaMax = 9999,
+    loose_dxyMax = 1000.,
+    loose_dzMax = 0.2,
+    loose_vetoLeptons = True,
+    loose_leptonVetoDR = 0.4,
+    loose_decayModeID = "decayModeFindingNewDMs", # ignored if not set or ""
+    loose_tauID = "byLooseCombinedIsolationDeltaBetaCorr3Hits",
+    loose_vetoLeptonsPOG = False, # If True, the following two IDs are required
+    loose_tauAntiMuonID = "againstMuonLoose3",
+    loose_tauAntiElectronID = "againstElectronLooseMVA5",
+    loose_tauLooseID = "decayModeFindingNewDMs"
   )
 )

VHbbAnalysis/Heppy/python/TTHtoTauTauAnalyzer.py

@@ -48,7 +48,7 @@ def process(self, event):
       self.readCollections( event.input )
 
       ##taus = list( self.handles['taus'].product() )
-      taus = event.selectedTaus
+      taus = event.inclusiveTaus
       taus_modified = []
       for idxTau in range(len(taus)):
          tau = Tau(taus[idxTau])
@@ -61,6 +61,6 @@ def process(self, event):
          #print " genMatchType = %i" % tau.genMatchType
          taus_modified.append(tau)
 
-      event.selectedTaus = taus_modified
+      event.inclusiveTaus = taus_modified
 
       return True

VHbbAnalysis/Heppy/python/TauGenJetAnalyzer.py

@@ -0,0 +1,54 @@
+
+from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
+from PhysicsTools.HeppyCore.framework.event import Event
+from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
+from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import GenParticle
+
+from PhysicsTools.Heppy.physicsutils.genutils import *
+import PhysicsTools.HeppyCore.framework.config as cfg
+
+from VHbbAnalysis.Heppy.genTauDecayMode import genTauDecayMode
+
+class TauGenJetAnalyzer( Analyzer ):
+    """Determine generator level tau decay mode."""
+
+    def __init__(self, cfg_ana, cfg_comp, looperName ):
+        super(TauGenJetAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
+
+    #---------------------------------------------
+    # DECLARATION OF HANDLES OF GEN LEVEL OBJECTS 
+    #---------------------------------------------
+
+    def declareHandles(self):
+        super(TauGenJetAnalyzer, self).declareHandles()
+
+        #mc information
+        self.mchandles['tauGenJetsSelectorAllHadrons'] = AutoHandle( 'tauGenJetsSelectorAllHadrons',
+                                                                     'vector<reco::GenJet>' )
+
+    def beginLoop(self,setup):
+        super(TauGenJetAnalyzer,self).beginLoop(setup)
+
+    def makeMCInfo(self, event):
+        event.tauGenJets = []
+        tauGenJets = list(self.mchandles['tauGenJetsSelectorAllHadrons'].product() )
+        for tauGenJet in tauGenJets:
+            tauGenJet.decayMode = genTauDecayMode(tauGenJet)[0]
+            event.tauGenJets.append(tauGenJet)
+
+    def process(self, event):
+        self.readCollections(event.input)
+
+        # if not MC, nothing to do
+        if not self.cfg_comp.isMC: 
+            return True
+
+        # do MC level analysis
+        self.makeMCInfo(event)
+
+        return True
+
+setattr(TauGenJetAnalyzer, "defaultConfig", cfg.Analyzer(
+    class_object = TauGenJetAnalyzer,
+    verbose = False,
+))

VHbbAnalysis/Heppy/python/VHbbAnalyzer.py

@@ -261,9 +261,19 @@ def classifyEvent(self,event):
 
     def fillTauIndices(self,event) :
         for j in event.cleanJetsAll :
-            j.tauIdxs = [event.selectedTaus.index(x) for x in j.taus if j.taus in event.selectedTaus]
-        for t in event.selectedTaus :
-            t.jetIdx =  event.cleanJetsAll.index(t.jet) if t.jet in event.cleanJetsAll else -1 
+            j.tauIdxs = [event.inclusiveTaus.index(x) for x in j.taus if j.taus in event.inclusiveTaus]
+        for t in event.inclusiveTaus :
+            dRmin = 1.e+3
+            t.jetIdx = -1
+            for jIdx, j in enumerate(event.cleanJetsAll) :
+                dR = None
+                if t.isPFTau():
+                   dR = deltaR(t.p4Jet().eta(),t.p4Jet().phi(),j.eta(),j.phi())
+                else:
+                   dR = deltaR(t.pfEssential().p4CorrJet_.eta(),t.pfEssential().p4CorrJet_.phi(),j.eta(),j.phi())
+                if dR < 0.3 and dR < dRmin :
+                    t.jetIdx = jIdx
+                    dRmin = dR
 
     def initOutputs (self,event) : 
         event.hJets = []

VHbbAnalysis/Heppy/python/genTauDecayMode.py

@@ -0,0 +1,60 @@
+
+#--------------------------------------------------------------------------------
+# Auxiliary function to determine generator level tau decay mode.
+#
+# The function corresponds to a translation of
+#    PhysicsTools/JetMCUtils/src/JetMCTag.cc
+# to python.
+#
+# author: Christian Veelken, NICPB Tallinn
+#--------------------------------------------------------------------------------
+
+def genTauDecayMode(tau):
+
+    numElectrons = 0
+    numMuons = 0
+    numChargedHadrons = 0
+    numNeutralHadrons = 0
+    numPhotons = 0
+
+    daughters = tau.daughterPtrVector();
+    for daughter in daughters:
+        pdg_id = abs(daughter.pdgId())
+        if pdg_id == 22:
+            numPhotons = numPhotons + 1
+        elif pdg_id == 11:
+            numElectrons = numElectrons + 1
+        elif pdg_id == 13:
+            numMuons = numMuons + 1
+        elif abs(daughter.charge()) > 0.5:
+            numChargedHadrons = numChargedHadrons + 1
+        else:
+            numNeutralHadrons = numNeutralHadrons + 1
+
+    ##print "<genTauDecayMode>: numChargedHadrons = %i, numNeutralHadrons = %i, numPhotons = %i" % (numChargedHadrons, numNeutralHadrons, numPhotons)
+
+    if numElectrons == 1:
+        return ( 0, "electron" )
+    elif numMuons == 1:
+        return ( 1, "muon" )
+    elif numChargedHadrons == 1:
+        if numNeutralHadrons != 0:
+            return ( 5, "oneProngOther" )
+        elif numPhotons == 0:
+            return ( 2, "oneProng0Pi0" )
+        elif numPhotons == 2:
+            return ( 3, "oneProng1Pi0" )
+        elif numPhotons == 4:
+            return ( 4, "oneProng2Pi0" )
+        else:
+            return ( 5, "oneProngOther" )
+    elif numChargedHadrons == 3:
+        if numNeutralHadrons != 0:
+            return ( 8, "threeProngOther" )
+        elif numPhotons == 0:
+            return ( 6, "threeProng0Pi0" )
+        elif numPhotons == 2:
+            return ( 7, "threeProng1Pi0" )
+        else:
+            return ( 8, "threeProngOther" )
+    return ( 9, "rare" )

VHbbAnalysis/Heppy/python/vhbbobj.py

@@ -32,11 +32,11 @@
     NTupleVariable("trkKink",      lambda lepton : lepton.combinedQuality().trkKink if abs(lepton.pdgId()) == 13 else 0, help="Tracker kink-finder"), 
     NTupleVariable("caloCompatibility",      lambda lepton : lepton.caloCompatibility() if abs(lepton.pdgId()) == 13 else 0, help="Calorimetric compatibility"), 
     NTupleVariable("globalTrackChi2",      lambda lepton : lepton.globalTrack().normalizedChi2() if abs(lepton.pdgId()) == 13 and lepton.globalTrack().isNonnull() else 0, help="Global track normalized chi2"), 
-    NTupleVariable("nChamberHits", lambda lepton: lepton.globalTrack().hitPattern().numberOfValidMuonHits() if abs(lepton.pdgId()) == 13 and lepton.globalTrack().isAvailable() else -1, help="Number of muon chamber hits (-1 for electrons)"),
+    NTupleVariable("nChamberHits", lambda lepton: lepton.globalTrack().hitPattern().numberOfValidMuonHits() if abs(lepton.pdgId()) == 13 and lepton.globalTrack().isNonnull() else -1, help="Number of muon chamber hits (-1 for electrons)"),
     NTupleVariable("isPFMuon", lambda lepton: lepton.isPFMuon() if abs(lepton.pdgId()) == 13 else 0, help="1 if muon passes particle flow ID"),
     NTupleVariable("isGlobalMuon", lambda lepton: lepton.isGlobalMuon() if abs(lepton.pdgId()) == 13 else 0, help="1 if muon is global muon"),
     NTupleVariable("isTrackerMuon", lambda lepton: lepton.isTrackerMuon() if abs(lepton.pdgId()) == 13 else 0, help="1 if muon is tracker muon"),
-    NTupleVariable("pixelHits", lambda lepton : lepton.innerTrack().hitPattern().numberOfValidPixelHits() if abs(lepton.pdgId()) == 13 else -1, help="Number of pixel hits (-1 for electrons)"),
+    NTupleVariable("pixelHits", lambda lepton : lepton.innerTrack().hitPattern().numberOfValidPixelHits() if abs(lepton.pdgId()) == 13 and lepton.innerTrack().isNonnull() else -1, help="Number of pixel hits (-1 for electrons)"),
     # Extra tracker-related id variables
     NTupleVariable("trackerLayers", lambda x : (x.track() if abs(x.pdgId())==13 else x.gsfTrack()).hitPattern().trackerLayersWithMeasurement(), int, help="Tracker Layers"),
     NTupleVariable("pixelLayers", lambda x : (x.track() if abs(x.pdgId())==13 else x.gsfTrack()).hitPattern().pixelLayersWithMeasurement(), int, help="Pixel Layers"),
@@ -183,6 +183,9 @@
 #    NTupleVariable("score",  lambda x : x.mass()), # to be added for 74X
 ])
 
+genTauJetType = NTupleObjectType("genTauJet", baseObjectTypes = [ genParticleType ], variables = [
+    NTupleVariable("decayMode", lambda x : x.decayMode, int, mcOnly=True, help="Generator level tau decay mode"),
+])
 
 def ptRel(p4,axis):
     a=ROOT.TVector3(axis.Vect().X(),axis.Vect().Y(),axis.Vect().Z())

VHbbAnalysis/Heppy/test/combined_cmssw.py

@@ -139,3 +139,12 @@
 process.OUT.outputCommands.append("keep *_ca15PFTrimmedJetsCHS_*_EX")
 process.OUT.outputCommands.append("keep *_looseMultiRHTT_*_EX")
 
+########################################
+# Generator level hadronic tau decays
+########################################
+
+process.load("PhysicsTools.JetMCAlgos.TauGenJets_cfi")
+process.tauGenJets.GenParticles = cms.InputTag('prunedGenParticles')
+process.load("PhysicsTools.JetMCAlgos.TauGenJetsDecayModeSelectorAllHadrons_cfi")
+
+process.OUT.outputCommands.append("keep *_tauGenJetsSelectorAllHadrons_*_EX")

VHbbAnalysis/Heppy/test/vhbb-combined.py

@@ -30,7 +30,6 @@
                                                              10,
                                                              help="MultiR HEPTopTagger Candidates")
 
-
 # Add b-Tagging Information
 
 btagana=cfg.Analyzer(
@@ -39,6 +38,16 @@
 )
 sequence.insert(sequence.index(VHbb),btagana)
 
+# Add Information on generator level hadronic tau decays
+
+from VHbbAnalysis.Heppy.TauGenJetAnalyzer import TauGenJetAnalyzer
+TauGenJet = cfg.Analyzer(
+    verbose = False,
+    class_object = TauGenJetAnalyzer,
+)
+sequence.insert(sequence.index(VHbb),TauGenJet)
+
+treeProducer.collections["tauGenJets"] = NTupleCollection("GenHadTaus", genTauJetType, 15, help="Generator level hadronic tau decays")
 
 # Run Everything