Logo 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;
59 np.setPayload(-mean);
60 return np._payload;
61 }
62 return TMath::Poisson(k,mean) ;
63}
64
65////////////////////////////////////////////////////////////////////////////////
66/// Compute multiple values of the Poisson distribution.
68 return RooBatchCompute::dispatch->computePoisson(this, evalData, x->getValues(evalData, normSet), mean->getValues(evalData, normSet), _protectNegative, _noRounding);
69}
70
71////////////////////////////////////////////////////////////////////////////////
72
73Int_t RooPoisson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
74{
75 if (matchArgs(allVars,analVars,x)) return 1 ;
76 if (matchArgs(allVars, analVars, mean)) return 2;
77 return 0 ;
78}
79
80////////////////////////////////////////////////////////////////////////////////
81
82Double_t RooPoisson::analyticalIntegral(Int_t code, const char* rangeName) const
83{
84 R__ASSERT(code == 1 || code == 2) ;
85
86 if(_protectNegative && mean<0)
87 return exp(-2*mean); // make it fall quickly
88
89 if (code == 1) {
90 // Implement integral over x as summation. Add special handling in case
91 // range boundaries are not on integer values of x
92 const double xmin = std::max(0., x.min(rangeName));
93 const double xmax = x.max(rangeName);
94
95 if (xmax < 0. || xmax < xmin) {
96 return 0.;
97 }
99 //Integrating the full Poisson distribution here
100 return 1.;
101 }
102
103 // The range as integers. ixmin is included, ixmax outside.
104 const unsigned int ixmin = xmin;
105 const unsigned int ixmax = std::min(xmax + 1.,
106 (double)std::numeric_limits<unsigned int>::max());
107
108 // Sum from 0 to just before the bin outside of the range.
109 if (ixmin == 0) {
110 return ROOT::Math::poisson_cdf(ixmax - 1, mean);
111 }
112 else {
113 // If necessary, subtract from 0 to the beginning of the range
114 if (ixmin <= mean) {
115 return ROOT::Math::poisson_cdf(ixmax - 1, mean) - ROOT::Math::poisson_cdf(ixmin - 1, mean);
116 }
117 else {
118 //Avoid catastrophic cancellation in the high tails:
119 return ROOT::Math::poisson_cdf_c(ixmin - 1, mean) - ROOT::Math::poisson_cdf_c(ixmax - 1, mean);
120 }
121
122 }
123
124 } else if(code == 2) {
125
126 // the integral with respect to the mean is the integral of a gamma distribution
127 Double_t mean_min = mean.min(rangeName);
128 Double_t mean_max = mean.max(rangeName);
129
130 Double_t ix;
131 if(_noRounding) ix = x + 1;
132 else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)
133
134 return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
135 }
136
137 return 0;
138
139}
140
141////////////////////////////////////////////////////////////////////////////////
142/// Advertise internal generator in x
143
144Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
145{
146 if (matchArgs(directVars,generateVars,x)) return 1 ;
147 return 0 ;
148}
149
150////////////////////////////////////////////////////////////////////////////////
151/// Implement internal generator using TRandom::Poisson
152
154{
155 R__ASSERT(code==1) ;
156 Double_t xgen ;
157 while(1) {
159 if (xgen<=x.max() && xgen>=x.min()) {
160 x = xgen ;
161 break;
162 }
163 }
164 return;
165}
int Int_t
Definition: RtypesCore.h:45
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
#define ClassImp(name)
Definition: Rtypes.h:364
#define R__ASSERT(e)
Definition: TError.h:118
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
double floor(double)
double exp(double)
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:61
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
by this change, please consult the release notes for ROOT 6.24 for guidance on how to make this trans...
Definition: RooAbsReal.cxx:311
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 RooSpan< double > computePoisson(const RooAbsReal *, RunContext &, RooSpan< const double > x, RooSpan< const double > mean, bool protectNegative, bool noRounding)=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:40
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:73
Bool_t _protectNegative
Definition: RooPoisson.h:43
RooSpan< double > evaluateSpan(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet) const override
Compute multiple values of the Poisson distribution.
Definition: RooPoisson.cxx:67
void generateEvent(Int_t code) override
Implement internal generator using TRandom::Poisson.
Definition: RooPoisson.cxx:153
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const override
Advertise internal generator in x.
Definition: RooPoisson.cxx:144
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooPoisson.cxx:82
Double_t evaluate() const override
Implementation in terms of the TMath::Poisson() function.
Definition: RooPoisson.cxx:54
Bool_t _noRounding
Definition: RooPoisson.h:42
RooRealProxy mean
Definition: RooPoisson.h:41
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:53
A simple container to hold a batch of data values.
Definition: RooSpan.h:34
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 * dispatch
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:703
Double_t Poisson(Double_t x, Double_t par)
Compute the Poisson distribution function for (x,par).
Definition: TMath.cxx:564
This struct enables passing computation data around between elements of a computation graph.
Definition: RunContext.h:31
Little struct that can pack a float into the unused bits of the mantissa of a NaN double.
Definition: RooNaNPacker.h:28
double _payload
Definition: RooNaNPacker.h:29
void setPayload(float payload)
Pack float into mantissa of NaN.
Definition: RooNaNPacker.h:46