diff --git a/CMGTools/RootTools/interface/ReclusterJets.h b/CMGTools/RootTools/interface/ReclusterJets.h new file mode 100644 index 0000000000000..e63f8f03bd250 --- /dev/null +++ b/CMGTools/RootTools/interface/ReclusterJets.h @@ -0,0 +1,52 @@ +#ifndef ReclusterJets_h +#define ReclusterJets_h + +#include +#include +#include +#include +#include +#include "DataFormats/Math/interface/LorentzVector.h" + +#include +#include +#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 & objects, double ktpower, double rparam); + + /// get grouping (inclusive jets) + std::vector getGrouping(double ptMin = 0.0); + + /// get grouping (exclusive jets, until n are left) + std::vector getGroupingExclusive(int njets); + /// get grouping (exclusive jets, up to cut dcut) + std::vector getGroupingExclusive(double dcut); + + private: + // pack the returns in a fwlite-friendly way + std::vector makeP4s(const std::vector &jets) ; + + // used to handle the inputs + std::vector fjInputs_; // fastjet inputs + + double ktpower_; + double rparam_; + + /// fastjet outputs + typedef boost::shared_ptr ClusterSequencePtr; + ClusterSequencePtr fjClusterSeq_; + std::vector inclusiveJets_; + std::vector exclusiveJets_; +}; + +#endif diff --git a/CMGTools/RootTools/src/ReclusterJets.cc b/CMGTools/RootTools/src/ReclusterJets.cc new file mode 100644 index 0000000000000..d47835d8f3a09 --- /dev/null +++ b/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 & 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 ReclusterJets::makeP4s(const std::vector &jets) { + std::vector JetObjectsAll; + for (const fastjet::PseudoJet & pj : jets) { + JetObjectsAll.push_back( LorentzVector( pj.px(), pj.py(), pj.pz(), pj.e() ) ); + } + return JetObjectsAll; +} +std::vector ReclusterJets::getGrouping(double ptMin) { + // recluster jet + inclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(ptMin)); + // return + return makeP4s(inclusiveJets_); +} + +std::vector ReclusterJets::getGroupingExclusive(double dcut) { + // recluster jet + exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(dcut)); + // return + return makeP4s(exclusiveJets_); +} + +std::vector ReclusterJets::getGroupingExclusive(int njets) { + // recluster jet + exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(njets)); + // return + return makeP4s(exclusiveJets_); +} diff --git a/CMGTools/RootTools/src/classes.h b/CMGTools/RootTools/src/classes.h index 917e88402ebc2..27cd589f7fc1b 100644 --- a/CMGTools/RootTools/src/classes.h +++ b/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 px, vector py, vector pz, vector E, int hemi_seed, int hemi_association); HemisphereViaKt hemisphere(vector px, vector py, vector pz, vector E, double ktpower); + ReclusterJets reclusterJets(vector px, vector py, vector pz, vector E, double ktpower, double rparam); Davismt2 mt2; mt2w_bisect::mt2w mt2wlept; AlphaT alphaT; diff --git a/CMGTools/RootTools/src/classes_def.xml b/CMGTools/RootTools/src/classes_def.xml index 3ed76fd895a15..6a6b1219ea2a7 100644 --- a/CMGTools/RootTools/src/classes_def.xml +++ b/CMGTools/RootTools/src/classes_def.xml @@ -6,6 +6,7 @@ + diff --git a/CMGTools/TTHAnalysis/cfg/run_susySinglelepton_cfg.py b/CMGTools/TTHAnalysis/cfg/run_susySinglelepton_cfg.py index 0fb4cb7c3c58b..562e779431cc8 100644 --- a/CMGTools/TTHAnalysis/cfg/run_susySinglelepton_cfg.py +++ b/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, ]) diff --git a/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusySingleLepton.py b/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusySingleLepton.py index e1dc45640b319..468d7be269f97 100644 --- a/CMGTools/TTHAnalysis/python/analyzers/treeProducerSusySingleLepton.py +++ b/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 diff --git a/CMGTools/TTHAnalysis/python/analyzers/ttHReclusterJetsAnalyzer.py b/CMGTools/TTHAnalysis/python/analyzers/ttHReclusterJetsAnalyzer.py new file mode 100644 index 0000000000000..a39224a49c8b2 --- /dev/null +++ b/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') + + 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