Merge pull request cms-sw#149 from gpetruc/pr143_withVectors

Merging.

You can port in the hemisphere code in a separate PR.

CMGTools/RootTools/interface/ReclusterJets.h

@@ -0,0 +1,52 @@
+#ifndef ReclusterJets_h
+#define ReclusterJets_h
+
+#include <vector>
+#include <iostream>
+#include <cmath>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include "DataFormats/Math/interface/LorentzVector.h"
+
+#include <boost/shared_ptr.hpp>
+#include <fastjet/internal/base.hh>
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/Selector.hh"
+#include "fastjet/PseudoJet.hh"
+
+
+class ReclusterJets {
+
+ public:
+  typedef math::XYZTLorentzVector LorentzVector;
+
+  ReclusterJets(const std::vector<LorentzVector> & objects, double ktpower, double rparam);
+
+  /// get grouping (inclusive jets)
+  std::vector<LorentzVector> getGrouping(double ptMin = 0.0);
+
+  /// 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);
+
+ private:
+  // pack the returns in a fwlite-friendly way
+  std::vector<LorentzVector> makeP4s(const std::vector<fastjet::PseudoJet> &jets) ;
+
+  // used to handle the inputs
+  std::vector<fastjet::PseudoJet> fjInputs_;        // fastjet inputs
+
+  double ktpower_;
+  double rparam_;
+
+  /// fastjet outputs
+  typedef boost::shared_ptr<fastjet::ClusterSequence>  ClusterSequencePtr;
+  ClusterSequencePtr fjClusterSeq_;    
+  std::vector<fastjet::PseudoJet> inclusiveJets_; 
+  std::vector<fastjet::PseudoJet> exclusiveJets_; 
+};
+
+#endif    

CMGTools/RootTools/src/ReclusterJets.cc

@@ -0,0 +1,69 @@
+#include "CMGTools/RootTools/interface/ReclusterJets.h"
+#include "FWCore/Utilities/interface/Exception.h"
+
+using namespace std;
+
+//using namespace std;
+using namespace fastjet;
+
+
+ReclusterJets::ReclusterJets(const std::vector<LorentzVector> & objects, double ktpower, double rparam) : 
+    ktpower_(ktpower), rparam_(rparam) 
+{
+  // define jet inputs
+  fjInputs_.clear();
+  int index=0;
+  for (const LorentzVector &o : objects) {
+    fastjet::PseudoJet j(o.Px(),o.Py(),o.Pz(),o.E());
+    j.set_user_index(index); index++; // in case we want to know which piece ended where
+    fjInputs_.push_back(j);
+  }
+
+  // choose a jet definition
+  fastjet::JetDefinition jet_def;
+
+  // prepare jet def 
+  if (ktpower_ == 1.0) {
+    jet_def = JetDefinition(kt_algorithm, rparam_);
+  }  else if (ktpower_ == 0.0) {
+    jet_def = JetDefinition(cambridge_algorithm, rparam_);
+  }  else if (ktpower_ == -1.0) {
+    jet_def = JetDefinition(antikt_algorithm, rparam_);
+  }  else {
+    throw cms::Exception("InvalidArgument", "Unsupported ktpower value");
+  }
+
+  // print out some infos
+  //  cout << "Clustering with " << jet_def.description() << endl;
+  ///
+  // define jet clustering sequence
+  fjClusterSeq_ = ClusterSequencePtr( new fastjet::ClusterSequence( fjInputs_, jet_def)); 
+}
+
+std::vector<math::XYZTLorentzVector> ReclusterJets::makeP4s(const std::vector<fastjet::PseudoJet> &jets) {
+  std::vector<math::XYZTLorentzVector> JetObjectsAll;
+  for (const fastjet::PseudoJet & pj : jets) {
+    JetObjectsAll.push_back( LorentzVector( pj.px(), pj.py(), pj.pz(), pj.e() ) );
+  }
+  return JetObjectsAll;
+}
+std::vector<math::XYZTLorentzVector> ReclusterJets::getGrouping(double ptMin) {
+  // recluster jet
+  inclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(ptMin));
+  // return
+  return makeP4s(inclusiveJets_);
+}
+
+std::vector<math::XYZTLorentzVector> ReclusterJets::getGroupingExclusive(double dcut) {
+  // recluster jet
+  exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(dcut));
+  // return
+  return makeP4s(exclusiveJets_);
+}
+
+std::vector<math::XYZTLorentzVector> ReclusterJets::getGroupingExclusive(int njets) {
+  // recluster jet
+  exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(njets));
+  // return
+  return makeP4s(exclusiveJets_);
+}

CMGTools/RootTools/src/classes.h

@@ -12,6 +12,7 @@
 #include "CMGTools/RootTools/interface/AlphaT.h"
 #include "CMGTools/RootTools/interface/HemisphereViaKt.h"
 #include "CMGTools/RootTools/interface/RecoilCorrector.h"
+#include "CMGTools/RootTools/interface/ReclusterJets.h"
 
 namespace {
   namespace {
@@ -24,6 +25,7 @@ namespace {
     EGammaMvaEleEstimatorFWLite egMVA;
     Hemisphere hemisphere(vector<float> px, vector<float> py, vector<float> pz, vector<float> E, int hemi_seed, int hemi_association);
     HemisphereViaKt hemisphere(vector<float> px, vector<float> py, vector<float> pz, vector<float> E, double ktpower);
+    ReclusterJets reclusterJets(vector<float> px, vector<float> py, vector<float> pz, vector<float> E, double ktpower, double rparam);
     Davismt2 mt2;
     mt2w_bisect::mt2w mt2wlept;
     AlphaT alphaT;

CMGTools/RootTools/src/classes_def.xml

@@ -6,6 +6,7 @@
   <class name="Davismt2" transient="true" />
   <class name="Hemisphere"/>
   <class name="HemisphereViaKt"/>
+  <class name="ReclusterJets"/>
   <class name="mt2w_bisect::mt2w"  transient="true" />
   <class name="MuScleFitCorrector" transient="true" />
   <class name="TriggerBitChecker" transient="true" />

CMGTools/TTHAnalysis/cfg/run_susySinglelepton_cfg.py

@@ -29,6 +29,10 @@
 ttHJetAna.minLepPt = 10 
 ttHJetMCAna.smearJets     = False # do we need to smear the jets?
 
+ttHReclusterJets = cfg.Analyzer(
+            'ttHReclusterJetsAnalyzer',
+            )
+
 # Event Analyzer for susy multi-lepton (at the moment, it's the TTH one)
 ttHEventAna = cfg.Analyzer(
     'ttHLepEventAnalyzer',
@@ -77,6 +81,7 @@
 sequence = cfg.Sequence(susyCoreSequence+[
     ttHIsoTrackAna,
     ttHEventAna,
+    ttHReclusterJets,
     treeProducer,
     ])
 

CMGTools/TTHAnalysis/python/analyzers/treeProducerSusySingleLepton.py

@@ -42,6 +42,7 @@ def __init__(self, cfg_ana, cfg_comp, looperName):
             "selectedIsoTrack"    : NTupleCollection("isoTrack", isoTrackType, 50, help="isoTrack, sorted by pt"),
             ##------------------------------------------------
             "cleanJetsAll"       : NTupleCollection("Jet",     jetTypeSusy, 25, help="Cental jets after full selection and cleaning, sorted by pt"),
+            "reclusteredFatJets" : NTupleCollection("FatJet",     fourVectorType,20, help="FatJets reclusterd from ak4 cleanJetsAll"),
         })
 
         ## Book the variables, but only if we're called explicitly and not through a base class

CMGTools/TTHAnalysis/python/analyzers/ttHReclusterJetsAnalyzer.py

@@ -0,0 +1,45 @@
+from CMGTools.RootTools.fwlite.Analyzer import Analyzer
+from CMGTools.RootTools.fwlite.AutoHandle import AutoHandle
+
+import ROOT
+
+class ttHReclusterJetsAnalyzer( Analyzer ):
+    def __init__(self, cfg_ana, cfg_comp, looperName ):
+        super(ttHReclusterJetsAnalyzer,self).__init__(cfg_ana,cfg_comp,looperName) 
+
+    def declareHandles(self):
+        super(ttHReclusterJetsAnalyzer, self).declareHandles()
+       #genJets                                                                                                                                                                     
+        self.handles['genJets'] = AutoHandle( 'slimmedGenJets','std::vector<reco::GenJet>')
+
+    def beginLoop(self):
+        super(ttHReclusterJetsAnalyzer,self).beginLoop()
+        self.counters.addCounter('pairs')
+        count = self.counters.counter('pairs')
+        count.register('all events')
+
+    def makeFatJets(self, event):
+        objects40jc = [ j for j in event.cleanJets if j.pt() > 40.0 and abs(j.eta())<2.4 ]
+
+        if len(objects40jc)>=1:
+
+           objects  = ROOT.std.vector(ROOT.reco.Particle.LorentzVector)()
+           for jet in objects40jc:
+                objects.push_back(jet.p4())
+
+           reclusterJets = ROOT.ReclusterJets(objects, 1.,1.2)
+           inclusiveJets = reclusterJets.getGrouping()
+
+           # maybe should dress them as Jet objects in the future
+           # for the moment, we just make LorentzVector
+           event.reclusteredFatJets = [ ROOT.reco.Particle.LorentzVector(p4) for p4 in inclusiveJets ]
+           # note 1: just taking inclusiveJets is not ok, since it's not a python list but a std::vector
+           # note 2: [p4 for p4 in inclusiveJets] is also bad, since these are references to values inside a temporary std::vector
+
+    def process(self, iEvent, event):
+        self.readCollections( iEvent )
+
+        event.reclusteredFatJets = []
+        self.makeFatJets(event)
+
+        return True