Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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 * *
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 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#include "RooAbsReal.h"
26
27#include "Math/SpecFunc.h"
28#include "TMath.h"
29
30#include <cmath>
31#include <string>
32#include <algorithm>
33
34using namespace std;
35
37
38////////////////////////////////////////////////////////////////////////////////
39
40namespace {
41 inline double a(int p, int l, int m) {
43 r /= pow(2.,m+2*p);
44 return p%2==0 ? r : -r ;
45 }
46
47 void checkCoeffs(int m1, int l1, int m2, int l2) {
48 if (m1 < 0 || m2 < 0) {
49 throw std::invalid_argument("RooLegendre: m coefficients need to be >= 0.");
50 }
51 if (l1 < m1 || l2 < m2) {
52 throw std::invalid_argument("RooLegendre: m coefficients need to be smaller than corresponding l.");
53 }
54 }
55}
56
57////////////////////////////////////////////////////////////////////////////////
58
60 _l1(1),_m1(1),_l2(0),_m2(0)
61{
62}
63
64////////////////////////////////////////////////////////////////////////////////
65///TODO: for now, we assume that ctheta has a range [-1,1]
66/// should map the ctheta range onto this interval, and adjust integrals...
67
68RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l, int m)
69 : RooAbsReal(name, title)
70 , _ctheta("ctheta", "ctheta", this, ctheta)
71 , _l1(l),_m1(m),_l2(0),_m2(0)
72{
73 checkCoeffs(_m1, _l1, _m2, _l2);
74}
75
76////////////////////////////////////////////////////////////////////////////////
77
78RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l1, int m1, int l2, int m2)
79 : RooAbsReal(name, title)
80 , _ctheta("ctheta", "ctheta", this, ctheta)
81 , _l1(l1),_m1(m1),_l2(l2),_m2(m2)
82{
83 checkCoeffs(_m1, _l1, _m2, _l2);
84}
85
86////////////////////////////////////////////////////////////////////////////////
87
88RooLegendre::RooLegendre(const RooLegendre& other, const char* name)
89 : RooAbsReal(other, name)
90 , _ctheta("ctheta", this, other._ctheta)
91 , _l1(other._l1), _m1(other._m1)
92 , _l2(other._l2), _m2(other._m2)
93{
94}
95
96////////////////////////////////////////////////////////////////////////////////
97/// Note: P_0^0 = 1, so P_l^m = P_l^m P_0^0
98
100{
101 double r = 1;
102 double ctheta = std::max(-1., std::min((double)_ctheta, +1.));
103 if (_l1!=0||_m1!=0) r *= ROOT::Math::assoc_legendre(_l1,_m1,ctheta);
104 if (_l2!=0||_m2!=0) r *= ROOT::Math::assoc_legendre(_l2,_m2,ctheta);
105 if ((_m1+_m2)%2==1) r = -r;
106 return r;
107}
108
109
110////////////////////////////////////////////////////////////////////////////////
111
112namespace {
113//Author: Emmanouil Michalainas, CERN 26 August 2019
114
115void compute( size_t batchSize, const int l1, const int m1, const int l2, const int m2,
116 double * __restrict output,
117 double const * __restrict TH)
118{
119 double legendre1=1.0, legendreMinus1=1.0;
120 if (l1+m1 > 0) {
121 legendre1 = ROOT::Math::internal::legendre(l1,m1,1.0);
122 legendreMinus1 = ROOT::Math::internal::legendre(l1,m1,-1.0);
123 }
124 if (l2+m2 > 0) {
125 legendre1 *= ROOT::Math::internal::legendre(l2,m2,1.0);
126 legendreMinus1 *= ROOT::Math::internal::legendre(l2,m2,-1.0);
127 }
128
129 for (size_t i=0; i<batchSize; i++) {
130 if (TH[i] <= -1.0) {
131 output[i] = legendreMinus1;
132 } else if (TH[i] >= 1.0) {
133 output[i] = legendre1;
134 }
135 else {
136 output[i] = 1.0;
137 if (l1+m1 > 0) {
138 output[i] *= ROOT::Math::internal::legendre(l1,m1,TH[i]);
139 }
140 if (l2+m2 > 0) {
141 output[i] *= ROOT::Math::internal::legendre(l2,m2,TH[i]);
142 }
143 }
144 }
145}
146};
147
148void RooLegendre::computeBatch(double* output, size_t size, RooFit::Detail::DataMap const& dataMap) const
149{
150 compute(size, _l1, _m1, _l2, _m2, output, dataMap.at(_ctheta).data());
151}
152
153
154////////////////////////////////////////////////////////////////////////////////
155
156namespace {
157 bool fullRange(const RooRealProxy& x ,const char* range)
158 {
159 return range == nullptr || strlen(range) == 0
160 ? std::abs(x.min() + 1.) < 1.e-8 && std::abs(x.max() - 1.) < 1.e-8
161 : std::abs(x.min(range) + 1.) < 1.e-8 && std::abs(x.max(range) - 1.) < 1.e-8;
162 }
163}
164
165////////////////////////////////////////////////////////////////////////////////
166
167Int_t RooLegendre::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
168{
169 // don't support indefinite integrals...
170 if (fullRange(_ctheta,rangeName) && matchArgs(allVars, analVars, _ctheta)) return 1;
171 return 0;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// this was verified to match mathematica for
176/// l1 in [0,2], m1 in [0,l1], l2 in [l1,4], m2 in [0,l2]
177
178double RooLegendre::analyticalIntegral(Int_t code, const char* ) const
179{
180 R__ASSERT(code==1) ;
181 if ( _m1==_m2 ) return ( _l1 == _l2) ? TMath::Factorial(_l1+_m2)/TMath::Factorial(_l1-_m1)*double(2)/(2*_l1+1) : 0.;
182 if ( (_l1+_l2-_m1-_m2)%2 != 0 ) return 0; // these combinations are odd under x -> -x
183
184 // from B.R. Wong, "On the overlap integral of associated Legendre Polynomials" 1998 J. Phys. A: Math. Gen. 31 1101
185 // TODO: update to the result of
186 // H. A. Mavromatis
187 // "A single-sum expression for the overlap integral of two associated Legendre polynomials"
188 // 1999 J. Phys. A: Math. Gen. 32 2601
189 // http://iopscience.iop.org/0305-4470/32/13/011/pdf/0305-4470_32_13_011.pdf
190 // For that we need Wigner 3-j, which Lorenzo has added for Root 5.28... (note: check Condon-Shortly convention in this paper!)
191 double r=0;
192 for (int p1=0; 2*p1 <= _l1-_m1 ;++p1) {
193 double a1 = a(p1,_l1,_m1);
194 for (int p2=0; 2*p2 <= _l2-_m2 ; ++p2) {
195 double a2 = a(p2,_l2,_m2);
196 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 );
197 }
198 }
199 r /= TMath::Gamma( double(_l1+_l2+3)/2 );
200
201 if ((_m1+_m2)%2==1) r = -r;
202 return r;
203}
204
205////////////////////////////////////////////////////////////////////////////////
206
207Int_t RooLegendre::getMaxVal( const RooArgSet& /*vars*/) const {
208 if (_m1==0&&_m2==0) return 1;
209 // does anyone know the analytical expression for the max values in case m!=0??
210 if (_l1<3&&_l2<3) return 1;
211 return 0;
212}
213
214namespace {
215 inline double maxSingle(int i, int j) {
216 R__ASSERT(j<=i);
217 // x0 : 1 (ordinary Legendre)
218 if (j==0) return 1;
219 R__ASSERT(i<3);
220 // 11: 1
221 if (i<2) return 1;
222 // 21: 3 22: 3
223 static const double m2[3] = { 3,3 };
224 return m2[j-1];
225 }
226}
227double RooLegendre::maxVal( Int_t /*code*/) const {
228 return maxSingle(_l1,_m1)*maxSingle(_l2,_m2);
229}
#define a(i)
Definition RSha256.hxx:99
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
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
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
Definition DataMap.cxx:22
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.
void computeBatch(double *output, size_t size, RooFit::Detail::DataMap const &) const override
Base function for computing multiple values of a RooAbsReal.
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.
double assoc_legendre(unsigned l, unsigned m, double x)
Computes the associated Legendre polynomials.
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition RVec.hxx:1809
Double_t x[n]
Definition legend1.C:17
double legendre(unsigned l, unsigned m, double x)
Double_t Factorial(Int_t i)
Computes factorial(n).
Definition TMath.cxx:252
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition TMath.cxx:353
TMarker m
Definition textangle.C:8
TLine l
Definition textangle.C:4
static void output()