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) 2000 Caltech, UCSB
11 // Module: EvtHelAmp.cc
13 // Description: Decay model for implementation of generic 2 body
14 // decay specified by the partial wave amplitudes
17 // Modification history:
19 // fkw February 2, 2001 changes to satisfy KCC
20 // RYD September 7, 2000 Module created
22 //------------------------------------------------------------------------
24 #include "EvtGenBase/EvtPatches.hh"
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"
38 EvtEvalHelAmp::~EvtEvalHelAmp() {
46 for(ib=0;ib<_nB;ib++){
53 for(ia=0;ia<_nA;ia++){
58 for(ib=0;ib<_nB;ib++){
63 for(ic=0;ic<_nC;ic++){
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];
87 EvtEvalHelAmp::EvtEvalHelAmp(EvtId idA,
90 EvtComplexPtrPtr HBC){
93 EvtSpinType::spintype typeA=EvtPDL::getSpinType(idA);
94 EvtSpinType::spintype typeB=EvtPDL::getSpinType(idB);
95 EvtSpinType::spintype typeC=EvtPDL::getSpinType(idC);
97 //find out how many states each particle have
98 _nA=EvtSpinType::getSpinStates(typeA);
99 _nB=EvtSpinType::getSpinStates(typeB);
100 _nC=EvtSpinType::getSpinStates(typeC);
102 //find out what 2 times the spin is
103 _JA2=EvtSpinType::getSpin2(typeA);
104 _JB2=EvtSpinType::getSpin2(typeB);
105 _JC2=EvtSpinType::getSpin2(typeC);
109 _lambdaA2=new int[_nA];
110 _lambdaB2=new int[_nB];
111 _lambdaC2=new int[_nC];
113 _HBC=new EvtComplexPtr[_nB];
115 for(ib=0;ib<_nB;ib++){
116 _HBC[ib]=new EvtComplex[_nC];
120 _RA=new EvtComplexPtr[_nA];
121 for(ia=0;ia<_nA;ia++){
122 _RA[ia]=new EvtComplex[_nA];
124 _RB=new EvtComplexPtr[_nB];
125 for(ib=0;ib<_nB;ib++){
126 _RB[ib]=new EvtComplex[_nB];
128 _RC=new EvtComplexPtr[_nC];
129 for(ic=0;ic<_nC;ic++){
130 _RC[ic]=new EvtComplex[_nC];
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];
147 //find the allowed helicities (actually 2*times the helicity!)
149 fillHelicity(_lambdaA2,_nA,_JA2,idA);
150 fillHelicity(_lambdaB2,_nB,_JB2,idB);
151 fillHelicity(_lambdaC2,_nC,_JC2,idC);
153 for(ib=0;ib<_nB;ib++){
154 for(ic=0;ic<_nC;ic++){
155 _HBC[ib][ic]=HBC[ib][ic];
165 double EvtEvalHelAmp::probMax(){
167 double c=1.0/sqrt(4*EvtConst::pi/(_JA2+1));
177 for(itheta=-10;itheta<=10;itheta++){
178 theta=acos(0.099999*itheta);
179 for(ia=0;ia<_nA;ia++){
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]));
195 if (prob>maxprob) maxprob=prob;
205 void EvtEvalHelAmp::evalAmp( EvtParticle *p, EvtAmp& amp){
207 //find theta and phi of the first daughter
209 EvtVector4R pB=p->getDaug(0)->getP4();
211 double theta=acos(pB.get(3)/pB.d3mag());
212 double phi=atan2(pB.get(2),pB.get(1));
214 double c=sqrt((_JA2+1)/(4*EvtConst::pi));
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);
228 _amp[ia][ib][ic]=c*_HBC[ib][ic]*
229 exp(EvtComplex(0.0,phi*0.5*(_lambdaA2[ia]-_lambdaB2[ib]+
230 _lambdaC2[ic])))*dfun;
232 prob1+=real(_amp[ia][ib][ic]*conj(_amp[ia][ib][ic]));
237 setUpRotationMatrices(p,theta,phi);
239 applyRotationMatrices();
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]));
250 amp.vertex(_amp[ia][ib][ic]);
253 amp.vertex(ic,_amp[ia][ib][ic]);
258 amp.vertex(ib,_amp[ia][ib][ic]);
261 amp.vertex(ib,ic,_amp[ia][ib][ic]);
267 amp.vertex(ia,_amp[ia][ib][ic]);
270 amp.vertex(ia,ic,_amp[ia][ib][ic]);
275 amp.vertex(ia,ib,_amp[ia][ib][ic]);
278 amp.vertex(ia,ib,ic,_amp[ia][ib][ic]);
286 if (fabs(prob1-prob2)>0.000001*prob1){
287 report(INFO,"EvtGen") << "prob1,prob2:"<<prob1<<" "<<prob2<<endl;
296 void EvtEvalHelAmp::fillHelicity(int* lambda2,int n,int J2, EvtId id){
300 //photon is special case!
307 //and so is the neutrino!
309 if (EvtPDL::getStdHep(id)>0){
310 //particle i.e. lefthanded
313 //anti particle i.e. righthanded
331 void EvtEvalHelAmp::setUpRotationMatrices(EvtParticle* p,double theta, double phi){
335 case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
339 EvtSpinDensity R=p->rotateToHelicityBasis();
351 _RA[i][j]=R.get(i,j);
360 report(ERROR,"EvtGen") << "Spin2(_JA2)="<<_JA2<<" not supported!"<<endl;
368 case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
374 EvtSpinDensity R=p->getDaug(0)->rotateToHelicityBasis(phi,theta,-phi);
382 _RB[i][j]=conj(R.get(i,j));
391 report(ERROR,"EvtGen") << "Spin2(_JB2)="<<_JB2<<" not supported!"<<endl;
397 case 0: case 1: case 2: case 3: case 4: case 5: case 6: case 7: case 8:
403 EvtSpinDensity R=p->getDaug(1)->rotateToHelicityBasis(phi,EvtConst::pi+theta,phi-EvtConst::pi);
411 _RC[i][j]=conj(R.get(i,j));
420 report(ERROR,"EvtGen") << "Spin2(_JC2)="<<_JC2<<" not supported!"<<endl;
429 void EvtEvalHelAmp::applyRotationMatrices(){
437 for(ia=0;ia<_nA;ia++){
438 for(ib=0;ib<_nB;ib++){
439 for(ic=0;ic<_nC;ic++){
442 temp+=_RC[i][ic]*_amp[ia][ib][i];
444 _amp1[ia][ib][ic]=temp;
451 for(ia=0;ia<_nA;ia++){
452 for(ic=0;ic<_nC;ic++){
453 for(ib=0;ib<_nB;ib++){
456 temp+=_RB[i][ib]*_amp1[ia][i][ic];
458 _amp3[ia][ib][ic]=temp;
465 for(ib=0;ib<_nB;ib++){
466 for(ic=0;ic<_nC;ic++){
467 for(ia=0;ia<_nA;ia++){
470 temp+=_RA[i][ia]*_amp3[i][ib][ic];
472 _amp[ia][ib][ic]=temp;