ROOT   Reference Guide
Searching...
No Matches
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,
30 RooAbsReal::Ref _mean,
31 bool 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 k = _noRounding ? x : floor(x);
57 if(_protectNegative && mean<0) {
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,
68 RooFit::Detail::DataMap const &dataMap) const
69{
71 RooBatchCompute::ArgVector extraArgs{static_cast<double>(_protectNegative), static_cast<double>(_noRounding)};
72 dispatch->compute(stream, RooBatchCompute::Poisson, output, nEvents, {dataMap.at(x), dataMap.at(mean)}, extraArgs);
73}
74
75////////////////////////////////////////////////////////////////////////////////
76
77Int_t RooPoisson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
78{
79 if (matchArgs(allVars,analVars,x)) return 1 ;
80 if (matchArgs(allVars, analVars, mean)) return 2;
81 return 0 ;
82}
83
84////////////////////////////////////////////////////////////////////////////////
85
86double RooPoisson::analyticalIntegral(Int_t code, const char* rangeName) const
87{
88 R__ASSERT(code == 1 || code == 2) ;
89
90 const double mu = mean; // evaluating the proxy once for less overhead
91
92 if(_protectNegative && mu < 0.0) {
93 return std::exp(-2.0 * mu); // make it fall quickly
94 }
95
96 if (code == 1) {
97 // Implement integral over x as summation. Add special handling in case
98 // range boundaries are not on integer values of x
99 const double xmin = std::max(0., x.min(rangeName));
100 const double xmax = x.max(rangeName);
101
102 if (xmax < 0. || xmax < xmin) {
103 return 0.;
104 }
105 const double delta = 100.0 * std::sqrt(mu);
106 // If the limits are more than many standard deviations away from the mean,
107 // we might as well return the integral of the full Poisson distribution to
108 // save computing time.
109 if (xmin < std::max(mu - delta, 0.0) && xmax > mu + delta) {
110 return 1.;
111 }
112
113 // The range as integers. ixmin is included, ixmax outside.
114 const unsigned int ixmin = xmin;
115 const unsigned int ixmax = std::min(xmax + 1.,
116 (double)std::numeric_limits<unsigned int>::max());
117
118 // Sum from 0 to just before the bin outside of the range.
119 if (ixmin == 0) {
120 return ROOT::Math::poisson_cdf(ixmax - 1, mu);
121 }
122 else {
123 // If necessary, subtract from 0 to the beginning of the range
124 if (ixmin <= mu) {
125 return ROOT::Math::poisson_cdf(ixmax - 1, mu) - ROOT::Math::poisson_cdf(ixmin - 1, mu);
126 }
127 else {
128 //Avoid catastrophic cancellation in the high tails:
129 return ROOT::Math::poisson_cdf_c(ixmin - 1, mu) - ROOT::Math::poisson_cdf_c(ixmax - 1, mu);
130 }
131
132 }
133
134 } else if(code == 2) {
135
136 // the integral with respect to the mean is the integral of a gamma distribution
137
138 // negative ix does not need protection (gamma returns 0.0)
139 const double ix = _noRounding ? x + 1 : int(TMath::Floor(x)) + 1.0;
140
142 return gamma_cdf(mean.max(rangeName), ix, 1.0) - gamma_cdf(mean.min(rangeName), ix, 1.0);
143 }
144
145 return 0.0;
146}
147
148////////////////////////////////////////////////////////////////////////////////
149/// Advertise internal generator in x
150
151Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
152{
153 if (matchArgs(directVars,generateVars,x)) return 1 ;
154 return 0 ;
155}
156
157////////////////////////////////////////////////////////////////////////////////
158/// Implement internal generator using TRandom::Poisson
159
161{
162 R__ASSERT(code==1) ;
163 double xgen ;
164 while(1) {
166 if (xgen<=x.max() && xgen>=x.min()) {
167 x = xgen ;
168 break;
169 }
170 }
171 return;
172}
#define ClassImp(name)
Definition Rtypes.h:377
#define R__ASSERT(e)
Definition TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
char name[80]
Definition TGX11.cxx:110
float xmin
float xmax
A RooAbsReal::Ref can be constructed from a RooAbsReal& or a double that will be implicitly converted...
Definition RooAbsReal.h:70
bool 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:55
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, ArgVector &)=0
RooSpan< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
Definition DataMap.cxx:21
Poisson pdf.
Definition RooPoisson.h:19
RooRealProxy x
Definition RooPoisson.h:49
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of the Poisson distribution.
bool _noRounding
Definition RooPoisson.h:51
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
bool _protectNegative
Definition RooPoisson.h:52
void generateEvent(Int_t code) override
Implement internal generator using TRandom::Poisson.
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
double evaluate() const override
Implementation in terms of the TMath::Poisson() function.
RooRealProxy mean
Definition RooPoisson.h:50
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:51
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower 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
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
std::vector< double > ArgVector
Double_t Floor(Double_t x)
Rounds x downward, returning the largest integral value that is not greater than x.
Definition TMath.h:680
Double_t Poisson(Double_t x, Double_t par)
Computes the Poisson distribution function for (x,par).
Definition TMath.cxx:587
Little struct that can pack a float into the unused bits of the mantissa of a NaN double.
static void output()