CLHEP/GenericFunctions/LegendreFit.icc

0001 // -*- C++ -*-
0002 // $Id:
0003 #include <sstream>
0004 #include <cmath>
0005 #include <gsl/gsl_sf_legendre.h>
0006 #include <complex>
0007 namespace Genfun {
0008
0009 FUNCTION_OBJECT_IMP(LegendreFit)
0010
0011 inline
0012 LegendreFit::LegendreFit(unsigned int N):
0013 N(N),coefA(N),coefASq(2*N)
0014 {
0015   for (unsigned int i=0;i<N;i++) {
0016     std::ostringstream stream;
0017     stream << "Fraction " << i;
0018     fraction.push_back(new Parameter(stream.str(), 0.5, 0.0, 1.0));
0019   }
0020   for (unsigned int i=0;i<N;i++) {
0021     std::ostringstream stream;
0022     stream << "Phase " << i;
0023     phase.push_back(new Parameter(stream.str(), M_PI, 0.0, 2.0*M_PI));
0024   }
0025 }
0026
0027 inline
0028 LegendreFit::~LegendreFit() {
0029   for (unsigned int i=0;i<N;i++) {
0030     delete fraction[i];
0031     delete phase[i];
0032   }
0033 }
0034
0035 inline
0036 LegendreFit::LegendreFit(const LegendreFit & right):
0037   AbsFunction(),
0038   N(right.N),coefA(right.N), coefASq(2*right.N)
0039 {
0040   for (unsigned int i=0;i<N;i++) {
0041     fraction.push_back(new Parameter(*right.fraction[i]));
0042     phase.push_back(new Parameter(*right.phase[i]));
0043   }
0044 }
0045
0046 inline
0047 double LegendreFit::operator() (double x) const {
0048   recomputeCoefficients();
0049   std::vector<double> Pk(N+1);
0050   gsl_sf_legendre_Pl_array(N, x, &Pk[0]);
0051   unsigned int n=N;
0052   std::complex<double> P=0.0;
0053   std::complex<double> I(0,1.0);
0054   while (1) {
0055     if (n==0) {
0056       double Pn=sqrt((2*n+1.0)/2.0)*Pk[n];
0057
0058       P+=coefA(n)*Pn;
0059       break;
0060     }
0061     else {
0062       double Pn=sqrt((2*n+1.0)/2.0)*Pk[n];
0063       P+=coefA(n)*Pn;
0064       n--;
0065     }
0066   }
0067   return std::norm(P);
0068 }
0069
0070 inline
0071 unsigned int LegendreFit::order() const{
0072   return N;
0073 }
0074 inline
0075 Parameter *LegendreFit::getFraction(unsigned int i) {
0076   return fraction[i];
0077 }
0078 inline
0079 const Parameter *LegendreFit::getFraction(unsigned int i) const{
0080   return fraction[i];
0081 }
0082 inline
0083 Parameter *LegendreFit::getPhase(unsigned int i) {
0084   return phase[i];
0085 }
0086 inline
0087 const Parameter *LegendreFit::getPhase(unsigned int i) const{
0088   return phase[i];
0089 }
0090 inline
0091 const LegendreCoefficientSet & LegendreFit::coefficientsA() const {
0092   recomputeCoefficients();
0093   return coefA;
0094 }
0095
0096 inline
0097 const LegendreCoefficientSet & LegendreFit::coefficientsASq() const {
0098   recomputeCoefficients();
0099   unsigned int LMAX=coefA.getLMax();
0100   for (unsigned int L=0;L<=2*LMAX;L++) {
0101     coefASq(L)=0.0;
0102     for (unsigned int l1=0;l1<=LMAX;l1++) {
0103       for (unsigned int l2=0;l2<=LMAX;l2++) {
0104     if (((l1+l2) >= L) && abs(l1-l2) <= int(L))  {
0105       coefASq(L) += (coefA(l1)*
0106              conj(coefA(l2))*
0107              sqrt((2*l1+1)*(2*l2+1)/4.0)*
0108              pow(ClebschGordan(l1,l2,0,0,L,0),2));
0109     }
0110       }
0111     }
0112   }
0113   return coefASq;
0114 }
0115
0116 inline
0117 void LegendreFit::recomputeCoefficients() const {
0118   unsigned int n=N;
0119   std::complex<double> P=0.0;
0120   std::complex<double> I(0,1.0);
0121   double f=1.0;
0122   while (1) {
0123     if (n==0) {
0124       double fn=1.0;
0125       coefA(n)=sqrt(f*fn);
0126       break;
0127     }
0128     else {
0129       double fn=getFraction(n-1)->getValue();
0130       double px=getPhase(n-1)->getValue();
0131       coefA(n)=exp(I*px)*sqrt(f*fn);
0132       f*=(1-fn);
0133       n--;
0134     }
0135   }
0136 }
0137
0138
0139 } // end namespace Genfun