ROOT   Reference Guide
RooGaussian.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/** \class RooGaussian
18 \ingroup Roofit
19
20Plain Gaussian p.d.f
21**/
22
23#include "RooGaussian.h"
24
25#include "RooFit.h"
26#include "BatchHelpers.h"
27#include "RooAbsReal.h"
28#include "RooRealVar.h"
29#include "RooRandom.h"
30#include "RooMath.h"
31
33
34using namespace BatchHelpers;
35using namespace std;
36
38
39////////////////////////////////////////////////////////////////////////////////
40
41RooGaussian::RooGaussian(const char *name, const char *title,
42 RooAbsReal& _x, RooAbsReal& _mean,
43 RooAbsReal& _sigma) :
44 RooAbsPdf(name,title),
45 x("x","Observable",this,_x),
46 mean("mean","Mean",this,_mean),
47 sigma("sigma","Width",this,_sigma)
48{
49}
50
51////////////////////////////////////////////////////////////////////////////////
52
53RooGaussian::RooGaussian(const RooGaussian& other, const char* name) :
54 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
55 sigma("sigma",this,other.sigma)
56{
57}
58
59////////////////////////////////////////////////////////////////////////////////
60
62{
63 const double arg = x - mean;
64 const double sig = sigma;
65 return exp(-0.5*arg*arg/(sig*sig));
66}
67
68
69namespace {
70
71///Actual computations for the batch evaluation of the Gaussian.
72///May vectorise over x, mean, sigma, depending on the types of the inputs.
73///\note The output and input spans are assumed to be non-overlapping. If they
74///overlap, results will likely be garbage.
75template<class Tx, class TMean, class TSig>
76void compute(RooSpan<double> output, Tx x, TMean mean, TSig sigma) {
77 const int n = output.size();
78
79 for (int i = 0; i < n; ++i) {
80 const double arg = x[i] - mean[i];
81 const double halfBySigmaSq = -0.5 / (sigma[i] * sigma[i]);
82
83 output[i] = _rf_fast_exp(arg*arg * halfBySigmaSq);
84 }
85}
86
87}
88
89////////////////////////////////////////////////////////////////////////////////
90/// Compute \f$\exp(-0.5 \cdot \frac{(x - \mu)^2}{\sigma^2} \f$ in batches.
91/// The local proxies {x, mean, sigma} will be searched for batch input data,
92/// and if found, the computation will be batched over their
93/// values. If batch data are not found for one of the proxies, the proxies value is assumed to
94/// be constant over the batch.
95/// \param[in] batchIndex Index of the batch to be computed.
96/// \param[in] batchSize Size of each batch. The last batch may be smaller.
97/// \return A span with the computed values.
98
99RooSpan<double> RooGaussian::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
100 auto xData = x.getValBatch(begin, batchSize);
101 auto meanData = mean.getValBatch(begin, batchSize);
102 auto sigmaData = sigma.getValBatch(begin, batchSize);
103
104 //Now explicitly write down all possible template instantiations of compute() above:
105 const bool batchX = !xData.empty();
106 const bool batchMean = !meanData.empty();
107 const bool batchSigma = !sigmaData.empty();
108
109 if (!(batchX || batchMean || batchSigma)) {
110 return {};
111 }
112
113 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
114
115 if (batchX && !batchMean && !batchSigma) {
117 }
118 else if (batchX && batchMean && !batchSigma) {
119 compute(output, xData, meanData, BracketAdapter<double>(sigma));
120 }
121 else if (batchX && !batchMean && batchSigma) {
123 }
124 else if (batchX && batchMean && batchSigma) {
125 compute(output, xData, meanData, sigmaData);
126 }
127 else if (!batchX && batchMean && !batchSigma) {
129 }
130 else if (!batchX && !batchMean && batchSigma) {
132 }
133 else if (!batchX && batchMean && batchSigma) {
135 }
136
137 return output;
138}
139
140////////////////////////////////////////////////////////////////////////////////
141
142Int_t RooGaussian::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
143{
144 if (matchArgs(allVars,analVars,x)) return 1 ;
145 if (matchArgs(allVars,analVars,mean)) return 2 ;
146 return 0 ;
147}
148
149////////////////////////////////////////////////////////////////////////////////
150
151Double_t RooGaussian::analyticalIntegral(Int_t code, const char* rangeName) const
152{
153 assert(code==1 || code==2);
154
155 //The normalisation constant 1./sqrt(2*pi*sigma^2) is left out in evaluate().
156 //Therefore, the integral is scaled up by that amount to make RooFit normalise
157 //correctly.
158 const double resultScale = sqrt(TMath::TwoPi()) * sigma;
159
160 //Here everything is scaled and shifted into a standard normal distribution:
161 const double xscale = TMath::Sqrt2() * sigma;
162 double max = 0.;
163 double min = 0.;
164 if (code == 1){
165 max = (x.max(rangeName)-mean)/xscale;
166 min = (x.min(rangeName)-mean)/xscale;
167 } else { //No == 2 test because of assert
168 max = (mean.max(rangeName)-x)/xscale;
169 min = (mean.min(rangeName)-x)/xscale;
170 }
171
172
173 //Here we go for maximum precision: We compute all integrals in the UPPER
174 //tail of the Gaussian, because erfc has the highest precision there.
175 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
176 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
177 const double ecmin = std::erfc(std::abs(min));
178 const double ecmax = std::erfc(std::abs(max));
179
180
181 return resultScale * 0.5 * (
182 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
183 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
184 );
185}
186
187////////////////////////////////////////////////////////////////////////////////
188
189Int_t RooGaussian::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
190{
191 if (matchArgs(directVars,generateVars,x)) return 1 ;
192 if (matchArgs(directVars,generateVars,mean)) return 2 ;
193 return 0 ;
194}
195
196////////////////////////////////////////////////////////////////////////////////
197
199{
200 assert(code==1 || code==2) ;
201 Double_t xgen ;
202 if(code==1){
203 while(1) {
205 if (xgen<x.max() && xgen>x.min()) {
206 x = xgen ;
207 break;
208 }
209 }
210 } else if(code==2){
211 while(1) {
213 if (xgen<mean.max() && xgen>mean.min()) {
214 mean = xgen ;
215 break;
216 }
217 }
218 } else {
219 cout << "error in RooGaussian generateEvent"<< endl;
220 }
221
222 return;
223}
double _rf_fast_exp(double x)
VDT headers for RooFit.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
double sqrt(double)
double exp(double)
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
Little adapter that gives a bracket operator to types that don't have one.
Definition: BatchHelpers.h:58
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Plain Gaussian p.d.f.
Definition: RooGaussian.h:25
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const override
Compute in batches.
Definition: RooGaussian.cxx:99
RooRealProxy mean
Definition: RooGaussian.h:45
void generateEvent(Int_t code) override
Interface for generation of an event using the algorithm corresponding to the specified code.
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const override
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
RooRealProxy sigma
Definition: RooGaussian.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.
Double_t evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooGaussian.cxx:61
RooRealProxy x
Definition: RooGaussian.h:44
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:54
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:263
double erfc(double x)
Complementary error function.
const Double_t sigma
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
constexpr Double_t Sqrt2()
Definition: TMath.h:89
constexpr Double_t TwoPi()
Definition: TMath.h:45
static void output(int code)
Definition: gifencode.c:226