CLHEP Reference Documentation

00001 // -*- C++ -*-
00002 // $Id: BasicVector3D.cc,v 1.3 2003/08/13 20:00:11 garren Exp $
00003 // ---------------------------------------------------------------------------
00004 
00005 #include <math.h>
00006 #include <float.h>
00007 #include <iostream>
00008 #include "CLHEP/Geometry/defs.h"
00009 #include "CLHEP/Geometry/BasicVector3D.h"
00010 
00011 namespace HepGeom {
00012   //--------------------------------------------------------------------------
00013   template<>
00014   float BasicVector3D<float>::pseudoRapidity() const {
00015     float ma = mag(), dz = z();
00016     if (ma ==  0)  return  0;
00017     if (ma ==  dz) return  FLT_MAX;
00018     if (ma == -dz) return -FLT_MAX;
00019     return 0.5*std::log((ma+dz)/(ma-dz));
00020   }
00021 
00022   //--------------------------------------------------------------------------
00023   template<>
00024   void BasicVector3D<float>::setEta(float a) {
00025     double ma = mag();
00026     if (ma == 0) return;
00027     double tanHalfTheta  = std::exp(-a);
00028     double tanHalfTheta2 = tanHalfTheta * tanHalfTheta;
00029     double cosTheta1      = (1 - tanHalfTheta2) / (1 + tanHalfTheta2);
00030     double rh            = ma * std::sqrt(1 - cosTheta1*cosTheta1);
00031     double ph            = phi();
00032     set(rh*std::cos(ph), rh*std::sin(ph), ma*cosTheta1);
00033   }
00034 
00035   //--------------------------------------------------------------------------
00036   template<>
00037   float BasicVector3D<float>::angle(const BasicVector3D<float> & v) const {
00038     double cosa = 0;
00039     double ptot = mag()*v.mag();
00040     if(ptot > 0) {
00041       cosa = dot(v)/ptot;
00042       if(cosa >  1) cosa =  1;
00043       if(cosa < -1) cosa = -1;
00044     }
00045     return std::acos(cosa);
00046   }
00047 
00048   //--------------------------------------------------------------------------
00049   template<>
00050   BasicVector3D<float> & BasicVector3D<float>::rotateX(float a) {
00051     double sina = std::sin(a), cosa = std::cos(a), dy = y(), dz = z();
00052     setY(dy*cosa-dz*sina);
00053     setZ(dz*cosa+dy*sina);
00054     return *this;
00055   }
00056 
00057   //--------------------------------------------------------------------------
00058   template<>
00059   BasicVector3D<float> & BasicVector3D<float>::rotateY(float a) {
00060     double sina = std::sin(a), cosa = std::cos(a), dz = z(), dx = x();
00061     setZ(dz*cosa-dx*sina);
00062     setX(dx*cosa+dz*sina);
00063     return *this;
00064   }
00065 
00066   //--------------------------------------------------------------------------
00067   template<>
00068   BasicVector3D<float> & BasicVector3D<float>::rotateZ(float a) {
00069     double sina = std::sin(a), cosa = std::cos(a), dx = x(), dy = y();
00070     setX(dx*cosa-dy*sina);
00071     setY(dy*cosa+dx*sina);
00072     return *this;
00073   }
00074 
00075   //--------------------------------------------------------------------------
00076   template<>
00077   BasicVector3D<float> &
00078   BasicVector3D<float>::rotate(float a, const BasicVector3D<float> & v) {
00079     if (a  == 0) return *this;
00080     double cx = v.x(), cy = v.y(), cz = v.z();
00081     double ll = std::sqrt(cx*cx + cy*cy + cz*cz);
00082     if (ll == 0) {
00083       std::cerr << "BasicVector<float>::rotate() : zero axis" << std::endl;
00084       return *this;
00085     }
00086     double cosa = std::cos(a), sina = std::sin(a);
00087     cx /= ll; cy /= ll; cz /= ll;   
00088 
00089     double xx = cosa + (1-cosa)*cx*cx;
00090     double xy =        (1-cosa)*cx*cy - sina*cz;
00091     double xz =        (1-cosa)*cx*cz + sina*cy;
00092     
00093     double yx =        (1-cosa)*cy*cx + sina*cz;
00094     double yy = cosa + (1-cosa)*cy*cy; 
00095     double yz =        (1-cosa)*cy*cz - sina*cx;
00096     
00097     double zx =        (1-cosa)*cz*cx - sina*cy;
00098     double zy =        (1-cosa)*cz*cy + sina*cx;
00099     double zz = cosa + (1-cosa)*cz*cz;
00100 
00101     cx = x(); cy = y(); cz = z();   
00102     set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz);
00103     return *this;
00104   }
00105 
00106   //--------------------------------------------------------------------------
00107   std::ostream &
00108   operator<<(std::ostream & os, const BasicVector3D<float> & a)
00109   {
00110     return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")";
00111   }
00112 
00113   //--------------------------------------------------------------------------
00114   std::istream &
00115   operator>> (std::istream & is, BasicVector3D<float> & a)
00116   {
00117     // Required format is ( a, b, c ) that is, three numbers, preceded by
00118     // (, followed by ), and separated by commas.  The three numbers are
00119     // taken as x, y, z.
00120 
00121     float x, y, z;
00122     char c;
00123 
00124     is >> std::ws >> c;
00125     // ws is defined to invoke eatwhite(istream & )
00126     // see (Stroustrup gray book) page 333 and 345.
00127     if (is.fail() || c != '(' ) {
00128       std::cerr
00129         << "Could not find required opening parenthesis "
00130         << "in input of a BasicVector3D<float>"
00131         << std::endl;
00132       return is;
00133     }
00134 
00135     is >> x >> std::ws >> c;
00136     if (is.fail() || c != ',' ) {
00137       std::cerr
00138         << "Could not find x value and required trailing comma "
00139         << "in input of a BasicVector3D<float>"
00140         << std::endl; 
00141       return is;
00142     }
00143 
00144     is >> y >> std::ws >> c;
00145     if (is.fail() || c != ',' ) {
00146       std::cerr
00147         << "Could not find y value and required trailing comma "
00148         <<  "in input of a BasicVector3D<float>"
00149         << std::endl;
00150       return is;
00151     }
00152 
00153     is >> z >> std::ws >> c;
00154     if (is.fail() || c != ')' ) {
00155       std::cerr
00156         << "Could not find z value and required close parenthesis "
00157         << "in input of a BasicVector3D<float>"
00158         << std::endl;
00159       return is;
00160     }
00161 
00162     a.setX(x);
00163     a.setY(y);
00164     a.setZ(z);
00165     return is;
00166   }
00167   
00168   //--------------------------------------------------------------------------
00169   template<>
00170   double BasicVector3D<double>::pseudoRapidity() const {
00171     double ma = mag(), dz = z();
00172     if (ma ==  0)  return  0;
00173     if (ma ==  dz) return  DBL_MAX;
00174     if (ma == -dz) return -DBL_MAX;
00175     return 0.5*std::log((ma+dz)/(ma-dz));
00176   }
00177 
00178   //--------------------------------------------------------------------------
00179   template<>
00180   void BasicVector3D<double>::setEta(double a) {
00181     double ma = mag();
00182     if (ma == 0) return;
00183     double tanHalfTheta  = std::exp(-a);
00184     double tanHalfTheta2 = tanHalfTheta * tanHalfTheta;
00185     double cosTheta1      = (1 - tanHalfTheta2) / (1 + tanHalfTheta2);
00186     double rh            = ma * std::sqrt(1 - cosTheta1*cosTheta1);
00187     double ph            = phi();
00188     set(rh*std::cos(ph), rh*std::sin(ph), ma*cosTheta1);
00189   }
00190 
00191   //--------------------------------------------------------------------------
00192   template<>
00193   double BasicVector3D<double>::angle(const BasicVector3D<double> & v) const {
00194     double cosa = 0;
00195     double ptot = mag()*v.mag();
00196     if(ptot > 0) {
00197       cosa = dot(v)/ptot;
00198       if(cosa >  1) cosa =  1;
00199       if(cosa < -1) cosa = -1;
00200     }
00201     return std::acos(cosa);
00202   }
00203 
00204   //--------------------------------------------------------------------------
00205   template<>
00206   BasicVector3D<double> & BasicVector3D<double>::rotateX(double a) {
00207     double sina = std::sin(a), cosa = std::cos(a), dy = y(), dz = z();
00208     setY(dy*cosa-dz*sina);
00209     setZ(dz*cosa+dy*sina);
00210     return *this;
00211   }
00212 
00213   //--------------------------------------------------------------------------
00214   template<>
00215   BasicVector3D<double> & BasicVector3D<double>::rotateY(double a) {
00216     double sina = std::sin(a), cosa = std::cos(a), dz = z(), dx = x();
00217     setZ(dz*cosa-dx*sina);
00218     setX(dx*cosa+dz*sina);
00219     return *this;
00220   }
00221 
00222   //--------------------------------------------------------------------------
00223   template<>
00224   BasicVector3D<double> & BasicVector3D<double>::rotateZ(double a) {
00225     double sina = std::sin(a), cosa = std::cos(a), dx = x(), dy = y();
00226     setX(dx*cosa-dy*sina);
00227     setY(dy*cosa+dx*sina);
00228     return *this;
00229   }
00230 
00231   //--------------------------------------------------------------------------
00232   template<>
00233   BasicVector3D<double> &
00234   BasicVector3D<double>::rotate(double a, const BasicVector3D<double> & v) {
00235     if (a  == 0) return *this;
00236     double cx = v.x(), cy = v.y(), cz = v.z();
00237     double ll = std::sqrt(cx*cx + cy*cy + cz*cz);
00238     if (ll == 0) {
00239       std::cerr << "BasicVector<double>::rotate() : zero axis" << std::endl;
00240       return *this;
00241     }
00242     double cosa = std::cos(a), sina = std::sin(a);
00243     cx /= ll; cy /= ll; cz /= ll;   
00244 
00245     double xx = cosa + (1-cosa)*cx*cx;
00246     double xy =        (1-cosa)*cx*cy - sina*cz;
00247     double xz =        (1-cosa)*cx*cz + sina*cy;
00248     
00249     double yx =        (1-cosa)*cy*cx + sina*cz;
00250     double yy = cosa + (1-cosa)*cy*cy; 
00251     double yz =        (1-cosa)*cy*cz - sina*cx;
00252     
00253     double zx =        (1-cosa)*cz*cx - sina*cy;
00254     double zy =        (1-cosa)*cz*cy + sina*cx;
00255     double zz = cosa + (1-cosa)*cz*cz;
00256 
00257     cx = x(); cy = y(); cz = z();   
00258     set(xx*cx+xy*cy+xz*cz, yx*cx+yy*cy+yz*cz, zx*cx+zy*cy+zz*cz);
00259     return *this;
00260   }
00261 
00262   //--------------------------------------------------------------------------
00263   std::ostream &
00264   operator<< (std::ostream & os, const BasicVector3D<double> & a)
00265   {
00266     return os << "(" << a.x() << "," << a.y() << "," << a.z() << ")";
00267   }
00268   
00269   //--------------------------------------------------------------------------
00270   std::istream &
00271   operator>> (std::istream & is, BasicVector3D<double> & a)
00272   {
00273     // Required format is ( a, b, c ) that is, three numbers, preceded by
00274     // (, followed by ), and separated by commas.  The three numbers are
00275     // taken as x, y, z.
00276     
00277     double x, y, z;
00278     char c;
00279     
00280     is >> std::ws >> c;
00281     // ws is defined to invoke eatwhite(istream & )
00282     // see (Stroustrup gray book) page 333 and 345.
00283     if (is.fail() || c != '(' ) {
00284       std::cerr
00285         << "Could not find required opening parenthesis "
00286         << "in input of a BasicVector3D<double>"
00287         << std::endl;
00288       return is;
00289     }
00290 
00291     is >> x >> std::ws >> c;
00292     if (is.fail() || c != ',' ) {
00293       std::cerr
00294         << "Could not find x value and required trailing comma "
00295         << "in input of a BasicVector3D<double>"
00296         << std::endl; 
00297       return is;
00298     }
00299 
00300     is >> y >> std::ws >> c;
00301     if (is.fail() || c != ',' ) {
00302       std::cerr
00303         << "Could not find y value and required trailing comma "
00304         <<  "in input of a BasicVector3D<double>"
00305         << std::endl;
00306       return is;
00307     }
00308 
00309     is >> z >> std::ws >> c;
00310     if (is.fail() || c != ')' ) {
00311       std::cerr
00312         << "Could not find z value and required close parenthesis "
00313         << "in input of a BasicVector3D<double>"
00314         << std::endl;
00315       return is;
00316     }
00317 
00318     a.setX(x);
00319     a.setY(y);
00320     a.setZ(z);
00321     return is;
00322   }
00323 } /* namespace HepGeom */