]>
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: EvtCGCoefSingle.cc | |
12 | // | |
13 | // Description: Evaluates Clebsch-Gordon coef for fixed j1 and j2. | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // fkw February 2, 2001 changes to satisfy KCC | |
18 | // RYD August 12, 2000 Module created | |
19 | // | |
20 | //------------------------------------------------------------------------ | |
21 | // | |
22 | #include "EvtGenBase/EvtPatches.hh" | |
23 | #include <stdlib.h> | |
24 | #include <math.h> | |
25 | #include <iostream> | |
26 | #include <assert.h> | |
27 | #include "EvtGenBase/EvtCGCoefSingle.hh" | |
28 | #include "EvtGenBase/EvtOrthogVector.hh" | |
29 | ||
30 | ||
31 | ||
32 | EvtCGCoefSingle::~EvtCGCoefSingle(){ | |
33 | } | |
34 | ||
35 | ||
36 | void EvtCGCoefSingle::init(int j1,int j2){ | |
37 | ||
38 | _j1=j1; | |
39 | _j2=j2; | |
40 | ||
41 | _Jmax=abs(j1+j2); | |
42 | _Jmin=abs(j1-j2); | |
43 | ||
44 | _table.resize((_Jmax-_Jmin)/2+1); | |
45 | ||
46 | int J,M; | |
47 | ||
48 | int lenmax=j1+1; | |
49 | if (j2<j1) lenmax=j2+1; | |
50 | ||
51 | //set vector sizes | |
52 | for(J=_Jmax;J>=_Jmin;J-=2){ | |
53 | _table[(J-_Jmin)/2].resize(J+1); | |
54 | for(M=J;J>=-M;M-=2){ | |
55 | int len=((_j1+_j2)-abs(M))/2+1; | |
56 | if (len>lenmax) len=lenmax; | |
57 | _table[(J-_Jmin)/2][(M+J)/2].resize(len); | |
58 | } | |
59 | } | |
60 | ||
61 | //now fill the vectors | |
62 | for(J=_Jmax;J>=_Jmin;J-=2){ | |
63 | //bootstrap with highest M(=J) as a special case | |
64 | if (J==_Jmax) { | |
65 | cg(J,J,_j1,_j2)=1.0; | |
66 | }else{ | |
67 | int n=(_Jmax-J)/2+1; | |
68 | std::vector<double>* vectors=new std::vector<double>[n-1]; | |
69 | int i,k; | |
70 | for(i=0;i<n-1;i++){ | |
71 | // i corresponds to J=Jmax-2*i | |
72 | vectors[i].resize(n); | |
73 | for(k=0;k<n;k++){ | |
74 | double tmp=_table[(_Jmax-_Jmin)/2-i][(J+_Jmax-2*i)/2][k]; | |
75 | vectors[i][k]=tmp; | |
76 | } | |
77 | } | |
78 | EvtOrthogVector getOrth(n,vectors); | |
79 | std::vector<double> orth=getOrth.getOrthogVector(); | |
80 | int sign=1; | |
81 | if (orth[n-1]<0.0) sign=-1; | |
82 | for(k=0;k<n;k++){ | |
83 | _table[(J-_Jmin)/2][J][k]=sign*orth[k]; | |
84 | } | |
85 | delete [] vectors ; | |
86 | } | |
87 | for(M=J-2;M>=-J;M-=2){ | |
88 | int len=((_j1+_j2)-abs(M))/2+1; | |
89 | if (len>lenmax) len=lenmax; | |
90 | int mmin=M-j2; | |
91 | if (mmin<-j1) mmin=-j1; | |
92 | int m1; | |
93 | for(m1=mmin;m1<mmin+len*2;m1+=2){ | |
94 | int m2=M-m1; | |
95 | double sum=0.0; | |
96 | float fkwTmp = _j1*(_j1+2)-(m1+2)*m1; | |
97 | //fkw 2/2/2001: changes needed to satisfy KCC | |
98 | //fkw if (m1+2<=_j1) sum+=0.5*sqrt(_j1*(_j1+2)-(m1+2)*m1)*cg(J,M+2,m1+2,m2); | |
99 | //fkw if (m2+2<=_j2) sum+=0.5*sqrt(_j2*(_j2+2)-(m2+2)*m2)*cg(J,M+2,m1,m2+2); | |
100 | //fkw sum/=(0.5*sqrt(J*(J+2)-(M+2)*M)); | |
101 | if (m1+2<=_j1) sum+=0.5*sqrt(fkwTmp)*cg(J,M+2,m1+2,m2); | |
102 | fkwTmp = _j2*(_j2+2)-(m2+2)*m2; | |
103 | if (m2+2<=_j2) sum+=0.5*sqrt(fkwTmp)*cg(J,M+2,m1,m2+2); | |
104 | fkwTmp = J*(J+2)-(M+2)*M; | |
105 | sum/=(0.5*sqrt(fkwTmp)); | |
106 | cg(J,M,m1,m2)=sum; | |
107 | } | |
108 | } | |
109 | } | |
110 | ||
111 | ||
112 | ||
113 | } | |
114 | ||
115 | ||
116 | double EvtCGCoefSingle::coef(int J,int M,int j1,int j2,int m1,int m2){ | |
117 | ||
118 | assert(j1==_j1); | |
119 | assert(j2==_j2); | |
120 | ||
121 | return cg(J,M,m1,m2); | |
122 | ||
123 | ||
124 | } | |
125 | ||
126 | ||
127 | double& EvtCGCoefSingle::cg(int J,int M, int m1, int m2){ | |
128 | ||
129 | assert(M==m1+m2); | |
130 | assert(abs(M)<=J); | |
131 | assert(J<=_Jmax); | |
132 | assert(J>=_Jmin); | |
133 | assert(abs(m1)<=_j1); | |
134 | assert(abs(m2)<=_j2); | |
135 | ||
136 | //find lowest m1 allowed for the given M | |
137 | ||
138 | int mmin=M-_j2; | |
139 | ||
140 | if (mmin<-_j1) mmin=-_j1; | |
141 | ||
142 | int n=m1-mmin; | |
143 | ||
144 | return _table[(J-_Jmin)/2][(M+J)/2][n/2]; | |
145 | ||
146 | } | |
147 | ||
148 | ||
149 |