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: EvtTensor4C.cc
12 : //
13 : // Description: Implementation of tensor particles.
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/EvtVector4C.hh"
27 : #include "EvtGenBase/EvtTensor4C.hh"
28 : using std::endl;
29 : using std::ostream;
30 :
31 :
32 :
33 0 : EvtTensor4C::EvtTensor4C( const EvtTensor4C& t1 ) {
34 :
35 : int i,j;
36 :
37 0 : for(i=0;i<4;i++) {
38 0 : for(j=0;j<4;j++) {
39 0 : t[i][j] = t1.t[i][j];
40 : }
41 : }
42 :
43 0 : }
44 :
45 0 : EvtTensor4C::~EvtTensor4C() { }
46 :
47 : const EvtTensor4C& EvtTensor4C::g(){
48 :
49 0 : static EvtTensor4C g_metric(1.0,-1.0,-1.0,-1.0);
50 :
51 0 : return g_metric;
52 :
53 0 : }
54 :
55 : EvtTensor4C& EvtTensor4C::operator=(const EvtTensor4C& t1) {
56 : int i,j;
57 :
58 0 : for(i=0;i<4;i++) {
59 0 : for(j=0;j<4;j++) {
60 0 : t[i][j] = t1.t[i][j];
61 : }
62 : }
63 0 : return *this;
64 : }
65 :
66 : EvtTensor4C EvtTensor4C::conj() const {
67 0 : EvtTensor4C temp;
68 :
69 : int i,j;
70 :
71 0 : for(i=0;i<4;i++) {
72 0 : for(j=0;j<4;j++) {
73 0 : temp.set(j,i,::conj(t[i][j]));
74 : }
75 : }
76 : return temp;
77 0 : }
78 :
79 :
80 : EvtTensor4C rotateEuler(const EvtTensor4C& rs,
81 : double alpha,double beta,double gamma){
82 :
83 0 : EvtTensor4C tmp(rs);
84 0 : tmp.applyRotateEuler(alpha,beta,gamma);
85 : return tmp;
86 :
87 0 : }
88 :
89 : EvtTensor4C boostTo(const EvtTensor4C& rs,
90 : const EvtVector4R p4){
91 :
92 0 : EvtTensor4C tmp(rs);
93 0 : tmp.applyBoostTo(p4);
94 : return tmp;
95 :
96 0 : }
97 :
98 : EvtTensor4C boostTo(const EvtTensor4C& rs,
99 : const EvtVector3R boost){
100 :
101 0 : EvtTensor4C tmp(rs);
102 0 : tmp.applyBoostTo(boost);
103 : return tmp;
104 :
105 0 : }
106 :
107 : void EvtTensor4C::applyBoostTo(const EvtVector4R& p4){
108 :
109 0 : double e=p4.get(0);
110 :
111 0 : EvtVector3R boost(p4.get(1)/e,p4.get(2)/e,p4.get(3)/e);
112 :
113 0 : applyBoostTo(boost);
114 :
115 : return;
116 :
117 0 : }
118 :
119 :
120 : void EvtTensor4C::applyBoostTo(const EvtVector3R& boost){
121 :
122 : double bx,by,bz,gamma,b2;
123 0 : double lambda[4][4];
124 0 : EvtComplex tt[4][4];
125 :
126 0 : bx=boost.get(0);
127 0 : by=boost.get(1);
128 0 : bz=boost.get(2);
129 :
130 0 : double bxx=bx*bx;
131 0 : double byy=by*by;
132 0 : double bzz=bz*bz;
133 :
134 0 : b2=bxx+byy+bzz;
135 :
136 :
137 0 : if (b2==0.0){
138 0 : return;
139 : }
140 :
141 0 : assert(b2<1.0);
142 :
143 0 : gamma=1.0/sqrt(1-b2);
144 :
145 :
146 : int i,j,k;
147 :
148 :
149 0 : if (b2==0.0){
150 0 : return ;
151 : }
152 :
153 0 : lambda[0][0]=gamma;
154 0 : lambda[0][1]=gamma*bx;
155 0 : lambda[1][0]=gamma*bx;
156 0 : lambda[0][2]=gamma*by;
157 0 : lambda[2][0]=gamma*by;
158 0 : lambda[0][3]=gamma*bz;
159 0 : lambda[3][0]=gamma*bz;
160 :
161 0 : lambda[1][1]=1.0+(gamma-1.0)*bx*bx/b2;
162 0 : lambda[2][2]=1.0+(gamma-1.0)*by*by/b2;
163 0 : lambda[3][3]=1.0+(gamma-1.0)*bz*bz/b2;
164 :
165 0 : lambda[1][2]=(gamma-1.0)*bx*by/b2;
166 0 : lambda[2][1]=(gamma-1.0)*bx*by/b2;
167 :
168 0 : lambda[1][3]=(gamma-1.0)*bx*bz/b2;
169 0 : lambda[3][1]=(gamma-1.0)*bx*bz/b2;
170 :
171 0 : lambda[3][2]=(gamma-1.0)*bz*by/b2;
172 0 : lambda[2][3]=(gamma-1.0)*bz*by/b2;
173 :
174 0 : for(i=0;i<4;i++){
175 0 : for(j=0;j<4;j++){
176 0 : tt[i][j] = EvtComplex(0.0);
177 0 : for(k=0;k<4;k++){
178 0 : tt[i][j]=tt[i][j]+lambda[j][k]*t[i][k];
179 : }
180 : }
181 : }
182 :
183 0 : for(i=0;i<4;i++){
184 0 : for(j=0;j<4;j++){
185 0 : t[i][j] = EvtComplex(0.0);
186 0 : for(k=0;k<4;k++){
187 0 : t[i][j]=t[i][j]+lambda[i][k]*tt[k][j];
188 : }
189 : }
190 : }
191 :
192 0 : }
193 :
194 : void EvtTensor4C::zero(){
195 : int i,j;
196 0 : for(i=0;i<4;i++){
197 0 : for(j=0;j<4;j++){
198 0 : t[i][j]=EvtComplex(0.0,0.0);
199 : }
200 : }
201 0 : }
202 :
203 :
204 :
205 : ostream& operator<<(ostream& s,const EvtTensor4C& t){
206 :
207 : int i,j;
208 0 : s<< endl;
209 0 : for(i=0;i<4;i++){
210 0 : for(j=0;j<4;j++){
211 0 : s << t.t[i][j];
212 : }
213 0 : s << endl;
214 : }
215 0 : return s;
216 : }
217 :
218 : void EvtTensor4C::setdiag(double g00, double g11, double g22, double g33){
219 0 : t[0][0]=EvtComplex(g00);
220 0 : t[1][1]=EvtComplex(g11);
221 0 : t[2][2]=EvtComplex(g22);
222 0 : t[3][3]=EvtComplex(g33);
223 0 : t[0][1] = EvtComplex(0.0);
224 0 : t[0][2] = EvtComplex(0.0);
225 0 : t[0][3] = EvtComplex(0.0);
226 0 : t[1][0] = EvtComplex(0.0);
227 0 : t[1][2] = EvtComplex(0.0);
228 0 : t[1][3] = EvtComplex(0.0);
229 0 : t[2][0] = EvtComplex(0.0);
230 0 : t[2][1] = EvtComplex(0.0);
231 0 : t[2][3] = EvtComplex(0.0);
232 0 : t[3][0] = EvtComplex(0.0);
233 0 : t[3][1] = EvtComplex(0.0);
234 0 : t[3][2] = EvtComplex(0.0);
235 0 : }
236 :
237 :
238 : EvtTensor4C& EvtTensor4C::operator+=(const EvtTensor4C& t2){
239 :
240 : int i,j;
241 :
242 0 : for (i=0;i<4;i++) {
243 0 : for (j=0;j<4;j++) {
244 0 : t[i][j]+=t2.get(i,j);
245 : }
246 : }
247 0 : return *this;
248 : }
249 :
250 : EvtTensor4C& EvtTensor4C::operator-=(const EvtTensor4C& t2){
251 :
252 : int i,j;
253 :
254 0 : for (i=0;i<4;i++) {
255 0 : for (j=0;j<4;j++) {
256 0 : t[i][j]-=t2.get(i,j);
257 : }
258 : }
259 0 : return *this;
260 : }
261 :
262 :
263 : EvtTensor4C& EvtTensor4C::operator*=(const EvtComplex& c) {
264 : int i,j;
265 :
266 0 : for (i=0;i<4;i++) {
267 0 : for (j=0;j<4;j++) {
268 0 : t[i][j]*=c;
269 : }
270 : }
271 0 : return *this;
272 : }
273 :
274 :
275 : EvtTensor4C operator*(const EvtTensor4C& t1,const EvtComplex& c){
276 :
277 0 : return EvtTensor4C(t1)*=c;
278 :
279 0 : }
280 :
281 : EvtTensor4C operator*(const EvtComplex& c,const EvtTensor4C& t1){
282 :
283 0 : return EvtTensor4C(t1)*=c;
284 :
285 0 : }
286 :
287 :
288 : EvtTensor4C& EvtTensor4C::operator*=(double d) {
289 : int i,j;
290 :
291 0 : for (i=0;i<4;i++) {
292 0 : for (j=0;j<4;j++) {
293 0 : t[i][j]*=EvtComplex(d,0.0);
294 : }
295 : }
296 0 : return *this;
297 : }
298 :
299 :
300 : EvtTensor4C operator*(const EvtTensor4C& t1, double d){
301 :
302 0 : return EvtTensor4C(t1)*=EvtComplex(d,0.0);
303 :
304 0 : }
305 :
306 : EvtTensor4C operator*(double d, const EvtTensor4C& t1){
307 :
308 0 : return EvtTensor4C(t1)*=EvtComplex(d,0.0);
309 :
310 0 : }
311 :
312 : EvtComplex cont(const EvtTensor4C& t1,const EvtTensor4C& t2){
313 :
314 0 : EvtComplex sum(0.0,0.0);
315 : int i,j;
316 :
317 0 : for (i=0;i<4;i++) {
318 0 : for (j=0;j<4;j++) {
319 0 : if ((i==0&&j!=0) || (j==0&&i!=0)) {
320 0 : sum -= t1.t[i][j]*t2.t[i][j];
321 0 : } else {
322 0 : sum += t1.t[i][j]*t2.t[i][j];
323 : }
324 : }
325 : }
326 :
327 : return sum;
328 0 : }
329 :
330 :
331 : EvtTensor4C EvtGenFunctions::directProd(const EvtVector4C& c1,
332 : const EvtVector4C& c2){
333 0 : EvtTensor4C temp;
334 : int i,j;
335 :
336 0 : for (i=0;i<4;i++) {
337 0 : for (j=0;j<4;j++) {
338 0 : temp.set(i,j,c1.get(i)*c2.get(j));
339 : }
340 : }
341 : return temp;
342 0 : }
343 :
344 :
345 : EvtTensor4C EvtGenFunctions::directProd(const EvtVector4C& c1,
346 : const EvtVector4R& c2){
347 0 : EvtTensor4C temp;
348 : int i,j;
349 :
350 0 : for (i=0;i<4;i++) {
351 0 : for (j=0;j<4;j++) {
352 0 : temp.set(i,j,c1.get(i)*c2.get(j));
353 : }
354 : }
355 : return temp;
356 0 : }
357 :
358 :
359 : EvtTensor4C EvtGenFunctions::directProd(const EvtVector4R& c1,
360 : const EvtVector4R& c2){
361 :
362 0 : EvtTensor4C temp;
363 : int i,j;
364 :
365 0 : for (i=0;i<4;i++) {
366 0 : for (j=0;j<4;j++) {
367 0 : temp.t[i][j]=EvtComplex(c1.get(i)*c2.get(j),0.0);
368 : }
369 : }
370 : return temp;
371 0 : }
372 :
373 : EvtTensor4C& EvtTensor4C::addDirProd(const EvtVector4R& p1,const EvtVector4R& p2){
374 :
375 : int i,j;
376 :
377 0 : for (i=0;i<4;i++) {
378 0 : for (j=0;j<4;j++) {
379 0 : t[i][j]+=p1.get(i)*p2.get(j);
380 : }
381 : }
382 0 : return *this;
383 : }
384 :
385 :
386 : EvtTensor4C dual(const EvtTensor4C& t2){
387 :
388 0 : EvtTensor4C temp;
389 :
390 0 : temp.set(0,0,EvtComplex(0.0,0.0));
391 0 : temp.set(1,1,EvtComplex(0.0,0.0));
392 0 : temp.set(2,2,EvtComplex(0.0,0.0));
393 0 : temp.set(3,3,EvtComplex(0.0,0.0));
394 :
395 0 : temp.set(0,1,t2.get(3,2)-t2.get(2,3));
396 0 : temp.set(0,2,-t2.get(3,1)+t2.get(1,3));
397 0 : temp.set(0,3,t2.get(2,1)-t2.get(1,2));
398 :
399 0 : temp.set(1,2,-t2.get(3,0)+t2.get(0,3));
400 0 : temp.set(1,3,t2.get(2,0)-t2.get(0,2));
401 :
402 0 : temp.set(2,3,-t2.get(1,0)+t2.get(0,1));
403 :
404 0 : temp.set(1,0,-temp.get(0,1));
405 0 : temp.set(2,0,-temp.get(0,2));
406 0 : temp.set(3,0,-temp.get(0,3));
407 :
408 0 : temp.set(2,1,-temp.get(1,2));
409 0 : temp.set(3,1,-temp.get(1,3));
410 :
411 0 : temp.set(3,2,-temp.get(2,3));
412 :
413 : return temp;
414 :
415 0 : }
416 :
417 :
418 : EvtTensor4C conj(const EvtTensor4C& t2) {
419 0 : EvtTensor4C temp;
420 :
421 : int i,j;
422 :
423 0 : for(i=0;i<4;i++){
424 0 : for(j=0;j<4;j++){
425 0 : temp.set(i,j,::conj((t2.get(i,j))));
426 : }
427 : }
428 :
429 : return temp;
430 0 : }
431 :
432 :
433 : EvtTensor4C cont22(const EvtTensor4C& t1,const EvtTensor4C& t2){
434 0 : EvtTensor4C temp;
435 :
436 : int i,j;
437 0 : EvtComplex c;
438 :
439 0 : for(i=0;i<4;i++){
440 0 : for(j=0;j<4;j++){
441 0 : c=t1.get(i,0)*t2.get(j,0)-t1.get(i,1)*t2.get(j,1)
442 0 : -t1.get(i,2)*t2.get(j,2)-t1.get(i,3)*t2.get(j,3);
443 0 : temp.set(i,j,c);
444 : }
445 : }
446 :
447 : return temp;
448 0 : }
449 :
450 : EvtTensor4C cont11(const EvtTensor4C& t1,const EvtTensor4C& t2){
451 0 : EvtTensor4C temp;
452 :
453 : int i,j;
454 0 : EvtComplex c;
455 :
456 0 : for(i=0;i<4;i++){
457 0 : for(j=0;j<4;j++){
458 0 : c=t1.get(0,i)*t2.get(0,j)-t1.get(1,i)*t2.get(1,j)
459 0 : -t1.get(2,i)*t2.get(2,j)-t1.get(3,i)*t2.get(3,j);
460 0 : temp.set(i,j,c);
461 : }
462 : }
463 :
464 : return temp;
465 0 : }
466 :
467 :
468 : EvtVector4C EvtTensor4C::cont1(const EvtVector4C& v4) const {
469 0 : EvtVector4C temp;
470 :
471 : int i;
472 :
473 0 : for(i=0;i<4;i++){
474 0 : temp.set(i,t[0][i]*v4.get(0)-t[1][i]*v4.get(1)
475 0 : -t[2][i]*v4.get(2)-t[3][i]*v4.get(3));
476 : }
477 :
478 : return temp;
479 0 : }
480 :
481 : EvtVector4C EvtTensor4C::cont2(const EvtVector4C& v4) const {
482 0 : EvtVector4C temp;
483 :
484 : int i;
485 :
486 0 : for(i=0;i<4;i++){
487 0 : temp.set(i,t[i][0]*v4.get(0)-t[i][1]*v4.get(1)
488 0 : -t[i][2]*v4.get(2)-t[i][3]*v4.get(3));
489 : }
490 :
491 : return temp;
492 0 : }
493 :
494 :
495 : EvtVector4C EvtTensor4C::cont1(const EvtVector4R& v4) const {
496 0 : EvtVector4C temp;
497 :
498 : int i;
499 :
500 0 : for(i=0;i<4;i++){
501 0 : temp.set(i,t[0][i]*v4.get(0)-t[1][i]*v4.get(1)
502 0 : -t[2][i]*v4.get(2)-t[3][i]*v4.get(3));
503 : }
504 :
505 : return temp;
506 0 : }
507 :
508 :
509 : EvtVector4C EvtTensor4C::cont2(const EvtVector4R& v4) const {
510 0 : EvtVector4C temp;
511 :
512 : int i;
513 :
514 0 : for(i=0;i<4;i++){
515 0 : temp.set(i,t[i][0]*v4.get(0)-t[i][1]*v4.get(1)
516 0 : -t[i][2]*v4.get(2)-t[i][3]*v4.get(3));
517 : }
518 :
519 : return temp;
520 0 : }
521 :
522 :
523 :
524 : void EvtTensor4C::applyRotateEuler(double phi,double theta,double ksi){
525 :
526 0 : EvtComplex tt[4][4];
527 : double sp,st,sk,cp,ct,ck;
528 0 : double lambda[4][4];
529 :
530 0 : sp=sin(phi);
531 0 : st=sin(theta);
532 0 : sk=sin(ksi);
533 0 : cp=cos(phi);
534 0 : ct=cos(theta);
535 0 : ck=cos(ksi);
536 :
537 :
538 0 : lambda[0][0]=1.0;
539 0 : lambda[0][1]=0.0;
540 0 : lambda[1][0]=0.0;
541 0 : lambda[0][2]=0.0;
542 0 : lambda[2][0]=0.0;
543 0 : lambda[0][3]=0.0;
544 0 : lambda[3][0]=0.0;
545 :
546 0 : lambda[1][1]= ck*ct*cp-sk*sp;
547 0 : lambda[1][2]=-sk*ct*cp-ck*sp;
548 0 : lambda[1][3]=st*cp;
549 :
550 0 : lambda[2][1]= ck*ct*sp+sk*cp;
551 0 : lambda[2][2]=-sk*ct*sp+ck*cp;
552 0 : lambda[2][3]=st*sp;
553 :
554 0 : lambda[3][1]=-ck*st;
555 0 : lambda[3][2]=sk*st;
556 0 : lambda[3][3]=ct;
557 :
558 :
559 : int i,j,k;
560 :
561 :
562 0 : for(i=0;i<4;i++){
563 0 : for(j=0;j<4;j++){
564 0 : tt[i][j] = EvtComplex(0.0);
565 0 : for(k=0;k<4;k++){
566 0 : tt[i][j]+=lambda[j][k]*t[i][k];
567 : }
568 : }
569 : }
570 :
571 0 : for(i=0;i<4;i++){
572 0 : for(j=0;j<4;j++){
573 0 : t[i][j] = EvtComplex(0.0);
574 0 : for(k=0;k<4;k++){
575 0 : t[i][j]+=lambda[i][k]*tt[k][j];
576 : }
577 : }
578 : }
579 :
580 0 : }
581 :
582 :
583 :
|
