``` /*****************************************************************************
* Project: RooFit                                                           *
*                                                                           *
* Simple Poisson PDF
* author: Kyle Cranmer <cranmer@cern.ch>
*                                                                           *
*****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// Poisson pdf
// END_HTML
//

#include <iostream>

#include "RooPoisson.h"
#include "RooAbsReal.h"
#include "RooAbsCategory.h"

#include "RooRandom.h"
#include "RooMath.h"
#include "TMath.h"
#include "Math/ProbFuncMathCore.h"

#include "TError.h"

using namespace std;

ClassImp(RooPoisson)

//_____________________________________________________________________________
RooPoisson::RooPoisson(const char *name, const char *title,
RooAbsReal& _x,
RooAbsReal& _mean,
Bool_t noRounding) :
RooAbsPdf(name,title),
x("x","x",this,_x),
mean("mean","mean",this,_mean),
_noRounding(noRounding),
_protectNegative(false)
{
// Constructor
}

//_____________________________________________________________________________
RooPoisson::RooPoisson(const RooPoisson& other, const char* name) :
RooAbsPdf(other,name),
x("x",this,other.x),
mean("mean",this,other.mean),
_noRounding(other._noRounding),
_protectNegative(other._protectNegative)
{
// Copy constructor
}

//_____________________________________________________________________________
Double_t RooPoisson::evaluate() const
{
// Implementation in terms of the TMath Poisson function

Double_t k = _noRounding ? x : floor(x);
if(_protectNegative && mean<0)
return 1e-3;
return TMath::Poisson(k,mean) ;
}

//_____________________________________________________________________________
Double_t RooPoisson::getLogVal(const RooArgSet* s) const
{
// calculate and return the negative log-likelihood of the Poisson
return RooAbsPdf::getLogVal(s) ;
//   Double_t prob = getVal(s) ;
//   return prob ;

// Make inputs to naming conventions of RooAbsPdf::extendedTerm
Double_t expected=mean ;
Double_t observed=x ;

// Explicitly handle case Nobs=Nexp=0
if (fabs(expected)<1e-10 && fabs(observed)<1e-10) {
return 0 ;
}

// Explicitly handle case Nexp=0
if (fabs(observed)<1e-10) {
return -1*expected;
}

// Michaels code for log(poisson) in RooAbsPdf::extendedTer with an approximated log(observed!) term
Double_t extra=0;
if(observed<1000000) {
extra = -observed*log(expected)+expected+TMath::LnGamma(observed+1.);
} else {
//if many observed events, use Gauss approximation
Double_t sigma_square=expected;
Double_t diff=observed-expected;
extra=-log(sigma_square)/2 + (diff*diff)/(2*sigma_square);
}

//   if (fabs(extra)>100 || log(prob)>100) {
//     cout << "RooPoisson::getLogVal(" << GetName() << ") mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
//   }

//   if (fabs(extra+log(prob))>1) {
//     cout << "RooPoisson::getLogVal(" << GetName() << ") WARNING mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
//   }

//return log(prob);
return -extra-analyticalIntegral(1,0) ; //log(prob);

}

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

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

if(_protectNegative && mean<0)
return exp(-2*mean); // make it fall quickly

if (code == 1) {
// Implement integral over x as summation. Add special handling in case
// range boundaries are not on integer values of x
Double_t xmin = x.min(rangeName) ;
Double_t xmax = x.max(rangeName) ;

// Protect against negative lower boundaries
if (xmin<0) xmin=0 ;

Int_t ixmin = Int_t (xmin) ;
Int_t ixmax = Int_t (xmax)+1 ;

Double_t fracLoBin = 1-(xmin-ixmin) ;
Double_t fracHiBin = 1-(ixmax-xmax) ;

if (!x.hasMax()) {
if (xmin<1e-6) {
return 1 ;
} else {

// Return 1 minus integral from 0 to x.min()

if(ixmin == 0){ // first bin
return TMath::Poisson(0, mean)*(xmin-0);
}
Double_t sum(0) ;
sum += TMath::Poisson(0,mean)*fracLoBin ;
sum+= ROOT::Math::poisson_cdf(ixmin-2, mean) - ROOT::Math::poisson_cdf(0,mean) ;
sum += TMath::Poisson(ixmin-1,mean)*fracHiBin ;
return 1-sum ;
}
}

if(ixmin == ixmax-1){ // first bin
return TMath::Poisson(ixmin, mean)*(xmax-xmin);
}

Double_t sum(0) ;
sum += TMath::Poisson(ixmin,mean)*fracLoBin ;
if (RooNumber::isInfinite(xmax)){
sum+= 1.-ROOT::Math::poisson_cdf(ixmin,mean) ;
}  else {
sum+= ROOT::Math::poisson_cdf(ixmax-2, mean) - ROOT::Math::poisson_cdf(ixmin,mean) ;
sum += TMath::Poisson(ixmax-1,mean)*fracHiBin ;
}

return sum ;

} else if(code == 2) {

// the integral with respect to the mean is the integral of a gamma distribution
Double_t mean_min = mean.min(rangeName);
Double_t mean_max = mean.max(rangeName);

Double_t ix;
if(_noRounding) ix = x + 1;
else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)

return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
}

return 0;

}

//_____________________________________________________________________________
Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
{
// Advertise internal generator in x

if (matchArgs(directVars,generateVars,x)) return 1 ;
return 0 ;
}

//_____________________________________________________________________________
void RooPoisson::generateEvent(Int_t code)
{
// Implement internal generator using TRandom::Poisson

R__ASSERT(code==1) ;
Double_t xgen ;
while(1) {
xgen = RooRandom::randomGenerator()->Poisson(mean);
if (xgen<=x.max() && xgen>=x.min()) {
x = xgen ;
break;
}
}
return;
}

```
RooPoisson.cxx:1
RooPoisson.cxx:2
RooPoisson.cxx:3
RooPoisson.cxx:4
RooPoisson.cxx:5
RooPoisson.cxx:6
RooPoisson.cxx:7
RooPoisson.cxx:8
RooPoisson.cxx:9
RooPoisson.cxx:10
RooPoisson.cxx:11
RooPoisson.cxx:12
RooPoisson.cxx:13
RooPoisson.cxx:14
RooPoisson.cxx:15
RooPoisson.cxx:16
RooPoisson.cxx:17
RooPoisson.cxx:18
RooPoisson.cxx:19
RooPoisson.cxx:20
RooPoisson.cxx:21
RooPoisson.cxx:22
RooPoisson.cxx:23
RooPoisson.cxx:24
RooPoisson.cxx:25
RooPoisson.cxx:26
RooPoisson.cxx:27
RooPoisson.cxx:28
RooPoisson.cxx:29
RooPoisson.cxx:30
RooPoisson.cxx:31
RooPoisson.cxx:32
RooPoisson.cxx:33
RooPoisson.cxx:34
RooPoisson.cxx:35
RooPoisson.cxx:36
RooPoisson.cxx:37
RooPoisson.cxx:38
RooPoisson.cxx:39
RooPoisson.cxx:40
RooPoisson.cxx:41
RooPoisson.cxx:42
RooPoisson.cxx:43
RooPoisson.cxx:44
RooPoisson.cxx:45
RooPoisson.cxx:46
RooPoisson.cxx:47
RooPoisson.cxx:48
RooPoisson.cxx:49
RooPoisson.cxx:50
RooPoisson.cxx:51
RooPoisson.cxx:52
RooPoisson.cxx:53
RooPoisson.cxx:54
RooPoisson.cxx:55
RooPoisson.cxx:56
RooPoisson.cxx:57
RooPoisson.cxx:58
RooPoisson.cxx:59
RooPoisson.cxx:60
RooPoisson.cxx:61
RooPoisson.cxx:62
RooPoisson.cxx:63
RooPoisson.cxx:64
RooPoisson.cxx:65
RooPoisson.cxx:66
RooPoisson.cxx:67
RooPoisson.cxx:68
RooPoisson.cxx:69
RooPoisson.cxx:70
RooPoisson.cxx:71
RooPoisson.cxx:72
RooPoisson.cxx:73
RooPoisson.cxx:74
RooPoisson.cxx:75
RooPoisson.cxx:76
RooPoisson.cxx:77
RooPoisson.cxx:78
RooPoisson.cxx:79
RooPoisson.cxx:80
RooPoisson.cxx:81
RooPoisson.cxx:82
RooPoisson.cxx:83
RooPoisson.cxx:84
RooPoisson.cxx:85
RooPoisson.cxx:86
RooPoisson.cxx:87
RooPoisson.cxx:88
RooPoisson.cxx:89
RooPoisson.cxx:90
RooPoisson.cxx:91
RooPoisson.cxx:92
RooPoisson.cxx:93
RooPoisson.cxx:94
RooPoisson.cxx:95
RooPoisson.cxx:96
RooPoisson.cxx:97
RooPoisson.cxx:98
RooPoisson.cxx:99
RooPoisson.cxx:100
RooPoisson.cxx:101
RooPoisson.cxx:102
RooPoisson.cxx:103
RooPoisson.cxx:104
RooPoisson.cxx:105
RooPoisson.cxx:106
RooPoisson.cxx:107
RooPoisson.cxx:108
RooPoisson.cxx:109
RooPoisson.cxx:110
RooPoisson.cxx:111
RooPoisson.cxx:112
RooPoisson.cxx:113
RooPoisson.cxx:114
RooPoisson.cxx:115
RooPoisson.cxx:116
RooPoisson.cxx:117
RooPoisson.cxx:118
RooPoisson.cxx:119
RooPoisson.cxx:120
RooPoisson.cxx:121
RooPoisson.cxx:122
RooPoisson.cxx:123
RooPoisson.cxx:124
RooPoisson.cxx:125
RooPoisson.cxx:126
RooPoisson.cxx:127
RooPoisson.cxx:128
RooPoisson.cxx:129
RooPoisson.cxx:130
RooPoisson.cxx:131
RooPoisson.cxx:132
RooPoisson.cxx:133
RooPoisson.cxx:134
RooPoisson.cxx:135
RooPoisson.cxx:136
RooPoisson.cxx:137
RooPoisson.cxx:138
RooPoisson.cxx:139
RooPoisson.cxx:140
RooPoisson.cxx:141
RooPoisson.cxx:142
RooPoisson.cxx:143
RooPoisson.cxx:144
RooPoisson.cxx:145
RooPoisson.cxx:146
RooPoisson.cxx:147
RooPoisson.cxx:148
RooPoisson.cxx:149
RooPoisson.cxx:150
RooPoisson.cxx:151
RooPoisson.cxx:152
RooPoisson.cxx:153
RooPoisson.cxx:154
RooPoisson.cxx:155
RooPoisson.cxx:156
RooPoisson.cxx:157
RooPoisson.cxx:158
RooPoisson.cxx:159
RooPoisson.cxx:160
RooPoisson.cxx:161
RooPoisson.cxx:162
RooPoisson.cxx:163
RooPoisson.cxx:164
RooPoisson.cxx:165
RooPoisson.cxx:166
RooPoisson.cxx:167
RooPoisson.cxx:168
RooPoisson.cxx:169
RooPoisson.cxx:170
RooPoisson.cxx:171
RooPoisson.cxx:172
RooPoisson.cxx:173
RooPoisson.cxx:174
RooPoisson.cxx:175
RooPoisson.cxx:176
RooPoisson.cxx:177
RooPoisson.cxx:178
RooPoisson.cxx:179
RooPoisson.cxx:180
RooPoisson.cxx:181
RooPoisson.cxx:182
RooPoisson.cxx:183
RooPoisson.cxx:184
RooPoisson.cxx:185
RooPoisson.cxx:186
RooPoisson.cxx:187
RooPoisson.cxx:188
RooPoisson.cxx:189
RooPoisson.cxx:190
RooPoisson.cxx:191
RooPoisson.cxx:192
RooPoisson.cxx:193
RooPoisson.cxx:194
RooPoisson.cxx:195
RooPoisson.cxx:196
RooPoisson.cxx:197
RooPoisson.cxx:198
RooPoisson.cxx:199
RooPoisson.cxx:200
RooPoisson.cxx:201
RooPoisson.cxx:202
RooPoisson.cxx:203
RooPoisson.cxx:204
RooPoisson.cxx:205
RooPoisson.cxx:206
RooPoisson.cxx:207
RooPoisson.cxx:208
RooPoisson.cxx:209
RooPoisson.cxx:210
RooPoisson.cxx:211
RooPoisson.cxx:212
RooPoisson.cxx:213
RooPoisson.cxx:214
RooPoisson.cxx:215
RooPoisson.cxx:216
RooPoisson.cxx:217
RooPoisson.cxx:218
RooPoisson.cxx:219
RooPoisson.cxx:220
RooPoisson.cxx:221
RooPoisson.cxx:222
RooPoisson.cxx:223
RooPoisson.cxx:224
RooPoisson.cxx:225
RooPoisson.cxx:226
RooPoisson.cxx:227
RooPoisson.cxx:228
RooPoisson.cxx:229
RooPoisson.cxx:230
RooPoisson.cxx:231
RooPoisson.cxx:232
RooPoisson.cxx:233
RooPoisson.cxx:234
RooPoisson.cxx:235
RooPoisson.cxx:236
RooPoisson.cxx:237
RooPoisson.cxx:238
RooPoisson.cxx:239
RooPoisson.cxx:240
RooPoisson.cxx:241
RooPoisson.cxx:242