Page MenuHomeHEPForge
Files F8309817

Recluster.cc
No OneTemporary
Actions

Recluster.cc
View Options

// $Id: Recluster.cc 682 2014-06-13 08:59:36Z gsoyez $
//
// Copyright (c) 2014-, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// 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 "Recluster.hh"
#include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
#include <sstream>
#include <typeinfo>
using namespace std;
// Comments:
//
// - If the jet comes from a C/A clustering (or is a composite jet
// made of C/A clusterings) and we recluster it with a C/A
// algorithm, we just use exclusive jets instead of doing the
// clustering explicitly. In this specific case, we need to make
// sure that all the pieces share the same cluster sequence.
//
// - If the recombiner has to be taken from the original jet and that
// jet is composite, we need to check that all the pieces were
// obtained with the same recombiner.
//
// TODO:
//
// - check this actually works!!!
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib{
LimitedWarning Recluster::_explicit_ghost_warning;
// class description
string Recluster::description() const {
ostringstream ostr;
ostr << "Recluster with subjet_def = ";
if (_use_full_def) {
ostr << _subjet_def.description();
} else {
if (_subjet_alg == kt_algorithm) {
ostr << "Longitudinally invariant kt algorithm with R = " << _subjet_radius;
} else if (_subjet_alg == cambridge_algorithm) {
ostr << "Longitudinally invariant Cambridge/Aachen algorithm with R = " << _subjet_radius;
} else if (_subjet_alg == antikt_algorithm) {
ostr << "Longitudinally invariant anti-kt algorithm with R = " << _subjet_radius;
} else if (_subjet_alg == genkt_algorithm) {
ostr << "Longitudinally invariant generalised kt algorithm with R = " << _subjet_radius
<< ", p = " << _subjet_extra;
} else if (_subjet_alg == cambridge_for_passive_algorithm) {
ostr << "Longitudinally invariant Cambridge/Aachen algorithm with R = " << _subjet_radius
<< " and a special hack whereby particles with kt < "
<< _subjet_extra << "are treated as passive ghosts";
} else if (_subjet_alg == ee_kt_algorithm) {
ostr << "e+e- kt (Durham) algorithm";
} else if (_subjet_alg == ee_genkt_algorithm) {
ostr << "e+e- generalised kt algorithm with R = " << _subjet_radius
<< ", p = " << _subjet_extra;
} else if (_subjet_alg == undefined_jet_algorithm) {
ostr << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
} else {
ostr << "unrecognized jet_algorithm";
}
ostr << ", a recombiner obtained from the jet being reclustered";
}
if (_single)
ostr << " and keeping the hardest subjet";
else
ostr << " and joining all subjets in a composite jet";
return ostr.str();
}
// return a vector of subjets, which are the ones that would be kept
// by the filtering
PseudoJet Recluster::result(const PseudoJet &jet) const {
// generic sanity checks
//-------------------------------------------------------------------
// make sure that the jet has constituents
if (! jet.has_constituents())
throw Error("Filter can only be applied on jets having constituents");
// tests particular to certain configurations
//-------------------------------------------------------------------
// for a whole variety of checks, we shall need the "recursive"
// pieces of the jet (the jet itself or recursing down to its most
// fundamental pieces). So we do compute these once and for all.
//
// Note that the pieces are always needed (either for C/A or for the
// area checks)
vector<PseudoJet> all_pieces; //.clear();
if ((!_get_all_pieces(jet, all_pieces)) || (all_pieces.size()==0)){
throw Error("Recluster: failed to retrieve all the pieces composing the jet.");
}
// decide which jet definition to use
//-------------------------------------------------------------------
JetDefinition subjet_def;
if (_use_full_def){
subjet_def = _subjet_def;
} else {
_build_jet_def_with_recombiner(all_pieces, subjet_def);
}
// the vector that will ultimately hold the subjets
vector<PseudoJet> subjets;
// check if we can apply the simplification for C/A jets reclustered
// with C/A
//
// we apply C/A clustering iff
// - the request subjet_def is C/A
// - the jet is either directly coming from C/A or if it is a
// superposition of C/A jets from the same cluster sequence
// - the pieces agree with the recombination scheme of subjet_def
//
// Note that in this case area support will be automatically
// inherted so we can only worry about this later
//-------------------------------------------------------------------
if (_check_ca(all_pieces, subjet_def)){
_recluster_cafilt(all_pieces, subjets, subjet_def.R());
subjets = sorted_by_pt(subjets);
return _single
? subjets[0]
: join(subjets, *(subjet_def.recombiner()));
}
// decide if area support has to be kept
//-------------------------------------------------------------------
bool include_area_support = jet.has_area();
if ((include_area_support) && (!_check_explicit_ghosts(all_pieces))){
_explicit_ghost_warning.warn("Recluster: the original cluster sequence is lacking explicit ghosts; area support will no longer be available after re-clustering");
include_area_support = false;
}
// extract the subjets
//-------------------------------------------------------------------
_recluster_generic(jet, subjets, subjet_def, include_area_support);
subjets = sorted_by_pt(subjets);
return _single
? subjets[0]
: join(subjets, *(subjet_def.recombiner()));
}
//----------------------------------------------------------------------
// the parts that really do the reclustering
//----------------------------------------------------------------------
// get the subjets in the simple case of C/A+C/A
void Recluster::_recluster_cafilt(const vector<PseudoJet> & all_pieces,
vector<PseudoJet> & subjets,
double Rfilt) const{
subjets.clear();
// each individual piece should have a C/A cluster sequence
// associated to it. So a simple loop would do the job
for (vector<PseudoJet>::const_iterator piece_it = all_pieces.begin();
piece_it!=all_pieces.end(); piece_it++){
// just extract the exclusive subjets of 'jet'
const ClusterSequence *cs = piece_it->associated_cluster_sequence();
vector<PseudoJet> local_subjets;
double dcut = Rfilt / cs->jet_def().R();
if (dcut>=1.0){
local_subjets.push_back(*piece_it);
} else {
local_subjets = piece_it->exclusive_subjets(dcut*dcut);
}
copy(local_subjets.begin(), local_subjets.end(), back_inserter(subjets));
}
}
// set the filtered elements in the generic re-clustering case (wo
// subtraction)
void Recluster::_recluster_generic(const PseudoJet & jet,
vector<PseudoJet> & subjets,
const JetDefinition & subjet_def,
bool do_areas) const{
// create a new, internal, ClusterSequence from the jet constituents
// get the subjets directly from there
//
// If the jet has area support then we separate the ghosts from the
// "regular" particles so the subjets will also have area
// support. Note that we do this regardless of whether rho is zero
// or not.
//
// Note that to be able to separate the ghosts, one needs explicit
// ghosts!!
// ---------------------------------------------------------------
if (do_areas){
vector<PseudoJet> all_constituents = jet.constituents();
vector<PseudoJet> regular_constituents, ghosts;
for (vector<PseudoJet>::iterator it = all_constituents.begin();
it != all_constituents.end(); it++){
if (it->is_pure_ghost())
ghosts.push_back(*it);
else
regular_constituents.push_back(*it);
}
// figure out the ghost area from the 1st ghost (if none, any value
// would probably do as the area will be 0 and subtraction will have
// no effect!)
double ghost_area = (ghosts.size()) ? ghosts[0].area() : 0.01;
ClusterSequenceActiveAreaExplicitGhosts * csa
= new ClusterSequenceActiveAreaExplicitGhosts(regular_constituents,
subjet_def,
ghosts, ghost_area);
subjets = csa->inclusive_jets();
// allow the cs to be deleted when it's no longer used
//
// Note that there is at least one constituent in the jet so there
// is in principle at least one subjet But one may have used a
// nasty recombiner that left an empty set of subjets, so we'd
// rather play it safe
if (subjets.size())
csa->delete_self_when_unused();
else
delete csa;
} else {
ClusterSequence * cs = new ClusterSequence(jet.constituents(), subjet_def);
subjets = cs->inclusive_jets();
// allow the cs to be deleted when it's no longer used (again, we
// add an extra safety check)
if (subjets.size())
cs->delete_self_when_unused();
else
delete cs;
}
}
//----------------------------------------------------------------------
// various checks and internal constructs
//----------------------------------------------------------------------
// fundamental info for CompositeJets
//----------------------------------------------------------------------
// get the pieces down to the fundamental pieces
//
// Note that this just checks that there is an associated CS to the
// fundamental pieces, not that it is still valid
bool Recluster::_get_all_pieces(const PseudoJet &jet, vector<PseudoJet> &all_pieces) const{
if (jet.has_associated_cluster_sequence()){
all_pieces.push_back(jet);
return true;
}
if (jet.has_pieces()){
const vector<PseudoJet> pieces = jet.pieces();
for (vector<PseudoJet>::const_iterator it=pieces.begin(); it!=pieces.end(); it++)
if (!_get_all_pieces(*it, all_pieces)) return false;
return true;
}
return false;
}
// treatment of recombiners
//----------------------------------------------------------------------
// get the common recombiner to all pieces (NULL if none)
//
// Note that if the jet has an associated cluster sequence that is no
// longer valid, an error will be thrown (needed since it could be the
// 1st check called after the enumeration of the pieces)
const JetDefinition::Recombiner* Recluster::_get_common_recombiner(const vector<PseudoJet> &all_pieces) const{
const JetDefinition & jd_ref = all_pieces[0].validated_cs()->jet_def();
for (unsigned int i=1; i<all_pieces.size(); i++)
if (!all_pieces[i].validated_cs()->jet_def().has_same_recombiner(jd_ref)) return NULL;
return jd_ref.recombiner();
}
void Recluster::_build_jet_def_with_recombiner(const vector<PseudoJet> &all_pieces,
JetDefinition &subjet_def) const{
// the recombiner has to be guessed from the pieces
const JetDefinition::Recombiner * common_recombiner = _get_common_recombiner(all_pieces);
if (common_recombiner) {
if (typeid(*common_recombiner) == typeid(JetDefinition::DefaultRecombiner)) {
RecombinationScheme scheme =
static_cast<const JetDefinition::DefaultRecombiner *>(common_recombiner)->scheme();
if (_has_subjet_extra)
subjet_def = JetDefinition(_subjet_alg, _subjet_radius, _subjet_extra, scheme);
else if (_has_subjet_radius)
subjet_def = JetDefinition(_subjet_alg, _subjet_radius, scheme);
else
subjet_def = JetDefinition(_subjet_alg, scheme);
} else {
if (_has_subjet_extra)
subjet_def = JetDefinition(_subjet_alg, _subjet_radius, _subjet_extra, common_recombiner);
else if (_has_subjet_radius)
subjet_def = JetDefinition(_subjet_alg, _subjet_radius, common_recombiner);
else
subjet_def = JetDefinition(_subjet_alg, common_recombiner);
}
} else {
throw Error("Recluster: requested to guess the recombination scheme (or recombiner) from the original jet but an inconsistency was found between the pieces constituing that jet.");
}
}
// area support
//----------------------------------------------------------------------
// check if the jet (or all its pieces) have explicit ghosts
// (assuming the jet has area support).
//
// Note that if the jet has an associated cluster sequence that is no
// longer valid, an error will be thrown (needed since it could be the
// 1st check called after the enumeration of the pieces)
bool Recluster::_check_explicit_ghosts(const vector<PseudoJet> &all_pieces) const{
for (vector<PseudoJet>::const_iterator it=all_pieces.begin(); it!=all_pieces.end(); it++)
if (! it->validated_csab()->has_explicit_ghosts()) return false;
return true;
}
// C/A specific tests
//----------------------------------------------------------------------
// check if one can apply the simplification for C/A subjets
//
// This includes:
// - the subjet definition asks for C/A subjets
// - all the pieces share the same CS
// - that CS is C/A with the same recombiner as the subjet def
// - the filtering radius is not larger than any of the pairwise
// distance between the pieces
//
// Note that if the jet has an associated cluster sequence that is no
// longer valid, an error will be thrown (needed since it could be the
// 1st check called after the enumeration of the pieces)
bool Recluster::_check_ca(const vector<PseudoJet> &all_pieces,
const JetDefinition &subjet_def) const{
if (subjet_def.jet_algorithm() != cambridge_algorithm) return false;
// check that the 1st of all the pieces (we're sure there is at
// least one) is coming from a C/A clustering. Then check that all
// the following pieces share the same ClusterSequence
const ClusterSequence * cs_ref = all_pieces[0].validated_cs();
if (cs_ref->jet_def().jet_algorithm() != cambridge_algorithm) return false;
for (unsigned int i=1; i<all_pieces.size(); i++)
if (all_pieces[i].validated_cs() != cs_ref) return false;
// check that the 1st peice has the same recombiner as the one used
// for the subjet clustering
// Note that since they share the same CS, checking the 1st one is enough
if (!cs_ref->jet_def().has_same_recombiner(subjet_def)) return false;
// we also have to make sure that the filtering radius is not larger
// than any of the inter-pieces distance
double Rfilt2 = subjet_def.R();
Rfilt2 *= Rfilt2;
for (unsigned int i=0; i<all_pieces.size()-1; i++){
for (unsigned int j=i+1; j<all_pieces.size(); j++){
if (all_pieces[i].squared_distance(all_pieces[j]) < Rfilt2) return false;
}
}
return true;
}
} // contrib namespace
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 4:37 PM (21 h, 39 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023475
Default Alt Text
Recluster.cc (15 KB)

Event Timeline