1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtDiracSpinor.cc
13 // Description: Class to describe (EvtDiracParticle) spinors.
15 // Modification history:
17 // DJL/RYD September 25, 1996 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
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"
33 EvtDiracSpinor::~EvtDiracSpinor(){}
35 EvtDiracSpinor::EvtDiracSpinor(const EvtComplex& sp0,const EvtComplex& sp1,
36 const EvtComplex& sp2,const EvtComplex& sp3){
40 void EvtDiracSpinor::set(const EvtComplex& sp0,const EvtComplex& sp1,
41 const EvtComplex& sp2,const EvtComplex& sp3){
43 spinor[0]=sp0;spinor[1]=sp1;spinor[2]=sp2;spinor[3]=sp3;
46 void EvtDiracSpinor::set_spinor(int i,const EvtComplex& sp){
51 ostream& operator<<(ostream& s, const EvtDiracSpinor& sp){
53 s <<"["<<sp.spinor[0]<<","<<sp.spinor[1]<<","
54 <<sp.spinor[2]<<","<<sp.spinor[3]<<"]";
60 const EvtComplex& EvtDiracSpinor::get_spinor(int i) const {
66 EvtDiracSpinor rotateEuler(const EvtDiracSpinor& sp,
67 double alpha,double beta,double gamma){
69 EvtDiracSpinor tmp(sp);
70 tmp.applyRotateEuler(alpha,beta,gamma);
75 EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
76 const EvtVector4R p4){
78 EvtDiracSpinor tmp(sp);
84 EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
85 const EvtVector3R boost){
87 EvtDiracSpinor tmp(sp);
88 tmp.applyBoostTo(boost);
93 void EvtDiracSpinor::applyBoostTo(const EvtVector4R& p4){
97 EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
107 void EvtDiracSpinor::applyBoostTo(const EvtVector3R& boost) {
109 double bx,by,bz,gamma,b2,f1,f2;
110 EvtComplex spinorp[4];
115 b2=bx*bx+by*by+bz*bz;
123 gamma=1.0/sqrt(1-b2);
125 f1=sqrt((gamma+1.0)/2.0);
126 f2=f1*gamma/(gamma+1.0);
128 spinorp[0]=f1*spinor[0]+f2*bz*spinor[2]+
129 f2*EvtComplex(bx,-by)*spinor[3];
130 spinorp[1]=f1*spinor[1]+f2*EvtComplex(bx,by)*spinor[2]-
132 spinorp[2]=f2*bz*spinor[0]+f2*EvtComplex(bx,-by)*spinor[1]+
134 spinorp[3]=f2*EvtComplex(bx,by)*spinor[0]-
135 f2*bz*spinor[1]+f1*spinor[3];
137 spinor[0]=spinorp[0];
138 spinor[1]=spinorp[1];
139 spinor[2]=spinorp[2];
140 spinor[3]=spinorp[3];
145 void EvtDiracSpinor::applyRotateEuler(double alpha,double beta,
148 EvtComplex retVal[4];
150 double cb2=cos(0.5*beta);
151 double sb2=sin(0.5*beta);
152 double capg2=cos(0.5*(alpha+gamma));
153 double camg2=cos(0.5*(alpha-gamma));
154 double sapg2=sin(0.5*(alpha+gamma));
155 double samg2=sin(0.5*(alpha-gamma));
157 EvtComplex m11(cb2*capg2,-cb2*sapg2);
158 EvtComplex m12(-sb2*camg2,sb2*samg2);
159 EvtComplex m21(sb2*camg2,sb2*samg2);
160 EvtComplex m22(cb2*capg2,cb2*sapg2);
162 retVal[0]=m11*spinor[0]+m12*spinor[1];
163 retVal[1]=m21*spinor[0]+m22*spinor[1];
164 retVal[2]=m11*spinor[2]+m12*spinor[3];
165 retVal[3]=m21*spinor[2]+m22*spinor[3];
177 EvtDiracSpinor EvtDiracSpinor::conj() const {
181 for ( int i=0; i<4; i++)
182 sp.set_spinor(i,::conj(spinor[i]));
187 EvtVector4C EvtLeptonVACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
189 //Old code; below is a new specialized code that does it more efficiently.
190 //EvtGammaMatrix mat;
193 //temp.set(0,d*(mat*dp));
195 //temp.set(1,d*(mat*dp));
197 //temp.set(2,d*(mat*dp));
199 //temp.set(3,d*(mat*dp));
203 EvtComplex u02=::conj(d.spinor[0]-d.spinor[2]);
204 EvtComplex u13=::conj(d.spinor[1]-d.spinor[3]);
206 EvtComplex v02=dp.spinor[0]-dp.spinor[2];
207 EvtComplex v13=dp.spinor[1]-dp.spinor[3];
209 EvtComplex a=u02*v02;
210 EvtComplex b=u13*v13;
212 EvtComplex c=u02*v13;
213 EvtComplex e=u13*v02;
215 return EvtVector4C(a+b,-(c+e),EvtComplex(0,1)*(c-e),b-a);
220 EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
224 // no conjugate here; done in the multiplication
225 // yes this is stupid and fooled me to for a long time (ryd)
227 temp.set(0,d*(EvtGammaMatrix::v0()*dp));
228 temp.set(1,d*(EvtGammaMatrix::v1()*dp));
229 temp.set(2,d*(EvtGammaMatrix::v2()*dp));
230 temp.set(3,d*(EvtGammaMatrix::v3()*dp));
236 EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
242 // no conjugate here; done in the multiplication
243 // yes this is stupid and fooled me to for a long time (ryd)
245 mat = EvtGammaMatrix::v0()-EvtGammaMatrix::va0();
246 temp.set(0,d*(mat*dp));
248 mat = EvtGammaMatrix::v1()-EvtGammaMatrix::va1();
249 temp.set(1,d*(mat*dp));
251 mat = EvtGammaMatrix::v2()-EvtGammaMatrix::va2();
252 temp.set(2,d*(mat*dp));
254 mat = EvtGammaMatrix::v3()-EvtGammaMatrix::va3();
255 temp.set(3,d*(mat*dp));
260 EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
264 // no conjugate here; done in the multiplication
265 // yes this is stupid and fooled me to for a long time (ryd)
267 temp=d*(EvtGammaMatrix::g0()*dp);
272 EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
276 // no conjugate here; done in the multiplication
277 // yes this is stupid and fooled me to for a long time (ryd)
278 static EvtGammaMatrix m=EvtGammaMatrix::g0()*EvtGammaMatrix::g5();
284 EvtTensor4C EvtLeptonTCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
288 EvtComplex i2(0,0.5);
290 static EvtGammaMatrix mat01=EvtGammaMatrix::g0()*
291 (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()-
292 EvtGammaMatrix::g1()*EvtGammaMatrix::g0());
293 static EvtGammaMatrix mat02=EvtGammaMatrix::g0()*
294 (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()-
295 EvtGammaMatrix::g2()*EvtGammaMatrix::g0());
296 static EvtGammaMatrix mat03=EvtGammaMatrix::g0()*
297 (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()-
298 EvtGammaMatrix::g3()*EvtGammaMatrix::g0());
299 static EvtGammaMatrix mat12=EvtGammaMatrix::g0()*
300 (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()-
301 EvtGammaMatrix::g2()*EvtGammaMatrix::g1());
302 static EvtGammaMatrix mat13=EvtGammaMatrix::g0()*
303 (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()-
304 EvtGammaMatrix::g3()*EvtGammaMatrix::g1());
305 static EvtGammaMatrix mat23=EvtGammaMatrix::g0()*
306 (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()-
307 EvtGammaMatrix::g3()*EvtGammaMatrix::g2());
310 temp.set(0,1,i2*(d*(mat01*dp)));
311 temp.set(1,0,-temp.get(0,1));
313 temp.set(0,2,i2*(d*(mat02*dp)));
314 temp.set(2,0,-temp.get(0,2));
316 temp.set(0,3,i2*(d*(mat03*dp)));
317 temp.set(3,0,-temp.get(0,3));
319 temp.set(1,2,i2*(d*(mat12*dp)));
320 temp.set(2,1,-temp.get(1,2));
322 temp.set(1,3,i2*(d*(mat13*dp)));
323 temp.set(3,1,-temp.get(1,3));
325 temp.set(2,3,i2*(d*(mat23*dp)));
326 temp.set(3,2,-temp.get(2,3));
332 EvtDiracSpinor operator*(const EvtComplex& c, const EvtDiracSpinor& d) {
333 EvtDiracSpinor result;
334 result.spinor[0] = c*d.spinor[0];
335 result.spinor[1] = c*d.spinor[1];
336 result.spinor[2] = c*d.spinor[2];
337 result.spinor[3] = c*d.spinor[3];
342 EvtComplex operator*(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
347 temp=EvtComplex(0.0,0.0);
350 temp+=::conj(d.get_spinor(i))*dp.get_spinor(i);
356 EvtDiracSpinor EvtDiracSpinor::adjoint() const
358 EvtDiracSpinor d = this->conj(); // first conjugate, then multiply with gamma0
359 EvtGammaMatrix g0 = EvtGammaMatrix::g0();
360 EvtDiracSpinor result; // automatically initialized to 0
362 for (int i=0; i<4; ++i)
363 for (int j=0; j<4; ++j)
364 result.spinor[i] += d.spinor[j] * g0._gamma[i][j];