Merge pull request #17 from gpetruc/heppy_7_2_2_patch2

Hi Giovanni, 
Thank you very much. I'll merge faster next time.

PhysicsTools/Heppy/interface/IsolationComputer.h

@@ -0,0 +1,64 @@
+#ifndef PhysicsTools_Heppy_IsolationComputer_h
+#define PhysicsTools_Heppy_IsolationComputer_h
+
+#include <vector>
+#include "DataFormats/PatCandidates/interface/PackedCandidate.h"
+
+namespace heppy {
+class IsolationComputer {
+    public:
+        /// Create the calculator; optionally specify a cone for computing deltaBeta weights
+        IsolationComputer(float weightCone=-1) : weightCone_(weightCone) {}
+
+        /// Self-veto policy
+        enum SelfVetoPolicy {
+            selfVetoNone=0, selfVetoAll=1, selfVetoFirst=2
+        };
+        /// Initialize with the list of packed candidates (note: clears also all vetos)
+        void setPackedCandidates(const std::vector<pat::PackedCandidate> & all, int fromPV_thresh=1, float dz_thresh=9999., bool also_leptons=false) ;
+
+
+        /// veto footprint from this candidate, for the isolation of all candidates and also for calculation of neutral weights (if used)
+        void addVetos(const reco::Candidate &cand) ;
+
+        /// clear all vetos
+        void clearVetos() ;
+
+        /// Isolation from charged from the PV 
+        float chargedAbsIso(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from charged from PU
+        float puAbsIso(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from all neutrals (uncorrected)
+        float neutralAbsIsoRaw(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from neutral hadrons (uncorrected)
+        float neutralHadAbsIsoRaw(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from photons (uncorrected)
+        float photonAbsIsoRaw(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from all neutrals (with weights)
+        float neutralAbsIsoWeighted(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from neutral hadrons (with weights)
+        float neutralHadAbsIsoWeighted(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+
+        /// Isolation from photons (with weights)
+        float photonAbsIsoWeighted(const reco::Candidate &cand, float dR, float innerR=0, float threshold=0, SelfVetoPolicy selfVeto=selfVetoAll) const ;
+    protected:
+        const std::vector<pat::PackedCandidate> * allcands_;
+        float weightCone_;
+        // collections of objects, sorted in eta
+        std::vector<const pat::PackedCandidate *> charged_, neutral_, pileup_;
+        mutable std::vector<float> weights_;
+        std::vector<const reco::Candidate *> vetos_;
+
+        float isoSumRaw(const std::vector<const pat::PackedCandidate *> & cands, const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto, int pdgId=-1) const ;
+        float isoSumNeutralsWeighted(const reco::Candidate &cand, float dR, float innerR, float threshold, SelfVetoPolicy selfVeto, int pdgId=-1) const ;
+};
+
+}
+
+#endif

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

@@ -24,6 +24,8 @@ def makeBranch(self,treeNumpy,isMC):
     def fillBranch(self,treeNumpy,object,isMC):
         if self.mcOnly and not isMC: return
         treeNumpy.fill(self.name, self(object))
+    def __repr__(self):
+        return "<NTupleVariable[%s]>" % self.name
 
 class NTupleObjectType:
     def __init__(self,name,baseObjectTypes=[],mcOnly=[],variables=[]):
@@ -60,6 +62,8 @@ def allBases(self):
         return ret
     def removeVariable(self,name):
         self.variables = [ v for v in self.variables if v.name != name]
+    def __repr__(self):
+        return "<NTupleObjectType[%s]>" % self.name
 
 class NTupleObject:
     def __init__(self, name, objectType, help="", mcOnly=False):
@@ -79,6 +83,8 @@ def fillBranches(self,treeNumpy,object,isMC):
         allvars = self.objectType.allVars(isMC)
         for v in allvars:
             treeNumpy.fill("%s_%s" % (self.name, v.name), v(object))
+    def __repr__(self):
+        return "<NTupleObject[%s]>" % self.name
 
 
 class NTupleCollection:
@@ -137,5 +143,7 @@ def fillBranchesVector(self,treeNumpy,collection,isMC):
         for v in allvars:
             name="%s_%s" % (self.name, v.name) if v.name != "" else self.name
             treeNumpy.vfill(name, [ v(collection[i]) for i in xrange(num) ])
+    def __repr__(self):
+        return "<NTupleCollection[%s]>" % self.name
 
 

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

@@ -64,7 +64,9 @@ def makeMCInfo(self, event):
         verbose = getattr(self.cfg_ana, 'verbose', False)
         rawGenParticles = self.mchandles['genParticles'].product() 
         good = []; keymap = {};
+        allGenParticles = []
         for rawIndex,p in enumerate(rawGenParticles):
+            if self.makeAllGenParticles: allGenParticles.append(p)
             id     = abs(p.pdgId())
             status = p.status()
             # particles must be status > 2, except for prompt leptons, photons, neutralinos
@@ -134,15 +136,15 @@ def makeMCInfo(self, event):
                         print "Error keying %d: motherIndex %d, ancestor.pdgId %d, good[gp.motherIndex].pdgId() %d " % (igp, gp.motherIndex, ancestor.pdgId(),  good[gp.motherIndex].pdgId())
                     break
                 ancestor = None if ancestor.numberOfMothers() == 0 else ancestor.motherRef(0)
-            if abs(gp.pdgId()) not in {1,2,3,4,5,11,12,13,14,15,16,21}:
+            if abs(gp.pdgId()) not in [1,2,3,4,5,11,12,13,14,15,16,21]:
                 gp.sourceId = gp.pdgId()
             if gp.motherIndex != -1:
                 ancestor = good[gp.motherIndex]
                 if ancestor.sourceId != 99 and (ancestor.mass() > gp.mass() or gp.sourceId == 99):
                     gp.sourceId = ancestor.sourceId
         event.generatorSummary = good
         # add the ID of the mother to be able to recreate last decay chains
-        for ip,p in enumerate(event.generatorSummary):
+        for ip,p in enumerate(good):
             moms = realGenMothers(p)
             if len(moms)==0:
                 p.motherId = 0
@@ -164,15 +166,8 @@ def makeMCInfo(self, event):
             print "\n\n"
 
         if self.makeAllGenParticles:
-            event.genParticles = []
-            for rawIndex,p in enumerate(rawGenParticles):
-                if rawIndex in keymap:
-                    gp = event.generatorSummary[keymap[rawIndex]]
-                else:
-                    gp = p
-                    gp.rawIndex = rawIndex 
-                    gp.genSummaryIndex = -1
-                event.genParticles.append(gp)
+            event.genParticles = allGenParticles
+
         if self.makeSplittedGenLists:
             event.genHiggsBosons = []
             event.genVBosons     = []
@@ -225,11 +220,11 @@ def process(self, event):
 setattr(GeneratorAnalyzer,"defaultConfig",
     cfg.Analyzer(GeneratorAnalyzer,
         # BSM particles that can appear with status <= 2 and should be kept
-        stableBSMParticleIds = { 1000022 }, 
+        stableBSMParticleIds = [ 1000022 ], 
         # Particles of which we want to save the pre-FSR momentum (a la status 3).
         # Note that for quarks and gluons the post-FSR doesn't make sense,
         # so those should always be in the list
-        savePreFSRParticleIds = { 1,2,3,4,5, 11,12,13,14,15,16, 21 },
+        savePreFSRParticleIds = [ 1,2,3,4,5, 11,12,13,14,15,16, 21 ],
         # Make also the list of all genParticles, for other analyzers to handle
         makeAllGenParticles = True,
         # Make also the splitted lists

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

@@ -17,6 +17,11 @@
 
 import PhysicsTools.HeppyCore.framework.config as cfg
 
+from ROOT import heppy
+
+
+
+
 def mtw(x1,x2):
     import math
     return math.sqrt(2*x1.pt()*x2.pt()*(1-math.cos(x1.phi()-x2.phi())))
@@ -42,14 +47,15 @@ class IsoTrackAnalyzer( Analyzer ):
 
     def __init__(self, cfg_ana, cfg_comp, looperName ):
         super(IsoTrackAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName)
+        self.IsoTrackIsolationComputer = heppy.IsolationComputer(self.cfg_ana.isoDR)
 
     #----------------------------------------
     # DECLARATION OF HANDLES OF LEPTONS STUFF   
     #----------------------------------------
     def declareHandles(self):
         super(IsoTrackAnalyzer, self).declareHandles()
-        self.handles['cmgCand'] = AutoHandle(self.cfg_ana.candidates,self.cfg_ana.candidatesTypes) 
         self.handles['met'] = AutoHandle( 'slimmedMETs', 'std::vector<pat::MET>' )
+        self.handles['packedCandidates'] = AutoHandle( 'packedPFCandidates', 'std::vector<pat::PackedCandidate>')
 
     def beginLoop(self, setup):
         super(IsoTrackAnalyzer,self).beginLoop(setup)
@@ -70,63 +76,40 @@ def makeIsoTrack(self, event):
         event.selectedIsoCleanTrack = []
         #event.preIsoTrack = []
 
-        pfcands = self.handles['cmgCand'].product()
+        patcands = self.handles['packedCandidates'].product()
 
-        charged = [ p for p in pfcands if ( p.charge() != 0 and abs(p.dz())<=self.cfg_ana.dzMax ) ]
+        charged = [ p for p in patcands if ( p.charge() != 0 and abs(p.dz())<=self.cfg_ana.dzMax ) ]
+
+        self.IsoTrackIsolationComputer.setPackedCandidates(patcands, -1, 0.1, True)
 
         alltrack = map( IsoTrack, charged )
 
 
         for track in alltrack:
 
-            foundNonIsoTrack = False
-
-## ===> require Track Candidate above some pt and charged
             if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.pt() < self.cfg_ana.ptMin) ): continue
             if ( track.pt() < self.cfg_ana.ptMinEMU ): continue
 
+            foundNonIsoTrack = False
 
 ## ===> require is not the leading lepton and opposite to the leading lepton 
             if( (self.cfg_ana.doSecondVeto) and len(event.selectedLeptons)>0) : 
                if( deltaR(event.selectedLeptons[0].eta(), event.selectedLeptons[0].phi(), track.eta(), track.phi()) <0.01) : continue
                if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.charge()*event.selectedLeptons[0].charge()) ): continue
 
+
 ## ===> Redundant:: require the Track Candidate with a  minimum dz
             track.associatedVertex = event.goodVertices[0]
 
 ## ===> compute the isolation and find the most isolated track
 
-            othertracks = [ p for p in charged if( deltaR(p.eta(), p.phi(), track.eta(), track.phi()) < self.cfg_ana.isoDR and p.pt()>self.cfg_ana.ptPartMin ) ]
-            #othertracks = alltrack
-
-            isoSum=0
-            for part in othertracks:
-                #### ===> skip pfcands with a pt min (this should be 0)
-                #if part.pt()<self.cfg_ana.ptPartMin : continue
-                #### ===> skip pfcands outside the cone (this should be 0.3)
-                #if deltaR(part.eta(), part.phi(), track.eta(), track.phi()) > self.cfg_ana.isoDR : continue
-                isoSum += part.pt()
-                ### break the loop to save time
-                if(isoSum > (self.cfg_ana.maxAbsIso + track.pt())):
-                    foundNonIsoTrack = True
-                    break
-
-            if foundNonIsoTrack: continue
-
-               ## reset
-               #isoSum=0
-               #for part in othertracks :
-               #### ===> skip pfcands with a pt min (this should be 0)
-               #    if part.pt()<self.cfg_ana.ptPartMin : continue
-               #### ===> skip pfcands outside the cone (this should be 0.3)
-               #    if deltaR(part.eta(), part.phi(), track.eta(), track.phi()) > self.cfg_ana.isoDR : continue
-               #    isoSum += part.pt()
+            isoSum = self.IsoTrackIsolationComputer.chargedAbsIso(track.physObj, self.cfg_ana.isoDR, 0., self.cfg_ana.ptPartMin)
 
-            #    ###            isoSum = isoSum/track.pt()  ## <--- this is for relIso
+            if(isoSum > (self.cfg_ana.maxAbsIso + track.pt())): continue
 
-            ### ===> the sum should not contain the track candidate
 
-            track.absIso = isoSum - track.pt()
+            #if abs(track.pdgId())==211 :
+            track.absIso = isoSum - track.pt() 
 
             #### store a preIso track
             #event.preIsoTrack.append(track)
@@ -151,6 +134,84 @@ def makeIsoTrack(self, event):
                                 else: 
                                     event.selectedIsoCleanTrack.append(track)
 
+
+
+##        alltrack = map( IsoTrack, charged )
+
+##        for track in alltrack:
+##
+##            foundNonIsoTrack = False
+##
+#### ===> require Track Candidate above some pt and charged
+##            if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.pt() < self.cfg_ana.ptMin) ): continue
+##            if ( track.pt() < self.cfg_ana.ptMinEMU ): continue
+##
+##
+#### ===> require is not the leading lepton and opposite to the leading lepton 
+##            if( (self.cfg_ana.doSecondVeto) and len(event.selectedLeptons)>0) : 
+##               if( deltaR(event.selectedLeptons[0].eta(), event.selectedLeptons[0].phi(), track.eta(), track.phi()) <0.01) : continue
+##               if ( (abs(track.pdgId())!=11) and (abs(track.pdgId())!=13) and (track.charge()*event.selectedLeptons[0].charge()) ): continue
+##
+#### ===> Redundant:: require the Track Candidate with a  minimum dz
+##            track.associatedVertex = event.goodVertices[0]
+##
+#### ===> compute the isolation and find the most isolated track
+##
+##            othertracks = [ p for p in charged if( deltaR(p.eta(), p.phi(), track.eta(), track.phi()) < self.cfg_ana.isoDR and p.pt()>self.cfg_ana.ptPartMin ) ]
+##            #othertracks = alltrack
+##
+##            isoSum=0
+##            for part in othertracks:
+##                #### ===> skip pfcands with a pt min (this should be 0)
+##                #if part.pt()<self.cfg_ana.ptPartMin : continue
+##                #### ===> skip pfcands outside the cone (this should be 0.3)
+##                #if deltaR(part.eta(), part.phi(), track.eta(), track.phi()) > self.cfg_ana.isoDR : continue
+##                isoSum += part.pt()
+##                ### break the loop to save time
+##                if(isoSum > (self.cfg_ana.maxAbsIso + track.pt())):
+##                    foundNonIsoTrack = True
+##                    break
+##
+##            if foundNonIsoTrack: continue
+##
+##               ## reset
+##               #isoSum=0
+##               #for part in othertracks :
+##               #### ===> skip pfcands with a pt min (this should be 0)
+##               #    if part.pt()<self.cfg_ana.ptPartMin : continue
+##               #### ===> skip pfcands outside the cone (this should be 0.3)
+##               #    if deltaR(part.eta(), part.phi(), track.eta(), track.phi()) > self.cfg_ana.isoDR : continue
+##               #    isoSum += part.pt()
+##
+##            #    ###            isoSum = isoSum/track.pt()  ## <--- this is for relIso
+##
+##            ### ===> the sum should not contain the track candidate
+##
+##            track.absIso = isoSum - track.pt()
+##
+##            #### store a preIso track
+##            #event.preIsoTrack.append(track)
+##            
+###            if (isoSum < minIsoSum ) :
+##            if(track.absIso < min(0.2*track.pt(), self.cfg_ana.maxAbsIso)): 
+##                event.selectedIsoTrack.append(track)
+##
+##                if self.cfg_ana.doPrune:
+##                    myMet = self.handles['met'].product()[0]
+##                    mtwIsoTrack = mtw(track, myMet)
+##                    if mtwIsoTrack < 100:
+##                        if abs(track.pdgId()) == 11 or abs(track.pdgId()) == 13:
+##                            if track.pt()>5 and track.absIso/track.pt()<0.2:
+##
+##                                myLeptons = [ l for l in event.selectedLeptons if l.pt() > 10 ] 
+##                                nearestSelectedLeptons = makeNearestLeptons(myLeptons,track, event)
+##                                if len(nearestSelectedLeptons) > 0:
+##                                    for lep in nearestSelectedLeptons:
+##                                        if deltaR(lep.eta(), lep.phi(), track.eta(), track.phi()) > 0.1:
+##                                            event.selectedIsoCleanTrack.append(track)
+##                                else: 
+##                                    event.selectedIsoCleanTrack.append(track)
+
         event.selectedIsoTrack.sort(key = lambda l : l.pt(), reverse = True)
         event.selectedIsoCleanTrack.sort(key = lambda l : l.pt(), reverse = True)