File indexing completed on 2023-03-17 11:18:22
0001 #include <memory>
0002
0003 #include "RecoJets/JetAlgorithms/interface/SubJetAlgorithm.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "fastjet/ClusterSequenceArea.hh"
0007
0008 using namespace std;
0009 using namespace edm;
0010
0011 void SubJetAlgorithm::set_zcut(double z) { zcut_ = z; }
0012
0013 void SubJetAlgorithm::set_rcut_factor(double r) { rcut_factor_ = r; }
0014
0015
0016
0017 void SubJetAlgorithm::run(const vector<fastjet::PseudoJet>& cell_particles, vector<CompoundPseudoJet>& hardjetsOutput) {
0018
0019
0020 fastjet::FastPrunePlugin PRplugin(*fjJetDefinition_, *fjJetDefinition_, zcut_, rcut_factor_);
0021 fastjet::JetDefinition pjetdef(&PRplugin);
0022
0023
0024
0025 std::shared_ptr<fastjet::ClusterSequence> fjClusterSeq;
0026 if (!doAreaFastjet_) {
0027 fjClusterSeq = std::make_shared<fastjet::ClusterSequence>(cell_particles, pjetdef);
0028 } else if (voronoiRfact_ <= 0) {
0029 fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>(
0030 new fastjet::ClusterSequenceActiveArea(cell_particles, pjetdef, *fjActiveArea_));
0031 } else {
0032 fjClusterSeq = std::shared_ptr<fastjet::ClusterSequence>(
0033 new fastjet::ClusterSequenceVoronoiArea(cell_particles, pjetdef, fastjet::VoronoiAreaSpec(voronoiRfact_)));
0034 }
0035
0036 vector<fastjet::PseudoJet> inclusiveJets = fjClusterSeq->inclusive_jets(ptMin_);
0037
0038
0039
0040 vector<vector<int> > indices(inclusiveJets.size());
0041
0042 vector<fastjet::PseudoJet>::iterator jetIt = inclusiveJets.begin();
0043 for (; jetIt != inclusiveJets.end(); ++jetIt) {
0044
0045 vector<fastjet::PseudoJet> subjets = fjClusterSeq->exclusive_subjets(*jetIt, nSubjets_);
0046
0047 vector<CompoundPseudoSubJet> subjetsOutput;
0048 std::vector<fastjet::PseudoJet>::const_iterator itSubJetBegin = subjets.begin(), itSubJet = itSubJetBegin,
0049 itSubJetEnd = subjets.end();
0050 for (; itSubJet != itSubJetEnd; ++itSubJet) {
0051
0052 vector<fastjet::PseudoJet> subjetFastjetConstituents = fjClusterSeq->constituents(*itSubJet);
0053
0054 vector<int> constituents;
0055 vector<fastjet::PseudoJet>::const_iterator fastSubIt = subjetFastjetConstituents.begin(),
0056 transConstEnd = subjetFastjetConstituents.end();
0057 for (; fastSubIt != transConstEnd; ++fastSubIt) {
0058 if (fastSubIt->user_index() >= 0) {
0059 constituents.push_back(fastSubIt->user_index());
0060 }
0061 }
0062
0063 double subJetArea =
0064 (doAreaFastjet_) ? dynamic_cast<fastjet::ClusterSequenceActiveArea&>(*fjClusterSeq).area(*itSubJet) : 0.0;
0065
0066
0067 subjetsOutput.push_back(CompoundPseudoSubJet(*itSubJet, subJetArea, constituents));
0068 }
0069
0070 double fatJetArea =
0071 (doAreaFastjet_) ? dynamic_cast<fastjet::ClusterSequenceActiveArea&>(*fjClusterSeq).area(*jetIt) : 0.0;
0072
0073
0074 hardjetsOutput.push_back(CompoundPseudoJet(*jetIt, fatJetArea, subjetsOutput));
0075 }
0076 }