30 #include "fastjet/Error.hh" 31 #include "fastjet/PseudoJet.hh" 32 #include "fastjet/ClusterSequence.hh" 34 #include "fastjet/ClusterSequenceAreaBase.hh" 36 #include "fastjet/CompositeJetStructure.hh" 44 FASTJET_BEGIN_NAMESPACE
51 PseudoJet::PseudoJet(
const double px_in,
const double py_in,
const double pz_in,
const double E_in) {
68 void PseudoJet::_finish_init () {
69 _kt2 = this->px()*this->px() + this->py()*this->py();
78 _rap = pseudojet_invalid_rap;
82 void PseudoJet::_set_rap_phi()
const {
87 _phi = atan2(this->py(),this->px());
89 if (_phi < 0.0) {_phi += twopi;}
90 if (_phi >= twopi) {_phi -= twopi;}
91 if (this->E() == abs(this->pz()) && _kt2 == 0) {
96 double MaxRapHere =
MaxRap + abs(this->pz());
97 if (this->pz() >= 0.0) {_rap = MaxRapHere;}
else {_rap = -MaxRapHere;}
102 double effective_m2 = max(0.0,m2());
103 double E_plus_pz = _E + abs(_pz);
105 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
106 if (_pz > 0) {_rap = - _rap;}
114 valarray<double> PseudoJet::four_mom()
const {
115 valarray<double> mom(4);
126 double PseudoJet::operator () (
int i)
const {
138 err <<
"PseudoJet subscripting: bad index (" << i <<
")";
139 throw Error(err.str());
146 double PseudoJet::pseudorapidity()
const {
147 if (px() == 0.0 && py() ==0.0)
return MaxRap;
148 if (pz() == 0.0)
return 0.0;
150 double theta = atan(perp()/pz());
151 if (theta < 0) theta += pi;
152 return -log(tan(theta/2));
162 jet1.E() +jet2.E() );
172 jet1.E() -jet2.E() );
181 coeff_times_jet *= coeff;
182 return coeff_times_jet;
194 return (1.0/coeff)*jet;
199 void PseudoJet::operator*=(
double coeff) {
210 void PseudoJet::operator/=(
double coeff) {
211 (*this) *= 1.0/coeff;
217 void PseudoJet::operator+=(
const PseudoJet & other_jet) {
218 _px += other_jet._px;
219 _py += other_jet._py;
220 _pz += other_jet._pz;
228 void PseudoJet::operator-=(
const PseudoJet & other_jet) {
229 _px -= other_jet._px;
230 _py -= other_jet._py;
231 _pz -= other_jet._pz;
238 if (a.px() != b.px())
return false;
239 if (a.py() != b.py())
return false;
240 if (a.pz() != b.pz())
return false;
241 if (a.E () != b.E ())
return false;
255 throw Error(
"comparing a PseudoJet with a non-zero constant (double) is not allowed.");
256 return (jet.px() == 0 && jet.py() == 0 &&
257 jet.pz() == 0 && jet.E() == 0);
270 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
273 double m_local = prest.
m();
274 assert(m_local != 0);
276 double pf4 = ( px()*prest.px() + py()*prest.py()
277 + pz()*prest.pz() + E()*prest.E() )/m_local;
278 double fn = (pf4 + E()) / (prest.E() + m_local);
279 _px += fn*prest.px();
280 _py += fn*prest.py();
281 _pz += fn*prest.pz();
297 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
300 double m_local = prest.
m();
301 assert(m_local != 0);
303 double pf4 = ( -px()*prest.px() - py()*prest.py()
304 - pz()*prest.pz() + E()*prest.E() )/m_local;
305 double fn = (pf4 + E()) / (prest.E() + m_local);
306 _px -= fn*prest.px();
307 _py -= fn*prest.py();
308 _pz -= fn*prest.pz();
319 return jeta.px() == jetb.px()
320 && jeta.py() == jetb.py()
321 && jeta.pz() == jetb.pz()
322 && jeta.E() == jetb.E();
326 void PseudoJet::set_cached_rap_phi(
double rap_in,
double phi_in) {
327 _rap = rap_in; _phi = phi_in;
328 if (_phi >= twopi) _phi -= twopi;
329 if (_phi < 0) _phi += twopi;
333 void PseudoJet::reset_momentum_PtYPhiM(
double pt_in,
double y_in,
double phi_in,
double m_in) {
334 assert(phi_in < 2*twopi && phi_in > -twopi);
335 double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
336 double exprap = exp(y_in);
337 double pminus = ptm/exprap;
338 double pplus = ptm*exprap;
339 double px_local = pt_in*cos(phi_in);
340 double py_local = pt_in*sin(phi_in);
341 reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
342 set_cached_rap_phi(y_in,phi_in);
348 assert(phi < 2*twopi && phi > -twopi);
349 double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
350 double exprap = exp(y);
351 double pminus = ptm/exprap;
352 double pplus = ptm*exprap;
353 double px = pt*cos(phi);
354 double py = pt*sin(phi);
355 PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
364 double PseudoJet::kt_distance(
const PseudoJet & other)
const {
366 double distance = min(_kt2, other._kt2);
367 double dphi = abs(phi() - other.
phi());
368 if (dphi > pi) {dphi = twopi - dphi;}
369 double drap = rap() - other.
rap();
370 distance = distance * (dphi*dphi + drap*drap);
377 double PseudoJet::plain_distance(
const PseudoJet & other)
const {
378 double dphi = abs(phi() - other.
phi());
379 if (dphi > pi) {dphi = twopi - dphi;}
380 double drap = rap() - other.
rap();
381 return (dphi*dphi + drap*drap);
387 double PseudoJet::delta_phi_to(
const PseudoJet & other)
const {
388 double dphi = other.
phi() - phi();
389 if (dphi > pi) dphi -= twopi;
390 if (dphi < -pi) dphi += twopi;
395 string PseudoJet::description()
const{
398 return "standard PseudoJet (with no associated clustering information)";
401 return _structure()->description();
416 bool PseudoJet::has_associated_cluster_sequence()
const{
417 return (_structure()) && (_structure->has_associated_cluster_sequence());
424 if (! has_associated_cluster_sequence())
return NULL;
426 return _structure->associated_cluster_sequence();
433 bool PseudoJet::has_valid_cluster_sequence()
const{
434 return (_structure()) && (_structure->has_valid_cluster_sequence());
444 return validated_structure_ptr()->validated_cs();
451 _structure = structure_in;
456 bool PseudoJet::has_structure()
const{
466 if (!_structure())
return NULL;
481 if (!_structure())
return NULL;
492 throw Error(
"Trying to access the structure of a PseudoJet which has no associated structure");
512 return validated_structure_ptr()->has_partner(*
this, partner);
523 return validated_structure_ptr()->has_child(*
this, child);
534 return validated_structure_ptr()->has_parents(*
this, parent1, parent2);
543 bool PseudoJet::contains(
const PseudoJet &constituent)
const{
544 return validated_structure_ptr()->object_in_jet(constituent, *
this);
554 return validated_structure_ptr()->object_in_jet(*
this, jet);
560 bool PseudoJet::has_constituents()
const{
561 return (_structure()) && (_structure->has_constituents());
566 vector<PseudoJet> PseudoJet::constituents()
const{
567 return validated_structure_ptr()->constituents(*
this);
573 bool PseudoJet::has_exclusive_subjets()
const{
574 return (_structure()) && (_structure->has_exclusive_subjets());
589 std::vector<PseudoJet> PseudoJet::exclusive_subjets (
const double & dcut)
const {
590 return validated_structure_ptr()->exclusive_subjets(*
this, dcut);
600 int PseudoJet::n_exclusive_subjets(
const double & dcut)
const {
601 return validated_structure_ptr()->n_exclusive_subjets(*
this, dcut);
613 std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (
int nsub)
const {
614 return validated_structure_ptr()->exclusive_subjets_up_to(*
this, nsub);
620 std::vector<PseudoJet> PseudoJet::exclusive_subjets (
int nsub)
const {
621 vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
622 if (
int(subjets.size()) < nsub) {
624 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only " 625 << subjets.size() <<
" particles in the jet";
626 throw Error(err.str());
637 double PseudoJet::exclusive_subdmerge(
int nsub)
const {
638 return validated_structure_ptr()->exclusive_subdmerge(*
this, nsub);
648 double PseudoJet::exclusive_subdmerge_max(
int nsub)
const {
649 return validated_structure_ptr()->exclusive_subdmerge_max(*
this, nsub);
657 bool PseudoJet::has_pieces()
const{
658 return ((_structure()) && (_structure->has_pieces(*
this)));
666 std::vector<PseudoJet> PseudoJet::pieces()
const{
667 return validated_structure_ptr()->pieces(*
this);
687 if (csab == NULL)
throw Error(
"you requested jet-area related information, but the PseudoJet does not have associated area information.");
694 bool PseudoJet::has_area()
const{
696 if (! has_structure())
return false;
697 return (validated_structure_ptr()->has_area() != 0);
703 double PseudoJet::area()
const{
704 return validated_structure_ptr()->
area(*
this);
711 double PseudoJet::area_error()
const{
712 return validated_structure_ptr()->area_error(*
this);
725 bool PseudoJet::is_pure_ghost()
const{
740 PseudoJet::InexistentUserInfo::InexistentUserInfo() :
Error(
"you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
747 void sort_indices(vector<int> & indices,
748 const vector<double> & values) {
749 IndexedSortHelper index_sort_helper(&values);
750 sort(indices.begin(), indices.end(), index_sort_helper);
760 const vector<T> & objects,
761 const vector<double> & values) {
763 assert(objects.size() == values.size());
766 vector<int> indices(values.size());
767 for (
size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
770 sort_indices(indices, values);
773 vector<T> objects_sorted(objects.size());
776 for (
size_t i = 0; i < indices.size(); i++) {
777 objects_sorted[i] = objects[indices[i]];
780 return objects_sorted;
786 vector<double> minus_kt2(jets.size());
787 for (
size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
794 vector<double> rapidities(jets.size());
795 for (
size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
802 vector<double> energies(jets.size());
803 for (
size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
810 vector<double> pz(jets.size());
811 for (
size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
829 for (
unsigned int i=0; i<pieces.size(); i++)
843 return join(vector<PseudoJet>(1,j1));
849 pieces.push_back(j1);
850 pieces.push_back(j2);
857 pieces.push_back(j1);
858 pieces.push_back(j2);
859 pieces.push_back(j3);
866 pieces.push_back(j1);
867 pieces.push_back(j2);
868 pieces.push_back(j3);
869 pieces.push_back(j4);
876 FASTJET_END_NAMESPACE
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
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...
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons, giving rapidity=infinity.
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
Selector operator*(const Selector &s1, const Selector &s2)
successive application of 2 selectors
vector< PseudoJet > sorted_by_pz(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing pz
vector< PseudoJet > sorted_by_rapidity(const vector< PseudoJet > &jets)
return a vector of jets sorted into increasing rapidity
The structure for a jet made of pieces.
virtual double area(const PseudoJet &) const
return the area associated with the given jet; this base class returns 0.
Contains any information related to the clustering that should be directly accessible to PseudoJet...
vector< T > objects_sorted_by_values(const vector< T > &objects, const vector< double > &values)
given a vector of values with a one-to-one correspondence with the vector of objects, sort objects into an order such that the associated values would be in increasing order
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...
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet...
base class that sets interface for extensions of ClusterSequence that provide information about the a...
vector< PseudoJet > sorted_by_E(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing energy
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
base class corresponding to errors that can be thrown by FastJet
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
an implementation of C++0x shared pointers (or boost's)
double phi() const
returns phi (in the range 0..2pi)
double rap() const
returns the rapidity or some large value when the rapidity is infinite
bool have_same_momentum(const PseudoJet &jeta, const PseudoJet &jetb)
returns true if the momenta of the two input jets are identical
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
const double pseudojet_invalid_phi
default value for phi, meaning it (and rapidity) have yet to be calculated)
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
int user_index() const
return the user_index,