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: EvtAmp.cc
12 : //
13 : // Description: Class to manipulate the amplitudes in the decays.
14 : //
15 : // Modification history:
16 : //
17 : // RYD May 29, 1997 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/EvtSpinDensity.hh"
27 : #include "EvtGenBase/EvtAmp.hh"
28 : #include "EvtGenBase/EvtReport.hh"
29 : #include "EvtGenBase/EvtId.hh"
30 : #include "EvtGenBase/EvtPDL.hh"
31 : #include "EvtGenBase/EvtParticle.hh"
32 : using std::endl;
33 :
34 :
35 :
36 0 : EvtAmp::EvtAmp(){
37 0 : _ndaug=0;
38 0 : _pstates=0;
39 0 : _nontrivial=0;
40 0 : }
41 :
42 :
43 0 : EvtAmp::EvtAmp(const EvtAmp& amp){
44 :
45 : int i;
46 :
47 0 : _ndaug=amp._ndaug;
48 0 : _pstates=amp._pstates;
49 0 : for(i=0;i<_ndaug;i++){
50 0 : dstates[i]=amp.dstates[i];
51 0 : _dnontrivial[i]=amp._dnontrivial[i];
52 : }
53 0 : _nontrivial=amp._nontrivial;
54 :
55 : int namp=1;
56 :
57 0 : for(i=0;i<_nontrivial;i++){
58 0 : _nstate[i]=amp._nstate[i];
59 0 : namp*=_nstate[i];
60 : }
61 :
62 0 : for(i=0;i<namp;i++){
63 0 : assert(i<125);
64 0 : _amp[i]=amp._amp[i];
65 : }
66 :
67 0 : }
68 :
69 :
70 :
71 : void EvtAmp::init(EvtId p,int ndaugs,EvtId *daug){
72 0 : setNDaug(ndaugs);
73 : int ichild;
74 0 : int daug_states[100],parstates;
75 0 : for(ichild=0;ichild<ndaugs;ichild++){
76 :
77 0 : daug_states[ichild]=
78 0 : EvtSpinType::getSpinStates(EvtPDL::getSpinType(daug[ichild]));
79 :
80 : }
81 :
82 0 : parstates=EvtSpinType::getSpinStates(EvtPDL::getSpinType(p));
83 :
84 0 : setNState(parstates,daug_states);
85 :
86 0 : }
87 :
88 :
89 :
90 :
91 : void EvtAmp::setNDaug(int n){
92 0 : _ndaug=n;
93 0 : }
94 :
95 : void EvtAmp::setNState(int parent_states,int *daug_states){
96 :
97 0 : _nontrivial=0;
98 0 : _pstates=parent_states;
99 :
100 0 : if(_pstates>1) {
101 0 : _nstate[_nontrivial]=_pstates;
102 0 : _nontrivial++;
103 0 : }
104 :
105 : int i;
106 :
107 0 : for(i=0;i<_ndaug;i++){
108 0 : dstates[i]=daug_states[i];
109 0 : _dnontrivial[i]=-1;
110 0 : if(daug_states[i]>1) {
111 0 : _nstate[_nontrivial]=daug_states[i];
112 0 : _dnontrivial[i]=_nontrivial;
113 0 : _nontrivial++;
114 0 : }
115 : }
116 :
117 0 : if (_nontrivial>5) {
118 0 : report(Severity::Error,"EvtGen") << "Too many nontrivial states in EvtAmp!"<<endl;
119 0 : }
120 :
121 0 : }
122 :
123 : void EvtAmp::setAmp(int *ind, const EvtComplex& a){
124 :
125 : int nstatepad = 1;
126 0 : int position = ind[0];
127 :
128 0 : for ( int i=1; i<_nontrivial; i++ ) {
129 0 : nstatepad *= _nstate[i-1];
130 0 : position += nstatepad*ind[i];
131 : }
132 0 : assert(position<125);
133 0 : _amp[position] = a;
134 :
135 0 : }
136 :
137 : const EvtComplex& EvtAmp::getAmp(int *ind)const{
138 :
139 : int nstatepad = 1;
140 0 : int position = ind[0];
141 :
142 0 : for ( int i=1; i<_nontrivial; i++ ) {
143 0 : nstatepad *= _nstate[i-1];
144 0 : position += nstatepad*ind[i];
145 : }
146 :
147 0 : return _amp[position];
148 : }
149 :
150 :
151 : EvtSpinDensity EvtAmp::getSpinDensity(){
152 :
153 0 : EvtSpinDensity rho;
154 0 : rho.setDim(_pstates);
155 :
156 0 : EvtComplex temp;
157 :
158 : int i,j,n;
159 :
160 0 : if (_pstates==1) {
161 :
162 0 : if (_nontrivial==0) {
163 :
164 0 : rho.set(0,0,_amp[0]*conj(_amp[0]));
165 0 : return rho;
166 :
167 : }
168 :
169 : n=1;
170 :
171 0 : temp = EvtComplex(0.0);
172 :
173 0 : for(i=0;i<_nontrivial;i++){
174 0 : n*=_nstate[i];
175 : }
176 :
177 0 : for(i=0;i<n;i++){
178 0 : temp+=_amp[i]*conj(_amp[i]);
179 : }
180 :
181 0 : rho.set(0,0,temp);;
182 :
183 0 : return rho;
184 :
185 : }
186 :
187 : else{
188 :
189 0 : for(i=0;i<_pstates;i++){
190 0 : for(j=0;j<_pstates;j++){
191 :
192 0 : temp = EvtComplex(0.0);
193 :
194 : int kk;
195 :
196 : int allloop = 1;
197 0 : for (kk=0;kk<_ndaug; kk++ ) {
198 0 : allloop *= dstates[kk];
199 : }
200 :
201 0 : for (kk=0; kk<allloop; kk++) {
202 0 : temp += _amp[_pstates*kk+i]*conj(_amp[_pstates*kk+j]);}
203 :
204 : // if (_nontrivial>3){
205 : //report(Severity::Error,"EvtGen") << "Can't handle so many states in EvtAmp!"<<endl;
206 : //}
207 :
208 0 : rho.set(i,j,temp);
209 :
210 : }
211 : }
212 0 : return rho;
213 : }
214 :
215 0 : }
216 :
217 :
218 : EvtSpinDensity EvtAmp::getBackwardSpinDensity(EvtSpinDensity *rho_list){
219 :
220 0 : EvtSpinDensity rho;
221 :
222 0 : rho.setDim(_pstates);
223 :
224 0 : if (_pstates==1){
225 0 : rho.set(0,0,EvtComplex(1.0,0.0));
226 0 : return rho;
227 : }
228 :
229 : int k;
230 :
231 0 : EvtAmp ampprime;
232 :
233 0 : ampprime=(*this);
234 :
235 0 : for(k=0;k<_ndaug;k++){
236 :
237 0 : if (dstates[k]!=1){
238 0 : ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
239 : }
240 : }
241 :
242 0 : return ampprime.contract(0,(*this));
243 :
244 0 : }
245 :
246 :
247 : EvtSpinDensity EvtAmp::getForwardSpinDensity(EvtSpinDensity *rho_list,int i){
248 :
249 0 : EvtSpinDensity rho;
250 :
251 0 : rho.setDim(dstates[i]);
252 :
253 : int k;
254 :
255 0 : if (dstates[i]==1) {
256 :
257 0 : rho.set(0,0,EvtComplex(1.0,0.0));
258 :
259 0 : return rho;
260 :
261 : }
262 :
263 0 : EvtAmp ampprime;
264 :
265 0 : ampprime=(*this);
266 :
267 0 : if (_pstates!=1){
268 0 : ampprime=ampprime.contract(0,rho_list[0]);
269 : }
270 :
271 0 : for(k=0;k<i;k++){
272 :
273 0 : if (dstates[k]!=1){
274 0 : ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
275 : }
276 :
277 : }
278 :
279 0 : return ampprime.contract(_dnontrivial[i],(*this));
280 :
281 0 : }
282 :
283 : EvtAmp EvtAmp::contract(int k,const EvtSpinDensity& rho){
284 :
285 0 : EvtAmp temp;
286 :
287 : int i,j;
288 0 : temp._ndaug=_ndaug;
289 0 : temp._pstates=_pstates;
290 0 : temp._nontrivial=_nontrivial;
291 :
292 0 : for(i=0;i<_ndaug;i++){
293 0 : temp.dstates[i]=dstates[i];
294 0 : temp._dnontrivial[i]=_dnontrivial[i];
295 : }
296 :
297 0 : if (_nontrivial==0) {
298 0 : report(Severity::Error,"EvtGen")<<"Should not be here EvtAmp!"<<endl;
299 0 : }
300 :
301 :
302 0 : for(i=0;i<_nontrivial;i++){
303 0 : temp._nstate[i]=_nstate[i];
304 : }
305 :
306 :
307 0 : EvtComplex c;
308 :
309 0 : int index[10];
310 0 : for (i=0;i<10;i++) {
311 0 : index[i] = 0;
312 : }
313 :
314 : int allloop = 1;
315 : int indflag,ii;
316 0 : for (i=0;i<_nontrivial;i++) {
317 0 : allloop *= _nstate[i];
318 : }
319 :
320 0 : for ( i=0;i<allloop;i++) {
321 :
322 0 : c = EvtComplex(0.0);
323 0 : int tempint = index[k];
324 0 : for (j=0;j<_nstate[k];j++) {
325 0 : index[k] = j;
326 0 : c+=rho.get(j,tempint)*getAmp(index);
327 : }
328 0 : index[k] = tempint;
329 :
330 0 : temp.setAmp(index,c);
331 :
332 : indflag = 0;
333 0 : for ( ii=0;ii<_nontrivial;ii++) {
334 0 : if ( indflag == 0 ) {
335 0 : if ( index[ii] == (_nstate[ii]-1) ) {
336 0 : index[ii] = 0;
337 0 : }
338 : else {
339 : indflag = 1;
340 0 : index[ii] += 1;
341 : }
342 : }
343 : }
344 :
345 : }
346 : return temp;
347 :
348 0 : }
349 :
350 :
351 : EvtSpinDensity EvtAmp::contract(int k,const EvtAmp& amp2){
352 :
353 : int i,j,l;
354 :
355 0 : EvtComplex temp;
356 0 : EvtSpinDensity rho;
357 :
358 0 : rho.setDim(_nstate[k]);
359 :
360 : int allloop = 1;
361 : int indflag,ii;
362 0 : for (i=0;i<_nontrivial;i++) {
363 0 : allloop *= _nstate[i];
364 : }
365 :
366 0 : int index[10];
367 0 : int index1[10];
368 : // int l;
369 0 : for(i=0;i<_nstate[k];i++){
370 :
371 0 : for(j=0;j<_nstate[k];j++){
372 0 : if (_nontrivial==0) {
373 0 : report(Severity::Error,"EvtGen")<<"Should not be here1 EvtAmp!"<<endl;
374 0 : rho.set(0,0,EvtComplex(1.0,0.0));
375 0 : return rho;
376 : }
377 :
378 0 : for (ii=0;ii<10;ii++) {
379 0 : index[ii] = 0;
380 0 : index1[ii] = 0;
381 : }
382 0 : index[k] = i;
383 0 : index1[k] = j;
384 :
385 0 : temp = EvtComplex(0.0);
386 :
387 0 : for ( l=0;l<int(allloop/_nstate[k]);l++) {
388 :
389 0 : temp+=getAmp(index)*conj(amp2.getAmp(index1));
390 : indflag = 0;
391 0 : for ( ii=0;ii<_nontrivial;ii++) {
392 0 : if ( ii!= k) {
393 0 : if ( indflag == 0 ) {
394 0 : if ( index[ii] == (_nstate[ii]-1) ) {
395 0 : index[ii] = 0;
396 0 : index1[ii] = 0;
397 0 : }
398 : else {
399 : indflag = 1;
400 0 : index[ii] += 1;
401 0 : index1[ii] += 1;
402 : }
403 : }
404 : }
405 : }
406 : }
407 0 : rho.set(i,j,temp);
408 :
409 : }
410 : }
411 :
412 0 : return rho;
413 0 : }
414 :
415 :
416 : EvtAmp EvtAmp::contract(int , const EvtAmp& ,const EvtAmp& ){
417 :
418 : //Do we need this method?
419 0 : EvtAmp tmp;
420 0 : report(Severity::Debug,"EvtGen") << "EvtAmp::contract not written yet" << endl;
421 0 : return tmp;
422 :
423 : }
424 :
425 :
426 : void EvtAmp::dump(){
427 :
428 0 : int i,list[10];
429 0 : for (i = 0; i < 10; i++) {list[i] = 0;}
430 :
431 0 : report(Severity::Debug,"EvtGen") << "Number of daugthers:"<<_ndaug<<endl;
432 0 : report(Severity::Debug,"EvtGen") << "Number of states of the parent:"<<_pstates<<endl;
433 0 : report(Severity::Debug,"EvtGen") << "Number of states on daughters:";
434 0 : for (i=0;i<_ndaug;i++){
435 0 : report(Severity::Debug,"EvtGen") <<dstates[i]<<" ";
436 : }
437 0 : report(Severity::Debug,"EvtGen") << endl;
438 0 : report(Severity::Debug,"EvtGen") << "Nontrivial index of daughters:";
439 0 : for (i=0;i<_ndaug;i++){
440 0 : report(Severity::Debug,"EvtGen") <<_dnontrivial[i]<<" ";
441 : }
442 0 : report(Severity::Debug,"EvtGen") <<endl;
443 0 : report(Severity::Debug,"EvtGen") <<"number of nontrivial states:"<<_nontrivial<<endl;
444 0 : report(Severity::Debug,"EvtGen") << "Nontrivial particles number of states:";
445 0 : for (i=0;i<_nontrivial;i++){
446 0 : report(Severity::Debug,"EvtGen") <<_nstate[i]<<" ";
447 : }
448 0 : report(Severity::Debug,"EvtGen") <<endl;
449 0 : report(Severity::Debug,"EvtGen") <<"Amplitudes:"<<endl;
450 0 : if (_nontrivial==0){
451 0 : list[0] = 0;
452 0 : report(Severity::Debug,"EvtGen") << getAmp(list) << endl;
453 0 : }
454 :
455 0 : int allloop[10];
456 0 : for (i = 0; i < 10; i++) {allloop[i] = 0;}
457 :
458 0 : allloop[0]=1;
459 0 : for (i=0;i<_nontrivial;i++) {
460 0 : if (i==0){
461 0 : allloop[i] *= _nstate[i];
462 0 : }
463 : else{
464 0 : allloop[i] = allloop[i-1]*_nstate[i];
465 : }
466 : }
467 : int index = 0;
468 0 : for (i=0;i<allloop[_nontrivial-1];i++) {
469 0 : report(Severity::Debug,"EvtGen") << getAmp(list) << " ";
470 0 : if ( i==allloop[index]-1 ) {
471 0 : index ++;
472 0 : report(Severity::Debug,"EvtGen") << endl;
473 0 : }
474 : }
475 :
476 0 : report(Severity::Debug,"EvtGen") << "-----------------------------------"<<endl;
477 :
478 0 : }
479 :
480 :
481 : void EvtAmp::vertex(const EvtComplex& c){
482 0 : int list[1];
483 0 : list[0] = 0;
484 0 : setAmp(list,c);
485 0 : }
486 :
487 : void EvtAmp::vertex(int i,const EvtComplex& c){
488 0 : int list[1];
489 0 : list[0] = i;
490 0 : setAmp(list,c);
491 0 : }
492 :
493 : void EvtAmp::vertex(int i,int j,const EvtComplex& c){
494 0 : int list[2];
495 0 : list[0] = i;
496 0 : list[1] = j;
497 0 : setAmp(list,c);
498 0 : }
499 :
500 : void EvtAmp::vertex(int i,int j,int k,const EvtComplex& c){
501 0 : int list[3];
502 0 : list[0] = i;
503 0 : list[1] = j;
504 0 : list[2] = k;
505 0 : setAmp(list,c);
506 0 : }
507 :
508 : void EvtAmp::vertex(int *i1,const EvtComplex& c){
509 :
510 0 : setAmp(i1,c);
511 0 : }
512 :
513 :
514 : EvtAmp& EvtAmp::operator=(const EvtAmp& amp){
515 :
516 : int i;
517 :
518 0 : _ndaug=amp._ndaug;
519 0 : _pstates=amp._pstates;
520 0 : for(i=0;i<_ndaug;i++){
521 0 : dstates[i]=amp.dstates[i];
522 0 : _dnontrivial[i]=amp._dnontrivial[i];
523 : }
524 0 : _nontrivial=amp._nontrivial;
525 :
526 : int namp=1;
527 :
528 0 : for(i=0;i<_nontrivial;i++){
529 0 : _nstate[i]=amp._nstate[i];
530 0 : namp*=_nstate[i];
531 : }
532 :
533 0 : for(i=0;i<namp;i++){
534 0 : assert(i<125);
535 0 : _amp[i]=amp._amp[i];
536 : }
537 :
538 0 : return *this;
539 : }
540 :
541 :
542 :
543 :
544 :
545 :
546 :
547 :
548 :
549 :
550 :
|
