File indexing completed on 2023-03-17 11:15:49
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
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
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++;
0023 fjInputs_.push_back(j);
0024 }
0025
0026
0027 fastjet::JetDefinition jet_def;
0028
0029
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
0041
0042
0043
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
0057 inclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(ptMin));
0058
0059 return makeP4s(inclusiveJets_);
0060 }
0061
0062 std::vector<math::XYZTLorentzVector> ReclusterJets::getGroupingExclusive(double dcut) {
0063
0064 exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(dcut));
0065
0066 return makeP4s(exclusiveJets_);
0067 }
0068
0069 std::vector<math::XYZTLorentzVector> ReclusterJets::getGroupingExclusive(int njets) {
0070
0071 exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(njets));
0072
0073 return makeP4s(exclusiveJets_);
0074 }
0075
0076 math::XYZTLorentzVector ReclusterJets::getPruned(double zcut, double rcutFactor) {
0077
0078 exclusiveJets_ = fastjet::sorted_by_pt(fjClusterSeq_->exclusive_jets(1));
0079
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
0102 fastjet::Pruner pruner(fastjet::cambridge_algorithm, zcut, rcutFactor);
0103
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 }