Line data Source code
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 0 : EvtDiracSpinor::~EvtDiracSpinor(){}
34 :
35 0 : EvtDiracSpinor::EvtDiracSpinor(const EvtComplex& sp0,const EvtComplex& sp1,
36 0 : const EvtComplex& sp2,const EvtComplex& sp3){
37 0 : set(sp0,sp1,sp2,sp3);
38 0 : }
39 :
40 : void EvtDiracSpinor::set(const EvtComplex& sp0,const EvtComplex& sp1,
41 : const EvtComplex& sp2,const EvtComplex& sp3){
42 :
43 0 : spinor[0]=sp0;spinor[1]=sp1;spinor[2]=sp2;spinor[3]=sp3;
44 0 : }
45 :
46 : void EvtDiracSpinor::set_spinor(int i,const EvtComplex& sp){
47 :
48 0 : spinor[i]=sp;
49 0 : }
50 :
51 : ostream& operator<<(ostream& s, const EvtDiracSpinor& sp){
52 :
53 0 : s <<"["<<sp.spinor[0]<<","<<sp.spinor[1]<<","
54 0 : <<sp.spinor[2]<<","<<sp.spinor[3]<<"]";
55 0 : return s;
56 :
57 : }
58 :
59 :
60 : const EvtComplex& EvtDiracSpinor::get_spinor(int i) const {
61 :
62 0 : return spinor[i];
63 :
64 : }
65 :
66 : EvtDiracSpinor rotateEuler(const EvtDiracSpinor& sp,
67 : double alpha,double beta,double gamma){
68 :
69 0 : EvtDiracSpinor tmp(sp);
70 0 : tmp.applyRotateEuler(alpha,beta,gamma);
71 : return tmp;
72 :
73 0 : }
74 :
75 : EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
76 : const EvtVector4R p4){
77 :
78 0 : EvtDiracSpinor tmp(sp);
79 0 : tmp.applyBoostTo(p4);
80 : return tmp;
81 :
82 0 : }
83 :
84 : EvtDiracSpinor boostTo(const EvtDiracSpinor& sp,
85 : const EvtVector3R boost){
86 :
87 0 : EvtDiracSpinor tmp(sp);
88 0 : tmp.applyBoostTo(boost);
89 : return tmp;
90 :
91 0 : }
92 :
93 : void EvtDiracSpinor::applyBoostTo(const EvtVector4R& p4){
94 :
95 0 : double e=p4.get(0);
96 :
97 0 : EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
98 :
99 0 : applyBoostTo(boost);
100 :
101 : return;
102 :
103 0 : }
104 :
105 :
106 :
107 : void EvtDiracSpinor::applyBoostTo(const EvtVector3R& boost) {
108 :
109 : double bx,by,bz,gamma,b2,f1,f2;
110 0 : EvtComplex spinorp[4];
111 :
112 0 : bx=boost.get(0);
113 0 : by=boost.get(1);
114 0 : bz=boost.get(2);
115 0 : b2=bx*bx+by*by+bz*bz;
116 :
117 0 : if (b2==0.0){
118 0 : return;
119 : }
120 :
121 : //assert(b2<1.0);
122 :
123 : gamma=1.0;
124 0 : if (b2 < 1.0) {gamma = 1.0/sqrt(1.0-b2);}
125 :
126 0 : f1=sqrt((gamma+1.0)/2.0);
127 0 : f2=f1*gamma/(gamma+1.0);
128 :
129 0 : spinorp[0]=f1*spinor[0]+f2*bz*spinor[2]+
130 0 : f2*EvtComplex(bx,-by)*spinor[3];
131 0 : spinorp[1]=f1*spinor[1]+f2*EvtComplex(bx,by)*spinor[2]-
132 0 : f2*bz*spinor[3];
133 0 : spinorp[2]=f2*bz*spinor[0]+f2*EvtComplex(bx,-by)*spinor[1]+
134 0 : f1*spinor[2];
135 0 : spinorp[3]=f2*EvtComplex(bx,by)*spinor[0]-
136 0 : f2*bz*spinor[1]+f1*spinor[3];
137 :
138 0 : spinor[0]=spinorp[0];
139 0 : spinor[1]=spinorp[1];
140 0 : spinor[2]=spinorp[2];
141 0 : spinor[3]=spinorp[3];
142 :
143 0 : return;
144 0 : }
145 :
146 : void EvtDiracSpinor::applyRotateEuler(double alpha,double beta,
147 : double gamma) {
148 :
149 0 : EvtComplex retVal[4];
150 :
151 0 : double cb2=cos(0.5*beta);
152 0 : double sb2=sin(0.5*beta);
153 0 : double capg2=cos(0.5*(alpha+gamma));
154 0 : double camg2=cos(0.5*(alpha-gamma));
155 0 : double sapg2=sin(0.5*(alpha+gamma));
156 0 : double samg2=sin(0.5*(alpha-gamma));
157 :
158 0 : EvtComplex m11(cb2*capg2,-cb2*sapg2);
159 0 : EvtComplex m12(-sb2*camg2,sb2*samg2);
160 0 : EvtComplex m21(sb2*camg2,sb2*samg2);
161 0 : EvtComplex m22(cb2*capg2,cb2*sapg2);
162 :
163 0 : retVal[0]=m11*spinor[0]+m12*spinor[1];
164 0 : retVal[1]=m21*spinor[0]+m22*spinor[1];
165 0 : retVal[2]=m11*spinor[2]+m12*spinor[3];
166 0 : retVal[3]=m21*spinor[2]+m22*spinor[3];
167 :
168 0 : spinor[0]=retVal[0];
169 0 : spinor[1]=retVal[1];
170 0 : spinor[2]=retVal[2];
171 0 : spinor[3]=retVal[3];
172 :
173 : return;
174 0 : }
175 :
176 :
177 :
178 : EvtDiracSpinor EvtDiracSpinor::conj() const {
179 :
180 0 : EvtDiracSpinor sp;
181 :
182 0 : for ( int i=0; i<4; i++)
183 0 : sp.set_spinor(i,::conj(spinor[i]));
184 :
185 : return sp;
186 0 : }
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 0 : EvtComplex u02=::conj(d.spinor[0]-d.spinor[2]);
205 0 : EvtComplex u13=::conj(d.spinor[1]-d.spinor[3]);
206 :
207 0 : EvtComplex v02=dp.spinor[0]-dp.spinor[2];
208 0 : EvtComplex v13=dp.spinor[1]-dp.spinor[3];
209 :
210 0 : EvtComplex a=u02*v02;
211 0 : EvtComplex b=u13*v13;
212 :
213 0 : EvtComplex c=u02*v13;
214 0 : EvtComplex e=u13*v02;
215 :
216 0 : return EvtVector4C(a+b,-(c+e),EvtComplex(0,1)*(c-e),b-a);
217 :
218 :
219 0 : }
220 :
221 : EvtVector4C EvtLeptonVCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
222 :
223 0 : 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 0 : temp.set(0,d*(EvtGammaMatrix::v0()*dp));
229 0 : temp.set(1,d*(EvtGammaMatrix::v1()*dp));
230 0 : temp.set(2,d*(EvtGammaMatrix::v2()*dp));
231 0 : temp.set(3,d*(EvtGammaMatrix::v3()*dp));
232 :
233 : return temp;
234 0 : }
235 :
236 :
237 : EvtVector4C EvtLeptonACurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
238 :
239 0 : EvtVector4C temp;
240 :
241 0 : 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 0 : mat = EvtGammaMatrix::v0()-EvtGammaMatrix::va0();
247 0 : temp.set(0,d*(mat*dp));
248 :
249 0 : mat = EvtGammaMatrix::v1()-EvtGammaMatrix::va1();
250 0 : temp.set(1,d*(mat*dp));
251 :
252 0 : mat = EvtGammaMatrix::v2()-EvtGammaMatrix::va2();
253 0 : temp.set(2,d*(mat*dp));
254 :
255 0 : mat = EvtGammaMatrix::v3()-EvtGammaMatrix::va3();
256 0 : temp.set(3,d*(mat*dp));
257 :
258 : return temp;
259 0 : }
260 :
261 : EvtComplex EvtLeptonSCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
262 :
263 0 : 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 0 : temp=d*(EvtGammaMatrix::g0()*dp);
269 :
270 0 : return temp;
271 0 : }
272 :
273 : EvtComplex EvtLeptonPCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
274 :
275 0 : 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 0 : static EvtGammaMatrix m=EvtGammaMatrix::g0()*EvtGammaMatrix::g5();
280 0 : temp=d*(m*dp);
281 :
282 0 : return temp;
283 0 : }
284 :
285 : EvtTensor4C EvtLeptonTCurrent(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
286 :
287 0 : EvtTensor4C temp;
288 0 : temp.zero();
289 0 : EvtComplex i2(0,0.5);
290 :
291 0 : static EvtGammaMatrix mat01=EvtGammaMatrix::g0()*
292 0 : (EvtGammaMatrix::g0()*EvtGammaMatrix::g1()-
293 0 : EvtGammaMatrix::g1()*EvtGammaMatrix::g0());
294 0 : static EvtGammaMatrix mat02=EvtGammaMatrix::g0()*
295 0 : (EvtGammaMatrix::g0()*EvtGammaMatrix::g2()-
296 0 : EvtGammaMatrix::g2()*EvtGammaMatrix::g0());
297 0 : static EvtGammaMatrix mat03=EvtGammaMatrix::g0()*
298 0 : (EvtGammaMatrix::g0()*EvtGammaMatrix::g3()-
299 0 : EvtGammaMatrix::g3()*EvtGammaMatrix::g0());
300 0 : static EvtGammaMatrix mat12=EvtGammaMatrix::g0()*
301 0 : (EvtGammaMatrix::g1()*EvtGammaMatrix::g2()-
302 0 : EvtGammaMatrix::g2()*EvtGammaMatrix::g1());
303 0 : static EvtGammaMatrix mat13=EvtGammaMatrix::g0()*
304 0 : (EvtGammaMatrix::g1()*EvtGammaMatrix::g3()-
305 0 : EvtGammaMatrix::g3()*EvtGammaMatrix::g1());
306 0 : static EvtGammaMatrix mat23=EvtGammaMatrix::g0()*
307 0 : (EvtGammaMatrix::g2()*EvtGammaMatrix::g3()-
308 0 : EvtGammaMatrix::g3()*EvtGammaMatrix::g2());
309 :
310 :
311 0 : temp.set(0,1,i2*(d*(mat01*dp)));
312 0 : temp.set(1,0,-temp.get(0,1));
313 :
314 0 : temp.set(0,2,i2*(d*(mat02*dp)));
315 0 : temp.set(2,0,-temp.get(0,2));
316 :
317 0 : temp.set(0,3,i2*(d*(mat03*dp)));
318 0 : temp.set(3,0,-temp.get(0,3));
319 :
320 0 : temp.set(1,2,i2*(d*(mat12*dp)));
321 0 : temp.set(2,1,-temp.get(1,2));
322 :
323 0 : temp.set(1,3,i2*(d*(mat13*dp)));
324 0 : temp.set(3,1,-temp.get(1,3));
325 :
326 0 : temp.set(2,3,i2*(d*(mat23*dp)));
327 0 : temp.set(3,2,-temp.get(2,3));
328 :
329 : return temp;
330 0 : }
331 :
332 :
333 : EvtDiracSpinor operator*(const EvtComplex& c, const EvtDiracSpinor& d) {
334 0 : EvtDiracSpinor result;
335 0 : result.spinor[0] = c*d.spinor[0];
336 0 : result.spinor[1] = c*d.spinor[1];
337 0 : result.spinor[2] = c*d.spinor[2];
338 0 : result.spinor[3] = c*d.spinor[3];
339 : return result;
340 0 : }
341 :
342 : EvtDiracSpinor EvtDiracSpinor::adjoint() const
343 : {
344 0 : EvtDiracSpinor d = this->conj(); // first conjugate, then multiply with gamma0
345 0 : EvtGammaMatrix g0 = EvtGammaMatrix::g0();
346 0 : EvtDiracSpinor result; // automatically initialized to 0
347 :
348 0 : for (int i=0; i<4; ++i)
349 0 : for (int j=0; j<4; ++j)
350 0 : result.spinor[i] += d.spinor[j] * g0._gamma[i][j];
351 :
352 : return result;
353 0 : }
354 :
355 : EvtComplex operator*(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
356 :
357 : int i;
358 0 : EvtComplex temp;
359 :
360 0 : temp=EvtComplex(0.0,0.0);
361 :
362 0 : for(i=0;i<4;i++){
363 0 : temp += conj( d.get_spinor(i) ) * dp.get_spinor( i ) ;
364 : }
365 : return temp;
366 0 : }
|
