JetAlgorithms/src/SubJetAlgorithm.cc

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 //  Run the algorithm
0016 //  ------------------
0017 void SubJetAlgorithm::run(const vector<fastjet::PseudoJet>& cell_particles, vector<CompoundPseudoJet>& hardjetsOutput) {
0018   //for actual jet clustering, either the pruned or the original version is used.
0019   //For the pruned version, a new jet definition using the PrunedRecombPlugin is required:
0020   fastjet::FastPrunePlugin PRplugin(*fjJetDefinition_, *fjJetDefinition_, zcut_, rcut_factor_);
0021   fastjet::JetDefinition pjetdef(&PRplugin);
0022
0023   // cluster the jets with the jet definition jetDef:
0024   // run algorithm
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   // These will store the indices of each subjet that
0039   // are present in each jet
0040   vector<vector<int> > indices(inclusiveJets.size());
0041   // Loop over inclusive jets, attempt to find substructure
0042   vector<fastjet::PseudoJet>::iterator jetIt = inclusiveJets.begin();
0043   for (; jetIt != inclusiveJets.end(); ++jetIt) {
0044     //decompose into requested number of subjets:
0045     vector<fastjet::PseudoJet> subjets = fjClusterSeq->exclusive_subjets(*jetIt, nSubjets_);
0046     //create the subjets objects to put into the "output" objects
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       // Get the transient subjet constituents from fastjet
0052       vector<fastjet::PseudoJet> subjetFastjetConstituents = fjClusterSeq->constituents(*itSubJet);
0053       // Get the indices of the constituents:
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       // Make a CompoundPseudoSubJet object to hold this subjet and the indices of its constituents
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     // Make a CompoundPseudoJet object to hold this hard jet, and the subjets that make it up
0074     hardjetsOutput.push_back(CompoundPseudoJet(*jetIt, fatJetArea, subjetsOutput));
0075   }
0076 }