Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/include/fastjet/PseudoJet.hh@ 227

Last change on this file since 227 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

File size: 11.7 KB
Line 
1//STARTHEADER
2// $Id: PseudoJet.hh,v 1.1 2008-11-06 14:32:08 ovyn Exp $
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 "Utilities/Fastjet/include/fastjet/internal/numconsts.hh"
41
42FASTJET_BEGIN_NAMESPACE // 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.
48const double MaxRap = 1e5;
49
50/// Class to contain pseudojets, including minimal information of use to
51/// to jet-clustering routines.
52class 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 const double phi() const {return phi_02pi();}
74
75 /// returns phi in the range -pi..pi
76 inline const double phi_std() const {
77 return _phi > pi ? _phi-twopi : _phi;}
78
79 /// returns phi in the range 0..2pi
80 inline const 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 /// returns component i, where X==0, Y==1, Z==2, E==3
110 double operator () (int i) const ;
111 /// returns component i, where X==0, Y==1, Z==2, E==3
112 inline double operator [] (int i) const { return (*this)(i); }; // this too
113
114 // taken from CLHEP
115 enum { X=0, Y=1, Z=2, T=3, NUM_COORDINATES=4, SIZE=NUM_COORDINATES };
116
117
118 /// transform this jet (given in the rest frame of prest) into a jet
119 /// in the lab frame [NOT FULLY TESTED]
120 PseudoJet & boost(const PseudoJet & prest);
121 /// transform this jet (given in lab) into a jet in the rest
122 /// frame of ps [NOT FULLY TESTED]
123 PseudoJet & unboost(const PseudoJet & prest);
124
125 /// return the cluster_hist_index, intended to be used by clustering
126 /// routines.
127 inline const int & cluster_hist_index() const {return _cluster_hist_index;}
128 /// set the cluster_hist_index, intended to be used by clustering routines.
129 inline void set_cluster_hist_index(const int index) {_cluster_hist_index = index;}
130
131 /// alternative name for cluster_hist_index() [perhaps more meaningful]
132 inline const int cluster_sequence_history_index() const {
133 return cluster_hist_index();}
134 /// alternative name for set_cluster_hist_index(...) [perhaps more
135 /// meaningful]
136 inline void set_cluster_sequence_history_index(const int index) {
137 set_cluster_hist_index(index);}
138
139
140 /// return the user_index, intended to allow the user to "add" information
141 inline const int & user_index() const {return _user_index;}
142 /// set the user_index, intended to allow the user to "add" information
143 inline void set_user_index(const int index) {_user_index = index;}
144
145 /// return a valarray containing the four-momentum (components 0-2
146 /// are 3-mom, component 3 is energy).
147 std::valarray<double> four_mom() const;
148
149 /// returns kt distance (R=1) between this jet and another
150 double kt_distance(const PseudoJet & other) const;
151
152 /// returns squared cylinder (rap-phi) distance between this jet and another
153 double plain_distance(const PseudoJet & other) const;
154 /// returns squared cylinder (rap-phi) distance between this jet and
155 /// another
156 inline double squared_distance(const PseudoJet & other) const {
157 return plain_distance(other);}
158
159 //// this seemed to compile except if it was used
160 //friend inline double
161 // kt_distance(const PseudoJet & jet1, const PseudoJet & jet2) {
162 // return jet1.kt_distance(jet2);}
163
164 /// returns distance between this jet and the beam
165 inline double beam_distance() const {return _kt2;}
166
167
168 void operator*=(double);
169 void operator/=(double);
170 void operator+=(const PseudoJet &);
171 void operator-=(const PseudoJet &);
172
173 /// reset the 4-momentum according to the supplied components and
174 /// put the user and history indices back to their default values
175 inline void reset(double px, double py, double pz, double E);
176
177 /// reset the PseudoJet to be equal to psjet (including its
178 /// indices); NB if the argument is derived from a PseudoJet then
179 /// the "reset" used will be the templated version (which does not
180 /// know about indices...)
181 inline void reset(const PseudoJet & psjet) {
182 (*this) = psjet;
183 }
184
185 /// reset the 4-momentum according to the supplied generic 4-vector
186 /// (accessible via indexing, [0]==px,...[3]==E) and put the user
187 /// and history indices back to their default values.
188 template <class L> inline void reset(const L & some_four_vector) {
189 reset(some_four_vector[0], some_four_vector[1],
190 some_four_vector[2], some_four_vector[3]);
191 }
192
193 private:
194 // NB: following order must be kept for things to behave sensibly...
195 double _px,_py,_pz,_E;
196 double _phi, _rap, _kt2;
197 int _cluster_hist_index, _user_index;
198 /// calculate phi, rap, kt2 based on the 4-momentum components
199 void _finish_init();
200 /// set the indices to default values
201 void _reset_indices();
202
203 //vertex_type * vertex0, vertex1;
204};
205
206
207//----------------------------------------------------------------------
208// routines for basic binary operations
209
210PseudoJet operator+(const PseudoJet &, const PseudoJet &);
211PseudoJet operator-(const PseudoJet &, const PseudoJet &);
212PseudoJet operator*(double, const PseudoJet &);
213PseudoJet operator*(const PseudoJet &, double);
214PseudoJet operator/(const PseudoJet &, double);
215
216/// returns true if the momenta of the two input jets are identical
217bool have_same_momentum(const PseudoJet &, const PseudoJet &);
218
219/// return a pseudojet with the given pt, y, phi and mass
220PseudoJet PtYPhiM(double pt, double y, double phi, double m = 0.0);
221
222//----------------------------------------------------------------------
223// Routines to do with providing sorted arrays of vectors.
224
225/// return a vector of jets sorted into decreasing transverse momentum
226std::vector<PseudoJet> sorted_by_pt(const std::vector<PseudoJet> & jets);
227
228/// return a vector of jets sorted into increasing rapidity
229std::vector<PseudoJet> sorted_by_rapidity(const std::vector<PseudoJet> & jets);
230
231/// return a vector of jets sorted into decreasing energy
232std::vector<PseudoJet> sorted_by_E(const std::vector<PseudoJet> & jets);
233
234//----------------------------------------------------------------------
235// some code to help sorting
236
237/// sort the indices so that values[indices[0->n-1]] is sorted
238/// into increasing order
239void sort_indices(std::vector<int> & indices,
240 const std::vector<double> & values);
241
242/// given a vector of values with a one-to-one correspondence with the
243/// vector of objects, sort objects into an order such that the
244/// associated values would be in increasing order (but don't actually
245/// touch the values vector in the process).
246template<class T> std::vector<T> objects_sorted_by_values(const std::vector<T> & objects,
247 const std::vector<double> & values);
248
249/// a class that helps us carry out indexed sorting.
250class IndexedSortHelper {
251public:
252 inline IndexedSortHelper (const std::vector<double> * reference_values) {
253 _ref_values = reference_values;
254 };
255 inline int operator() (const int & i1, const int & i2) const {
256 return (*_ref_values)[i1] < (*_ref_values)[i2];
257 };
258private:
259 const std::vector<double> * _ref_values;
260};
261
262
263//----------------------------------------------------------------------
264/// constructor from any object that has px,py,pz,E = some_four_vector[0--3],
265// NB: do not know if it really needs to be inline, but when it wasn't
266// linking failed with g++ (who knows what was wrong...)
267template <class L> inline PseudoJet::PseudoJet(const L & some_four_vector) {
268
269 _px = some_four_vector[0];
270 _py = some_four_vector[1];
271 _pz = some_four_vector[2];
272 _E = some_four_vector[3];
273 _finish_init();
274 // some default values for these two indices
275 _reset_indices();
276}
277
278
279//----------------------------------------------------------------------
280inline void PseudoJet::_reset_indices() {
281 set_cluster_hist_index(-1);
282 set_user_index(-1);
283}
284
285//----------------------------------------------------------------------
286/// specialization of the "reset" template for case where something
287/// is reset to a pseudojet -- it then takes the user and history
288/// indices from the psjet
289// template<> inline void PseudoJet::reset<PseudoJet>(const PseudoJet & psjet) {
290// (*this) = psjet;
291// }
292
293////// fun and games...
294////template<class L> class FJVector : public L {
295////// /** Default Constructor: create jet with no constituents */
296////// Vector<L>();
297////
298////};
299////
300
301// taken literally from CLHEP
302inline double PseudoJet::m() const {
303 double mm = m2();
304 return mm < 0.0 ? -std::sqrt(-mm) : std::sqrt(mm);
305}
306
307
308inline void PseudoJet::reset(double px, double py, double pz, double E) {
309 _px = px;
310 _py = py;
311 _pz = pz;
312 _E = E;
313 _finish_init();
314 _reset_indices();
315}
316
317
318FASTJET_END_NAMESPACE
319
320#endif // __FASTJET_PSEUDOJET_HH__
Note: See TracBrowser for help on using the repository browser.