SphericalHarmonicFit.icc

Go to the documentation of this file.

// -*- C++ -*-
// $Id: 
#include <sstream>
#include <cmath>
#include <gsl/gsl_sf_legendre.h>
#include <complex>
#include <cstdlib>
#include <stdexcept>
#include "CLHEP/GenericFunctions/ClebschGordanCoefficientSet.hh"
namespace Genfun {

FUNCTION_OBJECT_IMP(SphericalHarmonicFit)

class SphericalHarmonicFit::Clockwork {
  
public:
  
  Clockwork(unsigned int LMAX):LMAX(LMAX),coefficientsA(LMAX),coefficientsASq(2*LMAX) {}

  struct MStruct {
    unsigned int M;
    Genfun::Parameter *fractionAbsMOrHigher;
    Genfun::Parameter *fractionMPositive;
    Genfun::Parameter *phaseMPlus;
    Genfun::Parameter *phaseMMinus;
  };
  
  struct LStruct {
    unsigned  int      L;
    Genfun::Parameter *fractionLOrHigher;
    Genfun::Parameter *phaseLM0;
    std::vector<MStruct> mstruct;
  };

  std::vector<LStruct>   lstruct;
  const unsigned   int   LMAX;


  SphericalHarmonicCoefficientSet coefficientsA;
  SphericalHarmonicCoefficientSet coefficientsASq;
  ClebschGordanCoefficientSet     ClebschGordan;

  void recomputeCoefficients() {

    // 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 ampSq=0.0;
    
    std::complex<double> I(0,1.0);
    double f=1.0;
    double fThisSum=0.0;
    for (unsigned int l=0;l<=LMAX;l++) {
      
      // lStructThis is zero if l==0;
      // lStructNext is zero if l==LMAX;
      const LStruct *lStructThis= (l==0    ? NULL: & lstruct[l-1]);
      const LStruct *lStructNext= (l==LMAX ? NULL: & lstruct[l]);
      double fHigher = lStructNext ? lStructNext->fractionLOrHigher->getValue() : NULL;
      double fThis   = f*(1-fHigher);
      fThisSum+=fThis;
      
      double g=1.0;
      double gThisSum=0.0;
      for (int m=0;m<=int(l);m++) {
        
        // mStructThis is zero if m==0;
        // mStructNext is zero if m==l;
        const MStruct *mStructThis= ((m==0 || !lStructThis) ? NULL: & lStructThis->mstruct[m-1]);
        const MStruct *mStructNext= (m==int(l) ? NULL: & lStructThis->mstruct[m]);
        double gHigher = mStructNext ? mStructNext->fractionAbsMOrHigher->getValue() : NULL;
        double gThis   = g*(1-gHigher);
        gThisSum+=gThis;
        
        if (fThis<0) {
          std::cout << "L-fraction correction" << fThis << "-->0" << std::endl;
          fThis=0.0;
        }
        if (gThis<0) {
          std::cout << "M-fraction correction" << gThis << "-->0" << std::endl;
          gThis=0.0;
        }
        double px=0.0; // phase
        if (m==0) {
          if (lStructThis) {
            double amplitude = sqrt(fThis*gThis);
            px = lStructThis->phaseLM0->getValue();;
            coefficientsA(l,m)=exp(I*px)*amplitude;
          }
          // L=0 occurs here:
          else {
            double amplitude = sqrt(fThis*gThis);
            coefficientsA(l,m)=exp(I*px)*amplitude;
          }
        }
        // Split it between positive and negative:
        else {
          {
            double amplitude = sqrt(fThis*gThis*mStructThis->fractionMPositive->getValue());
            px = mStructThis->phaseMPlus->getValue();;
            coefficientsA(l,m)=exp(I*px)*amplitude;
          }
          {
            double amplitude = sqrt(fThis*gThis*(1-mStructThis->fractionMPositive->getValue()));
            px = mStructThis->phaseMMinus->getValue();;
            coefficientsA(l,-m)=exp(I*px)*amplitude;
          }
        }
        g*=gHigher;
      }
      f*=fHigher;
    }
  }
};


inline
SphericalHarmonicFit::SphericalHarmonicFit(unsigned int LMAX):
  c(new Clockwork(LMAX))
{
  for (unsigned int l=1;l<=LMAX;l++) {
    Clockwork::LStruct lstruct;
    lstruct.L=l;
    {
      std::ostringstream stream; 
      stream << "Fraction L>=" << l;
      lstruct.fractionLOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
    }
    {
      std::ostringstream stream; 
      stream << "Phase L=" << l << "; M=0";
      lstruct.phaseLM0= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
    }
    for (unsigned int m=1;m<=l;m++) {
      Clockwork::MStruct mstruct;
      mstruct.M=m;
      {
        std::ostringstream stream; 
        stream << "Fraction L= " << l << "; |M| >=" << m;
        mstruct.fractionAbsMOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
      }
      {
        std::ostringstream stream; 
        stream << "Fraction L=" << l << "; M=+" << m ;
        mstruct.fractionMPositive= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
      }
      {
        std::ostringstream stream; 
        stream << "Phase L=" << l << "; M=+" << m ;
        mstruct.phaseMPlus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
      }
      {
        std::ostringstream stream; 
        stream << "Phase L=" << l << "; M=-" << m ;
        mstruct.phaseMMinus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
      }
      
      lstruct.mstruct.push_back(mstruct);
    }
    c->lstruct.push_back(lstruct);
  }
}


inline
SphericalHarmonicFit::~SphericalHarmonicFit() {
  for (unsigned int i=0;i<c->lstruct.size();i++) {
    delete c->lstruct[i].fractionLOrHigher;
    delete c->lstruct[i].phaseLM0;
    for (unsigned int j=0;j<c->lstruct[i].mstruct.size();j++) {
      delete c->lstruct[i].mstruct[j].fractionAbsMOrHigher;
      delete c->lstruct[i].mstruct[j].fractionMPositive;
      delete c->lstruct[i].mstruct[j].phaseMPlus;
      delete c->lstruct[i].mstruct[j].phaseMMinus;
    }
  }
  delete c;
}

  inline
  SphericalHarmonicFit::SphericalHarmonicFit(const SphericalHarmonicFit & right):
    AbsFunction(),
    c(new Clockwork(right.c->LMAX))
  {
    for (unsigned int i=0;i<right.c->lstruct.size();i++) {
      Clockwork::LStruct lstruct;
      lstruct.L= right.c->lstruct[i].L;
      lstruct.fractionLOrHigher = new Parameter(*right.c->lstruct[i].fractionLOrHigher);
      lstruct.phaseLM0          = new Parameter(*right.c->lstruct[i].phaseLM0);
      for (unsigned int j=0;j<right.c->lstruct[i].mstruct.size();j++) {
        Clockwork::MStruct mstruct;
        mstruct.M=right.c->lstruct[i].mstruct[j].M;
        mstruct.fractionAbsMOrHigher=new Parameter(*right.c->lstruct[i].mstruct[j].fractionAbsMOrHigher);
        mstruct.fractionMPositive   =new Parameter(*right.c->lstruct[i].mstruct[j].fractionMPositive);
        mstruct.phaseMPlus          =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMPlus);
        mstruct.phaseMMinus         =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMMinus);
        lstruct.mstruct.push_back(mstruct);
      }
      c->lstruct.push_back(lstruct);
    }
  }
  
  inline
  double SphericalHarmonicFit::operator() (double ) const {
    throw std::runtime_error("Dimensionality error in SphericalHarmonicFit");
    return 0;
  }
  
  inline
  double SphericalHarmonicFit::operator() (const Argument & a ) const {
    unsigned int LMAX=c->LMAX;
    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];  
    std::vector< std::vector<double> > Plm(LMAX+1);
    for (int m=0;m<=int(LMAX);m++) {
      Plm[m].resize(LMAX+1);
      gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
    }

    c->recomputeCoefficients();
    std::complex<double> P=0.0;
    std::complex<double> I(0,1.0);
    for (unsigned int l=0;l<=LMAX;l++) {
      for (int m=0;m<=int(l);m++) {
        {
          int LP=l-abs(m);
          double Pn= Plm[abs(m)][LP];
          
          if (!finite(Pn)) return 0.0;

          // Once for positive m (in all cases):
          P+=(c->coefficientsA(l,m)*Pn*exp(I*(m*phi)));
          // Once for negative m (skip if m==0);
          if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficientsA(l,-m)*Pn*exp(-I*(m*phi)));
        }
      }
    }
    
    double retVal=std::norm(P);
    if (!finite(retVal)) {
      return 0.0;
    }
    return retVal;
  }
  
  inline
  unsigned int SphericalHarmonicFit::lMax() const {
    return c->LMAX;
  }
  
  inline
  unsigned int SphericalHarmonicFit::numComponents() const {
    return (c->LMAX+1)*(c->LMAX+1)-1;
  }
  
  // The fraction of Amplitude sq which is L or higher
  inline
  Parameter *SphericalHarmonicFit::getFractionLOrHigher(unsigned int L){
    return c->lstruct[L-1].fractionLOrHigher;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getFractionLOrHigher(unsigned int L) const {
    return c->lstruct[L-1].fractionLOrHigher;
  }

  // The phase of coefficient L, M=0;
  inline
  Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L){
    return c->lstruct[L-1].phaseLM0;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L) const{
    return c->lstruct[L-1].phaseLM0;
  }
  
  // The fraction of amplitude sq which is L which is +- M OR HIGHER
  inline
  Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
  }
  

  // The fraction of amplitude sq which is +- M, which is positive
  inline
  Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
  }

  inline
  const Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
  }


  // The phase of the positive M coefficient
  inline
  Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
  }
  
  
  // The phase of the negative M coefficient
  inline
  Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M){
    return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
  }
  
  inline
  const Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M) const{
    return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
  }
  
  inline
  const SphericalHarmonicCoefficientSet & SphericalHarmonicFit::coefficientsA() const {
    c->recomputeCoefficients();
    return c->coefficientsA;
  }

  inline
  const SphericalHarmonicCoefficientSet & SphericalHarmonicFit::coefficientsASq() const{
    c->recomputeCoefficients();
    unsigned int LMAX=c->coefficientsA.getLMax();
    for (unsigned int L=0;L<=2*LMAX;L++) {
      for (int M=-L; M<=int(L); M++) {
        c->coefficientsASq(L,M)=0.0;
        for (unsigned int l1=0;l1<=LMAX;l1++) {
          for (unsigned int l2=0;l2<=LMAX;l2++) {
            for (int m1=-l1;m1<=int(l1);m1++) {
              for (int m2=-l2;m2<=int(l2);m2++) {
                  if (m1-m2==M) {
                    if (((l1+l2) >= L) && abs(l1-l2) <=int(L))  {
                    c->coefficientsASq(L,M) += (c->coefficientsA(l1,m1)*
                                                conj(c->coefficientsA(l2,m2))*
                                                (m2%2 ? -1.0:1.0) * 
                                                sqrt((2*l1+1)*(2*l2+1)/(4*M_PI*(2*L+1)))*
                                                c->ClebschGordan(l1,l2,0,0,L,0)*c->ClebschGordan(l1,l2,m1,-m2,L,M));
                    }
                  }
              }
            }
          }
        }
      }
    }
    return c->coefficientsASq;
  }

  inline
  void SphericalHarmonicFit::recomputeCoefficients() const {
    c->recomputeCoefficients();
  }

} // end namespace Genfun