Commit 9b9f4b10 authored by Joseph Boudreau's avatar Joseph Boudreau
Browse files

Adding a round of work which includes multiparameter fitting functions:...

Adding a round of work which includes multiparameter fitting functions: legendre, sphericalHarmonic, and Fourier. These are normalized distributions built from orthogonal polynomials.
parent 22d73541
......@@ -154,15 +154,15 @@ private:
//
#define FUNCTION_OBJECT_IMP(classname) \
FunctionComposition classname::operator()(const AbsFunction & function) const\
inline FunctionComposition classname::operator()(const AbsFunction & function) const\
{ \
return AbsFunction::operator() (function); \
} \
ParameterComposition classname::operator()(const AbsParameter & p) const\
inline ParameterComposition classname::operator()(const AbsParameter & p) const\
{ \
return AbsFunction::operator() (p); \
} \
classname *classname::clone() const \
inline classname *classname::clone() const \
{ \
return new classname(*this); \
}
......
// -*- C++ -*-
// $Id:
//-------------------------------------------------------------//
// //
// This functional returns the Likelihood of a functoin //
// given some data //
// //
//-------------------------------------------------------------//
#ifndef _EfficiencyFunctional_h_
#define _EfficiencyFunctional_h_
#include "CLHEP/GenericFunctions/AbsFunctional.hh"
#include "CLHEP/GenericFunctions/ArgumentList.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class EfficiencyFunctional:public AbsFunctional {
public:
// Constructor:
EfficiencyFunctional(const ArgumentList & aList);
// Destructor:
~EfficiencyFunctional();
// Evaluate ChiSquared of a function w.r.t the data.
virtual double operator [] (const AbsFunction & function) const;
private:
const ArgumentList _aList;
};
} // namespace Genfun
#endif
// -*- C++ -*-
// $Id:
//---------------------FourierFit-------------------------------------------//
// //
// Class FourierFit. This is a fitting function consisting of a super //
// position of N legendre polynomials. Cascading fractions and phases are //
// the input parameters. Function is normalized to one (on [0,2PI]) //
// Joe Boudreau, Petar Maksimovic, January 2000 //
// //
//--------------------------------------------------------------------------//
#ifndef FourierFit_h
#define FourierFit_h 1
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class FourierFit : public AbsFunction {
FUNCTION_OBJECT_DEF(FourierFit)
public:
// Constructor
FourierFit(unsigned int N);
// Copy constructor
FourierFit(const FourierFit &right);
// Destructor
virtual ~FourierFit();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
unsigned int order() const;
Parameter *getFraction(unsigned int i);
const Parameter *getFraction(unsigned int i) const;
Parameter *getPhase(unsigned int i);
const Parameter *getPhase(unsigned int i) const;
private:
// It is illegal to assign an adjustable constant
const FourierFit & operator=(const FourierFit &right);
//
const unsigned int N;
std::vector <Genfun::Parameter *> fraction;
std::vector <Genfun::Parameter *> phase;
};
} // namespace Genfun
#include "CLHEP/GenericFunctions/FourierFit.icc"
#endif
// -*- C++ -*-
// $Id:
#include <sstream>
#include <cmath>
#include <complex>
namespace Genfun {
FUNCTION_OBJECT_IMP(FourierFit)
inline
FourierFit::FourierFit(unsigned int N):
N(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
FourierFit::~FourierFit() {
for (unsigned int i=0;i<N;i++) {
delete fraction[i];
delete phase[i];
}
}
inline
FourierFit::FourierFit(const FourierFit & right):
N(right.N)
{
for (int i=0;i<N;i++) {
fraction.push_back(new Parameter(*right.fraction[i]));
phase.push_back(new Parameter(*right.phase[i]));
}
}
inline
double FourierFit::operator() (double x) 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;
double Pn=sqrt(1/2.0/M_PI);
P+=(sqrt(f*fn)*Pn);
break;
}
else {
double fn=getFraction(n-1)->getValue();
double px=getPhase(n-1)->getValue();
double Pn=sqrt(1/M_PI)*sin(n*x/2.0);
P+=exp(I*px)*sqrt(f*fn)*Pn;
f*=(1-fn);
n--;
}
}
return std::norm(P);
}
inline
unsigned int FourierFit::order() const{
return N;
}
inline
Parameter *FourierFit::getFraction(unsigned int i) {
return fraction[i];
}
inline
const Parameter *FourierFit::getFraction(unsigned int i) const{
return fraction[i];
}
inline
Parameter *FourierFit::getPhase(unsigned int i) {
return phase[i];
}
inline
const Parameter *FourierFit::getPhase(unsigned int i) const{
return phase[i];
}
} // end namespace Genfun
// -*- C++ -*-
// $Id: LegendreFit.hh,v 1.2 2003/09/06 14:04:13 boudreau Exp $
//---------------------LegendreFit------------------------------------------//
// //
// Class LegendreFit. This is a fitting function consisting of a super //
// position of N legendre polynomials. Cascading fractions and phases are //
// the input parameters. Function is normalized to one (on [-1,1]) //
// Joe Boudreau, Petar Maksimovic, January 2000 //
// //
//--------------------------------------------------------------------------//
#ifndef LegendreFit_h
#define LegendreFit_h 1
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class LegendreFit : public AbsFunction {
FUNCTION_OBJECT_DEF(LegendreFit)
public:
// Constructor
LegendreFit(unsigned int N);
// Copy constructor
LegendreFit(const LegendreFit &right);
// Destructor
virtual ~LegendreFit();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
unsigned int order() const;
Parameter *getFraction(unsigned int i);
const Parameter *getFraction(unsigned int i) const;
Parameter *getPhase(unsigned int i);
const Parameter *getPhase(unsigned int i) const;
private:
// It is illegal to assign an adjustable constant
const LegendreFit & operator=(const LegendreFit &right);
//
const unsigned int N;
std::vector <Genfun::Parameter *> fraction;
std::vector <Genfun::Parameter *> phase;
};
} // namespace Genfun
#include "CLHEP/GenericFunctions/LegendreFit.icc"
#endif
// -*- 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)
{
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):
N(right.N)
{
for (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 {
double Pk[N+1];
gsl_sf_legendre_Pl_array(N, x, Pk);
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;
double Pn=sqrt((2*n+1.0))*Pk[n];
P+=(sqrt(f*fn)*Pn);
break;
}
else {
double fn=getFraction(n-1)->getValue();
double px=getPhase(n-1)->getValue();
double Pn=sqrt((2*n+1.0))*Pk[n];
P+=exp(I*px)*sqrt(f*fn)*Pn;
f*=(1-fn);
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];
}
} // end namespace Genfun
......@@ -44,6 +44,7 @@ pkginclude_HEADERS = \
DefiniteIntegral.hh \
DoubleParamToArgAdaptor.hh \
DoubleParamToArgAdaptor.icc \
EfficiencyFunctional.hh \
EllipticIntegral.hh \
EllipticIntegral.icc \
Erf.hh \
......@@ -51,6 +52,8 @@ pkginclude_HEADERS = \
Exponential.hh \
FixedConstant.hh \
FloatingConstant.hh \
FourierFit.hh \
FourierFit.icc \
FunctionComposition.hh \
FunctionConvolution.hh \
FunctionDifference.hh \
......@@ -70,6 +73,8 @@ pkginclude_HEADERS = \
IncompleteGamma.hh \
InterpolatingPolynomial.hh \
Landau.hh \
LegendreFit.hh \
LegendreFit.icc \
KroneckerDelta.hh \
Legendre.hh \
Legendre.icc \
......@@ -106,6 +111,8 @@ pkginclude_HEADERS = \
SphericalBessel.icc \
SphericalNeumann.hh \
SphericalNeumann.icc \
SphericalHarmonicFit.hh \
SphericalHarmonicFit.icc \
Sqrt.hh \
Square.hh \
SymToArgAdaptor.hh \
......
// -*- C++ -*-
// $Id:
//---------------------SphericalHarmonicFit------------------------------------------//
// //
// Class SphericalHarmonicFit. This is a fitting function consisting of a super //
// position of N legendre polynomials. Cascading fractions and phases are //
// the input parameters. Function is normalized to one (on [-1,1]) //
// Joe Boudreau, Petar Maksimovic, January 2000 //
// //
//--------------------------------------------------------------------------//
#ifndef SphericalHarmonicFit_h
#define SphericalHarmonicFit_h
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class SphericalHarmonicFit : public AbsFunction {
FUNCTION_OBJECT_DEF(SphericalHarmonicFit)
public:
// Constructor. Builds all the
SphericalHarmonicFit(unsigned int LMAX);
// Copy constructor
SphericalHarmonicFit(const SphericalHarmonicFit &right);
// Destructor
virtual ~SphericalHarmonicFit();
// Dimensionality=2. They are; cosTheta (not theta) and phi
virtual unsigned int dimensionality() const {return 2;}
// Retreive function value
virtual double operator ()(double argument) const; // Gives an error.
virtual double operator ()(const Argument & a) const; // Must use this one
unsigned int numComponents() const {
return (LMAX+1)*(LMAX+1)-1;
}
Parameter *getFraction(unsigned int i);
const Parameter *getFraction(unsigned int i) const;
Parameter *getPhase(unsigned int i);
const Parameter *getPhase(unsigned int i) const;
private:
// It is illegal to assign an adjustable constant
const SphericalHarmonicFit & operator=(const SphericalHarmonicFit &right);
//
const unsigned int LMAX;
std::vector <Genfun::Parameter *> fraction;
std::vector <Genfun::Parameter *> phase;
};
} // namespace Genfun
#include "CLHEP/GenericFunctions/SphericalHarmonicFit.icc"
#endif
// -*- C++ -*-
// $Id:
#include <sstream>
#include <cmath>
#include <gsl/gsl_sf_legendre.h>
#include <complex>
#include <cstdlib>
#include <stdexcept>
namespace Genfun {
FUNCTION_OBJECT_IMP(SphericalHarmonicFit)
inline
SphericalHarmonicFit::SphericalHarmonicFit(unsigned int LMAX):
LMAX(LMAX)
{
// SKIP L=1. We take the phase of this to be zero and the fraction
// to be whatever is left over after the other fraction take.
for (int l=1;l<=LMAX;l++) {
for (int m=-l;m<=l;m++) {
std::ostringstream stream;
stream << "Fraction l=" << l << " ;m=" << m << std::endl;
fraction.push_back(new Parameter(stream.str(), 0.5, 0.0, 1.0));
}
}
for (int l=1;l<=LMAX;l++) {
for (int m=-l;m<=l;m++) {
std::ostringstream stream;
stream << "Phase l=" << l << " ;m=" << m << std::endl;
phase.push_back(new Parameter(stream.str(), M_PI, 0.0, 2.0*M_PI));
}
}
}
inline
SphericalHarmonicFit::~SphericalHarmonicFit() {
for (unsigned int i=0;i<numComponents();i++) {
delete fraction[i];
delete phase[i];
}
}
inline
SphericalHarmonicFit::SphericalHarmonicFit(const SphericalHarmonicFit & right):
LMAX(right.LMAX)
{
for (int i=0;i<right.numComponents();i++) {
fraction.push_back(new Parameter(*right.fraction[i]));
phase.push_back(new Parameter(*right.phase[i]));
}
}
inline
double SphericalHarmonicFit::operator() (double ) const {
throw std::runtime_error("Dimensionality error in SphericalHarmonicFit");
return 0;
}
inline
double SphericalHarmonicFit::operator() (const Argument & a ) const {
double x = a[0];
double phi=a[1];
// Note, the calling sequence of the GSL Special Function forces us to
// transpose Plm from its "natural" order.. It is addressed as P[m][l].
double Plm[LMAX+1][LMAX+1];
for (int m=0;m<=LMAX;m++) {
gsl_sf_legendre_sphPlm_array (LMAX, m, x, Plm[m]);
}
unsigned int n=(LMAX+1)*(LMAX+1)-1;
int l=LMAX;
int m=LMAX;
std::complex<double> P=0.0;
std::complex<double> I(0,1.0);
double f=1.0;
double sumSq=0;
while (1) {
int LP=l-abs(m);
if (l==0) {
double fn=1.0;
double Pn=Plm[abs(m)][LP];
if (!finite(Pn)) {
throw std::runtime_error("Non-finite Pn l=0");
}
P+=(sqrt(std::max(0.0,f*fn))*Pn);
sumSq+= f*fn;
break;
}
else {
double fn=getFraction(n-1)->getValue();
double px=getPhase(n-1)->getValue();
int S= (m<0 && abs(m)%2) ? -1:1;
double Pn= S*Plm[abs(m)][LP];
if (!finite(Pn)) {
throw std::runtime_error("Non-finite Pn l!=0");
}
P+=(exp(I*(px+m*phi))*sqrt(std::max(0.0,f*fn))*Pn);
sumSq+= f*fn;
f*=(1-fn);
n--;
if (m>-l) {
m--;
}
else {
l--;
m=l;
}
}
}
double retVal=std::norm(P);
if (!finite(retVal)) {
throw std::runtime_error("Non-finite return value in SphericalHarmonicFit");
}
return retVal;
}
inline
Parameter *SphericalHarmonicFit::getFraction(unsigned int i) {
return fraction[i];
}
inline
const Parameter *SphericalHarmonicFit::getFraction(unsigned int i) const{
return fraction[i];
}
inline
Parameter *SphericalHarmonicFit::getPhase(unsigned int i) {
return phase[i];
}
inline
const Parameter *SphericalHarmonicFit::getPhase(unsigned int i) const{
return phase[i];
}