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#include "RooFit.h"
17
18/** \class RooBifurGauss
19 \ingroup Roofit
20
21Bifurcated Gaussian p.d.f with different widths on left and right
22side of maximum value
23**/
24
25#include "RooBifurGauss.h"
26
27#include "RooAbsReal.h"
28#include "RooMath.h"
29#include "RooBatchCompute.h"
30
31#include "TMath.h"
32
33#include <cmath>
34
35using namespace std;
36
38
39////////////////////////////////////////////////////////////////////////////////
40
41RooBifurGauss::RooBifurGauss(const char *name, const char *title,
42 RooAbsReal& _x, RooAbsReal& _mean,
43 RooAbsReal& _sigmaL, RooAbsReal& _sigmaR) :
44 RooAbsPdf(name, title),
45 x ("x" , "Dependent" , this, _x),
46 mean ("mean" , "Mean" , this, _mean),
47 sigmaL("sigmaL", "Left Sigma" , this, _sigmaL),
48 sigmaR("sigmaR", "Right Sigma", this, _sigmaR)
49
50{
51}
52
53////////////////////////////////////////////////////////////////////////////////
54
56 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
57 sigmaL("sigmaL",this,other.sigmaL), sigmaR("sigmaR", this, other.sigmaR)
58{
59}
60
61////////////////////////////////////////////////////////////////////////////////
62
64 Double_t arg = x - mean;
65
66 Double_t coef(0.0);
67
68 if (arg < 0.0){
69 if (TMath::Abs(sigmaL) > 1e-30) {
70 coef = -0.5/(sigmaL*sigmaL);
71 }
72 } else {
73 if (TMath::Abs(sigmaR) > 1e-30) {
74 coef = -0.5/(sigmaR*sigmaR);
75 }
76 }
77
78 return exp(coef*arg*arg);
79}
80
81////////////////////////////////////////////////////////////////////////////////
82/// Compute multiple values of BifurGauss distribution.
83void RooBifurGauss::computeBatch(cudaStream_t* stream, double* output, size_t nEvents, RooBatchCompute::DataMap& dataMap) const
84{
86 dispatch->compute(stream, RooBatchCompute::BifurGauss, output, nEvents, dataMap, {&*x,&*mean,&*sigmaL,&*sigmaR,&*_norm});
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_t RooBifurGauss::analyticalIntegral(Int_t code, const char* rangeName) const
100{
101 switch(code) {
102 case 1:
103 {
104 static Double_t root2 = sqrt(2.) ;
105 static Double_t rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
106
107// Double_t 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_t xscaleL = root2*sigmaL;
117 Double_t xscaleR = root2*sigmaR;
118
119 Double_t 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:364
char name[80]
Definition: TGX11.cxx:110
RooAbsReal * _norm
Definition: RooAbsPdf.h:364
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
Bool_t 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:35
virtual void compute(cudaStream_t *, Computer, RestrictArr, size_t, const DataMap &, 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
RooRealProxy mean
Definition: RooBifurGauss.h:41
RooRealProxy sigmaL
Definition: RooBifurGauss.h:42
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
void computeBatch(cudaStream_t *, double *output, size_t nEvents, RooBatchCompute::DataMap &) const
Compute multiple values of BifurGauss distribution.
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooRealProxy sigmaR
Definition: RooBifurGauss.h:43
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
RooRealProxy x
Definition: RooBifurGauss.h:40
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:575
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.
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
std::map< DataKey, RooSpan< const double > > DataMap
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)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226