TEvtGen/EvtGen/EvtGenBase/EvtCGCoefSingle.cpp

   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); _unused( j1 );
 119   assert(j2==_j2); _unused( 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); _unused( 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