Logo 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 "RooAbsReal.h"
27#include "RooMath.h"
28#include "RooBatchCompute.h"
29
30#include "TMath.h"
31
32#include <cmath>
33
34using namespace std;
35
37
38////////////////////////////////////////////////////////////////////////////////
39
40RooBifurGauss::RooBifurGauss(const char *name, const char *title,
41 RooAbsReal& _x, RooAbsReal& _mean,
42 RooAbsReal& _sigmaL, RooAbsReal& _sigmaR) :
43 RooAbsPdf(name, title),
44 x ("x" , "Dependent" , this, _x),
45 mean ("mean" , "Mean" , this, _mean),
46 sigmaL("sigmaL", "Left Sigma" , this, _sigmaL),
47 sigmaR("sigmaR", "Right Sigma", this, _sigmaR)
48
49{
50}
51
52////////////////////////////////////////////////////////////////////////////////
53
55 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
56 sigmaL("sigmaL",this,other.sigmaL), sigmaR("sigmaR", this, other.sigmaR)
57{
58}
59
60////////////////////////////////////////////////////////////////////////////////
61
63 double arg = x - mean;
64
65 double coef(0.0);
66
67 if (arg < 0.0){
68 if (TMath::Abs(sigmaL) > 1e-30) {
69 coef = -0.5/(sigmaL*sigmaL);
70 }
71 } else {
72 if (TMath::Abs(sigmaR) > 1e-30) {
73 coef = -0.5/(sigmaR*sigmaR);
74 }
75 }
76
77 return exp(coef*arg*arg);
78}
79
80////////////////////////////////////////////////////////////////////////////////
81/// Compute multiple values of BifurGauss distribution.
82void RooBifurGauss::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooFit::Detail::DataMap const& dataMap) const
83{
85 dispatch->compute(stream, RooBatchCompute::BifurGauss, output, nEvents,
86 {dataMap.at(x),dataMap.at(mean),dataMap.at(sigmaL),dataMap.at(sigmaR)});
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
91Int_t RooBifurGauss::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
92{
93 if (matchArgs(allVars,analVars,x)) return 1 ;
94 return 0 ;
95}
96
97////////////////////////////////////////////////////////////////////////////////
98
99double RooBifurGauss::analyticalIntegral(Int_t code, const char* rangeName) const
100{
101 switch(code) {
102 case 1:
103 {
104 static double root2 = sqrt(2.) ;
105 static double rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
106
107// double coefL(0.0), coefR(0.0);
108// if (TMath::Abs(sigmaL) > 1e-30) {
109// coefL = -0.5/(sigmaL*sigmaL);
110// }
111
112// if (TMath::Abs(sigmaR) > 1e-30) {
113// coefR = -0.5/(sigmaR*sigmaR);
114// }
115
116 double xscaleL = root2*sigmaL;
117 double xscaleR = root2*sigmaR;
118
119 double integral = 0.0;
120 if(x.max(rangeName) < mean)
121 {
122 integral = sigmaL * ( RooMath::erf((x.max(rangeName) - mean)/xscaleL) - RooMath::erf((x.min(rangeName) - mean)/xscaleL) );
123 }
124 else if (x.min(rangeName) > mean)
125 {
126 integral = sigmaR * ( RooMath::erf((x.max(rangeName) - mean)/xscaleR) - RooMath::erf((x.min(rangeName) - mean)/xscaleR) );
127 }
128 else
129 {
130 integral = sigmaR*RooMath::erf((x.max(rangeName) - mean)/xscaleR) - sigmaL*RooMath::erf((x.min(rangeName) - mean)/xscaleL);
131 }
132 // return rootPiBy2*(sigmaR*RooMath::erf((x.max(rangeName) - mean)/xscaleR) -
133 // sigmaL*RooMath::erf((x.min(rangeName) - mean)/xscaleL));
134 return integral*rootPiBy2;
135 }
136 }
137
138 assert(0) ;
139 return 0 ; // to prevent compiler warnings
140}
#define e(i)
Definition: RSha256.hxx:103
#define ClassImp(name)
Definition: Rtypes.h:375
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
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=0) const override
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
double analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
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
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 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.
RVec< PromoteTypes< T0, T1 > > atan2(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1764
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
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()