ROOT   Reference Guide
Searching...
No Matches
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 * *
9 * *
10 * Redistribution and use in source and binary forms, *
11 * with or without modification, are permitted according to the terms *
13 *****************************************************************************/
14
15/** \class RooSpHarmonic
16 \ingroup Roofit
17
18 Implementation of the so-called real spherical harmonics, using the orthonormal normalization,
19which 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
26which implies:
27\f[
28Y_{l0}(\cos\theta,\phi) = N_{l0} P_l^0 (\cos\theta) (m=0)
29Y_{lm}(\cos\theta,\phi) = \sqrt{2} N_{lm} P_l^m (\cos\theta) cos(|m|\phi) (m>0)
30Y_{lm}(\cos\theta,\phi) = \sqrt{2} N_{l|m|} P_l^{|m|}(\cos\theta) sin(|m|\phi) (m<0)
31\f]
32
33where
34\f[
35 N_{lm} = \sqrt{ \frac{2l+1}{4\pi} \frac{ (l-m)! }{ (l+m)! } }
36\f]
37
38Note that the normalization corresponds to the orthonormal case,
39and thus we have \f$\int d\cos\theta d\phi Y_{lm} Y_{l'm'} = \delta_{ll'} \delta{mm'}\f$
40
41Note 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$
43in order to represent a single spherical harmonics by multiplying it
44by \f$\sqrt{4\pi} Y_00\f$, as this makes it trivial to compute the analytical
45integrals, 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
57using namespace std;
58
60
61////////////////////////////////////////////////////////////////////////////////
62
63namespace {
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
73 _n(0),
74 _sgn1(0),
75 _sgn2(0)
76{
77}
78
79////////////////////////////////////////////////////////////////////////////////
80
81RooSpHarmonic::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
92RooSpHarmonic::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
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
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
124namespace {
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
143Int_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
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
171 return RooLegendre::getMaxVal(vars);
172}
173
175 double n = _n*N(_l1,_m1)*N(_l2,_m2);
176 return n*RooLegendre::maxVal(code);
177}
const Bool_t kFALSE
Definition RtypesCore.h:92
const Bool_t kTRUE
Definition RtypesCore.h:91
#define ClassImp(name)
Definition Rtypes.h:364
#define N
char name[80]
Definition TGX11.cxx:110
double cos(double)
double sqrt(double)
double sin(double)
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:61
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
Compute the associated Legendre polynomials using ROOT::Math::assoc_legendre().
Definition RooLegendre.h:20
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Double_t evaluate() const
Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0.
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
this was verified to match mathematica for l1 in [0,2], m1 in [0,l1], l2 in [l1,4],...
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise capability to determine maximum value of function for given set of observables.
RooRealProxy _ctheta
Definition RooLegendre.h:40
Implementation of the so-called real spherical harmonics, using the orthonormal normalization,...
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
TODO: check that phi.max - phi.min = 2 pi... ctheta.max = +1, and ctheta.min = -1 we don't support in...
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
this was verified to match mathematica for l1 in [0,2], m1 in [0,l1], l2 in [l1,4],...
RooRealProxy _phi
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
Double_t evaluate() const
Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0.
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise capability to determine maximum value of function for given set of observables.
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TMath.
Definition TMathBase.h:35
Double_t Factorial(Int_t i)
Compute factorial(n).
Definition TMath.cxx:247
constexpr Double_t Sqrt2()
Definition TMath.h:88
constexpr Double_t Pi()
Definition TMath.h:37
constexpr Double_t TwoPi()
Definition TMath.h:44
auto * m
Definition textangle.C:8
auto * l
Definition textangle.C:4