Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
RooSpHarmonic.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitModels *
4  * File: $Id$
5  * Authors: *
6  * GR, Gerhard Raven, Nikhef & VU, Gerhard.Raven@nikhef.nl
7  * *
8  * Copyright (c) 2010, Nikhef & VU. All rights reserved.
9  * *
10  * Redistribution and use in source and binary forms, *
11  * with or without modification, are permitted according to the terms *
12  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
13  *****************************************************************************/
14 
15 /** \class RooSpHarmonic
16  \ingroup Roofit
17 
18  Implementation of the so-called real spherical harmonics, using the orthonormal normalization,
19 which are related to spherical harmonics as:
20 \f[
21  Y_{l0} = Y_l^0 (m=0) \\
22  Y_{lm} = \frac{1}{\sqrt{2}} \left( Y_l^m + (-1)^m Y_l^{-m} \right) (m>0) \\
23  Y_{lm} = \frac{1}{i\sqrt{2}} \left( Y_l^{|m|} - (-1)^{|m|} Y_l^{-|m|} \right) (m<0)
24 \f]
25 
26 which implies:
27 \f[
28 Y_{l0}(\cos\theta,\phi) = N_{l0} P_l^0 (\cos\theta) (m=0)
29 Y_{lm}(\cos\theta,\phi) = \sqrt{2} N_{lm} P_l^m (\cos\theta) cos(|m|\phi) (m>0)
30 Y_{lm}(\cos\theta,\phi) = \sqrt{2} N_{l|m|} P_l^{|m|}(\cos\theta) sin(|m|\phi) (m<0)
31 \f]
32 
33 where
34 \f[
35  N_{lm} = \sqrt{ \frac{2l+1}{4\pi} \frac{ (l-m)! }{ (l+m)! } }
36 \f]
37 
38 Note that the normalization corresponds to the orthonormal case,
39 and thus we have \f$ \int d\cos\theta d\phi Y_{lm} Y_{l'm'} = \delta_{ll'} \delta{mm'}\f$
40 
41 Note that in addition, this code can also represent the product of two
42 (real) spherical harmonics -- it actually uses the fact that \f$Y_{00} = \sqrt{\frac{1}{4\pi}}\f$
43 in order to represent a single spherical harmonics by multiplying it
44 by \f$\sqrt{4\pi} Y_00\f$, as this makes it trivial to compute the analytical
45 integrals, using the orthogonality properties of \f$Y_l^m\f$...
46 
47 **/
48 
49 #include "RooFit.h"
50 #include "Riostream.h"
51 #include <math.h>
52 
53 #include "RooSpHarmonic.h"
54 #include "Math/SpecFunc.h"
55 #include "TMath.h"
56 
57 using namespace std;
58 
59 ClassImp(RooSpHarmonic);
60 
61 ////////////////////////////////////////////////////////////////////////////////
62 
63 namespace {
64  inline double N(int l, int m=0) {
65  double n = sqrt( double(2*l+1)/(4*TMath::Pi())*TMath::Factorial(l-m)/TMath::Factorial(l+m) );
66  return m==0 ? n : TMath::Sqrt2() * n;
67  }
68 }
69 
70 ////////////////////////////////////////////////////////////////////////////////
71 
72 RooSpHarmonic::RooSpHarmonic() :
73  _n(0),
74  _sgn1(0),
75  _sgn2(0)
76 {
77 }
78 
79 ////////////////////////////////////////////////////////////////////////////////
80 
81 RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l, int m)
82  : RooLegendre(name, title,ctheta,l,m<0?-m:m)
83  , _phi("phi", "phi", this, phi)
84  , _n( 2*sqrt(TMath::Pi()))
85  , _sgn1( m==0 ? 0 : m<0 ? -1 : +1 )
86  , _sgn2( 0 )
87 {
88 }
89 
90 ////////////////////////////////////////////////////////////////////////////////
91 
92 RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l1, int m1, int l2, int m2)
93  : RooLegendre(name, title,ctheta,l1, m1<0?-m1:m1,l2,m2<0?-m2:m2)
94  , _phi("phi", "phi", this, phi)
95  , _n(1)
96  , _sgn1( m1==0 ? 0 : m1<0 ? -1 : +1 )
97  , _sgn2( m2==0 ? 0 : m2<0 ? -1 : +1 )
98 {
99 }
100 
101 ////////////////////////////////////////////////////////////////////////////////
102 
103 RooSpHarmonic::RooSpHarmonic(const RooSpHarmonic& other, const char* name)
104  : RooLegendre(other, name)
105  , _phi("phi", this,other._phi)
106  , _n(other._n)
107  , _sgn1( other._sgn1 )
108  , _sgn2( other._sgn2 )
109 {
110 }
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 
114 Double_t RooSpHarmonic::evaluate() const
115 {
116  double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
117  if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
118  if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
119  return n;
120 }
121 
122 ////////////////////////////////////////////////////////////////////////////////
123 
124 namespace {
125  Bool_t fullRange(const RooRealProxy& x, const char* range, Bool_t phi)
126  {
127  if (phi) {
128  return range == 0 || strlen(range) == 0
129  ? std::fabs(x.max() - x.min() - TMath::TwoPi()) < 1.e-8
130  : std::fabs(x.max(range) - x.min(range) - TMath::TwoPi()) < 1.e-8;
131  }
132 
133  return range == 0 || strlen(range) == 0
134  ? std::fabs(x.min() + 1.) < 1.e-8 && std::fabs(x.max() - 1.) < 1.e-8
135  : std::fabs(x.min(range) + 1.) < 1.e-8 && std::fabs(x.max(range) - 1.) < 1.e-8;
136  }
137 }
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// TODO: check that phi.max - phi.min = 2 pi... ctheta.max = +1, and ctheta.min = -1
141 /// we don't support indefinite integrals... maybe one day, when there is a use for it.....
142 
143 Int_t RooSpHarmonic::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
144 {
145  // we don't support indefinite integrals... maybe one day, when there is a use for it.....
146  Bool_t cthetaOK = fullRange(_ctheta, rangeName, kFALSE);
147  Bool_t phiOK = fullRange(_phi, rangeName, kTRUE );
148  if (cthetaOK && phiOK && matchArgs(allVars, analVars, _ctheta, _phi)) return 3;
149  if ( phiOK && matchArgs(allVars, analVars, _phi)) return 2;
150  return RooLegendre::getAnalyticalIntegral(allVars, analVars, rangeName);
151 }
152 
153 ////////////////////////////////////////////////////////////////////////////////
154 
155 Double_t RooSpHarmonic::analyticalIntegral(Int_t code, const char* range) const
156 {
157  if (code==3) {
158  return (_l1==_l2 && _sgn1*_m1==_sgn2*_m2 ) ? _n : 0 ;
159  } else if (code == 2) {
160  if ( _sgn1*_m1 != _sgn2*_m2) return 0;
161  return ( _m1==0 ? 2 : 1 ) * TMath::Pi()*_n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
162  } else {
163  double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::analyticalIntegral(code,range);
164  if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
165  if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
166  return n;
167  }
168 }
169 
170 Int_t RooSpHarmonic::getMaxVal( const RooArgSet& vars) const {
171  return RooLegendre::getMaxVal(vars);
172 }
173 
174 Double_t RooSpHarmonic::maxVal( Int_t code) const {
175  double n = _n*N(_l1,_m1)*N(_l2,_m2);
176  return n*RooLegendre::maxVal(code);
177 }