/*****************************************************************************
* Project: RooFit                                                           *
*                                                                           *
* Simple Gamma distribution
* authors: Stefan A. Schmitz, Gregory Schott
*                                                                           *
*****************************************************************************/

//
// implementation of the Gamma PDF for RooFit/RooStats
// $f(x) = \frac{(x-\mu)^{\gamma-1} \cdot \exp^{(-(x-mu) / \beta}}{\Gamma(\gamma) \cdot \beta^{\gamma}}$
// defined for $x \geq 0$ if $\mu = 0$
//

// Notes from Kyle Cranmer
// Wikipedia and several sources refer to the Gamma distribution as
// G(mu|alpha,beta) = beta^alpha mu^(alpha-1) e^(-beta mu) / Gamma(alpha)
// Below is the correspondance
// Wikipedia | This Implementation
//---------------------------------
// alpha     | gamma
// beta      | 1/beta
// mu        | x
// 0         | mu
//
// Note, that for a model Pois(N|mu), a uniform prior on mu, and a measurement N
// the posterior is in the Wikipedia parametrization Gamma(mu, alpha=N+1, beta=1)
// thus for this implementation it is:
// RooGamma(_x=mu,_gamma=N+1,_beta=1,_mu=0)
// Also note, that in this case it is equivalent to
// RooPoison(N,mu) and treating the function as a PDF in mu.

#include "RooFit.h"

#include "Riostream.h"
#include "Riostream.h"
#include <math.h>

#include "RooGamma.h"
#include "RooAbsReal.h"
#include "RooRealVar.h"
#include "RooRandom.h"
#include "RooMath.h"

#include <iostream>
#include "TMath.h"

#include <Math/SpecFuncMathCore.h>
#include <Math/PdfFuncMathCore.h>
#include <Math/ProbFuncMathCore.h>

#include "TError.h"

using namespace std;

ClassImp(RooGamma)

//_____________________________________________________________________________
RooGamma::RooGamma(const char *name, const char *title,
RooAbsReal& _x, RooAbsReal& _gamma,
RooAbsReal& _beta, RooAbsReal& _mu) :
RooAbsPdf(name,title),
x("x","Observable",this,_x),
gamma("gamma","Mean",this,_gamma),
beta("beta","Width",this,_beta),
mu("mu","Para",this,_mu)
{
}

//_____________________________________________________________________________
RooGamma::RooGamma(const RooGamma& other, const char* name) :
RooAbsPdf(other,name), x("x",this,other.x), gamma("mean",this,other.gamma),
beta("beta",this,other.beta), mu("mu",this,other.mu)
{
}

//_____________________________________________________________________________
Double_t RooGamma::evaluate() const
{

Double_t arg= x ;
Double_t ret = TMath::GammaDist(arg, gamma, mu, beta) ;
return ret ;
}

//_____________________________________________________________________________
Int_t RooGamma::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
{
if (matchArgs(allVars,analVars,x)) return 1 ;
return 0 ;
}

//_____________________________________________________________________________
Double_t RooGamma::analyticalIntegral(Int_t code, const char* rangeName) const
{
R__ASSERT(code==1) ;

//integral of the Gamma distribution via ROOT::Math
Double_t integral = ROOT::Math::gamma_cdf(x.max(rangeName), gamma, beta, mu) - ROOT::Math::gamma_cdf(x.min(rangeName), gamma, beta, mu);
return integral ;

}

//_____________________________________________________________________________
Int_t RooGamma::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
{
if (matchArgs(directVars,generateVars,x)) return 1 ;
return 0 ;
}

//_____________________________________________________________________________
void RooGamma::generateEvent(Int_t code)
{
R__ASSERT(code==1) ;
//algorithm adapted from code example in:
//Marsaglia, G. and Tsang, W. W.
//A Simple Method for Generating Gamma Variables
//ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000
//
//The speed of this algorithm depends on the speed of generating normal variates.
//The algorithm is limited to $\gamma \geq 0$ !

while(1) {

double d = 0;
double c = 0;
double xgen = 0;
double v = 0;
double u = 0;
d = gamma -1./3.; c = 1./TMath::Sqrt(9.*d);

while(v <= 0.){
xgen = RooRandom::randomGenerator()->Gaus(); v = 1. + c*xgen;
}
v = v*v*v; u = RooRandom::randomGenerator()->Uniform();
if( u < 1.-.0331*(xgen*xgen)*(xgen*xgen) ) {
if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
x = ((d*v)* beta + mu) ;
break;
}
}
if( TMath::Log(u) < 0.5*xgen*xgen + d*(1.-v + TMath::Log(v)) ) {
if ( (((d*v)* beta + mu ) < x.max()) && (((d*v)* beta + mu) > x.min()) ) {
x = ((d*v)* beta + mu) ;
break;
}
}

}

return;
}


RooGamma.cxx:1
RooGamma.cxx:2
RooGamma.cxx:3
RooGamma.cxx:4
RooGamma.cxx:5
RooGamma.cxx:6
RooGamma.cxx:7
RooGamma.cxx:8
RooGamma.cxx:9
RooGamma.cxx:10
RooGamma.cxx:11
RooGamma.cxx:12
RooGamma.cxx:13
RooGamma.cxx:14
RooGamma.cxx:15
RooGamma.cxx:16
RooGamma.cxx:17
RooGamma.cxx:18
RooGamma.cxx:19
RooGamma.cxx:20
RooGamma.cxx:21
RooGamma.cxx:22
RooGamma.cxx:23
RooGamma.cxx:24
RooGamma.cxx:25
RooGamma.cxx:26
RooGamma.cxx:27
RooGamma.cxx:28
RooGamma.cxx:29
RooGamma.cxx:30
RooGamma.cxx:31
RooGamma.cxx:32
RooGamma.cxx:33
RooGamma.cxx:34
RooGamma.cxx:35
RooGamma.cxx:36
RooGamma.cxx:37
RooGamma.cxx:38
RooGamma.cxx:39
RooGamma.cxx:40
RooGamma.cxx:41
RooGamma.cxx:42
RooGamma.cxx:43
RooGamma.cxx:44
RooGamma.cxx:45
RooGamma.cxx:46
RooGamma.cxx:47
RooGamma.cxx:48
RooGamma.cxx:49
RooGamma.cxx:50
RooGamma.cxx:51
RooGamma.cxx:52
RooGamma.cxx:53
RooGamma.cxx:54
RooGamma.cxx:55
RooGamma.cxx:56
RooGamma.cxx:57
RooGamma.cxx:58
RooGamma.cxx:59
RooGamma.cxx:60
RooGamma.cxx:61
RooGamma.cxx:62
RooGamma.cxx:63
RooGamma.cxx:64
RooGamma.cxx:65
RooGamma.cxx:66
RooGamma.cxx:67
RooGamma.cxx:68
RooGamma.cxx:69
RooGamma.cxx:70
RooGamma.cxx:71
RooGamma.cxx:72
RooGamma.cxx:73
RooGamma.cxx:74
RooGamma.cxx:75
RooGamma.cxx:76
RooGamma.cxx:77
RooGamma.cxx:78
RooGamma.cxx:79
RooGamma.cxx:80
RooGamma.cxx:81
RooGamma.cxx:82
RooGamma.cxx:83
RooGamma.cxx:84
RooGamma.cxx:85
RooGamma.cxx:86
RooGamma.cxx:87
RooGamma.cxx:88
RooGamma.cxx:89
RooGamma.cxx:90
RooGamma.cxx:91
RooGamma.cxx:92
RooGamma.cxx:93
RooGamma.cxx:94
RooGamma.cxx:95
RooGamma.cxx:96
RooGamma.cxx:97
RooGamma.cxx:98
RooGamma.cxx:99
RooGamma.cxx:100
RooGamma.cxx:101
RooGamma.cxx:102
RooGamma.cxx:103
RooGamma.cxx:104
RooGamma.cxx:105
RooGamma.cxx:106
RooGamma.cxx:107
RooGamma.cxx:108
RooGamma.cxx:109
RooGamma.cxx:110
RooGamma.cxx:111
RooGamma.cxx:112
RooGamma.cxx:113
RooGamma.cxx:114
RooGamma.cxx:115
RooGamma.cxx:116
RooGamma.cxx:117
RooGamma.cxx:118
RooGamma.cxx:119
RooGamma.cxx:120
RooGamma.cxx:121
RooGamma.cxx:122
RooGamma.cxx:123
RooGamma.cxx:124
RooGamma.cxx:125
RooGamma.cxx:126
RooGamma.cxx:127
RooGamma.cxx:128
RooGamma.cxx:129
RooGamma.cxx:130
RooGamma.cxx:131
RooGamma.cxx:132
RooGamma.cxx:133
RooGamma.cxx:134
RooGamma.cxx:135
RooGamma.cxx:136
RooGamma.cxx:137
RooGamma.cxx:138
RooGamma.cxx:139
RooGamma.cxx:140
RooGamma.cxx:141
RooGamma.cxx:142
RooGamma.cxx:143
RooGamma.cxx:144
RooGamma.cxx:145
RooGamma.cxx:146
RooGamma.cxx:147
RooGamma.cxx:148
RooGamma.cxx:149
RooGamma.cxx:150
RooGamma.cxx:151
RooGamma.cxx:152
RooGamma.cxx:153
RooGamma.cxx:154
RooGamma.cxx:155
RooGamma.cxx:156
RooGamma.cxx:157
RooGamma.cxx:158
RooGamma.cxx:159
RooGamma.cxx:160
RooGamma.cxx:161
RooGamma.cxx:162
RooGamma.cxx:163
RooGamma.cxx:164
RooGamma.cxx:165
RooGamma.cxx:166
RooGamma.cxx:167
RooGamma.cxx:168
RooGamma.cxx:169
RooGamma.cxx:170
RooGamma.cxx:171