Page MenuHomeHEPForge
Files F8309823

BottomUpSoftDrop.cc
No OneTemporary
Actions

BottomUpSoftDrop.cc
View Options

// $Id: BottomUpSoftDrop.cc 1064 2017-09-08 09:19:57Z gsoyez $
//
// Copyright (c) 2017-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
// Kevin Zhou, Frederic Dreyer
//
//----------------------------------------------------------------------
// This file is part of FastJet contrib.
//
// It is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at
// your option) any later version.
//
// It is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this code. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#include "BottomUpSoftDrop.hh"
#include <cassert>
#include <algorithm>
#include <sstream>
#include <typeinfo>
#include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh"
#include "fastjet/Selector.hh"
#include "fastjet/config.h"
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib{
//----------------------------------------------------------------------
// BottomUpSoftDrop class
//----------------------------------------------------------------------
// action on a single jet
PseudoJet BottomUpSoftDrop::result(const PseudoJet &jet) const{
// soft drop can only be applied to jets that have constituents
if (!jet.has_constituents()){
throw Error("BottomUpSoftDrop: trying to apply the Soft Drop transformer to a jet that has no constituents");
}
// if the jet has area support and there are explicit ghosts, we can
// transfer that support to the internal re-clustering
//
// Note that this is just meant to maintain the information since
// all the jes will have a 0 area
bool do_areas = jet.has_area() && _check_explicit_ghosts(jet);
// build the soft drop plugin
BottomUpSoftDropPlugin * softdrop_plugin;
// for some constructors, we get the recombiner from the
// input jet -- some acrobatics are needed
if (_get_recombiner_from_jet) {
JetDefinition jet_def = _jet_def;
// if all the pieces have a shared recombiner, we'll use that
// one. Otherwise, use the one from _jet_def as a fallback.
JetDefinition jet_def_for_recombiner;
if (_check_common_recombiner(jet, jet_def_for_recombiner)){
#if FASTJET_VERSION_NUMBER >= 30100
// Note that this is better than the option directly passing the
// recombiner (for cases where th ejet def own its recombiner)
// but it's only available in FJ>=3.1
jet_def.set_recombiner(jet_def_for_recombiner);
#else
jet_def.set_recombiner(jet_def_for_recombiner.recombiner());
#endif
}
softdrop_plugin = new BottomUpSoftDropPlugin(jet_def, _beta, _symmetry_cut, _R0);
} else {
softdrop_plugin = new BottomUpSoftDropPlugin(_jet_def, _beta, _symmetry_cut, _R0);
}
// now recluster the constituents of the jet with that plugin
JetDefinition internal_jet_def(softdrop_plugin);
// flag the plugin for automatic deletion _before_ we make
// copies (so that as long as the copies are also present
// it doesn't get deleted).
internal_jet_def.delete_plugin_when_unused();
ClusterSequence * cs;
if (do_areas){
vector<PseudoJet> particles, ghosts;
SelectorIsPureGhost().sift(jet.constituents(), ghosts, particles);
// determine the ghost area from the 1st ghost (if none, any value
// will do, as the area will be 0 and subtraction will have
// no effect!)
double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
cs = new ClusterSequenceActiveAreaExplicitGhosts(particles, internal_jet_def,
ghosts, ghost_area);
} else {
cs = new ClusterSequence(jet.constituents(), internal_jet_def);
}
PseudoJet result_local = SelectorNHardest(1)(cs->inclusive_jets())[0];
BottomUpSoftDropStructure * s = new BottomUpSoftDropStructure(result_local);
s->_beta = _beta;
s->_symmetry_cut = _symmetry_cut;
s->_R0 = _R0;
result_local.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(s));
// make sure things remain persistent -- i.e. tell the jet definition
// and the cluster sequence that it is their responsibility to clean
// up memory once the "result" reaches the end of its life in the user's
// code. (The CS deletes itself when the result goes out of scope and
// that also triggers deletion of the plugin)
cs->delete_self_when_unused();
return result_local;
}
// global grooming on a full event
// note: does not support jet areas
vector<PseudoJet> BottomUpSoftDrop::global_grooming(const vector<PseudoJet> & event) const {
// start by reclustering the event into one very large jet
ClusterSequence cs(event, _jet_def);
std::vector<PseudoJet> global_jet = SelectorNHardest(1)(cs.inclusive_jets());
// if the event is empty, do nothing
if (global_jet.size() == 0) return vector<PseudoJet>();
// apply the groomer to the large jet
PseudoJet result = this->result(global_jet[0]);
return result.constituents();
}
// check if the jet has explicit_ghosts (knowing that there is an
// area support)
bool BottomUpSoftDrop::_check_explicit_ghosts(const PseudoJet &jet) const {
// if the jet comes from a Clustering check explicit ghosts in that
// clustering
if (jet.has_associated_cluster_sequence())
return jet.validated_csab()->has_explicit_ghosts();
// if the jet has pieces, recurse in the pieces
if (jet.has_pieces()){
vector<PseudoJet> pieces = jet.pieces();
for (unsigned int i=0;i<pieces.size(); i++)
if (!_check_explicit_ghosts(pieces[i])) return false;
// never returned false, so we're OK.
return true;
}
// return false for any other (unknown) structure
return false;
}
// see if there is a common recombiner among the pieces; if there
// is return true and set jet_def_for_recombiner so that the
// recombiner can be taken from that JetDefinition. Otherwise,
// return false. 'assigned' is initially false; when true, each
// time we meet a new jet definition, we'll check it shares the
// same recombiner as jet_def_for_recombiner.
bool BottomUpSoftDrop::_check_common_recombiner(const PseudoJet &jet,
JetDefinition &jet_def_for_recombiner,
bool assigned) const {
if (jet.has_associated_cluster_sequence()){
// if the jet def for recombination has already been assigned, check if we have the same
if (assigned)
return jet.validated_cs()->jet_def().has_same_recombiner(jet_def_for_recombiner);
// otherwise, assign it.
jet_def_for_recombiner = jet.validated_cs()->jet_def();
assigned = true;
return true;
}
// if the jet has pieces, recurse in the pieces
if (jet.has_pieces()){
vector<PseudoJet> pieces = jet.pieces();
if (pieces.size() == 0) return false;
for (unsigned int i=0;i<pieces.size(); i++)
if (!_check_common_recombiner(pieces[i], jet_def_for_recombiner, assigned)) return false;
// never returned false, so we're OK.
return true;
}
// return false for any other (unknown) structure
return false;
}
// description
string BottomUpSoftDrop::description() const{
ostringstream oss;
oss << "BottomUpSoftDrop with jet_definition = (" << _jet_def.description() << ")"
<< ", symmetry_cut = " << _symmetry_cut << ", beta = " << _beta << ", R0 = " << _R0;
return oss.str();
}
//----------------------------------------------------------------------
// BottomUpSoftDropRecombiner class
//----------------------------------------------------------------------
// perform a recombination taking into account the Soft Drop
// conditions
void BottomUpSoftDropRecombiner::recombine(const PseudoJet &pa,
const PseudoJet &pb,
PseudoJet &pab) const {
PseudoJet p;
_recombiner->recombine(pa, pb, p);
// calculate symmetry parameter
double symmetry_cut_fn =_symmetry_cut * pow(pa.squared_distance(pb)/_R0sqr, 0.5*_beta);
double pta = pa.pt();
double ptb = pb.pt();
// make sure denominator is non-zero
double sym = pta + ptb;
if (sym == 0){
pab = p;
return;
}
sym = min(pta, ptb) / sym;
// if the 2 particles pass the soft drop condition, do the
// recombination
if (sym > symmetry_cut_fn){
pab=p;
return;
}
// if the soft drop condition is not passed, return the particle
// with largest pt
if (pta < ptb) {
pab = pb;
_rejected.push_back(pa.cluster_hist_index());
} else {
pab = pa;
_rejected.push_back(pb.cluster_hist_index());
}
}
//----------------------------------------------------------------------
// BottomUpSoftDropPlugin class
//----------------------------------------------------------------------
// the actual clustering work for the plugin
void BottomUpSoftDropPlugin::run_clustering(ClusterSequence &input_cs) const {
// declare a soft drop recombiner
BottomUpSoftDropRecombiner softdrop_recombiner(_beta, _symmetry_cut, _R0, _jet_def.recombiner());
JetDefinition jet_def = _jet_def;
jet_def.set_recombiner(&softdrop_recombiner);
// cluster the particles using that recombiner
ClusterSequence internal_cs(input_cs.jets(), jet_def);
const vector<ClusterSequence::history_element> & internal_hist = internal_cs.history();
// transfer the list of "childless" elements into a bool vector
vector<bool> kept(internal_hist.size(), true);
const vector<unsigned int> &sd_rej = softdrop_recombiner.rejected();
for (unsigned int i=0;i<sd_rej.size(); i++) kept[sd_rej[i]]=false;
// browse the history, keeping only the elements that have not been
// vetoed.
//
// In the process we build a map for the history indices
vector<unsigned int> internal2input(internal_hist.size());
for (unsigned int i=0; i<input_cs.jets().size(); i++)
internal2input[i] = i;
for (unsigned int i=input_cs.jets().size(); i<internal_hist.size(); i++){
const ClusterSequence::history_element &he = internal_hist[i];
// deal with recombinations with the beam
if (he.parent2 == ClusterSequence::BeamJet){
int internal_jetp_index = internal_hist[he.parent1].jetp_index;
int internal_hist_index = internal_cs.jets()[internal_jetp_index].cluster_hist_index();
int input_jetp_index = input_cs.history()[internal2input[internal_hist_index]].jetp_index;
// cout << "Beam recomb for internal " << internal_hist_index
// << " (input jet index=" << input_jetp_index << endl;
input_cs.plugin_record_iB_recombination(input_jetp_index, he.dij);
continue;
}
// now, deal with two-body recombinations
if (!kept[he.parent1]){ // 1 is rejected, we keep only 2
internal2input[i]=internal2input[he.parent2];
// cout << "rejecting internal " << he.parent1
// << ", mapping internal " << i
// << " to internal " << he.parent2
// << " i.e. " << internal2input[i] << endl;
} else if (!kept[he.parent2]){ // 2 is rejected, we keep only 1
internal2input[i]=internal2input[he.parent1];
// cout << "rejecting internal " << he.parent2
// << ", mapping internal " << i
// << " to internal " << he.parent1
// << " i.e. " << internal2input[i] << endl;
} else { // do the recombination
int new_index;
input_cs.plugin_record_ij_recombination(input_cs.history()[internal2input[he.parent1]].jetp_index,
input_cs.history()[internal2input[he.parent2]].jetp_index,
he.dij, internal_cs.jets()[he.jetp_index], new_index);
internal2input[i]=input_cs.jets()[new_index].cluster_hist_index();
// cout << "merging " << internal2input[he.parent1] << " (int: " << he.parent1 << ")"
// << " and " << internal2input[he.parent2] << " (int: " << he.parent2 << ")"
// << " into " << internal2input[i] << " (int: " << i << ")" << endl;
}
}
}
// description of the plugin
string BottomUpSoftDropPlugin::description() const {
ostringstream oss;
oss << "BottomUpSoftDropPlugin with jet_definition = (" << _jet_def.description()
<<"), symmetry_cut = " << _symmetry_cut << ", beta = "
<< _beta << ", R0 = " << _R0;
return oss.str();
}
}
FASTJET_END_NAMESPACE

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 4:38 PM (21 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022624
Default Alt Text
BottomUpSoftDrop.cc (12 KB)

Event Timeline