Logo ROOT  
Reference Guide
RooJohnson.cxx
Go to the documentation of this file.
1// Author: Stephan Hageboeck, CERN, May 2019
2/*****************************************************************************
3 * Project: RooFit *
4 * Authors: *
5 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
6 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
7 * *
8 * Copyright (c) 2000-2019, Regents of the University of California *
9 * and Stanford University. All rights reserved. *
10 * *
11 * Redistribution and use in source and binary forms, *
12 * with or without modification, are permitted according to the terms *
13 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
14 *****************************************************************************/
15
16/** \class RooJohnson
17 \ingroup Roofit
18
19Johnson's \f$ S_{U} \f$ distribution.
20
21This PDF results from transforming a normally distributed variable \f$ x \f$ to this form:
22\f[
23 z = \gamma + \delta \sinh^{-1}\left( \frac{x - \mu}{\lambda} \right)
24\f]
25The resulting PDF is
26\f[
27 \mathrm{PDF}[\mathrm{Johnson}\ S_U] = \frac{\delta}{\lambda\sqrt{2\pi}}
28 \frac{1}{\sqrt{1 + \left( \frac{x-\mu}{\lambda} \right)^2}}
29 \;\exp\left[-\frac{1}{2} \left(\gamma + \delta \sinh^{-1}\left(\frac{x-\mu}{\lambda}\right) \right)^2\right].
30\f]
31
32It is often used to fit a mass difference for charm decays, and therefore the variable \f$ x \f$ is called
33"mass" in the implementation. A mass threshold allows to set the PDF to zero to the left of the threshold.
34
35###References:
36Johnson, N. L. (1949). *Systems of Frequency Curves Generated by Methods of Translation*. Biometrika **36(1/2)**, 149–176. [doi:10.2307/2332539](https://doi.org/10.2307%2F2332539)
37
38\image html RooJohnson_plot.png
39
40**/
41
42#include "RooJohnson.h"
43
44#include "RooRandom.h"
45#include "RooHelpers.h"
46#include "RooBatchCompute.h"
47
48#include <cmath>
49#include "TMath.h"
50
52
53////////////////////////////////////////////////////////////////////////////////
54/// Construct a new Johnson PDF.
55///
56/// \param name Name that identifies the PDF in computations
57/// \param title Title for plotting
58/// \param mass The variable of the PDF. Often this is a mass.
59/// \param mu Location parameter of the Gaussian component.
60/// \param lambda Width parameter (>0) of the Gaussian component.
61/// \param gamma Shape parameter that distorts distribution to left/right.
62/// \param delta Shape parameter (>0) that determines strength of Gaussian-like component.
63/// \param massThreshold Set PDF to zero below this threshold.
64RooJohnson::RooJohnson(const char *name, const char *title,
65 RooAbsReal& mass, RooAbsReal& mu, RooAbsReal& lambda,
67 double massThreshold) :
68 RooAbsPdf(name,title),
69 _mass("mass", "Mass observable", this, mass),
70 _mu("mu", "Location parameter of the underlying normal distribution.", this, mu),
71 _lambda("lambda", "Width parameter of the underlying normal distribution (=2 lambda)", this, lambda),
72 _gamma("gamma", "Shift of transformation", this, gamma),
73 _delta("delta", "Scale of transformation", this, delta),
74 _massThreshold(massThreshold)
75{
76 RooHelpers::checkRangeOfParameters(this, {&lambda, &delta}, 0.);
77}
78
79
80////////////////////////////////////////////////////////////////////////////////
81/// Copy a Johnson PDF.
82RooJohnson::RooJohnson(const RooJohnson& other, const char* newName) :
83 RooAbsPdf(other, newName),
84 _mass("Mass", this, other._mass),
85 _mu("mean", this, other._mu),
86 _lambda("lambda", this, other._lambda),
87 _gamma("gamma", this, other._gamma),
88 _delta("delta", this, other._delta),
89 _massThreshold(other._massThreshold)
90{
91
92}
93
94
95////////////////////////////////////////////////////////////////////////////////
96
98{
100 return 0.;
101
102 const double arg = (_mass-_mu)/_lambda;
103 const double expo = _gamma + _delta * asinh(arg);
104
105 const double result = _delta
106 / sqrt(TMath::TwoPi())
107 / (_lambda * sqrt(1. + arg*arg))
108 * exp(-0.5 * expo * expo);
109
110 return result;
111}
112
113////////////////////////////////////////////////////////////////////////////////
114/// Compute multiple values of the Johnson distribution.
115void RooJohnson::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
116{
118 dispatch->compute(stream, RooBatchCompute::Johnson, output, nEvents, dataMap, {&*_mass,&*_mu,&*_lambda,&*_gamma,&*_delta,&*_norm},{_massThreshold});
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
123int RooJohnson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
124{
125 if (matchArgs(allVars, analVars, _mass)) return kMass;
126 if (matchArgs(allVars, analVars, _mu)) return kMean;
127 if (matchArgs(allVars, analVars, _lambda)) return kLambda;
128 if (matchArgs(allVars, analVars, _gamma)) return kGamma;
129 if (matchArgs(allVars, analVars, _delta)) return kDelta;
130 //TODO write integral for others
131 return 0;
132}
133
134////////////////////////////////////////////////////////////////////////////////
135
136double RooJohnson::analyticalIntegral(Int_t code, const char* rangeName) const
137{
138 //The normalisation constant is left out in evaluate().
139 //Therefore, the integral is scaled up by that amount to make RooFit normalise
140 //correctly.
141 const double globalNorm = 1.;
142// const double globalNorm = sqrt(TMath::TwoPi());
143
144 //Here everything is scaled and shifted such that we only need to compute CDF(Gauss):
145 double min = -1.E300;
146 double max = 1.E300;
147 if (kMass <= code && code <= kLambda) {
148 double argMin, argMax;
149
150 if (code == kMass) {
151 argMin = (_mass.min(rangeName)-_mu)/_lambda;
152 argMax = (_mass.max(rangeName)-_mu)/_lambda;
153 } else if (code == kMean) {
154 argMin = (_mass-_mu.min(rangeName))/_lambda;
155 argMax = (_mass-_mu.max(rangeName))/_lambda;
156 } else {
157 assert(code == kLambda);
158 argMin = (_mass-_mu)/_lambda.min(rangeName);
159 argMax = (_mass-_mu)/_lambda.max(rangeName);
160 }
161
162 min = _gamma + _delta * asinh(argMin);
163 max = _gamma + _delta * asinh(argMax);
164 } else if (code == kGamma) {
165 const double arg = (_mass-_mu)/_lambda;
166 min = _gamma.min(rangeName) + _delta * asinh(arg);
167 max = _gamma.max(rangeName) + _delta * asinh(arg);
168 } else if (code == kDelta) {
169 const double arg = (_mass-_mu)/_lambda;
170 min = _gamma + _delta.min(rangeName) * asinh(arg);
171 max = _gamma + _delta.max(rangeName) * asinh(arg);
172 } else {
173 assert(false);
174 }
175
176
177
178 //Here we go for maximum precision: We compute all integrals in the UPPER
179 //tail of the Gaussian, because erfc has the highest precision there.
180 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
181 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
182 const double ecmin = std::erfc(std::abs(min/sqrt(2.)));
183 const double ecmax = std::erfc(std::abs(max/sqrt(2.)));
184
185 const double result = 0.5 * (
186 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
187 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
188 );
189
190 // Now, include the global norm that may be missing in evaluate and return
191 return globalNorm * (result != 0. ? result : 1.E-300);
192}
193
194
195
196////////////////////////////////////////////////////////////////////////////////
197/// Advertise which kind of direct event generation is supported.
198///
199/// So far, only generating mass values is supported.
200Int_t RooJohnson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
201{
202 if (matchArgs(directVars, generateVars, _mass)) return 1 ;
203// if (matchArgs(directVars, generateVars, _mu)) return 2 ;
204 return 0 ;
205}
206
207
208
209////////////////////////////////////////////////////////////////////////////////
210/// Generate events based on code obtained by getGenerator().
211///
212/// So far, only generating mass values is supported. Others will have to be generated
213/// by the slower accept/reject method.
215{
216 if (code == 1) {
217 while (true) {
218 const double gauss = RooRandom::randomGenerator()->Gaus(0., 1.);
219 const double mass = _lambda * sinh((gauss - _gamma)/_delta) + _mu;
220 if (_mass.min() <= mass && mass <= _mass.max() && _massThreshold <= mass) {
221 _mass = mass;
222 break;
223 }
224 }
225 } else {
226 throw std::logic_error("Generation in other variables not yet implemented.");
227 }
228}
#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
Johnson's distribution.
Definition: RooJohnson.h:24
Double_t evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooJohnson.cxx:97
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooJohnson.cxx:136
RooRealProxy _mass
Definition: RooJohnson.h:50
RooRealProxy _delta
Definition: RooJohnson.h:55
RooRealProxy _gamma
Definition: RooJohnson.h:54
void generateEvent(Int_t code) override
Generate events based on code obtained by getGenerator().
Definition: RooJohnson.cxx:214
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const override
Advertise which kind of direct event generation is supported.
Definition: RooJohnson.cxx:200
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const override
Compute multiple values of the Johnson distribution.
Definition: RooJohnson.cxx:115
RooRealProxy _mu
Definition: RooJohnson.h:51
RooRealProxy _lambda
Definition: RooJohnson.h:52
double _massThreshold
Definition: RooJohnson.h:57
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooJohnson.cxx:123
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.
double gamma(double x)
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
static constexpr double gauss
constexpr Double_t TwoPi()
Definition: TMath.h:44
static void output(int code)
Definition: gifencode.c:226