Merge pull request #40 from cbernet/heppy_7_2_2_patch2

Sync with central Heppy repository

PhysicsTools/Heppy/interface/Apc.h

@@ -0,0 +1,28 @@
+#ifndef PhysicsTools_Heppy_Apc_h
+#define PhysicsTools_Heppy_Apc_h
+
+#include <cmath>
+#include <numeric>
+#include <vector>
+
+
+/*  Apc
+ *
+ *  Calculates the Apc event variable for a given input jet collection
+ *
+ */
+
+namespace heppy {
+
+struct Apc {
+
+  static double getApcJetMetMin( const std::vector<double>& et,
+                                 const std::vector<double>& px,
+                                 const std::vector<double>& py,
+                                 const double metx, const double mety);
+
+};
+
+};
+
+#endif // Apc_h

PhysicsTools/Heppy/interface/ReclusterJets.h

@@ -29,13 +29,27 @@ class ReclusterJets {
 
   /// get grouping (exclusive jets, until n are left)
   std::vector<LorentzVector> getGroupingExclusive(int njets);
+
   /// get grouping (exclusive jets, up to cut dcut)
   std::vector<LorentzVector> getGroupingExclusive(double dcut);
 
+  /// get pruned 4-vector
+  LorentzVector getPruned(double zcut, double rcutFactor) ;
+
+  /// get pruned 4-vector for a given subject (must be called after getGroupingExclusive)
+  LorentzVector getPrunedSubjetExclusive(unsigned int isubjet, double zcut, double rcutFactor) ;
+
+  /// get pruned 4-vector for a given subject (must be called after getGroupingInclusive)
+  LorentzVector getPrunedSubjetInclusive(unsigned int isubjet, double zcut, double rcutFactor) ;
+
+
  private:
   // pack the returns in a fwlite-friendly way
   std::vector<LorentzVector> makeP4s(const std::vector<fastjet::PseudoJet> &jets) ;
 
+  // prune and return in fa fwlite-friendly way
+  LorentzVector getPruned(const fastjet::PseudoJet & jet, double zcut, double rcutFactor) ;
+
   // used to handle the inputs
   std::vector<fastjet::PseudoJet> fjInputs_;        // fastjet inputs
 

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

@@ -88,7 +88,7 @@ def makeMCInfo(self, event):
                     continue
             else:
                 # everything else, we want it after radiation, i.e. just before decay
-                if any((p.daughter(j).pdgId() == p.pdgId() and p.daughter(j).status > 2) for j in xrange(p.numberOfDaughters())):
+                if any((p.daughter(j).pdgId() == p.pdgId() and p.daughter(j).status() > 2) for j in xrange(p.numberOfDaughters())):
                     #print "    fail auto-decay"
                     continue
             # FIXME find a better criterion to discard there
@@ -112,6 +112,9 @@ def makeMCInfo(self, event):
                     if not any(mom.daughter(j2).pdgId() == mom.pdgId() for j2 in xrange(mom.numberOfDaughters())):
                         #print "         pass no-self-decay"
                         ok = True
+                    # Account for generator feature with Higgs decaying to itself with same four-vector but no daughters
+                    elif mom.pdgId() == 25 and any(mom.daughter(j2).pdgId() == 25 and mom.daughter(j2).numberOfDaughters()==0 for j2 in range(mom.numberOfDaughters())):
+                        ok = True
                 if abs(mom.pdgId()) == 15:
                     # if we're a tau daughter we're status 2
                     # if we passed all the previous steps, then we're a prompt lepton

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

@@ -99,7 +99,7 @@ def makeIsoTrack(self, event):
 
 
 ## ===> Redundant:: require the Track Candidate with a  minimum dz
-            track.associatedVertex = event.goodVertices[0]
+            track.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
 
 ## ===> compute the isolation and find the most isolated track
 

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

@@ -234,7 +234,7 @@ def makeAllMuons(self, event):
           else: raise RuntimeError,  "Unsupported value for mu_effectiveAreas: can only use Data2012 (rho: ?) and Phys14_v1 (rho: fixedGridRhoFastjetAll)"
         # Attach the vertex to them, for dxy/dz calculation
         for mu in allmuons:
-            mu.associatedVertex = event.goodVertices[0]
+            mu.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
 
         # Compute relIso in 0.3 and 0.4 cones
         for mu in allmuons:
@@ -308,7 +308,7 @@ def makeAllElectrons(self, event):
 
         # Attach the vertex
         for ele in allelectrons:
-            ele.associatedVertex = event.goodVertices[0]
+            ele.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
 
         # Compute relIso with R=0.3 and R=0.4 cones
         for ele in allelectrons:

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

@@ -24,14 +24,7 @@ def declareHandles(self):
         super(METAnalyzer, self).declareHandles()
         self.handles['met'] = AutoHandle( 'slimmedMETs', 'std::vector<pat::MET>' )
         self.handles['nopumet'] = AutoHandle( 'slimmedMETs', 'std::vector<pat::MET>' )
-        self.handles['cmgCand1'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand2'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand3'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand4'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand5'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand6'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand7'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
-        self.handles['cmgCand8'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
+        self.handles['cmgCand'] = AutoHandle( self.cfg_ana.candidates, self.cfg_ana.candidatesTypes )
         self.handles['vertices'] =  AutoHandle( "offlineSlimmedPrimaryVertices", 'std::vector<reco::Vertex>', fallbackLabel="offlinePrimaryVertices" )
 
     def beginLoop(self, setup):
@@ -44,18 +37,35 @@ def makeTkMETs(self, event):
         event.tkMet = 0  
 
         charged = []
+        chargedchs = []
+        chargedPVLoose = []
+        chargedPVTight = []
+
         pfcands = self.handles['cmgCand'].product()
 
         for i in xrange(pfcands.size()):
 
 ## ===> require the Track Candidate charge and with a  minimum dz 
 
-            if (pfcands.at(i).charge()!=0 and (abs(pfcands.at(i).dz())<=self.cfg_ana.dzMax )):
-
-                charged.append(pfcands.at(i))
-
+            if (pfcands.at(i).charge()!=0):
+
+                if abs(pfcands.at(i).dz())<=self.cfg_ana.dzMax:
+                    charged.append(pfcands.at(i))
+
+                if pfcands.at(i).fromPV()>0:
+                    chargedchs.append(pfcands.at(i))
+
+                if pfcands.at(i).fromPV()>1:
+                    chargedPVLoose.append(pfcands.at(i))
+
+                if pfcands.at(i).fromPV()>2:
+                    chargedPVTight.append(pfcands.at(i))
+
         import ROOT
-        event.tkMet = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in charged])) , -1.*(sum([x.py() for x in charged])), 0, 0 )
+        event.tkMet = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in charged])) , -1.*(sum([x.py() for x in charged])), 0, math.hypot((sum([x.px() for x in charged])),(sum([x.py() for x in charged]))) )
+        event.tkMetchs = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in chargedchs])) , -1.*(sum([x.py() for x in chargedchs])), 0, math.hypot((sum([x.px() for x in chargedchs])),(sum([x.py() for x in chargedchs]))) )
+        event.tkMetPVLoose = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in chargedPVLoose])) , -1.*(sum([x.py() for x in chargedPVLoose])), 0, math.hypot((sum([x.px() for x in chargedPVLoose])),(sum([x.py() for x in chargedPVLoose]))) )
+        event.tkMetPVTight = ROOT.reco.Particle.LorentzVector(-1.*(sum([x.px() for x in chargedPVTight])) , -1.*(sum([x.py() for x in chargedPVTight])), 0, math.hypot((sum([x.px() for x in chargedPVTight])),(sum([x.py() for x in chargedPVTight]))) )
 ##        print 'tkmet',event.tkMet.pt(),'tkmetphi',event.tkMet.phi()
 
 

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

@@ -105,7 +105,7 @@ def matchPhotons(self, event):
         event.genPhotons = [ x for x in event.genParticles if x.status() == 1 and abs(x.pdgId()) == 22 ]
         event.genPhotonsWithMom = [ x for x in event.genPhotons if x.numberOfMothers()>0 ]
         event.genPhotonsWithoutMom = [ x for x in event.genPhotons if x.numberOfMothers()==0 ]
-        event.genPhotonsMatched = [ x for x in event.genPhotonsWithMom if abs(x.mother(0).pdgId())<23 ]
+        event.genPhotonsMatched = [ x for x in event.genPhotonsWithMom if abs(x.mother(0).pdgId())<23 or x.mother(0).pdgId()==2212 ]
         match = matchObjectCollection3(event.allphotons, event.genPhotonsMatched, deltaRMax = 0.1)
         matchNoMom = matchObjectCollection3(event.allphotons, event.genPhotonsWithoutMom, deltaRMax = 0.1)
         packedGenParts = [ p for p in self.mchandles['packedGen'].product() if abs(p.eta()) < 3.1 ]

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

@@ -44,7 +44,7 @@ def makeTaus(self, event):
 
         foundTau = False
         for tau in alltaus:
-            tau.associatedVertex = event.goodVertices[0]
+            tau.associatedVertex = event.goodVertices[0] if len(event.goodVertices)>0 else event.vertices[0]
             tau.lepVeto = False
             tau.idDecayMode = tau.tauID("decayModeFinding")
             tau.idDecayModeNewDMs = tau.tauID("decayModeFindingNewDMs")

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

@@ -94,6 +94,7 @@
     NTupleVariable("dPhiScTrkIn",    lambda x : x.deltaPhiSuperClusterTrackAtVtx() if abs(x.pdgId())==11 else 0, help="Electron deltaPhiSuperClusterTrackAtVtx (without absolute value!)"),
     NTupleVariable("hadronicOverEm", lambda x : x.hadronicOverEm() if abs(x.pdgId())==11 else 0, help="Electron hadronicOverEm"),
     NTupleVariable("eInvMinusPInv",  lambda x : ((1.0/x.ecalEnergy() - x.eSuperClusterOverP()/x.ecalEnergy()) if x.ecalEnergy()>0. else 9e9) if abs(x.pdgId())==11 else 0, help="Electron 1/E - 1/p  (without absolute value!)"),
+    NTupleVariable("etaSc", lambda x : x.superCluster().eta() if abs(x.pdgId())==11 else -100, help="Electron supercluster pseudorapidity"),
 ])
 
 

PhysicsTools/Heppy/python/physicsobjects/Electron.py

@@ -147,24 +147,18 @@ def photonIsoR(self,R=0.4):
         elif R == 0.4: return self.physObj.photonIso()
         raise RuntimeError, "Electron photonIso missing for R=%s" % R
 
-    def chargedAllIso(self,R=0.4):
+    def chargedAllIsoR(self,R=0.4):
         if   R == 0.3: return self.physObj.pfIsolationVariables().sumChargedParticlePt 
         raise RuntimeError, "Electron chargedAllIso missing for R=%s" % R
 
+    def chargedAllIso(self):
+        raise RuntimeError, "Electron chargedAllIso missing"
+
     def puChargedHadronIsoR(self,R=0.4):
         if   R == 0.3: return self.physObj.pfIsolationVariables().sumPUPt 
         elif R == 0.4: return self.physObj.puChargedHadronIso()
         raise RuntimeError, "Electron chargedHadronIso missing for R=%s" % R
 
-
-
-
-    def chargedAllIso(self):
-        '''This function is used in the isolation, see Lepton class.
-        Here, we replace the all charged isolation by the all charged isolation with cone veto'''
-        return self.chargedAllIsoWithConeVeto()
-
-
     def dxy(self, vertex=None):
         '''Returns dxy.
         Computed using vertex (or self.associatedVertex if vertex not specified),