JetAlgorithms/src/FastPrunePlugin.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Implements the FastPrunePlugin class.  See FastPrunePlugin.hh for a
0004 //   description.
0005 //
0006 ///////////////////////////////////////////////////////////////////////////////
0007
0008 #include "RecoJets/JetAlgorithms/interface/FastPrunePlugin.hh"
0009
0010 #include <sstream>
0011 #include <cmath>
0012 #include <vector>
0013 #include <algorithm>
0014 using namespace std;
0015
0016 using namespace fastjet;
0017
0018 FastPrunePlugin::FastPrunePlugin(const JetDefinition& find_definition,
0019                                  const JetDefinition& prune_definition,
0020                                  const double& zcut,
0021                                  const double& Rcut_factor)
0022     : _find_definition(find_definition),
0023       _prune_definition(prune_definition),
0024       _minpT(20.),
0025       _unpruned_seq(nullptr),
0026       _pruned_recombiner(nullptr),
0027       _cut_setter(new DefaultCutSetter(zcut, Rcut_factor)) {
0028   // If the passed prune_definition (copied into _prune_def) has an external
0029   // recombiner, use it.  Otherwise, create our own DefaultRecombiner.  We
0030   // could just point to _prune_definition's DefaultRecombiner, but FastJet
0031   // actually sets this to use external_scheme later when we set
0032   // _prune_definition's Recombiner to be a PrunedRecombiner.  (Calling
0033   // JetDefinition::set_recombiner() also calls jetdef._default_recombiner =
0034   // DefaultRecombiner(external_scheme) as a sanity check to keep you from
0035   // using it.)
0036   RecombinationScheme scheme = _prune_definition.recombination_scheme();
0037   if (scheme == external_scheme)
0038     _pruned_recombiner = new PrunedRecombiner(_prune_definition.recombiner(), zcut, 0.0);
0039   else
0040     _pruned_recombiner = new PrunedRecombiner(scheme, zcut, 0.0);
0041 }
0042
0043 FastPrunePlugin::FastPrunePlugin(const JetDefinition& find_definition,
0044                                  const JetDefinition& prune_definition,
0045                                  const JetDefinition::Recombiner* recomb,
0046                                  const double& zcut,
0047                                  const double& Rcut_factor)
0048     : _find_definition(find_definition),
0049       _prune_definition(prune_definition),
0050       _minpT(20.),
0051       _unpruned_seq(nullptr),
0052       _pruned_recombiner(new PrunedRecombiner(recomb, zcut, 0.0)),
0053       _cut_setter(new DefaultCutSetter(zcut, Rcut_factor)) {}
0054
0055 FastPrunePlugin::FastPrunePlugin(const JetDefinition& find_definition,
0056                                  const JetDefinition& prune_definition,
0057                                  CutSetter* const cut_setter)
0058     : _find_definition(find_definition),
0059       _prune_definition(prune_definition),
0060       _minpT(20.),
0061       _unpruned_seq(nullptr),
0062       _pruned_recombiner(nullptr),
0063       _cut_setter(cut_setter) {
0064   // See comments in first constructor
0065   RecombinationScheme scheme = _prune_definition.recombination_scheme();
0066   if (scheme == external_scheme)
0067     _pruned_recombiner = new PrunedRecombiner(_prune_definition.recombiner());
0068   else
0069     _pruned_recombiner = new PrunedRecombiner(scheme);
0070 }
0071
0072 FastPrunePlugin::FastPrunePlugin(const JetDefinition& find_definition,
0073                                  const JetDefinition& prune_definition,
0074                                  CutSetter* const cut_setter,
0075                                  const JetDefinition::Recombiner* recomb)
0076     : _find_definition(find_definition),
0077       _prune_definition(prune_definition),
0078       _minpT(20.),
0079       _unpruned_seq(nullptr),
0080       _pruned_recombiner(new PrunedRecombiner(recomb)),
0081       _cut_setter(cut_setter) {}
0082
0083 string FastPrunePlugin::description() const {
0084   ostringstream desc;
0085
0086   desc << "Pruned jet algorithm \n"
0087        << "----------------------- \n"
0088        << "Jet finder: " << _find_definition.description() << "\n----------------------- \n"
0089        << "Prune jets with: " << _prune_definition.description() << "\n----------------------- \n"
0090        << "Pruning parameters: "
0091        << "zcut = " << _cut_setter->zcut << ", "
0092        << "Rcut_factor = ";
0093   if (DefaultCutSetter* cs = dynamic_cast<DefaultCutSetter*>(_cut_setter))
0094     desc << cs->Rcut_factor;
0095   else
0096     desc << "[dynamic]";
0097   desc << "\n"
0098        << "----------------------- \n";
0099
0100   return desc.str();
0101 }
0102
0103 // Meat of the plugin.  This function takes the input particles from input_seq
0104 //   and calls plugin_record_{ij, iB}_recombination repeatedly such that
0105 //   input_seq holds a cluster sequence corresponding to pruned jets
0106 // This function first creates unpruned jets using _find_definition, then uses
0107 //   the prune_definition, along with the PrunedRecombiner, on each jet to
0108 //   produce a pruned jet.  For each of these jets, output_mergings() is called,
0109 //   which reads the history of the pruned sequence and calls the
0110 //   plugin_record functions of input_seq to match this history.
0111 // Branches that are pruned appear in the history as PseudoJets with no
0112 //   children.  For this reason, only inclusive_jets() is sensible to use with
0113 //   the final ClusterSequence.  The substructure, e.g., constituents() of a
0114 //   pruned jet will not include the pruned away branchings.
0115 void FastPrunePlugin::run_clustering(ClusterSequence& input_seq) const {
0116   // this will be filled in the output_mergings() step
0117   _pruned_subjets.clear();
0118
0119   vector<PseudoJet> inputs = input_seq.jets();
0120   // Record user_index's so we can match PJ's in pruned jets to PJ's in
0121   //   input_seq.
0122   // Use i+1 to distinguish from the default, which in some places appears to
0123   //   be 0.
0124   // Note that we're working with a local copy of the input jets, so we don't
0125   //   change the indices in input_seq.
0126   for (unsigned int i = 0; i < inputs.size(); i++)
0127     inputs[i].set_user_index(i + 1);
0128
0129   // ClusterSequence for initial (unpruned) jet finding
0130   if (_unpruned_seq)
0131     delete _unpruned_seq;
0132   _unpruned_seq = new ClusterSequence(inputs, _find_definition);
0133
0134   // now, for each jet, construct a pruned version:
0135   vector<PseudoJet> unpruned_jets = sorted_by_pt(_unpruned_seq->inclusive_jets(_minpT));
0136
0137   for (unsigned int i = 0; i < unpruned_jets.size(); i++) {
0138     _cut_setter->SetCuts(unpruned_jets[i], *_unpruned_seq);
0139     _pruned_recombiner->reset(_cut_setter->zcut, _cut_setter->Rcut);
0140     _prune_definition.set_recombiner(_pruned_recombiner);
0141
0142     // temp way to get constituents, to compare to new version
0143     vector<PseudoJet> constituents;
0144     for (size_t j = 0; j < inputs.size(); ++j)
0145       if (_unpruned_seq->object_in_jet(inputs[j], unpruned_jets[i])) {
0146         constituents.push_back(inputs[j]);
0147       }
0148     ClusterSequence pruned_seq(constituents, _prune_definition);
0149
0150     vector<int> pruned_PJs = _pruned_recombiner->pruned_pseudojets();
0151     _output_mergings(pruned_seq, pruned_PJs, input_seq);
0152   }
0153 }
0154
0155 // Takes the merging structure of "in_seq" and feeds this into "out_seq" using
0156 //   _plugin_record_{ij,iB}_recombination().
0157 // PJ's in the input CS should have user_index's that are their (index+1) in
0158 //    _jets() in the output CS (the output CS should be the input CS to
0159 //    run_clustering()).
0160 // This allows us to build up the same jet in out_seq as already exists in
0161 //   in_seq.
0162 void FastPrunePlugin::_output_mergings(ClusterSequence& in_seq,
0163                                        vector<int>& pruned_pseudojets,
0164                                        ClusterSequence& out_seq) const {
0165   // vector to store the pruned subjets for this jet
0166   vector<PseudoJet> temp_pruned_subjets;
0167
0168   vector<PseudoJet> p = in_seq.jets();
0169
0170   // sort this vector so we can binary search it
0171   sort(pruned_pseudojets.begin(), pruned_pseudojets.end());
0172
0173   // get the history from in_seq
0174   const vector<ClusterSequence::history_element>& hist = in_seq.history();
0175   vector<ClusterSequence::history_element>::const_iterator iter = hist.begin();
0176
0177   // skip particle input elements
0178   while (iter->parent1 == ClusterSequence::InexistentParent)
0179     iter++;
0180
0181   // Walk through history.  PseudoJets in in_seq should have a user_index
0182   //   corresponding to their index in out_seq.  Note that when we create a
0183   //   new PJ via record_ij we need to set the user_index of our local copy.
0184   for (; iter != hist.end(); iter++) {
0185     int new_jet_index = -1;
0186     int jet_index = iter->jetp_index;
0187     int parent1 = iter->parent1;
0188     int parent2 = iter->parent2;
0189     int parent1_index = p[hist[iter->parent1].jetp_index].user_index() - 1;
0190
0191     if (parent2 == ClusterSequence::BeamJet) {
0192       out_seq.plugin_record_iB_recombination(parent1_index, iter->dij);
0193     } else {
0194       int parent2_index = p[hist[parent2].jetp_index].user_index() - 1;
0195
0196       // Check if either parent is stored in pruned_pseudojets
0197       //   Note that it is the history index that is stored
0198       if (binary_search(pruned_pseudojets.begin(), pruned_pseudojets.end(), parent2)) {
0199         // pruned away parent2 -- just give child parent1's index
0200         p[jet_index].set_user_index(p[hist[parent1].jetp_index].user_index());
0201         temp_pruned_subjets.push_back(out_seq.jets()[parent2_index]);
0202       } else if (binary_search(pruned_pseudojets.begin(), pruned_pseudojets.end(), parent1)) {
0203         // pruned away parent1 -- just give child parent2's index
0204         p[jet_index].set_user_index(p[hist[parent2].jetp_index].user_index());
0205         temp_pruned_subjets.push_back(out_seq.jets()[parent1_index]);
0206       } else {
0207         // no pruning -- record combination and index for child
0208         out_seq.plugin_record_ij_recombination(parent1_index, parent2_index, iter->dij, new_jet_index);
0209         p[jet_index].set_user_index(new_jet_index + 1);
0210       }
0211     }
0212   }
0213
0214   _pruned_subjets.push_back(temp_pruned_subjets);
0215 }