Logo ROOT  
Reference Guide
RooLognormal.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2 * Project: RooFit *
3 * @(#)root/roofit:$Id$ *
4 * *
5 * RooFit Lognormal PDF *
6 * *
7 * Author: Gregory Schott and Stefan Schmitz *
8 * *
9 *****************************************************************************/
10
11/** \class RooLognormal
12 \ingroup Roofit
13
14RooFit Lognormal PDF. The two parameters are:
15 - m0: the median of the distribution
16 - k=exp(sigma): sigma is called the shape parameter in the TMath parameterization
17
18\f[ Lognormal(x,m_0,k) = \frac{e^{(-ln^2(x/m_0))/(2ln^2(k))}}{\sqrt{2\pi \cdot ln(k)\cdot x}} \f]
19
20The parameterization here is physics driven and differs from the ROOT::Math::lognormal_pdf(x,m,s,x0) with:
21 - m = log(m0)
22 - s = log(k)
23 - x0 = 0
24**/
25
26#include "RooLognormal.h"
27#include "RooFit.h"
28#include "RooAbsReal.h"
29#include "RooRealVar.h"
30#include "RooRandom.h"
31#include "RooMath.h"
32#include "RooVDTHeaders.h"
33#include "RooHelpers.h"
34#include "BatchHelpers.h"
35
36#include "TMath.h"
40
41#include <cmath>
42using namespace std;
43
45
46////////////////////////////////////////////////////////////////////////////////
47
48RooLognormal::RooLognormal(const char *name, const char *title,
49 RooAbsReal& _x, RooAbsReal& _m0,
50 RooAbsReal& _k) :
51 RooAbsPdf(name,title),
52 x("x","Observable",this,_x),
53 m0("m0","m0",this,_m0),
54 k("k","k",this,_k)
55{
56 RooHelpers::checkRangeOfParameters(this, {&_x, &_m0, &_k}, 0.);
57
58 auto par = dynamic_cast<const RooAbsRealLValue*>(&_k);
59 if (par && par->getMin()<=1 && par->getMax()>=1 ) {
60 oocoutE(this, InputArguments) << "The parameter '" << par->GetName() << "' with range [" << par->getMin("") << ", "
61 << par->getMax() << "] of the " << this->IsA()->GetName() << " '" << this->GetName()
62 << "' can reach the unsafe value 1.0 " << ". Advise to limit its range." << std::endl;
63 }
64}
65
66////////////////////////////////////////////////////////////////////////////////
67
68RooLognormal::RooLognormal(const RooLognormal& other, const char* name) :
69 RooAbsPdf(other,name), x("x",this,other.x), m0("m0",this,other.m0),
70 k("k",this,other.k)
71{
72}
73
74////////////////////////////////////////////////////////////////////////////////
75/// ln(k)<1 would correspond to sigma < 0 in the parameterization
76/// resulting by transforming a normal random variable in its
77/// standard parameterization to a lognormal random variable
78/// => treat ln(k) as -ln(k) for k<1
79
81{
83 Double_t ln_m0 = TMath::Log(m0);
84
85 Double_t ret = ROOT::Math::lognormal_pdf(x,ln_m0,ln_k);
86 return ret ;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
91namespace {
92//Author: Emmanouil Michalainas, CERN 10 September 2019
93
94template<class Tx, class Tm0, class Tk>
95void compute( size_t batchSize,
96 double * __restrict output,
97 Tx X, Tm0 M0, Tk K)
98{
99 const double rootOf2pi = std::sqrt(2 * M_PI);
100 for (size_t i=0; i<batchSize; i++) {
101 double lnxOverM0 = _rf_fast_log(X[i]/M0[i]);
102 double lnk = _rf_fast_log(K[i]);
103 if (lnk<0) lnk = -lnk;
104 double arg = lnxOverM0/lnk;
105 arg *= -0.5*arg;
106 output[i] = _rf_fast_exp(arg) / (X[i]*lnk*rootOf2pi);
107 }
108}
109};
110
111RooSpan<double> RooLognormal::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
112 using namespace BatchHelpers;
113 auto xData = x.getValBatch(begin, batchSize);
114 auto m0Data = m0.getValBatch(begin, batchSize);
115 auto kData = k.getValBatch(begin, batchSize);
116 const bool batchX = !xData.empty();
117 const bool batchM0 = !m0Data.empty();
118 const bool batchK = !kData.empty();
119
120 if (!batchX && !batchM0 && !batchK) {
121 return {};
122 }
123 batchSize = findSize({ xData, m0Data, kData });
124 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
125
126 if (batchX && !batchM0 && !batchK ) {
127 compute(batchSize, output.data(), xData, BracketAdapter<double>(m0), BracketAdapter<double>(k));
128 }
129 else if (!batchX && batchM0 && !batchK ) {
130 compute(batchSize, output.data(), BracketAdapter<double>(x), m0Data, BracketAdapter<double>(k));
131 }
132 else if (batchX && batchM0 && !batchK ) {
133 compute(batchSize, output.data(), xData, m0Data, BracketAdapter<double>(k));
134 }
135 else if (!batchX && !batchM0 && batchK ) {
136 compute(batchSize, output.data(), BracketAdapter<double>(x), BracketAdapter<double>(m0), kData);
137 }
138 else if (batchX && !batchM0 && batchK ) {
139 compute(batchSize, output.data(), xData, BracketAdapter<double>(m0), kData);
140 }
141 else if (!batchX && batchM0 && batchK ) {
142 compute(batchSize, output.data(), BracketAdapter<double>(x), m0Data, kData);
143 }
144 else if (batchX && batchM0 && batchK ) {
145 compute(batchSize, output.data(), xData, m0Data, kData);
146 }
147 return output;
148}
149
150
151////////////////////////////////////////////////////////////////////////////////
152
153Int_t RooLognormal::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
154{
155 if (matchArgs(allVars,analVars,x)) return 1 ;
156 return 0 ;
157}
158
159////////////////////////////////////////////////////////////////////////////////
160
161Double_t RooLognormal::analyticalIntegral(Int_t code, const char* rangeName) const
162{
163 R__ASSERT(code==1) ;
164
165 static const Double_t root2 = sqrt(2.) ;
166
168 Double_t ret = 0.5*( RooMath::erf( TMath::Log(x.max(rangeName)/m0)/(root2*ln_k) ) - RooMath::erf( TMath::Log(x.min(rangeName)/m0)/(root2*ln_k) ) ) ;
169
170 return ret ;
171}
172
173////////////////////////////////////////////////////////////////////////////////
174
175Int_t RooLognormal::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
176{
177 if (matchArgs(directVars,generateVars,x)) return 1 ;
178 return 0 ;
179}
180
181////////////////////////////////////////////////////////////////////////////////
182
184{
185 R__ASSERT(code==1) ;
186
187 Double_t xgen ;
188 while(1) {
190 if (xgen<=x.max() && xgen>=x.min()) {
191 x = xgen ;
192 break;
193 }
194 }
195
196 return;
197}
#define oocoutE(o, a)
Definition: RooMsgService.h:49
double _rf_fast_exp(double x)
VDT headers for RooFit.
Definition: RooVDTHeaders.h:47
double _rf_fast_log(double x)
Definition: RooVDTHeaders.h:51
#define M_PI
Definition: Rotated.cxx:105
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
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
double sqrt(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
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
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
RooFit Lognormal PDF.
Definition: RooLognormal.h:19
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
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=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooRealProxy k
Definition: RooLognormal.h:38
Double_t evaluate() const
ln(k)<1 would correspond to sigma < 0 in the parameterization resulting by transforming a normal rand...
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
RooRealProxy x
Definition: RooLognormal.h:36
RooRealProxy m0
Definition: RooLognormal.h:37
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:580
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 const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
double lognormal_pdf(double x, double m, double s, double x0=0)
Probability density function of the lognormal distribution.
Double_t x[n]
Definition: legend1.C:17
size_t findSize(std::vector< RooSpan< const double > > parameters)
This function returns the minimum size of the non-zero-sized batches.
static const uint32_t K[64]
Definition: RSha256.hxx:148
Double_t Gaus(Double_t x, Double_t mean, Double_t sigma)
Gauss.
@ InputArguments
Definition: RooGlobalFunc.h:68
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 extraMessage="")
Check if the parameters have a range, and warn if the range extends below / above the set limits.
Definition: RooHelpers.cxx:75
Double_t Exp(Double_t x)
Definition: TMath.h:717
Double_t Log(Double_t x)
Definition: TMath.h:750
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226