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