// END_HTML
#include "RooFit.h"
#include "Riostream.h"
#include <math.h>
#include <string>
#include <algorithm>
#include "RooLegendre.h"
#include "RooAbsReal.h"
#include "Math/SpecFunc.h"
#include "TMath.h"
#include "TError.h"
using namespace std;
ClassImp(RooLegendre)
;
namespace {
inline double a(int p, int l, int m) {
double r = TMath::Factorial(l+m)/TMath::Factorial(m+p)/TMath::Factorial(p)/TMath::Factorial(l-m-2*p);
r /= pow(2.,m+2*p);
return p%2==0 ? r : -r ;
}
}
RooLegendre::RooLegendre() :
_l1(1),_m1(1),_l2(0),_m2(0)
{
}
RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l, int m)
: RooAbsReal(name, title)
, _ctheta("ctheta", "ctheta", this, ctheta)
, _l1(l),_m1(m),_l2(0),_m2(0)
{
}
RooLegendre::RooLegendre(const char* name, const char* title, RooAbsReal& ctheta, int l1, int m1, int l2, int m2)
: RooAbsReal(name, title)
, _ctheta("ctheta", "ctheta", this, ctheta)
, _l1(l1),_m1(m1),_l2(l2),_m2(m2)
{
}
RooLegendre::RooLegendre(const RooLegendre& other, const char* name)
: RooAbsReal(other, name)
, _ctheta("ctheta", this, other._ctheta)
, _l1(other._l1), _m1(other._m1)
, _l2(other._l2), _m2(other._m2)
{
}
Double_t RooLegendre::evaluate() const
{
#ifdef R__HAS_MATHMORE
double r = 1;
double ctheta = std::max(-1., std::min((double)_ctheta, +1.));
if (_l1!=0||_m1!=0) r *= ROOT::Math::assoc_legendre(_l1,_m1,ctheta);
if (_l2!=0||_m2!=0) r *= ROOT::Math::assoc_legendre(_l2,_m2,ctheta);
if ((_m1+_m2)%2==1) r = -r;
return r;
#else
throw std::string("RooLegendre: ERROR: This class require installation of the MathMore library") ;
return 0 ;
#endif
}
namespace {
bool fullRange(const RooRealProxy& x ,const char* range)
{ return range==0 || strlen(range)==0
|| ( x.min(range) == x.min() && x.max(range) == x.max() ) ;
}
}
Int_t RooLegendre::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
{
if (fullRange(_ctheta,rangeName) && matchArgs(allVars, analVars, _ctheta)) return 1;
return 0;
}
Double_t RooLegendre::analyticalIntegral(Int_t code, const char* ) const
{
R__ASSERT(code==1) ;
if ( _m1==_m2 ) return ( _l1 == _l2) ? TMath::Factorial(_l1+_m2)/TMath::Factorial(_l1-_m1)*double(2)/(2*_l1+1) : 0.;
if ( (_l1+_l2-_m1-_m2)%2 != 0 ) return 0;
double r=0;
for (int p1=0; 2*p1 <= _l1-_m1 ;++p1) {
double a1 = a(p1,_l1,_m1);
for (int p2=0; 2*p2 <= _l2-_m2 ; ++p2) {
double a2 = a(p2,_l2,_m2);
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 );
}
}
r /= TMath::Gamma( double(_l1+_l2+3)/2 );
if ((_m1+_m2)%2==1) r = -r;
return r;
}
Int_t RooLegendre::getMaxVal( const RooArgSet& ) const {
if (_m1==0&&_m2==0) return 1;
if (_l1<3&&_l2<3) return 1;
return 0;
}
namespace {
inline double maxSingle(int i, int j) {
R__ASSERT(j<=i);
if (j==0) return 1;
R__ASSERT(i<3);
if (i<2) return 1;
static const double m2[3] = { 3,3 };
return m2[j-1];
}
}
Double_t RooLegendre::maxVal( Int_t ) const {
return maxSingle(_l1,_m1)*maxSingle(_l2,_m2);
}