Skip to content

Commit

Permalink
Merge pull request cms-sw#77 from arizzi/vhbbQGL
Browse files Browse the repository at this point in the history 
Make QGL faster being on demand only for VBF jets
  • Loading branch information
@arizzi
arizzi committed Jun 23, 2015
2 parents 584ffeb + a70243a commit c1dd86a
Show file tree
Hide file tree
Showing 5 changed files with 104 additions and 75 deletions.
72 changes: 3 additions & 69 deletions PhysicsTools/Heppy/python/analyzers/objects/JetAnalyzer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ def declareHandles(self):

def beginLoop(self, setup):
super(JetAnalyzer,self).beginLoop(setup)

def process(self, event):
self.readCollections( event.input )
rho = float(self.handles['rho'].product()[0])
Expand Down Expand Up @@ -130,8 +130,8 @@ def process(self, event):
if self.testJetID (jet ):

if(self.cfg_ana.doQG):
self.computeQGvars(jet)
jet.qgl = self.qglcalc.computeQGLikelihood(jet, rho)
jet.qgl_calc = self.qglcalc.computeQGLikelihood
jet.qgl_rho = rho


self.jets.append(jet)
Expand Down Expand Up @@ -262,71 +262,6 @@ def testJetNoID( self, jet ):
return jet.pt() > self.cfg_ana.jetPt and \
abs( jet.eta() ) < self.cfg_ana.jetEta;

def computeQGvars(self, jet):

jet.mult = 0
sum_weight = 0.
sum_pt = 0.
sum_deta = 0.
sum_dphi = 0.
sum_deta2 = 0.
sum_detadphi = 0.
sum_dphi2 = 0.



for ii in range(0, jet.numberOfDaughters()) :

part = jet.daughter(ii)

if part.charge() == 0 : # neutral particles

if part.pt() < 1.: continue

else : # charged particles

if part.trackHighPurity()==False: continue
if part.fromPV()<=1: continue


jet.mult += 1

deta = part.eta() - jet.eta()
dphi = deltaPhi(part.phi(), jet.phi())
partPt = part.pt()
weight = partPt*partPt
sum_weight += weight
sum_pt += partPt
sum_deta += deta*weight
sum_dphi += dphi*weight
sum_deta2 += deta*deta*weight
sum_detadphi += deta*dphi*weight
sum_dphi2 += dphi*dphi*weight




a = 0.
b = 0.
c = 0.

if sum_weight > 0 :
jet.ptd = math.sqrt(sum_weight)/sum_pt
ave_deta = sum_deta/sum_weight
ave_dphi = sum_dphi/sum_weight
ave_deta2 = sum_deta2/sum_weight
ave_dphi2 = sum_dphi2/sum_weight
a = ave_deta2 - ave_deta*ave_deta
b = ave_dphi2 - ave_dphi*ave_dphi
c = -(sum_detadphi/sum_weight - ave_deta*ave_dphi)
else: jet.ptd = 0.

delta = math.sqrt(math.fabs((a-b)*(a-b)+4.*c*c))

if a+b-delta > 0: jet.axis2 = -math.log(math.sqrt(0.5*(a+b-delta)))
else: jet.axis2 = -1.


def jetFlavour(self,event):
def isFlavour(x,f):
id = abs(x.pdgId())
Expand Down Expand Up @@ -362,7 +297,6 @@ def isFlavour(x,f):

self.heaviestQCDFlavour = 5 if len(self.bqObjects) else (4 if len(self.cqObjects) else 1);


def matchJets(self, event, jets):
match = matchObjectCollection2(jets,
event.genbquarks + event.genwzquarks,
Expand Down
8 changes: 4 additions & 4 deletions PhysicsTools/Heppy/python/analyzers/objects/autophobj.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,10 +178,10 @@
jetTypeExtra = NTupleObjectType("jetExtra", baseObjectTypes = [ jetType ], variables = [
NTupleVariable("area", lambda x : x.jetArea(), help="Catchment area of jet"),
# QG variables:
NTupleVariable("qgl", lambda x : getattr(x,'qgl', 0) , float, mcOnly=False,help="QG Likelihood"),
NTupleVariable("ptd", lambda x : getattr(x,'ptd', 0), float, mcOnly=False,help="QG input variable: ptD"),
NTupleVariable("axis2", lambda x : getattr(x,'axis2', 0) , float, mcOnly=False,help="QG input variable: axis2"),
NTupleVariable("mult", lambda x : getattr(x,'mult', 0) , int, mcOnly=False,help="QG input variable: total multiplicity"),
NTupleVariable("qgl", lambda x :x.qgl() , float, mcOnly=False,help="QG Likelihood"),
NTupleVariable("ptd", lambda x : getattr(x.computeQGvars(),'ptd', 0), float, mcOnly=False,help="QG input variable: ptD"),
NTupleVariable("axis2", lambda x : getattr(x.computeQGvars(),'axis2', 0) , float, mcOnly=False,help="QG input variable: axis2"),
NTupleVariable("mult", lambda x : getattr(x.computeQGvars(),'mult', 0) , int, mcOnly=False,help="QG input variable: total multiplicity"),
NTupleVariable("partonId", lambda x : getattr(x,'partonId', 0), int, mcOnly=True, help="parton flavour (manually matching to status 23 particles)"),
NTupleVariable("partonMotherId", lambda x : getattr(x,'partonMotherId', 0), int, mcOnly=True, help="parton flavour (manually matching to status 23 particles)"),
NTupleVariable("nLeptons", lambda x : len(x.leptons) if hasattr(x,'leptons') else 0 , float, mcOnly=False,help="Number of associated leptons"),
Expand Down
83 changes: 83 additions & 0 deletions PhysicsTools/Heppy/python/physicsobjects/Jet.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
from PhysicsTools.Heppy.physicsobjects.PhysicsObject import *
from PhysicsTools.HeppyCore.utils.deltar import deltaPhi
import math

loose_WP = [
(0, 2.5, -0.8),
Expand Down Expand Up @@ -122,6 +124,87 @@ def leadTrackPt(self):
return lt.pt()
else :
return 0.
def qgl(self) :
if not hasattr(self,"qgl_value") :
if hasattr(self,"qgl_rho") : #check if qgl calculator is configured
self.computeQGvars()
self.qgl_value=self.qgl_calc(self,self.qgl_rho)
else :
self.qgl_value=0. #if no qgl calculator configured

return self.qgl_value

def computeQGvars(self):
#return immediately if qgvars already computed or if qgl is disabled
if not hasattr(self,"qgl_rho") or getattr(self,"hasQGVvars",False) :
return self
self.hasQGvars = True

jet = self
jet.mult = 0
sum_weight = 0.
sum_pt = 0.
sum_deta = 0.
sum_dphi = 0.
sum_deta2 = 0.
sum_detadphi = 0.
sum_dphi2 = 0.



for ii in range(0, jet.numberOfDaughters()) :

part = jet.daughter(ii)

if part.charge() == 0 : # neutral particles

if part.pt() < 1.: continue

else : # charged particles

if part.trackHighPurity()==False: continue
if part.fromPV()<=1: continue


jet.mult += 1

deta = part.eta() - jet.eta()
dphi = deltaPhi(part.phi(), jet.phi())
partPt = part.pt()
weight = partPt*partPt
sum_weight += weight
sum_pt += partPt
sum_deta += deta*weight
sum_dphi += dphi*weight
sum_deta2 += deta*deta*weight
sum_detadphi += deta*dphi*weight
sum_dphi2 += dphi*dphi*weight




a = 0.
b = 0.
c = 0.

if sum_weight > 0 :
jet.ptd = math.sqrt(sum_weight)/sum_pt
ave_deta = sum_deta/sum_weight
ave_dphi = sum_dphi/sum_weight
ave_deta2 = sum_deta2/sum_weight
ave_dphi2 = sum_dphi2/sum_weight
a = ave_deta2 - ave_deta*ave_deta
b = ave_dphi2 - ave_dphi*ave_dphi
c = -(sum_detadphi/sum_weight - ave_deta*ave_dphi)
else: jet.ptd = 0.

delta = math.sqrt(math.fabs((a-b)*(a-b)+4.*c*c))

if a+b-delta > 0: jet.axis2 = -math.log(math.sqrt(0.5*(a+b-delta)))
else: jet.axis2 = -1.
return jet



class GenJet( PhysicsObject):
pass
Expand Down
6 changes: 5 additions & 1 deletion VHbbAnalysis/Heppy/python/VHbbAnalyzer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,13 @@ def doSoftActivity(self,event) :
event.jetsForVBF = [x for x in event.cleanJetsAll if self.cfg_ana.higgsJetsPreSelection(x) ]
#compute SoftActivity only for events passing VBF selection
if len(event.jetsForVBF) < 4 or event.jetsForVBF[0] < 70 or event.jetsForVBF[1] < 55 or event.jetsForVBF[2] < 35 or event.jetsForVBF[3] < 20 :
return
return
event.jetsForVBF.sort(key=lambda x:x.pt(),reverse=True)
event.jetsForVBF=event.jetsForVBF[:4]

#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]
j1=event.bJetsForVBF[0]
j2=event.bJetsForVBF[1]
Expand Down
10 changes: 9 additions & 1 deletion VHbbAnalysis/Heppy/python/vhbbobj.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,10 +105,18 @@
NTupleVariable("vtxPosY", lambda x : x.userFloat("vtxPosY"), mcOnly=False, help="Y coord of vertex from btag"),
NTupleVariable("vtxPosZ", lambda x : x.userFloat("vtxPosZ"), mcOnly=False, help="Z coord of vertex from btag"),
# QG variables:
NTupleVariable("qgl", lambda x : getattr(x,'qgl', 0) , float, mcOnly=False,help="QG Likelihood"),
# this computes for all
# NTupleVariable("qgl", lambda x :x.qgl() , float, mcOnly=False,help="QG Likelihood"),
# NTupleVariable("ptd", lambda x : getattr(x.computeQGvars(),'ptd', 0), float, mcOnly=False,help="QG input variable: ptD"),
# NTupleVariable("axis2", lambda x : getattr(x.computeQGvars(),'axis2', 0) , float, mcOnly=False,help="QG input variable: axis2"),
# NTupleVariable("mult", lambda x : getattr(x.computeQGvars(),'mult', 0) , int, mcOnly=False,help="QG input variable: total multiplicity"),

# this only read qgl if it was explicitelly computed in the code
NTupleVariable("qgl", lambda x : getattr(x,'qgl_value',0) , float, mcOnly=False,help="QG Likelihood"),
NTupleVariable("ptd", lambda x : getattr(x,'ptd', 0), float, mcOnly=False,help="QG input variable: ptD"),
NTupleVariable("axis2", lambda x : getattr(x,'axis2', 0) , float, mcOnly=False,help="QG input variable: axis2"),
NTupleVariable("mult", lambda x : getattr(x,'mult', 0) , int, mcOnly=False,help="QG input variable: total multiplicity"),

NTupleVariable("numberOfDaughters", lambda x : x.numberOfDaughters(), int, mcOnly=False,help="number of daughters"),
])

Expand Down

0 comments on commit c1dd86a

Please sign in to comment.