EIC code displayed by LXR |
|
|
|
|||
File indexing completed on 2025-02-22 10:55:44
0001 // -*- C++ -*- 0002 /////////////////////////////////////////////////////////////////////////////// 0003 // File: momentum.h // 0004 // Description: header file for 4-momentum class Cmomentum // 0005 // This file is part of the SISCone project. // 0006 // WARNING: this is not the main SISCone trunk but // 0007 // an adaptation to spherical coordinates // 0008 // For more details, see http://projects.hepforge.org/siscone // 0009 // // 0010 // Copyright (c) 2006-2008 Gavin Salam and Gregory Soyez // 0011 // // 0012 // This program is free software; you can redistribute it and/or modify // 0013 // it under the terms of the GNU General Public License as published by // 0014 // the Free Software Foundation; either version 2 of the License, or // 0015 // (at your option) any later version. // 0016 // // 0017 // This program is distributed in the hope that it will be useful, // 0018 // but WITHOUT ANY WARRANTY; without even the implied warranty of // 0019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // 0020 // GNU General Public License for more details. // 0021 // // 0022 // You should have received a copy of the GNU General Public License // 0023 // along with this program; if not, write to the Free Software // 0024 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA // 0025 // // 0026 // $Revision:: $// 0027 // $Date:: $// 0028 /////////////////////////////////////////////////////////////////////////////// 0029 0030 #ifndef __SPH_VECTOR_H__ 0031 #define __SPH_VECTOR_H__ 0032 0033 #include <vector> 0034 #include <math.h> 0035 #include <siscone/reference.h> 0036 #include "geom_2d.h" 0037 #include <siscone/defines.h> 0038 0039 namespace siscone_spherical{ 0040 0041 /** 0042 * \class CSph3vector 0043 * \brief base class for managing the spatial part of Cmomentum (defined after) 0044 * 0045 * This class contains the information for particle or group of 0046 * particles management. 0047 * It is adapted to use spherical geometry, where, for our purposes, 0048 * the only time-consuming operation we need is the computation of 0049 * the norm. To compute it once-and-for-all and store it in a local 0050 * variable, you should call the 'build_norm' method. 0051 * On top of that, the angle phi is computed from the x-axis 0052 * and theta from the "north pole". 0053 */ 0054 class CSph3vector{ 0055 public: 0056 /// default ctor 0057 CSph3vector(); 0058 0059 /// ctor with initialisation 0060 CSph3vector(double _px, double _py, double _pz); 0061 0062 /// default dtor 0063 ~CSph3vector(); 0064 0065 /// assignment of vectors 0066 CSph3vector& operator = (const CSph3vector &v); 0067 0068 /// addition of vectors 0069 /// WARNING= norm is not updated 0070 const CSph3vector operator + (const CSph3vector &v); 0071 0072 /// subtraction of vectors 0073 /// WARNING= norm is not updated 0074 const CSph3vector operator - (const CSph3vector &v); 0075 0076 /// division by a constant 0077 /// WARNING= norm is not updated 0078 const CSph3vector operator / (const double &r); 0079 0080 /// incrementation of vectors 0081 /// WARNING= norm is not updated 0082 CSph3vector& operator += (const CSph3vector &v); 0083 0084 /// decrementation of vectors 0085 /// WARNING= norm is not updated 0086 CSph3vector& operator -= (const CSph3vector &v); 0087 0088 /// multiplication by a constant 0089 /// WARNING= norm is not updated 0090 CSph3vector& operator *= (const double &r); 0091 0092 /// division by a constant 0093 /// WARNING= norm is not updated 0094 CSph3vector& operator /= (const double &r); 0095 0096 /// computes pT 0097 inline double perp() const {return sqrt(perp2());} 0098 0099 /// computes pT^2 0100 inline double perp2() const {return px*px+py*py;} 0101 0102 /// 3-vect norm 0103 inline double norm() const {return sqrt(px*px+py*py+pz*pz);} 0104 0105 /// 3-vect norm squared 0106 inline double norm2() const {return px*px+py*py+pz*pz;} 0107 0108 /// 3-vect azimuthal angle 0109 inline double phi() const {return atan2(py, px);} 0110 0111 /// 3-vect polar angle 0112 inline double theta() const {return atan2(perp(),pz);} 0113 0114 /// build the spatial normfrom 4-momentum info 0115 /// !!! WARNING !!! 0116 /// !!! computing the norm is the only time-consuming !!! 0117 /// !!! information we need in all computations. !!! 0118 /// !!! use this whenever you need repeated access !!! 0119 /// !!! to the norm to store it in the local variable !!! 0120 void build_norm(); 0121 0122 /// just a useful tool to store theta and phi 0123 /// locally (in _theta and _phi) in case you need 0124 /// repeated access 0125 void build_thetaphi(); 0126 0127 /// for this direction, compute the two reference directions 0128 /// used to measure angles 0129 void get_angular_directions(CSph3vector &angular_dir1, CSph3vector &angular_dir2); 0130 0131 double px; ///< x-momentum 0132 double py; ///< y-momentum 0133 double pz; ///< z-momentum 0134 0135 double _norm; ///< particle spatial norm (available ONLY after a call to build_norm) 0136 double _theta; ///< particle theta angle (available ONLY after a call to build_thetaphi) 0137 double _phi; ///< particle phi angle (available ONLY after a call to build_thetaphi) 0138 0139 ////////////////////////////////////////////// 0140 // the following part is used for checksums // 0141 ////////////////////////////////////////////// 0142 siscone::Creference ref; ///< reference number for the vector 0143 }; 0144 0145 /** 0146 * \class CSphmomentum 0147 * \brief base class for dynamic coordinates management 0148 * 0149 * This class contains the information for particle or group of 0150 * particles management. 0151 * It is adapted to use spherical geometry, where, for our purposes, 0152 * the only time-consuming operation we need is the computation of 0153 * the norm. To compute it once-and-for-all and store it in a local 0154 * variable, you should call the 'build_norm' method. 0155 * On top of that, the angle phi is computed from the x-axis 0156 * and theta from the "north pole". 0157 */ 0158 class CSphmomentum : public CSph3vector{ 0159 public: 0160 /// default ctor 0161 CSphmomentum(); 0162 0163 /// init from a 3-vect 0164 CSphmomentum(CSph3vector &init, double E=0.0); 0165 0166 /// ctor with initialisation 0167 CSphmomentum(double _px, double _py, double _pz, double _E); 0168 0169 /// ctor with detailed initialisation 0170 //CSphmomentum(double _eta, double _phi, siscone::Creference _ref); 0171 0172 /// default dtor 0173 ~CSphmomentum(); 0174 0175 /// computes m 0176 inline double mass() const {return sqrt(mass2());} 0177 0178 /// computes m^2 0179 inline double mass2() const {return perpmass2()-perp2();} 0180 0181 /// transverse mass, mt = sqrt(pt^2+m^2) = sqrt(E^2 - pz^2) 0182 inline double perpmass() const {return sqrt((E-pz)*(E+pz));} 0183 0184 /// transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2 0185 inline double perpmass2() const {return (E-pz)*(E+pz);} 0186 0187 /// computes transverse energy 0188 inline double Et() const {return E/sqrt(1.0+pz*pz/perp2());} 0189 0190 /// computes transverse energy (squared) 0191 inline double Et2() const {return E*E/(1.0+pz*pz/perp2());} 0192 0193 /// assignment of vectors 0194 CSphmomentum& operator = (const CSphmomentum &v); 0195 0196 /// addition of vectors 0197 /// !!! WARNING !!! no updating of eta and phi !!! 0198 const CSphmomentum operator + (const CSphmomentum &v); 0199 0200 /// incrementation of vectors 0201 /// !!! WARNING !!! no updating of eta and phi !!! 0202 CSphmomentum& operator += (const CSphmomentum &v); 0203 0204 /// decrementation of vectors 0205 /// !!! WARNING !!! no updating of eta and phi !!! 0206 CSphmomentum& operator -= (const CSphmomentum &v); 0207 0208 double E; ///< energy 0209 0210 int parent_index; ///< particle number in the parent list 0211 int index; ///< internal particle number 0212 }; 0213 0214 /// ordering of two vectors 0215 /// this is by default done w.r.t. their references 0216 bool operator < (const CSphmomentum &v1, const CSphmomentum &v2); 0217 0218 /// ordering of vectors in eta (e.g. used in collinear tests) 0219 bool momentum_theta_less(const CSphmomentum &v1, const CSphmomentum &v2); 0220 0221 /// ordering of vectors in pt 0222 bool momentum_pt_less(const CSphmomentum &v1, const CSphmomentum &v2); 0223 0224 0225 ////////////////////////// 0226 // some handy utilities // 0227 ////////////////////////// 0228 0229 /// square 0230 inline double sqr(double x){return x*x;} 0231 0232 /// dot product for te spatial 3-vect 0233 /// \param v1 first 4-vect 0234 /// \param v2 second 4-vect 0235 inline double dot_product3(const CSph3vector &v1, const CSph3vector &v2){ 0236 //double tmp = v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz; 0237 //if (!isfinite(tmp)){ 0238 // std::cout << "dot_product inf: " << std::endl; 0239 // std::cout << " angles: " << v1._theta << " " << v1._phi << " and " << v2._theta << " " << v2._phi << std::endl; 0240 // std::cout << " moms : " << v1.px << " " << v1.py << " " << v1.pz 0241 // << " and " << v2.px << " " << v2.py << " " << v2.pz << std::endl; 0242 //} 0243 return v1.px*v2.px + v1.py*v2.py + v1.pz*v2.pz; 0244 } 0245 0246 /// cross product for the spatial 3-vect 0247 /// \param v1 first 4-vect 0248 /// \param v2 second 4-vect 0249 inline CSph3vector cross_product3(const CSph3vector &v1, const CSph3vector &v2){ 0250 //CSph3vector tmp; 0251 //tmp.px = v1.py*v2.pz-v1.pz*v2.py; 0252 //tmp.py = v1.pz*v2.px-v1.px*v2.pz; 0253 //tmp.pz = v1.px*v2.py-v1.py*v2.px; 0254 //return tmp; 0255 return CSph3vector(v1.py*v2.pz-v1.pz*v2.py, 0256 v1.pz*v2.px-v1.px*v2.pz, 0257 v1.px*v2.py-v1.py*v2.px); 0258 } 0259 0260 /// squared norm of the cross product for the spatial 3-vect (energy is set to 0) 0261 /// \param v1 first 4-vect 0262 /// \param v2 second 4-vect 0263 inline double norm2_cross_product3(const CSph3vector &v1, const CSph3vector &v2){ 0264 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); 0265 } 0266 0267 /// get tangent squared of the spherical distance between two vectors 0268 /// \param v1 vector defining the first point 0269 /// \param v2 vector defining the second point 0270 inline double get_tan2_distance(const CSphmomentum &v1, const CSphmomentum &v2){ 0271 return norm2_cross_product3(v1,v2)/sqr(dot_product3(v1,v2)); 0272 } 0273 0274 /// get spherical distance between to vectors 0275 /// \param v1 vector defining the first point 0276 /// \param v2 vector defining the second point 0277 inline double get_distance(const CSph3vector *v1, const CSph3vector *v2){ 0278 return atan2(sqrt(norm2_cross_product3(*v1,*v2)), dot_product3(*v1,*v2)); 0279 } 0280 0281 /// return true if the two points are distant by less than get spherical distance between two vectors 0282 /// \param v1 vector defining the first point 0283 /// \param v2 vector defining the second point 0284 /// \param tan2R tangent squared of the max distance 0285 /// WARNING: using the tangent here is dangerous for R>pi/2. 0286 /// this never happens per se for "regular R" but 0287 /// it may in the vicinity computation as we're using 0288 /// 2R there. 0289 inline bool is_closer(const CSph3vector *v1, const CSph3vector *v2, const double tan2R){ 0290 double dot = dot_product3(*v1,*v2); 0291 return (dot>=0) && (norm2_cross_product3(*v1,*v2)<=tan2R*dot*dot); 0292 } 0293 0294 /// return true if the two points are distant by less than get spherical distance between to vectors 0295 /// \param v1 vector defining the first point 0296 /// \param v2 vector defining the second point 0297 /// \param tan2R tangent squared of the max distance 0298 /// safer version but computes the norm 0299 inline bool is_closer_safer(const CSph3vector *v1, const CSph3vector *v2, const double cosR){ 0300 return dot_product3(*v1,*v2)>=cosR*sqrt(v1->norm2()*v2->norm2()); 0301 //double dot = dot_product3(*v1,*v2); 0302 //return (dot>=0) && (norm2_cross_product3(*v1,*v2)<tan2R*dot*dot); 0303 } 0304 0305 /// multiply a vector by a constant 0306 /// WARNING: norm not updated 0307 inline CSph3vector operator * (const double &r, const CSph3vector &v){ 0308 CSph3vector tmp = v; 0309 return tmp*=r; 0310 } 0311 } 0312 #endif
| [ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
|
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
|