CLHEP/GenericFunctions/SphericalHarmonicFit.icc

0001 // -*- C++ -*-
0002 // $Id:
0003 #include <sstream>
0004 #include <cmath>
0005 #include <gsl/gsl_sf_legendre.h>
0006 #include <complex>
0007 #include <cstdlib>
0008 #include <stdexcept>
0009 #include "CLHEP/GenericFunctions/ClebschGordanCoefficientSet.hh"
0010 namespace Genfun {
0011
0012 FUNCTION_OBJECT_IMP(SphericalHarmonicFit)
0013
0014 class SphericalHarmonicFit::Clockwork {
0015
0016 public:
0017
0018   Clockwork(unsigned int LMAX):LMAX(LMAX),coefficientsA(LMAX),coefficientsASq(2*LMAX) {}
0019
0020   struct MStruct {
0021     unsigned int M;
0022     Genfun::Parameter *fractionAbsMOrHigher;
0023     Genfun::Parameter *fractionMPositive;
0024     Genfun::Parameter *phaseMPlus;
0025     Genfun::Parameter *phaseMMinus;
0026   };
0027
0028   struct LStruct {
0029     unsigned  int      L;
0030     Genfun::Parameter *fractionLOrHigher;
0031     Genfun::Parameter *phaseLM0;
0032     std::vector<MStruct> mstruct;
0033   };
0034
0035   std::vector<LStruct>   lstruct;
0036   const unsigned   int   LMAX;
0037
0038
0039   SphericalHarmonicCoefficientSet coefficientsA;
0040   SphericalHarmonicCoefficientSet coefficientsASq;
0041   ClebschGordanCoefficientSet     ClebschGordan;
0042
0043   void recomputeCoefficients() {
0044
0045     // Note, the calling sequence of the GSL Special Function forces us to
0046     // transpose Plm from its "natural" order.. It is addressed as P[m][l].
0047
0048     //  double ampSq=0.0;
0049
0050     std::complex<double> I(0,1.0);
0051     double f=1.0;
0052     double fThisSum=0.0;
0053     for (unsigned int l=0;l<=LMAX;l++) {
0054
0055       // lStructThis is zero if l==0;
0056       // lStructNext is zero if l==LMAX;
0057       const LStruct *lStructThis= (l==0    ? NULL: & lstruct[l-1]);
0058       const LStruct *lStructNext= (l==LMAX ? NULL: & lstruct[l]);
0059       double fHigher = lStructNext ? lStructNext->fractionLOrHigher->getValue() : NULL;
0060       double fThis   = f*(1-fHigher);
0061       fThisSum+=fThis;
0062
0063       double g=1.0;
0064       double gThisSum=0.0;
0065       for (int m=0;m<=int(l);m++) {
0066
0067     // mStructThis is zero if m==0;
0068     // mStructNext is zero if m==l;
0069     const MStruct *mStructThis= ((m==0 || !lStructThis) ? NULL: & lStructThis->mstruct[m-1]);
0070     const MStruct *mStructNext= (m==int(l) ? NULL: & lStructThis->mstruct[m]);
0071     double gHigher = mStructNext ? mStructNext->fractionAbsMOrHigher->getValue() : NULL;
0072     double gThis   = g*(1-gHigher);
0073     gThisSum+=gThis;
0074
0075     if (fThis<0) {
0076       std::cout << "L-fraction correction" << fThis << "-->0" << std::endl;
0077       fThis=0.0;
0078     }
0079     if (gThis<0) {
0080       std::cout << "M-fraction correction" << gThis << "-->0" << std::endl;
0081       gThis=0.0;
0082     }
0083     double px=0.0; // phase
0084     if (m==0) {
0085       if (lStructThis) {
0086         double amplitude = sqrt(fThis*gThis);
0087         px = lStructThis->phaseLM0->getValue();;
0088         coefficientsA(l,m)=exp(I*px)*amplitude;
0089       }
0090       // L=0 occurs here:
0091       else {
0092         double amplitude = sqrt(fThis*gThis);
0093         coefficientsA(l,m)=exp(I*px)*amplitude;
0094       }
0095     }
0096     // Split it between positive and negative:
0097     else {
0098       {
0099         double amplitude = sqrt(fThis*gThis*mStructThis->fractionMPositive->getValue());
0100         px = mStructThis->phaseMPlus->getValue();;
0101         coefficientsA(l,m)=exp(I*px)*amplitude;
0102       }
0103       {
0104         double amplitude = sqrt(fThis*gThis*(1-mStructThis->fractionMPositive->getValue()));
0105         px = mStructThis->phaseMMinus->getValue();;
0106         coefficientsA(l,-m)=exp(I*px)*amplitude;
0107       }
0108     }
0109     g*=gHigher;
0110       }
0111       f*=fHigher;
0112     }
0113   }
0114 };
0115
0116
0117 inline
0118 SphericalHarmonicFit::SphericalHarmonicFit(unsigned int LMAX):
0119   c(new Clockwork(LMAX))
0120 {
0121   for (unsigned int l=1;l<=LMAX;l++) {
0122     Clockwork::LStruct lstruct;
0123     lstruct.L=l;
0124     {
0125       std::ostringstream stream;
0126       stream << "Fraction L>=" << l;
0127       lstruct.fractionLOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
0128     }
0129     {
0130       std::ostringstream stream;
0131       stream << "Phase L=" << l << "; M=0";
0132       lstruct.phaseLM0= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
0133     }
0134     for (unsigned int m=1;m<=l;m++) {
0135       Clockwork::MStruct mstruct;
0136       mstruct.M=m;
0137       {
0138     std::ostringstream stream;
0139     stream << "Fraction L= " << l << "; |M| >=" << m;
0140     mstruct.fractionAbsMOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
0141       }
0142       {
0143     std::ostringstream stream;
0144     stream << "Fraction L=" << l << "; M=+" << m ;
0145     mstruct.fractionMPositive= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
0146       }
0147       {
0148     std::ostringstream stream;
0149     stream << "Phase L=" << l << "; M=+" << m ;
0150     mstruct.phaseMPlus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
0151       }
0152       {
0153     std::ostringstream stream;
0154     stream << "Phase L=" << l << "; M=-" << m ;
0155     mstruct.phaseMMinus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
0156       }
0157
0158       lstruct.mstruct.push_back(mstruct);
0159     }
0160     c->lstruct.push_back(lstruct);
0161   }
0162 }
0163
0164
0165 inline
0166 SphericalHarmonicFit::~SphericalHarmonicFit() {
0167   for (unsigned int i=0;i<c->lstruct.size();i++) {
0168     delete c->lstruct[i].fractionLOrHigher;
0169     delete c->lstruct[i].phaseLM0;
0170     for (unsigned int j=0;j<c->lstruct[i].mstruct.size();j++) {
0171       delete c->lstruct[i].mstruct[j].fractionAbsMOrHigher;
0172       delete c->lstruct[i].mstruct[j].fractionMPositive;
0173       delete c->lstruct[i].mstruct[j].phaseMPlus;
0174       delete c->lstruct[i].mstruct[j].phaseMMinus;
0175     }
0176   }
0177   delete c;
0178 }
0179
0180   inline
0181   SphericalHarmonicFit::SphericalHarmonicFit(const SphericalHarmonicFit & right):
0182     AbsFunction(),
0183     c(new Clockwork(right.c->LMAX))
0184   {
0185     for (unsigned int i=0;i<right.c->lstruct.size();i++) {
0186       Clockwork::LStruct lstruct;
0187       lstruct.L= right.c->lstruct[i].L;
0188       lstruct.fractionLOrHigher = new Parameter(*right.c->lstruct[i].fractionLOrHigher);
0189       lstruct.phaseLM0          = new Parameter(*right.c->lstruct[i].phaseLM0);
0190       for (unsigned int j=0;j<right.c->lstruct[i].mstruct.size();j++) {
0191     Clockwork::MStruct mstruct;
0192     mstruct.M=right.c->lstruct[i].mstruct[j].M;
0193     mstruct.fractionAbsMOrHigher=new Parameter(*right.c->lstruct[i].mstruct[j].fractionAbsMOrHigher);
0194     mstruct.fractionMPositive   =new Parameter(*right.c->lstruct[i].mstruct[j].fractionMPositive);
0195     mstruct.phaseMPlus          =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMPlus);
0196     mstruct.phaseMMinus         =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMMinus);
0197     lstruct.mstruct.push_back(mstruct);
0198       }
0199       c->lstruct.push_back(lstruct);
0200     }
0201   }
0202
0203   inline
0204   double SphericalHarmonicFit::operator() (double ) const {
0205     throw std::runtime_error("Dimensionality error in SphericalHarmonicFit");
0206     return 0;
0207   }
0208
0209   inline
0210   double SphericalHarmonicFit::operator() (const Argument & a ) const {
0211     unsigned int LMAX=c->LMAX;
0212     double x = a[0];
0213     double phi=a[1];
0214
0215     // Note, the calling sequence of the GSL Special Function forces us to
0216     // transpose Plm from its "natural" order.. It is addressed as P[m][l].
0217
0218     //double Plm[LMAX+1][LMAX+1];
0219     std::vector< std::vector<double> > Plm(LMAX+1);
0220     for (int m=0;m<=int(LMAX);m++) {
0221       Plm[m].resize(LMAX+1);
0222       gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
0223     }
0224
0225     c->recomputeCoefficients();
0226     std::complex<double> P=0.0;
0227     std::complex<double> I(0,1.0);
0228     for (unsigned int l=0;l<=LMAX;l++) {
0229       for (int m=0;m<=int(l);m++) {
0230     {
0231       int LP=l-abs(m);
0232       double Pn= Plm[abs(m)][LP];
0233
0234       if (!finite(Pn)) return 0.0;
0235
0236       // Once for positive m (in all cases):
0237       P+=(c->coefficientsA(l,m)*Pn*exp(I*(m*phi)));
0238       // Once for negative m (skip if m==0);
0239       if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficientsA(l,-m)*Pn*exp(-I*(m*phi)));
0240     }
0241       }
0242     }
0243
0244     double retVal=std::norm(P);
0245     if (!finite(retVal)) {
0246       return 0.0;
0247     }
0248     return retVal;
0249   }
0250
0251   inline
0252   unsigned int SphericalHarmonicFit::lMax() const {
0253     return c->LMAX;
0254   }
0255
0256   inline
0257   unsigned int SphericalHarmonicFit::numComponents() const {
0258     return (c->LMAX+1)*(c->LMAX+1)-1;
0259   }
0260
0261   // The fraction of Amplitude sq which is L or higher
0262   inline
0263   Parameter *SphericalHarmonicFit::getFractionLOrHigher(unsigned int L){
0264     return c->lstruct[L-1].fractionLOrHigher;
0265   }
0266
0267   inline
0268   const Parameter *SphericalHarmonicFit::getFractionLOrHigher(unsigned int L) const {
0269     return c->lstruct[L-1].fractionLOrHigher;
0270   }
0271
0272   // The phase of coefficient L, M=0;
0273   inline
0274   Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L){
0275     return c->lstruct[L-1].phaseLM0;
0276   }
0277
0278   inline
0279   const Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L) const{
0280     return c->lstruct[L-1].phaseLM0;
0281   }
0282
0283   // The fraction of amplitude sq which is L which is +- M OR HIGHER
0284   inline
0285   Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M){
0286     return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
0287   }
0288
0289   inline
0290   const Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M) const{
0291     return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
0292   }
0293
0294
0295   // The fraction of amplitude sq which is +- M, which is positive
0296   inline
0297   Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M){
0298     return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
0299   }
0300
0301   inline
0302   const Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M) const{
0303     return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
0304   }
0305
0306
0307   // The phase of the positive M coefficient
0308   inline
0309   Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M){
0310     return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
0311   }
0312
0313   inline
0314   const Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M) const{
0315     return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
0316   }
0317
0318
0319   // The phase of the negative M coefficient
0320   inline
0321   Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M){
0322     return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
0323   }
0324
0325   inline
0326   const Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M) const{
0327     return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
0328   }
0329
0330   inline
0331   const SphericalHarmonicCoefficientSet & SphericalHarmonicFit::coefficientsA() const {
0332     c->recomputeCoefficients();
0333     return c->coefficientsA;
0334   }
0335
0336   inline
0337   const SphericalHarmonicCoefficientSet & SphericalHarmonicFit::coefficientsASq() const{
0338     c->recomputeCoefficients();
0339     unsigned int LMAX=c->coefficientsA.getLMax();
0340     for (unsigned int L=0;L<=2*LMAX;L++) {
0341       for (int M=-L; M<=int(L); M++) {
0342     c->coefficientsASq(L,M)=0.0;
0343     for (unsigned int l1=0;l1<=LMAX;l1++) {
0344       for (unsigned int l2=0;l2<=LMAX;l2++) {
0345         for (int m1=-l1;m1<=int(l1);m1++) {
0346           for (int m2=-l2;m2<=int(l2);m2++) {
0347           if (m1-m2==M) {
0348             if (((l1+l2) >= L) && abs(l1-l2) <=int(L))  {
0349             c->coefficientsASq(L,M) += (c->coefficientsA(l1,m1)*
0350                         conj(c->coefficientsA(l2,m2))*
0351                         (m2%2 ? -1.0:1.0) *
0352                         sqrt((2*l1+1)*(2*l2+1)/(4*M_PI*(2*L+1)))*
0353                         c->ClebschGordan(l1,l2,0,0,L,0)*c->ClebschGordan(l1,l2,m1,-m2,L,M));
0354             }
0355           }
0356           }
0357         }
0358       }
0359     }
0360       }
0361     }
0362     return c->coefficientsASq;
0363   }
0364
0365   inline
0366   void SphericalHarmonicFit::recomputeCoefficients() const {
0367     c->recomputeCoefficients();
0368   }
0369
0370 } // end namespace Genfun