/***************************************************************************** 
  * 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