Logo ROOT  
Reference Guide
 
Loading...
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 * *
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 <cmath>
51
52#include "RooSpHarmonic.h"
53#include "Math/SpecFunc.h"
54#include "TMath.h"
55
57
58////////////////////////////////////////////////////////////////////////////////
59
60namespace {
61 inline double N(int l, int m=0) {
62 double n = sqrt( double(2*l+1)/(4*TMath::Pi())*TMath::Factorial(l-m)/TMath::Factorial(l+m) );
63 return m==0 ? n : TMath::Sqrt2() * n;
64 }
65}
66
67////////////////////////////////////////////////////////////////////////////////
68
70 _n(0),
71 _sgn1(0),
72 _sgn2(0)
73{
74}
75
76////////////////////////////////////////////////////////////////////////////////
77
78RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l, int m)
79 : RooLegendre(name, title,ctheta,l,m<0?-m:m)
80 , _phi("phi", "phi", this, phi)
81 , _n( 2*sqrt(TMath::Pi()))
82 , _sgn1( m==0 ? 0 : m<0 ? -1 : +1 )
83 , _sgn2( 0 )
84{
85}
86
87////////////////////////////////////////////////////////////////////////////////
88
89RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l1, int m1, int l2, int m2)
90 : RooLegendre(name, title,ctheta,l1, m1<0?-m1:m1,l2,m2<0?-m2:m2)
91 , _phi("phi", "phi", this, phi)
92 , _n(1)
93 , _sgn1( m1==0 ? 0 : m1<0 ? -1 : +1 )
94 , _sgn2( m2==0 ? 0 : m2<0 ? -1 : +1 )
95{
96}
97
98////////////////////////////////////////////////////////////////////////////////
99
101 : RooLegendre(other, name)
102 , _phi("phi", this,other._phi)
103 , _n(other._n)
104 , _sgn1( other._sgn1 )
105 , _sgn2( other._sgn2 )
106{
107}
108
109////////////////////////////////////////////////////////////////////////////////
110
112{
113 double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
114 if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
115 if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
116 return n;
117}
118
119////////////////////////////////////////////////////////////////////////////////
120
121namespace {
122 bool fullRange(const RooRealProxy& x, const char* range, bool phi)
123 {
124 if (phi) {
125 return range == nullptr || strlen(range) == 0
126 ? std::abs(x.max() - x.min() - TMath::TwoPi()) < 1.e-8
127 : std::abs(x.max(range) - x.min(range) - TMath::TwoPi()) < 1.e-8;
128 }
129
130 return range == nullptr || strlen(range) == 0
131 ? std::abs(x.min() + 1.) < 1.e-8 && std::abs(x.max() - 1.) < 1.e-8
132 : std::abs(x.min(range) + 1.) < 1.e-8 && std::abs(x.max(range) - 1.) < 1.e-8;
133 }
134}
135
136////////////////////////////////////////////////////////////////////////////////
137/// TODO: check that phi.max - phi.min = 2 pi... ctheta.max = +1, and ctheta.min = -1
138/// we don't support indefinite integrals. maybe one day, when there is a use for it.
139
140Int_t RooSpHarmonic::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
141{
142 // we don't support indefinite integrals... maybe one day, when there is a use for it.....
143 bool cthetaOK = fullRange(_ctheta, rangeName, false);
144 bool phiOK = fullRange(_phi, rangeName, true );
145 if (cthetaOK && phiOK && matchArgs(allVars, analVars, _ctheta, _phi)) return 3;
146 if ( phiOK && matchArgs(allVars, analVars, _phi)) return 2;
147 return RooLegendre::getAnalyticalIntegral(allVars, analVars, rangeName);
148}
149
150////////////////////////////////////////////////////////////////////////////////
151
152double RooSpHarmonic::analyticalIntegral(Int_t code, const char* range) const
153{
154 if (code==3) {
155 return (_l1==_l2 && _sgn1*_m1==_sgn2*_m2 ) ? _n : 0 ;
156 } else if (code == 2) {
157 if ( _sgn1*_m1 != _sgn2*_m2) return 0;
158 return ( _m1==0 ? 2 : 1 ) * TMath::Pi()*_n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
159 } else {
160 double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::analyticalIntegral(code,range);
161 if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
162 if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
163 return n;
164 }
165}
166
168 return RooLegendre::getMaxVal(vars);
169}
170
171double RooSpHarmonic::maxVal( Int_t code) const {
172 double n = _n*N(_l1,_m1)*N(_l2,_m2);
173 return n*RooLegendre::maxVal(code);
174}
#define ClassImp(name)
Definition Rtypes.h:377
#define N
char name[80]
Definition TGX11.cxx:110
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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:55
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:39
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,...
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
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],...
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)
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