Logo ROOT  
Reference Guide
 
Loading...
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 *
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 "TMath.h"
49
50#include <array>
51#include <cmath>
52
54
55////////////////////////////////////////////////////////////////////////////////
56/// Construct a new Johnson PDF.
57///
58/// \param name Name that identifies the PDF in computations
59/// \param title Title for plotting
60/// \param mass The variable of the PDF. Often this is a mass.
61/// \param mu Location parameter of the Gaussian component.
62/// \param lambda Width parameter (>0) of the Gaussian component.
63/// \param gamma Shape parameter that distorts distribution to left/right.
64/// \param delta Shape parameter (>0) that determines strength of Gaussian-like component.
65/// \param massThreshold Set PDF to zero below this threshold.
66RooJohnson::RooJohnson(const char *name, const char *title,
67 RooAbsReal& mass, RooAbsReal& mu, RooAbsReal& lambda,
68 RooAbsReal& gamma, RooAbsReal& delta,
69 double massThreshold) :
70 RooAbsPdf(name,title),
71 _mass("mass", "Mass observable", this, mass),
72 _mu("mu", "Location parameter of the underlying normal distribution.", this, mu),
73 _lambda("lambda", "Width parameter of the underlying normal distribution (=2 lambda)", this, lambda),
74 _gamma("gamma", "Shift of transformation", this, gamma),
75 _delta("delta", "Scale of transformation", this, delta),
76 _massThreshold(massThreshold)
77{
78 RooHelpers::checkRangeOfParameters(this, {&lambda, &delta}, 0.);
79}
80
81
82////////////////////////////////////////////////////////////////////////////////
83/// Copy a Johnson PDF.
84RooJohnson::RooJohnson(const RooJohnson& other, const char* newName) :
85 RooAbsPdf(other, newName),
86 _mass("Mass", this, other._mass),
87 _mu("mean", this, other._mu),
88 _lambda("lambda", this, other._lambda),
89 _gamma("gamma", this, other._gamma),
90 _delta("delta", this, other._delta),
91 _massThreshold(other._massThreshold)
92{
93
94}
95
96
97////////////////////////////////////////////////////////////////////////////////
98
100{
101 if (_mass < _massThreshold)
102 return 0.;
103
104 const double arg = (_mass-_mu)/_lambda;
105 const double expo = _gamma + _delta * asinh(arg);
106
107 const double result = _delta
108 / sqrt(TMath::TwoPi())
109 / (_lambda * sqrt(1. + arg*arg))
110 * exp(-0.5 * expo * expo);
111
112 return result;
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// Compute multiple values of the Johnson distribution.
117void RooJohnson::computeBatch(double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
118{
119 std::array<double, 1> extraArgs{_massThreshold};
121 {dataMap.at(_mass), dataMap.at(_mu), dataMap.at(_lambda), dataMap.at(_gamma), dataMap.at(_delta)},
122 extraArgs);
123}
124
125////////////////////////////////////////////////////////////////////////////////
126
127int RooJohnson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
128{
129 if (matchArgs(allVars, analVars, _mass)) return kMass;
130 if (matchArgs(allVars, analVars, _mu)) return kMean;
131 if (matchArgs(allVars, analVars, _lambda)) return kLambda;
132 if (matchArgs(allVars, analVars, _gamma)) return kGamma;
133 if (matchArgs(allVars, analVars, _delta)) return kDelta;
134 //TODO write integral for others
135 return 0;
136}
137
138////////////////////////////////////////////////////////////////////////////////
139
140double RooJohnson::analyticalIntegral(Int_t code, const char* rangeName) const
141{
142 //The normalisation constant is left out in evaluate().
143 //Therefore, the integral is scaled up by that amount to make RooFit normalise
144 //correctly.
145 const double globalNorm = 1.;
146// const double globalNorm = sqrt(TMath::TwoPi());
147
148 //Here everything is scaled and shifted such that we only need to compute CDF(Gauss):
149 double min = -1.E300;
150 double max = 1.E300;
151 if (kMass <= code && code <= kLambda) {
152 double argMin;
153 double argMax;
154
155 if (code == kMass) {
156 argMin = (_mass.min(rangeName)-_mu)/_lambda;
157 argMax = (_mass.max(rangeName)-_mu)/_lambda;
158 } else if (code == kMean) {
159 argMin = (_mass-_mu.min(rangeName))/_lambda;
160 argMax = (_mass-_mu.max(rangeName))/_lambda;
161 } else {
162 assert(code == kLambda);
163 argMin = (_mass-_mu)/_lambda.min(rangeName);
164 argMax = (_mass-_mu)/_lambda.max(rangeName);
165 }
166
167 min = _gamma + _delta * asinh(argMin);
168 max = _gamma + _delta * asinh(argMax);
169 } else if (code == kGamma) {
170 const double arg = (_mass-_mu)/_lambda;
171 min = _gamma.min(rangeName) + _delta * asinh(arg);
172 max = _gamma.max(rangeName) + _delta * asinh(arg);
173 } else if (code == kDelta) {
174 const double arg = (_mass-_mu)/_lambda;
175 min = _gamma + _delta.min(rangeName) * asinh(arg);
176 max = _gamma + _delta.max(rangeName) * asinh(arg);
177 } else {
178 assert(false);
179 }
180
181
182
183 //Here we go for maximum precision: We compute all integrals in the UPPER
184 //tail of the Gaussian, because erfc has the highest precision there.
185 //Therefore, the different cases for range limits in the negative hemisphere are mapped onto
186 //the equivalent points in the upper hemisphere using erfc(-x) = 2. - erfc(x)
187 const double ecmin = std::erfc(std::abs(min/sqrt(2.)));
188 const double ecmax = std::erfc(std::abs(max/sqrt(2.)));
189
190 const double result = 0.5 * (
191 min*max < 0.0 ? 2.0 - (ecmin + ecmax)
192 : max <= 0. ? ecmax - ecmin : ecmin - ecmax
193 );
194
195 // Now, include the global norm that may be missing in evaluate and return
196 return globalNorm * (result != 0. ? result : 1.E-300);
197}
198
199
200
201////////////////////////////////////////////////////////////////////////////////
202/// Advertise which kind of direct event generation is supported.
203///
204/// So far, only generating mass values is supported.
205Int_t RooJohnson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
206{
207 if (matchArgs(directVars, generateVars, _mass)) return 1 ;
208// if (matchArgs(directVars, generateVars, _mu)) return 2 ;
209 return 0 ;
210}
211
212
213
214////////////////////////////////////////////////////////////////////////////////
215/// Generate events based on code obtained by getGenerator().
216///
217/// So far, only generating mass values is supported. Others will have to be generated
218/// by the slower accept/reject method.
220{
221 if (code == 1) {
222 while (true) {
223 const double gauss = RooRandom::randomGenerator()->Gaus(0., 1.);
224 const double mass = _lambda * sinh((gauss - _gamma)/_delta) + _mu;
225 if (_mass.min() <= mass && mass <= _mass.max() && _massThreshold <= mass) {
226 _mass = mass;
227 break;
228 }
229 }
230 } else {
231 throw std::logic_error("Generation in other variables not yet implemented.");
232 }
233}
#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:73
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:48
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
void compute(Config cfg, Computer comp, RestrictArr output, size_t size, VarSpan vars, ArgSpan 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()