00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <iostream>
00023 #include <math.h>
00024 #include <assert.h>
00025 #include "EvtGenBase/EvtVector4R.hh"
00026 #include "EvtGenBase/EvtVector3R.hh"
00027 #include "EvtGenBase/EvtVector4C.hh"
00028 #include "EvtGenBase/EvtTensor4C.hh"
00029
00030 using std::ostream;
00031
00032
00033
00034 EvtVector4R::EvtVector4R(double e,double p1,double p2, double p3){
00035
00036 v[0]=e; v[1]=p1; v[2]=p2; v[3]=p3;
00037 }
00038
00039 double EvtVector4R::mass() const{
00040
00041 double m2=v[0]*v[0]-v[1]*v[1]-v[2]*v[2]-v[3]*v[3];
00042
00043 if (m2>0.0) {
00044 return sqrt(m2);
00045 }
00046 else{
00047 return 0.0;
00048 }
00049 }
00050
00051
00052 EvtVector4R rotateEuler(const EvtVector4R& rs,
00053 double alpha,double beta,double gamma){
00054
00055 EvtVector4R tmp(rs);
00056 tmp.applyRotateEuler(alpha,beta,gamma);
00057 return tmp;
00058
00059 }
00060
00061 EvtVector4R boostTo(const EvtVector4R& rs,
00062 const EvtVector4R& p4){
00063
00064 EvtVector4R tmp(rs);
00065 tmp.applyBoostTo(p4);
00066 return tmp;
00067
00068 }
00069
00070 EvtVector4R boostTo(const EvtVector4R& rs,
00071 const EvtVector3R& boost){
00072
00073 EvtVector4R tmp(rs);
00074 tmp.applyBoostTo(boost);
00075 return tmp;
00076
00077 }
00078
00079
00080
00081 void EvtVector4R::applyRotateEuler(double phi,double theta,double ksi){
00082
00083 double sp=sin(phi);
00084 double st=sin(theta);
00085 double sk=sin(ksi);
00086 double cp=cos(phi);
00087 double ct=cos(theta);
00088 double ck=cos(ksi);
00089
00090 double x=( ck*ct*cp-sk*sp)*v[1]+( -sk*ct*cp-ck*sp)*v[2]+st*cp*v[3];
00091 double y=( ck*ct*sp+sk*cp)*v[1]+(-sk*ct*sp+ck*cp)*v[2]+st*sp*v[3];
00092 double z=-ck*st*v[1]+sk*st*v[2]+ct*v[3];
00093
00094 v[1]=x;
00095 v[2]=y;
00096 v[3]=z;
00097
00098 }
00099
00100 ostream& operator<<(ostream& s, const EvtVector4R& v){
00101
00102 s<<"("<<v.v[0]<<","<<v.v[1]<<","<<v.v[2]<<","<<v.v[3]<<")";
00103
00104 return s;
00105
00106 }
00107
00108 void EvtVector4R::applyBoostTo(const EvtVector4R& p4){
00109
00110 double e=p4.get(0);
00111
00112 EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
00113
00114 applyBoostTo(boost);
00115
00116 return;
00117
00118 }
00119
00120 void EvtVector4R::applyBoostTo(const EvtVector3R& boost){
00121
00122 double bx,by,bz,gamma,b2;
00123
00124 bx=boost.get(0);
00125 by=boost.get(1);
00126 bz=boost.get(2);
00127
00128 double bxx=bx*bx;
00129 double byy=by*by;
00130 double bzz=bz*bz;
00131
00132 b2=bxx+byy+bzz;
00133
00134
00135 if (b2==0.0){
00136 return;
00137 }
00138
00139 assert(b2<1.0);
00140
00141 gamma=1.0/sqrt(1-b2);
00142
00143
00144 double gb2=(gamma-1.0)/b2;
00145
00146 double gb2xy=gb2*bx*by;
00147 double gb2xz=gb2*bx*bz;
00148 double gb2yz=gb2*by*bz;
00149
00150 double gbx=gamma*bx;
00151 double gby=gamma*by;
00152 double gbz=gamma*bz;
00153
00154 double e2=v[0];
00155 double px2=v[1];
00156 double py2=v[2];
00157 double pz2=v[3];
00158
00159 v[0]=gamma*e2+gbx*px2+gby*py2+gbz*pz2;
00160
00161 v[1]=gbx*e2+gb2*bxx*px2+px2+gb2xy*py2+gb2xz*pz2;
00162
00163 v[2]=gby*e2+gb2*byy*py2+py2+gb2xy*px2+gb2yz*pz2;
00164
00165 v[3]=gbz*e2+gb2*bzz*pz2+pz2+gb2yz*py2+gb2xz*px2;
00166
00167 return;
00168
00169 }
00170
00171 EvtVector4R EvtVector4R::cross( const EvtVector4R& p2 ){
00172
00173
00174
00175
00176 EvtVector4R temp;
00177
00178 temp.v[0] = 0.0;
00179 temp.v[1] = v[2]*p2.v[3] - v[3]*p2.v[2];
00180 temp.v[2] = v[3]*p2.v[1] - v[1]*p2.v[3];
00181 temp.v[3] = v[1]*p2.v[2] - v[2]*p2.v[1];
00182
00183 return temp;
00184 }
00185
00186 double EvtVector4R::d3mag() const
00187
00188
00189 {
00190 double temp;
00191
00192 temp = v[1]*v[1]+v[2]*v[2]+v[3]*v[3];
00193
00194 temp = sqrt( temp );
00195
00196 return temp;
00197 }
00198
00199 double EvtVector4R::dot ( const EvtVector4R& p2 )const{
00200
00201
00202
00203
00204 double temp;
00205
00206 temp = v[1]*p2.v[1];
00207 temp += v[2]*p2.v[2];
00208 temp += v[3]*p2.v[3];
00209
00210 return temp;
00211
00212 }
00213
00214
00215
00216
00217 double EvtVector4R::scalartripler3( const EvtVector4R& p1,
00218 const EvtVector4R& p2, const EvtVector4R& p3 ) const
00219 {
00220 EvtVector4C lc=dual(directProd(*this, p1)).cont2(p2);
00221 EvtVector4R l(real(lc.get(0)), real(lc.get(1)), real(lc.get(2)),
00222 real(lc.get(3)));
00223
00224 return -1.0/mass() * (l * p3);
00225 }
00226
00227
00228
00229 double EvtVector4R::dotr3( const EvtVector4R& p1, const EvtVector4R& p2 ) const
00230 {
00231 return 1/mass2() * ((*this) * p1) * ((*this) * p2) - p1 * p2;
00232 }
00233
00234
00235
00236 double EvtVector4R::mag2r3( const EvtVector4R& p1 ) const
00237 {
00238 return Square((*this) * p1)/mass2() - p1.mass2();
00239 }
00240
00241
00242 double EvtVector4R::magr3( const EvtVector4R& p1 ) const
00243 {
00244 return sqrt(mag2r3(p1));
00245 }
00246
00247
00248
00249 double EvtVector4R::theta(){
00250 return (v[1]==0 && v[2]==0 && v[3]==0 )?0: atan2(sqrt(v[1]*v[1]+v[2]*v[2]),v[3]);
00251 }
00252
00253
00254 double EvtVector4R::phi(){
00255 return (v[1]==0 && v[2]==0)?0: atan2(v[2],v[1]);
00256 }
00257
00258
00259
00260
00261