CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
LegendreFit.icc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id:
3#include <sstream>
4#include <cmath>
5#include <gsl/gsl_sf_legendre.h>
6#include <complex>
7namespace Genfun {
8
9FUNCTION_OBJECT_IMP(LegendreFit)
10
11inline
13N(N),coefA(N),coefASq(2*N)
14{
15 for (unsigned int i=0;i<N;i++) {
16 std::ostringstream stream;
17 stream << "Fraction " << i;
18 fraction.push_back(new Parameter(stream.str(), 0.5, 0.0, 1.0));
19 }
20 for (unsigned int i=0;i<N;i++) {
21 std::ostringstream stream;
22 stream << "Phase " << i;
23 phase.push_back(new Parameter(stream.str(), M_PI, 0.0, 2.0*M_PI));
24 }
25}
26
27inline
29 for (unsigned int i=0;i<N;i++) {
30 delete fraction[i];
31 delete phase[i];
32 }
33}
34
35inline
38 N(right.N),coefA(right.N), coefASq(2*right.N)
39{
40 for (unsigned int i=0;i<N;i++) {
41 fraction.push_back(new Parameter(*right.fraction[i]));
42 phase.push_back(new Parameter(*right.phase[i]));
43 }
44}
45
46inline
47double LegendreFit::operator() (double x) const {
49 std::vector<double> Pk(N+1);
50 gsl_sf_legendre_Pl_array(N, x, &Pk[0]);
51 unsigned int n=N;
52 std::complex<double> P=0.0;
53 std::complex<double> I(0,1.0);
54 while (1) {
55 if (n==0) {
56 double Pn=sqrt((2*n+1.0)/2.0)*Pk[n];
57
58 P+=coefA(n)*Pn;
59 break;
60 }
61 else {
62 double Pn=sqrt((2*n+1.0)/2.0)*Pk[n];
63 P+=coefA(n)*Pn;
64 n--;
65 }
66 }
67 return std::norm(P);
68}
69
70inline
71unsigned int LegendreFit::order() const{
72 return N;
73}
74inline
76 return fraction[i];
77}
78inline
79const Parameter *LegendreFit::getFraction(unsigned int i) const{
80 return fraction[i];
81}
82inline
84 return phase[i];
85}
86inline
87const Parameter *LegendreFit::getPhase(unsigned int i) const{
88 return phase[i];
89}
90inline
93 return coefA;
94}
95
96inline
99 unsigned int LMAX=coefA.getLMax();
100 for (unsigned int L=0;L<=2*LMAX;L++) {
101 coefASq(L)=0.0;
102 for (unsigned int l1=0;l1<=LMAX;l1++) {
103 for (unsigned int l2=0;l2<=LMAX;l2++) {
104 if (((l1+l2) >= L) && abs(l1-l2) <= int(L)) {
105 coefASq(L) += (coefA(l1)*
106 conj(coefA(l2))*
107 sqrt((2*l1+1)*(2*l2+1)/4.0)*
108 pow(ClebschGordan(l1,l2,0,0,L,0),2));
109 }
110 }
111 }
112 }
113 return coefASq;
114}
115
116inline
118 unsigned int n=N;
119 std::complex<double> P=0.0;
120 std::complex<double> I(0,1.0);
121 double f=1.0;
122 while (1) {
123 if (n==0) {
124 double fn=1.0;
125 coefA(n)=sqrt(f*fn);
126 break;
127 }
128 else {
129 double fn=getFraction(n-1)->getValue();
130 double px=getPhase(n-1)->getValue();
131 coefA(n)=exp(I*px)*sqrt(f*fn);
132 f*=(1-fn);
133 n--;
134 }
135 }
136}
137
138
139} // end namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
Parameter * getFraction(unsigned int i)
LegendreFit(unsigned int N)
Parameter * getPhase(unsigned int i)
void recomputeCoefficients() const
unsigned int order() const
const LegendreCoefficientSet & coefficientsASq() const
virtual double operator()(double argument) const override
const LegendreCoefficientSet & coefficientsA() const
virtual double getValue() const
Definition Abs.hh:14