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#include "RooBatchCompute.h"
25#include "RooHelpers.h"
26#include "RooMath.h"
27#include "RooRandom.h"
28
29#include <vector>
30
32
33////////////////////////////////////////////////////////////////////////////////
34
35RooGaussian::RooGaussian(const char *name, const char *title,
37 RooAbsReal::Ref _sigma) :
38 RooAbsPdf(name,title),
39 x("x","Observable",this,_x),
40 mean("mean","Mean",this,_mean),
41 sigma("sigma","Width",this,_sigma)
42{
43 RooHelpers::checkRangeOfParameters(this, {&static_cast<RooAbsReal&>(_sigma)}, 0);
44}
45
46////////////////////////////////////////////////////////////////////////////////
47
48RooGaussian::RooGaussian(const RooGaussian& other, const char* name) :
49 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
50 sigma("sigma",this,other.sigma)
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55
57{
58 const double arg = x - mean;
59 const double sig = sigma;
60 return std::exp(-0.5*arg*arg/(sig*sig));
61}
62
63
64////////////////////////////////////////////////////////////////////////////////
65/// Compute multiple values of Gaussian distribution.
66void RooGaussian::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
67{
69 dispatch->compute(stream, RooBatchCompute::Gaussian, output, nEvents,
70 {dataMap.at(x), dataMap.at(mean), dataMap.at(sigma)});
71}
72
73////////////////////////////////////////////////////////////////////////////////
74
75Int_t RooGaussian::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
76{
77 if (matchArgs(allVars,analVars,x)) return 1 ;
78 if (matchArgs(allVars,analVars,mean)) return 2 ;
79 return 0 ;
80}
81
82////////////////////////////////////////////////////////////////////////////////
83
84double RooGaussian::analyticalIntegral(Int_t code, const char* rangeName) const
85{
86 assert(code==1 || code==2);
87
88 //The normalisation constant 1./sqrt(2*pi*sigma^2) is left out in evaluate().
89 //Therefore, the integral is scaled up by that amount to make RooFit normalise
90 //correctly.
91 const double resultScale = std::sqrt(TMath::TwoPi()) * sigma;
92
93 //Here everything is scaled and shifted into a standard normal distribution:
94 const double xscale = TMath::Sqrt2() * sigma;
95 double max = 0.;
96 double min = 0.;
97 if (code == 1){
98 max = (x.max(rangeName)-mean)/xscale;
99 min = (x.min(rangeName)-mean)/xscale;
100 } else { //No == 2 test because of assert
101 max = (mean.max(rangeName)-x)/xscale;
102 min = (mean.min(rangeName)-x)/xscale;
103 }
104
105
106 //Here we go for maximum precision: We compute all integrals in the UPPER
107 //tail of the Gaussian, because erfc has the highest precision there.
108 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
109 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
110 const double ecmin = std::erfc(std::abs(min));
111 const double ecmax = std::erfc(std::abs(max));
112
113
114 return resultScale * 0.5 * (
115 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
116 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
117 );
118}
119
120////////////////////////////////////////////////////////////////////////////////
121
122Int_t RooGaussian::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
123{
124 if (matchArgs(directVars,generateVars,x)) return 1 ;
125 if (matchArgs(directVars,generateVars,mean)) return 2 ;
126 return 0 ;
127}
128
129////////////////////////////////////////////////////////////////////////////////
130
132{
133 assert(code==1 || code==2) ;
134 double xgen ;
135 if(code==1){
136 while(1) {
138 if (xgen<x.max() && xgen>x.min()) {
139 x = xgen ;
140 break;
141 }
142 }
143 } else if(code==2){
144 while(1) {
146 if (xgen<mean.max() && xgen>mean.min()) {
147 mean = xgen ;
148 break;
149 }
150 }
151 } else {
152 std::cout << "error in RooGaussian generateEvent"<< std::endl;
153 }
154
155 return;
156}
#define ClassImp(name)
Definition: Rtypes.h:375
char name[80]
Definition: TGX11.cxx:110
A RooAbsReal::Ref can be constructed from a RooAbsReal& or a double that will be implicitly converted...
Definition: RooAbsReal.h:69
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
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:56
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
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
RooRealProxy mean
Definition: RooGaussian.h:57
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooGaussian.cxx:84
void generateEvent(Int_t code) override
Interface for generation of an event using the algorithm corresponding to the specified code.
void computeBatch(cudaStream_t *, double *output, size_t size, RooFit::Detail::DataMap const &) const override
Compute multiple values of Gaussian distribution.
Definition: RooGaussian.cxx:66
RooRealProxy sigma
Definition: RooGaussian.h:58
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooGaussian.cxx:75
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooGaussian.cxx:56
RooRealProxy x
Definition: RooGaussian.h:56
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 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:274
double erfc(double x)
Complementary error function.
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1778
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1783
const Double_t sigma
Double_t x[n]
Definition: legend1.C:17
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
R__EXTERN RooBatchComputeInterface * dispatchCUDA
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
void checkRangeOfParameters(const RooAbsReal *callingClass, std::initializer_list< const RooAbsReal * > pars, double min=-std::numeric_limits< double >::max(), double max=std::numeric_limits< double >::max(), bool limitsInAllowedRange=false, std::string const &extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.
Definition: RooHelpers.cxx:123
constexpr Double_t Sqrt2()
Definition: TMath.h:86
constexpr Double_t TwoPi()
Definition: TMath.h:44
static void output()