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 *
10 * *
11 * Redistribution and use in source and binary forms, *
12 * with or without modification, are permitted according to the terms *
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, RooFit::Detail::DataMap const& dataMap) const
116{
118 dispatch->compute(stream, RooBatchCompute::Johnson, output, nEvents,
119 {dataMap.at(_mass), dataMap.at(_mu), dataMap.at(_lambda), dataMap.at(_gamma), dataMap.at(_delta)},
121}
122
123////////////////////////////////////////////////////////////////////////////////
124
125int RooJohnson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
126{
127 if (matchArgs(allVars, analVars, _mass)) return kMass;
128 if (matchArgs(allVars, analVars, _mu)) return kMean;
129 if (matchArgs(allVars, analVars, _lambda)) return kLambda;
130 if (matchArgs(allVars, analVars, _gamma)) return kGamma;
131 if (matchArgs(allVars, analVars, _delta)) return kDelta;
132 //TODO write integral for others
133 return 0;
134}
135
136////////////////////////////////////////////////////////////////////////////////
137
138double RooJohnson::analyticalIntegral(Int_t code, const char* rangeName) const
139{
140 //The normalisation constant is left out in evaluate().
141 //Therefore, the integral is scaled up by that amount to make RooFit normalise
142 //correctly.
143 const double globalNorm = 1.;
144// const double globalNorm = sqrt(TMath::TwoPi());
145
146 //Here everything is scaled and shifted such that we only need to compute CDF(Gauss):
147 double min = -1.E300;
148 double max = 1.E300;
149 if (kMass <= code && code <= kLambda) {
150 double argMin, argMax;
151
152 if (code == kMass) {
153 argMin = (_mass.min(rangeName)-_mu)/_lambda;
154 argMax = (_mass.max(rangeName)-_mu)/_lambda;
155 } else if (code == kMean) {
156 argMin = (_mass-_mu.min(rangeName))/_lambda;
157 argMax = (_mass-_mu.max(rangeName))/_lambda;
158 } else {
159 assert(code == kLambda);
160 argMin = (_mass-_mu)/_lambda.min(rangeName);
161 argMax = (_mass-_mu)/_lambda.max(rangeName);
162 }
163
164 min = _gamma + _delta * asinh(argMin);
165 max = _gamma + _delta * asinh(argMax);
166 } else if (code == kGamma) {
167 const double arg = (_mass-_mu)/_lambda;
168 min = _gamma.min(rangeName) + _delta * asinh(arg);
169 max = _gamma.max(rangeName) + _delta * asinh(arg);
170 } else if (code == kDelta) {
171 const double arg = (_mass-_mu)/_lambda;
172 min = _gamma + _delta.min(rangeName) * asinh(arg);
173 max = _gamma + _delta.max(rangeName) * asinh(arg);
174 } else {
175 assert(false);
176 }
177
178
179
180 //Here we go for maximum precision: We compute all integrals in the UPPER
181 //tail of the Gaussian, because erfc has the highest precision there.
182 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
183 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
184 const double ecmin = std::erfc(std::abs(min/sqrt(2.)));
185 const double ecmax = std::erfc(std::abs(max/sqrt(2.)));
186
187 const double result = 0.5 * (
188 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
189 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
190 );
191
192 // Now, include the global norm that may be missing in evaluate and return
193 return globalNorm * (result != 0. ? result : 1.E-300);
194}
195
196
197
198////////////////////////////////////////////////////////////////////////////////
199/// Advertise which kind of direct event generation is supported.
200///
201/// So far, only generating mass values is supported.
202Int_t RooJohnson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
203{
204 if (matchArgs(directVars, generateVars, _mass)) return 1 ;
205// if (matchArgs(directVars, generateVars, _mu)) return 2 ;
206 return 0 ;
207}
208
209
210
211////////////////////////////////////////////////////////////////////////////////
212/// Generate events based on code obtained by getGenerator().
213///
214/// So far, only generating mass values is supported. Others will have to be generated
215/// by the slower accept/reject method.
217{
218 if (code == 1) {
219 while (true) {
220 const double gauss = RooRandom::randomGenerator()->Gaus(0., 1.);
221 const double mass = _lambda * sinh((gauss - _gamma)/_delta) + _mu;
222 if (_mass.min() <= mass && mass <= _mass.max() && _massThreshold <= mass) {
223 _mass = mass;
224 break;
225 }
226 }
227 } else {
228 throw std::logic_error("Generation in other variables not yet implemented.");
229 }
230}
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char name[80]
Definition: TGX11.cxx:110
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
bool 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:57
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
Johnson's distribution.
Definition: RooJohnson.h:24
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooJohnson.cxx:97
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooJohnson.cxx:138
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:216
RooRealProxy _mu
Definition: RooJohnson.h:51
RooRealProxy _lambda
Definition: RooJohnson.h:52
double _massThreshold
Definition: RooJohnson.h:57
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of the Johnson distribution.
Definition: RooJohnson.cxx:115
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Advertise which kind of direct event generation is supported.
Definition: RooJohnson.cxx:202
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:125
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition: RooRandom.cxx:51
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.
RVec< PromoteType< T > > asinh(const RVec< T > &v)
Definition: RVec.hxx:1769
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1739
RVec< PromoteType< T > > sinh(const RVec< T > &v)
Definition: RVec.hxx:1766
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
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
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:119
static constexpr double gauss
constexpr Double_t TwoPi()
Definition: TMath.h:44
static void output()