1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
13 // Description: Class to manipulate the amplitudes in the decays.
15 // Modification history:
17 // RYD May 29, 1997 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
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"
43 EvtAmp::EvtAmp(const EvtAmp& amp){
48 _pstates=amp._pstates;
49 for(i=0;i<_ndaug;i++){
50 dstates[i]=amp.dstates[i];
51 _dnontrivial[i]=amp._dnontrivial[i];
53 _nontrivial=amp._nontrivial;
57 for(i=0;i<_nontrivial;i++){
58 _nstate[i]=amp._nstate[i];
71 void EvtAmp::init(EvtId p,int ndaugs,EvtId *daug){
74 int daug_states[100],parstates;
75 for(ichild=0;ichild<ndaugs;ichild++){
78 EvtSpinType::getSpinStates(EvtPDL::getSpinType(daug[ichild]));
82 parstates=EvtSpinType::getSpinStates(EvtPDL::getSpinType(p));
84 setNState(parstates,daug_states);
91 void EvtAmp::setNDaug(int n){
95 void EvtAmp::setNState(int parent_states,int *daug_states){
98 _pstates=parent_states;
101 _nstate[_nontrivial]=_pstates;
107 for(i=0;i<_ndaug;i++){
108 dstates[i]=daug_states[i];
110 if(daug_states[i]>1) {
111 _nstate[_nontrivial]=daug_states[i];
112 _dnontrivial[i]=_nontrivial;
118 report(ERROR,"EvtGen") << "Too many nontrivial states in EvtAmp!"<<endl;
123 void EvtAmp::setAmp(int *ind, const EvtComplex& a){
126 int position = ind[0];
128 for ( int i=1; i<_nontrivial; i++ ) {
129 nstatepad *= _nstate[i-1];
130 position += nstatepad*ind[i];
132 assert(position<125);
137 const EvtComplex& EvtAmp::getAmp(int *ind)const{
140 int position = ind[0];
142 for ( int i=1; i<_nontrivial; i++ ) {
143 nstatepad *= _nstate[i-1];
144 position += nstatepad*ind[i];
147 return _amp[position];
151 EvtSpinDensity EvtAmp::getSpinDensity(){
154 rho.setDim(_pstates);
162 if (_nontrivial==0) {
164 rho.set(0,0,_amp[0]*conj(_amp[0]));
171 temp = EvtComplex(0.0);
173 for(i=0;i<_nontrivial;i++){
178 temp+=_amp[i]*conj(_amp[i]);
189 for(i=0;i<_pstates;i++){
190 for(j=0;j<_pstates;j++){
192 temp = EvtComplex(0.0);
197 for (kk=0;kk<_ndaug; kk++ ) {
198 allloop *= dstates[kk];
201 for (kk=0; kk<allloop; kk++) {
202 temp += _amp[_pstates*kk+i]*conj(_amp[_pstates*kk+j]);}
204 // if (_nontrivial>3){
205 //report(ERROR,"EvtGen") << "Can't handle so many states in EvtAmp!"<<endl;
218 EvtSpinDensity EvtAmp::getBackwardSpinDensity(EvtSpinDensity *rho_list){
222 rho.setDim(_pstates);
225 rho.set(0,0,EvtComplex(1.0,0.0));
235 for(k=0;k<_ndaug;k++){
238 ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
242 return ampprime.contract(0,(*this));
247 EvtSpinDensity EvtAmp::getForwardSpinDensity(EvtSpinDensity *rho_list,int i){
251 rho.setDim(dstates[i]);
257 rho.set(0,0,EvtComplex(1.0,0.0));
268 ampprime=ampprime.contract(0,rho_list[0]);
274 ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
279 return ampprime.contract(_dnontrivial[i],(*this));
283 EvtAmp EvtAmp::contract(int k,const EvtSpinDensity& rho){
289 temp._pstates=_pstates;
290 temp._nontrivial=_nontrivial;
292 for(i=0;i<_ndaug;i++){
293 temp.dstates[i]=dstates[i];
294 temp._dnontrivial[i]=_dnontrivial[i];
297 if (_nontrivial==0) {
298 report(ERROR,"EvtGen")<<"Should not be here EvtAmp!"<<endl;
302 for(i=0;i<_nontrivial;i++){
303 temp._nstate[i]=_nstate[i];
316 for (i=0;i<_nontrivial;i++) {
317 allloop *= _nstate[i];
320 for ( i=0;i<allloop;i++) {
323 int tempint = index[k];
324 for (j=0;j<_nstate[k];j++) {
326 c+=rho.get(j,tempint)*getAmp(index);
330 temp.setAmp(index,c);
333 for ( ii=0;ii<_nontrivial;ii++) {
334 if ( indflag == 0 ) {
335 if ( index[ii] == (_nstate[ii]-1) ) {
351 EvtSpinDensity EvtAmp::contract(int k,const EvtAmp& amp2){
358 rho.setDim(_nstate[k]);
362 for (i=0;i<_nontrivial;i++) {
363 allloop *= _nstate[i];
369 for(i=0;i<_nstate[k];i++){
371 for(j=0;j<_nstate[k];j++){
372 if (_nontrivial==0) {
373 report(ERROR,"EvtGen")<<"Should not be here1 EvtAmp!"<<endl;
374 rho.set(0,0,EvtComplex(1.0,0.0));
378 for (ii=0;ii<10;ii++) {
385 temp = EvtComplex(0.0);
387 for ( l=0;l<int(allloop/_nstate[k]);l++) {
389 temp+=getAmp(index)*conj(amp2.getAmp(index1));
391 for ( ii=0;ii<_nontrivial;ii++) {
393 if ( indflag == 0 ) {
394 if ( index[ii] == (_nstate[ii]-1) ) {
416 EvtAmp EvtAmp::contract(int , const EvtAmp& ,const EvtAmp& ){
418 //Do we need this method?
420 report(DEBUG,"EvtGen") << "EvtAmp::contract not written yet" << endl;
429 for (i = 0; i < 10; i++) {list[i] = 0;}
431 report(DEBUG,"EvtGen") << "Number of daugthers:"<<_ndaug<<endl;
432 report(DEBUG,"EvtGen") << "Number of states of the parent:"<<_pstates<<endl;
433 report(DEBUG,"EvtGen") << "Number of states on daughters:";
434 for (i=0;i<_ndaug;i++){
435 report(DEBUG,"EvtGen") <<dstates[i]<<" ";
437 report(DEBUG,"EvtGen") << endl;
438 report(DEBUG,"EvtGen") << "Nontrivial index of daughters:";
439 for (i=0;i<_ndaug;i++){
440 report(DEBUG,"EvtGen") <<_dnontrivial[i]<<" ";
442 report(DEBUG,"EvtGen") <<endl;
443 report(DEBUG,"EvtGen") <<"number of nontrivial states:"<<_nontrivial<<endl;
444 report(DEBUG,"EvtGen") << "Nontrivial particles number of states:";
445 for (i=0;i<_nontrivial;i++){
446 report(DEBUG,"EvtGen") <<_nstate[i]<<" ";
448 report(DEBUG,"EvtGen") <<endl;
449 report(DEBUG,"EvtGen") <<"Amplitudes:"<<endl;
452 report(DEBUG,"EvtGen") << getAmp(list) << endl;
456 for (i = 0; i < 10; i++) {allloop[i] = 0;}
459 for (i=0;i<_nontrivial;i++) {
461 allloop[i] *= _nstate[i];
464 allloop[i] = allloop[i-1]*_nstate[i];
468 for (i=0;i<allloop[_nontrivial-1];i++) {
469 report(DEBUG,"EvtGen") << getAmp(list) << " ";
470 if ( i==allloop[index]-1 ) {
472 report(DEBUG,"EvtGen") << endl;
476 report(DEBUG,"EvtGen") << "-----------------------------------"<<endl;
481 void EvtAmp::vertex(const EvtComplex& c){
487 void EvtAmp::vertex(int i,const EvtComplex& c){
493 void EvtAmp::vertex(int i,int j,const EvtComplex& c){
500 void EvtAmp::vertex(int i,int j,int k,const EvtComplex& c){
508 void EvtAmp::vertex(int *i1,const EvtComplex& c){
514 EvtAmp& EvtAmp::operator=(const EvtAmp& amp){
519 _pstates=amp._pstates;
520 for(i=0;i<_ndaug;i++){
521 dstates[i]=amp.dstates[i];
522 _dnontrivial[i]=amp._dnontrivial[i];
524 _nontrivial=amp._nontrivial;
528 for(i=0;i<_nontrivial;i++){
529 _nstate[i]=amp._nstate[i];