FastJet 3.0.2
|
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