30 #include "fastjet/Error.hh"
31 #include "fastjet/PseudoJet.hh"
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/ClusterSequenceAreaBase.hh"
34 #include "fastjet/CompositeJetStructure.hh"
42 FASTJET_BEGIN_NAMESPACE
49 PseudoJet::PseudoJet(
const double px_in,
const double py_in,
const double pz_in,
const double E_in) {
66 void PseudoJet::_finish_init () {
67 _kt2 = this->px()*this->px() + this->py()*this->py();
72 void PseudoJet::_set_rap_phi()
const {
77 _phi = atan2(this->py(),this->px());
79 if (_phi < 0.0) {_phi += twopi;}
80 if (_phi >= twopi) {_phi -= twopi;}
81 if (this->E() == abs(this->pz()) && _kt2 == 0) {
86 double MaxRapHere =
MaxRap + abs(this->pz());
87 if (this->pz() >= 0.0) {_rap = MaxRapHere;}
else {_rap = -MaxRapHere;}
92 double effective_m2 = max(0.0,m2());
93 double E_plus_pz = _E + abs(_pz);
95 _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz));
96 if (_pz > 0) {_rap = - _rap;}
104 valarray<double> PseudoJet::four_mom()
const {
105 valarray<double> mom(4);
116 double PseudoJet::operator () (
int i)
const {
128 err <<
"PseudoJet subscripting: bad index (" << i <<
")";
129 throw Error(err.str());
136 double PseudoJet::pseudorapidity()
const {
137 if (px() == 0.0 && py() ==0.0)
return MaxRap;
138 if (pz() == 0.0)
return 0.0;
140 double theta = atan(perp()/pz());
141 if (theta < 0) theta += pi;
142 return -log(tan(theta/2));
152 jet1.E() +jet2.E() );
157 PseudoJet operator- (
const PseudoJet & jet1,
const PseudoJet & jet2) {
159 return PseudoJet(jet1.px()-jet2.px(),
162 jet1.E() -jet2.E() );
167 PseudoJet
operator* (
double coeff,
const PseudoJet & jet) {
170 PseudoJet coeff_times_jet(jet);
171 coeff_times_jet *= coeff;
172 return coeff_times_jet;
177 PseudoJet
operator* (
const PseudoJet & jet,
double coeff) {
183 PseudoJet operator/ (
const PseudoJet & jet,
double coeff) {
184 return (1.0/coeff)*jet;
189 void PseudoJet::operator*=(
double coeff) {
200 void PseudoJet::operator/=(
double coeff) {
201 (*this) *= 1.0/coeff;
207 void PseudoJet::operator+=(
const PseudoJet & other_jet) {
208 _px += other_jet._px;
209 _py += other_jet._py;
210 _pz += other_jet._pz;
218 void PseudoJet::operator-=(
const PseudoJet & other_jet) {
219 _px -= other_jet._px;
220 _py -= other_jet._py;
221 _pz -= other_jet._pz;
228 if (a.px() != b.px())
return false;
229 if (a.py() != b.py())
return false;
230 if (a.pz() != b.pz())
return false;
231 if (a.E () != b.E ())
return false;
245 throw Error(
"comparing a PseudoJet with a non-zero constant (double) is not allowed.");
246 return (jet.px() == 0 && jet.py() == 0 &&
247 jet.pz() == 0 && jet.E() == 0);
260 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
263 double m_local = prest.
m();
264 assert(m_local != 0);
266 double pf4 = ( px()*prest.px() + py()*prest.py()
267 + pz()*prest.pz() + E()*prest.E() )/m_local;
268 double fn = (pf4 + E()) / (prest.E() + m_local);
269 _px += fn*prest.px();
270 _py += fn*prest.py();
271 _pz += fn*prest.pz();
287 if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0)
290 double m_local = prest.
m();
291 assert(m_local != 0);
293 double pf4 = ( -px()*prest.px() - py()*prest.py()
294 - pz()*prest.pz() + E()*prest.E() )/m_local;
295 double fn = (pf4 + E()) / (prest.E() + m_local);
296 _px -= fn*prest.px();
297 _py -= fn*prest.py();
298 _pz -= fn*prest.pz();
309 return jeta.px() == jetb.px()
310 && jeta.py() == jetb.py()
311 && jeta.pz() == jetb.pz()
312 && jeta.E() == jetb.E();
316 void PseudoJet::set_cached_rap_phi(
double rap_in,
double phi_in) {
317 _rap = rap_in; _phi = phi_in;
318 if (_phi >= twopi) _phi -= twopi;
319 if (_phi < 0) _phi += twopi;
323 void PseudoJet::reset_momentum_PtYPhiM(
double pt_in,
double y_in,
double phi_in,
double m_in) {
324 assert(phi_in < 2*twopi && phi_in > -twopi);
325 double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in);
326 double exprap = exp(y_in);
327 double pminus = ptm/exprap;
328 double pplus = ptm*exprap;
329 double px_local = pt_in*cos(phi_in);
330 double py_local = pt_in*sin(phi_in);
331 reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus));
332 set_cached_rap_phi(y_in,phi_in);
338 assert(phi < 2*twopi && phi > -twopi);
339 double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m);
340 double exprap = exp(y);
341 double pminus = ptm/exprap;
342 double pplus = ptm*exprap;
343 double px = pt*cos(phi);
344 double py = pt*sin(phi);
345 PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus));
354 double PseudoJet::kt_distance(
const PseudoJet & other)
const {
356 double distance = min(_kt2, other._kt2);
357 double dphi = abs(phi() - other.
phi());
358 if (dphi > pi) {dphi = twopi - dphi;}
359 double drap = rap() - other.
rap();
360 distance = distance * (dphi*dphi + drap*drap);
367 double PseudoJet::plain_distance(
const PseudoJet & other)
const {
368 double dphi = abs(phi() - other.
phi());
369 if (dphi > pi) {dphi = twopi - dphi;}
370 double drap = rap() - other.
rap();
371 return (dphi*dphi + drap*drap);
377 double PseudoJet::delta_phi_to(
const PseudoJet & other)
const {
378 double dphi = other.
phi() - phi();
379 if (dphi > pi) dphi -= twopi;
380 if (dphi < -pi) dphi += twopi;
385 string PseudoJet::description()
const{
388 return "standard PseudoJet (with no associated clustering information)";
391 return _structure()->description();
406 bool PseudoJet::has_associated_cluster_sequence()
const{
407 return (_structure()) && (_structure->has_associated_cluster_sequence());
414 if (! has_associated_cluster_sequence())
return NULL;
416 return _structure->associated_cluster_sequence();
423 bool PseudoJet::has_valid_cluster_sequence()
const{
424 return (_structure()) && (_structure->has_valid_cluster_sequence());
434 return validated_structure_ptr()->validated_cs();
441 _structure = structure_in;
446 bool PseudoJet::has_structure()
const{
456 if (!_structure())
return NULL;
471 if (!_structure())
return NULL;
482 throw Error(
"Trying to access the structure of a PseudoJet which has no associated structure");
502 return validated_structure_ptr()->
has_partner(*
this, partner);
513 return validated_structure_ptr()->has_child(*
this, child);
524 return validated_structure_ptr()->has_parents(*
this, parent1, parent2);
533 bool PseudoJet::contains(
const PseudoJet &constituent)
const{
534 return validated_structure_ptr()->object_in_jet(constituent, *
this);
544 return validated_structure_ptr()->object_in_jet(*
this, jet);
550 bool PseudoJet::has_constituents()
const{
551 return (_structure()) && (_structure->has_constituents());
556 vector<PseudoJet> PseudoJet::constituents()
const{
557 return validated_structure_ptr()->constituents(*
this);
563 bool PseudoJet::has_exclusive_subjets()
const{
564 return (_structure()) && (_structure->has_exclusive_subjets());
579 std::vector<PseudoJet> PseudoJet::exclusive_subjets (
const double & dcut)
const {
580 return validated_structure_ptr()->exclusive_subjets(*
this, dcut);
590 int PseudoJet::n_exclusive_subjets(
const double & dcut)
const {
591 return validated_structure_ptr()->n_exclusive_subjets(*
this, dcut);
603 std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (
int nsub)
const {
604 return validated_structure_ptr()->exclusive_subjets_up_to(*
this, nsub);
610 std::vector<PseudoJet> PseudoJet::exclusive_subjets (
int nsub)
const {
611 vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub);
612 if (
int(subjets.size()) < nsub) {
614 err <<
"Requested " << nsub <<
" exclusive subjets, but there were only "
615 << subjets.size() <<
" particles in the jet";
616 throw Error(err.str());
627 double PseudoJet::exclusive_subdmerge(
int nsub)
const {
628 return validated_structure_ptr()->exclusive_subdmerge(*
this, nsub);
638 double PseudoJet::exclusive_subdmerge_max(
int nsub)
const {
639 return validated_structure_ptr()->exclusive_subdmerge_max(*
this, nsub);
647 bool PseudoJet::has_pieces()
const{
648 return ((_structure()) && (_structure->has_pieces(*
this)));
656 std::vector<PseudoJet> PseudoJet::pieces()
const{
657 return validated_structure_ptr()->pieces(*
this);
675 if (csab == NULL)
throw Error(
"you requested jet-area related information, but the PseudoJet does not have associated area information.");
682 bool PseudoJet::has_area()
const{
684 if (! has_structure())
return false;
685 return (validated_structure_ptr()->has_area() != 0);
691 double PseudoJet::area()
const{
692 return validated_structure_ptr()->
area(*
this);
699 double PseudoJet::area_error()
const{
700 return validated_structure_ptr()->area_error(*
this);
713 bool PseudoJet::is_pure_ghost()
const{
727 PseudoJet::InexistentUserInfo::InexistentUserInfo() :
Error(
"you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
734 void sort_indices(vector<int> & indices,
735 const vector<double> & values) {
736 IndexedSortHelper index_sort_helper(&values);
737 sort(indices.begin(), indices.end(), index_sort_helper);
747 const vector<T> & objects,
748 const vector<double> & values) {
750 assert(objects.size() == values.size());
753 vector<int> indices(values.size());
754 for (
size_t i = 0; i < indices.size(); i++) {indices[i] = i;}
757 sort_indices(indices, values);
760 vector<T> objects_sorted(objects.size());
763 for (
size_t i = 0; i < indices.size(); i++) {
764 objects_sorted[i] = objects[indices[i]];
767 return objects_sorted;
773 vector<double> minus_kt2(jets.size());
774 for (
size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();}
781 vector<double> rapidities(jets.size());
782 for (
size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();}
789 vector<double> energies(jets.size());
790 for (
size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();}
797 vector<double> pz(jets.size());
798 for (
size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();}
816 for (
unsigned int i=0; i<pieces.size(); i++)
821 CompositeJetStructure *cj_struct =
new CompositeJetStructure(pieces);
823 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
830 return join(vector<PseudoJet>(1,j1));
836 pieces.push_back(j1);
837 pieces.push_back(j2);
844 pieces.push_back(j1);
845 pieces.push_back(j2);
846 pieces.push_back(j3);
853 pieces.push_back(j1);
854 pieces.push_back(j2);
855 pieces.push_back(j3);
856 pieces.push_back(j4);
863 FASTJET_END_NAMESPACE