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"
33#include "RooHelpers.h"
34#include "BatchHelpers.h"
35
36#include "TMath.h"
37#include "TClass.h"
41
42#include <cmath>
43using namespace std;
44
46
47////////////////////////////////////////////////////////////////////////////////
48
49RooLognormal::RooLognormal(const char *name, const char *title,
50 RooAbsReal& _x, RooAbsReal& _m0,
51 RooAbsReal& _k) :
52 RooAbsPdf(name,title),
53 x("x","Observable",this,_x),
54 m0("m0","m0",this,_m0),
55 k("k","k",this,_k)
56{
57 RooHelpers::checkRangeOfParameters(this, {&_x, &_m0, &_k}, 0.);
58
59 auto par = dynamic_cast<const RooAbsRealLValue*>(&_k);
60 if (par && par->getMin()<=1 && par->getMax()>=1 ) {
61 oocoutE(this, InputArguments) << "The parameter '" << par->GetName() << "' with range [" << par->getMin("") << ", "
62 << par->getMax() << "] of the " << this->IsA()->GetName() << " '" << this->GetName()
63 << "' can reach the unsafe value 1.0 " << ". Advise to limit its range." << std::endl;
64 }
65}
66
67////////////////////////////////////////////////////////////////////////////////
68
69RooLognormal::RooLognormal(const RooLognormal& other, const char* name) :
70 RooAbsPdf(other,name), x("x",this,other.x), m0("m0",this,other.m0),
71 k("k",this,other.k)
72{
73}
74
75////////////////////////////////////////////////////////////////////////////////
76/// ln(k)<1 would correspond to sigma < 0 in the parameterization
77/// resulting by transforming a normal random variable in its
78/// standard parameterization to a lognormal random variable
79/// => treat ln(k) as -ln(k) for k<1
80
82{
84 Double_t ln_m0 = TMath::Log(m0);
85
86 Double_t ret = ROOT::Math::lognormal_pdf(x,ln_m0,ln_k);
87 return ret ;
88}
89
90////////////////////////////////////////////////////////////////////////////////
91
92namespace {
93//Author: Emmanouil Michalainas, CERN 10 September 2019
94
95template<class Tx, class Tm0, class Tk>
96void compute( size_t batchSize,
97 double * __restrict output,
98 Tx X, Tm0 M0, Tk K)
99{
100 const double rootOf2pi = std::sqrt(2 * M_PI);
101 for (size_t i=0; i<batchSize; i++) {
102 double lnxOverM0 = _rf_fast_log(X[i]/M0[i]);
103 double lnk = _rf_fast_log(K[i]);
104 if (lnk<0) lnk = -lnk;
105 double arg = lnxOverM0/lnk;
106 arg *= -0.5*arg;
107 output[i] = _rf_fast_exp(arg) / (X[i]*lnk*rootOf2pi);
108 }
109}
110};
111
112RooSpan<double> RooLognormal::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
113 using namespace BatchHelpers;
114 auto xData = x.getValBatch(begin, batchSize);
115 auto m0Data = m0.getValBatch(begin, batchSize);
116 auto kData = k.getValBatch(begin, batchSize);
117 const bool batchX = !xData.empty();
118 const bool batchM0 = !m0Data.empty();
119 const bool batchK = !kData.empty();
120
121 if (!batchX && !batchM0 && !batchK) {
122 return {};
123 }
124 batchSize = findSize({ xData, m0Data, kData });
125 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
126
127 if (batchX && !batchM0 && !batchK ) {
129 }
130 else if (!batchX && batchM0 && !batchK ) {
132 }
133 else if (batchX && batchM0 && !batchK ) {
134 compute(batchSize, output.data(), xData, m0Data, BracketAdapter<double>(k));
135 }
136 else if (!batchX && !batchM0 && batchK ) {
138 }
139 else if (batchX && !batchM0 && batchK ) {
140 compute(batchSize, output.data(), xData, BracketAdapter<double>(m0), kData);
141 }
142 else if (!batchX && batchM0 && batchK ) {
143 compute(batchSize, output.data(), BracketAdapter<double>(x), m0Data, kData);
144 }
145 else if (batchX && batchM0 && batchK ) {
146 compute(batchSize, output.data(), xData, m0Data, kData);
147 }
148 return output;
149}
150
151
152////////////////////////////////////////////////////////////////////////////////
153
154Int_t RooLognormal::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
155{
156 if (matchArgs(allVars,analVars,x)) return 1 ;
157 return 0 ;
158}
159
160////////////////////////////////////////////////////////////////////////////////
161
162Double_t RooLognormal::analyticalIntegral(Int_t code, const char* rangeName) const
163{
164 R__ASSERT(code==1) ;
165
166 static const Double_t root2 = sqrt(2.) ;
167
169 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) ) ) ;
170
171 return ret ;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175
176Int_t RooLognormal::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
177{
178 if (matchArgs(directVars,generateVars,x)) return 1 ;
179 return 0 ;
180}
181
182////////////////////////////////////////////////////////////////////////////////
183
185{
186 R__ASSERT(code==1) ;
187
188 Double_t xgen ;
189 while(1) {
191 if (xgen<=x.max() && xgen>=x.min()) {
192 x = xgen ;
193 break;
194 }
195 }
196
197 return;
198}
#define oocoutE(o, a)
Definition: RooMsgService.h:48
double _rf_fast_exp(double x)
double _rf_fast_log(double x)
#define M_PI
Definition: Rotated.cxx:105
#define ClassImp(name)
Definition: Rtypes.h:361
#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, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:118
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:60
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:450
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:53
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
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.
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
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:118
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