| 1 | //FJSTARTHEADER
|
|---|
| 2 | // $Id: Pruner.cc 4354 2018-04-22 07:12:37Z salam $
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (c) 2005-2018, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
|
|---|
| 5 | //
|
|---|
| 6 | //----------------------------------------------------------------------
|
|---|
| 7 | // This file is part of FastJet.
|
|---|
| 8 | //
|
|---|
| 9 | // FastJet is free software; you can redistribute it and/or modify
|
|---|
| 10 | // it under the terms of the GNU General Public License as published by
|
|---|
| 11 | // the Free Software Foundation; either version 2 of the License, or
|
|---|
| 12 | // (at your option) any later version.
|
|---|
| 13 | //
|
|---|
| 14 | // The algorithms that underlie FastJet have required considerable
|
|---|
| 15 | // development. They are described in the original FastJet paper,
|
|---|
| 16 | // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
|
|---|
| 17 | // FastJet as part of work towards a scientific publication, please
|
|---|
| 18 | // quote the version you use and include a citation to the manual and
|
|---|
| 19 | // optionally also to hep-ph/0512210.
|
|---|
| 20 | //
|
|---|
| 21 | // FastJet is distributed in the hope that it will be useful,
|
|---|
| 22 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 23 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 24 | // GNU General Public License for more details.
|
|---|
| 25 | //
|
|---|
| 26 | // You should have received a copy of the GNU General Public License
|
|---|
| 27 | // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 28 | //----------------------------------------------------------------------
|
|---|
| 29 | //FJENDHEADER
|
|---|
| 30 |
|
|---|
| 31 | #include "fastjet/tools/Pruner.hh"
|
|---|
| 32 | #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
|
|---|
| 33 | #include "fastjet/Selector.hh"
|
|---|
| 34 | #include <cassert>
|
|---|
| 35 | #include <algorithm>
|
|---|
| 36 | #include <sstream>
|
|---|
| 37 | #include <typeinfo>
|
|---|
| 38 |
|
|---|
| 39 | using namespace std;
|
|---|
| 40 |
|
|---|
| 41 |
|
|---|
| 42 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 43 |
|
|---|
| 44 |
|
|---|
| 45 | //----------------------------------------------------------------------
|
|---|
| 46 | // class Pruner
|
|---|
| 47 | //----------------------------------------------------------------------
|
|---|
| 48 |
|
|---|
| 49 | //----------------------------------------------------------------------
|
|---|
| 50 | // alternative (dynamic) ctor
|
|---|
| 51 | // \param jet_def the jet definition for the internal clustering
|
|---|
| 52 | // \param zcut_dyn dynamic pt-fraction cut in the pruning
|
|---|
| 53 | // \param Rcut_dyn dynamic angular distance cut in the pruning
|
|---|
| 54 | Pruner::Pruner(const JetDefinition &jet_def,
|
|---|
| 55 | const FunctionOfPseudoJet<double> *zcut_dyn,
|
|---|
| 56 | const FunctionOfPseudoJet<double> *Rcut_dyn)
|
|---|
| 57 | : _jet_def(jet_def), _zcut(0), _Rcut_factor(0),
|
|---|
| 58 | _zcut_dyn(zcut_dyn), _Rcut_dyn(Rcut_dyn), _get_recombiner_from_jet(false) {
|
|---|
| 59 | assert(_zcut_dyn != 0 && _Rcut_dyn != 0);
|
|---|
| 60 | }
|
|---|
| 61 |
|
|---|
| 62 | //----------------------------------------------------------------------
|
|---|
| 63 | // action on a single jet
|
|---|
| 64 | PseudoJet Pruner::result(const PseudoJet &jet) const{
|
|---|
| 65 | // pruning can only be applied to jets that have constituents
|
|---|
| 66 | if (!jet.has_constituents()){
|
|---|
| 67 | throw Error("Pruner: trying to apply the Pruner transformer to a jet that has no constituents");
|
|---|
| 68 | }
|
|---|
| 69 |
|
|---|
| 70 | // if the jet has area support and there are explicit ghosts, we can
|
|---|
| 71 | // transfer that support to the internal re-clustering
|
|---|
| 72 | bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
|
|---|
| 73 |
|
|---|
| 74 | // build the pruning plugin
|
|---|
| 75 | double Rcut = (_Rcut_dyn) ? (*_Rcut_dyn)(jet) : _Rcut_factor * 2.0*jet.m()/jet.perp();
|
|---|
| 76 | double zcut = (_zcut_dyn) ? (*_zcut_dyn)(jet) : _zcut;
|
|---|
| 77 | PruningPlugin * pruning_plugin;
|
|---|
| 78 |
|
|---|
| 79 | // for some constructors, we get the recombiner from the
|
|---|
| 80 | // input jet -- some acrobatics are needed
|
|---|
| 81 | if (_get_recombiner_from_jet) {
|
|---|
| 82 | JetDefinition jet_def = _jet_def;
|
|---|
| 83 |
|
|---|
| 84 | // if all the pieces have a shared recombiner, we'll use that
|
|---|
| 85 | // one. Otherwise, use the one from _jet_def as a fallback.
|
|---|
| 86 | JetDefinition jet_def_for_recombiner;
|
|---|
| 87 | if (_check_common_recombiner(jet, jet_def_for_recombiner)){
|
|---|
| 88 | jet_def.set_recombiner(jet_def_for_recombiner);
|
|---|
| 89 | }
|
|---|
| 90 | pruning_plugin = new PruningPlugin(jet_def, zcut, Rcut);
|
|---|
| 91 | } else {
|
|---|
| 92 | pruning_plugin = new PruningPlugin(_jet_def, zcut, Rcut);
|
|---|
| 93 | }
|
|---|
| 94 |
|
|---|
| 95 | // now recluster the constituents of the jet with that plugin
|
|---|
| 96 | JetDefinition internal_jet_def(pruning_plugin);
|
|---|
| 97 | // flag the plugin for automatic deletion _before_ we make
|
|---|
| 98 | // copies (so that as long as the copies are also present
|
|---|
| 99 | // it doesn't get deleted).
|
|---|
| 100 | internal_jet_def.delete_plugin_when_unused();
|
|---|
| 101 |
|
|---|
| 102 | ClusterSequence * cs;
|
|---|
| 103 | if (do_areas){
|
|---|
| 104 | vector<PseudoJet> particles, ghosts;
|
|---|
| 105 | SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
|
|---|
| 106 | // determine the ghost area from the 1st ghost (if none, any value
|
|---|
| 107 | // will do, as the area will be 0 and subtraction will have
|
|---|
| 108 | // no effect!)
|
|---|
| 109 | double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
|
|---|
| 110 | cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def,
|
|---|
| 111 | ghosts, ghost_area);
|
|---|
| 112 | } else {
|
|---|
| 113 | cs = new ClusterSequence(jet.constituents(), internal_jet_def);
|
|---|
| 114 | }
|
|---|
| 115 |
|
|---|
| 116 | PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
|
|---|
| 117 | PrunerStructure * s = new PrunerStructure(result_local);
|
|---|
| 118 | s->_Rcut = Rcut;
|
|---|
| 119 | s->_zcut = zcut;
|
|---|
| 120 | result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
|
|---|
| 121 |
|
|---|
| 122 | // make sure things remain persistent -- i.e. tell the jet definition
|
|---|
| 123 | // and the cluster sequence that it is their responsibility to clean
|
|---|
| 124 | // up memory once the "result" reaches the end of its life in the user's
|
|---|
| 125 | // code. (The CS deletes itself when the result goes out of scope and
|
|---|
| 126 | // that also triggers deletion of the plugin)
|
|---|
| 127 | cs->delete_self_when_unused();
|
|---|
| 128 |
|
|---|
| 129 | return result_local;
|
|---|
| 130 | }
|
|---|
| 131 |
|
|---|
| 132 | // check if the jet has explicit_ghosts (knowing that there is an
|
|---|
| 133 | // area support)
|
|---|
| 134 | bool Pruner::_check_explicit_ghosts(const PseudoJet &jet) const{
|
|---|
| 135 | // if the jet comes from a Clustering check explicit ghosts in that
|
|---|
| 136 | // clustering
|
|---|
| 137 | if (jet.has_associated_cluster_sequence())
|
|---|
| 138 | return jet.validated_csab()->has_explicit_ghosts();
|
|---|
| 139 |
|
|---|
| 140 | // if the jet has pieces, recurse in the pieces
|
|---|
| 141 | if (jet.has_pieces()){
|
|---|
| 142 | vector<PseudoJet> pieces = jet.pieces();
|
|---|
| 143 | for (unsigned int i=0;i<pieces.size(); i++)
|
|---|
| 144 | if (!_check_explicit_ghosts(pieces[i])) return false;
|
|---|
| 145 | // never returned false, so we're OK.
|
|---|
| 146 | return true;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | // return false for any other (unknown) structure
|
|---|
| 150 | return false;
|
|---|
| 151 | }
|
|---|
| 152 |
|
|---|
| 153 | // see if there is a common recombiner among the pieces; if there is
|
|---|
| 154 | // return true and set jet_def_for_recombiner so that the recombiner
|
|---|
| 155 | // can be taken from that JetDefinition. Otherwise, return
|
|---|
| 156 | // false. 'assigned' is initially false; when true, each time we meet
|
|---|
| 157 | // a new jet definition, we'll check it shares the same recombiner as
|
|---|
| 158 | // jet_def_for_recombiner.
|
|---|
| 159 | bool Pruner::_check_common_recombiner(const PseudoJet &jet,
|
|---|
| 160 | JetDefinition &jet_def_for_recombiner,
|
|---|
| 161 | bool assigned) const{
|
|---|
| 162 | if (jet.has_associated_cluster_sequence()){
|
|---|
| 163 | // if the jet def for recombination has already been assigned, check if we have the same
|
|---|
| 164 | if (assigned)
|
|---|
| 165 | return jet.validated_cs()->jet_def().has_same_recombiner(jet_def_for_recombiner);
|
|---|
| 166 |
|
|---|
| 167 | // otherwise, assign it.
|
|---|
| 168 | jet_def_for_recombiner = jet.validated_cs()->jet_def();
|
|---|
| 169 | assigned = true;
|
|---|
| 170 | return true;
|
|---|
| 171 | }
|
|---|
| 172 |
|
|---|
| 173 | // if the jet has pieces, recurse in the pieces
|
|---|
| 174 | if (jet.has_pieces()){
|
|---|
| 175 | vector<PseudoJet> pieces = jet.pieces();
|
|---|
| 176 | if (pieces.size() == 0) return false;
|
|---|
| 177 | for (unsigned int i=0;i<pieces.size(); i++)
|
|---|
| 178 | if (!_check_common_recombiner(pieces[i], jet_def_for_recombiner, assigned)) return false;
|
|---|
| 179 | // never returned false, so we're OK.
|
|---|
| 180 | return true;
|
|---|
| 181 | }
|
|---|
| 182 |
|
|---|
| 183 | // return false for any other (unknown) structure
|
|---|
| 184 | return false;
|
|---|
| 185 | }
|
|---|
| 186 |
|
|---|
| 187 |
|
|---|
| 188 | // transformer description
|
|---|
| 189 | std::string Pruner::description() const{
|
|---|
| 190 | ostringstream oss;
|
|---|
| 191 | oss << "Pruner with jet_definition = (" << _jet_def.description() << ")";
|
|---|
| 192 | if (_zcut_dyn) {
|
|---|
| 193 | oss << ", dynamic zcut (" << _zcut_dyn->description() << ")"
|
|---|
| 194 | << ", dynamic Rcut (" << _Rcut_dyn->description() << ")";
|
|---|
| 195 | } else {
|
|---|
| 196 | oss << ", zcut = " << _zcut
|
|---|
| 197 | << ", Rcut_factor = " << _Rcut_factor;
|
|---|
| 198 | }
|
|---|
| 199 | return oss.str();
|
|---|
| 200 | }
|
|---|
| 201 |
|
|---|
| 202 |
|
|---|
| 203 |
|
|---|
| 204 | //----------------------------------------------------------------------
|
|---|
| 205 | // class PrunerStructure
|
|---|
| 206 | //----------------------------------------------------------------------
|
|---|
| 207 |
|
|---|
| 208 | // return the other jets that may have been found along with the
|
|---|
| 209 | // result of the pruning
|
|---|
| 210 | // The resulting vector is sorted in pt
|
|---|
| 211 | vector<PseudoJet> PrunerStructure::extra_jets() const{
|
|---|
| 212 | return sorted_by_pt((!SelectorNHardest(1))(validated_cs()->inclusive_jets()));;
|
|---|
| 213 | }
|
|---|
| 214 |
|
|---|
| 215 |
|
|---|
| 216 | //----------------------------------------------------------------------
|
|---|
| 217 | // class PruningRecombiner
|
|---|
| 218 | //----------------------------------------------------------------------
|
|---|
| 219 |
|
|---|
| 220 | // decide whether to recombine things or not
|
|---|
| 221 | void PruningRecombiner::recombine(const PseudoJet &pa,
|
|---|
| 222 | const PseudoJet &pb,
|
|---|
| 223 | PseudoJet &pab) const{
|
|---|
| 224 | PseudoJet p;
|
|---|
| 225 | _recombiner->recombine(pa, pb, p);
|
|---|
| 226 |
|
|---|
| 227 | // if the 2 particles are close enough, do the recombination
|
|---|
| 228 | if (pa.squared_distance(pb)<=_Rcut2){
|
|---|
| 229 | pab=p; return;
|
|---|
| 230 | }
|
|---|
| 231 |
|
|---|
| 232 | double pt2a = pa.perp2();
|
|---|
| 233 | double pt2b = pb.perp2();
|
|---|
| 234 |
|
|---|
| 235 | // check which is the softest
|
|---|
| 236 | if (pt2a < pt2b){
|
|---|
| 237 | if (pt2a<_zcut2*p.perp2()){
|
|---|
| 238 | pab = pb; _rejected.push_back(pa.cluster_hist_index());
|
|---|
| 239 | } else {
|
|---|
| 240 | pab = p;
|
|---|
| 241 | }
|
|---|
| 242 | } else {
|
|---|
| 243 | if (pt2b<_zcut2*p.perp2()) {
|
|---|
| 244 | pab = pa; _rejected.push_back(pb.cluster_hist_index());
|
|---|
| 245 | } else {
|
|---|
| 246 | pab = p;
|
|---|
| 247 | }
|
|---|
| 248 | }
|
|---|
| 249 | }
|
|---|
| 250 |
|
|---|
| 251 | // description
|
|---|
| 252 | string PruningRecombiner::description() const{
|
|---|
| 253 | ostringstream oss;
|
|---|
| 254 | oss << "Pruning recombiner with zcut = " << sqrt(_zcut2)
|
|---|
| 255 | << ", Rcut = " << sqrt(_Rcut2)
|
|---|
| 256 | << ", and underlying recombiner = " << _recombiner->description();
|
|---|
| 257 | return oss.str();
|
|---|
| 258 | }
|
|---|
| 259 |
|
|---|
| 260 |
|
|---|
| 261 |
|
|---|
| 262 |
|
|---|
| 263 | //----------------------------------------------------------------------
|
|---|
| 264 | // class PruningPlugin
|
|---|
| 265 | //----------------------------------------------------------------------
|
|---|
| 266 | // the actual clustering work for the plugin
|
|---|
| 267 | void PruningPlugin::run_clustering(ClusterSequence &input_cs) const{
|
|---|
| 268 | // declare a pruning recombiner
|
|---|
| 269 | PruningRecombiner pruning_recombiner(_zcut, _Rcut, _jet_def.recombiner());
|
|---|
| 270 | JetDefinition jet_def = _jet_def;
|
|---|
| 271 | jet_def.set_recombiner(&pruning_recombiner);
|
|---|
| 272 |
|
|---|
| 273 | // cluster the particles using that recombiner
|
|---|
| 274 | ClusterSequence internal_cs(input_cs.jets(), jet_def);
|
|---|
| 275 | const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
|
|---|
| 276 |
|
|---|
| 277 | // transfer the list of "childless" elements into a bool vector
|
|---|
| 278 | vector<bool> kept(internal_hist.size(), true);
|
|---|
| 279 | const vector<unsigned int> &pr_rej = pruning_recombiner.rejected();
|
|---|
| 280 | for (unsigned int i=0;i<pr_rej.size(); i++) kept[pr_rej[i]]=false;
|
|---|
| 281 |
|
|---|
| 282 | // browse the history, keeping only the elements that have not been
|
|---|
| 283 | // vetoed.
|
|---|
| 284 | //
|
|---|
| 285 | // In the process we build a map for the history indices
|
|---|
| 286 | vector<unsigned int> internal2input(internal_hist.size());
|
|---|
| 287 | for (unsigned int i=0; i<input_cs.jets().size(); i++)
|
|---|
| 288 | internal2input[i] = i;
|
|---|
| 289 |
|
|---|
| 290 | for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
|
|---|
| 291 | const ClusterSequence::history_element &he = internal_hist[i];
|
|---|
| 292 |
|
|---|
| 293 | // deal with recombinations with the beam
|
|---|
| 294 | if (he.parent2 == ClusterSequence::BeamJet){
|
|---|
| 295 | int internal_jetp_index = internal_hist[he.parent1].jetp_index;
|
|---|
| 296 | int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
|
|---|
| 297 |
|
|---|
| 298 | int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
|
|---|
| 299 |
|
|---|
| 300 | // cout << "Beam recomb for internal " << internal_hist_index
|
|---|
| 301 | // << " (input jet index=" << input_jetp_index << endl;
|
|---|
| 302 |
|
|---|
| 303 | input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
|
|---|
| 304 | continue;
|
|---|
| 305 | }
|
|---|
| 306 |
|
|---|
| 307 | // now, deal with two-body recombinations
|
|---|
| 308 | if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
|
|---|
| 309 | internal2input[i]=internal2input[he.parent2];
|
|---|
| 310 | // cout << "rejecting internal " << he.parent1
|
|---|
| 311 | // << ", mapping internal " << i
|
|---|
| 312 | // << " to internal " << he.parent2
|
|---|
| 313 | // << " i.e. " << internal2input[i] << endl;
|
|---|
| 314 | } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
|
|---|
| 315 | internal2input[i]=internal2input[he.parent1];
|
|---|
| 316 | // cout << "rejecting internal " << he.parent2
|
|---|
| 317 | // << ", mapping internal " << i
|
|---|
| 318 | // << " to internal " << he.parent1
|
|---|
| 319 | // << " i.e. " << internal2input[i] << endl;
|
|---|
| 320 | } else { // do the recombination
|
|---|
| 321 | int new_index;
|
|---|
| 322 | input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
|
|---|
| 323 | input_cs.history()[internal2input[he.parent2]].jetp_index,
|
|---|
| 324 | he.dij, internal_cs.jets()[he.jetp_index], new_index);
|
|---|
| 325 | internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
|
|---|
| 326 | // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
|
|---|
| 327 | // << " and " << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
|
|---|
| 328 | // << " into " << internal2input[i] << " (int: " << i << ")" << endl;
|
|---|
| 329 | }
|
|---|
| 330 | }
|
|---|
| 331 | }
|
|---|
| 332 |
|
|---|
| 333 | // returns the plugin description
|
|---|
| 334 | string PruningPlugin::description() const{
|
|---|
| 335 | ostringstream oss;
|
|---|
| 336 | oss << "Pruning plugin with jet_definition = (" << _jet_def.description()
|
|---|
| 337 | <<"), zcut = " << _zcut
|
|---|
| 338 | << ", Rcut = " << _Rcut;
|
|---|
| 339 | return oss.str();
|
|---|
| 340 | }
|
|---|
| 341 |
|
|---|
| 342 |
|
|---|
| 343 | FASTJET_END_NAMESPACE
|
|---|
