32#include "fastjet/Error.hh"
33#include "fastjet/PseudoJet.hh"
34#include "fastjet/ClusterSequence.hh"
36#include "fastjet/ClusterSequenceAreaBase.hh"
38#include "fastjet/CompositeJetStructure.hh"
46FASTJET_BEGIN_NAMESPACE
53PseudoJet::PseudoJet(
const double px_in,
const double py_in,
const double pz_in,
const double E_in) {
68#ifdef FASTJET_HAVE_THREAD_SAFETY
73 _structure = other_pj._structure;
74 _user_info = other_pj._user_info;
77 _cluster_hist_index = other_pj._cluster_hist_index;
78 _user_index = other_pj._user_index;
88 _init_status.store(other_pj._init_status);
121void PseudoJet::_finish_init () {
122 _kt2 = this->px()*this->px() + this->py()*this->py();
123 _phi = pseudojet_invalid_phi;
131 _rap = pseudojet_invalid_rap;
133#ifdef FASTJET_HAVE_THREAD_SAFETY
134 _init_status = Init_NotDone;
139#ifdef FASTJET_HAVE_THREAD_SAFETY
140void PseudoJet::_ensure_valid_rap_phi()
const{
144 if (_init_status!=Init_Done){
145 int expected = Init_NotDone;
156 if (_init_status.compare_exchange_strong(expected, Init_InProgress,
157 std::memory_order_seq_cst,
158 std::memory_order_relaxed)){
160 _init_status = Init_Done;
167 expected = Init_Done;
175 }
while (!_init_status.compare_exchange_weak(expected, Init_Done,
176 std::memory_order_relaxed,
177 std::memory_order_relaxed));
185void PseudoJet::_set_rap_phi()
const {
190 _phi = atan2(this->py(),this->px());
192 if (_phi < 0.0) {_phi += twopi;}
193 if (_phi >= twopi) {_phi -= twopi;}
194 if (this->E() == abs(this->pz()) && _kt2 == 0) {
199 double MaxRapHere =
MaxRap + abs(this->pz());
200 if (this->pz() >= 0.0) {_rap = MaxRapHere;}
else {_rap = -MaxRapHere;}
205 double effective_m2 = max(0.0,m2());
206 double E_plus_pz = _E + abs(_pz);
208 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
209 if (_pz > 0) {_rap = - _rap;}
217valarray<double> PseudoJet::four_mom()
const {
218 valarray<double> mom(4);
229double PseudoJet::operator () (
int i)
const {
241 err <<
"PseudoJet subscripting: bad index (" << i <<
")";
242 throw Error(err.str());
249double PseudoJet::pseudorapidity()
const {
250 if (px() == 0.0 && py() ==0.0)
return MaxRap;
251 if (pz() == 0.0)
return 0.0;
253 double theta = atan(perp()/pz());
255 return -log(tan(
theta/2));
265 jet1.E() +jet2.E() );
270PseudoJet operator- (
const PseudoJet & jet1,
const PseudoJet & jet2) {
272 return PseudoJet(jet1.px()-jet2.px(),
275 jet1.E() -jet2.E() );
280PseudoJet operator* (
double coeff,
const PseudoJet & jet) {
284 jet._ensure_valid_rap_phi();
287 PseudoJet coeff_times_jet(jet);
288 coeff_times_jet *= coeff;
289 return coeff_times_jet;
294PseudoJet operator* (
const PseudoJet & jet,
double coeff) {
300PseudoJet operator/ (
const PseudoJet & jet,
double coeff) {
301 return (1.0/coeff)*jet;
314 _ensure_valid_rap_phi();
327 (*this) *= 1.0/coeff;
335 _px += other_jet._px;
336 _py += other_jet._py;
337 _pz += other_jet._pz;
347 _px -= other_jet._px;
348 _py -= other_jet._py;
349 _pz -= other_jet._pz;
357 if (a.px() != b.px())
return false;
358 if (a.py() != b.py())
return false;
359 if (a.pz() != b.pz())
return false;
360 if (a.E () != b.E ())
return false;
374 throw Error(
"comparing a PseudoJet with a non-zero constant (double) is not allowed.");
375 return (jet.px() == 0 && jet.py() == 0 &&
376 jet.pz() == 0 && jet.E() == 0);
389 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
392 double m_local = prest.
m();
393 assert(m_local != 0);
395 double pf4 = ( px()*prest.px() + py()*prest.py()
396 + pz()*prest.pz() + E()*prest.E() )/m_local;
397 double fn = (pf4 + E()) / (prest.E() + m_local);
398 _px += fn*prest.px();
399 _py += fn*prest.py();
400 _pz += fn*prest.pz();
416 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
419 double m_local = prest.
m();
420 assert(m_local != 0);
422 double pf4 = ( -px()*prest.px() - py()*prest.py()
423 - pz()*prest.pz() + E()*prest.E() )/m_local;
424 double fn = (pf4 + E()) / (prest.E() + m_local);
425 _px -= fn*prest.px();
426 _py -= fn*prest.py();
427 _pz -= fn*prest.pz();
438 return jeta.px() == jetb.px()
439 && jeta.py() == jetb.py()
440 && jeta.pz() == jetb.pz()
441 && jeta.E() == jetb.E();
445void PseudoJet::set_cached_rap_phi(
double rap_in,
double phi_in) {
446 _rap = rap_in; _phi = phi_in;
447 if (_phi >= twopi) _phi -= twopi;
448 if (_phi < 0) _phi += twopi;
449#ifdef FASTJET_HAVE_THREAD_SAFETY
450 _init_status = Init_Done;
455void PseudoJet::reset_momentum_PtYPhiM(
double pt_in,
double y_in,
double phi_in,
double m_in) {
456 assert(phi_in < 2*twopi && phi_in > -twopi);
457 double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
458 double exprap = exp(y_in);
459 double pminus = ptm/exprap;
460 double pplus = ptm*exprap;
461 double px_local = pt_in*cos(phi_in);
462 double py_local = pt_in*sin(phi_in);
463 reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
464 set_cached_rap_phi(y_in,phi_in);
470 assert(phi < 2*twopi && phi > -twopi);
471 double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
472 double exprap = exp(y);
473 double pminus = ptm/exprap;
474 double pplus = ptm*exprap;
475 double px = pt*cos(phi);
476 double py = pt*sin(phi);
477 PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
486double PseudoJet::kt_distance(
const PseudoJet & other)
const {
488 double distance = min(_kt2, other._kt2);
489 double dphi = abs(phi() - other.
phi());
490 if (dphi > pi) {dphi = twopi - dphi;}
491 double drap = rap() - other.
rap();
492 distance = distance * (dphi*dphi + drap*drap);
499double PseudoJet::plain_distance(
const PseudoJet & other)
const {
500 double dphi = abs(phi() - other.
phi());
501 if (dphi > pi) {dphi = twopi - dphi;}
502 double drap = rap() - other.
rap();
503 return (dphi*dphi + drap*drap);
509double PseudoJet::delta_phi_to(
const PseudoJet & other)
const {
510 double dphi = other.
phi() - phi();
511 if (dphi > pi) dphi -= twopi;
512 if (dphi < -pi) dphi += twopi;
517string PseudoJet::description()
const{
520 return "standard PseudoJet (with no associated clustering information)";
523 return _structure->description();
538bool PseudoJet::has_associated_cluster_sequence()
const{
539 return (_structure) && (_structure->has_associated_cluster_sequence());
546 if (! has_associated_cluster_sequence())
return NULL;
548 return _structure->associated_cluster_sequence();
555bool PseudoJet::has_valid_cluster_sequence()
const{
556 return (_structure) && (_structure->has_valid_cluster_sequence());
566 return validated_structure_ptr()->validated_cs();
577 _structure = structure_in;
582bool PseudoJet::has_structure()
const{
583 return (
bool) _structure;
593 return _structure.get();
608 return _structure.get();
618 throw Error(
"Trying to access the structure of a PseudoJet which has no associated structure");
619 return _structure.get();
638 return validated_structure_ptr()->has_partner(*
this, partner);
649 return validated_structure_ptr()->has_child(*
this, child);
660 return validated_structure_ptr()->has_parents(*
this, parent1, parent2);
669bool PseudoJet::contains(
const PseudoJet &constituent)
const{
670 return validated_structure_ptr()->object_in_jet(constituent, *
this);
680 return validated_structure_ptr()->object_in_jet(*
this, jet);
686bool PseudoJet::has_constituents()
const{
687 return (_structure) && (_structure->has_constituents());
692vector<PseudoJet> PseudoJet::constituents()
const{
693 return validated_structure_ptr()->constituents(*
this);
699bool PseudoJet::has_exclusive_subjets()
const{
700 return (_structure) && (_structure->has_exclusive_subjets());
715std::vector<PseudoJet> PseudoJet::exclusive_subjets (
const double dcut)
const {
716 return validated_structure_ptr()->exclusive_subjets(*
this, dcut);
726int PseudoJet::n_exclusive_subjets(
const double dcut)
const {
727 return validated_structure_ptr()->n_exclusive_subjets(*
this, dcut);
739std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (
int nsub)
const {
740 return validated_structure_ptr()->exclusive_subjets_up_to(*
this, nsub);
746std::vector<PseudoJet> PseudoJet::exclusive_subjets (
int nsub)
const {
747 vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
748 if (
int(subjets.size()) < nsub) {
750 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
751 << subjets.size() <<
" particles in the jet";
752 throw Error(err.str());
763double PseudoJet::exclusive_subdmerge(
int nsub)
const {
764 return validated_structure_ptr()->exclusive_subdmerge(*
this, nsub);
774double PseudoJet::exclusive_subdmerge_max(
int nsub)
const {
775 return validated_structure_ptr()->exclusive_subdmerge_max(*
this, nsub);
783bool PseudoJet::has_pieces()
const{
784 return ((_structure) && (_structure->has_pieces(*
this)));
792std::vector<PseudoJet> PseudoJet::pieces()
const{
793 return validated_structure_ptr()->pieces(*
this);
813 if (csab == NULL)
throw Error(
"you requested jet-area related information, but the PseudoJet does not have associated area information.");
820bool PseudoJet::has_area()
const{
822 if (! has_structure())
return false;
823 return (validated_structure_ptr()->has_area() != 0);
829double PseudoJet::area()
const{
830 return validated_structure_ptr()->
area(*
this);
837double PseudoJet::area_error()
const{
838 return validated_structure_ptr()->area_error(*
this);
851bool PseudoJet::is_pure_ghost()
const{
866PseudoJet::InexistentUserInfo::InexistentUserInfo() :
Error(
"you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
873void sort_indices(vector<int> & indices,
874 const vector<double> & values) {
876 sort(indices.begin(), indices.end(), index_sort_helper);
883 vector<double> minus_kt2(jets.size());
884 for (
size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
891 vector<double> rapidities(jets.size());
892 for (
size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
899 vector<double> energies(jets.size());
900 for (
size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
907 vector<double> pz(jets.size());
908 for (
size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
926 for (
unsigned int i=0; i<
pieces.size(); i++)
931 CompositeJetStructure *cj_struct =
new CompositeJetStructure(
pieces);
933 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
940 return join(vector<PseudoJet>(1,j1));
base class that sets interface for extensions of ClusterSequence that provide information about the a...
virtual double area(const PseudoJet &) const
return the area associated with the given jet; this base class returns 0.
base class corresponding to errors that can be thrown by FastJet
Contains any information related to the clustering that should be directly accessible to PseudoJet.
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
PseudoJet()
default constructor, which as of FJ3.0 provides an object for which all operations are now valid and ...
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...
double rap() const
returns the rapidity or some large value when the rapidity is infinite
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
double phi() const
returns phi (in the range 0..2pi)
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet.
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
int user_index() const
return the user_index,
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons,...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
bool have_same_momentum(const PseudoJet &jeta, const PseudoJet &jetb)
returns true if the momenta of the two input jets are identical
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...
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,...
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
double theta(const PseudoJet &a, const PseudoJet &b)
returns the angle between a and b
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
vector< PseudoJet > sorted_by_pz(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing pz