]>
Commit | Line | Data |
---|---|---|
370be031 | 1 | //STARTHEADER |
2 | // $Id: PseudoJet.hh 1510 2009-04-13 08:48:41Z salam $ | |
3 | // | |
4 | // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam | |
5 | // | |
6 | //---------------------------------------------------------------------- | |
7 | // This file is part of FastJet. | |
8 | // | |
9 | // FastJet is free software; you can redistribute it and/or modify | |
10 | // it under the terms of the GNU General Public License as published by | |
11 | // the Free Software Foundation; either version 2 of the License, or | |
12 | // (at your option) any later version. | |
13 | // | |
14 | // The algorithms that underlie FastJet have required considerable | |
15 | // development and are described in hep-ph/0512210. If you use | |
16 | // FastJet as part of work towards a scientific publication, please | |
17 | // include a citation to the FastJet paper. | |
18 | // | |
19 | // FastJet is distributed in the hope that it will be useful, | |
20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
22 | // GNU General Public License for more details. | |
23 | // | |
24 | // You should have received a copy of the GNU General Public License | |
25 | // along with FastJet; if not, write to the Free Software | |
26 | // Foundation, Inc.: | |
27 | // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA | |
28 | //---------------------------------------------------------------------- | |
29 | //ENDHEADER | |
30 | ||
31 | ||
32 | #ifndef __FASTJET_PSEUDOJET_HH__ | |
33 | #define __FASTJET_PSEUDOJET_HH__ | |
34 | ||
35 | #include<valarray> | |
36 | #include<vector> | |
37 | #include<cassert> | |
38 | #include<cmath> | |
39 | #include<iostream> | |
40 | #include "fastjet/internal/numconsts.hh" | |
41 | ||
42 | namespace fastjet { // defined in fastjet/internal/base.hh | |
43 | ||
44 | //using namespace std; | |
45 | ||
46 | /// Used to protect against parton-level events where pt can be zero | |
47 | /// for some partons, giving rapidity=infinity. KtJet fails in those cases. | |
48 | const double MaxRap = 1e5; | |
49 | ||
50 | /// Class to contain pseudojets, including minimal information of use to | |
51 | /// to jet-clustering routines. | |
52 | class PseudoJet { | |
53 | ||
54 | public: | |
55 | PseudoJet() {}; | |
56 | /// construct a pseudojet from explicit components | |
57 | PseudoJet(const double px, const double py, const double pz, const double E); | |
58 | /// constructor from any object that has px,py,pz,E = some_four_vector[0--3], | |
59 | template <class L> PseudoJet(const L & some_four_vector) ; | |
60 | ||
61 | // first "const double &" says that result is a reference to the | |
62 | // stored value and that we will not change that stored value. | |
63 | // | |
64 | // second "const" says that "this" will not be modified by these | |
65 | // functions. | |
66 | inline double E() const {return _E;} | |
67 | inline double e() const {return _E;} // like CLHEP | |
68 | inline double px() const {return _px;} | |
69 | inline double py() const {return _py;} | |
70 | inline double pz() const {return _pz;} | |
71 | ||
72 | /// returns phi (in the range 0..2pi) | |
73 | inline double phi() const {return phi_02pi();} | |
74 | ||
75 | /// returns phi in the range -pi..pi | |
76 | inline double phi_std() const { | |
77 | return _phi > pi ? _phi-twopi : _phi;} | |
78 | ||
79 | /// returns phi in the range 0..2pi | |
80 | inline double phi_02pi() const {return _phi;} | |
81 | ||
82 | /// returns the rapidity or some large value when the rapidity | |
83 | /// is infinite | |
84 | inline double rap() const {return _rap;} | |
85 | ||
86 | /// the same as rap() | |
87 | inline double rapidity() const {return _rap;} // like CLHEP | |
88 | ||
89 | /// returns the pseudo-rapidity or some large value when the | |
90 | /// rapidity is infinite | |
91 | double pseudorapidity() const; | |
92 | double eta() const {return pseudorapidity();} | |
93 | ||
94 | /// returns the squared transverse momentum | |
95 | inline double kt2() const {return _kt2;} | |
96 | /// returns the squared transverse momentum | |
97 | inline double perp2() const {return _kt2;} // like CLHEP | |
98 | /// returns the scalar transverse momentum | |
99 | inline double perp() const {return sqrt(_kt2);} // like CLHEP | |
100 | /// returns the squared invariant mass // like CLHEP | |
101 | inline double m2() const {return (_E+_pz)*(_E-_pz)-_kt2;} | |
102 | /// returns the squared transverse mass = kt^2+m^2 | |
103 | inline double mperp2() const {return (_E+_pz)*(_E-_pz);} | |
104 | /// returns the transverse mass = sqrt(kt^2+m^2) | |
105 | inline double mperp() const {return sqrt(std::abs(mperp2()));} | |
106 | /// returns the invariant mass | |
107 | /// (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) | |
108 | inline double m() const; | |
109 | /// return px^2+py^2+pz^2 | |
110 | inline double modp2() const {return _kt2+_pz*_pz;} | |
111 | /// return the transverse energy | |
112 | inline double Et() const {return (_kt2==0) ? 0.0 : _E/sqrt(1.0+_pz*_pz/_kt2);} | |
113 | /// return the transverse energy squared | |
114 | inline double Et2() const {return (_kt2==0) ? 0.0 : _E*_E/(1.0+_pz*_pz/_kt2);} | |
115 | ||
116 | /// returns component i, where X==0, Y==1, Z==2, E==3 | |
117 | double operator () (int i) const ; | |
118 | /// returns component i, where X==0, Y==1, Z==2, E==3 | |
119 | inline double operator [] (int i) const { return (*this)(i); }; // this too | |
120 | ||
121 | ||
122 | // taken from CLHEP | |
123 | enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES }; | |
124 | ||
125 | ||
126 | /// transform this jet (given in the rest frame of prest) into a jet | |
127 | /// in the lab frame [NOT FULLY TESTED] | |
128 | PseudoJet & boost(const PseudoJet & prest); | |
129 | /// transform this jet (given in lab) into a jet in the rest | |
130 | /// frame of prest [NOT FULLY TESTED] | |
131 | PseudoJet & unboost(const PseudoJet & prest); | |
132 | ||
133 | /// return the cluster_hist_index, intended to be used by clustering | |
134 | /// routines. | |
135 | inline int cluster_hist_index() const {return _cluster_hist_index;} | |
136 | /// set the cluster_hist_index, intended to be used by clustering routines. | |
137 | inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;} | |
138 | ||
139 | /// alternative name for cluster_hist_index() [perhaps more meaningful] | |
140 | inline int cluster_sequence_history_index() const { | |
141 | return cluster_hist_index();} | |
142 | /// alternative name for set_cluster_hist_index(...) [perhaps more | |
143 | /// meaningful] | |
144 | inline void set_cluster_sequence_history_index(const int index) { | |
145 | set_cluster_hist_index(index);} | |
146 | ||
147 | ||
148 | /// return the user_index, intended to allow the user to "add" information | |
149 | inline int user_index() const {return _user_index;} | |
150 | /// set the user_index, intended to allow the user to "add" information | |
151 | inline void set_user_index(const int index) {_user_index = index;} | |
152 | ||
153 | /// return a valarray containing the four-momentum (components 0-2 | |
154 | /// are 3-mom, component 3 is energy). | |
155 | std::valarray<double> four_mom() const; | |
156 | ||
157 | /// returns kt distance (R=1) between this jet and another | |
158 | double kt_distance(const PseudoJet & other) const; | |
159 | ||
160 | /// returns squared cylinder (rap-phi) distance between this jet and another | |
161 | double plain_distance(const PseudoJet & other) const; | |
162 | /// returns squared cylinder (rap-phi) distance between this jet and | |
163 | /// another | |
164 | inline double squared_distance(const PseudoJet & other) const { | |
165 | return plain_distance(other);} | |
166 | ||
167 | /// returns other.phi() - this.phi(), constrained to be in | |
168 | /// range -pi .. pi | |
169 | double delta_phi_to(const PseudoJet & other) const; | |
170 | ||
171 | //// this seemed to compile except if it was used | |
172 | //friend inline double | |
173 | // kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) { | |
174 | // return jet1.kt_distance(jet2);} | |
175 | ||
176 | /// returns distance between this jet and the beam | |
177 | inline double beam_distance() const {return _kt2;} | |
178 | ||
179 | ||
180 | void operator*=(double); | |
181 | void operator/=(double); | |
182 | void operator+=(const PseudoJet &); | |
183 | void operator-=(const PseudoJet &); | |
184 | ||
185 | /// reset the 4-momentum according to the supplied components and | |
186 | /// put the user and history indices back to their default values | |
187 | inline void reset(double px, double py, double pz, double E); | |
188 | ||
189 | /// reset the PseudoJet to be equal to psjet (including its | |
190 | /// indices); NB if the argument is derived from a PseudoJet then | |
191 | /// the "reset" used will be the templated version (which does not | |
192 | /// know about indices...) | |
193 | inline void reset(const PseudoJet & psjet) { | |
194 | (*this) = psjet; | |
195 | } | |
196 | ||
197 | /// reset the 4-momentum according to the supplied generic 4-vector | |
198 | /// (accessible via indexing, [0]==px,...[3]==E) and put the user | |
199 | /// and history indices back to their default values. | |
200 | template <class L> inline void reset(const L & some_four_vector) { | |
201 | reset(some_four_vector[0], some_four_vector[1], | |
202 | some_four_vector[2], some_four_vector[3]); | |
203 | } | |
204 | ||
205 | private: | |
206 | // NB: following order must be kept for things to behave sensibly... | |
207 | double _px,_py,_pz,_E; | |
208 | double _phi, _rap, _kt2; | |
209 | int _cluster_hist_index, _user_index; | |
210 | /// calculate phi, rap, kt2 based on the 4-momentum components | |
211 | void _finish_init(); | |
212 | /// set the indices to default values | |
213 | void _reset_indices(); | |
214 | ||
215 | //vertex_type * vertex0, vertex1; | |
216 | }; | |
217 | ||
218 | ||
219 | //---------------------------------------------------------------------- | |
220 | // routines for basic binary operations | |
221 | ||
222 | PseudoJet operator+(const PseudoJet &, const PseudoJet &); | |
223 | PseudoJet operator-(const PseudoJet &, const PseudoJet &); | |
224 | PseudoJet operator*(double, const PseudoJet &); | |
225 | PseudoJet operator*(const PseudoJet &, double); | |
226 | PseudoJet operator/(const PseudoJet &, double); | |
227 | ||
228 | inline double dot_product(const PseudoJet & a, const PseudoJet & b) { | |
229 | return a.E()*b.E() - a.px()*b.px() - a.py()*b.py() - a.pz()*b.pz(); | |
230 | } | |
231 | ||
232 | /// returns true if the momenta of the two input jets are identical | |
233 | bool have_same_momentum(const PseudoJet &, const PseudoJet &); | |
234 | ||
235 | /// return a pseudojet with the given pt, y, phi and mass | |
236 | PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0); | |
237 | ||
238 | //---------------------------------------------------------------------- | |
239 | // Routines to do with providing sorted arrays of vectors. | |
240 | ||
241 | /// return a vector of jets sorted into decreasing transverse momentum | |
242 | std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets); | |
243 | ||
244 | /// return a vector of jets sorted into increasing rapidity | |
245 | std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets); | |
246 | ||
247 | /// return a vector of jets sorted into decreasing energy | |
248 | std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets); | |
249 | ||
250 | /// return a vector of jets sorted into increasing pz | |
251 | std::vector<PseudoJet> sorted_by_pz(const std::vector<PseudoJet> & jets); | |
252 | ||
253 | //---------------------------------------------------------------------- | |
254 | // some code to help sorting | |
255 | ||
256 | /// sort the indices so that values[indices[0->n-1]] is sorted | |
257 | /// into increasing order | |
258 | void sort_indices(std::vector<int> & indices, | |
259 | const std::vector<double> & values); | |
260 | ||
261 | /// given a vector of values with a one-to-one correspondence with the | |
262 | /// vector of objects, sort objects into an order such that the | |
263 | /// associated values would be in increasing order (but don't actually | |
264 | /// touch the values vector in the process). | |
265 | template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects, | |
266 | const std::vector<double> & values); | |
267 | ||
268 | /// a class that helps us carry out indexed sorting. | |
269 | class IndexedSortHelper { | |
270 | public: | |
271 | inline IndexedSortHelper (const std::vector<double> * reference_values) { | |
272 | _ref_values = reference_values; | |
273 | }; | |
274 | inline int operator() (const int & i1, const int & i2) const { | |
275 | return (*_ref_values)[i1] < (*_ref_values)[i2]; | |
276 | }; | |
277 | private: | |
278 | const std::vector<double> * _ref_values; | |
279 | }; | |
280 | ||
281 | ||
282 | //---------------------------------------------------------------------- | |
283 | /// constructor from any object that has px,py,pz,E = some_four_vector[0--3], | |
284 | // NB: do not know if it really needs to be inline, but when it wasn't | |
285 | // linking failed with g++ (who knows what was wrong...) | |
286 | template <class L> inline PseudoJet::PseudoJet(const L & some_four_vector) { | |
287 | ||
288 | _px = some_four_vector[0]; | |
289 | _py = some_four_vector[1]; | |
290 | _pz = some_four_vector[2]; | |
291 | _E = some_four_vector[3]; | |
292 | _finish_init(); | |
293 | // some default values for these two indices | |
294 | _reset_indices(); | |
295 | } | |
296 | ||
297 | ||
298 | //---------------------------------------------------------------------- | |
299 | inline void PseudoJet::_reset_indices() { | |
300 | set_cluster_hist_index(-1); | |
301 | set_user_index(-1); | |
302 | } | |
303 | ||
304 | //---------------------------------------------------------------------- | |
305 | /// specialization of the "reset" template for case where something | |
306 | /// is reset to a pseudojet -- it then takes the user and history | |
307 | /// indices from the psjet | |
308 | // template<> inline void PseudoJet::reset<PseudoJet>(const PseudoJet & psjet) { | |
309 | // (*this) = psjet; | |
310 | // } | |
311 | ||
312 | ////// fun and games... | |
313 | ////template<class L> class FJVector : public L { | |
314 | ////// /** Default Constructor: create jet with no constituents */ | |
315 | ////// Vector<L>(); | |
316 | //// | |
317 | ////}; | |
318 | //// | |
319 | ||
320 | // taken literally from CLHEP | |
321 | inline double PseudoJet::m() const { | |
322 | double mm = m2(); | |
323 | return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm); | |
324 | } | |
325 | ||
326 | ||
327 | inline void PseudoJet::reset(double px, double py, double pz, double E) { | |
328 | _px = px; | |
329 | _py = py; | |
330 | _pz = pz; | |
331 | _E = E; | |
332 | _finish_init(); | |
333 | _reset_indices(); | |
334 | } | |
335 | ||
336 | ||
337 | } // fastjet namespace | |
338 | ||
339 | #endif // __FASTJET_PSEUDOJET_HH__ |