Loading...
Searching...
No Matches
PseudoJet.cc
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet 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.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet 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 FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
#include "fastjet/Error.hh"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
#ifndef __FJCORE__
#include "fastjet/ClusterSequenceAreaBase.hh"
#endif // __FJCORE__
#include "fastjet/CompositeJetStructure.hh"
#include<valarray>
#include<iostream>
#include<sstream>
#include<cmath>
#include<algorithm>
#include <cstdarg>
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
using namespace std;
//----------------------------------------------------------------------
// another constructor...
PseudoJet::PseudoJet(const double px_in, const double py_in, const double pz_in, const double E_in) {
_E = E_in ;
_px = px_in;
_py = py_in;
_pz = pz_in;
this->_finish_init();
// some default values for the history and user indices
// note: no reset of shared pointers needed
_reset_indices();
}
#ifdef FASTJET_HAVE_THREAD_SAFETY
/// copy-assignmemt
///
/// this has to be explicitly specified since atomic does not support it.
PseudoJet & PseudoJet::operator=(const PseudoJet & other_pj){
_structure = other_pj._structure;
_user_info = other_pj._user_info;
_kt2 = other_pj._kt2;
_cluster_hist_index = other_pj._cluster_hist_index;
_user_index = other_pj._user_index;
_px = other_pj._px;
_py = other_pj._py;
_pz = other_pj._pz;
_E = other_pj._E;
_phi = other_pj._phi;
_rap = other_pj._rap;
_init_status.store(other_pj._init_status);
return *this;
}
#endif // FASTJET_HAVE_THREAD_SAFETY
//std::shared_ptr: //----------------------------------------------------------------------
//std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
//std::shared_ptr: PseudoJet::~PseudoJet(){
//std::shared_ptr: _release_jet_from_cs();
//std::shared_ptr: }
//std::shared_ptr:
//std::shared_ptr: // this has to be called everytime one tries to alter the jet
//std::shared_ptr: // structural info
//std::shared_ptr: void PseudoJet::_release_jet_from_cs(){
//std::shared_ptr: // check if the jet has the structure of type CSstruct in which case
//std::shared_ptr: // we have to check if there is a need for self-deletion of the CS
//std::shared_ptr: ClusterSequenceStructure * assoc_css = dynamic_cast<ClusterSequenceStructure *>(_structure.get());
//std::shared_ptr: if (assoc_css) {
//std::shared_ptr: const ClusterSequence * assoc_cs = assoc_css->associated_cluster_sequence();
//std::shared_ptr: if (assoc_cs) assoc_cs->release_pseudojet(*this);
//std::shared_ptr: }
//std::shared_ptr:
//std::shared_ptr: // slightly less efficient version
//std::shared_ptr: // if ((has_structure_of<ClusterSequence>()) && (has_valid_cluster_sequence())){
//std::shared_ptr: // associated_cs()->release_pseudojet(*this);
//std::shared_ptr: // }
//std::shared_ptr: }
//std::shared_ptr:
//std::shared_ptr: #endif // FASTJET_HAVE_THREAD_SAFETY
//----------------------------------------------------------------------
/// do standard end of initialisation
void PseudoJet::_finish_init () {
_kt2 = this->px()*this->px() + this->py()*this->py();
_phi = pseudojet_invalid_phi;
// strictly speaking, _rap does not need initialising, because
// it's never used as long as _phi == pseudojet_invalid_phi
// (and gets set when _phi is requested). However ATLAS
// 2013-03-28 complained that they sometimes have NaN's in
// _rap and this interferes with some of their internal validation.
// So we initialise it; penalty is about 0.3ns per PseudoJet out of about
// 10ns total initialisation time (on a intel Core i7 2.7GHz)
_rap = pseudojet_invalid_rap;
#ifdef FASTJET_HAVE_THREAD_SAFETY
_init_status = Init_NotDone;
#endif
}
//----------------------------------------------------------------------
#ifdef FASTJET_HAVE_THREAD_SAFETY
void PseudoJet::_ensure_valid_rap_phi() const{
//TODO: for _init_status we can use memory_order_release when
//writing and memory_order_acquire when reading. Check if that has
//an impact on timing
if (_init_status!=Init_Done){
int expected = Init_NotDone;
// if it's 0 (not done), set thing to "operation in progress"
// (-1) and do the init
//
// Note:
// - this has to be the strong version so we can make a single
// test
// - we need to make sure that the -1 is loaded properly
// (success), hence a std::memory_order_seq_cst model
// - the failure model can be anything since we're not going
// to use the value anyway
if (_init_status.compare_exchange_strong(expected, Init_InProgress,
std::memory_order_seq_cst,
std::memory_order_relaxed)){
_set_rap_phi();
_init_status = Init_Done; // can safely be done after all physics varlables are set
} else {
// wait until done
do{
// the operation below will reset expected to whatever is in
// init_state if the test fails, so we need to reset it to
// the 1 (aka init_done) we want!
expected = Init_Done;
// the next line
// - here we could potentially use the weak form
// - on success the value is unchanged so I think we can use relaxed ordering
// - expected will be reinitialised anyway so again, relaxed ordering should be fi
//} while (!_init_state.compare_exchange_strong(expected, 1));
} while (!_init_status.compare_exchange_weak(expected, Init_Done,
std::memory_order_relaxed,
std::memory_order_relaxed));
}
}
}
#endif
//----------------------------------------------------------------------
void PseudoJet::_set_rap_phi() const {
if (_kt2 == 0.0) {
_phi = 0.0; }
else {
_phi = atan2(this->py(),this->px());
}
if (_phi < 0.0) {_phi += twopi;}
if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|?
if (this->E() == abs(this->pz()) && _kt2 == 0) {
// Point has infinite rapidity -- convert that into a very large
// number, but in such a way that different 0-pt momenta will have
// different rapidities (so as to lift the degeneracy between
// them) [this can be relevant at parton-level]
double MaxRapHere = MaxRap + abs(this->pz());
if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;}
} else {
// get the rapidity in a way that's modestly insensitive to roundoff
// error when things pz,E are large (actually the best we can do without
// explicit knowledge of mass)
double effective_m2 = max(0.0,m2()); // force non tachyonic mass
double E_plus_pz = _E + abs(_pz); // the safer of p+, p-
// p+/p- = (p+ p-) / (p-)^2 = (kt^2+m^2)/(p-)^2
_rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
if (_pz > 0) {_rap = - _rap;}
}
}
//----------------------------------------------------------------------
// return a valarray four-momentum
valarray<double> PseudoJet::four_mom() const {
valarray<double> mom(4);
mom[0] = _px;
mom[1] = _py;
mom[2] = _pz;
mom[3] = _E ;
return mom;
}
//----------------------------------------------------------------------
// Return the component corresponding to the specified index.
// taken from CLHEP
double PseudoJet::operator () (int i) const {
switch(i) {
case X:
return px();
case Y:
return py();
case Z:
return pz();
case T:
return e();
default:
ostringstream err;
err << "PseudoJet subscripting: bad index (" << i << ")";
throw Error(err.str());
}
return 0.;
}
//----------------------------------------------------------------------
// return the pseudorapidity
double PseudoJet::pseudorapidity() const {
if (px() == 0.0 && py() ==0.0) return MaxRap;
if (pz() == 0.0) return 0.0;
double theta = atan(perp()/pz());
if (theta < 0) theta += pi;
return -log(tan(theta/2));
}
//----------------------------------------------------------------------
// return "sum" of two pseudojets
PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) {
//return PseudoJet(jet1.four_mom()+jet2.four_mom());
return PseudoJet(jet1.px()+jet2.px(),
jet1.py()+jet2.py(),
jet1.pz()+jet2.pz(),
jet1.E() +jet2.E() );
}
//----------------------------------------------------------------------
// return difference of two pseudojets
PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) {
//return PseudoJet(jet1.four_mom()-jet2.four_mom());
return PseudoJet(jet1.px()-jet2.px(),
jet1.py()-jet2.py(),
jet1.pz()-jet2.pz(),
jet1.E() -jet2.E() );
}
//----------------------------------------------------------------------
// return the product, coeff * jet
PseudoJet operator* (double coeff, const PseudoJet & jet) {
// see the comment in operator*= about ensuring valid rap phi
// before a multiplication to handle case of multiplication by
// zero, while maintaining rapidity and phi
jet._ensure_valid_rap_phi();
//return PseudoJet(coeff*jet.four_mom());
// the following code is hopefully more efficient
PseudoJet coeff_times_jet(jet);
coeff_times_jet *= coeff;
return coeff_times_jet;
}
//----------------------------------------------------------------------
// return the product, coeff * jet
PseudoJet operator* (const PseudoJet & jet, double coeff) {
return coeff*jet;
}
//----------------------------------------------------------------------
// return the ratio, jet / coeff
PseudoJet operator/ (const PseudoJet & jet, double coeff) {
return (1.0/coeff)*jet;
}
//----------------------------------------------------------------------
/// multiply the jet's momentum by the coefficient
PseudoJet & PseudoJet::operator*=(double coeff) {
// operator*= aims to maintain the rapidity and azimuth
// for the PseudoJet; if they have already been evaluated
// this is fine, but if they haven't and coeff is sufficiently
// small as to cause a zero or underflow result, then a subsequent
// invocation of rap or phi will lead to a non-sensical result.
// So, here, we preemptively ensure that rapidity and phi
// are correctly cached
_ensure_valid_rap_phi();
_px *= coeff;
_py *= coeff;
_pz *= coeff;
_E *= coeff;
_kt2*= coeff*coeff;
// phi and rap are unchanged
return *this;
}
//----------------------------------------------------------------------
/// divide the jet's momentum by the coefficient
PseudoJet & PseudoJet::operator/=(double coeff) {
(*this) *= 1.0/coeff;
return *this;
}
//----------------------------------------------------------------------
/// add the other jet's momentum to this jet
PseudoJet & PseudoJet::operator+=(const PseudoJet & other_jet) {
_px += other_jet._px;
_py += other_jet._py;
_pz += other_jet._pz;
_E += other_jet._E ;
_finish_init(); // we need to recalculate phi,rap,kt2
return *this;
}
//----------------------------------------------------------------------
/// subtract the other jet's momentum from this jet
PseudoJet & PseudoJet::operator-=(const PseudoJet & other_jet) {
_px -= other_jet._px;
_py -= other_jet._py;
_pz -= other_jet._pz;
_E -= other_jet._E ;
_finish_init(); // we need to recalculate phi,rap,kt2
return *this;
}
//----------------------------------------------------------------------
bool operator==(const PseudoJet & a, const PseudoJet & b) {
if (a.px() != b.px()) return false;
if (a.py() != b.py()) return false;
if (a.pz() != b.pz()) return false;
if (a.E () != b.E ()) return false;
if (a.user_index() != b.user_index()) return false;
if (a.cluster_hist_index() != b.cluster_hist_index()) return false;
if (a.user_info_ptr() != b.user_info_ptr()) return false;
if (a.structure_ptr() != b.structure_ptr()) return false;
return true;
}
//----------------------------------------------------------------------
// check if the jet has zero momentum
bool operator==(const PseudoJet & jet, const double val) {
if (val != 0)
throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed.");
return (jet.px() == 0 && jet.py() == 0 &&
jet.pz() == 0 && jet.E() == 0);
}
//----------------------------------------------------------------------
/// transform this jet (given in the rest frame of prest) into a jet
/// in the lab frame
//
// NB: code adapted from that in herwig f77 (checked how it worked
// long ago)
PseudoJet & PseudoJet::boost(const PseudoJet & prest) {
if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
return *this;
double m_local = prest.m();
assert(m_local != 0);
double pf4 = ( px()*prest.px() + py()*prest.py()
+ pz()*prest.pz() + E()*prest.E() )/m_local;
double fn = (pf4 + E()) / (prest.E() + m_local);
_px += fn*prest.px();
_py += fn*prest.py();
_pz += fn*prest.pz();
_E = pf4;
_finish_init(); // we need to recalculate phi,rap,kt2
return *this;
}
//----------------------------------------------------------------------
/// transform this jet (given in lab) into a jet in the rest
/// frame of prest
//
// NB: code adapted from that in herwig f77 (checked how it worked
// long ago)
PseudoJet & PseudoJet::unboost(const PseudoJet & prest) {
if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
return *this;
double m_local = prest.m();
assert(m_local != 0);
double pf4 = ( -px()*prest.px() - py()*prest.py()
- pz()*prest.pz() + E()*prest.E() )/m_local;
double fn = (pf4 + E()) / (prest.E() + m_local);
_px -= fn*prest.px();
_py -= fn*prest.py();
_pz -= fn*prest.pz();
_E = pf4;
_finish_init(); // we need to recalculate phi,rap,kt2
return *this;
}
//----------------------------------------------------------------------
/// returns true if the momenta of the two input jets are identical
bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) {
return jeta.px() == jetb.px()
&& jeta.py() == jetb.py()
&& jeta.pz() == jetb.pz()
&& jeta.E() == jetb.E();
}
//----------------------------------------------------------------------
void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) {
_rap = rap_in; _phi = phi_in;
if (_phi >= twopi) _phi -= twopi;
if (_phi < 0) _phi += twopi;
#ifdef FASTJET_HAVE_THREAD_SAFETY
_init_status = Init_Done;
#endif
}
//----------------------------------------------------------------------
void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) {
assert(phi_in < 2*twopi && phi_in > -twopi);
double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
double exprap = exp(y_in);
double pminus = ptm/exprap;
double pplus = ptm*exprap;
double px_local = pt_in*cos(phi_in);
double py_local = pt_in*sin(phi_in);
reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
set_cached_rap_phi(y_in,phi_in);
}
//----------------------------------------------------------------------
/// return a pseudojet with the given pt, y, phi and mass
PseudoJet PtYPhiM(double pt, double y, double phi, double m) {
assert(phi < 2*twopi && phi > -twopi);
double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
double exprap = exp(y);
double pminus = ptm/exprap;
double pplus = ptm*exprap;
double px = pt*cos(phi);
double py = pt*sin(phi);
PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
mom.set_cached_rap_phi(y,phi);
return mom;
//return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y));
}
//----------------------------------------------------------------------
// return kt-distance between this jet and another one
double PseudoJet::kt_distance(const PseudoJet & other) const {
//double distance = min(this->kt2(), other.kt2());
double distance = min(_kt2, other._kt2);
double dphi = abs(phi() - other.phi());
if (dphi > pi) {dphi = twopi - dphi;}
double drap = rap() - other.rap();
distance = distance * (dphi*dphi + drap*drap);
return distance;
}
//----------------------------------------------------------------------
// return squared cylinder (eta-phi) distance between this jet and another one
double PseudoJet::plain_distance(const PseudoJet & other) const {
double dphi = abs(phi() - other.phi());
if (dphi > pi) {dphi = twopi - dphi;}
double drap = rap() - other.rap();
return (dphi*dphi + drap*drap);
}
//----------------------------------------------------------------------
/// returns other.phi() - this.phi(), i.e. the phi distance to
/// other, constrained to be in range -pi .. pi
double PseudoJet::delta_phi_to(const PseudoJet & other) const {
double dphi = other.phi() - phi();
if (dphi > pi) dphi -= twopi;
if (dphi < -pi) dphi += twopi;
return dphi;
}
string PseudoJet::description() const{
// the "default" case of a PJ which does not belong to any cluster sequence
if (!_structure)
return "standard PseudoJet (with no associated clustering information)";
// for all the other cases, the description comes from the structure
return _structure->description();
}
//----------------------------------------------------------------------
//
// The following methods access the associated jet structure (if any)
//
//----------------------------------------------------------------------
//----------------------------------------------------------------------
// check whether this PseudoJet has an associated parent
// ClusterSequence
bool PseudoJet::has_associated_cluster_sequence() const{
return (_structure) && (_structure->has_associated_cluster_sequence());
}
//----------------------------------------------------------------------
// get a (const) pointer to the associated ClusterSequence (NULL if
// inexistent)
const ClusterSequence* PseudoJet::associated_cluster_sequence() const{
if (! has_associated_cluster_sequence()) return NULL;
return _structure->associated_cluster_sequence();
}
//----------------------------------------------------------------------
// check whether this PseudoJet has an associated parent
// ClusterSequence that is still valid
bool PseudoJet::has_valid_cluster_sequence() const{
return (_structure) && (_structure->has_valid_cluster_sequence());
}
//----------------------------------------------------------------------
// If there is a valid cluster sequence associated with this jet,
// returns a pointer to it; otherwise throws an Error.
//
// Open question: should these errors be upgraded to classes of their
// own so that they can be caught? [Maybe, but later]
const ClusterSequence * PseudoJet::validated_cs() const {
return validated_structure_ptr()->validated_cs();
}
//----------------------------------------------------------------------
// set the associated structure
void PseudoJet::set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in){
//std::shared_ptr: #ifdef FASTJET_HAVE_THREAD_SAFETY
//std::shared_ptr: // if the jet currently belongs to a cs, we need to release it before any chenge
//std::shared_ptr: _release_jet_from_cs();
//std::shared_ptr: #endif //FASTJET_HAVE_THREAD_SAFETY
_structure = structure_in;
}
//----------------------------------------------------------------------
// return true if there is some structure associated with this PseudoJet
bool PseudoJet::has_structure() const{
return (bool) _structure; // cast to bool has to be made explicit
}
//----------------------------------------------------------------------
// return a pointer to the structure (of type
// PseudoJetStructureBase*) associated with this PseudoJet.
//
// return NULL if there is no associated structure
const PseudoJetStructureBase* PseudoJet::structure_ptr() const {
//if (!_structure) return NULL;
return _structure.get();
}
//----------------------------------------------------------------------
// return a non-const pointer to the structure (of type
// PseudoJetStructureBase*) associated with this PseudoJet.
//
// return NULL if there is no associated structure
//
// Only use this if you know what you are doing. In any case,
// prefer the 'structure_ptr()' (the const version) to this method,
// unless you really need a write access to the PseudoJet's
// underlying structure.
PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){
//if (!_structure) return NULL;
return _structure.get();
}
//----------------------------------------------------------------------
// return a pointer to the structure (of type
// PseudoJetStructureBase*) associated with this PseudoJet.
//
// throw an error if there is no associated structure
const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const {
if (!_structure)
throw Error("Trying to access the structure of a PseudoJet which has no associated structure");
return _structure.get();
}
//----------------------------------------------------------------------
// return a reference to the shared pointer to the
// PseudoJetStructureBase associated with this PseudoJet
const SharedPtr<PseudoJetStructureBase> & PseudoJet::structure_shared_ptr() const {
return _structure;
}
//----------------------------------------------------------------------
// check if it has been recombined with another PseudoJet in which
// case, return its partner through the argument. Otherwise,
// 'partner' is set to 0.
//
// false is also returned if this PseudoJet has no associated
// ClusterSequence
bool PseudoJet::has_partner(PseudoJet &partner) const{
return validated_structure_ptr()->has_partner(*this, partner);
}
//----------------------------------------------------------------------
// check if it has been recombined with another PseudoJet in which
// case, return its child through the argument. Otherwise, 'child'
// is set to 0.
//
// false is also returned if this PseudoJet has no associated
// ClusterSequence, with the child set to 0
bool PseudoJet::has_child(PseudoJet &child) const{
return validated_structure_ptr()->has_child(*this, child);
}
//----------------------------------------------------------------------
// check if it is the product of a recombination, in which case
// return the 2 parents through the 'parent1' and 'parent2'
// arguments. Otherwise, set these to 0.
//
// false is also returned if this PseudoJet has no parent
// ClusterSequence
bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{
return validated_structure_ptr()->has_parents(*this, parent1, parent2);
}
//----------------------------------------------------------------------
// check if the current PseudoJet contains the one passed as
// argument
//
// false is also returned if this PseudoJet has no associated
// ClusterSequence.
bool PseudoJet::contains(const PseudoJet &constituent) const{
return validated_structure_ptr()->object_in_jet(constituent, *this);
}
//----------------------------------------------------------------------
// check if the current PseudoJet is contained the one passed as
// argument
//
// false is also returned if this PseudoJet has no associated
// ClusterSequence
bool PseudoJet::is_inside(const PseudoJet &jet) const{
return validated_structure_ptr()->object_in_jet(*this, jet);
}
//----------------------------------------------------------------------
// returns true if the PseudoJet has constituents
bool PseudoJet::has_constituents() const{
return (_structure) && (_structure->has_constituents());
}
//----------------------------------------------------------------------
// retrieve the constituents.
vector<PseudoJet> PseudoJet::constituents() const{
return validated_structure_ptr()->constituents(*this);
}
//----------------------------------------------------------------------
// returns true if the PseudoJet has support for exclusive subjets
bool PseudoJet::has_exclusive_subjets() const{
return (_structure) && (_structure->has_exclusive_subjets());
}
//----------------------------------------------------------------------
// return a vector of all subjets of the current jet (in the sense
// of the exclusive algorithm) that would be obtained when running
// the algorithm with the given dcut.
//
// Time taken is O(m ln m), where m is the number of subjets that
// are found. If m gets to be of order of the total number of
// constituents in the jet, this could be substantially slower than
// just getting that list of constituents.
//
// an Error is thrown if this PseudoJet has no currently valid
// associated ClusterSequence
std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double dcut) const {
return validated_structure_ptr()->exclusive_subjets(*this, dcut);
}
//----------------------------------------------------------------------
// return the size of exclusive_subjets(...); still n ln n with same
// coefficient, but marginally more efficient than manually taking
// exclusive_subjets.size()
//
// an Error is thrown if this PseudoJet has no currently valid
// associated ClusterSequence
int PseudoJet::n_exclusive_subjets(const double dcut) const {
return validated_structure_ptr()->n_exclusive_subjets(*this, dcut);
}
//----------------------------------------------------------------------
// return the list of subjets obtained by unclustering the supplied
// jet down to n subjets (or all constituents if there are fewer
// than n).
//
// requires n ln n time
//
// an Error is thrown if this PseudoJet has no currently valid
// associated ClusterSequence
std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const {
return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub);
}
//----------------------------------------------------------------------
// Same as exclusive_subjets_up_to but throws an error if there are
// fewer than nsub particles in the jet
std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const {
vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
if (int(subjets.size()) < nsub) {
ostringstream err;
err << "Requested " << nsub << " exclusive subjets, but there were only "
<< subjets.size() << " particles in the jet";
throw Error(err.str());
}
return subjets;
}
//----------------------------------------------------------------------
// return the dij that was present in the merging nsub+1 -> nsub
// subjets inside this jet.
//
// an Error is thrown if this PseudoJet has no currently valid
// associated ClusterSequence
double PseudoJet::exclusive_subdmerge(int nsub) const {
return validated_structure_ptr()->exclusive_subdmerge(*this, nsub);
}
//----------------------------------------------------------------------
// return the maximum dij that occurred in the whole event at the
// stage that the nsub+1 -> nsub merge of subjets occurred inside
// this jet.
//
// an Error is thrown if this PseudoJet has no currently valid
// associated ClusterSequence
double PseudoJet::exclusive_subdmerge_max(int nsub) const {
return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub);
}
// returns true if a jet has pieces
//
// By default a single particle or a jet coming from a
// ClusterSequence have no pieces and this methos will return false.
bool PseudoJet::has_pieces() const{
return ((_structure) && (_structure->has_pieces(*this)));
}
// retrieve the pieces that make up the jet.
//
// By default a jet does not have pieces.
// If the underlying interface supports "pieces" retrieve the
// pieces from there.
std::vector<PseudoJet> PseudoJet::pieces() const{
return validated_structure_ptr()->pieces(*this);
// if (!has_pieces())
// throw Error("Trying to retrieve the pieces of a PseudoJet that has no support for pieces.");
//
// return _structure->pieces(*this);
}
//----------------------------------------------------------------------
// the following ones require a computation of the area in the
// associated ClusterSequence (See ClusterSequenceAreaBase for details)
//----------------------------------------------------------------------
#ifndef __FJCORE__
//----------------------------------------------------------------------
// if possible, return a valid ClusterSequenceAreaBase pointer; otherwise
// throw an error
const ClusterSequenceAreaBase * PseudoJet::validated_csab() const {
const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase*>(validated_cs());
if (csab == NULL) throw Error("you requested jet-area related information, but the PseudoJet does not have associated area information.");
return csab;
}
//----------------------------------------------------------------------
// check if it has a defined area
bool PseudoJet::has_area() const{
//if (! has_associated_cluster_sequence()) return false;
if (! has_structure()) return false;
return (validated_structure_ptr()->has_area() != 0);
}
//----------------------------------------------------------------------
// return the jet (scalar) area.
// throw an Error if there is no support for area in the associated CS
double PseudoJet::area() const{
return validated_structure_ptr()->area(*this);
}
//----------------------------------------------------------------------
// return the error (uncertainty) associated with the determination
// of the area of this jet.
// throws an Error if there is no support for area in the associated CS
double PseudoJet::area_error() const{
return validated_structure_ptr()->area_error(*this);
}
//----------------------------------------------------------------------
// return the jet 4-vector area
// throws an Error if there is no support for area in the associated CS
PseudoJet PseudoJet::area_4vector() const{
return validated_structure_ptr()->area_4vector(*this);
}
//----------------------------------------------------------------------
// true if this jet is made exclusively of ghosts
// throws an Error if there is no support for area in the associated CS
bool PseudoJet::is_pure_ghost() const{
return validated_structure_ptr()->is_pure_ghost(*this);
}
#endif // __FJCORE__
//----------------------------------------------------------------------
//
// end of the methods accessing the information in the associated
// Cluster Sequence
//
//----------------------------------------------------------------------
//----------------------------------------------------------------------
/// provide a meaningful error message for InexistentUserInfo
PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
{}
//----------------------------------------------------------------------
// sort the indices so that values[indices[0..n-1]] is sorted
// into increasing order
void sort_indices(vector<int> & indices,
const vector<double> & values) {
IndexedSortHelper index_sort_helper(&values);
sort(indices.begin(), indices.end(), index_sort_helper);
}
//----------------------------------------------------------------------
/// return a vector of jets sorted into decreasing kt2
vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) {
vector<double> minus_kt2(jets.size());
for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
return objects_sorted_by_values(jets, minus_kt2);
}
//----------------------------------------------------------------------
/// return a vector of jets sorted into increasing rapidity
vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) {
vector<double> rapidities(jets.size());
for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
return objects_sorted_by_values(jets, rapidities);
}
//----------------------------------------------------------------------
/// return a vector of jets sorted into decreasing energy
vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) {
vector<double> energies(jets.size());
for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
return objects_sorted_by_values(jets, energies);
}
//----------------------------------------------------------------------
/// return a vector of jets sorted into increasing pz
vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) {
vector<double> pz(jets.size());
for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
return objects_sorted_by_values(jets, pz);
}
//-------------------------------------------------------------------------------
// helper functions to build a jet made of pieces
//-------------------------------------------------------------------------------
// build a "CompositeJet" from the vector of its pieces
//
// In this case, E-scheme recombination is assumed to compute the
// total momentum
PseudoJet join(const vector<PseudoJet> & pieces){
// compute the total momentum
//--------------------------------------------------
PseudoJet result; // automatically initialised to 0
for (unsigned int i=0; i<pieces.size(); i++)
result += pieces[i];
// attach a CompositeJetStructure to the result
//--------------------------------------------------
CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces);
result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
return result;
}
// build a "CompositeJet" from a single PseudoJet
PseudoJet join(const PseudoJet & j1){
return join(vector<PseudoJet>(1,j1));
}
// build a "CompositeJet" from two PseudoJet
PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){
vector<PseudoJet> pieces;
pieces.reserve(2);
pieces.push_back(j1);
pieces.push_back(j2);
return join(pieces);
}
// build a "CompositeJet" from 3 PseudoJet
PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){
vector<PseudoJet> pieces;
pieces.reserve(3);
pieces.push_back(j1);
pieces.push_back(j2);
pieces.push_back(j3);
return join(pieces);
}
// build a "CompositeJet" from 4 PseudoJet
PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){
vector<PseudoJet> pieces;
pieces.reserve(4);
pieces.push_back(j1);
pieces.push_back(j2);
pieces.push_back(j3);
pieces.push_back(j4);
return join(pieces);
}
FASTJET_END_NAMESPACE
fastjet::ClusterSequenceAreaBase
base class that sets interface for extensions of ClusterSequence that provide information about the a...
Definition ClusterSequenceAreaBase.hh:48
fastjet::ClusterSequence
deals with clustering
Definition ClusterSequence.hh:63
fastjet::Error
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
fastjet::Error::Error
Error()
default constructors
Definition Error.hh:55
fastjet::IndexedSortHelper
Definition PseudoJet.hh:1014
fastjet::PseudoJetStructureBase
Contains any information related to the clustering that should be directly accessible to PseudoJet.
Definition PseudoJetStructureBase.hh:59
fastjet::PseudoJet::InexistentUserInfo::InexistentUserInfo
InexistentUserInfo()
provide a meaningful error message for InexistentUserInfo
Definition PseudoJet.cc:866
fastjet::PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
fastjet::PseudoJet::structure_non_const_ptr
PseudoJetStructureBase * structure_non_const_ptr()
return a non-const pointer to the structure (of type PseudoJetStructureBase*) associated with this Ps...
Definition PseudoJet.cc:606
fastjet::PseudoJet::boost
PseudoJet & boost(const PseudoJet &prest)
transform this jet (given in the rest frame of prest) into a jet in the lab frame
Definition PseudoJet.cc:387
fastjet::PseudoJet::PseudoJet
PseudoJet()
default constructor, which as of FJ3.0 provides an object for which all operations are now valid and ...
Definition PseudoJet.hh:82
fastjet::PseudoJet::exclusive_subjets_up_to
std::vector< PseudoJet > exclusive_subjets_up_to(int nsub) const
return the list of subjets obtained by unclustering the supplied jet down to nsub subjets (or all con...
Definition PseudoJet.cc:739
fastjet::PseudoJet::has_child
virtual bool has_child(PseudoJet &child) const
check if it has been recombined with another PseudoJet in which case, return its child through the ar...
Definition PseudoJet.cc:648
fastjet::PseudoJet::has_structure
bool has_structure() const
return true if there is some structure associated with this PseudoJet
Definition PseudoJet.cc:582
fastjet::PseudoJet::reset_momentum_PtYPhiM
void reset_momentum_PtYPhiM(double pt, double y, double phi, double m=0.0)
reset the 4-momentum according to the specified pt, rapidity, azimuth and mass (phi should satisfy -2...
Definition PseudoJet.cc:455
fastjet::PseudoJet::operator*=
PseudoJet & operator*=(double)
multiply the jet's momentum by the coefficient
Definition PseudoJet.cc:306
fastjet::PseudoJet::exclusive_subdmerge
double exclusive_subdmerge(int nsub) const
Returns the dij that was present in the merging nsub+1 -> nsub subjets inside this jet.
Definition PseudoJet.cc:763
fastjet::PseudoJet::constituents
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition PseudoJet.cc:692
fastjet::PseudoJet::set_cached_rap_phi
void set_cached_rap_phi(double rap, double phi)
in some cases when setting a 4-momentum, the user/program knows what rapidity and azimuth are associa...
Definition PseudoJet.cc:445
fastjet::PseudoJet::n_exclusive_subjets
int n_exclusive_subjets(const double dcut) const
return the size of exclusive_subjets(...); still n ln n with same coefficient, but marginally more ef...
Definition PseudoJet.cc:726
fastjet::PseudoJet::rap
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition PseudoJet.hh:138
fastjet::PseudoJet::is_inside
virtual bool is_inside(const PseudoJet &jet) const
check if the current PseudoJet is contained the one passed as argument.
Definition PseudoJet.cc:679
fastjet::PseudoJet::contains
virtual bool contains(const PseudoJet &constituent) const
check if the current PseudoJet contains the one passed as argument.
Definition PseudoJet.cc:669
fastjet::PseudoJet::plain_distance
double plain_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
Definition PseudoJet.cc:499
fastjet::PseudoJet::user_info_ptr
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
Definition PseudoJet.hh:496
fastjet::PseudoJet::operator-=
PseudoJet & operator-=(const PseudoJet &)
subtract the other jet's momentum from this jet
Definition PseudoJet.cc:346
fastjet::PseudoJet::reset_momentum
void reset_momentum(double px, double py, double pz, double E)
reset the 4-momentum according to the supplied components but leave all other information (indices,...
Definition PseudoJet.hh:1076
fastjet::PseudoJet::has_partner
virtual bool has_partner(PseudoJet &partner) const
check if it has been recombined with another PseudoJet in which case, return its partner through the ...
Definition PseudoJet.cc:637
fastjet::PseudoJet::is_pure_ghost
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
Definition PseudoJet.cc:851
fastjet::PseudoJet::validated_csab
const ClusterSequenceAreaBase * validated_csab() const
shorthand for validated_cluster_sequence_area_base()
Definition PseudoJet.cc:811
fastjet::PseudoJet::phi
double phi() const
returns phi (in the range 0..2pi)
Definition PseudoJet.hh:123
fastjet::PseudoJet::area_error
virtual double area_error() const
return the error (uncertainty) associated with the determination of the area of this jet.
Definition PseudoJet.cc:837
fastjet::PseudoJet::four_mom
std::valarray< double > four_mom() const
return a valarray containing the four-momentum (components 0-2 are 3-mom, component 3 is energy).
Definition PseudoJet.cc:217
fastjet::PseudoJet::has_area
virtual bool has_area() const
check if it has a defined area
Definition PseudoJet.cc:820
fastjet::PseudoJet::perp
double perp() const
returns the scalar transverse momentum
Definition PseudoJet.hh:158
fastjet::PseudoJet::exclusive_subdmerge_max
double exclusive_subdmerge_max(int nsub) const
Returns the maximum dij that occurred in the whole event at the stage that the nsub+1 -> nsub merge o...
Definition PseudoJet.cc:774
fastjet::PseudoJet::structure_ptr
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet.
Definition PseudoJet.cc:591
fastjet::PseudoJet::operator+=
PseudoJet & operator+=(const PseudoJet &)
add the other jet's momentum to this jet
Definition PseudoJet.cc:334
fastjet::PseudoJet::kt_distance
double kt_distance(const PseudoJet &other) const
returns kt distance (R=1) between this jet and another
Definition PseudoJet.cc:486
fastjet::PseudoJet::has_parents
virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const
check if it is the product of a recombination, in which case return the 2 parents through the 'parent...
Definition PseudoJet.cc:659
fastjet::PseudoJet::theta
double theta() const
polar angle
Definition PseudoJet.hh:193
fastjet::PseudoJet::has_valid_cluster_sequence
bool has_valid_cluster_sequence() const
returns true if this PseudoJet has an associated and still valid(ated) ClusterSequence.
Definition PseudoJet.cc:555
fastjet::PseudoJet::has_exclusive_subjets
virtual bool has_exclusive_subjets() const
returns true if the PseudoJet has support for exclusive subjets
Definition PseudoJet.cc:699
fastjet::PseudoJet::area
virtual double area() const
return the jet (scalar) area.
Definition PseudoJet.cc:829
fastjet::PseudoJet::has_pieces
virtual bool has_pieces() const
returns true if a jet has pieces
Definition PseudoJet.cc:783
fastjet::PseudoJet::cluster_hist_index
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
Definition PseudoJet.hh:830
fastjet::PseudoJet::has_constituents
virtual bool has_constituents() const
returns true if the PseudoJet has constituents
Definition PseudoJet.cc:686
fastjet::PseudoJet::set_structure_shared_ptr
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure_in)
set the associated structure
Definition PseudoJet.cc:572
fastjet::PseudoJet::validated_structure_ptr
const PseudoJetStructureBase * validated_structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet.
Definition PseudoJet.cc:616
fastjet::PseudoJet::exclusive_subjets
std::vector< PseudoJet > exclusive_subjets(const double dcut) const
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that woul...
Definition PseudoJet.cc:715
fastjet::PseudoJet::m
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
Definition PseudoJet.hh:1060
fastjet::PseudoJet::user_index
int user_index() const
return the user_index,
Definition PseudoJet.hh:389
fastjet::PseudoJet::description
std::string description() const
return a string describing what kind of PseudoJet we are dealing with
Definition PseudoJet.cc:517
fastjet::PseudoJet::associated_cluster_sequence
const ClusterSequence * associated_cluster_sequence() const
get a (const) pointer to the parent ClusterSequence (NULL if inexistent)
Definition PseudoJet.cc:545
fastjet::PseudoJet::operator/=
PseudoJet & operator/=(double)
divide the jet's momentum by the coefficient
Definition PseudoJet.cc:326
fastjet::PseudoJet::validated_cs
const ClusterSequence * validated_cs() const
shorthand for validated_cluster_sequence()
Definition PseudoJet.cc:565
fastjet::PseudoJet::pieces
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
Definition PseudoJet.cc:792
fastjet::PseudoJet::area_4vector
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition PseudoJet.cc:844
fastjet::PseudoJet::delta_phi_to
double delta_phi_to(const PseudoJet &other) const
returns other.phi() - this.phi(), constrained to be in range -pi .
Definition PseudoJet.cc:509
fastjet::PseudoJet::structure_shared_ptr
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
return a reference to the shared pointer to the PseudoJetStructureBase associated with this PseudoJet
Definition PseudoJet.cc:625
fastjet::PseudoJet::has_associated_cluster_sequence
bool has_associated_cluster_sequence() const
returns true if this PseudoJet has an associated ClusterSequence.
Definition PseudoJet.cc:538
fastjet::PseudoJet::pseudorapidity
double pseudorapidity() const
returns the pseudo-rapidity or some large value when the rapidity is infinite
Definition PseudoJet.cc:249
fastjet::PseudoJet::unboost
PseudoJet & unboost(const PseudoJet &prest)
transform this jet (given in lab) into a jet in the rest frame of prest
Definition PseudoJet.cc:414
fastjet::SharedPtr
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition SharedPtr.hh:341
fastjet::sorted_by_rapidity
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
Definition PseudoJet.cc:890
fastjet::MaxRap
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons,...
Definition PseudoJet.hh:53
fastjet::sorted_by_E
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
Definition PseudoJet.cc:898
fastjet::have_same_momentum
bool have_same_momentum(const PseudoJet &jeta, const PseudoJet &jetb)
returns true if the momenta of the two input jets are identical
Definition PseudoJet.cc:437
fastjet::operator==
bool operator==(const PseudoJet &a, const PseudoJet &b)
returns true if the 4 momentum components of the two PseudoJets are identical and all the internal in...
Definition PseudoJet.cc:356
fastjet::objects_sorted_by_values
std::vector< T > objects_sorted_by_values(const std::vector< T > &objects, const std::vector< double > &values)
given a vector of values with a one-to-one correspondence with the vector of objects,...
Definition PseudoJet.hh:984
fastjet::sorted_by_pt
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition PseudoJet.cc:882
fastjet::PtYPhiM
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition PseudoJet.cc:469
fastjet::sorted_by_pz
vector< PseudoJet > sorted_by_pz(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing pz
Definition PseudoJet.cc:906