Merge pull request cms-sw#87 from arizzi/severalFixes

Several fixes

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

@@ -2,7 +2,7 @@
 from PhysicsTools.Heppy.analyzers.core.Analyzer import Analyzer
 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
 from PhysicsTools.Heppy.physicsobjects.PhysicsObjects import Jet
-from PhysicsTools.HeppyCore.utils.deltar import deltaR2, deltaPhi, matchObjectCollection, matchObjectCollection2, bestMatch
+from PhysicsTools.HeppyCore.utils.deltar import deltaR2, deltaPhi, matchObjectCollection, matchObjectCollection2, bestMatch,matchObjectCollection3
 from PhysicsTools.Heppy.physicsutils.JetReCalibrator import JetReCalibrator
 import PhysicsTools.HeppyCore.framework.config as cfg
 
@@ -208,7 +208,16 @@ def process(self, event):
 
             if self.cfg_ana.cleanGenJetsFromPhoton:
                 self.cleanGenJets = cleanNearestJetOnly(self.cleanGenJets, photons, self.jetLepDR)
-
+	    if hasattr(self.cfg_ana,"genNuSelection") :
+		jetNus=[x for x in event.genParticles if abs(x.pdgId()) in [12,14,16] and self.cfg_ana.genNuSelection(x) ]
+	 	pairs= matchObjectCollection3 ( self.cleanGenJets, jetNus, 0.4)
+		for (genJet,nu) in pairs.iteritems() :
+		     if nu is not None :
+			if not hasattr(genJet,"nu") :
+				genJet.nu=nu.p4()
+			else :
+				genJet.nu+=nu.p4()
+
 
             #event.nGenJets25 = 0
             #event.nGenJets25Cen = 0
@@ -366,6 +375,7 @@ def smearJets(self, event, jets):
     alwaysCleanPhotons = False,
     jecPath = "",
     cleanGenJetsFromPhoton = False,
+    genNuSelection = lambda nu : True, #FIXME: add here check for ispromptfinalstate
     collectionPostFix = ""
     )
 )

VHbbAnalysis/Heppy/python/VHbbAnalyzer.py

@@ -61,15 +61,18 @@ def doVBF(self,event) :
 	#compute QGL here for VBF jets if passing VBF pre-selection 
         map(lambda x :x.qgl(),event.jetsForVBF)
 
-        event.bJetsForVBF=sorted(event.jetsForVBF,key = lambda jet : jet.btag('combinedInclusiveSecondaryVertexV2BJetTags'), reverse=True)[:2]
+        event.bJetsForVBF=sorted(event.jetsForVBF,key = lambda jet : jet.btag('pfCombinedInclusiveSecondaryVertexV2BJetTags'), reverse=True)[:2]
         j1=event.bJetsForVBF[0]
         j2=event.bJetsForVBF[1]
-	event.softActivityJets=self.softActivity(event,j1,j2,event.jetsForVBF)
+#print "VBF"
+#print event.selectedElectrons,event.selectedMuons,
+	event.softActivityJets=self.softActivity(event,j1,j2,event.jetsForVBF+event.selectedElectrons+event.selectedMuons)
 
     def doSoftActivityVH(self,event) :
         j1=event.hJetsCSV[0]
         j2=event.hJetsCSV[1]
-        event.softActivityVHJets=self.softActivity(event,j1,j2,event.hJetsCSV+event.vLeptons)
+#print "VH"
+        event.softActivityVHJets=self.softActivity(event,j1,j2,event.hJetsCSV+event.selectedElectrons+event.selectedMuons)
 
 
     def softActivity(self,event,j1,j2,excludedJets) :
@@ -138,22 +141,22 @@ def doFakeMET(self,event):
 
     def doHiggsHighCSV(self,event) :
         #leading csv interpretation
-        event.hJetsCSV=sorted(event.jetsForHiggs,key = lambda jet : jet.btag('combinedInclusiveSecondaryVertexV2BJetTags'), reverse=True)[0:2]
+        event.hJetsCSV=sorted(event.jetsForHiggs,key = lambda jet : jet.btag('pfCombinedInclusiveSecondaryVertexV2BJetTags'), reverse=True)[0:2]
         event.aJetsCSV = [x for x in event.cleanJets if x not in event.hJetsCSV]
-        event.hjidxCSV=[event.cleanJets.index(x) for x in event.hJetsCSV ]
-        event.ajidxCSV=[event.cleanJets.index(x) for x in event.aJetsCSV ]
+        event.hjidxCSV=[event.cleanJetsAll.index(x) for x in event.hJetsCSV ]
+        event.ajidxCSV=[event.cleanJetsAll.index(x) for x in event.aJetsCSV ]
         event.aJetsCSV+=event.cleanJetsFwd
         event.HCSV = event.hJetsCSV[0].p4()+event.hJetsCSV[1].p4()
 
     def doHiggsHighPt(self,event) :
         #highest pair interpretations
         event.hJets=list(max(itertools.combinations(event.jetsForHiggs,2), key = lambda x : (x[0].p4()+x[1].p4()).pt() ))
         event.aJets = [x for x in event.cleanJets if x not in event.hJets]
-        event.hjidx=[event.cleanJets.index(x) for x in event.hJets ]
-        event.ajidx=[event.cleanJets.index(x) for x in event.aJets ]
+        event.hjidx=[event.cleanJetsAll.index(x) for x in event.hJets ]
+        event.ajidx=[event.cleanJetsAll.index(x) for x in event.aJets ]
         event.aJets+=event.cleanJetsFwd
-        hJetsByCSV = sorted(event.hJets , key =  lambda jet : jet.btag('combinedInclusiveSecondaryVertexV2BJetTags'), reverse=True)
-        event.hjidxDiJetPtByCSV = [event.cleanJets.index(x) for x in hJetsByCSV]
+        hJetsByCSV = sorted(event.hJets , key =  lambda jet : jet.btag('pfCombinedInclusiveSecondaryVertexV2BJetTags'), reverse=True)
+        event.hjidxDiJetPtByCSV = [event.cleanJetsAll.index(x) for x in hJetsByCSV]
         event.H = event.hJets[0].p4()+event.hJets[1].p4()
 
 
@@ -172,8 +175,8 @@ def doHiggs3cj(self,event) :
            event.hJets3cj.append(j3)
 #          event.hJets3cj=list(min(itertools.combinations(event.jetsForHiggs,3), key = lambda x : ( x[2]>15 and (deltaR( x[0].eta(), x[0].phi(), x[2].eta(), x[2].phi()) +  deltaR( x[1].eta(), x[1].phi(), x[2].eta(), x[2].phi()) ) ) ))
            event.aJets3cj = [x for x in event.cleanJets if x not in event.hJets3cj]
-           event.hjidx3cj=[event.cleanJets.index(x) for x in event.hJets3cj ]
-           event.ajidx3cj=[event.cleanJets.index(x) for x in event.aJets3cj ]
+           event.hjidx3cj=[event.cleanJetsAll.index(x) for x in event.hJets3cj ]
+           event.ajidx3cj=[event.cleanJetsAll.index(x) for x in event.aJets3cj ]
            event.aJets3cj+=event.cleanJetsFwd
            event.H3cj = event.hJets3cj[0].p4()+event.hJets3cj[1].p4()+event.hJets3cj[2].p4()
            event.minDr3 = min(deltaR( event.hJets3cj[0].eta(), event.hJets3cj[0].phi(), event.hJets3cj[2].eta(), event.hJets3cj[2].phi()) ,  deltaR( event.hJets3cj[1].eta(), event.hJets3cj[1].phi(), event.hJets3cj[2].eta(), event.hJets3cj[2].phi()))  
@@ -315,6 +318,7 @@ def initOutputs (self,event) :
         event.hjidxDiJetPtByCSV = []
         event.softActivityJets=[]
     def process(self, event):
+	#print "Event number",event.iEv
         self.readCollections( event.input )
         self.inputCounter.Fill(1)
         if self.cfg_comp.isMC:
@@ -340,6 +344,7 @@ def process(self, event):
         if event.Vtype < 0 and not ( sum(x.pt() > 30 for x in event.jetsForHiggs) >= 4 or sum(x.pt() for x in event.jetsForHiggs[:4]) > 160 ):
                 return self.cfg_ana.passall
 
+        map(lambda x :x.qgl(),event.jetsForHiggs[:4])
 	self.doHiggsHighCSV(event)
 	self.doHiggsHighPt(event)
         self.doHiggs3cj(event)

VHbbAnalysis/Heppy/python/vhbbobj.py

@@ -43,7 +43,7 @@
     # TTH-id related variables
     NTupleVariable("mvaTTH",     lambda lepton : lepton.mvaValue if hasattr(lepton,'mvaValue') else -1, help="Lepton MVA (ttH version)"),
     NTupleVariable("jetPtRatio", lambda lepton : lepton.pt()/lepton.jet.pt() if hasattr(lepton,'jet') else -1, help="pt(lepton)/pt(nearest jet)"),
-    NTupleVariable("jetBTagCSV", lambda lepton : lepton.jet.btag('combinedSecondaryVertexBJetTags') if hasattr(lepton,'jet') and hasattr(lepton.jet, 'btag') else -99, help="btag of nearest jet"),
+    NTupleVariable("jetBTagCSV", lambda lepton : lepton.jet.btag('pfCombinedInclusiveSecondaryVertexV2BJetTags') if hasattr(lepton,'jet') and hasattr(lepton.jet, 'btag') else -99, help="btag of nearest jet"),
     NTupleVariable("jetDR",      lambda lepton : deltaR(lepton.eta(),lepton.phi(),lepton.jet.eta(),lepton.jet.phi()) if hasattr(lepton,'jet') else -1, help="deltaR(lepton, nearest jet)"),
     NTupleVariable("pfRelIso03",      lambda ele : (ele.pfIsolationVariables().sumChargedHadronPt + max(ele.pfIsolationVariables().sumNeutralHadronEt + ele.pfIsolationVariables().sumPhotonEt - 0.5 * ele.pfIsolationVariables().sumPUPt,0.0)) / ele.pt()  if abs(ele.pdgId()) == 11 else -1, help="0.3 particle based iso"),
     NTupleVariable("pfRelIso04",      lambda mu : (mu.pfIsolationR04().sumChargedHadronPt + max( mu.pfIsolationR04().sumNeutralHadronEt + mu.pfIsolationR04().sumPhotonEt - 0.5 * mu.pfIsolationR04().sumPUPt,0.0)) / mu.pt() if abs(mu.pdgId()) == 13 else -1, help="0.4 particle based iso"),
@@ -198,7 +198,7 @@
     NTupleVariable("sip3d", lambda x : x.d3d.significance(), help="S_{ip3d} with respect to PV (absolute value)"),
     NTupleVariable("cosTheta", lambda x : x.cosTheta, help="Cosine of the angle between the 3D displacement and the momentum"),
     NTupleVariable("jetPt",  lambda x : x.jet.pt() if x.jet != None else 0, help="pT of associated jet"),
-    NTupleVariable("jetBTag",  lambda x : x.jet.btag('combinedSecondaryVertexBJetTags') if x.jet != None else -99, help="CSV b-tag of associated jet"),
+    NTupleVariable("jetBTag",  lambda x : x.jet.btag('pfCombinedInclusiveSecondaryVertexV2BJetTags') if x.jet != None else -99, help="CSV b-tag of associated jet"),
     NTupleVariable("mcMatchNTracks", lambda x : x.mcMatchNTracks, int, mcOnly=True, help="Number of mc-matched tracks in SV"),
     NTupleVariable("mcMatchNTracksHF", lambda x : x.mcMatchNTracksHF, int, mcOnly=True, help="Number of mc-matched tracks from b/c in SV"),
     NTupleVariable("mcMatchFraction", lambda x : x.mcMatchFraction, mcOnly=True, help="Fraction of mc-matched tracks from b/c matched to a single hadron (or -1 if mcMatchNTracksHF < 2)"),
@@ -222,7 +222,7 @@
     NTupleVariable("svSip3d", lambda x : x.sv.d3d.significance() if x.sv else -99., help="SV: S_{ip3d} with respect to PV (absolute value)"),
     NTupleVariable("svCosTheta", lambda x : x.sv.cosTheta if x.sv else -99., help="SV: Cosine of the angle between the 3D displacement and the momentum"),
     NTupleVariable("jetPt",  lambda x : x.jet.pt() if x.jet != None else 0, help="Jet: pT"),
-    NTupleVariable("jetBTag",  lambda x : x.jet.btag('combinedSecondaryVertexBJetTags') if x.jet != None else -99, help="CSV b-tag of associated jet"),
+    NTupleVariable("jetBTag",  lambda x : x.jet.btag('pfCombinedInclusiveSecondaryVertexV2BJetTags') if x.jet != None else -99, help="CSV b-tag of associated jet"),
 ])
 shiftedMetType= NTupleObjectType("shiftedMetType", baseObjectTypes=[twoVectorType], variables=[
     NTupleVariable("sumEt", lambda x : x.sumEt() ),
@@ -249,6 +249,10 @@
     NTupleVariable("numCHadronsFromTop", lambda x : x.numCHadronsFromTop, int, mcOnly=True, help="number of matched c hadrons from top quark decay"),
     NTupleVariable("numBHadronsAfterTop", lambda x : x.numBHadronsAfterTop, int, mcOnly=True, help="number of matched b hadrons after top quark decay"),
     NTupleVariable("numCHadronsAfterTop", lambda x : x.numCHadronsAfterTop, int, mcOnly=True, help="number of matched c hadrons after top quark decay"),
+    NTupleVariable("wNuPt", lambda x : (x.p4()+x.nu).pt() if hasattr(x,"nu") else x.p4().pt() ,float, mcOnly=True, help="pt of jet adding back the neutrinos"),
+    NTupleVariable("wNuEta", lambda x : (x.p4()+x.nu).eta() if hasattr(x,"nu") else x.p4().eta() ,float, mcOnly=True, help="eta of jet adding back the neutrinos"),
+    NTupleVariable("wNuPhi", lambda x : (x.p4()+x.nu).phi() if hasattr(x,"nu") else x.p4().phi() ,float, mcOnly=True, help="phi of jet adding back the neutrinos"),
+
 ])
 
 softActivityType = NTupleObjectType("softActivity", baseObjectTypes = [  ], variables = [

VHbbAnalysis/Heppy/src/FastSoftActivity.cc

@@ -52,6 +52,12 @@ FastSoftActivity::FastSoftActivity(const std::vector<LorentzVector> & objects, d
 std::vector<math::XYZTLorentzVector> FastSoftActivity::makeP4s(const std::vector<fastjet::PseudoJet> &jets) {
   std::vector<math::XYZTLorentzVector> JetObjectsAll;
   for (const fastjet::PseudoJet & pj : jets) {
+/*    std::vector<fastjet::PseudoJet> constituents = pj.constituents();
+    std::cout << "Constituents for " << pj.pt() << " " ;
+    for (unsigned j = 0; j < constituents.size(); j++) {
+		std::cout << constituents[j].user_index() << " ";
+    }
+    std::cout << std::endl;*/
     JetObjectsAll.push_back( LorentzVector( pj.px(), pj.py(), pj.pz(), pj.e() ) );
   }
   return JetObjectsAll;

VHbbAnalysis/Heppy/test/crab/heppy_crab_script.sh

@@ -34,8 +34,10 @@ edmProvDump
 
 # Update library path
 # Needed so recompiled modules are found
-export LD_LIBRARY_PATH=./lib/slc6_amd64_gcc481:$LD_LIBRARY_PATH 
-
+#export LD_LIBRARY_PATH=./lib/slc6_amd64_gcc481:$LD_LIBRARY_PATH 
+eval `scram runtime -sh`
+echo "LD LIBRARY PATH IS"
+echo $LD_LIBRARY_PATH
 
 # Move JEC files into flace
 mkdir jec
@@ -45,3 +47,5 @@ mv PHYS14_V4_MC_L3Absolute_AK4PFchs.txt jec/
 mkdir csv
 mv csv*root csv/
 python heppy_crab_script.py $1
+echo "======================== CMSRUN LOG ============================"
+cat Output/cmsRun.log