Logo ROOT   6.08/07
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 /**
16 \file RooSpHarmonic.cxx
17 \class RooSpHarmonic
18 \ingroup Roofit
19 
20  Implementation of the so-called real spherical harmonics, using the orthonormal normalization,
21 which are related to spherical harmonics as:
22 \f[
23  Y_{l0} = Y_l^0 (m=0) \\
24  Y_{lm} = \frac{1}{\sqrt{2}} \left( Y_l^m + (-1)^m Y_l^{-m} \right) (m>0) \\
25  Y_{lm} = \frac{1}{i\sqrt{2}} \left( Y_l^{|m|} - (-1)^{|m|} Y_l^{-|m|} \right) (m<0)
26 \f]
27 
28 which implies:
29 \f[
30 Y_{l0}(\cos\theta,\phi) = N_{l0} P_l^0 (\cos\theta) (m=0)
31 Y_{lm}(\cos\theta,\phi) = \sqrt{2} N_{lm} P_l^m (\cos\theta) cos(|m|\phi) (m>0)
32 Y_{lm}(\cos\theta,\phi) = \sqrt{2} N_{l|m|} P_l^{|m|}(\cos\theta) sin(|m|\phi) (m<0)
33 \f]
34 
35 where
36 \f[
37  N_{lm} = \sqrt{ \frac{2l+1}{4\pi} \frac{ (l-m)! }{ (l+m)! } }
38 \f]
39 
40 Note that the normalization corresponds to the orthonormal case,
41 and thus we have \f$ \int d\cos\theta d\phi Y_{lm} Y_{l'm'} = \delta_{ll'} \delta{mm'}\f$
42 
43 Note that in addition, this code can also represent the product of two
44 (real) spherical harmonics -- it actually uses the fact that \f$Y_{00} = \sqrt{\frac{1}{4\pi}}\f$
45 in order to represent a single spherical harmonics by multiplying it
46 by \f$\sqrt{4\pi} Y_00\f$, as this makes it trivial to compute the analytical
47 integrals, using the orthogonality properties of \f$Y_l^m\f$...
48 
49 **/
50 
51 #include "RooFit.h"
52 #include "Riostream.h"
53 #include <math.h>
54 
55 #include "RooSpHarmonic.h"
56 #include "Math/SpecFunc.h"
57 #include "TMath.h"
58 
59 using namespace std;
60 
62 ;
63 
64 
65 ////////////////////////////////////////////////////////////////////////////////
66 
67 namespace {
68  inline double N(int l, int m=0) {
69  double n = sqrt( double(2*l+1)/(4*TMath::Pi())*TMath::Factorial(l-m)/TMath::Factorial(l+m) );
70  return m==0 ? n : TMath::Sqrt2() * n;
71  }
72 }
73 
74 ////////////////////////////////////////////////////////////////////////////////
75 
77  _n(0),
78  _sgn1(0),
79  _sgn2(0)
80 {
81 }
82 
83 ////////////////////////////////////////////////////////////////////////////////
84 
85 RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l, int m)
86  : RooLegendre(name, title,ctheta,l,m<0?-m:m)
87  , _phi("phi", "phi", this, phi)
88  , _n( 2*sqrt(TMath::Pi()))
89  , _sgn1( m==0 ? 0 : m<0 ? -1 : +1 )
90  , _sgn2( 0 )
91 {
92 }
93 
94 ////////////////////////////////////////////////////////////////////////////////
95 
96 RooSpHarmonic::RooSpHarmonic(const char* name, const char* title, RooAbsReal& ctheta, RooAbsReal& phi, int l1, int m1, int l2, int m2)
97  : RooLegendre(name, title,ctheta,l1, m1<0?-m1:m1,l2,m2<0?-m2:m2)
98  , _phi("phi", "phi", this, phi)
99  , _n(1)
100  , _sgn1( m1==0 ? 0 : m1<0 ? -1 : +1 )
101  , _sgn2( m2==0 ? 0 : m2<0 ? -1 : +1 )
102 {
103 }
104 
105 ////////////////////////////////////////////////////////////////////////////////
106 
108  : RooLegendre(other, name)
109  , _phi("phi", this,other._phi)
110  , _n(other._n)
111  , _sgn1( other._sgn1 )
112  , _sgn2( other._sgn2 )
113 {
114 }
115 
116 ////////////////////////////////////////////////////////////////////////////////
117 
119 {
120  double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
121  if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
122  if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
123  return n;
124 }
125 
126 //_____________________________________________________________________________
127 namespace {
128  Bool_t fullRange(const RooRealProxy& x, const char* range, Bool_t phi)
129  {
130  if (phi) {
131  return range == 0 || strlen(range) == 0
132  ? std::fabs(x.max() - x.min() - TMath::TwoPi()) < 1.e-8
133  : std::fabs(x.max(range) - x.min(range) - TMath::TwoPi()) < 1.e-8;
134  }
135 
136  return range == 0 || strlen(range) == 0
137  ? std::fabs(x.min() + 1.) < 1.e-8 && std::fabs(x.max() - 1.) < 1.e-8
138  : std::fabs(x.min(range) + 1.) < 1.e-8 && std::fabs(x.max(range) - 1.) < 1.e-8;
139  }
140 }
141 
142 ////////////////////////////////////////////////////////////////////////////////
143 /// TODO: check that phi.max - phi.min = 2 pi... ctheta.max = +1, and ctheta.min = -1
144 /// we don't support indefinite integrals... maybe one day, when there is a use for it.....
145 
146 Int_t RooSpHarmonic::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
147 {
148  // we don't support indefinite integrals... maybe one day, when there is a use for it.....
149  Bool_t cthetaOK = fullRange(_ctheta, rangeName, kFALSE);
150  Bool_t phiOK = fullRange(_phi, rangeName, kTRUE );
151  if (cthetaOK && phiOK && matchArgs(allVars, analVars, _ctheta, _phi)) return 3;
152  if ( phiOK && matchArgs(allVars, analVars, _phi)) return 2;
153  return RooLegendre::getAnalyticalIntegral(allVars, analVars, rangeName);
154 }
155 
156 ////////////////////////////////////////////////////////////////////////////////
157 
158 Double_t RooSpHarmonic::analyticalIntegral(Int_t code, const char* range) const
159 {
160  if (code==3) {
161  return (_l1==_l2 && _sgn1*_m1==_sgn2*_m2 ) ? _n : 0 ;
162  } else if (code == 2) {
163  if ( _sgn1*_m1 != _sgn2*_m2) return 0;
164  return ( _m1==0 ? 2 : 1 ) * TMath::Pi()*_n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::evaluate();
165  } else {
166  double n = _n*N(_l1,_m1)*N(_l2,_m2)*RooLegendre::analyticalIntegral(code,range);
167  if (_sgn1!=0) n *= (_sgn1<0 ? sin(_m1*_phi) : cos(_m1*_phi) );
168  if (_sgn2!=0) n *= (_sgn2<0 ? sin(_m2*_phi) : cos(_m2*_phi) );
169  return n;
170  }
171 }
172 
174  return RooLegendre::getMaxVal(vars);
175 }
176 
178  double n = _n*N(_l1,_m1)*N(_l2,_m2);
179  return n*RooLegendre::maxVal(code);
180 }
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise capability to determine maximum value of function for given set of observables.
Double_t Sqrt2()
Definition: TMath.h:51
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
#define N
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], m2 in [0,l2]
RooRealProxy _ctheta
Definition: RooLegendre.h:40
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
STL namespace.
double cos(double)
Double_t evaluate() const
TODO: check that 0<=m_i<=l_i; on the other hand, assoc_legendre already does that ;-) Note: P_0^0 = 1...
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
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], m2 in [0,l2]
RooRealProxy _phi
Definition: RooSpHarmonic.h:37
Double_t TwoPi()
Definition: TMath.h:45
double sin(double)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Implementation of the so-called real spherical harmonics, using the orthonormal normalization, which are related to spherical harmonics as: .
Definition: RooSpHarmonic.h:20
double Pi()
Mathematical constants.
Definition: Math.h:68
TMarker * m
Definition: textangle.C:8
TLine * l
Definition: textangle.C:4
Double_t Pi()
Definition: TMath.h:44
#define ClassImp(name)
Definition: Rtypes.h:279
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise capability to determine maximum value of function for given set of observables.
Double_t Factorial(Int_t i)
Compute factorial(n).
Definition: TMath.cxx:250
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
RooRealProxy is the concrete proxy for RooAbsReal objects A RooRealProxy is the general mechanism to ...
Definition: RooRealProxy.h:23
const Bool_t kTRUE
Definition: Rtypes.h:91
Double_t evaluate() const
TODO: check that 0<=m_i<=l_i; on the other hand, assoc_legendre already does that ;-) Note: P_0^0 = 1...
Definition: RooLegendre.cxx:94
const Int_t n
Definition: legend1.C:16
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
char name[80]
Definition: TGX11.cxx:109
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
TODO: check that phi.max - phi.min = 2 pi...