ROOT   Reference Guide
Searching...
No Matches
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,
66 RooAbsReal& gamma, RooAbsReal& delta,
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(double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
116{
119 {dataMap.at(_mass), dataMap.at(_mu), dataMap.at(_lambda), dataMap.at(_gamma), dataMap.at(_delta)},
120 extraArgs);
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;
151 double argMax;
152
153 if (code == kMass) {
154 argMin = (_mass.min(rangeName)-_mu)/_lambda;
155 argMax = (_mass.max(rangeName)-_mu)/_lambda;
156 } else if (code == kMean) {
157 argMin = (_mass-_mu.min(rangeName))/_lambda;
158 argMax = (_mass-_mu.max(rangeName))/_lambda;
159 } else {
160 assert(code == kLambda);
161 argMin = (_mass-_mu)/_lambda.min(rangeName);
162 argMax = (_mass-_mu)/_lambda.max(rangeName);
163 }
164
165 min = _gamma + _delta * asinh(argMin);
166 max = _gamma + _delta * asinh(argMax);
167 } else if (code == kGamma) {
168 const double arg = (_mass-_mu)/_lambda;
169 min = _gamma.min(rangeName) + _delta * asinh(arg);
170 max = _gamma.max(rangeName) + _delta * asinh(arg);
171 } else if (code == kDelta) {
172 const double arg = (_mass-_mu)/_lambda;
173 min = _gamma + _delta.min(rangeName) * asinh(arg);
174 max = _gamma + _delta.max(rangeName) * asinh(arg);
175 } else {
176 assert(false);
177 }
178
179
180
181 //Here we go for maximum precision: We compute all integrals in the UPPER
182 //tail of the Gaussian, because erfc has the highest precision there.
183 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
184 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
185 const double ecmin = std::erfc(std::abs(min/sqrt(2.)));
186 const double ecmax = std::erfc(std::abs(max/sqrt(2.)));
187
188 const double result = 0.5 * (
189 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
190 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
191 );
192
193 // Now, include the global norm that may be missing in evaluate and return
194 return globalNorm * (result != 0. ? result : 1.E-300);
195}
196
197
198
199////////////////////////////////////////////////////////////////////////////////
200/// Advertise which kind of direct event generation is supported.
201///
202/// So far, only generating mass values is supported.
203Int_t RooJohnson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
204{
205 if (matchArgs(directVars, generateVars, _mass)) return 1 ;
206// if (matchArgs(directVars, generateVars, _mu)) return 2 ;
207 return 0 ;
208}
209
210
211
212////////////////////////////////////////////////////////////////////////////////
213/// Generate events based on code obtained by getGenerator().
214///
215/// So far, only generating mass values is supported. Others will have to be generated
216/// by the slower accept/reject method.
218{
219 if (code == 1) {
220 while (true) {
221 const double gauss = RooRandom::randomGenerator()->Gaus(0., 1.);
222 const double mass = _lambda * sinh((gauss - _gamma)/_delta) + _mu;
223 if (_mass.min() <= mass && mass <= _mass.max() && _massThreshold <= mass) {
224 _mass = mass;
225 break;
226 }
227 }
228 } else {
229 throw std::logic_error("Generation in other variables not yet implemented.");
230 }
231}
#define ClassImp(name)
Definition Rtypes.h:377
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
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
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:55
RooBatchCompute::Config config(RooAbsArg const *arg) const
Definition DataMap.cxx:40
Johnson's distribution.
Definition RooJohnson.h:24
void computeBatch(double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of the Johnson distribution.
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
RooRealProxy _mass
Definition RooJohnson.h:48
RooRealProxy _delta
Definition RooJohnson.h:53
RooRealProxy _gamma
Definition RooJohnson.h:52
void generateEvent(Int_t code) override
Generate events based on code obtained by getGenerator().
RooRealProxy _mu
Definition RooJohnson.h:49
RooRealProxy _lambda
Definition RooJohnson.h:50
double _massThreshold
Definition RooJohnson.h:55
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Advertise which kind of direct event generation is supported.
static TRandom * randomGenerator()
Return a pointer to a singleton random-number generator implementation.
Definition RooRandom.cxx:51
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower 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:275
std::vector< double > ArgVector
void compute(Config cfg, Computer comp, RestrictArr output, size_t size, const VarVector &vars, ArgVector &extraArgs)
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.
constexpr Double_t TwoPi()
Definition TMath.h:44
static void output()