[u/mrichter/AliRoot.git] / JETAN / fastjet / siscone / spherical / momentum.h

// -*- C++ -*-
///////////////////////////////////////////////////////////////////////////////
// File: momentum.h                                                          //
// Description: header file for 4-momentum class Cmomentum                   //
// This file is part of the SISCone project.                                 //
// WARNING: this is not the main SISCone trunk but                           //
//          an adaptation to spherical coordinates                           //
// For more details, see http://projects.hepforge.org/siscone                //
//                                                                           //
// Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez                     //
//                                                                           //
// This program is free software; you can redistribute it and/or modify      //
// it under the terms of the GNU General Public License as published by      //
// the Free Software Foundation; either version 2 of the License, or         //
// (at your option) any later version.                                       //
//                                                                           //
// This program is distributed in the hope that it will be useful,           //
// but WITHOUT ANY WARRANTY; without even the implied warranty of            //
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             //
// GNU General Public License for more details.                              //
//                                                                           //
// You should have received a copy of the GNU General Public License         //
// along with this program; if not, write to the Free Software               //
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
//                                                                           //
// $Revision:: 256                                                          $//
// $Date:: 2008-07-14 13:52:16 +0200 (Mon, 14 Jul 2008)                     $//
///////////////////////////////////////////////////////////////////////////////

#ifndef __SPH_VECTOR_H__
#define __SPH_VECTOR_H__

#include <vector>
#include <math.h>
#include <siscone/reference.h>
#include "geom_2d.h"
#include <siscone/defines.h>

namespace siscone_spherical{

/**
 * \class CSph3vector
 * \brief base class for managing the spatial part of Cmomentum (defined after)
 *
 * This class contains the information for particle or group of
 * particles management.
 * It is adapted to use spherical geometry, where, for our purposes,
 * the only time-consuming operation we need is the computation of 
 * the norm. To compute it once-and-for-all and store it in a local 
 * variable, you should call the 'build_norm' method.
 * On top of that, the angle phi is computed from the x-axis
 * and theta from the "north pole". 
 */
class CSph3vector{
 public:
  /// default ctor
  CSph3vector();

  /// ctor with initialisation
  CSph3vector(double _px, double _py, double _pz);

  /// default dtor
  ~CSph3vector();

  /// assignment of vectors
  CSph3vector& operator = (const CSph3vector &v);

  /// addition of vectors
  /// WARNING= norm is not updated
  const CSph3vector operator + (const CSph3vector &v);

  /// subtraction of vectors
  /// WARNING= norm is not updated
  const CSph3vector operator - (const CSph3vector &v);

  /// division by a constant
  /// WARNING= norm is not updated
  const CSph3vector operator / (const double &r);

  /// incrementation of vectors
  /// WARNING= norm is not updated
  CSph3vector& operator += (const CSph3vector &v);

  /// decrementation of vectors
  /// WARNING= norm is not updated
  CSph3vector& operator -= (const CSph3vector &v);

  /// multiplication by a constant
  /// WARNING= norm is not updated
  CSph3vector& operator *= (const double &r);

  /// division by a constant
  /// WARNING= norm is not updated
  CSph3vector& operator /= (const double &r);

  /// computes pT
  inline double perp() const {return sqrt(perp2());}

  /// computes pT^2
  inline double perp2() const {return px*px+py*py;}

  /// 3-vect norm
  inline double norm() const {return sqrt(px*px+py*py+pz*pz);}

  /// 3-vect norm squared
  inline double norm2() const {return px*px+py*py+pz*pz;}

  /// 3-vect azimuthal angle
  inline double phi() const {return atan2(py, px);}

  /// 3-vect polar angle
  inline double theta() const {return atan2(perp(),pz);}

  /// build the spatial normfrom 4-momentum info
  /// !!!                  WARNING                       !!!
  /// !!! computing the norm is the only time-consuming  !!!
  /// !!! information we need in all computations.       !!!
  /// !!! use this whenever you need repeated access     !!!
  /// !!! to the norm to store it in the local variable  !!!
  void build_norm();

  /// just a useful tool to store theta and phi
  /// locally (in _theta and _phi) in case you need 
  /// repeated access
  void build_thetaphi();

  /// for this direction, compute the two reference directions
  /// used to measure angles
  void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2);

  double px;        ///< x-momentum
  double py;        ///< y-momentum
  double pz;        ///< z-momentum

  double _norm;     ///< particle spatial norm (available ONLY after a call to build_norm)
  double _theta;    ///< particle theta angle (available ONLY after a call to build_thetaphi)
  double _phi;      ///< particle phi angle   (available ONLY after a call to build_thetaphi)

  //////////////////////////////////////////////
  // the following part is used for checksums //
  //////////////////////////////////////////////
  siscone::Creference ref;   ///< reference number for the vector
};

/**
 * \class CSphmomentum
 * \brief base class for dynamic coordinates management
 *
 * This class contains the information for particle or group of
 * particles management.
 * It is adapted to use spherical geometry, where, for our purposes,
 * the only time-consuming operation we need is the computation of 
 * the norm. To compute it once-and-for-all and store it in a local 
 * variable, you should call the 'build_norm' method.
 * On top of that, the angle phi is computed from the x-axis
 * and theta from the "north pole". 
 */
class CSphmomentum : public CSph3vector{
 public:
  /// default ctor
  CSphmomentum();

  /// init from a 3-vect
  CSphmomentum(CSph3vector &init, double E=0.0);

  /// ctor with initialisation
  CSphmomentum(double _px, double _py, double _pz, double _E);

  /// ctor with detailed initialisation
  //CSphmomentum(double _eta, double _phi, siscone::Creference _ref);

  /// default dtor
  ~CSphmomentum();

  /// computes m
  inline double mass() const {return sqrt(mass2());}

  /// computes m^2
  inline double mass2() const {return perpmass2()-perp2();}

  /// transverse mass, mt = sqrt(pt^2+m^2) = sqrt(E^2 - pz^2)
  inline double perpmass() const {return sqrt((E-pz)*(E+pz));}

  /// transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2
  inline double perpmass2() const {return (E-pz)*(E+pz);}

  /// computes transverse energy
  inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());}

  /// computes transverse energy (squared)
  inline double Et2() const {return E*E/(1.0+pz*pz/perp2());}

  /// assignment of vectors
  CSphmomentum& operator = (const CSphmomentum &v);

  /// addition of vectors
  /// !!! WARNING !!! no updating of eta and phi !!!
  const CSphmomentum operator + (const CSphmomentum &v);

  /// incrementation of vectors
  /// !!! WARNING !!! no updating of eta and phi !!!
  CSphmomentum& operator += (const CSphmomentum &v);

  /// decrementation of vectors
  /// !!! WARNING !!! no updating of eta and phi !!!
  CSphmomentum& operator -= (const CSphmomentum &v);

  double E;         ///< energy

  int parent_index; ///< particle number in the parent list
  int index;        ///< internal particle number
};

/// ordering of two vectors
/// this is by default done w.r.t. their references
bool operator < (const CSphmomentum &v1, const CSphmomentum &v2);

/// ordering of vectors in eta (e.g. used in collinear tests)
bool momentum_theta_less(const CSphmomentum &v1, const CSphmomentum &v2);

/// ordering of vectors in pt
bool momentum_pt_less(const CSphmomentum &v1, const CSphmomentum &v2);


//////////////////////////
// some handy utilities //
//////////////////////////

/// square
inline double sqr(double x){return x*x;}

/// dot product for te spatial 3-vect
/// \param v1    first 4-vect
/// \param v2    second 4-vect
inline double dot_product3(const CSph3vector &v1, const CSph3vector &v2){
  //double tmp = v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz;
  //if (!isfinite(tmp)){
  //  std::cout << "dot_product inf: " << std::endl;
  //  std::cout << "  angles: " << v1._theta << " " << v1._phi << " and " << v2._theta << " " << v2._phi << std::endl; 
  //  std::cout << "  moms  : " << v1.px << " " << v1.py << " " << v1.pz 
  //	      << " and "      << v2.px << " " << v2.py << " " << v2.pz << std::endl;
  //}
  return v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz;
}

/// cross product for the spatial 3-vect
/// \param v1    first 4-vect
/// \param v2    second 4-vect
inline CSph3vector cross_product3(const CSph3vector &v1, const CSph3vector &v2){
  //CSph3vector tmp;
  //tmp.px = v1.py*v2.pz-v1.pz*v2.py;
  //tmp.py = v1.pz*v2.px-v1.px*v2.pz;
  //tmp.pz = v1.px*v2.py-v1.py*v2.px;
  //return tmp;
  return CSph3vector(v1.py*v2.pz-v1.pz*v2.py,
		  v1.pz*v2.px-v1.px*v2.pz,
		  v1.px*v2.py-v1.py*v2.px);
}

/// squared norm of the cross product for the spatial 3-vect (energy is set to 0)
/// \param v1    first 4-vect
/// \param v2    second 4-vect
inline double norm2_cross_product3(const CSph3vector &v1, const CSph3vector &v2){
  return sqr(v1.py*v2.pz-v1.pz*v2.py) + sqr(v1.pz*v2.px-v1.px*v2.pz) + sqr(v1.px*v2.py-v1.py*v2.px);
}

/// get tangent squared of the spherical distance between two vectors
/// \param v1    vector defining the first point
/// \param v2    vector defining the second point
inline double get_tan2_distance(const CSphmomentum &v1, const CSphmomentum &v2){
  return norm2_cross_product3(v1,v2)/sqr(dot_product3(v1,v2));
}

/// get spherical distance between to vectors
/// \param v1    vector defining the first point
/// \param v2    vector defining the second point
inline double get_distance(const CSph3vector *v1, const CSph3vector *v2){
  return atan2(sqrt(norm2_cross_product3(*v1,*v2)), dot_product3(*v1,*v2));
}

/// return true if the two points are distant by less than get spherical distance between two vectors
/// \param v1      vector defining the first point
/// \param v2      vector defining the second point
/// \param tan2R   tangent squared of the max distance
/// WARNING: using the tangent here is dangerous for R>pi/2.
///          this never happens per se for "regular R" but 
///          it may in the vicinity computation as we're using
///          2R there. 
inline bool is_closer(const CSph3vector *v1, const CSph3vector *v2, const double tan2R){
  double dot = dot_product3(*v1,*v2);
  return (dot>=0) && (norm2_cross_product3(*v1,*v2)<=tan2R*dot*dot);
}

/// return true if the two points are distant by less than  get spherical distance between to vectors
/// \param v1      vector defining the first point
/// \param v2      vector defining the second point
/// \param tan2R   tangent squared of the max distance
/// safer version but computes the norm
inline bool is_closer_safer(const CSph3vector *v1, const CSph3vector *v2, const double cosR){
  return dot_product3(*v1,*v2)>=cosR*sqrt(v1->norm2()*v2->norm2());
  //double dot = dot_product3(*v1,*v2);
  //return (dot>=0) && (norm2_cross_product3(*v1,*v2)<tan2R*dot*dot);
}

/// multiply a vector by a constant
/// WARNING: norm not updated
inline CSph3vector operator * (const double &r, const CSph3vector &v){
  CSph3vector tmp = v;
  return tmp*=r;
}
}
#endif
Commit	Line	Data
370be031	1	// -- C++ --
	2	///////////////////////////////////////////////////////////////////////////////
	3	// File: momentum.h //
	4	// Description: header file for 4-momentum class Cmomentum //
	5	// This file is part of the SISCone project. //
	6	// WARNING: this is not the main SISCone trunk but //
	7	// an adaptation to spherical coordinates //
	8	// For more details, see http://projects.hepforge.org/siscone //
	9	// //
	10	// Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez //
	11	// //
	12	// This program is free software; you can redistribute it and/or modify //
	13	// it under the terms of the GNU General Public License as published by //
	14	// the Free Software Foundation; either version 2 of the License, or //
	15	// (at your option) any later version. //
	16	// //
	17	// This program is distributed in the hope that it will be useful, //
	18	// but WITHOUT ANY WARRANTY; without even the implied warranty of //
	19	// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
	20	// GNU General Public License for more details. //
	21	// //
	22	// You should have received a copy of the GNU General Public License //
	23	// along with this program; if not, write to the Free Software //
	24	// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA //
	25	// //
	26	// $Revision:: 256 $//
	27	// $Date:: 2008-07-14 13:52:16 +0200 (Mon, 14 Jul 2008) $//
	28	///////////////////////////////////////////////////////////////////////////////
	29
	30	#ifndef __SPH_VECTOR_H__
	31	#define __SPH_VECTOR_H__
	32
	33	#include <vector>
	34	#include <math.h>
	35	#include <siscone/reference.h>
	36	#include "geom_2d.h"
	37	#include <siscone/defines.h>
	38
	39	namespace siscone_spherical{
	40
	41	/**
	42	* \class CSph3vector
	43	* \brief base class for managing the spatial part of Cmomentum (defined after)
	44	*
	45	* This class contains the information for particle or group of
	46	* particles management.
	47	* It is adapted to use spherical geometry, where, for our purposes,
	48	* the only time-consuming operation we need is the computation of
	49	* the norm. To compute it once-and-for-all and store it in a local
	50	* variable, you should call the 'build_norm' method.
	51	* On top of that, the angle phi is computed from the x-axis
	52	* and theta from the "north pole".
	53	*/
	54	class CSph3vector{
	55	public:
	56	/// default ctor
	57	CSph3vector();
	58
	59	/// ctor with initialisation
	60	CSph3vector(double _px, double _py, double _pz);
	61
	62	/// default dtor
	63	~CSph3vector();
	64
65	/// assignment of vectors
66	CSph3vector& operator = (const CSph3vector &v);
67
68	/// addition of vectors
69	/// WARNING= norm is not updated
70	const CSph3vector operator + (const CSph3vector &v);
71
72	/// subtraction of vectors
73	/// WARNING= norm is not updated
74	const CSph3vector operator - (const CSph3vector &v);
75
76	/// division by a constant
77	/// WARNING= norm is not updated
78	const CSph3vector operator / (const double &r);
79
80	/// incrementation of vectors
81	/// WARNING= norm is not updated
82	CSph3vector& operator += (const CSph3vector &v);
83
84	/// decrementation of vectors
85	/// WARNING= norm is not updated
86	CSph3vector& operator -= (const CSph3vector &v);
87
88	/// multiplication by a constant
89	/// WARNING= norm is not updated
90	CSph3vector& operator *= (const double &r);
91
92	/// division by a constant
93	/// WARNING= norm is not updated
94	CSph3vector& operator /= (const double &r);
95
96	/// computes pT
97	inline double perp() const {return sqrt(perp2());}
98
99	/// computes pT^2
100	inline double perp2() const {return pxpx+pypy;}
101
102	/// 3-vect norm
103	inline double norm() const {return sqrt(pxpx+pypy+pz*pz);}
104
105	/// 3-vect norm squared
106	inline double norm2() const {return pxpx+pypy+pz*pz;}
107
108	/// 3-vect azimuthal angle
109	inline double phi() const {return atan2(py, px);}
110
111	/// 3-vect polar angle
112	inline double theta() const {return atan2(perp(),pz);}
113
114	/// build the spatial normfrom 4-momentum info
115	/// !!! WARNING !!!
116	/// !!! computing the norm is the only time-consuming !!!
117	/// !!! information we need in all computations. !!!
118	/// !!! use this whenever you need repeated access !!!
119	/// !!! to the norm to store it in the local variable !!!
120	void build_norm();
121
122	/// just a useful tool to store theta and phi
123	/// locally (in _theta and _phi) in case you need
124	/// repeated access
125	void build_thetaphi();
126
127	/// for this direction, compute the two reference directions
128	/// used to measure angles
129	void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2);
130
131	double px; ///< x-momentum
132	double py; ///< y-momentum
133	double pz; ///< z-momentum
134
135	double _norm; ///< particle spatial norm (available ONLY after a call to build_norm)
136	double _theta; ///< particle theta angle (available ONLY after a call to build_thetaphi)
137	double _phi; ///< particle phi angle (available ONLY after a call to build_thetaphi)
138
139	//////////////////////////////////////////////
140	// the following part is used for checksums //
141	//////////////////////////////////////////////
142	siscone::Creference ref; ///< reference number for the vector
143	};
144
145	/**
146	* \class CSphmomentum
147	* \brief base class for dynamic coordinates management
148	*
149	* This class contains the information for particle or group of
150	* particles management.
151	* It is adapted to use spherical geometry, where, for our purposes,
152	* the only time-consuming operation we need is the computation of
153	* the norm. To compute it once-and-for-all and store it in a local
154	* variable, you should call the 'build_norm' method.
155	* On top of that, the angle phi is computed from the x-axis
156	* and theta from the "north pole".
157	*/
158	class CSphmomentum : public CSph3vector{
159	public:
160	/// default ctor
161	CSphmomentum();
162
163	/// init from a 3-vect
164	CSphmomentum(CSph3vector &init, double E=0.0);
165
166	/// ctor with initialisation
167	CSphmomentum(double _px, double _py, double _pz, double _E);
168
169	/// ctor with detailed initialisation
170	//CSphmomentum(double _eta, double _phi, siscone::Creference _ref);
171
172	/// default dtor
173	~CSphmomentum();
174
175	/// computes m
176	inline double mass() const {return sqrt(mass2());}
177
178	/// computes m^2
179	inline double mass2() const {return perpmass2()-perp2();}
180
181	/// transverse mass, mt = sqrt(pt^2+m^2) = sqrt(E^2 - pz^2)
182	inline double perpmass() const {return sqrt((E-pz)*(E+pz));}
183
184	/// transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2
185	inline double perpmass2() const {return (E-pz)*(E+pz);}
186
187	/// computes transverse energy
188	inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());}
189
190	/// computes transverse energy (squared)
191	inline double Et2() const {return EE/(1.0+pzpz/perp2());}
192
193	/// assignment of vectors
194	CSphmomentum& operator = (const CSphmomentum &v);
195
196	/// addition of vectors
197	/// !!! WARNING !!! no updating of eta and phi !!!
198	const CSphmomentum operator + (const CSphmomentum &v);
199
200	/// incrementation of vectors
201	/// !!! WARNING !!! no updating of eta and phi !!!
202	CSphmomentum& operator += (const CSphmomentum &v);
203
204	/// decrementation of vectors
205	/// !!! WARNING !!! no updating of eta and phi !!!
206	CSphmomentum& operator -= (const CSphmomentum &v);
207
208	double E; ///< energy
209
210	int parent_index; ///< particle number in the parent list
211	int index; ///< internal particle number
212	};
213
214	/// ordering of two vectors
215	/// this is by default done w.r.t. their references
216	bool operator < (const CSphmomentum &v1, const CSphmomentum &v2);
217
218	/// ordering of vectors in eta (e.g. used in collinear tests)
219	bool momentum_theta_less(const CSphmomentum &v1, const CSphmomentum &v2);
220
221	/// ordering of vectors in pt
222	bool momentum_pt_less(const CSphmomentum &v1, const CSphmomentum &v2);
223
224
225	//////////////////////////
226	// some handy utilities //
227	//////////////////////////
228
229	/// square
230	inline double sqr(double x){return x*x;}
231
232	/// dot product for te spatial 3-vect
233	/// \param v1 first 4-vect
234	/// \param v2 second 4-vect
235	inline double dot_product3(const CSph3vector &v1, const CSph3vector &v2){
236	//double tmp = v1.pxv2.px + v1.pyv2.py + v1.pz*v2.pz;
237	//if (!isfinite(tmp)){
238	// std::cout << "dot_product inf: " << std::endl;
239	// std::cout << " angles: " << v1._theta << " " << v1._phi << " and " << v2._theta << " " << v2._phi << std::endl;
240	// std::cout << " moms : " << v1.px << " " << v1.py << " " << v1.pz
241	// << " and " << v2.px << " " << v2.py << " " << v2.pz << std::endl;
242	//}
243	return v1.pxv2.px + v1.pyv2.py + v1.pz*v2.pz;
244	}
245
246	/// cross product for the spatial 3-vect
247	/// \param v1 first 4-vect
248	/// \param v2 second 4-vect
249	inline CSph3vector cross_product3(const CSph3vector &v1, const CSph3vector &v2){
250	//CSph3vector tmp;
251	//tmp.px = v1.pyv2.pz-v1.pzv2.py;
252	//tmp.py = v1.pzv2.px-v1.pxv2.pz;
253	//tmp.pz = v1.pxv2.py-v1.pyv2.px;
254	//return tmp;
255	return CSph3vector(v1.pyv2.pz-v1.pzv2.py,
256	v1.pzv2.px-v1.pxv2.pz,
257	v1.pxv2.py-v1.pyv2.px);
258	}
259
260	/// squared norm of the cross product for the spatial 3-vect (energy is set to 0)
261	/// \param v1 first 4-vect
262	/// \param v2 second 4-vect
263	inline double norm2_cross_product3(const CSph3vector &v1, const CSph3vector &v2){
264	return sqr(v1.pyv2.pz-v1.pzv2.py) + sqr(v1.pzv2.px-v1.pxv2.pz) + sqr(v1.pxv2.py-v1.pyv2.px);
265	}
266
267	/// get tangent squared of the spherical distance between two vectors
268	/// \param v1 vector defining the first point
269	/// \param v2 vector defining the second point
270	inline double get_tan2_distance(const CSphmomentum &v1, const CSphmomentum &v2){
271	return norm2_cross_product3(v1,v2)/sqr(dot_product3(v1,v2));
272	}
273
274	/// get spherical distance between to vectors
275	/// \param v1 vector defining the first point
276	/// \param v2 vector defining the second point
277	inline double get_distance(const CSph3vector v1, const CSph3vector v2){
278	return atan2(sqrt(norm2_cross_product3(v1,v2)), dot_product3(v1,v2));
279	}
280
281	/// return true if the two points are distant by less than get spherical distance between two vectors
282	/// \param v1 vector defining the first point
283	/// \param v2 vector defining the second point
284	/// \param tan2R tangent squared of the max distance
285	/// WARNING: using the tangent here is dangerous for R>pi/2.
286	/// this never happens per se for "regular R" but
287	/// it may in the vicinity computation as we're using
288	/// 2R there.
289	inline bool is_closer(const CSph3vector v1, const CSph3vector v2, const double tan2R){
290	double dot = dot_product3(v1,v2);
291	return (dot>=0) && (norm2_cross_product3(v1,v2)<=tan2Rdotdot);
292	}
293
294	/// return true if the two points are distant by less than get spherical distance between to vectors
295	/// \param v1 vector defining the first point
296	/// \param v2 vector defining the second point
297	/// \param tan2R tangent squared of the max distance
298	/// safer version but computes the norm
299	inline bool is_closer_safer(const CSph3vector v1, const CSph3vector v2, const double cosR){
300	return dot_product3(v1,v2)>=cosRsqrt(v1->norm2()v2->norm2());
301	//double dot = dot_product3(v1,v2);
302	//return (dot>=0) && (norm2_cross_product3(v1,v2)<tan2Rdotdot);
303	}
304
305	/// multiply a vector by a constant
306	/// WARNING: norm not updated
307	inline CSph3vector operator * (const double &r, const CSph3vector &v){
308	CSph3vector tmp = v;
309	return tmp*=r;
310	}
311	}
312	#endif