TEvtGen/EvtGenBase/EvtVector4R.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtVector4R.cc
  12 //
  13 // Description: Real implementation of 4-vectors
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD  September 25, 1996            Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <math.h>
  24 #include <assert.h>
  25 #include "EvtGenBase/EvtVector4R.hh"
  26 #include "EvtGenBase/EvtVector3R.hh"
  27 #include "EvtGenBase/EvtVector4C.hh"
  28 #include "EvtGenBase/EvtTensor4C.hh"
  29
  30 using std::ostream;
  31
  32
  33
  34 EvtVector4R::EvtVector4R(double e,double p1,double p2, double p3){
  35
  36   v[0]=e; v[1]=p1; v[2]=p2; v[3]=p3;
  37 }
  38
  39 double EvtVector4R::mass() const{
  40
  41   double m2=v[0]*v[0]-v[1]*v[1]-v[2]*v[2]-v[3]*v[3];
  42
  43   if (m2>0.0) {
  44     return sqrt(m2);
  45   }
  46   else{
  47     return 0.0;
  48   }
  49 }
  50
  51
  52 EvtVector4R rotateEuler(const EvtVector4R& rs,
  53                         double alpha,double beta,double gamma){
  54
  55   EvtVector4R tmp(rs);
  56   tmp.applyRotateEuler(alpha,beta,gamma);
  57   return tmp;
  58
  59 }
  60
  61 EvtVector4R boostTo(const EvtVector4R& rs,
  62                     const EvtVector4R& p4){
  63
  64   EvtVector4R tmp(rs);
  65   tmp.applyBoostTo(p4);
  66   return tmp;
  67
  68 }
  69
  70 EvtVector4R boostTo(const EvtVector4R& rs,
  71                     const EvtVector3R& boost){
  72
  73   EvtVector4R tmp(rs);
  74   tmp.applyBoostTo(boost);
  75   return tmp;
  76
  77 }
  78
  79
  80
  81 void EvtVector4R::applyRotateEuler(double phi,double theta,double ksi){
  82
  83   double sp=sin(phi);
  84   double st=sin(theta);
  85   double sk=sin(ksi);
  86   double cp=cos(phi);
  87   double ct=cos(theta);
  88   double ck=cos(ksi);
  89
  90   double x=( ck*ct*cp-sk*sp)*v[1]+( -sk*ct*cp-ck*sp)*v[2]+st*cp*v[3];
  91   double y=( ck*ct*sp+sk*cp)*v[1]+(-sk*ct*sp+ck*cp)*v[2]+st*sp*v[3];
  92   double z=-ck*st*v[1]+sk*st*v[2]+ct*v[3];
  93
  94   v[1]=x;
  95   v[2]=y;
  96   v[3]=z;
  97
  98 }
  99
 100 ostream& operator<<(ostream& s, const EvtVector4R& v){
 101
 102   s<<"("<<v.v[0]<<","<<v.v[1]<<","<<v.v[2]<<","<<v.v[3]<<")";
 103
 104   return s;
 105
 106 }
 107
 108 void EvtVector4R::applyBoostTo(const EvtVector4R& p4){
 109
 110   double e=p4.get(0);
 111
 112   EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
 113
 114   applyBoostTo(boost);
 115
 116   return;
 117
 118 }
 119
 120 void EvtVector4R::applyBoostTo(const EvtVector3R& boost){
 121
 122   double bx,by,bz,gamma,b2;
 123
 124   bx=boost.get(0);
 125   by=boost.get(1);
 126   bz=boost.get(2);
 127
 128   double bxx=bx*bx;
 129   double byy=by*by;
 130   double bzz=bz*bz;
 131
 132   b2=bxx+byy+bzz;
 133
 134
 135   if (b2==0.0){
 136     return;
 137   }
 138
 139   assert(b2<1.0);
 140
 141   gamma=1.0/sqrt(1-b2);
 142
 143
 144   double gb2=(gamma-1.0)/b2;
 145
 146   double gb2xy=gb2*bx*by;
 147   double gb2xz=gb2*bx*bz;
 148   double gb2yz=gb2*by*bz;
 149
 150   double gbx=gamma*bx;
 151   double gby=gamma*by;
 152   double gbz=gamma*bz;
 153
 154   double e2=v[0];
 155   double px2=v[1];
 156   double py2=v[2];
 157   double pz2=v[3];
 158
 159   v[0]=gamma*e2+gbx*px2+gby*py2+gbz*pz2;
 160
 161   v[1]=gbx*e2+gb2*bxx*px2+px2+gb2xy*py2+gb2xz*pz2;
 162
 163   v[2]=gby*e2+gb2*byy*py2+py2+gb2xy*px2+gb2yz*pz2;
 164
 165   v[3]=gbz*e2+gb2*bzz*pz2+pz2+gb2yz*py2+gb2xz*px2;
 166
 167   return;
 168
 169 }
 170
 171 EvtVector4R EvtVector4R::cross( const EvtVector4R& p2 ){
 172
 173   //Calcs the cross product.  Added by djl on July 27, 1995.
 174   //Modified for real vectros by ryd Aug 28-96
 175
 176   EvtVector4R temp;
 177
 178   temp.v[0] = 0.0;
 179   temp.v[1] = v[2]*p2.v[3] - v[3]*p2.v[2];
 180   temp.v[2] = v[3]*p2.v[1] - v[1]*p2.v[3];
 181   temp.v[3] = v[1]*p2.v[2] - v[2]*p2.v[1];
 182
 183   return temp;
 184 }
 185
 186 double EvtVector4R::d3mag() const
 187
 188 // returns the 3 momentum mag.
 189 {
 190   double temp;
 191
 192   temp = v[1]*v[1]+v[2]*v[2]+v[3]*v[3];
 193
 194   temp = sqrt( temp );
 195
 196   return temp;
 197 } // r3mag
 198
 199 double EvtVector4R::dot ( const EvtVector4R& p2 )const{
 200
 201   //Returns the dot product of the 3 momentum.  Added by
 202   //djl on July 27, 1995.  for real!!!
 203
 204   double temp;
 205
 206   temp = v[1]*p2.v[1];
 207   temp += v[2]*p2.v[2];
 208   temp += v[3]*p2.v[3];
 209
 210   return temp;
 211
 212 } //dot
 213
 214 // Functions below added by AJB
 215
 216 // Calculate ( \vec{p1} cross \vec{p2} ) \cdot \vec{p3} in rest frame of object
 217 double EvtVector4R::scalartripler3( const EvtVector4R& p1,
 218         const EvtVector4R& p2, const EvtVector4R& p3 ) const
 219 {
 220     EvtVector4C lc=dual(directProd(*this, p1)).cont2(p2);
 221     EvtVector4R  l(real(lc.get(0)), real(lc.get(1)), real(lc.get(2)),
 222             real(lc.get(3)));
 223
 224     return -1.0/mass() * (l * p3);
 225 }
 226
 227 // Calculate the 3-d dot product of 4-vectors p1 and p2 in the rest frame of
 228 // 4-vector p0
 229 double EvtVector4R::dotr3( const EvtVector4R& p1, const EvtVector4R& p2 ) const
 230 {
 231     return 1/mass2() * ((*this) * p1) * ((*this) * p2) - p1 * p2;
 232 }
 233
 234 // Calculate the 3-d magnitude squared of 4-vector p1 in the rest frame of
 235 // 4-vector p0
 236 double EvtVector4R::mag2r3( const EvtVector4R& p1 ) const
 237 {
 238     return Square((*this) * p1)/mass2() - p1.mass2();
 239 }
 240
 241 // Calculate the 3-d magnitude 4-vector p1 in the rest frame of 4-vector p0.
 242 double EvtVector4R::magr3( const EvtVector4R& p1 ) const
 243 {
 244     return sqrt(mag2r3(p1));
 245 }
 246
 247
 248
 249
 250
 251
 252