Line data Source code
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 0 : EvtCGCoefSingle::~EvtCGCoefSingle(){
33 0 : }
34 :
35 :
36 : void EvtCGCoefSingle::init(int j1,int j2){
37 :
38 0 : _j1=j1;
39 0 : _j2=j2;
40 :
41 0 : _Jmax=abs(j1+j2);
42 0 : _Jmin=abs(j1-j2);
43 :
44 0 : _table.resize((_Jmax-_Jmin)/2+1);
45 :
46 : int J,M;
47 :
48 0 : int lenmax=j1+1;
49 0 : if (j2<j1) lenmax=j2+1;
50 :
51 : //set vector sizes
52 0 : for(J=_Jmax;J>=_Jmin;J-=2){
53 0 : _table[(J-_Jmin)/2].resize(J+1);
54 0 : for(M=J;J>=-M;M-=2){
55 0 : int len=((_j1+_j2)-abs(M))/2+1;
56 0 : if (len>lenmax) len=lenmax;
57 0 : _table[(J-_Jmin)/2][(M+J)/2].resize(len);
58 : }
59 : }
60 :
61 : //now fill the vectors
62 0 : for(J=_Jmax;J>=_Jmin;J-=2){
63 : //bootstrap with highest M(=J) as a special case
64 0 : if (J==_Jmax) {
65 0 : cg(J,J,_j1,_j2)=1.0;
66 0 : }else{
67 0 : int n=(_Jmax-J)/2+1;
68 0 : std::vector<double>* vectors=new std::vector<double>[n-1];
69 : int i,k;
70 0 : for(i=0;i<n-1;i++){
71 : // i corresponds to J=Jmax-2*i
72 0 : vectors[i].resize(n);
73 0 : for(k=0;k<n;k++){
74 0 : double tmp=_table[(_Jmax-_Jmin)/2-i][(J+_Jmax-2*i)/2][k];
75 0 : vectors[i][k]=tmp;
76 : }
77 : }
78 0 : EvtOrthogVector getOrth(n,vectors);
79 0 : std::vector<double> orth=getOrth.getOrthogVector();
80 : int sign=1;
81 0 : if (orth[n-1]<0.0) sign=-1;
82 0 : for(k=0;k<n;k++){
83 0 : _table[(J-_Jmin)/2][J][k]=sign*orth[k];
84 : }
85 0 : delete [] vectors ;
86 0 : }
87 0 : for(M=J-2;M>=-J;M-=2){
88 0 : int len=((_j1+_j2)-abs(M))/2+1;
89 0 : if (len>lenmax) len=lenmax;
90 0 : int mmin=M-j2;
91 0 : if (mmin<-j1) mmin=-j1;
92 : int m1;
93 0 : for(m1=mmin;m1<mmin+len*2;m1+=2){
94 0 : int m2=M-m1;
95 : double sum=0.0;
96 0 : 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 0 : if (m1+2<=_j1) sum+=0.5*sqrt(fkwTmp)*cg(J,M+2,m1+2,m2);
102 0 : fkwTmp = _j2*(_j2+2)-(m2+2)*m2;
103 0 : if (m2+2<=_j2) sum+=0.5*sqrt(fkwTmp)*cg(J,M+2,m1,m2+2);
104 0 : fkwTmp = J*(J+2)-(M+2)*M;
105 0 : sum/=(0.5*sqrt(fkwTmp));
106 0 : cg(J,M,m1,m2)=sum;
107 : }
108 : }
109 : }
110 :
111 :
112 :
113 0 : }
114 :
115 :
116 : double EvtCGCoefSingle::coef(int J,int M,int j1,int j2,int m1,int m2){
117 :
118 0 : assert(j1==_j1); _unused( j1 );
119 0 : assert(j2==_j2); _unused( j2 );
120 :
121 0 : 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 0 : assert(M==m1+m2); _unused( m2 );
130 0 : assert(abs(M)<=J);
131 0 : assert(J<=_Jmax);
132 0 : assert(J>=_Jmin);
133 0 : assert(abs(m1)<=_j1);
134 0 : assert(abs(m2)<=_j2);
135 :
136 : //find lowest m1 allowed for the given M
137 :
138 0 : int mmin=M-_j2;
139 :
140 0 : if (mmin<-_j1) mmin=-_j1;
141 :
142 0 : int n=m1-mmin;
143 :
144 0 : return _table[(J-_Jmin)/2][(M+J)/2][n/2];
145 :
146 : }
147 :
148 :
149 :
|
