Merge pull request #17 from arizzi/vhbbHeppy722patch2-taustuff

Vhbb heppy722patch2 taustuff

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

@@ -165,7 +165,16 @@ def process(self, event):
                 lep.jet = lep
             else:
                 lep.jet = jet
+        ## Associate jets to taus 
+        taus = getattr(event,'selectedTaus',[])
+        jtaupairs = matchObjectCollection( taus, allJets, self.jetLepDR**2)
 
+        for jet in allJets:
+            jet.taus = [l for l in jtaupairs if jtaupairs[l] == jet ]
+        for tau in taus:
+            tau.jet = jtaupairs[tau]
+
+        #MC stuff
         if self.cfg_comp.isMC:
             event.deltaMetFromJetSmearing = [0, 0]
             for j in event.cleanJetsAll:

VHbbAnalysis/Heppy/python/TTHtoTauTauAnalyzer.py

@@ -0,0 +1,66 @@
+from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
+from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
+from PhysicsTools.Heppy.physicsobjects.Tau import Tau
+
+from PhysicsTools.HeppyCore.utils.deltar import deltaR, matchObjectCollection3
+
+import PhysicsTools.HeppyCore.framework.config as cfg
+
+class TTHtoTauTauAnalyzer( Analyzer ):
+   '''Analyze ttH, H -> tautau events'''
+
+   def declareHandles(self):
+      super(TTHtoTauTauAnalyzer, self).declareHandles()
+
+      ##self.handles['taus'] = AutoHandle( ('slimmedTaus',''), 'std::vector<pat::Tau>' )
+
+      #mc information      
+      self.mchandles['genParticles'] = AutoHandle( 'prunedGenParticles',
+                                                   'std::vector<reco::GenParticle>' )
+
+   def addTau_genMatchType(self, event, tau):      
+      '''Determine if given tau matched gen level hadronic tau decay or is due to a misidentified jet, electron or muon
+             tau.genMatchType = 0 for matched to gen level hadronic tau decay
+                              = 1 for matched to gen level jet
+                              = 2 for matched to gen level electron
+                              = 3 for matched to gen level muon
+      '''
+
+      genParticles = list(self.mchandles['genParticles'].product() )
+
+      genMatchType = 1 # assume hadronic tau to be due to misidentified jet per default
+      if tau.genJet():
+         genMatchType = 0
+      if genMatchType == 1:
+         match = matchObjectCollection3([ tau ], genParticles, deltaRMax = 0.4, filter = lambda x,y : True if (y.pt() > 0.5*x.pt() and abs(y.pdgId()) == 11) else False)
+         if match[tau]:
+            genMatchType = 2
+      if genMatchType == 1:
+         match = matchObjectCollection3([ tau ], genParticles, deltaRMax = 0.4, filter = lambda x,y : True if (y.pt() > 0.5*x.pt() and abs(y.pdgId()) == 13) else False)
+         if match[tau]:
+            genMatchType = 3
+
+      return genMatchType
+
+   def process(self, event):
+      #print "<TTHtoTauTauAnalyzer::process>:"
+
+      self.readCollections( event.input )
+
+      ##taus = list( self.handles['taus'].product() )
+      taus = event.selectedTaus
+      taus_modified = []
+      for idxTau in range(len(taus)):
+         tau = Tau(taus[idxTau])
+         #print "processing tau #%i: Pt = %1.2f, eta = %1.2f, phi = %1.2f" % (idxTau, tau.pt(), tau.eta(), tau.phi())
+         # if not MC, nothing to do
+         if self.cfg_comp.isMC:
+            tau.genMatchType = self.addTau_genMatchType(event, tau)
+         else:
+            tau.genMatchType = -1
+         #print " genMatchType = %i" % tau.genMatchType
+         taus_modified.append(tau)
+
+      event.selectedTaus = taus_modified
+
+      return True

VHbbAnalysis/Heppy/python/TTHtoTauTauGeneratorAnalyzer.py

@@ -0,0 +1,148 @@
+
+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
+
+class TTHtoTauTauGeneratorAnalyzer( Analyzer ):
+    """Do generator-level analysis of ttH, H -> tautau decay:
+
+       Creates in the event:
+         event.genTTHtoTauTauDecayMode = 0 for 2b_2taul_2wl
+                                         1 for 2b_2taulh_2wl
+                                         2 for 2b_2tauh_2wl
+                                         3 for 2b_2taul_2wlj
+                                         4 for 2b_2taulh_2wlj
+                                         5 for 2b_2tauh_2wlj
+                                         6 for 2b_2taul_2wj
+                                         7 for 2b_2taulh_2wj
+                                         8 for 2b_2tauh_2wj
+                                        -1 other
+    """
+    def __init__(self, cfg_ana, cfg_comp, looperName ):
+        super(TTHtoTauTauGeneratorAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
+
+    #---------------------------------------------
+    # DECLARATION OF HANDLES OF GEN LEVEL OBJECTS 
+    #---------------------------------------------
+
+    def declareHandles(self):
+        super(TTHtoTauTauGeneratorAnalyzer, self).declareHandles()
+
+        #mc information
+        self.mchandles['genParticles'] = AutoHandle( 'prunedGenParticles',
+                                                     'std::vector<reco::GenParticle>' )
+
+    def beginLoop(self,setup):
+        super(TTHtoTauTauGeneratorAnalyzer,self).beginLoop(setup)
+
+    def findFirstDaughterGivenPdgId(self, givenMother, givenPdgIds):
+        """Get the first gen level daughter particle with given pdgId"""
+
+        mothers = []
+        if givenMother:
+            mothers.append(givenMother)
+        while len(mothers) > 0:
+            daughters = []
+            for mother in mothers:
+                for idxDaughter in range(mother.numberOfDaughters()):
+                    daughter = GenParticle(mother.daughter(idxDaughter))
+                    daughters.append(daughter)
+            for daughter in daughters:
+                if daughter.pdgId() in givenPdgIds:
+                    return daughter
+            mothers = daughters
+
+        # no daughter particle with given pdgId found
+        return None
+
+    def fillGenTTHtoTauTauDecayMode(self, event, genParticles):
+        """Determine gen level ttH, H -> tautau decay mode"""
+
+        t = None
+        tbar = None
+        H = None
+
+        for genParticle in genParticles:
+            pdgId = genParticle.pdgId()
+            if pdgId == +6 and not t:
+                t = genParticle
+            if pdgId == -6 and not tbar:
+                tbar = genParticle
+            if pdgId in [ 25, 35, 36 ] and not H:
+                H = genParticle
+
+        numElectrons_or_Muons_fromH = 0
+        numHadTaus_fromH = 0
+        if self.findFirstDaughterGivenPdgId(H, [ -11, -13 ]):
+            numElectrons_or_Muons_fromH = numElectrons_or_Muons_fromH + 1
+        elif self.findFirstDaughterGivenPdgId(H, [ -15 ]):
+            numHadTaus_fromH = numHadTaus_fromH + 1
+        if self.findFirstDaughterGivenPdgId(H, [ +11, +13 ]):
+            numElectrons_or_Muons_fromH = numElectrons_or_Muons_fromH + 1
+        elif self.findFirstDaughterGivenPdgId(H, [ +15 ]):
+            numHadTaus_fromH = numHadTaus_fromH + 1
+
+        numBs = 0
+        if self.findFirstDaughterGivenPdgId(t, [ +5 ]):
+            numBs = numBs + 1
+        if self.findFirstDaughterGivenPdgId(tbar, [ -5 ]):
+            numBs = numBs + 1
+
+        Wplus = self.findFirstDaughterGivenPdgId(t, [ +24 ])
+        Wminus = self.findFirstDaughterGivenPdgId(tbar, [ -24 ])
+
+        numElectrons_or_Muons_fromW = 0
+        numHadTaus_fromW = 0
+        if self.findFirstDaughterGivenPdgId(Wplus, [ -11, -13 ]):
+            numElectrons_or_Muons_fromW = numElectrons_or_Muons_fromW + 1
+        elif self.findFirstDaughterGivenPdgId(Wplus, [ -15 ]):
+            numHadTaus_fromW = numHadTaus_fromW + 1
+        if self.findFirstDaughterGivenPdgId(Wminus, [ +11, +13 ]):
+            numElectrons_or_Muons_fromW = numElectrons_or_Muons_fromW + 1
+        elif self.findFirstDaughterGivenPdgId(Wminus, [ +15 ]):
+            numHadTaus_fromW = numHadTaus_fromW + 1
+
+        event.genTTHtoTauTauDecayMode = -1
+        if numElectrons_or_Muons_fromH == 2 and (numElectrons_or_Muons_fromW + numHadTaus_fromW) == 2:
+            event.genTTHtoTauTauDecayMode = 0
+        elif numElectrons_or_Muons_fromH == 1 and numHadTaus_fromH == 1 and (numElectrons_or_Muons_fromW + numHadTaus_fromW) == 2:
+            event.genTTHtoTauTauDecayMode = 1
+        elif numHadTaus_fromH == 2 and numElectrons_or_Muons_fromW == 2:
+            event.genTTHtoTauTauDecayMode = 2
+        elif numElectrons_or_Muons_fromH == 2 and (numElectrons_or_Muons_fromW + numHadTaus_fromW) == 1:
+            event.genTTHtoTauTauDecayMode = 3
+        elif numElectrons_or_Muons_fromH == 1 and numHadTaus_fromH == 1 and (numElectrons_or_Muons_fromW + numHadTaus_fromW) == 1:
+            event.genTTHtoTauTauDecayMode = 4
+        elif numHadTaus_fromH == 2 and (numElectrons_or_Muons_fromW + numHadTaus_fromW) == 1:
+            event.genTTHtoTauTauDecayMode = 5
+        elif numElectrons_or_Muons_fromH == 2 and (numElectrons_or_Muons_fromW + numHadTaus_fromW) == 0:
+            event.genTTHtoTauTauDecayMode = 6
+        elif numElectrons_or_Muons_fromH == 1 and numHadTaus_fromH == 1 and numElectrons_or_Muons_fromW == 0:
+            event.genTTHtoTauTauDecayMode = 7
+        elif numHadTaus_fromH == 2 and numElectrons_or_Muons_fromW == 0:
+            event.genTTHtoTauTauDecayMode = 8
+
+    def makeMCInfo(self, event):
+        genParticles = list(self.mchandles['genParticles'].product() )
+        self.fillGenTTHtoTauTauDecayMode(event, genParticles)
+
+    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(TTHtoTauTauGeneratorAnalyzer, "defaultConfig", cfg.Analyzer(
+    class_object = TTHtoTauTauGeneratorAnalyzer,
+    verbose = False,
+))

VHbbAnalysis/Heppy/python/VHGeneratorAnalyzer.py

@@ -294,7 +294,7 @@ def process(self, event):
         return True
 
 setattr(GeneratorAnalyzer,"defaultConfig",cfg.Analyzer(
-    class_object=GeneratorAnalyzer,
+    class_object = GeneratorAnalyzer,
     filterHiggsDecays = False, 
     verbose = False,
     PDFWeights = []

VHbbAnalysis/Heppy/python/VHbbAnalyzer.py

@@ -228,6 +228,13 @@ def classifyEvent(self,event):
 	event.aLeptons = [x for x in event.inclusiveLeptons if x not in event.vLeptons]
 
 	return True
+
+    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 
+
     def initOutputs (self,event) : 
         event.hJets = []
         event.aJets = []
@@ -274,7 +281,7 @@ def process(self, event):
 	self.doHiggsHighCSV(event)
 	self.doHiggsHighPt(event)
         self.doHiggs3cj(event)
-
+        self.fillTauIndices(event)
 
 
     #    event.jee = list(self.handles['jee'].product())

VHbbAnalysis/Heppy/python/vhbbobj.py

@@ -6,7 +6,6 @@
 from PhysicsTools.Heppy.analyzers.objects.autophobj import *
 import copy
 
-
 leptonTypeVHbb = NTupleObjectType("leptonTypeVHbb", baseObjectTypes = [ leptonType ], variables = [
     # Loose id 
     NTupleVariable("looseIdSusy", lambda x : x.looseIdSusy if hasattr(x, 'looseIdSusy') else -1, int, help="Loose ID for Susy ntuples (always true on selected leptons)"),
@@ -47,11 +46,21 @@
 #    NTupleVariable("mcMatchTau",  lambda x : x.mcMatchTau, int, mcOnly=True, help="True if the leptons comes from a tau"),
 ])
 
+##------------------------------------------  
+## TAU
+##------------------------------------------  
+
+tauTypeVHbb = NTupleObjectType("tauTypeVHbb", baseObjectTypes = [ tauType ], variables = [
+    NTupleVariable("idxJetMatch", lambda x : x.jetIdx, int, help="index of the matching jet"),
+    NTupleVariable("genMatchType", lambda x : x.genMatchType, int,mcOnly=True, help="..FILLME PLEASE..")
+])
+
 ##------------------------------------------  
 ## JET
 ##------------------------------------------  
 
 jetTypeVHbb = NTupleObjectType("jet",  baseObjectTypes = [ jetType ], variables = [
+    NTupleVariable("idxFirstTauMatch", lambda x : x.tauIdxs[0] if len(x.tauIdxs) > 0 else -1, int,help='index of the first matching tau'),
     NTupleVariable("hadronFlavour", lambda x : x.mcFlavour, int,     mcOnly=True, help="match to heavy hadrons"),
     NTupleVariable("btagBDT", lambda x : getattr(x,"btagBDT",-99), help="btag"),
     NTupleVariable("btagProb", lambda x : x.btag('jetProbabilityBJetTags') , help="btag"),

VHbbAnalysis/Heppy/test/vhbb-btag.py

@@ -10,7 +10,7 @@
 config.preprocessor=preprocessor
 if __name__ == '__main__':
     from PhysicsTools.HeppyCore.framework.looper import Looper 
-    looper = Looper( 'Loop', config, nPrint = 20, nEvents = 300)
+    looper = Looper( 'Loop', config, nPrint = 20, nEvents = 3000)
     import time
     import cProfile
     p = cProfile.Profile(time.clock)