Logo 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 *
10 * and Stanford University. All rights reserved. *
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#include "RunContext.h"
29
30#include <vector>
31
33
34////////////////////////////////////////////////////////////////////////////////
35
36RooGaussian::RooGaussian(const char *name, const char *title,
37 RooAbsReal& _x, RooAbsReal& _mean,
38 RooAbsReal& _sigma) :
39 RooAbsPdf(name,title),
40 x("x","Observable",this,_x),
41 mean("mean","Mean",this,_mean),
42 sigma("sigma","Width",this,_sigma)
43{
44 RooHelpers::checkRangeOfParameters(this, {&_sigma}, 0);
45}
46
47////////////////////////////////////////////////////////////////////////////////
48
49RooGaussian::RooGaussian(const RooGaussian& other, const char* name) :
50 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
51 sigma("sigma",this,other.sigma)
52{
53}
54
55////////////////////////////////////////////////////////////////////////////////
56
58{
59 const double arg = x - mean;
60 const double sig = sigma;
61 return std::exp(-0.5*arg*arg/(sig*sig));
62}
63
64
65////////////////////////////////////////////////////////////////////////////////
66/// Compute multiple values of Gaussian distribution.
67void RooGaussian::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
68{
70 dispatch->compute(stream, RooBatchCompute::Gaussian, output, nEvents, dataMap, {&*x,&*mean,&*sigma,&*_norm});
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_t 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_t /*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_t 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}
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
char name[80]
Definition: TGX11.cxx:110
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
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 void compute(cudaStream_t *, Computer, RestrictArr, size_t, const DataMap &, const VarVector &, const ArgVector &={})=0
Plain Gaussian p.d.f.
Definition: RooGaussian.h:24
RooRealProxy mean
Definition: RooGaussian.h:44
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...
void computeBatch(cudaStream_t *, double *output, size_t size, RooBatchCompute::DataMap &) const override
Compute multiple values of Gaussian distribution.
Definition: RooGaussian.cxx:67
RooRealProxy sigma
Definition: RooGaussian.h:45
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooGaussian.cxx:75
Double_t evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooGaussian.cxx:57
RooRealProxy x
Definition: RooGaussian.h:43
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooGaussian.cxx:84
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:53
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 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.
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
std::map< DataKey, RooSpan< const double > > DataMap
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:118
constexpr Double_t Sqrt2()
Definition: TMath.h:88
constexpr Double_t TwoPi()
Definition: TMath.h:44
static void output(int code)
Definition: gifencode.c:226