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"
46 FASTJET_BEGIN_NAMESPACE
53 PseudoJet::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);
92 #endif // FASTJET_HAVE_THREAD_SAFETY
121 void 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
140 void 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));
185 void 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;}
217 valarray<double> PseudoJet::four_mom()
const {
218 valarray<double> mom(4);
229 double PseudoJet::operator () (
int i)
const {
241 err <<
"PseudoJet subscripting: bad index (" << i <<
")";
242 throw Error(err.str());
249 double 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() );
270 PseudoJet operator- (
const PseudoJet & jet1,
const PseudoJet & jet2) {
272 return PseudoJet(jet1.px()-jet2.px(),
275 jet1.E() -jet2.E() );
280 PseudoJet
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;
294 PseudoJet
operator* (
const PseudoJet & jet,
double coeff) {
300 PseudoJet 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();
445 void 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;
455 void 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));
486 double 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);
499 double 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);
509 double 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;
517 string PseudoJet::description()
const{
520 return "standard PseudoJet (with no associated clustering information)";
523 return _structure->description();
538 bool 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();
555 bool 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;
582 bool 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);
669 bool 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);
686 bool PseudoJet::has_constituents()
const{
687 return (_structure) && (_structure->has_constituents());
692 vector<PseudoJet> PseudoJet::constituents()
const{
693 return validated_structure_ptr()->constituents(*
this);
699 bool PseudoJet::has_exclusive_subjets()
const{
700 return (_structure) && (_structure->has_exclusive_subjets());
715 std::vector<PseudoJet> PseudoJet::exclusive_subjets (
const double dcut)
const {
716 return validated_structure_ptr()->exclusive_subjets(*
this, dcut);
726 int PseudoJet::n_exclusive_subjets(
const double dcut)
const {
727 return validated_structure_ptr()->n_exclusive_subjets(*
this, dcut);
739 std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (
int nsub)
const {
740 return validated_structure_ptr()->exclusive_subjets_up_to(*
this, nsub);
746 std::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());
763 double PseudoJet::exclusive_subdmerge(
int nsub)
const {
764 return validated_structure_ptr()->exclusive_subdmerge(*
this, nsub);
774 double PseudoJet::exclusive_subdmerge_max(
int nsub)
const {
775 return validated_structure_ptr()->exclusive_subdmerge_max(*
this, nsub);
783 bool PseudoJet::has_pieces()
const{
784 return ((_structure) && (_structure->has_pieces(*
this)));
792 std::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.");
820 bool PseudoJet::has_area()
const{
822 if (! has_structure())
return false;
823 return (validated_structure_ptr()->has_area() != 0);
829 double PseudoJet::area()
const{
830 return validated_structure_ptr()->
area(*
this);
837 double PseudoJet::area_error()
const{
838 return validated_structure_ptr()->area_error(*
this);
851 bool PseudoJet::is_pure_ghost()
const{
866 PseudoJet::InexistentUserInfo::InexistentUserInfo() :
Error(
"you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null")
873 void sort_indices(vector<int> & indices,
874 const vector<double> & values) {
875 IndexedSortHelper index_sort_helper(&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));
976 FASTJET_END_NAMESPACE