FastJet 3.0.2
Pruner.cc
00001 //STARTHEADER
00002 // $Id: Pruner.cc 2689 2011-11-14 14:51:06Z soyez $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 #include "fastjet/tools/Pruner.hh"
00030 #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
00031 #include "fastjet/Selector.hh"
00032 #include <cassert>
00033 #include <algorithm>
00034 #include <sstream>
00035 #include <typeinfo>
00036 
00037 using namespace std;
00038 
00039 
00040 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00041 
00042 
00043 //----------------------------------------------------------------------
00044 // class Pruner
00045 //----------------------------------------------------------------------
00046 
00047 //----------------------------------------------------------------------
00048 // alternative (dynamic) ctor
00049 //  \param jet_def the jet definition for the internal clustering
00050 //  \param zcut_dyn    dynamic pt-fraction cut in the pruning
00051 //  \param Rcut_dyn    dynamic angular distance cut in the pruning
00052 Pruner::Pruner(const JetDefinition &jet_def, 
00053          FunctionOfPseudoJet<double> *zcut_dyn,
00054          FunctionOfPseudoJet<double> *Rcut_dyn)
00055   : _jet_def(jet_def), _zcut(0), _Rcut_factor(0),
00056     _zcut_dyn(zcut_dyn), _Rcut_dyn(Rcut_dyn), _get_recombiner_from_jet(false)  {
00057   assert(_zcut_dyn != 0 && _Rcut_dyn != 0);
00058 }
00059 
00060 //----------------------------------------------------------------------
00061 // action on a single jet
00062 PseudoJet Pruner::result(const PseudoJet &jet) const{
00063   // pruning can only be applied to jets that have constituents
00064   if (!jet.has_constituents()){
00065     throw Error("Pruner: trying to apply the Pruner transformer to a jet that has no constituents");
00066   }
00067 
00068   // if the jet has area support and there are explicit ghosts, we can
00069   // transfer that support to the internal re-clustering
00070   bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
00071 
00072   // build the pruning plugin
00073   double Rcut = (_Rcut_dyn) ? (*_Rcut_dyn)(jet) : _Rcut_factor * 2.0*jet.m()/jet.perp();
00074   double zcut = (_zcut_dyn) ? (*_zcut_dyn)(jet) : _zcut;
00075   PruningPlugin * pruning_plugin;
00076   // for some constructors, we get the recombiner from the 
00077   // input jet -- some acrobatics are needed (see plans for FJ3.1
00078   // for a hopefully better solution).
00079   if (_get_recombiner_from_jet) {
00080     const JetDefinition::Recombiner * common_recombiner = 
00081                                               _get_common_recombiner(jet);
00082     if (common_recombiner) {
00083       JetDefinition jet_def = _jet_def;
00084       if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
00085         RecombinationScheme scheme = 
00086           static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
00087         jet_def.set_recombination_scheme(scheme);
00088       } else {
00089         jet_def.set_recombiner(common_recombiner);
00090       }
00091       pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
00092     } else {
00093       // if there wasn't a common recombiner, we just use the default
00094       // recombiner that was in _jet_def
00095       pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
00096     }
00097   } else {
00098     pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
00099   }
00100 
00101   // now recluster the constituents of the jet with that plugin
00102   JetDefinition internal_jet_def(pruning_plugin);
00103   // flag the plugin for automatic deletion _before_ we make
00104   // copies (so that as long as the copies are also present
00105   // it doesn't get deleted).
00106   internal_jet_def.delete_plugin_when_unused();
00107 
00108   ClusterSequence * cs;
00109   if (do_areas){
00110     vector<PseudoJet> particles, ghosts;
00111     SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
00112     // determine the ghost area from the 1st ghost (if none, any value
00113     // will do, as the area will be 0 and subtraction will have
00114     // no effect!)
00115     double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
00116     cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def, 
00117                                                      ghosts, ghost_area);
00118   } else {
00119     cs = new ClusterSequence(jet.constituents(), internal_jet_def);
00120   }
00121   
00122   PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
00123   PrunerStructure * s = new PrunerStructure(result_local);
00124   result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
00125   
00126   // make sure things remain persistent -- i.e. tell the jet definition
00127   // and the cluster sequence that it is their responsibility to clean 
00128   // up memory once the "result" reaches the end of its life in the user's
00129   // code. (The CS deletes itself when the result goes out of scope and
00130   // that also triggers deletion of the plugin)
00131   cs->delete_self_when_unused();
00132 
00133   return result_local;  
00134 }
00135 
00136 // check if the jet has explicit_ghosts (knowing that there is an
00137 // area support)
00138 bool Pruner::_check_explicit_ghosts(const PseudoJet &jet) const{
00139   // if the jet comes from a Clustering check explicit ghosts in that
00140   // clustering
00141   if (jet.has_associated_cluster_sequence())
00142     return jet.validated_csab()->has_explicit_ghosts();
00143 
00144   // if the jet has pieces, recurse in the pieces
00145   if (jet.has_pieces()){
00146     vector<PseudoJet> pieces = jet.pieces();
00147     for (unsigned int i=0;i<pieces.size(); i++)
00148       if (!_check_explicit_ghosts(pieces[i])) return false;
00149     // never returned false, so we're OK.
00150     return true;
00151   }
00152 
00153   // return false for any other (unknown) structure
00154   return false;
00155 }
00156 
00157 // see if there is a common recombiner among the pieces; if there
00158 // is return a pointer to it; otherwise, return NULL.
00159 // 
00160 // NB: this way of doing things is not ideal, because quite some work
00161 //     is needed to get a correct handling of the final recombiner
00162 //     (e.g. default v. non-default). In future add
00163 //     set_recombiner(jet_def) to JetDefinition, maybe also add
00164 //     an invalid_scheme to the default recombiner and then 
00165 //     do all the work below directly with a JetDefinition directly
00166 //     together with JD::has_same_recombiner(...)
00167 const JetDefinition::Recombiner * Pruner::_get_common_recombiner(const PseudoJet &jet) const{
00168   if (jet.has_associated_cluster_sequence())
00169     return jet.validated_cs()->jet_def().recombiner();
00170 
00171   // if the jet has pieces, recurse in the pieces
00172   if (jet.has_pieces()){
00173     vector<PseudoJet> pieces = jet.pieces();
00174     if (pieces.size() == 0) return 0;
00175     const JetDefinition::Recombiner * reco = _get_common_recombiner(pieces[0]);
00176     for (unsigned int i=1;i<pieces.size(); i++)
00177       if (_get_common_recombiner(pieces[i]) != reco) return 0;
00178     // never returned false, so we're OK.
00179     return reco;
00180   }
00181 
00182   // return false for any other (unknown) structure
00183   return 0;
00184 }
00185 
00186 
00187 // transformer description
00188 std::string Pruner::description() const{
00189   ostringstream oss;
00190   oss << "Pruner with jet_definition = (" << _jet_def.description() << ")";
00191   if (_zcut_dyn) {
00192     oss << ", dynamic zcut (" << _zcut_dyn->description() << ")"
00193         << ", dynamic Rcut (" << _Rcut_dyn->description() << ")";
00194   } else {
00195     oss << ", zcut = " << _zcut
00196         << ", Rcut_factor = " << _Rcut_factor;
00197   }
00198   return oss.str();
00199 }
00200 
00201 
00202 
00203 //----------------------------------------------------------------------
00204 // class PrunerStructure
00205 //----------------------------------------------------------------------
00206 
00207 // return the other jets that may have been found along with the
00208 // result of the pruning
00209 // The resulting vector is sorted in pt
00210 vector<PseudoJet> PrunerStructure::extra_jets() const{ 
00211   return sorted_by_pt((!SelectorNHardest(1))(validated_cs()->inclusive_jets()));;
00212 }
00213 
00214 
00215 //----------------------------------------------------------------------
00216 // class PruningRecombiner
00217 //----------------------------------------------------------------------
00218 
00219 // decide whether to recombine things or not
00220 void PruningRecombiner::recombine(const PseudoJet &pa, 
00221                                   const PseudoJet &pb,
00222                                   PseudoJet &pab) const{
00223   PseudoJet p;
00224   _recombiner->recombine(pa, pb, p);
00225 
00226   // if the 2 particles are close enough, do the recombination
00227   if (pa.squared_distance(pb)<=_Rcut2){
00228     pab=p; return;
00229   }
00230 
00231   double pt2a = pa.perp2();
00232   double pt2b = pb.perp2();
00233 
00234   // check which is the softest
00235   if (pt2a < pt2b){
00236     if (pt2a<_zcut2*p.perp2()){
00237       pab = pb; _rejected.push_back(pa.cluster_hist_index());
00238     } else {
00239       pab = p;
00240     }
00241   } else {
00242     if (pt2b<_zcut2*p.perp2()) {
00243       pab = pa; _rejected.push_back(pb.cluster_hist_index());
00244     } else {
00245       pab = p;
00246     }
00247   }
00248 }
00249 
00250 // description
00251 string PruningRecombiner::description() const{
00252   ostringstream oss;
00253   oss << "Pruning recombiner with zcut = " << sqrt(_zcut2)
00254       << ", Rcut = " << sqrt(_Rcut2)
00255       << ", and underlying recombiner = " << _recombiner->description();
00256   return oss.str();
00257 }
00258 
00259 
00260 
00261 
00262 //----------------------------------------------------------------------
00263 // class PruningPlugin
00264 //----------------------------------------------------------------------
00265 // the actual clustering work for the plugin
00266 void PruningPlugin::run_clustering(ClusterSequence &input_cs) const{
00267   // declare a pruning recombiner
00268   PruningRecombiner pruning_recombiner(_zcut, _Rcut, _jet_def.recombiner());
00269   JetDefinition jet_def = _jet_def;
00270   jet_def.set_recombiner(&pruning_recombiner);
00271 
00272   // cluster the particles using that recombiner
00273   ClusterSequence internal_cs(input_cs.jets(), jet_def);
00274   const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
00275 
00276   // transfer the list of "childless" elements into a bool vector
00277   vector<bool> kept(internal_hist.size(), true);
00278   const vector<unsigned int> &pr_rej = pruning_recombiner.rejected();
00279   for (unsigned int i=0;i<pr_rej.size(); i++) kept[pr_rej[i]]=false;
00280 
00281   // browse the history, keeping only the elements that have not been
00282   // vetoed.
00283   //
00284   // In the process we build a map for the history indices
00285   vector<unsigned int> internal2input(internal_hist.size());
00286   for (unsigned int i=0; i<input_cs.jets().size(); i++)
00287     internal2input[i] = i;
00288 
00289   for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
00290     const ClusterSequence::history_element &he = internal_hist[i];
00291 
00292     // deal with recombinations with the beam
00293     if (he.parent2 == ClusterSequence::BeamJet){
00294       int internal_jetp_index = internal_hist[he.parent1].jetp_index;
00295       int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
00296 
00297       int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
00298 
00299       // cout << "Beam recomb for internal " << internal_hist_index
00300       //            << " (input jet index=" << input_jetp_index << endl;
00301 
00302       input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
00303       continue;
00304     }
00305 
00306     // now, deal with two-body recombinations
00307     if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
00308       internal2input[i]=internal2input[he.parent2];
00309       // cout << "rejecting internal " << he.parent1
00310       //            << ", mapping internal " << i 
00311       //            << " to internal " << he.parent2
00312       //            << " i.e. " << internal2input[i] << endl;
00313     } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
00314       internal2input[i]=internal2input[he.parent1];
00315       // cout << "rejecting internal " << he.parent2 
00316       //            << ", mapping internal " << i 
00317       //            << " to internal " << he.parent1
00318       //            << " i.e. " << internal2input[i] << endl;
00319     } else { // do the recombination
00320       int new_index;
00321       input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
00322                                               input_cs.history()[internal2input[he.parent2]].jetp_index,
00323                                               he.dij, internal_cs.jets()[he.jetp_index], new_index);
00324       internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
00325       // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
00326       //      << " and "    << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
00327       //      << " into "   << internal2input[i] << " (int: " << i << ")" << endl;
00328     }
00329   }
00330 }
00331 
00332 // returns the plugin description
00333 string PruningPlugin::description() const{
00334   ostringstream oss;
00335   oss << "Pruning plugin with jet_definition = (" << _jet_def.description()
00336       <<"), zcut = " << _zcut
00337       << ", Rcut = " << _Rcut;
00338   return oss.str();
00339 }
00340 
00341 
00342 FASTJET_END_NAMESPACE
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends