2d84f9f6fa60a1c3dadbcef3bee96faaf039ba2b
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtGammaMatrix.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: 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>
32 using std::endl;
33 using std::ostream;
34
35 EvtGammaMatrix::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
47 EvtGammaMatrix operator*(const EvtGammaMatrix& g, const EvtComplex& c)
48 {
49     return c*g;
50 }
51
52
53 EvtGammaMatrix 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
69 ostream& 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
83 EvtGammaMatrix::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
93 EvtGammaMatrix::~EvtGammaMatrix() {}
94
95 EvtGammaMatrix& 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
106 void 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
118 const 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
148 const 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
179 const 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
211 const 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
244 const 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
274 const 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
302 const 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
331 const 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
359 const 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
387 const 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
407 const 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
436 const 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
464 const 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
492 const 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
521 const 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
549 EvtGammaMatrix& 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
565 EvtGammaMatrix& 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
579 EvtGammaMatrix& 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
603 EvtDiracSpinor 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 // upper index
619 const EvtGammaMatrix& EvtGammaMatrix::sigmaUpper(unsigned int mu, unsigned int nu)
620 {
621     EvtGammaMatrix a, b;
622     static const EvtTensor4C eta = EvtTensor4C::g(); //metric
623     static EvtGammaMatrix sigma[4][4];
624     static bool hasBeenCalled = false;
625     if (!hasBeenCalled)
626     {
627         EvtComplex I(0, 1);
628         for (int i=0; i<4; ++i)
629             sigma[i][i].init(); // set to 0
630         
631         EvtGammaMatrix s01 = I/2 * (g0()*g1() - g1()*g0());
632         EvtGammaMatrix s02 = I/2 * (g0()*g2() - g2()*g0());
633         EvtGammaMatrix s03 = I/2 * (g0()*g3() - g3()*g0());
634         EvtGammaMatrix s12 = I/2 * (g1()*g2() - g2()*g1());
635         EvtGammaMatrix s13 = I/2 * (g1()*g3() - g3()*g1());
636         EvtGammaMatrix s23 = I/2 * (g2()*g3() - g3()*g2());
637         sigma[0][1] = s01;
638         sigma[1][0] = -1*s01;
639         sigma[0][2] = s02;
640         sigma[2][0] = -1*s02;
641         sigma[0][3] = s03;
642         sigma[3][0] = -1*s03;
643         sigma[1][2] = s12;
644         sigma[2][1] = -1*s12;
645         sigma[1][3] = s13;
646         sigma[3][1] = -1*s13;
647         sigma[2][3] = s23;
648         sigma[3][2] = -1*s23;
649     }
650     hasBeenCalled = true;
651         
652     if (mu > 3 || nu > 3)
653     {
654         report(ERROR, "EvtSigmaTensor") << "Expected index between 0 and 3, but found " << nu << "!" << endl;
655         assert(0);
656     }
657     return sigma[mu][nu];
658     
659 }
660
661 const EvtGammaMatrix& EvtGammaMatrix::sigmaLower(unsigned int mu, unsigned int nu)
662 {
663     const EvtComplex I(0, 1);
664     EvtGammaMatrix a, b;
665     static EvtGammaMatrix sigma[4][4];
666     static bool hasBeenCalled = false;
667     static const EvtTensor4C eta = EvtTensor4C::g();
668     
669     if (!hasBeenCalled) // has to be initialized only at the first call
670     {
671         // lower index
672         for (int i=0; i<4; ++i)
673         {
674             a = eta.get(i, 0)*g0() + eta.get(i, 1)*g1() + eta.get(i, 2)*g2() + eta.get(i, 3)*g3();
675             for (int j=0; j<4; ++j)
676             {
677                 b = eta.get(j, 0)*g0() + eta.get(j, 1)*g1() + eta.get(j, 2)*g2() + eta.get(j, 3)*g3();
678                 sigma[i][j] = I/2 * (a*b - b*a);
679             }
680         }
681     }
682     return sigma[mu][nu];    
683 }
684
685
686 EvtGammaMatrix EvtGenFunctions::slash(const EvtVector4C& p)
687 {
688     return EvtGammaMatrix::g0()*p.get(0) + 
689       EvtGammaMatrix::g1()*p.get(1) + 
690       EvtGammaMatrix::g2()*p.get(2) + 
691       EvtGammaMatrix::g3()*p.get(3);
692 }
693
694 EvtGammaMatrix EvtGenFunctions::slash(const EvtVector4R& p)
695 {
696   return EvtGammaMatrix::g0()*p.get(0) + 
697     EvtGammaMatrix::g1()*p.get(1) + 
698     EvtGammaMatrix::g2()*p.get(2) + 
699     EvtGammaMatrix::g3()*p.get(3);
700 }