Merge pull request cms-sw#285 from mmarionncern/MM-CMG022515-v2

Addition of Jet variables/gen event histograms to the susy multilepton analysis/heppy updates (jet-tau, LHE analyzer and PDF weights)

CMGTools/TTHAnalysis/cfg/run_susyMultilepton_cfg.py

@@ -34,6 +34,13 @@
     minJets25 = 0,
     )
 
+##==== tau jet analyzer, to be called (for the moment) once bjetsMedium are produced
+from CMGTools.TTHAnalysis.analyzers.ttHJetTauAnalyzer import ttHJetTauAnalyzer
+ttHJetTauAna = cfg.Analyzer(
+    ttHJetTauAnalyzer, name="ttHJetTauAnalyzer",
+    )
+
+
 ## Insert the SV analyzer in the sequence
 susyCoreSequence.insert(susyCoreSequence.index(ttHCoreEventAna), 
                         ttHFatJetAna)
@@ -96,6 +103,12 @@
      collections = susyMultilepton_collections,
 )
 
+## histo counter
+TFileServiceMode=False
+if TFileServiceMode:
+    susyCoreSequence.insert(susyCoreSequence.index(skimAnalyzer),
+                            susyCounter)
+
 #-------- SAMPLES AND TRIGGERS -----------
 
 #-------- SEQUENCE
@@ -117,14 +130,16 @@
 #selectedComponents = [T5ttttDeg_mGo1000_mStop300_mCh285_mChi280, T5ttttDeg_mGo1000_mStop300_mCh285_mChi280_dil, TTWJets, TTZJets, WZJetsTo3LNu]
 
 sequence = cfg.Sequence(susyCoreSequence+[
-    ttHEventAna,
-    treeProducer,
+        ttHJetTauAna,
+        ttHEventAna,
+        treeProducer,
     ])
 
 # -- fine splitting, for some private MC samples with a single file
 #for comp in selectedComponents:
 #    comp.splitFactor = 1
 #    comp.fineSplitFactor = 40
+
 
 from PhysicsTools.HeppyCore.framework.heppy import getHeppyOption
 test = getHeppyOption('test')
@@ -176,8 +191,18 @@
     comp.fineSplitFactor = 10
     selectedComponents = [ comp ]
 
-
-
+## output histogram
+outputService=[]
+if TFileServiceMode:
+    from PhysicsTools.HeppyCore.framework.services.tfile import TFileService
+    output_service = cfg.Service(
+        TFileService,
+        'outputfile',
+        name="outputfile",
+        fname='treeProducerSusyMultilepton/tree.root',
+        option='recreate'
+        )    
+    outputService.append(output_service)
 
 # the following is declared in case this cfg is used in input to the heppy.py script
 from PhysicsTools.HeppyCore.framework.eventsfwlite import Events
@@ -187,5 +212,5 @@
     event_class = Events 
 config = cfg.Config( components = selectedComponents,
                      sequence = sequence,
-                     services = [],  
+                     services = outputService,  
                      events_class = event_class)

CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyCore.py

@@ -4,16 +4,31 @@
             NTupleVariable("rho",  lambda ev: ev.rho, float, help="kt6PFJets rho"),
             NTupleVariable("nVert",  lambda ev: len(ev.goodVertices), int, help="Number of good vertices"),
             NTupleVariable("nJet25", lambda ev: len(ev.cleanJets), int, help="Number of jets with pt > 25"),
+            NTupleVariable("nJet25NoTau", lambda ev: sum([ (j.pt() > 25 and not j.taus ) for j in ev.cleanJets]), int, help="Number of jets with pt > 25, not matched with taus"),
             NTupleVariable("nBJetLoose25", lambda ev: len(ev.bjetsLoose), int, help="Number of jets with pt > 25 passing CSV loose"),
             NTupleVariable("nBJetMedium25", lambda ev: len(ev.bjetsMedium), int, help="Number of jets with pt > 25 passing CSV medium"),
             NTupleVariable("nBJetTight25", lambda ev: sum([j.btagWP("CSVv2IVFT") for j in ev.bjetsMedium]), int, help="Number of jets with pt > 25 passing CSV tight"),
             NTupleVariable("nBJetCMVALoose25", lambda ev: sum([j.btagWP("CMVAL") for j in ev.cleanJets]), int, help="Number of jets with pt > 25 passing CMVA Loose"),
             NTupleVariable("nBJetCMVAMedium25", lambda ev: sum([j.btagWP("CMVAM") for j in ev.cleanJets]), int, help="Number of jets with pt > 25 passing CMVA Medium"),
             NTupleVariable("nBJetCMVATight25", lambda ev: sum([j.btagWP("CMVAT") for j in ev.cleanJets]), int, help="Number of jets with pt > 25 passing CMVA Tight"),
+            NTupleVariable("nBJetLoose25NoTau", lambda ev: sum([(j.btagWP("CSVv2IVFL") and not j.taus) for j in ev.cleanJets]), int, help="Number of jets with pt > 25 passing CSV loose, not matched with taus"),
+            NTupleVariable("nBJetMedium25NoTau", lambda ev: sum([(j.btagWP("CSVv2IVFM") and not j.taus) for j in ev.bjetsMedium]), int, help="Number of jets with pt > 25 passing CSV medium, not matched with taus"),
+
+            NTupleVariable("nJet30", lambda ev: sum([j.pt() > 30 for j in ev.cleanJets]), int, help="Number of jets with pt > 30, |eta|<2.4"),
+            NTupleVariable("nJet30a", lambda ev: sum([j.pt() > 30 for j in ev.cleanJetsAll]), int, help="Number of jets with pt > 30, |eta|<4.7"),
+            NTupleVariable("nJet30NoTau", lambda ev: sum([ (j.pt() > 30 and not j.taus ) for j in ev.cleanJets]), int, help="Number of jets with pt > 30, not matched with taus"),
+            NTupleVariable("nBJetLoose30", lambda ev: sum([j.btagWP("CSVv2IVFL") for j in ev.cleanJets if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV loose"),
+            NTupleVariable("nBJetLoose30NoTau", lambda ev: sum([(j.btagWP("CSVv2IVFL") and not j.taus) for j in ev.cleanJets if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV loose, not matched with taus"),
+            NTupleVariable("nBJetMedium30", lambda ev: sum([j.btagWP("CSVv2IVFM") for j in ev.bjetsMedium if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV medium"),
+            NTupleVariable("nBJetMedium30NoTau", lambda ev: sum([(j.btagWP("CSVv2IVFM") and not j.taus) for j in ev.bjetsMedium if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV medium, not matched with taus"),
+            NTupleVariable("nBJetTight30", lambda ev: sum([j.btagWP("CSVv2IVFT") for j in ev.bjetsMedium if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV tight"),
             NTupleVariable("nJet40", lambda ev: sum([j.pt() > 40 for j in ev.cleanJets]), int, help="Number of jets with pt > 40, |eta|<2.4"),
             NTupleVariable("nJet40a", lambda ev: sum([j.pt() > 40 for j in ev.cleanJetsAll]), int, help="Number of jets with pt > 40, |eta|<4.7"),
+            NTupleVariable("nJet40NoTau", lambda ev: sum([ (j.pt() > 40 and not j.taus ) for j in ev.cleanJets]), int, help="Number of jets with pt > 40, not matched with taus"),
             NTupleVariable("nBJetLoose40", lambda ev: sum([j.btagWP("CSVv2IVFL") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV loose"),
+            NTupleVariable("nBJetLoose40NoTau", lambda ev: sum([(j.btagWP("CSVv2IVFL") and not j.taus) for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV loose, not matched with taus"),
             NTupleVariable("nBJetMedium40", lambda ev: sum([j.btagWP("CSVv2IVFM") for j in ev.bjetsMedium if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV medium"),
+            NTupleVariable("nBJetMedium40NoTau", lambda ev: sum([(j.btagWP("CSVv2IVFM") and not j.taus) for j in ev.bjetsMedium if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV medium, not matched with taus"),
             NTupleVariable("nBJetTight40", lambda ev: sum([j.btagWP("CSVv2IVFT") for j in ev.bjetsMedium if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV tight"),
             NTupleVariable("nBJetCMVALoose40", lambda ev: sum([j.btagWP("CMVAL") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CMVA Loose"),
             NTupleVariable("nBJetCMVAMedium40", lambda ev: sum([j.btagWP("CMVAM") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CMVA Medium"),

CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyGenLevelOnly.py

@@ -18,6 +18,8 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
             NTupleVariable("nJet40a", lambda ev: sum([j.pt() > 40 for j in ev.cleanJetsAll]), int, help="Number of jets with pt > 40, |eta|<4.7"),
             NTupleVariable("nBJetLoose25", lambda ev: sum([j.btagWP("CSVv2IVFL") for j in ev.cleanJets if j.pt() > 25]), int, help="Number of jets with pt > 25 passing CSV loose"),
             NTupleVariable("nBJetMedium25", lambda ev: sum([j.btagWP("CSVv2IVFM") for j in ev.cleanJets if j.pt() > 25]), int, help="Number of jets with pt > 25 passing CSV medium"),
+            NTupleVariable("nBJetLoose30", lambda ev: sum([j.btagWP("CSVv2IVFL") for j in ev.cleanJets if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV loose"),
+            NTupleVariable("nBJetMedium30", lambda ev: sum([j.btagWP("CSVv2IVFM") for j in ev.cleanJets if j.pt() > 30]), int, help="Number of jets with pt > 30 passing CSV medium"),
             NTupleVariable("nBJetLoose40", lambda ev: sum([j.btagWP("CSVv2IVFL") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV loose"),
             NTupleVariable("nBJetMedium40", lambda ev: sum([j.btagWP("CSVv2IVFM") for j in ev.cleanJets if j.pt() > 40]), int, help="Number of jets with pt > 40 passing CSV medium"),
             ##--------------------------------------------------

CMGTools/TTHAnalysis/python/analyzers/treeProducerSusyMultilepton.py

@@ -14,6 +14,7 @@
             NTupleVariable("nSoftBLoose25",  lambda ev: sum([(sv.mva>0.3 and (sv.jet == None or sv.jet.pt() < 25)) for sv in ev.ivf]), int, help="SV from ivf with loose sv mva not matched to a jet with pt > 25 GeV"),
             NTupleVariable("nSoftBMedium25", lambda ev: sum([(sv.mva>0.7 and (sv.jet == None or sv.jet.pt() < 25)) for sv in ev.ivf]), int, help="SV from ivf with medium sv mva not matched to a jet with pt > 25 GeV"),
             NTupleVariable("nSoftBTight25",  lambda ev: sum([(sv.mva>0.9 and (sv.jet == None or sv.jet.pt() < 25)) for sv in ev.ivf]), int, help="SV from ivf with tight sv mva not matched to a jet with pt > 25 GeV"),
+
             ##--------------------------------------------------            
             NTupleVariable("mZ1", lambda ev : ev.bestZ1[0], help="Best m(ll) SF/OS"),
             NTupleVariable("mZ1SFSS", lambda ev : ev.bestZ1sfss[0], help="Best m(ll) SF/SS"),
@@ -116,7 +117,7 @@
             "selectedLeptons" : NTupleCollection("LepGood",  leptonTypeSusyExtra, 8, help="Leptons after the preselection"),
             "otherLeptons"    : NTupleCollection("LepOther", leptonTypeSusyExtra, 8, help="Leptons after the preselection"),
             ##------------------------------------------------
-            "cleanJets"       : NTupleCollection("Jet",     jetTypeSusy, 15, help="Cental jets after full selection and cleaning, sorted by pt"),
+            "cleanJets"       : NTupleCollection("Jet",     jetTypeSusyExtra, 15, help="Cental jets after full selection and cleaning, sorted by pt"),
             "cleanJetsFwd"    : NTupleCollection("JetFwd",  jetTypeSusy,  6, help="Forward jets after full selection and cleaning, sorted by pt"),            
             "fatJets"         : NTupleCollection("FatJet",  fatJetType,  15, help="AK8 jets, sorted by pt"),
             ##------------------------------------------------
@@ -126,4 +127,8 @@
             "ivf"       : NTupleCollection("SV",     svType, 20, help="SVs from IVF"),
             "genBHadrons"  : NTupleCollection("GenBHad", heavyFlavourHadronType, 20, mcOnly=True, help="Gen-level B hadrons"),
             "genDHadrons"  : NTupleCollection("GenDHad", heavyFlavourHadronType, 20, mcOnly=True, help="Gen-level D hadrons"),
+            ##------------------------------------------------
+            "jetsNonTauIdx" : NTupleCollection("JetNoTauIdx",objectInt, 10, help="index of jets not associated to taus"),
+
+
 })

CMGTools/TTHAnalysis/python/analyzers/ttHJetTauAnalyzer.py

@@ -0,0 +1,28 @@
+from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
+#from CMGTools.RootTools.fwlite.AutoHandle import AutoHandle
+
+
+class ttHJetTauAnalyzer( Analyzer ):
+    def __init__(self, cfg_ana, cfg_comp, looperName ):
+        super(ttHJetTauAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
+
+    def declareHandles(self):
+        super(ttHJetTauAnalyzer, self).declareHandles()
+
+    def beginLoop(self, setup):
+        super(ttHJetTauAnalyzer,self).beginLoop(setup)
+
+    def findNonTauJets(self, jets):
+        iJetNoTau =  []
+        for ij in xrange(len(jets)):
+            if not jets[ij].taus:
+                iJetNoTau.append(ij)
+
+        return iJetNoTau
+
+    def process(self, event):
+        self.readCollections(event.input)
+
+        event.jetsNonTauIdx = self.findNonTauJets(event.cleanJets)
+
+        return True

PhysicsTools/Heppy/python/analyzers/core/AutoFillTreeProducer.py

@@ -170,4 +170,33 @@ def fillTree(self, event, resetFirst=True):
                 c.fillBranchesVector(self.tree, getattr(event, cn), isMC)
 
         self.tree.tree.Fill()      
+
+    def getPythonWrapper(self):
+        """
+        This function produces a string that contains a Python wrapper for the event.
+        The wrapper is automatically generated based on the collections and allows the full
+        event contents to be accessed from subsequent Analyzers using e.g.
+        
+        leps = event.selLeptons #is of type selLeptons
+        pt0 = leps[0].pt
+
+        One just needs to add the EventAnalyzer to the sequence.
+        """
+
+        isMC = self.cfg_comp.isMC 
+
+        classes = ""
+        anclass = ""
+        anclass += "from PhysicsTools.HeppyCore.framework.analyzer import Analyzer\n"
+        anclass += "class EventAnalyzer(Analyzer):\n"
+        anclass += "    def __init__(self, cfg_ana, cfg_comp, looperName):\n"
+        anclass += "        super(EventAnalyzer, self).__init__(cfg_ana, cfg_comp, looperName)\n"
+
+        anclass += "    def process(self, event):\n"
+
+        for cname, coll in self.collections.items():
+            classes += coll.get_py_wrapper_class(isMC)
+            anclass += "        event.{0} = {0}.make_array(event)\n".format(coll.name)
+
+        return classes + "\n" + anclass
 

PhysicsTools/Heppy/python/analyzers/core/TreeAnalyzerNumpy.py

@@ -27,8 +27,8 @@ def beginLoop(self, setup) :
             print 'Compression', isCompressed
             self.file = TFile( fileName, 'recreate', '', isCompressed )
         self.file.cd()
-	if self.file.Get(self.treename) :
-		raise RuntimeError, "You are booking two Trees with the same name in the same file"
+        if self.file.Get(self.treename) :
+            raise RuntimeError, "You are booking two Trees with the same name in the same file"
         self.tree = Tree(self.treename, self.name)
         self.tree.setDefaultFloatType(getattr(self.cfg_ana, 'defaultFloatType','D')); # or 'F'
         self.declareVariables(setup)

PhysicsTools/Heppy/python/analyzers/core/autovars.py

@@ -146,4 +146,32 @@ def fillBranchesVector(self,treeNumpy,collection,isMC):
     def __repr__(self):
         return "<NTupleCollection[%s]>" % self.name
 
+    def get_cpp_declaration(self, isMC):
+        s = []
+        for v in self.objectType.allVars(isMC):
+            s += ["{0} {1}__{2}[{3}];".format(v.type.__name__, self.name, v.name, self.maxlen)]
+        return "\n".join(s)
+
+    def get_cpp_wrapper_class(self, isMC):
+        s = "class %s {\n" % self.name
+        s += "public:\n"
+        for v in self.objectType.allVars(isMC):
+            s += "    {0} {1};\n".format(v.type.__name__, v.name)
+        s += "};\n"
+        return s
+
+    def get_py_wrapper_class(self, isMC):
+        s = "class %s:\n" % self.name
+        s += "    def __init__(self, tree, n):\n"
+        for v in self.objectType.allVars(isMC):
+            if len(v.name)>0:
+                s += "        self.{0} = tree.{1}_{2}[n];\n".format(v.name, self.name, v.name)
+            else:
+                s += "        self.{0} = tree.{0}[n];\n".format(self.name)
+
+        s += "    @staticmethod\n"
+        s += "    def make_array(event):\n"
+        s += "        return [{0}(event.input, i) for i in range(event.input.n{0})]\n".format(self.name)
+        return s
+
 

PhysicsTools/Heppy/python/analyzers/gen/GeneratorAnalyzer.py

@@ -40,6 +40,7 @@ class GeneratorAnalyzer( Analyzer ):
             event.genwzquarks   = [] # quarks from W,Z decays
             event.genbquarksFromTop = []
             event.genbquarksFromH   = []
+            event.genlepsFromTop = [] #mu/ele that have a t->W chain as ancestor, also contained in event.genleps
        event.genwzquarks and event.genbquarks, might have overlaps 
        event.genbquarksFromTop and event.genbquarksFromH are all contained in event.genbquarks
        
@@ -183,6 +184,7 @@ def makeMCInfo(self, event):
             event.genwzquarks    = []
             event.genbquarksFromTop = []
             event.genbquarksFromH   = []
+            event.genlepsFromTop = []
             for p in event.generatorSummary:
                 id = abs(p.pdgId())
                 if id == 25: 
@@ -192,10 +194,23 @@ def makeMCInfo(self, event):
                 elif id in {12,14,16}:
                     event.gennus.append(p)
                 elif id in {11,13}:
+                    #taus to separate vector
                     if abs(p.motherId) == 15:
                         event.gentauleps.append(p)
+                    #all muons and electrons
                     else:
                         event.genleps.append(p)
+                        momids = [(m, abs(m.pdgId())) for m in realGenMothers(p)]
+
+                        #have a look at the lepton mothers
+                        for mom, momid in momids:
+                            #lepton from W
+                            if momid == 24:
+                                wmomids = [abs(m.pdgId()) for m in realGenMothers(mom)]
+                                #W from t
+                                if 6 in wmomids:
+                                    #save mu,e from t->W->mu/e
+                                    event.genlepsFromTop.append(p)
                 elif id == 15:
                     if self.allGenTaus or not any([abs(d.pdgId()) in {11,13} for d in realGenDaughters(p)]):
                         event.gentaus.append(p)

PhysicsTools/Heppy/python/analyzers/gen/LHEAnalyzer.py

@@ -24,6 +24,10 @@ def process(self, event):
             return True
         event.lheHT=0
         event.lheNj=0
+        event.lheNb=0
+        event.lheNc=0
+        event.lheNl=0
+        event.lheNg=0
         event.lheV_pt = 0
         try:
           event.input.getByLabel( 'externalLHEProducer',self.lheh)
@@ -46,6 +50,14 @@ def process(self, event):
           if status == 1 and ( ( idabs == 21 ) or (idabs > 0 and idabs < 7) ) : # gluons and quarks
               event.lheHT += sqrt( pup[i][0]**2 + pup[i][1]**2 ) # first entry is px, second py
               event.lheNj +=1
+              if idabs==5:
+                event.lheNb += 1
+              if idabs==4:
+                event.lheNc += 1
+              if idabs in [1,2,3]:
+                event.lheNl += 1
+              if idabs==21:
+                event.lheNg += 1
           if idabs in [12,14,16] :  
               if id > 0 :
                 nu = i