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() );
167 PseudoJet operator- (
const PseudoJet & jet1,
const PseudoJet & jet2) {
169 return PseudoJet(jet1.px()-jet2.px(),
172 jet1.E() -jet2.E() );
177 PseudoJet
operator* (
double coeff,
const PseudoJet & jet) {
180 PseudoJet coeff_times_jet(jet);
181 coeff_times_jet *= coeff;
182 return coeff_times_jet;
187 PseudoJet
operator* (
const PseudoJet & jet,
double coeff) {
193 PseudoJet operator/ (
const PseudoJet & jet,
double coeff) {
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++)
834 CompositeJetStructure *cj_struct =
new CompositeJetStructure(pieces);
836 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
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
double rap() const
returns the rapidity or some large value when the rapidity is infinite
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
const UserInfoBase * user_info_ptr() const
retrieve a pointer to the (const) user information
const double MaxRap
Used to protect against parton-level events where pt can be zero for some partons, giving rapidity=infinity.
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
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
Contains any information related to the clustering that should be directly accessible to PseudoJet...
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
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 double area(const PseudoJet &) const
return the area associated with the given jet; this base class returns 0.
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
base class corresponding to errors that can be thrown by FastJet
an implementation of C++0x shared pointers (or boost's)
virtual bool is_pure_ghost() const
true if this jet is made exclusively of ghosts.
bool have_same_momentum(const PseudoJet &jeta, const PseudoJet &jetb)
returns true if the momenta of the two input jets are identical
double phi() const
returns phi (in the range 0..2pi)
int user_index() const
return the user_index,
PseudoJet()
default constructor, which as of FJ3.0 provides an object for which all operations are now valid and ...
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
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...
const PseudoJetStructureBase * structure_ptr() const
return a pointer to the structure (of type PseudoJetStructureBase*) associated with this PseudoJet...
virtual PseudoJet area_4vector() const
return the jet 4-vector area.