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