ROOT   6.10/09 Reference Guide
RooLegendre.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 RooLegendre
16  \ingroup Roofit
17
18 **/
19
20 #include "RooFit.h"
21 #include "Riostream.h"
22 #include <math.h>
23 #include <string>
24 #include <algorithm>
25
26 #include "RooLegendre.h"
27 #include "RooAbsReal.h"
28 #include "Math/SpecFunc.h"
29 #include "TMath.h"
30
31 #include "TError.h"
32
33 using namespace std;
34
36
37 ////////////////////////////////////////////////////////////////////////////////
38
39 namespace {
40  inline double a(int p, int l, int m) {
42  r /= pow(2.,m+2*p);
43  return p%2==0 ? r : -r ;
44  }
45 }
46
47 ////////////////////////////////////////////////////////////////////////////////
48
50  _l1(1),_m1(1),_l2(0),_m2(0)
51 {
52 }
53
54 ////////////////////////////////////////////////////////////////////////////////
55 ///TODO: for now, we assume that ctheta has a range [-1,1]
56 /// should map the ctheta range onto this interval, and adjust integrals...
57
58 RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l, int m)
59  : RooAbsReal(name, title)
60  , _ctheta("ctheta", "ctheta", this, ctheta)
61  , _l1(l),_m1(m),_l2(0),_m2(0)
62 {
63  //TODO: we assume m>=0
64  // should map m<0 back to m>=0...
65 }
66
67 ////////////////////////////////////////////////////////////////////////////////
68
69 RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l1, int m1, int l2, int m2)
70  : RooAbsReal(name, title)
71  , _ctheta("ctheta", "ctheta", this, ctheta)
72  , _l1(l1),_m1(m1),_l2(l2),_m2(m2)
73 {
74 }
75
76 ////////////////////////////////////////////////////////////////////////////////
77
78 RooLegendre::RooLegendre(const RooLegendre& other, const char* name)
79  : RooAbsReal(other, name)
80  , _ctheta("ctheta", this, other._ctheta)
81  , _l1(other._l1), _m1(other._m1)
82  , _l2(other._l2), _m2(other._m2)
83 {
84 }
85
86 ////////////////////////////////////////////////////////////////////////////////
87 /// TODO: check that 0<=m_i<=l_i; on the other hand, assoc_legendre already does that ;-)
88 /// Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0
89
91 {
92 #ifdef R__HAS_MATHMORE
93  double r = 1;
94  double ctheta = std::max(-1., std::min((double)_ctheta, +1.));
95  if (_l1!=0||_m1!=0) r *= ROOT::Math::assoc_legendre(_l1,_m1,ctheta);
96  if (_l2!=0||_m2!=0) r *= ROOT::Math::assoc_legendre(_l2,_m2,ctheta);
97  if ((_m1+_m2)%2==1) r = -r;
98  return r;
99 #else
100  throw std::string("RooLegendre: ERROR: This class require installation of the MathMore library") ;
101  return 0 ;
102 #endif
103 }
104
105 ////////////////////////////////////////////////////////////////////////////////
106
107 namespace {
108  Bool_t fullRange(const RooRealProxy& x ,const char* range)
109  {
110  return range == 0 || strlen(range) == 0
111  ? std::fabs(x.min() + 1.) < 1.e-8 && std::fabs(x.max() - 1.) < 1.e-8
112  : std::fabs(x.min(range) + 1.) < 1.e-8 && std::fabs(x.max(range) - 1.) < 1.e-8;
113  }
114 }
115
116 ////////////////////////////////////////////////////////////////////////////////
117
118 Int_t RooLegendre::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
119 {
120  // don't support indefinite integrals...
121  if (fullRange(_ctheta,rangeName) && matchArgs(allVars, analVars, _ctheta)) return 1;
122  return 0;
123 }
124
125 ////////////////////////////////////////////////////////////////////////////////
126 /// this was verified to match mathematica for
127 /// l1 in [0,2], m1 in [0,l1], l2 in [l1,4], m2 in [0,l2]
128
130 {
131  R__ASSERT(code==1) ;
132  if ( _m1==_m2 ) return ( _l1 == _l2) ? TMath::Factorial(_l1+_m2)/TMath::Factorial(_l1-_m1)*double(2)/(2*_l1+1) : 0.;
133  if ( (_l1+_l2-_m1-_m2)%2 != 0 ) return 0; // these combinations are odd under x -> -x
134
135  // from B.R. Wong, "On the overlap integral of associated Legendre Polynomials" 1998 J. Phys. A: Math. Gen. 31 1101
136  // TODO: update to the result of
137  // H. A. Mavromatis
138  // "A single-sum expression for the overlap integral of two associated Legendre polynomials"
139  // 1999 J. Phys. A: Math. Gen. 32 2601
140  // http://iopscience.iop.org/0305-4470/32/13/011/pdf/0305-4470_32_13_011.pdf
141  // For that we need Wigner 3-j, which Lorenzo has added for Root 5.28... (note: check Condon-Shortly convention in this paper!)
142  double r=0;
143  for (int p1=0; 2*p1 <= _l1-_m1 ;++p1) {
144  double a1 = a(p1,_l1,_m1);
145  for (int p2=0; 2*p2 <= _l2-_m2 ; ++p2) {
146  double a2 = a(p2,_l2,_m2);
147  r+= a1*a2*TMath::Gamma( double(_l1+_l2-_m1-_m2-2*p1-2*p2+1)/2 )*TMath::Gamma( double(_m1+_m2+2*p1+2*p2+2)/2 );
148  }
149  }
150  r /= TMath::Gamma( double(_l1+_l2+3)/2 );
151
152  if ((_m1+_m2)%2==1) r = -r;
153  return r;
154 }
155
156 ////////////////////////////////////////////////////////////////////////////////
157
158 Int_t RooLegendre::getMaxVal( const RooArgSet& /*vars*/) const {
159  if (_m1==0&&_m2==0) return 1;
160  // does anyone know the analytical expression for the max values in case m!=0??
161  if (_l1<3&&_l2<3) return 1;
162  return 0;
163 }
164
165 namespace {
166  inline double maxSingle(int i, int j) {
167  R__ASSERT(j<=i);
168  // x0 : 1 (ordinary Legendre)
169  if (j==0) return 1;
170  R__ASSERT(i<3);
171  // 11: 1
172  if (i<2) return 1;
173  // 21: 3 22: 3
174  static const double m2[3] = { 3,3 };
175  return m2[j-1];
176  }
177 }
179  return maxSingle(_l1,_m1)*maxSingle(_l2,_m2);
180 }
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
#define R__ASSERT(e)
Definition: TError.h:96
RooRealProxy _ctheta
Definition: RooLegendre.h:40
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:352
STL namespace.
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]
static double p2(double t, double a, double b, double c)
double pow(double, double)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
TRandom2 r(17)
TMarker * m
Definition: textangle.C:8
double assoc_legendre(unsigned l, unsigned m, double x)
Computes the associated Legendre polynomials.
TLine * l
Definition: textangle.C:4
static double p1(double t, double a, double b)
#define ClassImp(name)
Definition: Rtypes.h:336
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
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
RooRealProxy is the concrete proxy for RooAbsReal objects A RooRealProxy is the general mechanism to ...
Definition: RooRealProxy.h:23
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:90
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.