Heppy/src/ReclusterJets.cc

0001 #include <memory>
0002
0003 #include "PhysicsTools/Heppy/interface/ReclusterJets.h"
0004 #include "FWCore/Utilities/interface/Exception.h"
0005 #include "fastjet/tools/Pruner.hh"
0006
0007 using namespace std;
0008
0009 //using namespace std;
0010 using namespace fastjet;
0011
0012 namespace heppy {
0013
0014   ReclusterJets::ReclusterJets(const std::vector<LorentzVector> &objects, double ktpower, double rparam)
0015       : ktpower_(ktpower), rparam_(rparam) {
0016     // define jet inputs
0017     fjInputs_.clear();
0018     int index = 0;
0019     for (const LorentzVector &o : objects) {
0020       fastjet::PseudoJet j(o.Px(), o.Py(), o.Pz(), o.E());
0021       j.set_user_index(index);
0022       index++;  // in case we want to know which piece ended where
0023       fjInputs_.push_back(j);
0024     }
0025
0026     // choose a jet definition
0027     fastjet::JetDefinition jet_def;
0028
0029     // prepare jet def
0030     if (ktpower_ == 1.0) {
0031       jet_def = JetDefinition(kt_algorithm, rparam_);
0032     } else if (ktpower_ == 0.0) {
0033       jet_def = JetDefinition(cambridge_algorithm, rparam_);
0034     } else if (ktpower_ == -1.0) {
0035       jet_def = JetDefinition(antikt_algorithm, rparam_);
0036     } else {
0037       throw cms::Exception("InvalidArgument", "Unsupported ktpower value");
0038     }
0039
0040     // print out some infos
0041     //  cout << "Clustering with " << jet_def.description() << endl;
0042     ///
0043     // define jet clustering sequence
0044     fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, jet_def);
0045   }
0046
0047   std::vector<math::XYZTLorentzVector> ReclusterJets::makeP4s(const std::vector<fastjet::PseudoJet> &jets) {
0048     std::vector<math::XYZTLorentzVector> JetObjectsAll;
0049     JetObjectsAll.reserve(jets.size());
0050     for (const fastjet::PseudoJet &pj : jets) {
0051       JetObjectsAll.push_back(LorentzVector(pj.px(), pj.py(), pj.pz(), pj.e()));
0052     }
0053     return JetObjectsAll;
0054   }
0055   std::vector<math::XYZTLorentzVector> ReclusterJets::getGrouping(double ptMin) {
0056     // recluster jet
0057     inclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(ptMin));
0058     // return
0059     return makeP4s(inclusiveJets_);
0060   }
0061
0062   std::vector<math::XYZTLorentzVector> ReclusterJets::getGroupingExclusive(double dcut) {
0063     // recluster jet
0064     exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(dcut));
0065     // return
0066     return makeP4s(exclusiveJets_);
0067   }
0068
0069   std::vector<math::XYZTLorentzVector> ReclusterJets::getGroupingExclusive(int njets) {
0070     // recluster jet
0071     exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(njets));
0072     // return
0073     return makeP4s(exclusiveJets_);
0074   }
0075
0076   math::XYZTLorentzVector ReclusterJets::getPruned(double zcut, double rcutFactor) {
0077     // cluster everything first
0078     exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(1));
0079     // get pruned exclusive
0080     return getPrunedSubjetExclusive(0, zcut, rcutFactor);
0081   }
0082
0083   math::XYZTLorentzVector ReclusterJets::getPrunedSubjetExclusive(unsigned int isubjet,
0084                                                                   double zcut,
0085                                                                   double rcutFactor) {
0086     if (isubjet >= exclusiveJets_.size()) {
0087       throw cms::Exception("InvalidArgument", "getPrunedSubjetExclusive called for non-existing exclusive subjet");
0088     }
0089     return getPruned(exclusiveJets_[isubjet], zcut, rcutFactor);
0090   }
0091   math::XYZTLorentzVector ReclusterJets::getPrunedSubjetInclusive(unsigned int isubjet,
0092                                                                   double zcut,
0093                                                                   double rcutFactor) {
0094     if (isubjet >= inclusiveJets_.size()) {
0095       throw cms::Exception("InvalidArgument", "getPrunedSubjetInclusive called for non-existing inclusive subjet");
0096     }
0097     return getPruned(inclusiveJets_[isubjet], zcut, rcutFactor);
0098   }
0099
0100   math::XYZTLorentzVector ReclusterJets::getPruned(const fastjet::PseudoJet &jet, double zcut, double rcutFactor) {
0101     // create pruner
0102     fastjet::Pruner pruner(fastjet::cambridge_algorithm, zcut, rcutFactor);
0103     // Prune the jet
0104     fastjet::PseudoJet pruned_jet = pruner(jet);
0105     return LorentzVector(pruned_jet.px(), pruned_jet.py(), pruned_jet.pz(), pruned_jet.e());
0106   }
0107
0108 }  // namespace heppy