CLHEP 2.4.7.1
C++ Class Library for High Energy Physics
SphericalHarmonicFit.icc
Go to the documentation of this file.
1// -*- C++ -*-
2// $Id:
3#include <sstream>
4#include <cmath>
5#include <gsl/gsl_sf_legendre.h>
6#include <complex>
7#include <cstdlib>
8#include <stdexcept>
10namespace Genfun {
11
12FUNCTION_OBJECT_IMP(SphericalHarmonicFit)
13
15
16public:
17
18 Clockwork(unsigned int LMAX):LMAX(LMAX),coefficientsA(LMAX),coefficientsASq(2*LMAX) {}
19
27
28 struct LStruct {
29 unsigned int L;
32 std::vector<MStruct> mstruct;
33 };
34
35 std::vector<LStruct> lstruct;
36 const unsigned int LMAX;
37
38
42
44
45 // Note, the calling sequence of the GSL Special Function forces us to
46 // transpose Plm from its "natural" order.. It is addressed as P[m][l].
47
48 // double ampSq=0.0;
49
50 std::complex<double> I(0,1.0);
51 double f=1.0;
52 double fThisSum=0.0;
53 for (unsigned int l=0;l<=LMAX;l++) {
54
55 // lStructThis is zero if l==0;
56 // lStructNext is zero if l==LMAX;
57 const LStruct *lStructThis= (l==0 ? NULL: & lstruct[l-1]);
58 const LStruct *lStructNext= (l==LMAX ? NULL: & lstruct[l]);
59 double fHigher = lStructNext ? lStructNext->fractionLOrHigher->getValue() : NULL;
60 double fThis = f*(1-fHigher);
62
63 double g=1.0;
64 double gThisSum=0.0;
65 for (int m=0;m<=int(l);m++) {
66
67 // mStructThis is zero if m==0;
68 // mStructNext is zero if m==l;
69 const MStruct *mStructThis= ((m==0 || !lStructThis) ? NULL: & lStructThis->mstruct[m-1]);
70 const MStruct *mStructNext= (m==int(l) ? NULL: & lStructThis->mstruct[m]);
71 double gHigher = mStructNext ? mStructNext->fractionAbsMOrHigher->getValue() : NULL;
72 double gThis = g*(1-gHigher);
74
75 if (fThis<0) {
76 std::cout << "L-fraction correction" << fThis << "-->0" << std::endl;
77 fThis=0.0;
78 }
79 if (gThis<0) {
80 std::cout << "M-fraction correction" << gThis << "-->0" << std::endl;
81 gThis=0.0;
82 }
83 double px=0.0; // phase
84 if (m==0) {
85 if (lStructThis) {
86 double amplitude = sqrt(fThis*gThis);
87 px = lStructThis->phaseLM0->getValue();;
88 coefficientsA(l,m)=exp(I*px)*amplitude;
89 }
90 // L=0 occurs here:
91 else {
92 double amplitude = sqrt(fThis*gThis);
93 coefficientsA(l,m)=exp(I*px)*amplitude;
94 }
95 }
96 // Split it between positive and negative:
97 else {
98 {
99 double amplitude = sqrt(fThis*gThis*mStructThis->fractionMPositive->getValue());
100 px = mStructThis->phaseMPlus->getValue();;
101 coefficientsA(l,m)=exp(I*px)*amplitude;
102 }
103 {
104 double amplitude = sqrt(fThis*gThis*(1-mStructThis->fractionMPositive->getValue()));
105 px = mStructThis->phaseMMinus->getValue();;
106 coefficientsA(l,-m)=exp(I*px)*amplitude;
107 }
108 }
109 g*=gHigher;
110 }
111 f*=fHigher;
112 }
113 }
114};
115
116
117inline
119 c(new Clockwork(LMAX))
120{
121 for (unsigned int l=1;l<=LMAX;l++) {
122 Clockwork::LStruct lstruct;
123 lstruct.L=l;
124 {
125 std::ostringstream stream;
126 stream << "Fraction L>=" << l;
127 lstruct.fractionLOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
128 }
129 {
130 std::ostringstream stream;
131 stream << "Phase L=" << l << "; M=0";
132 lstruct.phaseLM0= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
133 }
134 for (unsigned int m=1;m<=l;m++) {
135 Clockwork::MStruct mstruct;
136 mstruct.M=m;
137 {
138 std::ostringstream stream;
139 stream << "Fraction L= " << l << "; |M| >=" << m;
140 mstruct.fractionAbsMOrHigher= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
141 }
142 {
143 std::ostringstream stream;
144 stream << "Fraction L=" << l << "; M=+" << m ;
145 mstruct.fractionMPositive= new Genfun::Parameter(stream.str(), 0.5, 0, 1);
146 }
147 {
148 std::ostringstream stream;
149 stream << "Phase L=" << l << "; M=+" << m ;
150 mstruct.phaseMPlus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
151 }
152 {
153 std::ostringstream stream;
154 stream << "Phase L=" << l << "; M=-" << m ;
155 mstruct.phaseMMinus= new Genfun::Parameter(stream.str(), M_PI, -2*M_PI, 2*M_PI);
156 }
157
158 lstruct.mstruct.push_back(mstruct);
159 }
160 c->lstruct.push_back(lstruct);
161 }
162}
163
164
165inline
167 for (unsigned int i=0;i<c->lstruct.size();i++) {
168 delete c->lstruct[i].fractionLOrHigher;
169 delete c->lstruct[i].phaseLM0;
170 for (unsigned int j=0;j<c->lstruct[i].mstruct.size();j++) {
171 delete c->lstruct[i].mstruct[j].fractionAbsMOrHigher;
172 delete c->lstruct[i].mstruct[j].fractionMPositive;
173 delete c->lstruct[i].mstruct[j].phaseMPlus;
174 delete c->lstruct[i].mstruct[j].phaseMMinus;
175 }
176 }
177 delete c;
178}
179
180 inline
182 AbsFunction(),
183 c(new Clockwork(right.c->LMAX))
184 {
185 for (unsigned int i=0;i<right.c->lstruct.size();i++) {
186 Clockwork::LStruct lstruct;
187 lstruct.L= right.c->lstruct[i].L;
188 lstruct.fractionLOrHigher = new Parameter(*right.c->lstruct[i].fractionLOrHigher);
189 lstruct.phaseLM0 = new Parameter(*right.c->lstruct[i].phaseLM0);
190 for (unsigned int j=0;j<right.c->lstruct[i].mstruct.size();j++) {
191 Clockwork::MStruct mstruct;
192 mstruct.M=right.c->lstruct[i].mstruct[j].M;
193 mstruct.fractionAbsMOrHigher=new Parameter(*right.c->lstruct[i].mstruct[j].fractionAbsMOrHigher);
194 mstruct.fractionMPositive =new Parameter(*right.c->lstruct[i].mstruct[j].fractionMPositive);
195 mstruct.phaseMPlus =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMPlus);
196 mstruct.phaseMMinus =new Parameter(*right.c->lstruct[i].mstruct[j].phaseMMinus);
197 lstruct.mstruct.push_back(mstruct);
198 }
199 c->lstruct.push_back(lstruct);
200 }
201 }
202
203 inline
204 double SphericalHarmonicFit::operator() (double ) const {
205 throw std::runtime_error("Dimensionality error in SphericalHarmonicFit");
206 return 0;
207 }
208
209 inline
211 unsigned int LMAX=c->LMAX;
212 double x = a[0];
213 double phi=a[1];
214
215 // Note, the calling sequence of the GSL Special Function forces us to
216 // transpose Plm from its "natural" order.. It is addressed as P[m][l].
217
218 //double Plm[LMAX+1][LMAX+1];
219 std::vector< std::vector<double> > Plm(LMAX+1);
220 for (int m=0;m<=int(LMAX);m++) {
221 Plm[m].resize(LMAX+1);
222 gsl_sf_legendre_sphPlm_array (LMAX, m, x, &*Plm[m].begin());
223 }
224
226 std::complex<double> P=0.0;
227 std::complex<double> I(0,1.0);
228 for (unsigned int l=0;l<=LMAX;l++) {
229 for (int m=0;m<=int(l);m++) {
230 {
231 int LP=l-abs(m);
232 double Pn= Plm[abs(m)][LP];
233
234 if (!finite(Pn)) return 0.0;
235
236 // Once for positive m (in all cases):
237 P+=(c->coefficientsA(l,m)*Pn*exp(I*(m*phi)));
238 // Once for negative m (skip if m==0);
239 if (m!=0) P+= ( (m%2 ?-1.0:1.0)*c->coefficientsA(l,-m)*Pn*exp(-I*(m*phi)));
240 }
241 }
242 }
243
244 double retVal=std::norm(P);
245 if (!finite(retVal)) {
246 return 0.0;
247 }
248 return retVal;
249 }
250
251 inline
252 unsigned int SphericalHarmonicFit::lMax() const {
253 return c->LMAX;
254 }
255
256 inline
258 return (c->LMAX+1)*(c->LMAX+1)-1;
259 }
260
261 // The fraction of Amplitude sq which is L or higher
262 inline
264 return c->lstruct[L-1].fractionLOrHigher;
265 }
266
267 inline
269 return c->lstruct[L-1].fractionLOrHigher;
270 }
271
272 // The phase of coefficient L, M=0;
273 inline
275 return c->lstruct[L-1].phaseLM0;
276 }
277
278 inline
279 const Parameter *SphericalHarmonicFit::getPhaseLM0(unsigned int L) const{
280 return c->lstruct[L-1].phaseLM0;
281 }
282
283 // The fraction of amplitude sq which is L which is +- M OR HIGHER
284 inline
286 return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
287 }
288
289 inline
290 const Parameter *SphericalHarmonicFit::getFractionAbsMOrHigher(unsigned int L, unsigned int M) const{
291 return c->lstruct[L-1].mstruct[M-1].fractionAbsMOrHigher;
292 }
293
294
295 // The fraction of amplitude sq which is +- M, which is positive
296 inline
298 return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
299 }
300
301 inline
302 const Parameter *SphericalHarmonicFit::getFractionMPositive(unsigned int L, unsigned int M) const{
303 return c->lstruct[L-1].mstruct[M-1].fractionMPositive;
304 }
305
306
307 // The phase of the positive M coefficient
308 inline
309 Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M){
310 return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
311 }
312
313 inline
314 const Parameter *SphericalHarmonicFit::getPhaseMPlus(unsigned int L, unsigned int M) const{
315 return c->lstruct[L-1].mstruct[M-1].phaseMPlus;
316 }
317
318
319 // The phase of the negative M coefficient
320 inline
321 Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M){
322 return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
323 }
324
325 inline
326 const Parameter *SphericalHarmonicFit::getPhaseMMinus(unsigned int L, unsigned int M) const{
327 return c->lstruct[L-1].mstruct[M-1].phaseMMinus;
328 }
329
330 inline
335
336 inline
339 unsigned int LMAX=c->coefficientsA.getLMax();
340 for (unsigned int L=0;L<=2*LMAX;L++) {
341 for (int M=-L; M<=int(L); M++) {
342 c->coefficientsASq(L,M)=0.0;
343 for (unsigned int l1=0;l1<=LMAX;l1++) {
344 for (unsigned int l2=0;l2<=LMAX;l2++) {
345 for (int m1=-l1;m1<=int(l1);m1++) {
346 for (int m2=-l2;m2<=int(l2);m2++) {
347 if (m1-m2==M) {
348 if (((l1+l2) >= L) && abs(l1-l2) <=int(L)) {
349 c->coefficientsASq(L,M) += (c->coefficientsA(l1,m1)*
350 conj(c->coefficientsA(l2,m2))*
351 (m2%2 ? -1.0:1.0) *
352 sqrt((2*l1+1)*(2*l2+1)/(4*M_PI*(2*L+1)))*
353 c->ClebschGordan(l1,l2,0,0,L,0)*c->ClebschGordan(l1,l2,m1,-m2,L,M));
354 }
355 }
356 }
357 }
358 }
359 }
360 }
361 }
362 return c->coefficientsASq;
363 }
364
365 inline
369
370} // end namespace Genfun
#define FUNCTION_OBJECT_IMP(classname)
SphericalHarmonicCoefficientSet coefficientsASq
SphericalHarmonicCoefficientSet coefficientsA
const SphericalHarmonicCoefficientSet & coefficientsASq() const
Parameter * getFractionAbsMOrHigher(unsigned int L, unsigned int M)
Parameter * getPhaseMMinus(unsigned int L, unsigned int M)
virtual double operator()(double argument) const override
const SphericalHarmonicCoefficientSet & coefficientsA() const
Parameter * getPhaseMPlus(unsigned int L, unsigned int M)
Parameter * getFractionMPositive(unsigned int L, unsigned int M)
Parameter * getPhaseLM0(unsigned int L)
Parameter * getFractionLOrHigher(unsigned int L)
Definition Abs.hh:14