ROOT   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 Compute the associated Legendre polynomials using ROOT::Math::assoc_legendre().
19
20 Since the Legendre polynomials have a value range of [-1, 1], these cannot be implemented as a PDF.
21 They can be used in sums, though, for example using a RooRealSumFunc of RooLegendre plus an offset.
22**/
23
24#include "RooLegendre.h"
25
26#include "RooFit.h"
27
28#include "RooAbsReal.h"
29#include "Math/SpecFunc.h"
30#include "TMath.h"
31
32#include <cmath>
33#include <string>
34#include <algorithm>
35
36using namespace std;
37
39
40////////////////////////////////////////////////////////////////////////////////
41
42namespace {
43 inline double a(int p, int l, int m) {
45 r /= pow(2.,m+2*p);
46 return p%2==0 ? r : -r ;
47 }
48
49 void throwIfNoMathMore() {
50#ifndef R__HAS_MATHMORE
51 throw std::runtime_error("RooLegendre needs functions from MathMore. It is not available in this root build.");
52#endif
53 }
54
55 void checkCoeffs(int m1, int l1, int m2, int l2) {
56 if (m1 < 0 || m2 < 0) {
57 throw std::invalid_argument("RooLegendre: m coefficients need to be >= 0.");
58 }
59 if (l1 < m1 || l2 < m2) {
60 throw std::invalid_argument("RooLegendre: m coefficients need to be smaller than corresponding l.");
61 }
62 }
63}
64
65////////////////////////////////////////////////////////////////////////////////
66
68 _l1(1),_m1(1),_l2(0),_m2(0)
69{
70 throwIfNoMathMore();
71}
72
73////////////////////////////////////////////////////////////////////////////////
74///TODO: for now, we assume that ctheta has a range [-1,1]
75/// should map the ctheta range onto this interval, and adjust integrals...
76
77RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l, int m)
78 : RooAbsReal(name, title)
79 , _ctheta("ctheta", "ctheta", this, ctheta)
80 , _l1(l),_m1(m),_l2(0),_m2(0)
81{
82 checkCoeffs(_m1, _l1, _m2, _l2);
83
84 throwIfNoMathMore();
85}
86
87////////////////////////////////////////////////////////////////////////////////
88
89RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l1, int m1, int l2, int m2)
90 : RooAbsReal(name, title)
91 , _ctheta("ctheta", "ctheta", this, ctheta)
92 , _l1(l1),_m1(m1),_l2(l2),_m2(m2)
93{
94 checkCoeffs(_m1, _l1, _m2, _l2);
95
96 throwIfNoMathMore();
97}
98
99////////////////////////////////////////////////////////////////////////////////
100
101RooLegendre::RooLegendre(const RooLegendre& other, const char* name)
102 : RooAbsReal(other, name)
103 , _ctheta("ctheta", this, other._ctheta)
104 , _l1(other._l1), _m1(other._m1)
105 , _l2(other._l2), _m2(other._m2)
106{
107}
108
109////////////////////////////////////////////////////////////////////////////////
110/// Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0
111
113{
114#ifdef R__HAS_MATHMORE
115 double r = 1;
116 double ctheta = std::max(-1., std::min((double)_ctheta, +1.));
117 if (_l1!=0||_m1!=0) r *= ROOT::Math::assoc_legendre(_l1,_m1,ctheta);
118 if (_l2!=0||_m2!=0) r *= ROOT::Math::assoc_legendre(_l2,_m2,ctheta);
119 if ((_m1+_m2)%2==1) r = -r;
120 return r;
121#else
122 throwIfNoMathMore();
123 return 0.;
124#endif
125}
126
127
128////////////////////////////////////////////////////////////////////////////////
129
130namespace {
131//Author: Emmanouil Michalainas, CERN 26 August 2019
132
133void compute( size_t batchSize, const int l1, const int m1, const int l2, const int m2,
134 double * __restrict output,
135 double const * __restrict TH)
136{
137#ifdef R__HAS_MATHMORE
138 double legendre1=1.0, legendreMinus1=1.0;
139 if (l1+m1 > 0) {
140 legendre1 = ROOT::Math::internal::legendre(l1,m1,1.0);
141 legendreMinus1 = ROOT::Math::internal::legendre(l1,m1,-1.0);
142 }
143 if (l2+m2 > 0) {
144 legendre1 *= ROOT::Math::internal::legendre(l2,m2,1.0);
145 legendreMinus1 *= ROOT::Math::internal::legendre(l2,m2,-1.0);
146 }
147
148 for (size_t i=0; i<batchSize; i++) {
149 if (TH[i] <= -1.0) {
150 output[i] = legendreMinus1;
151 } else if (TH[i] >= 1.0) {
152 output[i] = legendre1;
153 }
154 else {
155 output[i] = 1.0;
156 if (l1+m1 > 0) {
157 output[i] *= ROOT::Math::internal::legendre(l1,m1,TH[i]);
158 }
159 if (l2+m2 > 0) {
161 }
162 }
163 }
164
165#else
166 (void) batchSize, (void) l1, (void)m1, (void)l2, (void)m2, (void)output, (void)TH;
167 throwIfNoMathMore();
168#endif
169}
170};
171
172RooSpan<double> RooLegendre::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
173 auto cthetaData = _ctheta.getValBatch(begin, batchSize);
174
176 return {};
177 }
178
180 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
181
182 compute(batchSize, _l1, _m1, _l2, _m2, output.data(), cthetaData.data());
183
184 return output;
185}
186
187
188////////////////////////////////////////////////////////////////////////////////
189
190namespace {
191 Bool_t fullRange(const RooRealProxy& x ,const char* range)
192 {
193 return range == 0 || strlen(range) == 0
194 ? std::fabs(x.min() + 1.) < 1.e-8 && std::fabs(x.max() - 1.) < 1.e-8
195 : std::fabs(x.min(range) + 1.) < 1.e-8 && std::fabs(x.max(range) - 1.) < 1.e-8;
196 }
197}
198
199////////////////////////////////////////////////////////////////////////////////
200
201Int_t RooLegendre::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
202{
203 // don't support indefinite integrals...
204 if (fullRange(_ctheta,rangeName) && matchArgs(allVars, analVars, _ctheta)) return 1;
205 return 0;
206}
207
208////////////////////////////////////////////////////////////////////////////////
209/// this was verified to match mathematica for
210/// l1 in [0,2], m1 in [0,l1], l2 in [l1,4], m2 in [0,l2]
211
213{
214 R__ASSERT(code==1) ;
215 if ( _m1==_m2 ) return ( _l1 == _l2) ? TMath::Factorial(_l1+_m2)/TMath::Factorial(_l1-_m1)*double(2)/(2*_l1+1) : 0.;
216 if ( (_l1+_l2-_m1-_m2)%2 != 0 ) return 0; // these combinations are odd under x -> -x
217
218 // from B.R. Wong, "On the overlap integral of associated Legendre Polynomials" 1998 J. Phys. A: Math. Gen. 31 1101
219 // TODO: update to the result of
220 // H. A. Mavromatis
221 // "A single-sum expression for the overlap integral of two associated Legendre polynomials"
222 // 1999 J. Phys. A: Math. Gen. 32 2601
223 // http://iopscience.iop.org/0305-4470/32/13/011/pdf/0305-4470_32_13_011.pdf
224 // For that we need Wigner 3-j, which Lorenzo has added for Root 5.28... (note: check Condon-Shortly convention in this paper!)
225 double r=0;
226 for (int p1=0; 2*p1 <= _l1-_m1 ;++p1) {
227 double a1 = a(p1,_l1,_m1);
228 for (int p2=0; 2*p2 <= _l2-_m2 ; ++p2) {
229 double a2 = a(p2,_l2,_m2);
230 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 );
231 }
232 }
233 r /= TMath::Gamma( double(_l1+_l2+3)/2 );
234
235 if ((_m1+_m2)%2==1) r = -r;
236 return r;
237}
238
239////////////////////////////////////////////////////////////////////////////////
240
241Int_t RooLegendre::getMaxVal( const RooArgSet& /*vars*/) const {
242 if (_m1==0&&_m2==0) return 1;
243 // does anyone know the analytical expression for the max values in case m!=0??
244 if (_l1<3&&_l2<3) return 1;
245 return 0;
246}
247
248namespace {
249 inline double maxSingle(int i, int j) {
250 R__ASSERT(j<=i);
251 // x0 : 1 (ordinary Legendre)
252 if (j==0) return 1;
253 R__ASSERT(i<3);
254 // 11: 1
255 if (i<2) return 1;
256 // 21: 3 22: 3
257 static const double m2[3] = { 3,3 };
258 return m2[j-1];
259 }
260}
262 return maxSingle(_l1,_m1)*maxSingle(_l2,_m2);
263}
double
Definition: Converters.cxx:921
ROOT::R::TRInterface & r
Definition: Object.C:4
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
double pow(double, double)
typedef void((*Func_t)())
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:118
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:450
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Compute the associated Legendre polynomials using ROOT::Math::assoc_legendre().
Definition: RooLegendre.h:20
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
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
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
double assoc_legendre(unsigned l, unsigned m, double x)
Computes the associated Legendre polynomials.
Double_t x[n]
Definition: legend1.C:17
double legendre(unsigned l, unsigned m, double x)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static constexpr double m2
Double_t Factorial(Int_t i)
Compute factorial(n).
Definition: TMath.cxx:247
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:348
auto * m
Definition: textangle.C:8
auto * l
Definition: textangle.C:4
auto * a
Definition: textangle.C:12
static void output(int code)
Definition: gifencode.c:226