| 1 | // $Id: BottomUpSoftDrop.cc 1064 2017-09-08 09:19:57Z gsoyez $
|
|---|
| 2 | //
|
|---|
| 3 | // Copyright (c) 2017-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
|
|---|
| 4 | // Kevin Zhou, Frederic Dreyer
|
|---|
| 5 | //
|
|---|
| 6 | //----------------------------------------------------------------------
|
|---|
| 7 | // This file is part of FastJet contrib.
|
|---|
| 8 | //
|
|---|
| 9 | // It is free software; you can redistribute it and/or modify it under
|
|---|
| 10 | // the terms of the GNU General Public License as published by the
|
|---|
| 11 | // Free Software Foundation; either version 2 of the License, or (at
|
|---|
| 12 | // your option) any later version.
|
|---|
| 13 | //
|
|---|
| 14 | // It is distributed in the hope that it will be useful, but WITHOUT
|
|---|
| 15 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|---|
| 16 | // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
|
|---|
| 17 | // License for more details.
|
|---|
| 18 | //
|
|---|
| 19 | // You should have received a copy of the GNU General Public License
|
|---|
| 20 | // along with this code. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 21 | //----------------------------------------------------------------------
|
|---|
| 22 |
|
|---|
| 23 | #include "BottomUpSoftDrop.hh"
|
|---|
| 24 | #include <cassert>
|
|---|
| 25 | #include <algorithm>
|
|---|
| 26 | #include <sstream>
|
|---|
| 27 | #include <typeinfo>
|
|---|
| 28 | #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
|
|---|
| 29 | #include "fastjet/Selector.hh"
|
|---|
| 30 | #include "fastjet/config.h"
|
|---|
| 31 |
|
|---|
| 32 |
|
|---|
| 33 | using namespace std;
|
|---|
| 34 |
|
|---|
| 35 |
|
|---|
| 36 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 37 |
|
|---|
| 38 | namespace contrib{
|
|---|
| 39 |
|
|---|
| 40 | //----------------------------------------------------------------------
|
|---|
| 41 | // BottomUpSoftDrop class
|
|---|
| 42 | //----------------------------------------------------------------------
|
|---|
| 43 |
|
|---|
| 44 | // action on a single jet
|
|---|
| 45 | PseudoJet BottomUpSoftDrop::result(const PseudoJet &jet) const{
|
|---|
| 46 | // soft drop can only be applied to jets that have constituents
|
|---|
| 47 | if (!jet.has_constituents()){
|
|---|
| 48 | throw Error("BottomUpSoftDrop: trying to apply the Soft Drop transformer to a jet that has no constituents");
|
|---|
| 49 | }
|
|---|
| 50 |
|
|---|
| 51 | // if the jet has area support and there are explicit ghosts, we can
|
|---|
| 52 | // transfer that support to the internal re-clustering
|
|---|
| 53 | //
|
|---|
| 54 | // Note that this is just meant to maintain the information since
|
|---|
| 55 | // all the jes will have a 0 area
|
|---|
| 56 | bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
|
|---|
| 57 |
|
|---|
| 58 | // build the soft drop plugin
|
|---|
| 59 | BottomUpSoftDropPlugin * softdrop_plugin;
|
|---|
| 60 |
|
|---|
| 61 | // for some constructors, we get the recombiner from the
|
|---|
| 62 | // input jet -- some acrobatics are needed
|
|---|
| 63 | if (_get_recombiner_from_jet) {
|
|---|
| 64 | JetDefinition jet_def = _jet_def;
|
|---|
| 65 |
|
|---|
| 66 | // if all the pieces have a shared recombiner, we'll use that
|
|---|
| 67 | // one. Otherwise, use the one from _jet_def as a fallback.
|
|---|
| 68 | JetDefinition jet_def_for_recombiner;
|
|---|
| 69 | if (_check_common_recombiner(jet, jet_def_for_recombiner)){
|
|---|
| 70 | #if FASTJET_VERSION_NUMBER >= 30100
|
|---|
| 71 | // Note that this is better than the option directly passing the
|
|---|
| 72 | // recombiner (for cases where th ejet def own its recombiner)
|
|---|
| 73 | // but it's only available in FJ>=3.1
|
|---|
| 74 | jet_def.set_recombiner(jet_def_for_recombiner);
|
|---|
| 75 | #else
|
|---|
| 76 | jet_def.set_recombiner(jet_def_for_recombiner.recombiner());
|
|---|
| 77 | #endif
|
|---|
| 78 | }
|
|---|
| 79 | softdrop_plugin = new BottomUpSoftDropPlugin(jet_def, _beta, _symmetry_cut, _R0);
|
|---|
| 80 | } else {
|
|---|
| 81 | softdrop_plugin = new BottomUpSoftDropPlugin(_jet_def, _beta, _symmetry_cut, _R0);
|
|---|
| 82 | }
|
|---|
| 83 |
|
|---|
| 84 | // now recluster the constituents of the jet with that plugin
|
|---|
| 85 | JetDefinition internal_jet_def(softdrop_plugin);
|
|---|
| 86 | // flag the plugin for automatic deletion _before_ we make
|
|---|
| 87 | // copies (so that as long as the copies are also present
|
|---|
| 88 | // it doesn't get deleted).
|
|---|
| 89 | internal_jet_def.delete_plugin_when_unused();
|
|---|
| 90 |
|
|---|
| 91 | ClusterSequence * cs;
|
|---|
| 92 | if (do_areas){
|
|---|
| 93 | vector<PseudoJet> particles, ghosts;
|
|---|
| 94 | SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
|
|---|
| 95 | // determine the ghost area from the 1st ghost (if none, any value
|
|---|
| 96 | // will do, as the area will be 0 and subtraction will have
|
|---|
| 97 | // no effect!)
|
|---|
| 98 | double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
|
|---|
| 99 | cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def,
|
|---|
| 100 | ghosts, ghost_area);
|
|---|
| 101 | } else {
|
|---|
| 102 | cs = new ClusterSequence(jet.constituents(), internal_jet_def);
|
|---|
| 103 | }
|
|---|
| 104 |
|
|---|
| 105 | PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
|
|---|
| 106 | BottomUpSoftDropStructure * s = new BottomUpSoftDropStructure(result_local);
|
|---|
| 107 | s->_beta = _beta;
|
|---|
| 108 | s->_symmetry_cut = _symmetry_cut;
|
|---|
| 109 | s->_R0 = _R0;
|
|---|
| 110 | result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
|
|---|
| 111 |
|
|---|
| 112 | // make sure things remain persistent -- i.e. tell the jet definition
|
|---|
| 113 | // and the cluster sequence that it is their responsibility to clean
|
|---|
| 114 | // up memory once the "result" reaches the end of its life in the user's
|
|---|
| 115 | // code. (The CS deletes itself when the result goes out of scope and
|
|---|
| 116 | // that also triggers deletion of the plugin)
|
|---|
| 117 | cs->delete_self_when_unused();
|
|---|
| 118 |
|
|---|
| 119 | return result_local;
|
|---|
| 120 | }
|
|---|
| 121 |
|
|---|
| 122 | // global grooming on a full event
|
|---|
| 123 | // note: does not support jet areas
|
|---|
| 124 | vector<PseudoJet> BottomUpSoftDrop::global_grooming(const vector<PseudoJet> & event) const {
|
|---|
| 125 | // start by reclustering the event into one very large jet
|
|---|
| 126 | ClusterSequence cs(event, _jet_def);
|
|---|
| 127 | std::vector<PseudoJet> global_jet = SelectorNHardest(1)(cs.inclusive_jets());
|
|---|
| 128 | // if the event is empty, do nothing
|
|---|
| 129 | if (global_jet.size() == 0) return vector<PseudoJet>();
|
|---|
| 130 |
|
|---|
| 131 | // apply the groomer to the large jet
|
|---|
| 132 | PseudoJet result = this->result(global_jet[0]);
|
|---|
| 133 | return result.constituents();
|
|---|
| 134 | }
|
|---|
| 135 |
|
|---|
| 136 | // check if the jet has explicit_ghosts (knowing that there is an
|
|---|
| 137 | // area support)
|
|---|
| 138 | bool BottomUpSoftDrop::_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 true and set jet_def_for_recombiner so that the
|
|---|
| 159 | // recombiner can be taken from that JetDefinition. Otherwise,
|
|---|
| 160 | // return false. 'assigned' is initially false; when true, each
|
|---|
| 161 | // time we meet a new jet definition, we'll check it shares the
|
|---|
| 162 | // same recombiner as jet_def_for_recombiner.
|
|---|
| 163 | bool BottomUpSoftDrop::_check_common_recombiner(const PseudoJet &jet,
|
|---|
| 164 | JetDefinition &jet_def_for_recombiner,
|
|---|
| 165 | bool assigned) const {
|
|---|
| 166 | if (jet.has_associated_cluster_sequence()){
|
|---|
| 167 | // if the jet def for recombination has already been assigned, check if we have the same
|
|---|
| 168 | if (assigned)
|
|---|
| 169 | return jet.validated_cs()->jet_def().has_same_recombiner(jet_def_for_recombiner);
|
|---|
| 170 |
|
|---|
| 171 | // otherwise, assign it.
|
|---|
| 172 | jet_def_for_recombiner = jet.validated_cs()->jet_def();
|
|---|
| 173 | assigned = true;
|
|---|
| 174 | return true;
|
|---|
| 175 | }
|
|---|
| 176 |
|
|---|
| 177 | // if the jet has pieces, recurse in the pieces
|
|---|
| 178 | if (jet.has_pieces()){
|
|---|
| 179 | vector<PseudoJet> pieces = jet.pieces();
|
|---|
| 180 | if (pieces.size() == 0) return false;
|
|---|
| 181 | for (unsigned int i=0;i<pieces.size(); i++)
|
|---|
| 182 | if (!_check_common_recombiner(pieces[i], jet_def_for_recombiner, assigned)) return false;
|
|---|
| 183 | // never returned false, so we're OK.
|
|---|
| 184 | return true;
|
|---|
| 185 | }
|
|---|
| 186 | // return false for any other (unknown) structure
|
|---|
| 187 | return false;
|
|---|
| 188 | }
|
|---|
| 189 |
|
|---|
| 190 |
|
|---|
| 191 | // description
|
|---|
| 192 | string BottomUpSoftDrop::description() const{
|
|---|
| 193 | ostringstream oss;
|
|---|
| 194 | oss << "BottomUpSoftDrop with jet_definition = (" << _jet_def.description() << ")"
|
|---|
| 195 | << ", symmetry_cut = " << _symmetry_cut << ", beta = " << _beta << ", R0 = " << _R0;
|
|---|
| 196 | return oss.str();
|
|---|
| 197 | }
|
|---|
| 198 |
|
|---|
| 199 | //----------------------------------------------------------------------
|
|---|
| 200 | // BottomUpSoftDropRecombiner class
|
|---|
| 201 | //----------------------------------------------------------------------
|
|---|
| 202 |
|
|---|
| 203 | // perform a recombination taking into account the Soft Drop
|
|---|
| 204 | // conditions
|
|---|
| 205 | void BottomUpSoftDropRecombiner::recombine(const PseudoJet &pa,
|
|---|
| 206 | const PseudoJet &pb,
|
|---|
| 207 | PseudoJet &pab) const {
|
|---|
| 208 | PseudoJet p;
|
|---|
| 209 | _recombiner->recombine(pa, pb, p);
|
|---|
| 210 |
|
|---|
| 211 | // calculate symmetry parameter
|
|---|
| 212 | double symmetry_cut_fn =_symmetry_cut * pow(pa.squared_distance(pb)/_R0sqr, 0.5*_beta);
|
|---|
| 213 | double pta = pa.pt();
|
|---|
| 214 | double ptb = pb.pt();
|
|---|
| 215 | // make sure denominator is non-zero
|
|---|
| 216 | double sym = pta + ptb;
|
|---|
| 217 | if (sym == 0){
|
|---|
| 218 | pab = p;
|
|---|
| 219 | return;
|
|---|
| 220 | }
|
|---|
| 221 | sym = min(pta, ptb) / sym;
|
|---|
| 222 |
|
|---|
| 223 | // if the 2 particles pass the soft drop condition, do the
|
|---|
| 224 | // recombination
|
|---|
| 225 | if (sym > symmetry_cut_fn){
|
|---|
| 226 | pab=p;
|
|---|
| 227 | return;
|
|---|
| 228 | }
|
|---|
| 229 |
|
|---|
| 230 | // if the soft drop condition is not passed, return the particle
|
|---|
| 231 | // with largest pt
|
|---|
| 232 | if (pta < ptb) {
|
|---|
| 233 | pab = pb;
|
|---|
| 234 | _rejected.push_back(pa.cluster_hist_index());
|
|---|
| 235 | } else {
|
|---|
| 236 | pab = pa;
|
|---|
| 237 | _rejected.push_back(pb.cluster_hist_index());
|
|---|
| 238 | }
|
|---|
| 239 | }
|
|---|
| 240 |
|
|---|
| 241 | //----------------------------------------------------------------------
|
|---|
| 242 | // BottomUpSoftDropPlugin class
|
|---|
| 243 | //----------------------------------------------------------------------
|
|---|
| 244 |
|
|---|
| 245 | // the actual clustering work for the plugin
|
|---|
| 246 | void BottomUpSoftDropPlugin::run_clustering(ClusterSequence &input_cs) const {
|
|---|
| 247 | // declare a soft drop recombiner
|
|---|
| 248 | BottomUpSoftDropRecombiner softdrop_recombiner(_beta, _symmetry_cut, _R0, _jet_def.recombiner());
|
|---|
| 249 | JetDefinition jet_def = _jet_def;
|
|---|
| 250 | jet_def.set_recombiner(&softdrop_recombiner);
|
|---|
| 251 |
|
|---|
| 252 | // cluster the particles using that recombiner
|
|---|
| 253 | ClusterSequence internal_cs(input_cs.jets(), jet_def);
|
|---|
| 254 | const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
|
|---|
| 255 |
|
|---|
| 256 | // transfer the list of "childless" elements into a bool vector
|
|---|
| 257 | vector<bool> kept(internal_hist.size(), true);
|
|---|
| 258 | const vector<unsigned int> &sd_rej = softdrop_recombiner.rejected();
|
|---|
| 259 | for (unsigned int i=0;i<sd_rej.size(); i++) kept[sd_rej[i]]=false;
|
|---|
| 260 |
|
|---|
| 261 | // browse the history, keeping only the elements that have not been
|
|---|
| 262 | // vetoed.
|
|---|
| 263 | //
|
|---|
| 264 | // In the process we build a map for the history indices
|
|---|
| 265 | vector<unsigned int> internal2input(internal_hist.size());
|
|---|
| 266 | for (unsigned int i=0; i<input_cs.jets().size(); i++)
|
|---|
| 267 | internal2input[i] = i;
|
|---|
| 268 |
|
|---|
| 269 | for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
|
|---|
| 270 | const ClusterSequence::history_element &he = internal_hist[i];
|
|---|
| 271 |
|
|---|
| 272 | // deal with recombinations with the beam
|
|---|
| 273 | if (he.parent2 == ClusterSequence::BeamJet){
|
|---|
| 274 | int internal_jetp_index = internal_hist[he.parent1].jetp_index;
|
|---|
| 275 | int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
|
|---|
| 276 |
|
|---|
| 277 | int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
|
|---|
| 278 |
|
|---|
| 279 | // cout << "Beam recomb for internal " << internal_hist_index
|
|---|
| 280 | // << " (input jet index=" << input_jetp_index << endl;
|
|---|
| 281 |
|
|---|
| 282 | input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
|
|---|
| 283 | continue;
|
|---|
| 284 | }
|
|---|
| 285 |
|
|---|
| 286 | // now, deal with two-body recombinations
|
|---|
| 287 | if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
|
|---|
| 288 | internal2input[i]=internal2input[he.parent2];
|
|---|
| 289 | // cout << "rejecting internal " << he.parent1
|
|---|
| 290 | // << ", mapping internal " << i
|
|---|
| 291 | // << " to internal " << he.parent2
|
|---|
| 292 | // << " i.e. " << internal2input[i] << endl;
|
|---|
| 293 | } else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
|
|---|
| 294 | internal2input[i]=internal2input[he.parent1];
|
|---|
| 295 | // cout << "rejecting internal " << he.parent2
|
|---|
| 296 | // << ", mapping internal " << i
|
|---|
| 297 | // << " to internal " << he.parent1
|
|---|
| 298 | // << " i.e. " << internal2input[i] << endl;
|
|---|
| 299 | } else { // do the recombination
|
|---|
| 300 | int new_index;
|
|---|
| 301 | input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
|
|---|
| 302 | input_cs.history()[internal2input[he.parent2]].jetp_index,
|
|---|
| 303 | he.dij, internal_cs.jets()[he.jetp_index], new_index);
|
|---|
| 304 | internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
|
|---|
| 305 | // cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
|
|---|
| 306 | // << " and " << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
|
|---|
| 307 | // << " into " << internal2input[i] << " (int: " << i << ")" << endl;
|
|---|
| 308 | }
|
|---|
| 309 | }
|
|---|
| 310 | }
|
|---|
| 311 |
|
|---|
| 312 | // description of the plugin
|
|---|
| 313 | string BottomUpSoftDropPlugin::description() const {
|
|---|
| 314 | ostringstream oss;
|
|---|
| 315 | oss << "BottomUpSoftDropPlugin with jet_definition = (" << _jet_def.description()
|
|---|
| 316 | <<"), symmetry_cut = " << _symmetry_cut << ", beta = "
|
|---|
| 317 | << _beta << ", R0 = " << _R0;
|
|---|
| 318 | return oss.str();
|
|---|
| 319 | }
|
|---|
| 320 |
|
|---|
| 321 |
|
|---|
| 322 | }
|
|---|
| 323 |
|
|---|
| 324 | FASTJET_END_NAMESPACE
|
|---|
