00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <stdlib.h>
00024 #include <math.h>
00025 #include <iostream>
00026 #include <assert.h>
00027 #include "EvtGenBase/EvtCGCoefSingle.hh"
00028 #include "EvtGenBase/EvtOrthogVector.hh"
00029
00030
00031
00032 EvtCGCoefSingle::~EvtCGCoefSingle(){
00033 }
00034
00035
00036 void EvtCGCoefSingle::init(int j1,int j2){
00037
00038 _j1=j1;
00039 _j2=j2;
00040
00041 _Jmax=abs(j1+j2);
00042 _Jmin=abs(j1-j2);
00043
00044 _table.resize((_Jmax-_Jmin)/2+1);
00045
00046 int J,M;
00047
00048 int lenmax=j1+1;
00049 if (j2<j1) lenmax=j2+1;
00050
00051
00052 for(J=_Jmax;J>=_Jmin;J-=2){
00053 _table[(J-_Jmin)/2].resize(J+1);
00054 for(M=J;J>=-M;M-=2){
00055 int len=((_j1+_j2)-abs(M))/2+1;
00056 if (len>lenmax) len=lenmax;
00057 _table[(J-_Jmin)/2][(M+J)/2].resize(len);
00058 }
00059 }
00060
00061
00062 for(J=_Jmax;J>=_Jmin;J-=2){
00063
00064 if (J==_Jmax) {
00065 cg(J,J,_j1,_j2)=1.0;
00066 }else{
00067 int n=(_Jmax-J)/2+1;
00068 std::vector<double>* vectors=new std::vector<double>[n-1];
00069 int i,k;
00070 for(i=0;i<n-1;i++){
00071
00072 vectors[i].resize(n);
00073 for(k=0;k<n;k++){
00074 double tmp=_table[(_Jmax-_Jmin)/2-i][(J+_Jmax-2*i)/2][k];
00075 vectors[i][k]=tmp;
00076 }
00077 }
00078 EvtOrthogVector getOrth(n,vectors);
00079 std::vector<double> orth=getOrth.getOrthogVector();
00080 int sign=1;
00081 if (orth[n-1]<0.0) sign=-1;
00082 for(k=0;k<n;k++){
00083 _table[(J-_Jmin)/2][J][k]=sign*orth[k];
00084 }
00085
00086 }
00087 for(M=J-2;M>=-J;M-=2){
00088 int len=((_j1+_j2)-abs(M))/2+1;
00089 if (len>lenmax) len=lenmax;
00090 int mmin=M-j2;
00091 if (mmin<-j1) mmin=-j1;
00092 int m1;
00093 for(m1=mmin;m1<mmin+len*2;m1+=2){
00094 int m2=M-m1;
00095 double sum=0.0;
00096 float fkwTmp = _j1*(_j1+2)-(m1+2)*m1;
00097
00098
00099
00100
00101 if (m1+2<=_j1) sum+=0.5*sqrt(fkwTmp)*cg(J,M+2,m1+2,m2);
00102 fkwTmp = _j2*(_j2+2)-(m2+2)*m2;
00103 if (m2+2<=_j2) sum+=0.5*sqrt(fkwTmp)*cg(J,M+2,m1,m2+2);
00104 fkwTmp = J*(J+2)-(M+2)*M;
00105 sum/=(0.5*sqrt(fkwTmp));
00106 cg(J,M,m1,m2)=sum;
00107 }
00108 }
00109 }
00110
00111
00112
00113 }
00114
00115
00116 double EvtCGCoefSingle::coef(int J,int M,int j1,int j2,int m1,int m2){
00117
00118 assert(j1==_j1);
00119 assert(j2==_j2);
00120
00121 return cg(J,M,m1,m2);
00122
00123
00124 }
00125
00126
00127 double& EvtCGCoefSingle::cg(int J,int M, int m1, int m2){
00128
00129 assert(M==m1+m2);
00130 assert(abs(M)<=J);
00131 assert(J<=_Jmax);
00132 assert(J>=_Jmin);
00133 assert(abs(m1)<=_j1);
00134 assert(abs(m2)<=_j2);
00135
00136
00137
00138 int mmin=M-_j2;
00139
00140 if (mmin<-_j1) mmin=-_j1;
00141
00142 int n=m1-mmin;
00143
00144 return _table[(J-_Jmin)/2][(M+J)/2][n/2];
00145
00146 }
00147
00148
00149