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 "BatchHelpers.h"
47#include "RooVDTHeaders.h"
48
49#include <cmath>
50#include "TMath.h"
51
52using namespace BatchHelpers;
53
55
56////////////////////////////////////////////////////////////////////////////////
57/// Construct a new Johnson PDF.
58///
59/// \param name Name that identifies the PDF in computations
60/// \param title Title for plotting
61/// \param mass The variable of the PDF. Often this is a mass.
62/// \param mu Location parameter of the Gaussian component.
63/// \param lambda Width parameter (>0) of the Gaussian component.
64/// \param gamma Shape parameter that distorts distribution to left/right.
65/// \param delta Shape parameter (>0) that determines strength of Gaussian-like component.
66/// \param massThreshold Set PDF to zero below this threshold.
67RooJohnson::RooJohnson(const char *name, const char *title,
68 RooAbsReal& mass, RooAbsReal& mu, RooAbsReal& lambda,
70 double massThreshold) :
71 RooAbsPdf(name,title),
72 _mass("mass", "Mass observable", this, mass),
73 _mu("mu", "Location parameter of the underlying normal distribution.", this, mu),
74 _lambda("lambda", "Width parameter of the underlying normal distribution (=2 lambda)", this, lambda),
75 _gamma("gamma", "Shift of transformation", this, gamma),
76 _delta("delta", "Scale of transformation", this, delta),
77 _massThreshold(massThreshold)
78{
79 RooHelpers::checkRangeOfParameters(this, {&lambda, &delta}, 0.);
80}
81
82
83////////////////////////////////////////////////////////////////////////////////
84/// Copy a Johnson PDF.
85RooJohnson::RooJohnson(const RooJohnson& other, const char* newName) :
86 RooAbsPdf(other, newName),
87 _mass("Mass", this, other._mass),
88 _mu("mean", this, other._mu),
89 _lambda("lambda", this, other._lambda),
90 _gamma("gamma", this, other._gamma),
91 _delta("delta", this, other._delta),
92 _massThreshold(other._massThreshold)
93{
94
95}
96
97
98////////////////////////////////////////////////////////////////////////////////
99
101{
102 if (_mass < _massThreshold)
103 return 0.;
104
105 const double arg = (_mass-_mu)/_lambda;
106 const double expo = _gamma + _delta * asinh(arg);
107
108 const double result = _delta
109 / sqrt(TMath::TwoPi())
110 / (_lambda * sqrt(1. + arg*arg))
111 * exp(-0.5 * expo * expo);
112
113 return result;
114}
115
116////////////////////////////////////////////////////////////////////////////////
117
118namespace {
119
120///Actual computations for the batch evaluation of the Johnson.
121///May vectorise over observables depending on types of inputs.
122///\note The output and input spans are assumed to be non-overlapping. If they
123///overlap, results will likely be garbage.
124template<class TMass, class TMu, class TLambda, class TGamma, class TDelta>
125void compute(RooSpan<double> output, TMass mass, TMu mu, TLambda lambda, TGamma gamma,
126 TDelta delta, double massThreshold) {
127 const int n = output.size();
128
129 const double sqrt_twoPi = sqrt(TMath::TwoPi());
130
131 for (int i = 0; i < n; ++i) { //CHECK_VECTORISE
132 const double arg = (mass[i] - mu[i]) / lambda[i];
133#ifdef R__HAS_VDT
134 const double asinh_arg = _rf_fast_log(arg + std::sqrt(arg*arg+1));
135#else
136 const double asinh_arg = asinh(arg);
137#endif
138 const double expo = gamma[i] + delta[i] * asinh_arg;
139 const double result = delta[i] / sqrt_twoPi
140 / (lambda[i] * std::sqrt(1. + arg*arg))
141 * _rf_fast_exp(-0.5 * expo * expo);
142
143 const double passThrough = mass[i] >= massThreshold;
144 output[i] = result * passThrough;
145 }
146}
147
148}
149
150////////////////////////////////////////////////////////////////////////////////
151/// Compute \f$ \exp(-0.5 \cdot \frac{(x - \mu)^2}{\sigma^2} \f$ in batches.
152/// The local proxies {x, mean, sigma} will be searched for batch input data,
153/// and if found, the computation will be batched over their
154/// values. If batch data are not found for one of the proxies, the proxies value is assumed to
155/// be constant over the batch.
156/// \param[in] batchIndex Index of the batch to be computed.
157/// \param[in] maxSize Maximal size of the batches. May return smaller batches depending on inputs.
158/// \return A span with the computed values.
159
160RooSpan<double> RooJohnson::evaluateBatch(std::size_t begin, std::size_t maxSize) const {
161 auto massData = _mass.getValBatch(begin, maxSize);
162 auto muData = _mu.getValBatch(begin, maxSize);
163 auto lambdaData = _lambda.getValBatch(begin, maxSize);
164 auto gammaData = _gamma.getValBatch(begin, maxSize);
165 auto deltaData = _delta.getValBatch(begin, maxSize);
166
167 maxSize = std::min({massData, muData, lambdaData, gammaData, deltaData},
169 return l.size() != 0 && l.size() < r.size();
170 }).size();
171
172 if (maxSize == 0) {
173 return {};
174 }
175
176 auto output = _batchData.makeWritableBatchUnInit(begin, maxSize);
177
178 if (!massData.empty()
179 && (muData.empty() && lambdaData.empty() && gammaData.empty() && deltaData.empty())) {
180 compute(output, massData, BracketAdapter<double>(_mu),
183 }
184 else {
185 compute(output,
186 BracketAdapterWithMask(_mass, massData),
188 BracketAdapterWithMask(_lambda, lambdaData),
189 BracketAdapterWithMask(_gamma, gammaData),
191 }
192
193 return output;
194}
195
196////////////////////////////////////////////////////////////////////////////////
197
198int RooJohnson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
199{
200 if (matchArgs(allVars, analVars, _mass)) return kMass;
201 if (matchArgs(allVars, analVars, _mu)) return kMean;
202 if (matchArgs(allVars, analVars, _lambda)) return kLambda;
203 if (matchArgs(allVars, analVars, _gamma)) return kGamma;
204 if (matchArgs(allVars, analVars, _delta)) return kDelta;
205 //TODO write integral for others
206 return 0;
207}
208
209////////////////////////////////////////////////////////////////////////////////
210
211double RooJohnson::analyticalIntegral(Int_t code, const char* rangeName) const
212{
213 //The normalisation constant is left out in evaluate().
214 //Therefore, the integral is scaled up by that amount to make RooFit normalise
215 //correctly.
216 const double globalNorm = 1.;
217// const double globalNorm = sqrt(TMath::TwoPi());
218
219 //Here everything is scaled and shifted such that we only need to compute CDF(Gauss):
220 double min = -1.E300;
221 double max = 1.E300;
222 if (kMass <= code && code <= kLambda) {
223 double argMin, argMax;
224
225 if (code == kMass) {
226 argMin = (_mass.min(rangeName)-_mu)/_lambda;
227 argMax = (_mass.max(rangeName)-_mu)/_lambda;
228 } else if (code == kMean) {
229 argMin = (_mass-_mu.min(rangeName))/_lambda;
230 argMax = (_mass-_mu.max(rangeName))/_lambda;
231 } else {
232 assert(code == kLambda);
233 argMin = (_mass-_mu)/_lambda.min(rangeName);
234 argMax = (_mass-_mu)/_lambda.max(rangeName);
235 }
236
237 min = _gamma + _delta * asinh(argMin);
238 max = _gamma + _delta * asinh(argMax);
239 } else if (code == kGamma) {
240 const double arg = (_mass-_mu)/_lambda;
241 min = _gamma.min(rangeName) + _delta * asinh(arg);
242 max = _gamma.max(rangeName) + _delta * asinh(arg);
243 } else if (code == kDelta) {
244 const double arg = (_mass-_mu)/_lambda;
245 min = _gamma + _delta.min(rangeName) * asinh(arg);
246 max = _gamma + _delta.max(rangeName) * asinh(arg);
247 } else {
248 assert(false);
249 }
250
251
252
253 //Here we go for maximum precision: We compute all integrals in the UPPER
254 //tail of the Gaussian, because erfc has the highest precision there.
255 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
256 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
257 const double ecmin = std::erfc(std::abs(min/sqrt(2.)));
258 const double ecmax = std::erfc(std::abs(max/sqrt(2.)));
259
260 const double result = 0.5 * (
261 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
262 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
263 );
264
265 // Now, include the global norm that may be missing in evaluate and return
266 return globalNorm * (result != 0. ? result : 1.E-300);
267}
268
269
270
271////////////////////////////////////////////////////////////////////////////////
272/// Advertise which kind of direct event generation is supported.
273///
274/// So far, only generating mass values is supported.
275Int_t RooJohnson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
276{
277 if (matchArgs(directVars, generateVars, _mass)) return 1 ;
278// if (matchArgs(directVars, generateVars, _mu)) return 2 ;
279 return 0 ;
280}
281
282
283
284////////////////////////////////////////////////////////////////////////////////
285/// Generate events based on code obtained by getGenerator().
286///
287/// So far, only generating mass values is supported. Others will have to be generated
288/// by the slower accept/reject method.
290{
291 if (code == 1) {
292 while (true) {
293 const double gauss = RooRandom::randomGenerator()->Gaus(0., 1.);
294 const double mass = _lambda * sinh((gauss - _gamma)/_delta) + _mu;
295 if (_mass.min() <= mass && mass <= _mass.max() && _massThreshold <= mass) {
296 _mass = mass;
297 break;
298 }
299 }
300 } else {
301 throw std::logic_error("Generation in other variables not yet implemented.");
302 }
303}
ROOT::R::TRInterface & r
Definition: Object.C:4
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 ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
double sinh(double)
double sqrt(double)
double exp(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
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
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:100
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooJohnson.cxx:211
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:289
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:275
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t end) const override
Compute in batches.
Definition: RooJohnson.cxx:160
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:198
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 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:263
double erfc(double x)
Complementary error function.
const Int_t n
Definition: legend1.C:16
double gamma(double x)
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
static constexpr double gauss
constexpr Double_t TwoPi()
Definition: TMath.h:45
auto * l
Definition: textangle.C:4
static void output(int code)
Definition: gifencode.c:226