]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenBase/EvtGammaMatrix.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coding conventions
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtGammaMatrix.cxx
CommitLineData
da0e9ce3 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: EvtGammaMatrix.cc
12//
13// Description: Make gamma matrices availible for the calc. of amplitudes, etc.
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/EvtGammaMatrix.hh"
27#include "EvtGenBase/EvtDiracSpinor.hh"
28#include "EvtGenBase/EvtReport.hh"
29#include "EvtGenBase/EvtTensor4C.hh"
30#include "EvtGenBase/EvtVector4C.hh"
31#include <stdlib.h>
32using std::endl;
33using std::ostream;
34
35EvtGammaMatrix::EvtGammaMatrix(){
36 int i,j;
37
38 static EvtComplex zero(0.0,0.0);
39
40 for(i=0;i<4;i++){
41 for(j=0;j<4;j++){
42 _gamma[i][j]=zero;
43 }
44 }
45}
46
47EvtGammaMatrix operator*(const EvtGammaMatrix& g, const EvtComplex& c)
48{
49 return c*g;
50}
51
52
53EvtGammaMatrix operator*(const EvtComplex& c,const EvtGammaMatrix& g){
54 int i,j;
55
56 EvtGammaMatrix temp;
57
58 for(i=0;i<4;i++){
59 for(j=0;j<4;j++){
60 temp._gamma[i][j]=g._gamma[i][j]*c;
61 }
62 }
63
64 return temp;
65
66}
67
68
69ostream& operator<<(ostream& s, const EvtGammaMatrix& g){
70
71
72 s<<"["<<g._gamma[0][0]<<","<<g._gamma[0][1]<<","<<g._gamma[0][2]<<","<<g._gamma[0][3]<<"]"<<endl;
73 s<<"["<<g._gamma[1][0]<<","<<g._gamma[1][1]<<","<<g._gamma[1][2]<<","<<g._gamma[1][3]<<"]"<<endl;
74 s<<"["<<g._gamma[2][0]<<","<<g._gamma[2][1]<<","<<g._gamma[2][2]<<","<<g._gamma[2][3]<<"]"<<endl;
75 s<<"["<<g._gamma[3][0]<<","<<g._gamma[3][1]<<","<<g._gamma[3][2]<<","<<g._gamma[3][3]<<"]"<<endl;
76
77 return s;
78
79}
80
81
82
83EvtGammaMatrix::EvtGammaMatrix(const EvtGammaMatrix& gm){
84 int i,j;
85
86 for(i=0;i<4;i++){
87 for(j=0;j<4;j++){
88 _gamma[i][j]=gm._gamma[i][j];
89 }
90 }
91}
92
93EvtGammaMatrix::~EvtGammaMatrix() {}
94
95EvtGammaMatrix& EvtGammaMatrix::operator=(const EvtGammaMatrix& gm){
96 int i,j;
97
98 for(i=0;i<4;i++){
99 for(j=0;j<4;j++){
100 _gamma[i][j]=gm._gamma[i][j];
101 }
102 }
103 return *this;
104}
105
106void EvtGammaMatrix::init(){
107 int i,j;
108
109 static EvtComplex zero(0.0,0.0);
110
111 for(i=0;i<4;i++){
112 for(j=0;j<4;j++){
113 _gamma[i][j]=zero;
114 }
115 }
116}
117
118const EvtGammaMatrix& EvtGammaMatrix::va0(){
119
120 static EvtGammaMatrix g;
121 static int first=1;
122
123 if (first){
124 first = 0;
125 g._gamma[0][0]=EvtComplex(1.0,0.0);
126 g._gamma[0][1]=EvtComplex(0.0,0.0);
127 g._gamma[0][2]=EvtComplex(-1.0,0.0);
128 g._gamma[0][3]=EvtComplex(0.0,0.0);
129 g._gamma[1][0]=EvtComplex(0.0,0.0);
130 g._gamma[1][1]=EvtComplex(1.0,0.0);
131 g._gamma[1][2]=EvtComplex(0.0,0.0);
132 g._gamma[1][3]=EvtComplex(-1.0,0.0);
133 g._gamma[2][0]=EvtComplex(-1.0,0.0);
134 g._gamma[2][1]=EvtComplex(0.0,0.0);
135 g._gamma[2][2]=EvtComplex(1.0,0.0);
136 g._gamma[2][3]=EvtComplex(0.0,0.0);
137 g._gamma[3][0]=EvtComplex(0.0,0.0);
138 g._gamma[3][1]=EvtComplex(-1.0,0.0);
139 g._gamma[3][2]=EvtComplex(0.0,0.0);
140 g._gamma[3][3]=EvtComplex(1.0,0.0);
141 }
142
143 return g;
144
145}
146
147
148const EvtGammaMatrix& EvtGammaMatrix::va1(){
149
150 static EvtGammaMatrix g;
151 static int first=1;
152
153 if (first){
154 first = 0;
155 g._gamma[0][0]=EvtComplex(0.0,0.0);
156 g._gamma[0][1]=EvtComplex(-1.0,0.0);
157 g._gamma[0][2]=EvtComplex(0.0,0.0);
158 g._gamma[0][3]=EvtComplex(1.0,0.0);
159 g._gamma[1][0]=EvtComplex(-1.0,0.0);
160 g._gamma[1][1]=EvtComplex(0.0,0.0);
161 g._gamma[1][2]=EvtComplex(1.0,0.0);
162 g._gamma[1][3]=EvtComplex(0.0,0.0);
163 g._gamma[2][0]=EvtComplex(0.0,0.0);
164 g._gamma[2][1]=EvtComplex(1.0,0.0);
165 g._gamma[2][2]=EvtComplex(0.0,0.0);
166 g._gamma[2][3]=EvtComplex(-1.0,0.0);
167 g._gamma[3][0]=EvtComplex(1.0,0.0);
168 g._gamma[3][1]=EvtComplex(0.0,0.0);
169 g._gamma[3][2]=EvtComplex(-1.0,0.0);
170 g._gamma[3][3]=EvtComplex(0.0,0.0);
171 }
172
173 return g;
174
175}
176
177
178
179const EvtGammaMatrix& EvtGammaMatrix::va2(){
180
181 static EvtGammaMatrix g;
182 static int first=1;
183
184 if (first){
185 first = 0;
186 g._gamma[0][0]=EvtComplex(0.0,0.0);
187 g._gamma[0][1]=EvtComplex(0.0,1.0);
188 g._gamma[0][2]=EvtComplex(0.0,0.0);
189 g._gamma[0][3]=EvtComplex(0.0,-1.0);
190 g._gamma[1][0]=EvtComplex(0.0,-1.0);
191 g._gamma[1][1]=EvtComplex(0.0,0.0);
192 g._gamma[1][2]=EvtComplex(0.0,1.0);
193 g._gamma[1][3]=EvtComplex(0.0,0.0);
194 g._gamma[2][0]=EvtComplex(0.0,0.0);
195 g._gamma[2][1]=EvtComplex(0.0,-1.0);
196 g._gamma[2][2]=EvtComplex(0.0,0.0);
197 g._gamma[2][3]=EvtComplex(0.0,1.0);
198 g._gamma[3][0]=EvtComplex(0.0,1.0);
199 g._gamma[3][1]=EvtComplex(0.0,0.0);
200 g._gamma[3][2]=EvtComplex(0.0,-1.0);
201 g._gamma[3][3]=EvtComplex(0.0,0.0);
202 }
203
204 return g;
205
206}
207
208
209
210
211const EvtGammaMatrix& EvtGammaMatrix::va3(){
212
213 static EvtGammaMatrix g;
214 static int first=1;
215
216 if (first){
217 first = 0;
218 g._gamma[0][0]=EvtComplex(-1.0,0.0);
219 g._gamma[0][1]=EvtComplex(0.0,0.0);
220 g._gamma[0][2]=EvtComplex(1.0,0.0);
221 g._gamma[0][3]=EvtComplex(0.0,0.0);
222 g._gamma[1][0]=EvtComplex(0.0,0.0);
223 g._gamma[1][1]=EvtComplex(1.0,0.0);
224 g._gamma[1][2]=EvtComplex(0.0,0.0);
225 g._gamma[1][3]=EvtComplex(-1.0,0.0);
226 g._gamma[2][0]=EvtComplex(1.0,0.0);
227 g._gamma[2][1]=EvtComplex(0.0,0.0);
228 g._gamma[2][2]=EvtComplex(-1.0,0.0);
229 g._gamma[2][3]=EvtComplex(0.0,0.0);
230 g._gamma[3][0]=EvtComplex(0.0,0.0);
231 g._gamma[3][1]=EvtComplex(-1.0,0.0);
232 g._gamma[3][2]=EvtComplex(0.0,0.0);
233 g._gamma[3][3]=EvtComplex(1.0,0.0);
234 }
235
236 return g;
237
238}
239
240
241
242
243
244const EvtGammaMatrix& EvtGammaMatrix::g0(){
245
246 static EvtGammaMatrix g;
247 static int first=1;
248
249 if (first){
250
251 first=0;
252
253 int i,j;
254
255 for(i=0;i<4;i++){
256 for(j=0;j<4;j++){
257 g._gamma[i][j]=EvtComplex(0.0,0.0);
258 }
259 }
260
261 g._gamma[0][0]=EvtComplex(1.0,0.0);
262 g._gamma[1][1]=EvtComplex(1.0,0.0);
263 g._gamma[2][2]=EvtComplex(-1.0,0.0);
264 g._gamma[3][3]=EvtComplex(-1.0,0.0);
265 }
266
267 return g;
268
269}
270
271
272
273
274const EvtGammaMatrix& EvtGammaMatrix::g1(){
275
276 static EvtGammaMatrix g;
277 static int first=1;
278
279 if (first){
280 first=0;
281 int i,j;
282
283 for(i=0;i<4;i++){
284 for(j=0;j<4;j++){
285 g._gamma[i][j]=EvtComplex(0.0,0.0);
286 }
287 }
288
289 g._gamma[0][3]=EvtComplex(1.0,0.0);
290 g._gamma[1][2]=EvtComplex(1.0,0.0);
291 g._gamma[2][1]=EvtComplex(-1.0,0.0);
292 g._gamma[3][0]=EvtComplex(-1.0,0.0);
293 }
294
295 return g;
296
297}
298
299
300
301
302const EvtGammaMatrix& EvtGammaMatrix::g2(){
303
304 static EvtGammaMatrix g;
305 static int first=1;
306
307 if (first){
308 first=0;
309 int i,j;
310
311 for(i=0;i<4;i++){
312 for(j=0;j<4;j++){
313 g._gamma[i][j]=EvtComplex(0.0,0.0);
314 }
315 }
316
317 g._gamma[0][3]=EvtComplex(0.0,-1.0);
318 g._gamma[1][2]=EvtComplex(0.0,1.0);
319 g._gamma[2][1]=EvtComplex(0.0,1.0);
320 g._gamma[3][0]=EvtComplex(0.0,-1.0);
321 }
322
323 return g;
324
325}
326
327
328
329
330
331const EvtGammaMatrix& EvtGammaMatrix::g3(){
332
333 static EvtGammaMatrix g;
334 static int first=1;
335
336 if (first){
337 first=0;
338 int i,j;
339
340 for(i=0;i<4;i++){
341 for(j=0;j<4;j++){
342 g._gamma[i][j]=EvtComplex(0.0,0.0);
343 }
344 }
345
346 g._gamma[0][2]=EvtComplex(1.0,0.0);
347 g._gamma[1][3]=EvtComplex(-1.0,0.0);
348 g._gamma[2][0]=EvtComplex(-1.0,0.0);
349 g._gamma[3][1]=EvtComplex(1.0,0.0);
350 }
351
352 return g;
353
354}
355
356
357
358
359const EvtGammaMatrix& EvtGammaMatrix::g5(){
360
361 static EvtGammaMatrix g;
362 static int first=1;
363
364 if (first){
365 first = 0;
366 int i,j;
367
368 for(i=0;i<4;i++){
369 for(j=0;j<4;j++){
370 g._gamma[i][j]=EvtComplex(0.0,0.0);
371 }
372 }
373
374 g._gamma[0][2]=EvtComplex(1.0,0.0);
375 g._gamma[1][3]=EvtComplex(1.0,0.0);
376 g._gamma[2][0]=EvtComplex(1.0,0.0);
377 g._gamma[3][1]=EvtComplex(1.0,0.0);
378 }
379
380 return g;
381
382}
383
384
385
386
387const EvtGammaMatrix& EvtGammaMatrix::g(int index) {
388 switch (index) {
389 case 0:
390 return g0();
391 case 1:
392 return g1();
393 case 2:
394 return g2();
395 case 3:
396 return g3();
397 case 5:
398 return g5();
399 default:
400 report(ERROR, "EvtGen") << "Invalid index for four vector: " << index << endl;
401 exit(-2);
402 }
403}
404
405
406
407const EvtGammaMatrix& EvtGammaMatrix::v0(){
408
409 static EvtGammaMatrix g;
410 static int first=1;
411
412 if (first){
413 first = 0;
414 int i,j;
415
416 for(i=0;i<4;i++){
417 for(j=0;j<4;j++){
418 g._gamma[i][j]=EvtComplex(0.0,0.0);
419 }
420 }
421
422 g._gamma[0][0]=EvtComplex(1.0,0.0);
423 g._gamma[1][1]=EvtComplex(1.0,0.0);
424 g._gamma[2][2]=EvtComplex(1.0,0.0);
425 g._gamma[3][3]=EvtComplex(1.0,0.0);
426 }
427
428 return g;
429
430}
431
432
433
434
435
436const EvtGammaMatrix& EvtGammaMatrix::v1(){
437
438 static EvtGammaMatrix g;
439 static int first=1;
440
441 if (first){
442 first = 0;
443 int i,j;
444
445 for(i=0;i<4;i++){
446 for(j=0;j<4;j++){
447 g._gamma[i][j]=EvtComplex(0.0,0.0);
448 }
449 }
450
451 g._gamma[0][3]=EvtComplex(1.0,0.0);
452 g._gamma[1][2]=EvtComplex(1.0,0.0);
453 g._gamma[2][1]=EvtComplex(1.0,0.0);
454 g._gamma[3][0]=EvtComplex(1.0,0.0);
455 }
456
457 return g;
458
459}
460
461
462
463
464const EvtGammaMatrix& EvtGammaMatrix::v2(){
465
466 static EvtGammaMatrix g;
467 static int first=1;
468
469 if (first){
470 first = 0;
471 int i,j;
472
473 for(i=0;i<4;i++){
474 for(j=0;j<4;j++){
475 g._gamma[i][j]=EvtComplex(0.0,0.0);
476 }
477 }
478
479 g._gamma[0][3]=EvtComplex(0.0,-1.0);
480 g._gamma[1][2]=EvtComplex(0.0,1.0);
481 g._gamma[2][1]=EvtComplex(0.0,-1.0);
482 g._gamma[3][0]=EvtComplex(0.0,1.0);
483 }
484
485 return g;
486
487}
488
489
490
491
492const EvtGammaMatrix& EvtGammaMatrix::v3(){
493
494 static EvtGammaMatrix g;
495 static int first=1;
496
497 if (first){
498 first = 0;
499 int i,j;
500
501 for(i=0;i<4;i++){
502 for(j=0;j<4;j++){
503 g._gamma[i][j]=EvtComplex(0.0,0.0);
504 }
505 }
506
507 g._gamma[0][2]=EvtComplex(1.0,0.0);
508 g._gamma[1][3]=EvtComplex(-1.0,0.0);
509 g._gamma[2][0]=EvtComplex(1.0,0.0);
510 g._gamma[3][1]=EvtComplex(-1.0,0.0);
511 }
512
513 return g;
514
515}
516
517
518
519
520
521const EvtGammaMatrix& EvtGammaMatrix::id(){
522
523 static EvtGammaMatrix g;
524 static int first=1;
525
526 if (first){
527 first = 0;
528 int i,j;
529
530 for(i=0;i<4;i++){
531 for(j=0;j<4;j++){
532 g._gamma[i][j]=EvtComplex(0.0,0.0);
533 }
534 }
535
536 g._gamma[0][0]=EvtComplex(1.0,0.0);
537 g._gamma[1][1]=EvtComplex(1.0,0.0);
538 g._gamma[2][2]=EvtComplex(1.0,0.0);
539 g._gamma[3][3]=EvtComplex(1.0,0.0);
540 }
541
542 return g;
543
544}
545
546
547
548
549EvtGammaMatrix& EvtGammaMatrix::operator+=(const EvtGammaMatrix &g){
550
551 int i,j;
552
553 for(i=0;i<4;i++){
554 for(j=0;j<4;j++){
555 _gamma[i][j]+=g._gamma[i][j];
556 }
557 }
558 return *this;
559}
560
561
562
563
564
565EvtGammaMatrix& EvtGammaMatrix::operator-=(const EvtGammaMatrix &g){
566
567 int i,j;
568
569 for(i=0;i<4;i++){
570 for(j=0;j<4;j++){
571 _gamma[i][j]-=g._gamma[i][j];
572 }
573 }
574 return *this;
575}
576
577
578
579EvtGammaMatrix& EvtGammaMatrix::operator*=(const EvtGammaMatrix &g){
580
581 int i,j,k;
582 EvtGammaMatrix temp;
583
584 for(i=0;i<4;i++){
585 for(j=0;j<4;j++){
586 temp._gamma[i][j]=EvtComplex(0.0,0.0);
587 for(k=0;k<4;k++){
588 temp._gamma[i][j]+=_gamma[i][k]*g._gamma[k][j];
589 }
590 }
591 }
592
593 for(i=0;i<4;i++){
594 for(j=0;j<4;j++){
595 _gamma[i][j]=temp._gamma[i][j];
596 }
597 }
598
599 return *this;
600}
601
602
603EvtDiracSpinor operator*(const EvtGammaMatrix& g,const EvtDiracSpinor& d){
604
605 int i,j;
606 EvtDiracSpinor temp;
607
608 for(i=0;i<4;i++){
609 temp.set_spinor(i,EvtComplex(0.0,0.0));
610 for(j=0;j<4;j++){
611 temp.set_spinor(i,temp.get_spinor(i)+g._gamma[i][j]*d.get_spinor(j));
612 }
613 }
614
615 return temp;
616}
617
618/*
619EvtComplex operator*(const EvtDiracSpinor& d,const EvtDiracSpinor& dp){
620
621 int i;
622 EvtComplex temp;
623
624 temp=EvtComplex(0.0,0.0);
625
626 for(i=0;i<4;i++){
627 temp+=::conj(d.get_spinor(i))*dp.get_spinor(i);
628 }
629 return temp;
630}
631*/
632
633// upper index
634const EvtGammaMatrix& EvtGammaMatrix::sigmaUpper(unsigned int mu, unsigned int nu)
635{
636 EvtGammaMatrix a, b;
637 static const EvtTensor4C eta = EvtTensor4C::g(); //metric
638 static EvtGammaMatrix sigma[4][4];
639 static bool hasBeenCalled = false;
640 if (!hasBeenCalled)
641 {
642 EvtComplex I(0, 1);
643 for (int i=0; i<4; ++i)
644 sigma[i][i].init(); // set to 0
645
646 EvtGammaMatrix s01 = I/2 * (g0()*g1() - g1()*g0());
647 EvtGammaMatrix s02 = I/2 * (g0()*g2() - g2()*g0());
648 EvtGammaMatrix s03 = I/2 * (g0()*g3() - g3()*g0());
649 EvtGammaMatrix s12 = I/2 * (g1()*g2() - g2()*g1());
650 EvtGammaMatrix s13 = I/2 * (g1()*g3() - g3()*g1());
651 EvtGammaMatrix s23 = I/2 * (g2()*g3() - g3()*g2());
652 sigma[0][1] = s01;
653 sigma[1][0] = -1*s01;
654 sigma[0][2] = s02;
655 sigma[2][0] = -1*s02;
656 sigma[0][3] = s03;
657 sigma[3][0] = -1*s03;
658 sigma[1][2] = s12;
659 sigma[2][1] = -1*s12;
660 sigma[1][3] = s13;
661 sigma[3][1] = -1*s13;
662 sigma[2][3] = s23;
663 sigma[3][2] = -1*s23;
664 }
665 hasBeenCalled = true;
666
667 if (mu > 3 || nu > 3)
668 {
669 report(ERROR, "EvtSigmaTensor") << "Expected index between 0 and 3, but found " << nu << "!" << endl;
670 assert(0);
671 }
672 return sigma[mu][nu];
673
674}
675
676const EvtGammaMatrix& EvtGammaMatrix::sigmaLower(unsigned int mu, unsigned int nu)
677{
678 const EvtComplex I(0, 1);
679 EvtGammaMatrix a, b;
680 static EvtGammaMatrix sigma[4][4];
681 static bool hasBeenCalled = false;
682 static const EvtTensor4C eta = EvtTensor4C::g();
683
684 if (!hasBeenCalled) // has to be initialized only at the first call
685 {
686 // lower index
687 for (int i=0; i<4; ++i)
688 {
689 a = eta.get(i, 0)*g0() + eta.get(i, 1)*g1() + eta.get(i, 2)*g2() + eta.get(i, 3)*g3();
690 for (int j=0; j<4; ++j)
691 {
692 b = eta.get(j, 0)*g0() + eta.get(j, 1)*g1() + eta.get(j, 2)*g2() + eta.get(j, 3)*g3();
693 sigma[i][j] = I/2 * (a*b - b*a);
694 }
695 }
696 }
697 return sigma[mu][nu];
698}
699
700
701EvtGammaMatrix slash(const EvtVector4C& p)
702{
703 return EvtGammaMatrix::g0()*p.get(0) + EvtGammaMatrix::g1()*p.get(1) + EvtGammaMatrix::g2()*p.get(2) + EvtGammaMatrix::g3()*p.get(3);
704}