 ROOT   Reference Guide RooPoisson.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2 * Project: RooFit *
3 * *
4 * Simple Poisson PDF
5 * author: Kyle Cranmer <cranmer@cern.ch>
6 * *
7 *****************************************************************************/
8
9/** \class RooPoisson
10 \ingroup Roofit
11
12Poisson pdf
13**/
14
15#include "RooPoisson.h"
16#include "RooRandom.h"
17#include "RooMath.h"
18#include "RooNaNPacker.h"
19#include "RooBatchCompute.h"
20
22
24
25////////////////////////////////////////////////////////////////////////////////
26/// Constructor
27
28RooPoisson::RooPoisson(const char *name, const char *title,
29 RooAbsReal& _x,
30 RooAbsReal& _mean,
31 Bool_t noRounding) :
32 RooAbsPdf(name,title),
33 x("x","x",this,_x),
34 mean("mean","mean",this,_mean),
35 _noRounding(noRounding)
36{
37}
38
39////////////////////////////////////////////////////////////////////////////////
40/// Copy constructor
41
42 RooPoisson::RooPoisson(const RooPoisson& other, const char* name) :
43 RooAbsPdf(other,name),
44 x("x",this,other.x),
45 mean("mean",this,other.mean),
46 _noRounding(other._noRounding),
47 _protectNegative(other._protectNegative)
48{
49}
50
51////////////////////////////////////////////////////////////////////////////////
52/// Implementation in terms of the TMath::Poisson() function.
53
55{
56 Double_t k = _noRounding ? x : floor(x);
57 if(_protectNegative && mean<0) {
58 RooNaNPacker np;
61 }
62 return TMath::Poisson(k,mean) ;
63}
64
65////////////////////////////////////////////////////////////////////////////////
66/// Compute multiple values of the Poisson distribution.
67void RooPoisson::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
68{
70 dispatch->compute(stream, RooBatchCompute::Poisson, output, nEvents, dataMap, {&*x,&*mean,&*_norm},
71 {static_cast<double>(_protectNegative), static_cast<double>(_noRounding)});
72}
73
74////////////////////////////////////////////////////////////////////////////////
75
76Int_t RooPoisson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
77{
78 if (matchArgs(allVars,analVars,x)) return 1 ;
79 if (matchArgs(allVars, analVars, mean)) return 2;
80 return 0 ;
81}
82
83////////////////////////////////////////////////////////////////////////////////
84
85Double_t RooPoisson::analyticalIntegral(Int_t code, const char* rangeName) const
86{
87 R__ASSERT(code == 1 || code == 2) ;
88
89 if(_protectNegative && mean<0)
90 return exp(-2*mean); // make it fall quickly
91
92 if (code == 1) {
93 // Implement integral over x as summation. Add special handling in case
94 // range boundaries are not on integer values of x
95 const double xmin = std::max(0., x.min(rangeName));
96 const double xmax = x.max(rangeName);
97
98 if (xmax < 0. || xmax < xmin) {
99 return 0.;
100 }
101 if (!x.hasMax() || RooNumber::isInfinite(xmax)) {
102 //Integrating the full Poisson distribution here
103 return 1.;
104 }
105
106 // The range as integers. ixmin is included, ixmax outside.
107 const unsigned int ixmin = xmin;
108 const unsigned int ixmax = std::min(xmax + 1.,
109 (double)std::numeric_limits<unsigned int>::max());
110
111 // Sum from 0 to just before the bin outside of the range.
112 if (ixmin == 0) {
113 return ROOT::Math::poisson_cdf(ixmax - 1, mean);
114 }
115 else {
116 // If necessary, subtract from 0 to the beginning of the range
117 if (ixmin <= mean) {
118 return ROOT::Math::poisson_cdf(ixmax - 1, mean) - ROOT::Math::poisson_cdf(ixmin - 1, mean);
119 }
120 else {
121 //Avoid catastrophic cancellation in the high tails:
122 return ROOT::Math::poisson_cdf_c(ixmin - 1, mean) - ROOT::Math::poisson_cdf_c(ixmax - 1, mean);
123 }
124
125 }
126
127 } else if(code == 2) {
128
129 // the integral with respect to the mean is the integral of a gamma distribution
130 Double_t mean_min = mean.min(rangeName);
131 Double_t mean_max = mean.max(rangeName);
132
133 Double_t ix;
134 if(_noRounding) ix = x + 1;
135 else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)
136
137 return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
138 }
139
140 return 0;
141
142}
143
144////////////////////////////////////////////////////////////////////////////////
145/// Advertise internal generator in x
146
147Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
148{
149 if (matchArgs(directVars,generateVars,x)) return 1 ;
150 return 0 ;
151}
152
153////////////////////////////////////////////////////////////////////////////////
154/// Implement internal generator using TRandom::Poisson
155
157{
158 R__ASSERT(code==1) ;
159 Double_t xgen ;
160 while(1) {
162 if (xgen<=x.max() && xgen>=x.min()) {
163 x = xgen ;
164 break;
165 }
166 }
167 return;
168}
int Int_t
Definition: RtypesCore.h:45
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
#define R__ASSERT(e)
Definition: TError.h:118
char name
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
Bool_t 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:35
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const DataMap &, const VarVector &, const ArgVector &={})=0
static Int_t isInfinite(Double_t x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:58
Poisson pdf.
Definition: RooPoisson.h:19
RooRealProxy x
Definition: RooPoisson.h:46
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooPoisson.cxx:76
Bool_t _protectNegative
Definition: RooPoisson.h:49
void generateEvent(Int_t code) override
Implement internal generator using TRandom::Poisson.
Definition: RooPoisson.cxx:156
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const override
Definition: RooPoisson.cxx:147
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooPoisson.cxx:85
Double_t evaluate() const override
Implementation in terms of the TMath::Poisson() function.
Definition: RooPoisson.cxx:54
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const override
Compute multiple values of the Poisson distribution.
Definition: RooPoisson.cxx:67
Bool_t _noRounding
Definition: RooPoisson.h:48
RooRealProxy mean
Definition: RooPoisson.h:47
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:53
bool hasMax(const char *rname=0) const
Check if the range has a upper bound. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
double max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:402
double poisson_cdf(unsigned int n, double mu)
Cumulative distribution function of the Poisson distribution Lower tail of the integral of the poisso...
double poisson_cdf_c(unsigned int n, double mu)
Complement of the cumulative distribution function of the Poisson distribution.
double gamma_cdf(double x, double alpha, double theta, double x0=0)
Cumulative distribution function of the gamma distribution (lower tail).
Double_t x[n]
Definition: legend1.C:17
R__EXTERN RooBatchComputeInterface * dispatchCUDA
std::map< DataKey, RooSpan< const double > > DataMap
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
Double_t Floor(Double_t x)
Definition: TMath.h:653
Double_t Poisson(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par).
Definition: TMath.cxx:564
Little struct that can pack a float into the unused bits of the mantissa of a NaN double.
Definition: RooNaNPacker.h:28