Merge pull request cms-sw#572 from bachtis/VV_7_4_12

VV: Improvements and Statistics tools

CMGTools/VVResonances/cfg/runVV_cfg.py

@@ -38,36 +38,36 @@
     "HT800":triggers_HT800,
     "HT900":triggers_HT900,
     "JJ":triggers_dijet_fat,  
-    "MET90":triggers_Jet80MET90,
-    "MET120":triggers_Jet80MET120
+    "MET90":triggers_met90_mht90+triggers_metNoMu90_mhtNoMu90,
+    "MET120":triggers_metNoMu120_mhtNoMu120
 }
 
 
 #-------- HOW TO RUN
-test = 0
+test = 1
 if test==1:
     # test a single component, using a single thread.
-    selectedComponents = [ZJetsToNuNu_HT100to200]
+    selectedComponents = [VBF_RadionToZZ_narrow_4500]
     for c in selectedComponents:
+        c.files = c.files[:1]
         c.splitFactor = 1
-
 elif test==2:    
     # test all components (1 thread per component).
-    selectedComponents = [dataSamples[0]]
+    selectedComponents = [BulkGravToWW_narrow_2500]
     for comp in selectedComponents:
         comp.splitFactor = 1
 #        comp.files = comp.files[:1]
+elif test==3:    
+    # test all components (1 thread per component).
+    selectedComponents = [DYJetsToLL_M50_HT600toInf]
+    for comp in selectedComponents:
+        comp.splitFactor = 1
 
-elif test==3:
-    # test a single component, using a single thread.
-    selectedComponents = [signalSamples[0]]
-    for c in selectedComponents:
-        c.files = c.files[:1]
-        c.splitFactor = 1
-
-
-
-
+elif test==4:    
+    # test all components (1 thread per component).
+    selectedComponents = [RSGravToWWToLNQQ_kMpl01_4500]
+    for comp in selectedComponents:
+        comp.splitFactor = 20
 
 ## output histogram
 outputService=[]

CMGTools/VVResonances/interface/FastJetInterface.h

@@ -30,12 +30,12 @@ class FastJetInterface {
   void makeExclusiveJets( int);
   void makeExclusiveJets( double);
   void makeExclusiveJetsUpTo(int);
-  void makeSubJets(unsigned int);
-  void makeSubJets(unsigned int, double);
-  void makeSubJets(unsigned int, int);
-  void makeSubJetsUpTo(unsigned int, int);
-  void prune(bool,double zcut,double rcutfactor);
-  void softDrop(bool,double beta,double zcut,double R0);
+  void makeSubJets(bool,unsigned int);
+  void makeSubJets(bool,unsigned int, double);
+  void makeSubJets(bool,unsigned int, int);
+  void makeSubJetsUpTo(bool,unsigned int, int);
+  void prune(bool,unsigned int ,double zcut,double rcutfactor);
+  void softDrop(bool,unsigned int,double beta,double zcut,double R0);
   bool  massDropTag( unsigned int, double&,double&);
 
 

CMGTools/VVResonances/python/analyzers/EventInterpretationBase.py

@@ -6,7 +6,7 @@
 from PhysicsTools.HeppyCore.utils.deltar import *
 import PhysicsTools.HeppyCore.framework.config as cfg
 from CMGTools.VVResonances.tools.PyJetToolbox import *
-
+from copy import copy
 
 
 
@@ -36,13 +36,18 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
           if self.cfg_comp.isMC:
               if hasattr(self.cfg_comp,'globalTag'):
                   self.jetReCalibrator = JetReCalibrator(self.cfg_comp.globalTag,self.cfg_ana.recalibrationType, False,cfg_ana.jecPath)
+                  self.jetReCalibratorFAT = JetReCalibrator(self.cfg_comp.globalTag,self.cfg_ana.recalibrationTypeFAT, False,cfg_ana.jecPath)
+
               else:
                   self.jetReCalibrator = JetReCalibrator(mcGT,self.cfg_ana.recalibrationType, False,cfg_ana.jecPath)
+                  self.jetReCalibratorFAT = JetReCalibrator(mcGT,self.cfg_ana.recalibrationTypeFAT, False,cfg_ana.jecPath)
           else:
               if hasattr(self.cfg_comp,'globalTag'):
                   self.jetReCalibrator = JetReCalibrator(self.cfg_comp.globalTag,self.cfg_ana.recalibrationType, True,cfg_ana.jecPath)
+                  self.jetReCalibratorFAT = JetReCalibrator(self.cfg_comp.globalTag,self.cfg_ana.recalibrationTypeFAT, True,cfg_ana.jecPath)
               else:    
                   self.jetReCalibrator = JetReCalibrator(dataGT,self.cfg_ana.recalibrationType, True,cfg_ana.jecPath)
+                  self.jetReCalibratorFAT = JetReCalibrator(dataGT,self.cfg_ana.recalibrationTypeFAT, True,cfg_ana.jecPath)          
 
         self.attachBTag = cfg_ana.attachBTag    
         if self.attachBTag:
@@ -69,16 +74,22 @@ def declareHandles(self):
 
     def removeLeptonFootPrint(self,leptons,cands):
         toRemove=[]
+        cList=list(cands)
         for lepton in leptons:
+
             for p in range(0,lepton.numberOfSourceCandidatePtrs()):
                 index = lepton.sourceCandidatePtr(p).key()
-                toRemove.append(cands[index])
-        return list(set(cands)-set(toRemove))
+                if not cands[index] in toRemove: 
+                    toRemove.append(cands[index])
+
+        for r in toRemove:
+            cList.remove(r)
+        return cList
 
 
     def matchSubJets(self,jets,genquarks):
         for j in jets:
-            for s in j.subjets:
+            for s in j.subjets_SD+j.subjets:
                 s.mcquark = None
                 s.matched=0
                 for g in genquarks:
@@ -118,21 +129,34 @@ def makeFatJets(self,cands):
         else:
             toolboxFat.setSoftDrop(self.cfg_ana.softdrop)
         # Lets cluster !! Fat jets first
-        fatJets=toolboxFat.inclusiveJets(100.0,True)
+        fatJets=toolboxFat.inclusiveJets(200.0,True)
         filtered = filter(self.selectFat,fatJets)
 
+        ##Apply JECS in SoftDrop and Pruned:
+        if self.jetReCalibratorFAT is not None:
+
+            prunedJets=[]
+            for j in fatJets:
+                prunedJets.append(j.softDropJet)
+                prunedJets.append(j.prunedJet)
+
+            self.jetReCalibratorFAT.correctAll(prunedJets, self.rho, self.shiftJEC,True,False,[0.,0.],[0.,0.,0.])
+            for p in prunedJets:
+                p.setRawFactor(1.0/p.corr)
         standardFatJets = self.handles['fatjets'].product()
 
         if self.attachBTag:
             for fat in filtered:
                 for standardFat in standardFatJets:
                     fat.btag = standardFat.bDiscriminator(self.btagDiscriminator)
-                for j in fat.subjets:
+                for j in fat.subjets_SD+fat.subjets:
                     for standard in self.handles['subjets'].product():
                         if deltaR(j.eta(),j.phi(),standard.eta(),standard.phi())<0.1:
                             j.btag = standard.bDiscriminator(self.btagDiscriminator)
                             break
 
+
+
         return filtered
 
     def makeSatelliteJets(self,cands):
@@ -142,9 +166,12 @@ def makeSatelliteJets(self,cands):
         toolbox.setSubjets(False,'inc',2)
         toolbox.setPruning(False)
         toolbox.setNtau(False)
-        toolbox.setSoftDrop(False)
-
+        toolbox.setSoftDrop(False)       
         unfiltered =  toolbox.inclusiveJets(30.0,False)
+
+
+
+
         if self.attachBTag:
             for j in unfiltered:
                     for standard in self.handles['jets'].product():
@@ -154,18 +181,63 @@ def makeSatelliteJets(self,cands):
 
         if self.jetReCalibrator is not None:
             self.jetReCalibrator.correctAll(unfiltered, self.rho, self.shiftJEC,True,False,[0.,0.],[0.,0.,0.])
-            return filter(lambda x: x.pt()>30, unfiltered)
+            for p in unfiltered:
+                p.setRawFactor(1.0/p.corr)
+
+
+            filtered=filter(lambda x: x.pt()>30, unfiltered)    
+            return filtered
         else:
             return unfiltered
 
     def removeJetFootPrint(self,jets,cands):
+
+
         toRemove=[]
         for j in jets:
             toRemove.extend(j.constituents)
-        return list(set(cands)-set(toRemove))
+        out=list(set(cands)-set(toRemove))
+        return out
 
 
     def topology(self,obj):
+
+        #first clean nicely the jets from leptons and subjets
+        objToClean=[]
+        if abs(obj['pair'].leg1.pdgId())==23:
+            objToClean.extend([obj['pair'].leg1.leg1,obj['pair'].leg1.leg2])
+#            print 'Z->ll', obj['pair'].leg1.leg1.pt(),obj['pair'].leg1.leg1.eta(),obj['pair'].leg1.leg1.phi(),obj['pair'].leg1.leg2.pt(),obj['pair'].leg1.leg2.eta(),obj['pair'].leg1.leg2.phi()
+        if abs(obj['pair'].leg1.pdgId())==24:
+            objToClean.extend([obj['pair'].leg1.leg1])
+#            print 'W->lnu', obj['pair'].leg1.leg1.pt(),obj['pair'].leg1.leg1.eta(),obj['pair'].leg1.leg1.phi()
+
+        if abs(obj['pair'].leg1.pdgId())==1:
+            objToClean.extend(obj['pair'].leg1.subjets)
+
+        if abs(obj['pair'].leg2.pdgId())==23:
+            objToClean.extend([obj['pair'].leg2.leg1,obj['pair'].leg2.leg2])
+        if abs(obj['pair'].leg2.pdgId())==24:
+            objToClean.extend([obj['pair'].leg2.leg1])
+        if abs(obj['pair'].leg2.pdgId())==1:
+            objToClean.extend(obj['pair'].leg2.subjets)
+#            print 'jet', obj['pair'].leg2.pt(),obj['pair'].leg2.eta(),obj['pair'].leg2.phi()
+
+        newJets=[]    
+        for j in  obj['satelliteJets']:
+            keep=True
+            for c in objToClean:
+                dr=deltaR(j.eta(),j.phi(),c.eta(),c.phi())
+#                print 'dr=',dr
+                if dr<0.4:
+                    keep=False
+                    break
+            if keep:
+#                print 'filtered satellite jet',j.pt(),j.eta(),j.phi() 
+                newJets.append(j)
+
+
+        obj['satelliteJets']=newJets            
+
         if len(obj['satelliteJets'])<2:
             obj['vbfDEta'] = -1.0
             obj['vbfMJJ'] = -1.0

CMGTools/VVResonances/python/analyzers/LeptonicVMaker.py

@@ -24,13 +24,110 @@ def makeDiLeptons(self,leptonList):
         for l1,l2 in combinations(leptonList,2):
             if  (l1.pdgId() == -l2.pdgId()):
                 pair = Pair(l1,l2,23)
+                if abs(l1.pdgId())==11:                  
+                    self.eeIsolationFootPrint(pair)
+                elif abs(l1.pdgId())==13:                      
+                    self.mumuIsolationFootPrint(pair)
                 m=pair.p4().mass()
                 if self.selectLLPair(pair):
 #                    print 'New Z with mass ',m
                     output.append(pair)
         return output            
 
 
+    def muIsolationFootPrint(self,z):#does nothing / for common interface
+        footPrintLeg1=0.0
+        z.leg1.cleanedChargedIso = max(z.leg1.pfIsolationR04().sumChargedHadronPt,0.0)                    
+        z.leg1.cleanedNeutralIsoDB=max( z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-z.leg1.puChargedHadronIsoR(0.4)/2, 0.0)
+        z.leg1.cleanedNeutralIsoRho=max(z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-z.leg1.rho*z.leg1.EffectiveArea04,0)        
+
+
+
+    def mumuIsolationFootPrint(self,z):
+        footPrintLeg1=0.0
+        footPrintLeg2=0.0
+        dr=deltaR(z.leg1.eta(),z.leg1.phi(),z.leg2.eta(),z.leg2.phi())
+        if dr<0.4 and not z.leg2.isPFMuon():
+            footPrintLeg1=z.leg2.pt()
+        if dr<0.4 and not z.leg1.isPFMuon():
+            footPrintLeg2=z.leg2.pt()
+
+        z.leg1.cleanedChargedIso = max(z.leg1.pfIsolationR04().sumChargedHadronPt-footPrintLeg1,0.0)                    
+        z.leg2.cleanedChargedIso = max(z.leg2.pfIsolationR04().sumChargedHadronPt-footPrintLeg2,0.0)                    
+        z.leg1.cleanedNeutralIsoDB=max( z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-z.leg1.puChargedHadronIsoR(0.4)/2, 0.0)
+        z.leg1.cleanedNeutralIsoRho=max(z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-z.leg1.rho*z.leg1.EffectiveArea04,0)        
+        z.leg2.cleanedNeutralIsoDB=max( z.leg2.neutralHadronIsoR(0.4)+z.leg2.photonIsoR(0.4)-z.leg2.puChargedHadronIsoR(0.4)/2, 0.0)
+        z.leg2.cleanedNeutralIsoRho=max(z.leg2.neutralHadronIsoR(0.4)+z.leg2.photonIsoR(0.4)-z.leg2.rho*z.leg2.EffectiveArea04,0)        
+
+
+
+    def eIsolationFootPrint(self,z):
+        footPrintChargedLeg1=0.0
+        footPrintNeutralLeg1=0.0
+
+        if not z.leg1.isPF():
+
+            for i in range(0,z.leg1.associatedPackedPFCandidates().size()):
+                c=z.leg1.associatedPackedPFCandidates()[i]
+                dr = deltaR(z.leg1.eta(),z.leg1.phi(),c.eta(),c.phi())
+                if z.leg1.isEB() or (dr>0.08 and dr<0.4):
+                    if c.charge()>0:
+                        footPrintChargedLeg1=footPrintChargedLeg1+c.pt()
+                    else:    
+                        footPrintNeutralLeg1=footPrintNeutralLeg1+c.pt()
+
+        z.leg1.cleanedChargedIso = max(z.leg1.chargedHadronIsoR(0.4)-footPrintChargedLeg1,0.0)                    
+        z.leg1.cleanedNeutralIsoDB=max( z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-footPrintNeutralLeg1-z.leg1.puChargedHadronIsoR(0.4)/2, 0.0)
+        z.leg1.cleanedNeutralIsoRho=max(z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-footPrintNeutralLeg1-z.leg1.rho*z.leg1.EffectiveArea04,0)        
+
+    def eeIsolationFootPrint(self,z):
+        footPrintChargedLeg1=0.0
+        footPrintChargedLeg2=0.0
+        footPrintNeutralLeg1=0.0
+        footPrintNeutralLeg2=0.0
+
+        if not z.leg1.isPF():
+
+            for i in range(0,z.leg1.associatedPackedPFCandidates().size()):
+                c=z.leg1.associatedPackedPFCandidates()[i]
+                dr = deltaR(z.leg1.eta(),z.leg1.phi(),c.eta(),c.phi())
+                if z.leg1.isEB() or (dr>0.08 and dr<0.4):
+                    if c.charge()>0:
+                        footPrintChargedLeg1=footPrintChargedLeg1+c.pt()
+                    else:    
+                        footPrintNeutralLeg1=footPrintNeutralLeg1+c.pt()
+                dr = deltaR(z.leg2.eta(),z.leg2.phi(),c.eta(),c.phi())
+                if z.leg2.isEB() or (dr>0.08 and dr<0.4):
+                    if c.charge()>0:
+                        footPrintChargedLeg2=footPrintChargedLeg2+c.pt()
+                    else:    
+                        footPrintNeutralLeg2=footPrintNeutralLeg2+c.pt()
+
+        if not z.leg2.isPF():
+            for i in range(0,z.leg2.associatedPackedPFCandidates().size()):
+                c=z.leg2.associatedPackedPFCandidates()[i]
+                dr = deltaR(z.leg2.eta(),z.leg2.phi(),c.eta(),c.phi())
+                if z.leg2.isEB() or (dr>0.08 and dr<0.4):
+                    if c.charge()>0:
+                        footPrintChargedLeg2=footPrintChargedLeg2+c.pt()
+                    else:    
+                        footPrintNeutralLeg2=footPrintNeutralLeg2+c.pt()
+                dr = deltaR(z.leg1.eta(),z.leg1.phi(),c.eta(),c.phi())
+                if z.leg1.isEB() or (dr>0.08 and dr<0.4):
+                    if c.charge()>0:
+                        footPrintChargedLeg1=footPrintChargedLeg1+c.pt()
+                    else:    
+                        footPrintNeutralLeg1=footPrintNeutralLeg1+c.pt()
+
+
+        z.leg1.cleanedChargedIso = max(z.leg1.chargedHadronIsoR(0.4)-footPrintChargedLeg1,0.0)                    
+        z.leg2.cleanedChargedIso = max(z.leg1.chargedHadronIsoR(0.4)-footPrintChargedLeg2,0.0)                    
+        z.leg1.cleanedNeutralIsoDB=max( z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-footPrintNeutralLeg1-z.leg1.puChargedHadronIsoR(0.4)/2, 0.0)
+        z.leg1.cleanedNeutralIsoRho=max(z.leg1.neutralHadronIsoR(0.4)+z.leg1.photonIsoR(0.4)-footPrintNeutralLeg1-z.leg1.rho*z.leg1.EffectiveArea04,0)        
+        z.leg2.cleanedNeutralIsoDB=max( z.leg2.neutralHadronIsoR(0.4)+z.leg2.photonIsoR(0.4)-footPrintNeutralLeg2-z.leg2.puChargedHadronIsoR(0.4)/2, 0.0)
+        z.leg2.cleanedNeutralIsoRho=max(z.leg2.neutralHadronIsoR(0.4)+z.leg2.photonIsoR(0.4)--footPrintNeutralLeg2-z.leg2.rho*z.leg2.EffectiveArea04,0)        
+
+
     def simpleWKinematicFit(self,pair):
         MW=80.390
 
@@ -84,7 +181,7 @@ def simpleWKinematicFit(self,pair):
     def defaultWKinematicFit(self,pair):
         MW=80.390
 
-        muonLV = ROOT.TLorentzVector(pair.leg1.px(),pair.leg1.py(),pair.leg1.pz(),pair.leg1.energy())
+        muLonLV = ROOT.TLorentzVector(pair.leg1.px(),pair.leg1.py(),pair.leg1.pz(),pair.leg1.energy())
         metLV = ROOT.TLorentzVector(pair.leg2.px(),pair.leg2.py(),pair.leg2.pz(),pair.leg2.energy())
 
 
@@ -124,6 +221,10 @@ def makeLeptonsMET(self,leptonList,MET):
         for l1 in leptonList:
             pair = Pair(l1,MET,l1.charge()*24)
             mt=pair.mt()
+            if abs(l1.pdgId())==11:
+                self.eIsolationFootPrint(pair)
+            elif     abs(l1.pdgId())==13:
+                self.muIsolationFootPrint(pair)
             if self.selectLNuPair(pair):
 #                self.defaultWKinematicFit(pair)
                 self.simpleWKinematicFit(pair)
@@ -145,14 +246,18 @@ def process(self, event):
 
 
         #now make Z first . for the remaining leptons after Z make W
-        leptonsSet = set(event.selectedLeptons)
+        leptons=list(event.selectedLeptons)
         used = []
+
+
         for z in event.LL:
             used.extend([z.leg1,z.leg2])
-        usedSet = set(used)
 
-        remaining = leptonsSet-usedSet
-        event.LNu = self.makeLeptonsMET(list(remaining),event.met)
+        for u in used:
+            leptons.remove(u)
+
+
+        event.LNu = self.makeLeptonsMET(leptons,event.met)
         return True