Logo ROOT  
Reference Guide
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,
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 "Riostream.h"
50#include <math.h>
51
52#include "RooSpHarmonic.h"
53#include "Math/SpecFunc.h"
54#include "TMath.h"
55
56using namespace std;
57
59
60////////////////////////////////////////////////////////////////////////////////
61
62namespace {
63 inline double N(int l, int m=0) {
64 double n = sqrt( double(2*l+1)/(4*TMath::Pi())*TMath::Factorial(l-m)/TMath::Factorial(l+m) );
65 return m==0 ? n : TMath::Sqrt2() * n;
66 }
67}
68
69////////////////////////////////////////////////////////////////////////////////
70
72 _n(0),
73 _sgn1(0),
74 _sgn2(0)
75{
76}
77
78////////////////////////////////////////////////////////////////////////////////
79
80RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l, int m)
81 : RooLegendre(name, title,ctheta,l,m<0?-m:m)
82 , _phi("phi", "phi", this, phi)
83 , _n( 2*sqrt(TMath::Pi()))
84 , _sgn1( m==0 ? 0 : m<0 ? -1 : +1 )
85 , _sgn2( 0 )
86{
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
91RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l1, int m1, int l2, int m2)
92 : RooLegendre(name, title,ctheta,l1, m1<0?-m1:m1,l2,m2<0?-m2:m2)
93 , _phi("phi", "phi", this, phi)
94 , _n(1)
95 , _sgn1( m1==0 ? 0 : m1<0 ? -1 : +1 )
96 , _sgn2( m2==0 ? 0 : m2<0 ? -1 : +1 )
97{
98}
99
100////////////////////////////////////////////////////////////////////////////////
101
103 : RooLegendre(other, name)
104 , _phi("phi", this,other._phi)
105 , _n(other._n)
106 , _sgn1( other._sgn1 )
107 , _sgn2( other._sgn2 )
108{
109}
110
111////////////////////////////////////////////////////////////////////////////////
112
114{
115 double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
116 if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
117 if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
118 return n;
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
123namespace {
124 bool fullRange(const RooRealProxy& x, const char* range, bool phi)
125 {
126 if (phi) {
127 return range == 0 || strlen(range) == 0
128 ? std::fabs(x.max() - x.min() - TMath::TwoPi()) < 1.e-8
129 : std::fabs(x.max(range) - x.min(range) - TMath::TwoPi()) < 1.e-8;
130 }
131
132 return range == 0 || strlen(range) == 0
133 ? std::fabs(x.min() + 1.) < 1.e-8 && std::fabs(x.max() - 1.) < 1.e-8
134 : std::fabs(x.min(range) + 1.) < 1.e-8 && std::fabs(x.max(range) - 1.) < 1.e-8;
135 }
136}
137
138////////////////////////////////////////////////////////////////////////////////
139/// TODO: check that phi.max - phi.min = 2 pi... ctheta.max = +1, and ctheta.min = -1
140/// we don't support indefinite integrals. maybe one day, when there is a use for it.
141
142Int_t RooSpHarmonic::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
143{
144 // we don't support indefinite integrals... maybe one day, when there is a use for it.....
145 bool cthetaOK = fullRange(_ctheta, rangeName, false);
146 bool phiOK = fullRange(_phi, rangeName, true );
147 if (cthetaOK && phiOK && matchArgs(allVars, analVars, _ctheta, _phi)) return 3;
148 if ( phiOK && matchArgs(allVars, analVars, _phi)) return 2;
149 return RooLegendre::getAnalyticalIntegral(allVars, analVars, rangeName);
150}
151
152////////////////////////////////////////////////////////////////////////////////
153
154double RooSpHarmonic::analyticalIntegral(Int_t code, const char* range) const
155{
156 if (code==3) {
157 return (_l1==_l2 && _sgn1*_m1==_sgn2*_m2 ) ? _n : 0 ;
158 } else if (code == 2) {
159 if ( _sgn1*_m1 != _sgn2*_m2) return 0;
160 return ( _m1==0 ? 2 : 1 ) * TMath::Pi()*_n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
161 } else {
162 double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::analyticalIntegral(code,range);
163 if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
164 if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
165 return n;
166 }
167}
168
170 return RooLegendre::getMaxVal(vars);
171}
172
173double RooSpHarmonic::maxVal( Int_t code) const {
174 double n = _n*N(_l1,_m1)*N(_l2,_m2);
175 return n*RooLegendre::maxVal(code);
176}
ClassImp(RooSpHarmonic)
#define N
char name[80]
Definition: TGX11.cxx:110
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
bool 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:56
Compute the associated Legendre polynomials using ROOT::Math::assoc_legendre().
Definition: RooLegendre.h:20
double evaluate() const override
Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0.
Int_t getMaxVal(const RooArgSet &vars) const override
Advertise capability to determine maximum value of function for given set of observables.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
this was verified to match mathematica for l1 in [0,2], m1 in [0,l1], l2 in [l1,4],...
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooRealProxy _ctheta
Definition: RooLegendre.h:40
double maxVal(Int_t code) const override
Return maximum value for set of observables identified by code assigned in getMaxVal.
Implementation of the so-called real spherical harmonics, using the orthonormal normalization,...
Definition: RooSpHarmonic.h:20
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
TODO: check that phi.max - phi.min = 2 pi... ctheta.max = +1, and ctheta.min = -1 we don't support in...
RooRealProxy _phi
Definition: RooSpHarmonic.h:37
Int_t getMaxVal(const RooArgSet &vars) const override
Advertise capability to determine maximum value of function for given set of observables.
double maxVal(Int_t code) const override
Return maximum value for set of observables identified by code assigned in getMaxVal.
double evaluate() const override
Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
this was verified to match mathematica for l1 in [0,2], m1 in [0,l1], l2 in [l1,4],...
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition: RVec.hxx:1776
RVec< PromoteType< T > > sin(const RVec< T > &v)
Definition: RVec.hxx:1775
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
double Pi()
Mathematical constants.
Definition: Math.h:88
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static constexpr double m2
TMath.
Definition: TMathBase.h:35
Double_t Factorial(Int_t i)
Computes factorial(n).
Definition: TMath.cxx:252
constexpr Double_t Sqrt2()
Definition: TMath.h:86
constexpr Double_t Pi()
Definition: TMath.h:37
constexpr Double_t TwoPi()
Definition: TMath.h:44
TMarker m
Definition: textangle.C:8
TLine l
Definition: textangle.C:4