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 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) {
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.
67void RooPoisson::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
68{
70 dispatch->compute(stream, RooBatchCompute::Poisson, output, nEvents,
71 {dataMap.at(x), dataMap.at(mean)},
72 {static_cast<double>(_protectNegative), static_cast<double>(_noRounding)});
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 if(_protectNegative && mean<0)
91 return exp(-2*mean); // make it fall quickly
92
93 if (code == 1) {
94 // Implement integral over x as summation. Add special handling in case
95 // range boundaries are not on integer values of x
96 const double xmin = std::max(0., x.min(rangeName));
97 const double xmax = x.max(rangeName);
98
99 if (xmax < 0. || xmax < xmin) {
100 return 0.;
101 }
102 if (!x.hasMax() || RooNumber::isInfinite(xmax)) {
103 //Integrating the full Poisson distribution here
104 return 1.;
105 }
106
107 // The range as integers. ixmin is included, ixmax outside.
108 const unsigned int ixmin = xmin;
109 const unsigned int ixmax = std::min(xmax + 1.,
110 (double)std::numeric_limits<unsigned int>::max());
111
112 // Sum from 0 to just before the bin outside of the range.
113 if (ixmin == 0) {
114 return ROOT::Math::poisson_cdf(ixmax - 1, mean);
115 }
116 else {
117 // If necessary, subtract from 0 to the beginning of the range
118 if (ixmin <= mean) {
119 return ROOT::Math::poisson_cdf(ixmax - 1, mean) - ROOT::Math::poisson_cdf(ixmin - 1, mean);
120 }
121 else {
122 //Avoid catastrophic cancellation in the high tails:
123 return ROOT::Math::poisson_cdf_c(ixmin - 1, mean) - ROOT::Math::poisson_cdf_c(ixmax - 1, mean);
124 }
125
126 }
127
128 } else if(code == 2) {
129
130 // the integral with respect to the mean is the integral of a gamma distribution
131 double mean_min = mean.min(rangeName);
132 double mean_max = mean.max(rangeName);
133
134 double ix;
135 if(_noRounding) ix = x + 1;
136 else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)
137
138 return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
139 }
140
141 return 0;
142
143}
144
145////////////////////////////////////////////////////////////////////////////////
146/// Advertise internal generator in x
147
148Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
149{
150 if (matchArgs(directVars,generateVars,x)) return 1 ;
151 return 0 ;
152}
153
154////////////////////////////////////////////////////////////////////////////////
155/// Implement internal generator using TRandom::Poisson
156
158{
159 R__ASSERT(code==1) ;
160 double xgen ;
161 while(1) {
163 if (xgen<=x.max() && xgen>=x.min()) {
164 x = xgen ;
165 break;
166 }
167 }
168 return;
169}
int Int_t
Definition: RtypesCore.h:45
#define ClassImp(name)
Definition: Rtypes.h:375
#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
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
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:57
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
static Int_t isInfinite(double x)
Return true if x is infinite by RooNumBer internal specification.
Definition: RooNumber.cxx:57
Poisson pdf.
Definition: RooPoisson.h:19
RooRealProxy x
Definition: RooPoisson.h:46
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of the Poisson distribution.
Definition: RooPoisson.cxx:67
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:77
bool _noRounding
Definition: RooPoisson.h:48
bool _protectNegative
Definition: RooPoisson.h:49
void generateEvent(Int_t code) override
Implement internal generator using TRandom::Poisson.
Definition: RooPoisson.cxx:157
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Advertise internal generator in x.
Definition: RooPoisson.cxx:148
double evaluate() const override
Implementation in terms of the TMath::Poisson() function.
Definition: RooPoisson.cxx:54
RooRealProxy mean
Definition: RooPoisson.h:47
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooPoisson.cxx:86
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:51
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).
RVec< PromoteType< T > > floor(const RVec< T > &v)
Definition: RVec.hxx:1773
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
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...
Double_t Floor(Double_t x)
Rounds x downward, returning the largest integral value that is not greater than x.
Definition: TMath.h:677
Double_t Poisson(Double_t x, Double_t par)
Computes the Poisson distribution function for (x,par).
Definition: TMath.cxx:569
Little struct that can pack a float into the unused bits of the mantissa of a NaN double.
Definition: RooNaNPacker.h:28
static void output()