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/EvtComplex.hh"
00026 #include "EvtGenBase/EvtVector4C.hh"
00027 #include "EvtGenBase/EvtTensor4C.hh"
00028 using std::endl;
00029 using std::ostream;
00030
00031
00032
00033 EvtTensor4C::EvtTensor4C( const EvtTensor4C& t1 ) {
00034
00035 int i,j;
00036
00037 for(i=0;i<4;i++) {
00038 for(j=0;j<4;j++) {
00039 t[i][j] = t1.t[i][j];
00040 }
00041 }
00042
00043 }
00044
00045 EvtTensor4C::~EvtTensor4C() { }
00046
00047 const EvtTensor4C& EvtTensor4C::g(){
00048
00049 static EvtTensor4C g_metric(1.0,-1.0,-1.0,-1.0);
00050
00051 return g_metric;
00052
00053 }
00054
00055 EvtTensor4C& EvtTensor4C::operator=(const EvtTensor4C& t1) {
00056 int i,j;
00057
00058 for(i=0;i<4;i++) {
00059 for(j=0;j<4;j++) {
00060 t[i][j] = t1.t[i][j];
00061 }
00062 }
00063 return *this;
00064 }
00065
00066 EvtTensor4C EvtTensor4C::conj() const {
00067 EvtTensor4C temp;
00068
00069 int i,j;
00070
00071 for(i=0;i<4;i++) {
00072 for(j=0;j<4;j++) {
00073 temp.set(j,i,::conj(t[i][j]));
00074 }
00075 }
00076 return temp;
00077 }
00078
00079
00080 EvtTensor4C rotateEuler(const EvtTensor4C& rs,
00081 double alpha,double beta,double gamma){
00082
00083 EvtTensor4C tmp(rs);
00084 tmp.applyRotateEuler(alpha,beta,gamma);
00085 return tmp;
00086
00087 }
00088
00089 EvtTensor4C boostTo(const EvtTensor4C& rs,
00090 const EvtVector4R p4){
00091
00092 EvtTensor4C tmp(rs);
00093 tmp.applyBoostTo(p4);
00094 return tmp;
00095
00096 }
00097
00098 EvtTensor4C boostTo(const EvtTensor4C& rs,
00099 const EvtVector3R boost){
00100
00101 EvtTensor4C tmp(rs);
00102 tmp.applyBoostTo(boost);
00103 return tmp;
00104
00105 }
00106
00107 void EvtTensor4C::applyBoostTo(const EvtVector4R& p4){
00108
00109 double e=p4.get(0);
00110
00111 EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
00112
00113 applyBoostTo(boost);
00114
00115 return;
00116
00117 }
00118
00119
00120 void EvtTensor4C::applyBoostTo(const EvtVector3R& boost){
00121
00122 double bx,by,bz,gamma,b2;
00123 double lambda[4][4];
00124 EvtComplex tt[4][4];
00125
00126 bx=boost.get(0);
00127 by=boost.get(1);
00128 bz=boost.get(2);
00129
00130 double bxx=bx*bx;
00131 double byy=by*by;
00132 double bzz=bz*bz;
00133
00134 b2=bxx+byy+bzz;
00135
00136
00137 if (b2==0.0){
00138 return;
00139 }
00140
00141 assert(b2<1.0);
00142
00143 gamma=1.0/sqrt(1-b2);
00144
00145
00146 int i,j,k;
00147
00148
00149 if (b2==0.0){
00150 return ;
00151 }
00152
00153 lambda[0][0]=gamma;
00154 lambda[0][1]=gamma*bx;
00155 lambda[1][0]=gamma*bx;
00156 lambda[0][2]=gamma*by;
00157 lambda[2][0]=gamma*by;
00158 lambda[0][3]=gamma*bz;
00159 lambda[3][0]=gamma*bz;
00160
00161 lambda[1][1]=1.0+(gamma-1.0)*bx*bx/b2;
00162 lambda[2][2]=1.0+(gamma-1.0)*by*by/b2;
00163 lambda[3][3]=1.0+(gamma-1.0)*bz*bz/b2;
00164
00165 lambda[1][2]=(gamma-1.0)*bx*by/b2;
00166 lambda[2][1]=(gamma-1.0)*bx*by/b2;
00167
00168 lambda[1][3]=(gamma-1.0)*bx*bz/b2;
00169 lambda[3][1]=(gamma-1.0)*bx*bz/b2;
00170
00171 lambda[3][2]=(gamma-1.0)*bz*by/b2;
00172 lambda[2][3]=(gamma-1.0)*bz*by/b2;
00173
00174 for(i=0;i<4;i++){
00175 for(j=0;j<4;j++){
00176 tt[i][j] = EvtComplex(0.0);
00177 for(k=0;k<4;k++){
00178 tt[i][j]=tt[i][j]+lambda[j][k]*t[i][k];
00179 }
00180 }
00181 }
00182
00183 for(i=0;i<4;i++){
00184 for(j=0;j<4;j++){
00185 t[i][j] = EvtComplex(0.0);
00186 for(k=0;k<4;k++){
00187 t[i][j]=t[i][j]+lambda[i][k]*tt[k][j];
00188 }
00189 }
00190 }
00191
00192 }
00193
00194 void EvtTensor4C::zero(){
00195 int i,j;
00196 for(i=0;i<4;i++){
00197 for(j=0;j<4;j++){
00198 t[i][j]=EvtComplex(0.0,0.0);
00199 }
00200 }
00201 }
00202
00203
00204
00205 ostream& operator<<(ostream& s,const EvtTensor4C& t){
00206
00207 int i,j;
00208 s<< endl;
00209 for(i=0;i<4;i++){
00210 for(j=0;j<4;j++){
00211 s << t.t[i][j];
00212 }
00213 s << endl;
00214 }
00215 return s;
00216 }
00217
00218 void EvtTensor4C::setdiag(double g00, double g11, double g22, double g33){
00219 t[0][0]=EvtComplex(g00);
00220 t[1][1]=EvtComplex(g11);
00221 t[2][2]=EvtComplex(g22);
00222 t[3][3]=EvtComplex(g33);
00223 t[0][1] = EvtComplex(0.0);
00224 t[0][2] = EvtComplex(0.0);
00225 t[0][3] = EvtComplex(0.0);
00226 t[1][0] = EvtComplex(0.0);
00227 t[1][2] = EvtComplex(0.0);
00228 t[1][3] = EvtComplex(0.0);
00229 t[2][0] = EvtComplex(0.0);
00230 t[2][1] = EvtComplex(0.0);
00231 t[2][3] = EvtComplex(0.0);
00232 t[3][0] = EvtComplex(0.0);
00233 t[3][1] = EvtComplex(0.0);
00234 t[3][2] = EvtComplex(0.0);
00235 }
00236
00237
00238 EvtTensor4C& EvtTensor4C::operator+=(const EvtTensor4C& t2){
00239
00240 int i,j;
00241
00242 for (i=0;i<4;i++) {
00243 for (j=0;j<4;j++) {
00244 t[i][j]+=t2.get(i,j);
00245 }
00246 }
00247 return *this;
00248 }
00249
00250 EvtTensor4C& EvtTensor4C::operator-=(const EvtTensor4C& t2){
00251
00252 int i,j;
00253
00254 for (i=0;i<4;i++) {
00255 for (j=0;j<4;j++) {
00256 t[i][j]-=t2.get(i,j);
00257 }
00258 }
00259 return *this;
00260 }
00261
00262
00263 EvtTensor4C& EvtTensor4C::operator*=(const EvtComplex& c) {
00264 int i,j;
00265
00266 for (i=0;i<4;i++) {
00267 for (j=0;j<4;j++) {
00268 t[i][j]*=c;
00269 }
00270 }
00271 return *this;
00272 }
00273
00274
00275 EvtTensor4C operator*(const EvtTensor4C& t1,const EvtComplex& c){
00276
00277 return EvtTensor4C(t1)*=c;
00278
00279 }
00280
00281 EvtTensor4C operator*(const EvtComplex& c,const EvtTensor4C& t1){
00282
00283 return EvtTensor4C(t1)*=c;
00284
00285 }
00286
00287
00288 EvtTensor4C& EvtTensor4C::operator*=(double d) {
00289 int i,j;
00290
00291 for (i=0;i<4;i++) {
00292 for (j=0;j<4;j++) {
00293 t[i][j]*=EvtComplex(d,0.0);
00294 }
00295 }
00296 return *this;
00297 }
00298
00299
00300 EvtTensor4C operator*(const EvtTensor4C& t1, double d){
00301
00302 return EvtTensor4C(t1)*=EvtComplex(d,0.0);
00303
00304 }
00305
00306 EvtTensor4C operator*(double d, const EvtTensor4C& t1){
00307
00308 return EvtTensor4C(t1)*=EvtComplex(d,0.0);
00309
00310 }
00311
00312 EvtComplex cont(const EvtTensor4C& t1,const EvtTensor4C& t2){
00313
00314 EvtComplex sum(0.0,0.0);
00315 int i,j;
00316
00317 for (i=0;i<4;i++) {
00318 for (j=0;j<4;j++) {
00319 sum+=t1.t[i][j]*t2.t[i][j];
00320 }
00321 }
00322
00323 return sum;
00324 }
00325
00326
00327 EvtTensor4C directProd(const EvtVector4C& c1,const EvtVector4C& c2){
00328 EvtTensor4C temp;
00329 int i,j;
00330
00331 for (i=0;i<4;i++) {
00332 for (j=0;j<4;j++) {
00333 temp.set(i,j,c1.get(i)*c2.get(j));
00334 }
00335 }
00336 return temp;
00337 }
00338
00339
00340 EvtTensor4C directProd(const EvtVector4C& c1,const EvtVector4R& c2){
00341 EvtTensor4C temp;
00342 int i,j;
00343
00344 for (i=0;i<4;i++) {
00345 for (j=0;j<4;j++) {
00346 temp.set(i,j,c1.get(i)*c2.get(j));
00347 }
00348 }
00349 return temp;
00350 }
00351
00352
00353 EvtTensor4C directProd(const EvtVector4R& c1,const EvtVector4R& c2){
00354
00355 EvtTensor4C temp;
00356 int i,j;
00357
00358 for (i=0;i<4;i++) {
00359 for (j=0;j<4;j++) {
00360 temp.t[i][j]=EvtComplex(c1.get(i)*c2.get(j),0.0);
00361 }
00362 }
00363 return temp;
00364 }
00365
00366 EvtTensor4C& EvtTensor4C::addDirProd(const EvtVector4R& p1,const EvtVector4R& p2){
00367
00368 int i,j;
00369
00370 for (i=0;i<4;i++) {
00371 for (j=0;j<4;j++) {
00372 t[i][j]+=p1.get(i)*p2.get(j);
00373 }
00374 }
00375 return *this;
00376 }
00377
00378
00379 EvtTensor4C dual(const EvtTensor4C& t2){
00380
00381 EvtTensor4C temp;
00382
00383 temp.set(0,0,EvtComplex(0.0,0.0));
00384 temp.set(1,1,EvtComplex(0.0,0.0));
00385 temp.set(2,2,EvtComplex(0.0,0.0));
00386 temp.set(3,3,EvtComplex(0.0,0.0));
00387
00388 temp.set(0,1,t2.get(3,2)-t2.get(2,3));
00389 temp.set(0,2,-t2.get(3,1)+t2.get(1,3));
00390 temp.set(0,3,t2.get(2,1)-t2.get(1,2));
00391
00392 temp.set(1,2,-t2.get(3,0)+t2.get(0,3));
00393 temp.set(1,3,t2.get(2,0)-t2.get(0,2));
00394
00395 temp.set(2,3,-t2.get(1,0)+t2.get(0,1));
00396
00397 temp.set(1,0,-temp.get(0,1));
00398 temp.set(2,0,-temp.get(0,2));
00399 temp.set(3,0,-temp.get(0,3));
00400
00401 temp.set(2,1,-temp.get(1,2));
00402 temp.set(3,1,-temp.get(1,3));
00403
00404 temp.set(3,2,-temp.get(2,3));
00405
00406 return temp;
00407
00408 }
00409
00410
00411 EvtTensor4C conj(const EvtTensor4C& t2) {
00412 EvtTensor4C temp;
00413
00414 int i,j;
00415
00416 for(i=0;i<4;i++){
00417 for(j=0;j<4;j++){
00418 temp.set(i,j,::conj((t2.get(i,j))));
00419 }
00420 }
00421
00422 return temp;
00423 }
00424
00425
00426 EvtTensor4C cont22(const EvtTensor4C& t1,const EvtTensor4C& t2){
00427 EvtTensor4C temp;
00428
00429 int i,j;
00430 EvtComplex c;
00431
00432 for(i=0;i<4;i++){
00433 for(j=0;j<4;j++){
00434 c=t1.get(i,0)*t2.get(j,0)-t1.get(i,1)*t2.get(j,1)
00435 -t1.get(i,2)*t2.get(j,2)-t1.get(i,3)*t2.get(j,3);
00436 temp.set(i,j,c);
00437 }
00438 }
00439
00440 return temp;
00441 }
00442
00443 EvtTensor4C cont11(const EvtTensor4C& t1,const EvtTensor4C& t2){
00444 EvtTensor4C temp;
00445
00446 int i,j;
00447 EvtComplex c;
00448
00449 for(i=0;i<4;i++){
00450 for(j=0;j<4;j++){
00451 c=t1.get(0,i)*t2.get(0,j)-t1.get(1,i)*t2.get(1,j)
00452 -t1.get(2,i)*t2.get(2,j)-t1.get(3,i)*t2.get(3,j);
00453 temp.set(i,j,c);
00454 }
00455 }
00456
00457 return temp;
00458 }
00459
00460
00461 EvtVector4C EvtTensor4C::cont1(const EvtVector4C& v4) const {
00462 EvtVector4C temp;
00463
00464 int i;
00465
00466 for(i=0;i<4;i++){
00467 temp.set(i,t[0][i]*v4.get(0)-t[1][i]*v4.get(1)
00468 -t[2][i]*v4.get(2)-t[3][i]*v4.get(3));
00469 }
00470
00471 return temp;
00472 }
00473
00474 EvtVector4C EvtTensor4C::cont2(const EvtVector4C& v4) const {
00475 EvtVector4C temp;
00476
00477 int i;
00478
00479 for(i=0;i<4;i++){
00480 temp.set(i,t[i][0]*v4.get(0)-t[i][1]*v4.get(1)
00481 -t[i][2]*v4.get(2)-t[i][3]*v4.get(3));
00482 }
00483
00484 return temp;
00485 }
00486
00487
00488 EvtVector4C EvtTensor4C::cont1(const EvtVector4R& v4) const {
00489 EvtVector4C temp;
00490
00491 int i;
00492
00493 for(i=0;i<4;i++){
00494 temp.set(i,t[0][i]*v4.get(0)-t[1][i]*v4.get(1)
00495 -t[2][i]*v4.get(2)-t[3][i]*v4.get(3));
00496 }
00497
00498 return temp;
00499 }
00500
00501
00502 EvtVector4C EvtTensor4C::cont2(const EvtVector4R& v4) const {
00503 EvtVector4C temp;
00504
00505 int i;
00506
00507 for(i=0;i<4;i++){
00508 temp.set(i,t[i][0]*v4.get(0)-t[i][1]*v4.get(1)
00509 -t[i][2]*v4.get(2)-t[i][3]*v4.get(3));
00510 }
00511
00512 return temp;
00513 }
00514
00515
00516
00517 void EvtTensor4C::applyRotateEuler(double phi,double theta,double ksi){
00518
00519 EvtComplex tt[4][4];
00520 double sp,st,sk,cp,ct,ck;
00521 double lambda[4][4];
00522
00523 sp=sin(phi);
00524 st=sin(theta);
00525 sk=sin(ksi);
00526 cp=cos(phi);
00527 ct=cos(theta);
00528 ck=cos(ksi);
00529
00530
00531 lambda[0][0]=1.0;
00532 lambda[0][1]=0.0;
00533 lambda[1][0]=0.0;
00534 lambda[0][2]=0.0;
00535 lambda[2][0]=0.0;
00536 lambda[0][3]=0.0;
00537 lambda[3][0]=0.0;
00538
00539 lambda[1][1]= ck*ct*cp-sk*sp;
00540 lambda[1][2]=-sk*ct*cp-ck*sp;
00541 lambda[1][3]=st*cp;
00542
00543 lambda[2][1]= ck*ct*sp+sk*cp;
00544 lambda[2][2]=-sk*ct*sp+ck*cp;
00545 lambda[2][3]=st*sp;
00546
00547 lambda[3][1]=-ck*st;
00548 lambda[3][2]=sk*st;
00549 lambda[3][3]=ct;
00550
00551
00552 int i,j,k;
00553
00554
00555 for(i=0;i<4;i++){
00556 for(j=0;j<4;j++){
00557 tt[i][j] = EvtComplex(0.0);
00558 for(k=0;k<4;k++){
00559 tt[i][j]+=lambda[j][k]*t[i][k];
00560 }
00561 }
00562 }
00563
00564 for(i=0;i<4;i++){
00565 for(j=0;j<4;j++){
00566 t[i][j] = EvtComplex(0.0);
00567 for(k=0;k<4;k++){
00568 t[i][j]+=lambda[i][k]*tt[k][j];
00569 }
00570 }
00571 }
00572
00573 }
00574
00575
00576