00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "EvtGenBase/EvtPatches.hh"
00025 #include <stdlib.h>
00026 #include "EvtGenBase/EvtParticle.hh"
00027 #include "EvtGenBase/EvtGenKine.hh"
00028 #include "EvtGenBase/EvtPDL.hh"
00029 #include "EvtGenBase/EvtVector4C.hh"
00030 #include "EvtGenBase/EvtTensor4C.hh"
00031 #include "EvtGenBase/EvtReport.hh"
00032 #include "EvtGenModels/EvtPartWave.hh"
00033 #include "EvtGenBase/EvtEvalHelAmp.hh"
00034 #include "EvtGenBase/EvtId.hh"
00035 #include <string>
00036 #include "EvtGenBase/EvtConst.hh"
00037 #include "EvtGenBase/EvtKine.hh"
00038 #include "EvtGenBase/EvtCGCoefSingle.hh"
00039 #include <algorithm>
00040 using std::endl;
00041 EvtPartWave::~EvtPartWave() {}
00042
00043 void EvtPartWave::getName(std::string& model_name){
00044
00045 model_name="PARTWAVE";
00046
00047 }
00048
00049
00050 EvtDecayBase* EvtPartWave::clone(){
00051
00052 return new EvtPartWave;
00053
00054 }
00055
00056 void EvtPartWave::init(){
00057
00058 checkNDaug(2);
00059
00060
00061 int _nA=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getParentId()));
00062 int _nB=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(0)));
00063 int _nC=EvtSpinType::getSpinStates(EvtPDL::getSpinType(getDaug(1)));
00064
00065 if (verbose()){
00066 report(INFO,"EvtGen")<<"_nA,_nB,_nC:"
00067 <<_nA<<","<<_nB<<","<<_nC<<endl;
00068 }
00069
00070
00071 int _JA2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getParentId()));
00072 int _JB2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(0)));
00073 int _JC2=EvtSpinType::getSpin2(EvtPDL::getSpinType(getDaug(1)));
00074
00075 if (verbose()){
00076 report(INFO,"EvtGen")<<"_JA2,_JB2,_JC2:"
00077 <<_JA2<<","<<_JB2<<","<<_JC2<<endl;
00078 }
00079
00080
00081
00082 int* _lambdaA2=new int[_nA];
00083 int* _lambdaB2=new int[_nB];
00084 int* _lambdaC2=new int[_nC];
00085
00086 EvtComplexPtr* _HBC=new EvtComplexPtr[_nB];
00087 int ib,ic;
00088 for(ib=0;ib<_nB;ib++){
00089 _HBC[ib]=new EvtComplex[_nC];
00090 }
00091
00092
00093 int i;
00094
00095
00096 fillHelicity(_lambdaA2,_nA,_JA2);
00097 fillHelicity(_lambdaB2,_nB,_JB2);
00098 fillHelicity(_lambdaC2,_nC,_JC2);
00099
00100 if (verbose()){
00101 report(INFO,"EvtGen")<<"Helicity states of particle A:"<<endl;
00102 for(i=0;i<_nA;i++){
00103 report(INFO,"EvtGen")<<_lambdaA2[i]<<endl;
00104 }
00105
00106 report(INFO,"EvtGen")<<"Helicity states of particle B:"<<endl;
00107 for(i=0;i<_nB;i++){
00108 report(INFO,"EvtGen")<<_lambdaB2[i]<<endl;
00109 }
00110
00111 report(INFO,"EvtGen")<<"Helicity states of particle C:"<<endl;
00112 for(i=0;i<_nC;i++){
00113 report(INFO,"EvtGen")<<_lambdaC2[i]<<endl;
00114 }
00115
00116 report(INFO,"EvtGen")<<"Will now figure out the valid (M_LS) states:"<<endl;
00117
00118 }
00119
00120 int Lmin=std::max(_JA2-_JB2-_JC2,std::max(_JB2-_JA2-_JC2,_JC2-_JA2-_JB2));
00121 if (Lmin<0) Lmin=0;
00122
00123 int Lmax=_JA2+_JB2+_JC2;
00124
00125 int L;
00126
00127 int _nPartialWaveAmp=0;
00128
00129 int _nL[50];
00130 int _nS[50];
00131
00132 for (L=Lmin;L<=Lmax;L+=2){
00133 int Smin=abs(L-_JA2);
00134 if (Smin<abs(_JB2-_JC2)) Smin=abs(_JB2-_JC2);
00135 int Smax=L+_JA2;
00136 if (Smax>abs(_JB2+_JC2)) Smax=abs(_JB2+_JC2);
00137 int S;
00138 for (S=Smin;S<=Smax;S+=2){
00139 _nL[_nPartialWaveAmp]=L;
00140 _nS[_nPartialWaveAmp]=S;
00141
00142 _nPartialWaveAmp++;
00143 if (verbose()){
00144 report(INFO,"EvtGen")<<"M["<<L<<"]["<<S<<"]"<<endl;
00145 }
00146 }
00147 }
00148
00149 checkNArg(_nPartialWaveAmp*2);
00150
00151 int argcounter=0;
00152
00153 EvtComplex _M[50];
00154
00155 for(i=0;i<_nPartialWaveAmp;i++){
00156 _M[i]=getArg(argcounter)*exp(EvtComplex(0.0,getArg(argcounter+1)));;
00157 argcounter+=2;
00158 if (verbose()){
00159 report(INFO,"EvtGen")<<"M["<<_nL[i]<<"]["<<_nS[i]<<"]="<<_M[i]<<endl;
00160 }
00161 }
00162
00163
00164
00165
00166
00167 for(ib=0;ib<_nB;ib++){
00168 for(ic=0;ic<_nC;ic++){
00169 _HBC[ib][ic]=0.0;
00170 if (abs(_lambdaB2[ib]-_lambdaC2[ic])<=_JA2){
00171 for(i=0;i<_nPartialWaveAmp;i++){
00172 int L=_nL[i];
00173 int S=_nS[i];
00174 int lambda2=_lambdaB2[ib];
00175 int lambda3=_lambdaC2[ic];
00176 int s1=_JA2;
00177 int s2=_JB2;
00178 int s3=_JC2;
00179 int m1=lambda2-lambda3;
00180 EvtCGCoefSingle c1(s2,s3);
00181 EvtCGCoefSingle c2(L,S);
00182
00183 if (verbose()){
00184 report(INFO,"EvtGen") << "s2,lambda2:"<<s2<<" "<<lambda2<<endl;
00185 }
00186
00187 double fkwTmp = (L+1.0)/(s1+1.0);
00188
00189 if (S>=abs(m1)){
00190
00191 EvtComplex tmp=sqrt(fkwTmp)
00192 *c1.coef(S,m1,s2,s3,lambda2,-lambda3)
00193 *c2.coef(s1,m1,L,S,0,m1)*_M[i];
00194
00195 _HBC[ib][ic]+=tmp;
00196 }
00197 }
00198 if (verbose()){
00199 report(INFO,"EvtGen")<<"_HBC["<<ib<<"]["<<ic<<"]="<<_HBC[ib][ic]<<endl;
00200 }
00201 }
00202 }
00203 }
00204
00205 _evalHelAmp=new EvtEvalHelAmp(EvtPDL::getSpinType(getParentId()),
00206 EvtPDL::getSpinType(getDaug(0)),
00207 EvtPDL::getSpinType(getDaug(1)),
00208 _HBC);
00209
00210
00211 }
00212
00213
00214 void EvtPartWave::initProbMax(){
00215
00216 double maxprob=_evalHelAmp->probMax();
00217
00218 if (verbose()){
00219 report(INFO,"EvtGen")<<"Calculated probmax"<<maxprob<<endl;
00220 }
00221
00222 setProbMax(maxprob);
00223
00224 }
00225
00226
00227 void EvtPartWave::decay( EvtParticle *p){
00228
00229
00230 p->initializePhaseSpace(getNDaug(),getDaugs());
00231
00232 _evalHelAmp->evalAmp(p,_amp2);
00233
00234 return;
00235
00236 }
00237
00238
00239
00240 void EvtPartWave::fillHelicity(int* lambda2,int n,int J2){
00241
00242 int i;
00243
00244
00245 if (n==2&&J2==2) {
00246 lambda2[0]=2;
00247 lambda2[1]=-2;
00248 return;
00249 }
00250
00251 assert(n==J2+1);
00252
00253 for(i=0;i<n;i++){
00254 lambda2[i]=n-i*2-1;
00255 }
00256
00257 return;
00258
00259 }
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270