CLHEP/GenericFunctions/FourierFit.icc

0001 // -*- C++ -*-
0002 // $Id:
0003 #include <sstream>
0004 #include <cmath>
0005 #include <complex>
0006 namespace Genfun {
0007
0008 FUNCTION_OBJECT_IMP(FourierFit)
0009
0010 inline
0011 FourierFit::FourierFit(unsigned int N):
0012    N(N)
0013 {
0014   for (unsigned int i=0;i<N;i++) {
0015     std::ostringstream stream;
0016     stream << "Fraction " << i;
0017     fraction.push_back(new Parameter(stream.str(), 0.5, 0.0, 1.0));
0018   }
0019   for (unsigned int i=0;i<N;i++) {
0020     std::ostringstream stream;
0021     stream << "Phase " << i;
0022     phase.push_back(new Parameter(stream.str(), M_PI, 0.0, 2.0*M_PI));
0023   }
0024 }
0025
0026 inline
0027 FourierFit::~FourierFit() {
0028   for (unsigned int i=0;i<N;i++) {
0029     delete fraction[i];
0030     delete phase[i];
0031   }
0032 }
0033
0034 inline
0035 FourierFit::FourierFit(const FourierFit & right):
0036   N(right.N)
0037 {
0038   for (int i=0;i<N;i++) {
0039     fraction.push_back(new Parameter(*right.fraction[i]));
0040     phase.push_back(new Parameter(*right.phase[i]));
0041   }
0042 }
0043
0044 inline
0045 double FourierFit::operator() (double x) const {
0046
0047     unsigned int n=N;
0048     std::complex<double> P=0.0;
0049     std::complex<double> I(0,1.0);
0050     double f=1.0;
0051     while (1) {
0052       if (n==0) {
0053     double fn=1.0;
0054     double Pn=sqrt(1/2.0/M_PI);
0055
0056     P+=(sqrt(f*fn)*Pn);
0057     break;
0058       }
0059       else {
0060     double fn=getFraction(n-1)->getValue();
0061     double px=getPhase(n-1)->getValue();
0062     double Pn=sqrt(1/M_PI)*sin(n*x/2.0);
0063     P+=exp(I*px)*sqrt(f*fn)*Pn;
0064     f*=(1-fn);
0065     n--;
0066       }
0067     }
0068     return std::norm(P);
0069 }
0070 inline
0071 unsigned int FourierFit::order() const{
0072   return N;
0073 }
0074 inline
0075 Parameter *FourierFit::getFraction(unsigned int i) {
0076   return fraction[i];
0077 }
0078 inline
0079 const Parameter *FourierFit::getFraction(unsigned int i) const{
0080   return fraction[i];
0081 }
0082 inline
0083 Parameter *FourierFit::getPhase(unsigned int i) {
0084   return phase[i];
0085 }
0086 inline
0087 const Parameter *FourierFit::getPhase(unsigned int i) const{
0088   return phase[i];
0089 }
0090
0091
0092 } // end namespace Genfun