[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtEvalHelAmp.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) 2000      Caltech, UCSB
//
// Module: EvtHelAmp.cc
//
// Description: Decay model for implementation of generic 2 body
//              decay specified by the partial wave amplitudes
//
//
// Modification history:
//
//    fkw        February 2, 2001     changes to satisfy KCC
//    RYD       September 7, 2000       Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtEvalHelAmp.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtdFunction.hh"
#include "EvtGenBase/EvtAmp.hh"
using std::endl;


EvtEvalHelAmp::~EvtEvalHelAmp() {

  //deallocate memory
  delete [] _lambdaA2;
  delete [] _lambdaB2;
  delete [] _lambdaC2;

  int ia,ib,ic;
  for(ib=0;ib<_nB;ib++){
    delete [] _HBC[ib];
  }

  delete [] _HBC;


  for(ia=0;ia<_nA;ia++){
    delete [] _RA[ia];
  }
  delete [] _RA;

  for(ib=0;ib<_nB;ib++){
    delete [] _RB[ib];
  }
  delete [] _RB;

  for(ic=0;ic<_nC;ic++){
    delete [] _RC[ic];
  }
  delete [] _RC;

 
  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      delete [] _amp[ia][ib];
      delete [] _amp1[ia][ib];
      delete [] _amp3[ia][ib];
    }
    delete [] _amp[ia];
    delete [] _amp1[ia];
    delete [] _amp3[ia];
  }

  delete [] _amp;
  delete [] _amp1;
  delete [] _amp3;

}


EvtEvalHelAmp::EvtEvalHelAmp(EvtId idA,
			     EvtId idB,
			     EvtId idC,
			     EvtComplexPtrPtr HBC){


  EvtSpinType::spintype typeA=EvtPDL::getSpinType(idA);
  EvtSpinType::spintype typeB=EvtPDL::getSpinType(idB);
  EvtSpinType::spintype typeC=EvtPDL::getSpinType(idC);

  //find out how many states each particle have
  _nA=EvtSpinType::getSpinStates(typeA);
  _nB=EvtSpinType::getSpinStates(typeB);
  _nC=EvtSpinType::getSpinStates(typeC);

  //find out what 2 times the spin is
  _JA2=EvtSpinType::getSpin2(typeA);
  _JB2=EvtSpinType::getSpin2(typeB);
  _JC2=EvtSpinType::getSpin2(typeC);


  //allocate memory
  _lambdaA2=new int[_nA];
  _lambdaB2=new int[_nB];
  _lambdaC2=new int[_nC];

  _HBC=new EvtComplexPtr[_nB];
  int ia,ib,ic;
  for(ib=0;ib<_nB;ib++){
    _HBC[ib]=new EvtComplex[_nC];
  }


  _RA=new EvtComplexPtr[_nA];
  for(ia=0;ia<_nA;ia++){
    _RA[ia]=new EvtComplex[_nA];
  }
  _RB=new EvtComplexPtr[_nB];
  for(ib=0;ib<_nB;ib++){
    _RB[ib]=new EvtComplex[_nB];
  }
  _RC=new EvtComplexPtr[_nC];
  for(ic=0;ic<_nC;ic++){
    _RC[ic]=new EvtComplex[_nC];
  }
  
  _amp=new EvtComplexPtrPtr[_nA];
  _amp1=new EvtComplexPtrPtr[_nA];
  _amp3=new EvtComplexPtrPtr[_nA];
  for(ia=0;ia<_nA;ia++){
    _amp[ia]=new EvtComplexPtr[_nB];
    _amp1[ia]=new EvtComplexPtr[_nB];
    _amp3[ia]=new EvtComplexPtr[_nB];
    for(ib=0;ib<_nB;ib++){
      _amp[ia][ib]=new EvtComplex[_nC];
      _amp1[ia][ib]=new EvtComplex[_nC];
      _amp3[ia][ib]=new EvtComplex[_nC];
    }
  }

  //find the allowed helicities (actually 2*times the helicity!)

  fillHelicity(_lambdaA2,_nA,_JA2,idA);
  fillHelicity(_lambdaB2,_nB,_JB2,idB);
  fillHelicity(_lambdaC2,_nC,_JC2,idC);

  for(ib=0;ib<_nB;ib++){
    for(ic=0;ic<_nC;ic++){
      _HBC[ib][ic]=HBC[ib][ic];
    }
  }
}


double EvtEvalHelAmp::probMax(){

  double c=1.0/sqrt(4*EvtConst::pi/(_JA2+1));

  int ia,ib,ic;


  double theta;
  int itheta;

  double maxprob=0.0;

  for(itheta=-10;itheta<=10;itheta++){
    theta=acos(0.099999*itheta);
    for(ia=0;ia<_nA;ia++){
      double prob=0.0;
      for(ib=0;ib<_nB;ib++){
	for(ic=0;ic<_nC;ic++){
	  _amp[ia][ib][ic]=0.0;
	  if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
	    _amp[ia][ib][ic]=c*_HBC[ib][ic]*
	      EvtdFunction::d(_JA2,_lambdaA2[ia],
			      _lambdaB2[ib]-_lambdaC2[ic],theta);
	    prob+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
	  }
	}
      }
      
      prob*=sqrt(1.0*_nA);
      
      if (prob>maxprob) maxprob=prob;

    }
  }

  return maxprob;

}


void EvtEvalHelAmp::evalAmp( EvtParticle *p, EvtAmp& amp){

  //find theta and phi of the first daughter
  
  EvtVector4R pB=p->getDaug(0)->getP4();

  double theta=acos(pB.get(3)/pB.d3mag());
  double phi=atan2(pB.get(2),pB.get(1));

  double c=sqrt((_JA2+1)/(4*EvtConst::pi));

  int ia,ib,ic;

  double prob1=0.0;

  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      for(ic=0;ic<_nC;ic++){
	_amp[ia][ib][ic]=0.0;
	if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
	  double dfun=EvtdFunction::d(_JA2,_lambdaA2[ia],
				      _lambdaB2[ib]-_lambdaC2[ic],theta);

	  _amp[ia][ib][ic]=c*_HBC[ib][ic]*
	    exp(EvtComplex(0.0,phi*0.5*(_lambdaA2[ia]-_lambdaB2[ib]+
	  				_lambdaC2[ic])))*dfun;
	}
	prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
      }
    }
  }

  setUpRotationMatrices(p,theta,phi);

  applyRotationMatrices();

  double prob2=0.0;

  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      for(ic=0;ic<_nC;ic++){
	prob2+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
	if (_nA==1){
	  if (_nB==1){
	    if (_nC==1){
	      amp.vertex(_amp[ia][ib][ic]);
	    }
	    else{
	      amp.vertex(ic,_amp[ia][ib][ic]);
	    }
	  }
	  else{
	    if (_nC==1){
	      amp.vertex(ib,_amp[ia][ib][ic]);
	    }
	    else{
	      amp.vertex(ib,ic,_amp[ia][ib][ic]);
	    }
	  }
	}else{
	  if (_nB==1){
	    if (_nC==1){
	      amp.vertex(ia,_amp[ia][ib][ic]);
	    }
	    else{
	      amp.vertex(ia,ic,_amp[ia][ib][ic]);
	    }
	  }
	  else{
	    if (_nC==1){
	      amp.vertex(ia,ib,_amp[ia][ib][ic]);
	    }
	    else{
	      amp.vertex(ia,ib,ic,_amp[ia][ib][ic]);
	    }
	  }
	}
      }
    }
  }

  if (fabs(prob1-prob2)>0.000001*prob1){
    report(INFO,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl;
    ::abort();
  }
    
  return ;

}


void EvtEvalHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){
  
  int i;
  
  //photon is special case!
  if (n==2&&J2==2) {
    lambda2[0]=2;
    lambda2[1]=-2;
    return;
  }
  
  //and so is the neutrino!
  if (n==1&&J2==1) {
    if (EvtPDL::getStdHep(id)>0){
	//particle i.e. lefthanded
        lambda2[0]=-1;
    }else{
	//anti particle i.e. righthanded
        lambda2[0]=1;
    }
    return;
  }


  assert(n==J2+1);

  for(i=0;i<n;i++){
    lambda2[i]=n-i*2-1;
  }

  return;

}


void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){

  switch(_JA2){

  case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:

    {

      EvtSpinDensity R=p->rotateToHelicityBasis();

      
      int i,j,n;
      
      n=R.getDim();
      
      assert(n==_nA);
	
      
      for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	  _RA[i][j]=R.get(i,j);
	}
      }

    }

    break;

  default:
    report(ERROR,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl;
    ::abort();
  }
  
  
  switch(_JB2){


  case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:

    {
      
      int i,j,n;

      EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi);
      
      n=R.getDim();
      
      assert(n==_nB);
	
      for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	  _RB[i][j]=conj(R.get(i,j));
	}
      }

    }

    break;

  default:
    report(ERROR,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl;
    ::abort();
  }
  
  switch(_JC2){

  case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:

    {

      int i,j,n;

      EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi);
            
      n=R.getDim();

      assert(n==_nC);

      for(i=0;i<n;i++){
	for(j=0;j<n;j++){
	  _RC[i][j]=conj(R.get(i,j));
	}
      }

    }

    break;

  default:
    report(ERROR,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl;
    ::abort();
  }
  
  
}


void EvtEvalHelAmp::applyRotationMatrices(){

  int ia,ib,ic,i;
  
  EvtComplex temp;


  for(ia=0;ia<_nA;ia++){
    for(ib=0;ib<_nB;ib++){
      for(ic=0;ic<_nC;ic++){
	temp=0;
	for(i=0;i<_nC;i++){
	  temp+=_RC[i][ic]*_amp[ia][ib][i];
	}
	_amp1[ia][ib][ic]=temp;
      }
    }
  }


  for(ia=0;ia<_nA;ia++){
    for(ic=0;ic<_nC;ic++){
      for(ib=0;ib<_nB;ib++){
  	temp=0;
  	for(i=0;i<_nB;i++){
  	  temp+=_RB[i][ib]*_amp1[ia][i][ic];
  	}
  	_amp3[ia][ib][ic]=temp;
      }
    }
  }
  

  for(ib=0;ib<_nB;ib++){
    for(ic=0;ic<_nC;ic++){
      for(ia=0;ia<_nA;ia++){
	temp=0;
	for(i=0;i<_nA;i++){
	  temp+=_RA[i][ia]*_amp3[i][ib][ic];
	}
	_amp[ia][ib][ic]=temp;
      }
    }
  }


}
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) 2000 Caltech, UCSB
	10	//
	11	// Module: EvtHelAmp.cc
	12	//
	13	// Description: Decay model for implementation of generic 2 body
	14	// decay specified by the partial wave amplitudes
	15	//
	16	//
	17	// Modification history:
	18	//
	19	// fkw February 2, 2001 changes to satisfy KCC
	20	// RYD September 7, 2000 Module created
	21	//
	22	//------------------------------------------------------------------------
	23	//
	24	#include "EvtGenBase/EvtPatches.hh"
	25	#include <stdlib.h>
	26	#include "EvtGenBase/EvtParticle.hh"
	27	#include "EvtGenBase/EvtPDL.hh"
	28	#include "EvtGenBase/EvtVector4C.hh"
	29	#include "EvtGenBase/EvtReport.hh"
	30	#include "EvtGenBase/EvtEvalHelAmp.hh"
	31	#include "EvtGenBase/EvtId.hh"
	32	#include "EvtGenBase/EvtConst.hh"
	33	#include "EvtGenBase/EvtdFunction.hh"
	34	#include "EvtGenBase/EvtAmp.hh"
	35	using std::endl;
	36
	37
	38	EvtEvalHelAmp::~EvtEvalHelAmp() {
	39
	40	//deallocate memory
	41	delete [] _lambdaA2;
	42	delete [] _lambdaB2;
	43	delete [] _lambdaC2;
	44
	45	int ia,ib,ic;
	46	for(ib=0;ib<_nB;ib++){
	47	delete [] _HBC[ib];
	48	}
	49
	50	delete [] _HBC;
	51
	52
	53	for(ia=0;ia<_nA;ia++){
	54	delete [] _RA[ia];
	55	}
	56	delete [] _RA;
	57
	58	for(ib=0;ib<_nB;ib++){
	59	delete [] _RB[ib];
	60	}
	61	delete [] _RB;
	62
	63	for(ic=0;ic<_nC;ic++){
	64	delete [] _RC[ic];
65	}
66	delete [] _RC;
67
68
69	for(ia=0;ia<_nA;ia++){
70	for(ib=0;ib<_nB;ib++){
71	delete [] _amp[ia][ib];
72	delete [] _amp1[ia][ib];
73	delete [] _amp3[ia][ib];
74	}
75	delete [] _amp[ia];
76	delete [] _amp1[ia];
77	delete [] _amp3[ia];
78	}
79
80	delete [] _amp;
81	delete [] _amp1;
82	delete [] _amp3;
83
84	}
85
86
87	EvtEvalHelAmp::EvtEvalHelAmp(EvtId idA,
88	EvtId idB,
89	EvtId idC,
90	EvtComplexPtrPtr HBC){
91
92
93	EvtSpinType::spintype typeA=EvtPDL::getSpinType(idA);
94	EvtSpinType::spintype typeB=EvtPDL::getSpinType(idB);
95	EvtSpinType::spintype typeC=EvtPDL::getSpinType(idC);
96
97	//find out how many states each particle have
98	_nA=EvtSpinType::getSpinStates(typeA);
99	_nB=EvtSpinType::getSpinStates(typeB);
100	_nC=EvtSpinType::getSpinStates(typeC);
101
102	//find out what 2 times the spin is
103	_JA2=EvtSpinType::getSpin2(typeA);
104	_JB2=EvtSpinType::getSpin2(typeB);
105	_JC2=EvtSpinType::getSpin2(typeC);
106
107
108	//allocate memory
109	_lambdaA2=new int[_nA];
110	_lambdaB2=new int[_nB];
111	_lambdaC2=new int[_nC];
112
113	_HBC=new EvtComplexPtr[_nB];
114	int ia,ib,ic;
115	for(ib=0;ib<_nB;ib++){
116	_HBC[ib]=new EvtComplex[_nC];
117	}
118
119
120	_RA=new EvtComplexPtr[_nA];
121	for(ia=0;ia<_nA;ia++){
122	_RA[ia]=new EvtComplex[_nA];
123	}
124	_RB=new EvtComplexPtr[_nB];
125	for(ib=0;ib<_nB;ib++){
126	_RB[ib]=new EvtComplex[_nB];
127	}
128	_RC=new EvtComplexPtr[_nC];
129	for(ic=0;ic<_nC;ic++){
130	_RC[ic]=new EvtComplex[_nC];
131	}
132
133	_amp=new EvtComplexPtrPtr[_nA];
134	_amp1=new EvtComplexPtrPtr[_nA];
135	_amp3=new EvtComplexPtrPtr[_nA];
136	for(ia=0;ia<_nA;ia++){
137	_amp[ia]=new EvtComplexPtr[_nB];
138	_amp1[ia]=new EvtComplexPtr[_nB];
139	_amp3[ia]=new EvtComplexPtr[_nB];
140	for(ib=0;ib<_nB;ib++){
141	_amp[ia][ib]=new EvtComplex[_nC];
142	_amp1[ia][ib]=new EvtComplex[_nC];
143	_amp3[ia][ib]=new EvtComplex[_nC];
144	}
145	}
146
147	//find the allowed helicities (actually 2*times the helicity!)
148
149	fillHelicity(_lambdaA2,_nA,_JA2,idA);
150	fillHelicity(_lambdaB2,_nB,_JB2,idB);
151	fillHelicity(_lambdaC2,_nC,_JC2,idC);
152
153	for(ib=0;ib<_nB;ib++){
154	for(ic=0;ic<_nC;ic++){
155	_HBC[ib][ic]=HBC[ib][ic];
156	}
157	}
158	}
159
160
161
162
163
164
165	double EvtEvalHelAmp::probMax(){
166
167	double c=1.0/sqrt(4*EvtConst::pi/(_JA2+1));
168
169	int ia,ib,ic;
170
171
172	double theta;
173	int itheta;
174
175	double maxprob=0.0;
176
177	for(itheta=-10;itheta<=10;itheta++){
178	theta=acos(0.099999*itheta);
179	for(ia=0;ia<_nA;ia++){
180	double prob=0.0;
181	for(ib=0;ib<_nB;ib++){
182	for(ic=0;ic<_nC;ic++){
183	_amp[ia][ib][ic]=0.0;
184	if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
185	_amp[ia][ib][ic]=c_HBC[ib][ic]
186	EvtdFunction::d(_JA2,_lambdaA2[ia],
187	_lambdaB2[ib]-_lambdaC2[ic],theta);
188	prob+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
189	}
190	}
191	}
192
193	prob=sqrt(1.0_nA);
194
195	if (prob>maxprob) maxprob=prob;
196
197	}
198	}
199
200	return maxprob;
201
202	}
203
204
205	void EvtEvalHelAmp::evalAmp( EvtParticle *p, EvtAmp& amp){
206
207	//find theta and phi of the first daughter
208
209	EvtVector4R pB=p->getDaug(0)->getP4();
210
211	double theta=acos(pB.get(3)/pB.d3mag());
212	double phi=atan2(pB.get(2),pB.get(1));
213
214	double c=sqrt((_JA2+1)/(4*EvtConst::pi));
215
216	int ia,ib,ic;
217
218	double prob1=0.0;
219
220	for(ia=0;ia<_nA;ia++){
221	for(ib=0;ib<_nB;ib++){
222	for(ic=0;ic<_nC;ic++){
223	_amp[ia][ib][ic]=0.0;
224	if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2) {
225	double dfun=EvtdFunction::d(_JA2,_lambdaA2[ia],
226	_lambdaB2[ib]-_lambdaC2[ic],theta);
227
228	_amp[ia][ib][ic]=c_HBC[ib][ic]
229	exp(EvtComplex(0.0,phi0.5(_lambdaA2[ia]-_lambdaB2[ib]+
230	_lambdaC2[ic])))*dfun;
231	}
232	prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
233	}
234	}
235	}
236
237	setUpRotationMatrices(p,theta,phi);
238
239	applyRotationMatrices();
240
241	double prob2=0.0;
242
243	for(ia=0;ia<_nA;ia++){
244	for(ib=0;ib<_nB;ib++){
245	for(ic=0;ic<_nC;ic++){
246	prob2+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
247	if (_nA==1){
248	if (_nB==1){
249	if (_nC==1){
250	amp.vertex(_amp[ia][ib][ic]);
251	}
252	else{
253	amp.vertex(ic,_amp[ia][ib][ic]);
254	}
255	}
256	else{
257	if (_nC==1){
258	amp.vertex(ib,_amp[ia][ib][ic]);
259	}
260	else{
261	amp.vertex(ib,ic,_amp[ia][ib][ic]);
262	}
263	}
264	}else{
265	if (_nB==1){
266	if (_nC==1){
267	amp.vertex(ia,_amp[ia][ib][ic]);
268	}
269	else{
270	amp.vertex(ia,ic,_amp[ia][ib][ic]);
271	}
272	}
273	else{
274	if (_nC==1){
275	amp.vertex(ia,ib,_amp[ia][ib][ic]);
276	}
277	else{
278	amp.vertex(ia,ib,ic,_amp[ia][ib][ic]);
279	}
280	}
281	}
282	}
283	}
284	}
285
286	if (fabs(prob1-prob2)>0.000001*prob1){
287	report(INFO,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl;
288	::abort();
289	}
290
291	return ;
292
293	}
294
295
296	void EvtEvalHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){
297
298	int i;
299
300	//photon is special case!
301	if (n==2&&J2==2) {
302	lambda2[0]=2;
303	lambda2[1]=-2;
304	return;
305	}
306
307	//and so is the neutrino!
308	if (n==1&&J2==1) {
309	if (EvtPDL::getStdHep(id)>0){
310	//particle i.e. lefthanded
311	lambda2[0]=-1;
312	}else{
313	//anti particle i.e. righthanded
314	lambda2[0]=1;
315	}
316	return;
317	}
318
319
320	assert(n==J2+1);
321
322	for(i=0;i<n;i++){
323	lambda2[i]=n-i*2-1;
324	}
325
326	return;
327
328	}
329
330
331	void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){
332
333	switch(_JA2){
334
335	case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
336
337	{
338
339	EvtSpinDensity R=p->rotateToHelicityBasis();
340
341
342	int i,j,n;
343
344	n=R.getDim();
345
346	assert(n==_nA);
347
348
349	for(i=0;i<n;i++){
350	for(j=0;j<n;j++){
351	_RA[i][j]=R.get(i,j);
352	}
353	}
354
355	}
356
357	break;
358
359	default:
360	report(ERROR,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl;
361	::abort();
362	}
363
364
365	switch(_JB2){
366
367
368	case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
369
370	{
371
372	int i,j,n;
373
374	EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi);
375
376	n=R.getDim();
377
378	assert(n==_nB);
379
380	for(i=0;i<n;i++){
381	for(j=0;j<n;j++){
382	_RB[i][j]=conj(R.get(i,j));
383	}
384	}
385
386	}
387
388	break;
389
390	default:
391	report(ERROR,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl;
392	::abort();
393	}
394
395	switch(_JC2){
396
397	case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
398
399	{
400
401	int i,j,n;
402
403	EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi);
404
405	n=R.getDim();
406
407	assert(n==_nC);
408
409	for(i=0;i<n;i++){
410	for(j=0;j<n;j++){
411	_RC[i][j]=conj(R.get(i,j));
412	}
413	}
414
415	}
416
417	break;
418
419	default:
420	report(ERROR,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl;
421	::abort();
422	}
423
424
425
426	}
427
428
429	void EvtEvalHelAmp::applyRotationMatrices(){
430
431	int ia,ib,ic,i;
432
433	EvtComplex temp;
434
435
436
437	for(ia=0;ia<_nA;ia++){
438	for(ib=0;ib<_nB;ib++){
439	for(ic=0;ic<_nC;ic++){
440	temp=0;
441	for(i=0;i<_nC;i++){
442	temp+=_RC[i][ic]*_amp[ia][ib][i];
443	}
444	_amp1[ia][ib][ic]=temp;
445	}
446	}
447	}
448
449
450
451	for(ia=0;ia<_nA;ia++){
452	for(ic=0;ic<_nC;ic++){
453	for(ib=0;ib<_nB;ib++){
454	temp=0;
455	for(i=0;i<_nB;i++){
456	temp+=_RB[i][ib]*_amp1[ia][i][ic];
457	}
458	_amp3[ia][ib][ic]=temp;
459	}
460	}
461	}
462
463
464
465	for(ib=0;ib<_nB;ib++){
466	for(ic=0;ic<_nC;ic++){
467	for(ia=0;ia<_nA;ia++){
468	temp=0;
469	for(i=0;i<_nA;i++){
470	temp+=_RA[i][ia]*_amp3[i][ib][ic];
471	}
472	_amp[ia][ib][ic]=temp;
473	}
474	}
475	}
476
477
478	}
479
480
481
482
483
484
485
486
487
488
489
490