TEvtGen/EvtGenBase/EvtTensor4C.cpp

   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: EvtTensor4C.cc
  12 //
  13 // Description: Implementation of tensor particles.
  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/EvtComplex.hh"
  26 #include "EvtGenBase/EvtVector4C.hh"
  27 #include "EvtGenBase/EvtTensor4C.hh"
  28 using std::endl;
  29 using std::ostream;
  30
  31
  32
  33 EvtTensor4C::EvtTensor4C( const EvtTensor4C& t1 ) {
  34
  35   int i,j;
  36
  37   for(i=0;i<4;i++) {
  38     for(j=0;j<4;j++) {
  39       t[i][j] = t1.t[i][j];
  40     }
  41   }
  42
  43 }
  44
  45 EvtTensor4C::~EvtTensor4C() { }
  46
  47 const EvtTensor4C& EvtTensor4C::g(){
  48
  49   static EvtTensor4C g_metric(1.0,-1.0,-1.0,-1.0);
  50
  51   return g_metric;
  52
  53 }
  54
  55 EvtTensor4C& EvtTensor4C::operator=(const EvtTensor4C& t1) {
  56   int i,j;
  57
  58   for(i=0;i<4;i++) {
  59     for(j=0;j<4;j++) {
  60       t[i][j] = t1.t[i][j];
  61     }
  62   }
  63   return *this;
  64 }
  65
  66 EvtTensor4C EvtTensor4C::conj() const {
  67   EvtTensor4C temp;
  68
  69   int i,j;
  70
  71   for(i=0;i<4;i++) {
  72     for(j=0;j<4;j++) {
  73       temp.set(j,i,::conj(t[i][j]));
  74     }
  75   }
  76   return temp;
  77 }
  78
  79
  80 EvtTensor4C rotateEuler(const EvtTensor4C& rs,
  81                         double alpha,double beta,double gamma){
  82
  83   EvtTensor4C tmp(rs);
  84   tmp.applyRotateEuler(alpha,beta,gamma);
  85   return tmp;
  86
  87 }
  88
  89 EvtTensor4C boostTo(const EvtTensor4C& rs,
  90                     const EvtVector4R p4){
  91
  92   EvtTensor4C tmp(rs);
  93   tmp.applyBoostTo(p4);
  94   return tmp;
  95
  96 }
  97
  98 EvtTensor4C boostTo(const EvtTensor4C& rs,
  99                     const EvtVector3R boost){
 100
 101   EvtTensor4C tmp(rs);
 102   tmp.applyBoostTo(boost);
 103   return tmp;
 104
 105 }
 106
 107 void EvtTensor4C::applyBoostTo(const EvtVector4R& p4){
 108
 109   double e=p4.get(0);
 110
 111   EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
 112
 113   applyBoostTo(boost);
 114
 115   return;
 116
 117 }
 118
 119
 120 void EvtTensor4C::applyBoostTo(const EvtVector3R& boost){
 121
 122   double bx,by,bz,gamma,b2;
 123   double lambda[4][4];
 124   EvtComplex tt[4][4];
 125
 126   bx=boost.get(0);
 127   by=boost.get(1);
 128   bz=boost.get(2);
 129
 130   double bxx=bx*bx;
 131   double byy=by*by;
 132   double bzz=bz*bz;
 133
 134   b2=bxx+byy+bzz;
 135
 136
 137   if (b2==0.0){
 138     return;
 139   }
 140
 141   assert(b2<1.0);
 142
 143   gamma=1.0/sqrt(1-b2);
 144
 145
 146   int i,j,k;
 147
 148
 149   if (b2==0.0){
 150     return ;
 151   }
 152
 153   lambda[0][0]=gamma;
 154   lambda[0][1]=gamma*bx;
 155   lambda[1][0]=gamma*bx;
 156   lambda[0][2]=gamma*by;
 157   lambda[2][0]=gamma*by;
 158   lambda[0][3]=gamma*bz;
 159   lambda[3][0]=gamma*bz;
 160
 161   lambda[1][1]=1.0+(gamma-1.0)*bx*bx/b2;
 162   lambda[2][2]=1.0+(gamma-1.0)*by*by/b2;
 163   lambda[3][3]=1.0+(gamma-1.0)*bz*bz/b2;
 164
 165   lambda[1][2]=(gamma-1.0)*bx*by/b2;
 166   lambda[2][1]=(gamma-1.0)*bx*by/b2;
 167
 168   lambda[1][3]=(gamma-1.0)*bx*bz/b2;
 169   lambda[3][1]=(gamma-1.0)*bx*bz/b2;
 170
 171   lambda[3][2]=(gamma-1.0)*bz*by/b2;
 172   lambda[2][3]=(gamma-1.0)*bz*by/b2;
 173
 174   for(i=0;i<4;i++){
 175     for(j=0;j<4;j++){
 176       tt[i][j] = EvtComplex(0.0);
 177       for(k=0;k<4;k++){
 178         tt[i][j]=tt[i][j]+lambda[j][k]*t[i][k];
 179       }
 180     }
 181   }
 182
 183   for(i=0;i<4;i++){
 184     for(j=0;j<4;j++){
 185       t[i][j] = EvtComplex(0.0);
 186       for(k=0;k<4;k++){
 187         t[i][j]=t[i][j]+lambda[i][k]*tt[k][j];
 188       }
 189     }
 190   }
 191
 192 }
 193
 194 void EvtTensor4C::zero(){
 195   int i,j;
 196   for(i=0;i<4;i++){
 197     for(j=0;j<4;j++){
 198       t[i][j]=EvtComplex(0.0,0.0);
 199     }
 200   }
 201 }
 202
 203
 204
 205 ostream& operator<<(ostream& s,const EvtTensor4C& t){
 206
 207   int i,j;
 208   s<< endl;
 209   for(i=0;i<4;i++){
 210     for(j=0;j<4;j++){
 211       s << t.t[i][j];
 212     }
 213    s << endl;
 214   }
 215   return s;
 216 }
 217
 218 void EvtTensor4C::setdiag(double g00, double g11, double g22, double g33){
 219   t[0][0]=EvtComplex(g00);
 220   t[1][1]=EvtComplex(g11);
 221   t[2][2]=EvtComplex(g22);
 222   t[3][3]=EvtComplex(g33);
 223   t[0][1] = EvtComplex(0.0);
 224   t[0][2] = EvtComplex(0.0);
 225   t[0][3] = EvtComplex(0.0);
 226   t[1][0] = EvtComplex(0.0);
 227   t[1][2] = EvtComplex(0.0);
 228   t[1][3] = EvtComplex(0.0);
 229   t[2][0] = EvtComplex(0.0);
 230   t[2][1] = EvtComplex(0.0);
 231   t[2][3] = EvtComplex(0.0);
 232   t[3][0] = EvtComplex(0.0);
 233   t[3][1] = EvtComplex(0.0);
 234   t[3][2] = EvtComplex(0.0);
 235 }
 236
 237
 238 EvtTensor4C& EvtTensor4C::operator+=(const EvtTensor4C& t2){
 239
 240   int i,j;
 241
 242   for (i=0;i<4;i++) {
 243     for (j=0;j<4;j++) {
 244       t[i][j]+=t2.get(i,j);
 245     }
 246   }
 247   return *this;
 248 }
 249
 250 EvtTensor4C& EvtTensor4C::operator-=(const EvtTensor4C& t2){
 251
 252   int i,j;
 253
 254   for (i=0;i<4;i++) {
 255     for (j=0;j<4;j++) {
 256       t[i][j]-=t2.get(i,j);
 257     }
 258   }
 259   return *this;
 260 }
 261
 262
 263 EvtTensor4C& EvtTensor4C::operator*=(const EvtComplex& c) {
 264   int i,j;
 265
 266   for (i=0;i<4;i++) {
 267     for (j=0;j<4;j++) {
 268       t[i][j]*=c;
 269     }
 270   }
 271   return *this;
 272 }
 273
 274
 275 EvtTensor4C operator*(const EvtTensor4C& t1,const EvtComplex& c){
 276
 277   return EvtTensor4C(t1)*=c;
 278
 279 }
 280
 281 EvtTensor4C operator*(const EvtComplex& c,const EvtTensor4C& t1){
 282
 283   return EvtTensor4C(t1)*=c;
 284
 285 }
 286
 287
 288 EvtTensor4C& EvtTensor4C::operator*=(double d) {
 289   int i,j;
 290
 291   for (i=0;i<4;i++) {
 292     for (j=0;j<4;j++) {
 293       t[i][j]*=EvtComplex(d,0.0);
 294     }
 295   }
 296   return *this;
 297 }
 298
 299
 300 EvtTensor4C operator*(const EvtTensor4C& t1, double d){
 301
 302   return EvtTensor4C(t1)*=EvtComplex(d,0.0);
 303
 304 }
 305
 306 EvtTensor4C operator*(double d, const EvtTensor4C& t1){
 307
 308   return EvtTensor4C(t1)*=EvtComplex(d,0.0);
 309
 310 }
 311
 312 EvtComplex cont(const EvtTensor4C& t1,const EvtTensor4C& t2){
 313
 314   EvtComplex sum(0.0,0.0);
 315   int i,j;
 316
 317   for (i=0;i<4;i++) {
 318     for (j=0;j<4;j++) {
 319       if ((i==0&&j!=0) || (j==0&&i!=0)) {
 320         sum -= t1.t[i][j]*t2.t[i][j];
 321       } else {
 322         sum += t1.t[i][j]*t2.t[i][j];
 323       }
 324     }
 325   }
 326
 327   return sum;
 328 }
 329
 330
 331 EvtTensor4C EvtGenFunctions::directProd(const EvtVector4C& c1,
 332                                         const EvtVector4C& c2){
 333   EvtTensor4C temp;
 334   int i,j;
 335
 336   for (i=0;i<4;i++) {
 337     for (j=0;j<4;j++) {
 338       temp.set(i,j,c1.get(i)*c2.get(j));
 339     }
 340   }
 341   return temp;
 342 }
 343
 344
 345 EvtTensor4C EvtGenFunctions::directProd(const EvtVector4C& c1,
 346                                         const EvtVector4R& c2){
 347   EvtTensor4C temp;
 348   int i,j;
 349
 350   for (i=0;i<4;i++) {
 351     for (j=0;j<4;j++) {
 352       temp.set(i,j,c1.get(i)*c2.get(j));
 353     }
 354   }
 355   return temp;
 356 }
 357
 358
 359 EvtTensor4C EvtGenFunctions::directProd(const EvtVector4R& c1,
 360                                         const EvtVector4R& c2){
 361
 362   EvtTensor4C temp;
 363   int i,j;
 364
 365   for (i=0;i<4;i++) {
 366     for (j=0;j<4;j++) {
 367       temp.t[i][j]=EvtComplex(c1.get(i)*c2.get(j),0.0);
 368     }
 369   }
 370   return temp;
 371 }
 372
 373 EvtTensor4C& EvtTensor4C::addDirProd(const EvtVector4R& p1,const EvtVector4R& p2){
 374
 375   int i,j;
 376
 377   for (i=0;i<4;i++) {
 378     for (j=0;j<4;j++) {
 379       t[i][j]+=p1.get(i)*p2.get(j);
 380     }
 381   }
 382   return *this;
 383 }
 384
 385
 386 EvtTensor4C dual(const EvtTensor4C& t2){
 387
 388   EvtTensor4C temp;
 389
 390   temp.set(0,0,EvtComplex(0.0,0.0));
 391   temp.set(1,1,EvtComplex(0.0,0.0));
 392   temp.set(2,2,EvtComplex(0.0,0.0));
 393   temp.set(3,3,EvtComplex(0.0,0.0));
 394
 395   temp.set(0,1,t2.get(3,2)-t2.get(2,3));
 396   temp.set(0,2,-t2.get(3,1)+t2.get(1,3));
 397   temp.set(0,3,t2.get(2,1)-t2.get(1,2));
 398
 399   temp.set(1,2,-t2.get(3,0)+t2.get(0,3));
 400   temp.set(1,3,t2.get(2,0)-t2.get(0,2));
 401
 402   temp.set(2,3,-t2.get(1,0)+t2.get(0,1));
 403
 404   temp.set(1,0,-temp.get(0,1));
 405   temp.set(2,0,-temp.get(0,2));
 406   temp.set(3,0,-temp.get(0,3));
 407
 408   temp.set(2,1,-temp.get(1,2));
 409   temp.set(3,1,-temp.get(1,3));
 410
 411   temp.set(3,2,-temp.get(2,3));
 412
 413   return temp;
 414
 415 }
 416
 417
 418 EvtTensor4C conj(const EvtTensor4C& t2) {
 419   EvtTensor4C temp;
 420
 421   int i,j;
 422
 423   for(i=0;i<4;i++){
 424     for(j=0;j<4;j++){
 425       temp.set(i,j,::conj((t2.get(i,j))));
 426     }
 427   }
 428
 429   return temp;
 430 }
 431
 432
 433 EvtTensor4C cont22(const EvtTensor4C& t1,const EvtTensor4C& t2){
 434   EvtTensor4C temp;
 435
 436   int i,j;
 437   EvtComplex c;
 438
 439   for(i=0;i<4;i++){
 440     for(j=0;j<4;j++){
 441       c=t1.get(i,0)*t2.get(j,0)-t1.get(i,1)*t2.get(j,1)
 442         -t1.get(i,2)*t2.get(j,2)-t1.get(i,3)*t2.get(j,3);
 443       temp.set(i,j,c);
 444     }
 445   }
 446
 447   return temp;
 448 }
 449
 450 EvtTensor4C cont11(const EvtTensor4C& t1,const EvtTensor4C& t2){
 451   EvtTensor4C temp;
 452
 453   int i,j;
 454   EvtComplex c;
 455
 456   for(i=0;i<4;i++){
 457     for(j=0;j<4;j++){
 458         c=t1.get(0,i)*t2.get(0,j)-t1.get(1,i)*t2.get(1,j)
 459           -t1.get(2,i)*t2.get(2,j)-t1.get(3,i)*t2.get(3,j);
 460         temp.set(i,j,c);
 461     }
 462   }
 463
 464   return temp;
 465 }
 466
 467
 468 EvtVector4C EvtTensor4C::cont1(const EvtVector4C& v4) const {
 469   EvtVector4C temp;
 470
 471   int i;
 472
 473   for(i=0;i<4;i++){
 474     temp.set(i,t[0][i]*v4.get(0)-t[1][i]*v4.get(1)
 475              -t[2][i]*v4.get(2)-t[3][i]*v4.get(3));
 476   }
 477
 478   return temp;
 479 }
 480
 481 EvtVector4C EvtTensor4C::cont2(const EvtVector4C& v4) const {
 482   EvtVector4C temp;
 483
 484   int i;
 485
 486   for(i=0;i<4;i++){
 487     temp.set(i,t[i][0]*v4.get(0)-t[i][1]*v4.get(1)
 488              -t[i][2]*v4.get(2)-t[i][3]*v4.get(3));
 489   }
 490
 491   return temp;
 492 }
 493
 494
 495 EvtVector4C EvtTensor4C::cont1(const EvtVector4R& v4) const {
 496   EvtVector4C temp;
 497
 498   int i;
 499
 500   for(i=0;i<4;i++){
 501     temp.set(i,t[0][i]*v4.get(0)-t[1][i]*v4.get(1)
 502              -t[2][i]*v4.get(2)-t[3][i]*v4.get(3));
 503   }
 504
 505   return temp;
 506 }
 507
 508
 509 EvtVector4C EvtTensor4C::cont2(const EvtVector4R& v4) const {
 510   EvtVector4C temp;
 511
 512   int i;
 513
 514   for(i=0;i<4;i++){
 515     temp.set(i,t[i][0]*v4.get(0)-t[i][1]*v4.get(1)
 516              -t[i][2]*v4.get(2)-t[i][3]*v4.get(3));
 517   }
 518
 519   return temp;
 520 }
 521
 522
 523
 524 void EvtTensor4C::applyRotateEuler(double phi,double theta,double ksi){
 525
 526   EvtComplex tt[4][4];
 527   double sp,st,sk,cp,ct,ck;
 528   double lambda[4][4];
 529
 530   sp=sin(phi);
 531   st=sin(theta);
 532   sk=sin(ksi);
 533   cp=cos(phi);
 534   ct=cos(theta);
 535   ck=cos(ksi);
 536
 537
 538   lambda[0][0]=1.0;
 539   lambda[0][1]=0.0;
 540   lambda[1][0]=0.0;
 541   lambda[0][2]=0.0;
 542   lambda[2][0]=0.0;
 543   lambda[0][3]=0.0;
 544   lambda[3][0]=0.0;
 545
 546   lambda[1][1]= ck*ct*cp-sk*sp;
 547   lambda[1][2]=-sk*ct*cp-ck*sp;
 548   lambda[1][3]=st*cp;
 549
 550   lambda[2][1]= ck*ct*sp+sk*cp;
 551   lambda[2][2]=-sk*ct*sp+ck*cp;
 552   lambda[2][3]=st*sp;
 553
 554   lambda[3][1]=-ck*st;
 555   lambda[3][2]=sk*st;
 556   lambda[3][3]=ct;
 557
 558
 559   int i,j,k;
 560
 561
 562   for(i=0;i<4;i++){
 563     for(j=0;j<4;j++){
 564       tt[i][j] = EvtComplex(0.0);
 565       for(k=0;k<4;k++){
 566         tt[i][j]+=lambda[j][k]*t[i][k];
 567       }
 568     }
 569   }
 570
 571   for(i=0;i<4;i++){
 572     for(j=0;j<4;j++){
 573       t[i][j] = EvtComplex(0.0);
 574       for(k=0;k<4;k++){
 575         t[i][j]+=lambda[i][k]*tt[k][j];
 576       }
 577     }
 578   }
 579
 580 }
 581
 582
 583