New version of heppy

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/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/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/core/AutoFillTreeProducer.py

@@ -82,7 +82,7 @@ def declareVariables(self,setup):
         tree = self.tree
         self.declareCoreVariables(tree, isMC)
 
-        if not hasattr(self.cfg_ana,"ignoreAnalyzerBookings") or not self.cfg_ana.ignoreAnalyzerBooking :
+        if not hasattr(self.cfg_ana,"ignoreAnalyzerBookings") or not self.cfg_ana.ignoreAnalyzerBookings :
             #import variables declared by the analyzers
             if hasattr(setup,"globalVariables"):
                 self.globalVariables+=setup.globalVariables
@@ -137,6 +137,9 @@ def fillCoreVariables(self, tr, event, isMC):
                     tr.vfill('pdfWeight_%s' % pdf, event.pdfWeights[pdf])
 
     def process(self, event):
+	if hasattr(self.cfg_ana,"filter") :	
+		if not self.cfg_ana.filter(event) :
+			return True #do not stop processing, just filter myself
         self.readCollections( event.input)
         self.fillTree(event)
 
@@ -167,4 +170,33 @@ def fillTree(self, event, resetFirst=True):
                 c.fillBranchesVector(self.tree, getattr(event, cn), isMC)
 
         self.tree.tree.Fill()      
+
+    def getPythonWrapper(self):
+        """
+        This function produces a string that contains a Python wrapper for the event.
+        The wrapper is automatically generated based on the collections and allows the full
+        event contents to be accessed from subsequent Analyzers using e.g.
+
+        leps = event.selLeptons #is of type selLeptons
+        pt0 = leps[0].pt
+
+        One just needs to add the EventAnalyzer to the sequence.
+        """
+
+        isMC = self.cfg_comp.isMC 
+
+        classes = ""
+        anclass = ""
+        anclass += "from PhysicsTools.HeppyCore.framework.analyzer import Analyzer\n"
+        anclass += "class EventAnalyzer(Analyzer):\n"
+        anclass += "    def __init__(self, cfg_ana, cfg_comp, looperName):\n"
+        anclass += "        super(EventAnalyzer, self).__init__(cfg_ana, cfg_comp, looperName)\n"
+
+        anclass += "    def process(self, event):\n"
+
+        for cname, coll in self.collections.items():
+            classes += coll.get_py_wrapper_class(isMC)
+            anclass += "        event.{0} = {0}.make_array(event)\n".format(coll.name)
+
+        return classes + "\n" + anclass
 

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

@@ -4,6 +4,8 @@
 from PhysicsTools.Heppy.analyzers.core.AutoHandle import AutoHandle
 from PhysicsTools.HeppyCore.statistics.average import Average
 from PhysicsTools.Heppy.physicsutils.PileUpSummaryInfo import PileUpSummaryInfo
+import PhysicsTools.HeppyCore.framework.config as cfg
+
 from ROOT import TFile, TH1F
 
 class PileUpAnalyzer( Analyzer ):
@@ -57,7 +59,7 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
           self.cfg_comp.puFileData = None
 
         if self.cfg_comp.isMC or self.cfg_comp.isEmbed:
-            if self.cfg_comp.puFileMC is None and self.cfg_comp.puFileData is None:
+            if not hasattr(self.cfg_comp,"puFileMC") or (self.cfg_comp.puFileMC is None and self.cfg_comp.puFileData is None):
                 self.enable = False
             else:
                 assert( os.path.isfile(self.cfg_comp.puFileMC) )
@@ -145,3 +147,12 @@ def write(self, setup):
         super(PileUpAnalyzer, self).write(setup)
         if self.cfg_comp.isMC and self.doHists:
             self.rawmcpileup.write()
+
+
+setattr(PileUpAnalyzer,"defaultConfig", cfg.Analyzer(
+    class_object = PileUpAnalyzer,
+    true = True,  # use number of true interactions for reweighting
+    makeHists=False
+)
+)
+

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

@@ -11,14 +11,12 @@ class TreeAnalyzerNumpy( Analyzer ):
 
     def __init__(self, cfg_ana, cfg_comp, looperName):
         super(TreeAnalyzerNumpy,self).__init__(cfg_ana, cfg_comp, looperName)
-        self.outservicename = "outputfile"
-        if hasattr(cfg_ana,"outservicename") :
-            self.outservicename = cfg_ana.outservicename
+        self.outservicename = getattr(cfg_ana,"outservicename","outputfile")
+        self.treename = getattr(cfg_ana,"treename","tree")
 
 
     def beginLoop(self, setup) :
         super(TreeAnalyzerNumpy, self).beginLoop(setup)
-        print setup.services
         if self.outservicename in setup.services:
             print "Using outputfile given in", self.outservicename
             self.file = setup.services[self.outservicename].file
@@ -28,7 +26,10 @@ def beginLoop(self, setup) :
             isCompressed = self.cfg_ana.isCompressed if hasattr(self.cfg_ana,'isCompressed') else 1
             print 'Compression', isCompressed
             self.file = TFile( fileName, 'recreate', '', isCompressed )
-        self.tree = Tree('tree', self.name)
+        self.file.cd()
+        if self.file.Get(self.treename) :
+            raise RuntimeError, "You are booking two Trees with the same name in the same file"
+        self.tree = Tree(self.treename, self.name)
         self.tree.setDefaultFloatType(getattr(self.cfg_ana, 'defaultFloatType','D')); # or 'F'
         self.declareVariables(setup)
 

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,35 @@ 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
+
+    def get_cpp_declaration(self, isMC):
+        s = []
+        for v in self.objectType.allVars(isMC):
+            s += ["{0} {1}__{2}[{3}];".format(v.type.__name__, self.name, v.name, self.maxlen)]
+        return "\n".join(s)
+
+    def get_cpp_wrapper_class(self, isMC):
+        s = "class %s {\n" % self.name
+        s += "public:\n"
+        for v in self.objectType.allVars(isMC):
+            s += "    {0} {1};\n".format(v.type.__name__, v.name)
+        s += "};\n"
+        return s
+
+    def get_py_wrapper_class(self, isMC):
+        s = "class %s:\n" % self.name
+        s += "    def __init__(self, tree, n):\n"
+        for v in self.objectType.allVars(isMC):
+            if len(v.name)>0:
+                s += "        self.{0} = tree.{1}_{2}[n];\n".format(v.name, self.name, v.name)
+            else:
+                s += "        self.{0} = tree.{0}[n];\n".format(self.name)
+
+        s += "    @staticmethod\n"
+        s += "    def make_array(event):\n"
+        s += "        return [{0}(event.input, i) for i in range(event.input.n{0})]\n".format(self.name)
+        return s