TEvtGen/EvtGenBase/EvtDiracSpinor.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: EvtDiracSpinor.cc
  12 //
  13 // Description:  Class to describe (EvtDiracParticle) spinors.
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     September 25, 1996         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <math.h>
  23 #include <assert.h>
  24 #include "EvtGenBase/EvtDiracSpinor.hh"
  25 #include "EvtGenBase/EvtGammaMatrix.hh"
  26 #include "EvtGenBase/EvtComplex.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenBase/EvtVector4C.hh"
  29 #include "EvtGenBase/EvtTensor4C.hh"
  30 using std::ostream;
  31
  32
  33 EvtDiracSpinor::~EvtDiracSpinor(){}
  34
  35 EvtDiracSpinor::EvtDiracSpinor(const EvtComplex& sp0,const EvtComplex& sp1,
  36                                     const EvtComplex& sp2,const EvtComplex& sp3){
  37   set(sp0,sp1,sp2,sp3);
  38 }
  39
  40 void EvtDiracSpinor::set(const EvtComplex& sp0,const EvtComplex& sp1,
  41                          const EvtComplex& sp2,const EvtComplex& sp3){
  42
  43   spinor[0]=sp0;spinor[1]=sp1;spinor[2]=sp2;spinor[3]=sp3;
  44 }
  45
  46 void EvtDiracSpinor::set_spinor(int i,const EvtComplex& sp){
  47
  48   spinor[i]=sp;
  49 }
  50
  51 ostream& operator<<(ostream& s, const EvtDiracSpinor& sp){
  52
  53   s <<"["<<sp.spinor[0]<<","<<sp.spinor[1]<<","
  54     <<sp.spinor[2]<<","<<sp.spinor[3]<<"]";
  55   return s;
  56
  57 }
  58
  59
  60 const EvtComplex& EvtDiracSpinor::get_spinor(int i) const {
  61
  62   return spinor[i];
  63
  64 }
  65
  66 EvtDiracSpinor rotateEuler(const EvtDiracSpinor& sp,
  67                            double alpha,double beta,double gamma){
  68
  69   EvtDiracSpinor tmp(sp);
  70   tmp.applyRotateEuler(alpha,beta,gamma);
  71   return tmp;
  72
  73 }
  74
  75 EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
  76                        const EvtVector4R p4){
  77
  78   EvtDiracSpinor tmp(sp);
  79   tmp.applyBoostTo(p4);
  80   return tmp;
  81
  82 }
  83
  84 EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
  85                        const EvtVector3R boost){
  86
  87   EvtDiracSpinor tmp(sp);
  88   tmp.applyBoostTo(boost);
  89   return tmp;
  90
  91 }
  92
  93 void EvtDiracSpinor::applyBoostTo(const EvtVector4R& p4){
  94
  95   double e=p4.get(0);
  96
  97   EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
  98
  99   applyBoostTo(boost);
 100
 101   return;
 102
 103 }
 104
 105
 106
 107 void EvtDiracSpinor::applyBoostTo(const EvtVector3R& boost) {
 108
 109   double bx,by,bz,gamma,b2,f1,f2;
 110   EvtComplex spinorp[4];
 111
 112   bx=boost.get(0);
 113   by=boost.get(1);
 114   bz=boost.get(2);
 115   b2=bx*bx+by*by+bz*bz;
 116
 117   if (b2==0.0){
 118     return;
 119   }
 120
 121   //assert(b2<1.0);
 122
 123   gamma=1.0;
 124   if (b2 < 1.0) {gamma = 1.0/sqrt(1.0-b2);}
 125
 126   f1=sqrt((gamma+1.0)/2.0);
 127   f2=f1*gamma/(gamma+1.0);
 128
 129   spinorp[0]=f1*spinor[0]+f2*bz*spinor[2]+
 130     f2*EvtComplex(bx,-by)*spinor[3];
 131   spinorp[1]=f1*spinor[1]+f2*EvtComplex(bx,by)*spinor[2]-
 132     f2*bz*spinor[3];
 133   spinorp[2]=f2*bz*spinor[0]+f2*EvtComplex(bx,-by)*spinor[1]+
 134     f1*spinor[2];
 135   spinorp[3]=f2*EvtComplex(bx,by)*spinor[0]-
 136     f2*bz*spinor[1]+f1*spinor[3];
 137
 138   spinor[0]=spinorp[0];
 139   spinor[1]=spinorp[1];
 140   spinor[2]=spinorp[2];
 141   spinor[3]=spinorp[3];
 142
 143   return;
 144 }
 145
 146 void EvtDiracSpinor::applyRotateEuler(double alpha,double beta,
 147                                       double gamma) {
 148
 149   EvtComplex retVal[4];
 150
 151   double cb2=cos(0.5*beta);
 152   double sb2=sin(0.5*beta);
 153   double capg2=cos(0.5*(alpha+gamma));
 154   double camg2=cos(0.5*(alpha-gamma));
 155   double sapg2=sin(0.5*(alpha+gamma));
 156   double samg2=sin(0.5*(alpha-gamma));
 157
 158   EvtComplex m11(cb2*capg2,-cb2*sapg2);
 159   EvtComplex m12(-sb2*camg2,sb2*samg2);
 160   EvtComplex m21(sb2*camg2,sb2*samg2);
 161   EvtComplex m22(cb2*capg2,cb2*sapg2);
 162
 163   retVal[0]=m11*spinor[0]+m12*spinor[1];
 164   retVal[1]=m21*spinor[0]+m22*spinor[1];
 165   retVal[2]=m11*spinor[2]+m12*spinor[3];
 166   retVal[3]=m21*spinor[2]+m22*spinor[3];
 167
 168   spinor[0]=retVal[0];
 169   spinor[1]=retVal[1];
 170   spinor[2]=retVal[2];
 171   spinor[3]=retVal[3];
 172
 173   return;
 174 }
 175
 176
 177
 178 EvtDiracSpinor EvtDiracSpinor::conj() const {
 179
 180   EvtDiracSpinor sp;
 181
 182   for ( int i=0; i<4; i++)
 183     sp.set_spinor(i,::conj(spinor[i]));
 184
 185   return sp;
 186 }
 187
 188 EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 189
 190   //Old code; below is a new specialized code that does it more efficiently.
 191   //EvtGammaMatrix mat;
 192   //EvtVector4C temp;
 193   //mat.va0();
 194   //temp.set(0,d*(mat*dp));
 195   //mat.va1();
 196   //temp.set(1,d*(mat*dp));
 197   //mat.va2();
 198   //temp.set(2,d*(mat*dp));
 199   //mat.va3();
 200   //temp.set(3,d*(mat*dp));
 201   //return temp;
 202
 203
 204   EvtComplex u02=::conj(d.spinor[0]-d.spinor[2]);
 205   EvtComplex u13=::conj(d.spinor[1]-d.spinor[3]);
 206
 207   EvtComplex v02=dp.spinor[0]-dp.spinor[2];
 208   EvtComplex v13=dp.spinor[1]-dp.spinor[3];
 209
 210   EvtComplex a=u02*v02;
 211   EvtComplex b=u13*v13;
 212
 213   EvtComplex c=u02*v13;
 214   EvtComplex e=u13*v02;
 215
 216   return EvtVector4C(a+b,-(c+e),EvtComplex(0,1)*(c-e),b-a);
 217
 218
 219 }
 220
 221 EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 222
 223   EvtVector4C temp;
 224
 225   // no conjugate here; done in the multiplication
 226   // yes this is stupid and fooled me to for a long time (ryd)
 227
 228   temp.set(0,d*(EvtGammaMatrix::v0()*dp));
 229   temp.set(1,d*(EvtGammaMatrix::v1()*dp));
 230   temp.set(2,d*(EvtGammaMatrix::v2()*dp));
 231   temp.set(3,d*(EvtGammaMatrix::v3()*dp));
 232
 233   return temp;
 234 }
 235
 236
 237 EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 238
 239   EvtVector4C temp;
 240
 241   EvtGammaMatrix mat;
 242
 243   // no conjugate here; done in the multiplication
 244   // yes this is stupid and fooled me to for a long time (ryd)
 245
 246   mat = EvtGammaMatrix::v0()-EvtGammaMatrix::va0();
 247   temp.set(0,d*(mat*dp));
 248
 249   mat = EvtGammaMatrix::v1()-EvtGammaMatrix::va1();
 250   temp.set(1,d*(mat*dp));
 251
 252   mat = EvtGammaMatrix::v2()-EvtGammaMatrix::va2();
 253   temp.set(2,d*(mat*dp));
 254
 255   mat = EvtGammaMatrix::v3()-EvtGammaMatrix::va3();
 256   temp.set(3,d*(mat*dp));
 257
 258   return temp;
 259 }
 260
 261 EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 262
 263   EvtComplex temp;
 264
 265   // no conjugate here; done in the multiplication
 266   // yes this is stupid and fooled me to for a long time (ryd)
 267
 268   temp=d*(EvtGammaMatrix::g0()*dp);
 269
 270   return temp;
 271 }
 272
 273 EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 274
 275   EvtComplex temp;
 276
 277   // no conjugate here; done in the multiplication
 278   // yes this is stupid and fooled me to for a long time (ryd)
 279   static EvtGammaMatrix m=EvtGammaMatrix::g0()*EvtGammaMatrix::g5();
 280   temp=d*(m*dp);
 281
 282   return temp;
 283 }
 284
 285 EvtTensor4C EvtLeptonTCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 286
 287   EvtTensor4C temp;
 288   temp.zero();
 289   EvtComplex i2(0,0.5);
 290
 291   static EvtGammaMatrix mat01=EvtGammaMatrix::g0()*
 292     (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()-
 293      EvtGammaMatrix::g1()*EvtGammaMatrix::g0());
 294   static EvtGammaMatrix mat02=EvtGammaMatrix::g0()*
 295     (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()-
 296      EvtGammaMatrix::g2()*EvtGammaMatrix::g0());
 297   static EvtGammaMatrix mat03=EvtGammaMatrix::g0()*
 298     (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()-
 299     EvtGammaMatrix::g3()*EvtGammaMatrix::g0());
 300   static EvtGammaMatrix mat12=EvtGammaMatrix::g0()*
 301     (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()-
 302     EvtGammaMatrix::g2()*EvtGammaMatrix::g1());
 303   static EvtGammaMatrix mat13=EvtGammaMatrix::g0()*
 304     (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()-
 305      EvtGammaMatrix::g3()*EvtGammaMatrix::g1());
 306   static EvtGammaMatrix mat23=EvtGammaMatrix::g0()*
 307     (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()-
 308      EvtGammaMatrix::g3()*EvtGammaMatrix::g2());
 309
 310
 311   temp.set(0,1,i2*(d*(mat01*dp)));
 312   temp.set(1,0,-temp.get(0,1));
 313
 314   temp.set(0,2,i2*(d*(mat02*dp)));
 315   temp.set(2,0,-temp.get(0,2));
 316
 317   temp.set(0,3,i2*(d*(mat03*dp)));
 318   temp.set(3,0,-temp.get(0,3));
 319
 320   temp.set(1,2,i2*(d*(mat12*dp)));
 321   temp.set(2,1,-temp.get(1,2));
 322
 323   temp.set(1,3,i2*(d*(mat13*dp)));
 324   temp.set(3,1,-temp.get(1,3));
 325
 326   temp.set(2,3,i2*(d*(mat23*dp)));
 327   temp.set(3,2,-temp.get(2,3));
 328
 329   return temp;
 330 }
 331
 332
 333 EvtDiracSpinor operator*(const EvtComplex& c, const EvtDiracSpinor& d) {
 334      EvtDiracSpinor result;
 335      result.spinor[0] = c*d.spinor[0];
 336      result.spinor[1] = c*d.spinor[1];
 337      result.spinor[2] = c*d.spinor[2];
 338      result.spinor[3] = c*d.spinor[3];
 339      return result;
 340  }
 341
 342 EvtDiracSpinor EvtDiracSpinor::adjoint() const
 343 {
 344     EvtDiracSpinor d = this->conj(); // first conjugate, then multiply with gamma0
 345     EvtGammaMatrix g0 = EvtGammaMatrix::g0();
 346     EvtDiracSpinor result; // automatically initialized to 0
 347
 348     for (int i=0; i<4; ++i)
 349         for (int j=0; j<4; ++j)
 350             result.spinor[i] += d.spinor[j] * g0._gamma[i][j];
 351
 352     return result;
 353 }
 354
 355 EvtComplex operator*(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
 356
 357   int i;
 358   EvtComplex temp;
 359
 360   temp=EvtComplex(0.0,0.0);
 361
 362   for(i=0;i<4;i++){
 363     temp += conj( d.get_spinor(i) ) * dp.get_spinor( i ) ;
 364   }
 365   return temp;
 366 }