ROOT   Reference Guide
RooBifurGauss.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * Abi Soffer, Colorado State University, abi@slac.stanford.edu *
7 * *
8 * Copyright (c) 2000-2005, Regents of the University of California, *
9 * Colorado State University *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/** \class RooBifurGauss
18 \ingroup Roofit
19
20Bifurcated Gaussian p.d.f with different widths on left and right
21side of maximum value
22**/
23
24#include "RooBifurGauss.h"
25
26#include "RooMath.h"
27#include "RooBatchCompute.h"
28
29#include "TMath.h"
30
31#include <cmath>
32
33using namespace std;
34
36
37////////////////////////////////////////////////////////////////////////////////
38
39RooBifurGauss::RooBifurGauss(const char *name, const char *title,
40 RooAbsReal& _x, RooAbsReal& _mean,
41 RooAbsReal& _sigmaL, RooAbsReal& _sigmaR) :
42 RooAbsPdf(name, title),
43 x ("x" , "Dependent" , this, _x),
44 mean ("mean" , "Mean" , this, _mean),
45 sigmaL("sigmaL", "Left Sigma" , this, _sigmaL),
46 sigmaR("sigmaR", "Right Sigma", this, _sigmaR)
47
48{
49}
50
51////////////////////////////////////////////////////////////////////////////////
52
54 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
55 sigmaL("sigmaL",this,other.sigmaL), sigmaR("sigmaR", this, other.sigmaR)
56{
57}
58
59////////////////////////////////////////////////////////////////////////////////
60
62 double arg = x - mean;
63
64 double coef(0.0);
65
66 if (arg < 0.0){
67 if (TMath::Abs(sigmaL) > 1e-30) {
68 coef = -0.5/(sigmaL*sigmaL);
69 }
70 } else {
71 if (TMath::Abs(sigmaR) > 1e-30) {
72 coef = -0.5/(sigmaR*sigmaR);
73 }
74 }
75
76 return exp(coef*arg*arg);
77}
78
79////////////////////////////////////////////////////////////////////////////////
80/// Compute multiple values of BifurGauss distribution.
81void RooBifurGauss::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
82{
84 dispatch->compute(stream, RooBatchCompute::BifurGauss, output, nEvents,
85 {dataMap.at(x),dataMap.at(mean),dataMap.at(sigmaL),dataMap.at(sigmaR)});
86}
87
88////////////////////////////////////////////////////////////////////////////////
89
90Int_t RooBifurGauss::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
91{
92 if (matchArgs(allVars,analVars,x)) return 1 ;
93 return 0 ;
94}
95
96////////////////////////////////////////////////////////////////////////////////
97
98double RooBifurGauss::analyticalIntegral(Int_t code, const char* rangeName) const
99{
100 switch(code) {
101 case 1:
102 {
103 static double root2 = sqrt(2.) ;
104 static double rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
105
106// double coefL(0.0), coefR(0.0);
107// if (TMath::Abs(sigmaL) > 1e-30) {
108// coefL = -0.5/(sigmaL*sigmaL);
109// }
110
111// if (TMath::Abs(sigmaR) > 1e-30) {
112// coefR = -0.5/(sigmaR*sigmaR);
113// }
114
115 double xscaleL = root2*sigmaL;
116 double xscaleR = root2*sigmaR;
117
118 double integral = 0.0;
119 if(x.max(rangeName) < mean)
120 {
121 integral = sigmaL * ( RooMath::erf((x.max(rangeName) - mean)/xscaleL) - RooMath::erf((x.min(rangeName) - mean)/xscaleL) );
122 }
123 else if (x.min(rangeName) > mean)
124 {
125 integral = sigmaR * ( RooMath::erf((x.max(rangeName) - mean)/xscaleR) - RooMath::erf((x.min(rangeName) - mean)/xscaleR) );
126 }
127 else
128 {
129 integral = sigmaR*RooMath::erf((x.max(rangeName) - mean)/xscaleR) - sigmaL*RooMath::erf((x.min(rangeName) - mean)/xscaleL);
130 }
131 // return rootPiBy2*(sigmaR*RooMath::erf((x.max(rangeName) - mean)/xscaleR) -
132 // sigmaL*RooMath::erf((x.min(rangeName) - mean)/xscaleL));
133 return integral*rootPiBy2;
134 }
135 }
136
137 assert(0) ;
138 return 0 ; // to prevent compiler warnings
139}
#define e(i)
Definition: RSha256.hxx:103
ClassImp(RooBifurGauss)
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:62
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:56
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const VarVector &, const ArgVector &={})=0
Bifurcated Gaussian p.d.f with different widths on left and right side of maximum value.
Definition: RooBifurGauss.h:24
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
RooRealProxy mean
Definition: RooBifurGauss.h:41
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooFit::Detail::DataMap const &) const override
Compute multiple values of BifurGauss distribution.
RooRealProxy sigmaL
Definition: RooBifurGauss.h:42
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooRealProxy sigmaR
Definition: RooBifurGauss.h:43
RooRealProxy x
Definition: RooBifurGauss.h:40
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:60
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.
RVec< PromoteTypes< T0, T1 > > atan2(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1781
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1761
Double_t x[n]
Definition: legend1.C:17
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...
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
static void output()