[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtAmp.cxx

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtAmp.cc
//
// Description: Class to manipulate the amplitudes in the decays.
//
// Modification history:
//
//    RYD     May 29, 1997         Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include <math.h>
#include <assert.h>
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtSpinDensity.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
using std::endl;


EvtAmp::EvtAmp(){
  _ndaug=0;
  _pstates=0;
  _nontrivial=0;
}


EvtAmp::EvtAmp(const EvtAmp& amp){

  int i;

  _ndaug=amp._ndaug;
  _pstates=amp._pstates;
  for(i=0;i<_ndaug;i++){  
    dstates[i]=amp.dstates[i];
    _dnontrivial[i]=amp._dnontrivial[i];
  }
  _nontrivial=amp._nontrivial;

  int namp=1;

  for(i=0;i<_nontrivial;i++){    
    _nstate[i]=amp._nstate[i];
    namp*=_nstate[i];
  }

  for(i=0;i<namp;i++){ 
    assert(i<125);
    _amp[i]=amp._amp[i];
  }
  
}


void EvtAmp::init(EvtId p,int ndaugs,EvtId *daug){
  setNDaug(ndaugs);
  int ichild;
  int daug_states[100],parstates;
  for(ichild=0;ichild<ndaugs;ichild++){

    daug_states[ichild]=
      EvtSpinType::getSpinStates(EvtPDL::getSpinType(daug[ichild]));
    
  }
  
  parstates=EvtSpinType::getSpinStates(EvtPDL::getSpinType(p));

  setNState(parstates,daug_states);

}


void EvtAmp::setNDaug(int n){
  _ndaug=n;
}

void EvtAmp::setNState(int parent_states,int *daug_states){

  _nontrivial=0;
  _pstates=parent_states;
  
  if(_pstates>1) {
     _nstate[_nontrivial]=_pstates;
     _nontrivial++;
  }

  int i;

  for(i=0;i<_ndaug;i++){
    dstates[i]=daug_states[i];
    _dnontrivial[i]=-1;
    if(daug_states[i]>1) {
      _nstate[_nontrivial]=daug_states[i];
      _dnontrivial[i]=_nontrivial;
      _nontrivial++;
    }
  }

  if (_nontrivial>5) {
    report(ERROR,"EvtGen") << "Too many nontrivial states in EvtAmp!"<<endl;
  }

}

void EvtAmp::setAmp(int *ind, const EvtComplex& a){

  int nstatepad = 1;
  int position = ind[0];

  for ( int i=1; i<_nontrivial; i++ ) {
    nstatepad *= _nstate[i-1];
    position += nstatepad*ind[i];
  }
  assert(position<125);
  _amp[position] = a;

}

const EvtComplex& EvtAmp::getAmp(int *ind)const{

  int nstatepad = 1;
  int position = ind[0];

  for ( int i=1; i<_nontrivial; i++ ) {
    nstatepad *= _nstate[i-1];
    position += nstatepad*ind[i];
  }

  return _amp[position];
}


EvtSpinDensity EvtAmp::getSpinDensity(){

  EvtSpinDensity rho;
  rho.setDim(_pstates);

  EvtComplex temp;

  int i,j,n;

  if (_pstates==1) {

    if (_nontrivial==0) {

       rho.set(0,0,_amp[0]*conj(_amp[0]));
       return rho;

    }
    
    n=1;

    temp = EvtComplex(0.0); 

    for(i=0;i<_nontrivial;i++){
      n*=_nstate[i];
    }

    for(i=0;i<n;i++){
      temp+=_amp[i]*conj(_amp[i]);
    }

    rho.set(0,0,temp);;

    return rho;

  }

  else{

    for(i=0;i<_pstates;i++){
      for(j=0;j<_pstates;j++){

        temp = EvtComplex(0.0);

	int kk;

        int allloop = 1;
        for (kk=0;kk<(_nontrivial-1); kk++ ) {
	  allloop *= dstates[kk];
	}
        
        for (kk=0; kk<allloop; kk++) {
	  temp += _amp[_pstates*kk+i]*conj(_amp[_pstates*kk+j]);}

	//        if (_nontrivial>3){
	//report(ERROR,"EvtGen") << "Can't handle so many states in EvtAmp!"<<endl;
	//}
        
        rho.set(i,j,temp);

      }
    }
    return rho; 
  }

} 


EvtSpinDensity EvtAmp::getBackwardSpinDensity(EvtSpinDensity *rho_list){

  EvtSpinDensity rho;

  rho.setDim(_pstates);

  if (_pstates==1){
    rho.set(0,0,EvtComplex(1.0,0.0));
    return rho;
  }

  int k;

  EvtAmp ampprime;

  ampprime=(*this);

  for(k=0;k<_ndaug;k++){
   
    if (dstates[k]!=1){
      ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
    }
  }

  return ampprime.contract(0,(*this));

}


EvtSpinDensity EvtAmp::getForwardSpinDensity(EvtSpinDensity *rho_list,int i){

  EvtSpinDensity rho;

  rho.setDim(dstates[i]);

  int k;

  if (dstates[i]==1) {

    rho.set(0,0,EvtComplex(1.0,0.0));

    return rho;

  }

  EvtAmp ampprime;

  ampprime=(*this);

  if (_pstates!=1){
    ampprime=ampprime.contract(0,rho_list[0]);
  }

  for(k=0;k<i;k++){

    if (dstates[k]!=1){
      ampprime=ampprime.contract(_dnontrivial[k],rho_list[k+1]);
    }
      
  }

  return ampprime.contract(_dnontrivial[i],(*this));

}

EvtAmp EvtAmp::contract(int k,const EvtSpinDensity& rho){

  EvtAmp temp;
  
  int i,j;
  temp._ndaug=_ndaug;
  temp._pstates=_pstates;
  temp._nontrivial=_nontrivial;

  for(i=0;i<_ndaug;i++){
    temp.dstates[i]=dstates[i];
    temp._dnontrivial[i]=_dnontrivial[i];
  }

  if (_nontrivial==0) {
    report(ERROR,"EvtGen")<<"Should not be here EvtAmp!"<<endl;
  }


  for(i=0;i<_nontrivial;i++){
    temp._nstate[i]=_nstate[i];
  }


  EvtComplex c;

  int index[10];
  for (i=0;i<10;i++) {
     index[i] = 0;
  }

  int allloop = 1;
  int indflag,ii;
  for (i=0;i<_nontrivial;i++) {
     allloop *= _nstate[i];
  }

  for ( i=0;i<allloop;i++) {

     c = EvtComplex(0.0);
     int tempint = index[k];
     for (j=0;j<_nstate[k];j++) {
       index[k] = j;
       c+=rho.get(j,tempint)*getAmp(index);
     }
     index[k] = tempint;
       
     temp.setAmp(index,c);

     indflag = 0;
     for ( ii=0;ii<_nontrivial;ii++) {
       if ( indflag == 0 ) {
	 if ( index[ii] == (_nstate[ii]-1) ) {
	   index[ii] = 0;
	 }
	 else {
	   indflag = 1;
	   index[ii] += 1;
	 }
       }
     }

  }
  return temp;

}


EvtSpinDensity EvtAmp::contract(int k,const EvtAmp& amp2){

  int i,j,l;

  EvtComplex temp;
  EvtSpinDensity rho;

  rho.setDim(_nstate[k]);

  int allloop = 1;
  int indflag,ii;
  for (i=0;i<_nontrivial;i++) {
     allloop *= _nstate[i];
  }

  int index[10];
  int index1[10];
  //  int l;
  for(i=0;i<_nstate[k];i++){

    for(j=0;j<_nstate[k];j++){
      if (_nontrivial==0) {
	report(ERROR,"EvtGen")<<"Should not be here1 EvtAmp!"<<endl;
        rho.set(0,0,EvtComplex(1.0,0.0)); 
        return rho;
      }

      for (ii=0;ii<10;ii++) {
	index[ii] = 0;
	index1[ii] = 0;
      }
      index[k] = i;
      index1[k] = j;

      temp = EvtComplex(0.0);

      for ( l=0;l<int(allloop/_nstate[k]);l++) {

	temp+=getAmp(index)*conj(amp2.getAmp(index1));
	indflag = 0;
	for ( ii=0;ii<_nontrivial;ii++) {
          if ( ii!= k) {
	    if ( indflag == 0 ) {
	      if ( index[ii] == (_nstate[ii]-1) ) {
		index[ii] = 0;
		index1[ii] = 0;
	      }
	      else {
		indflag = 1;
		index[ii] += 1;
		index1[ii] += 1;
	      }
	    }
	  }
	}
      }
      rho.set(i,j,temp);
      
    }
  }

  return rho;
}


EvtAmp EvtAmp::contract(int i, const EvtAmp& a1,const EvtAmp& a2){
  
  //Do we need this method?

  assert(a2._pstates>1&&a2._nontrivial==1);
  assert(i<=a1._nontrivial);

  EvtAmp tmp;
  report(DEBUG,"EvtGen") << "EvtAmp::contract not written yet" << endl;
  return tmp;

}


void EvtAmp::dump(){

  int i,list[10];

  report(DEBUG,"EvtGen") << "Number of daugthers:"<<_ndaug<<endl;
  report(DEBUG,"EvtGen") << "Number of states of the parent:"<<_pstates<<endl;
  report(DEBUG,"EvtGen") << "Number of states on daughters:";
  for (i=0;i<_ndaug;i++){
    report(DEBUG,"EvtGen") <<dstates[i]<<" ";
  }
  report(DEBUG,"EvtGen") << endl;
  report(DEBUG,"EvtGen") << "Nontrivial index of  daughters:";
  for (i=0;i<_ndaug;i++){
    report(DEBUG,"EvtGen") <<_dnontrivial[i]<<" ";
  }
  report(DEBUG,"EvtGen") <<endl;
  report(DEBUG,"EvtGen") <<"number of nontrivial states:"<<_nontrivial<<endl;
  report(DEBUG,"EvtGen") << "Nontrivial particles number of states:";
  for (i=0;i<_nontrivial;i++){
    report(DEBUG,"EvtGen") <<_nstate[i]<<" ";
  }
  report(DEBUG,"EvtGen") <<endl;
  report(DEBUG,"EvtGen") <<"Amplitudes:"<<endl;
  if (_nontrivial==0){
    list[0] = 0;
    report(DEBUG,"EvtGen") << getAmp(list) << endl;
  }

  int allloop[10];
  allloop[0]=1;
  for (i=0;i<_nontrivial;i++) {
    if (i==0){
      allloop[i] *= _nstate[i];
    }
    else{
      allloop[i] = allloop[i-1]*_nstate[i];
    }
  }
  int index = 0;
  for (i=0;i<allloop[_nontrivial-1];i++) {
    report(DEBUG,"EvtGen") << getAmp(list) << " ";
    if ( i==allloop[index]-1 ) {
      index ++;
      report(DEBUG,"EvtGen") << endl;
    }
  }

  report(DEBUG,"EvtGen") << "-----------------------------------"<<endl;

}


void EvtAmp::vertex(const EvtComplex& c){
   int list[1];
   list[0] = 0;
   setAmp(list,c);
}

void EvtAmp::vertex(int i,const EvtComplex& c){
   int list[1];
   list[0] = i;
   setAmp(list,c);
}

void EvtAmp::vertex(int i,int j,const EvtComplex& c){
   int list[2];
   list[0] = i;
   list[1] = j;
   setAmp(list,c);
}

void EvtAmp::vertex(int i,int j,int k,const EvtComplex& c){
   int list[3];
   list[0] = i;
   list[1] = j;
   list[2] = k;
   setAmp(list,c);
}

void EvtAmp::vertex(int *i1,const EvtComplex& c){

   setAmp(i1,c);
}


EvtAmp& EvtAmp::operator=(const EvtAmp& amp){

  int i;

  _ndaug=amp._ndaug;
  _pstates=amp._pstates;
  for(i=0;i<_ndaug;i++){  
    dstates[i]=amp.dstates[i];
    _dnontrivial[i]=amp._dnontrivial[i];
  }
  _nontrivial=amp._nontrivial;

  int namp=1;

  for(i=0;i<_nontrivial;i++){    
    _nstate[i]=amp._nstate[i];
    namp*=_nstate[i];
  }

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