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 *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
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 "BatchHelpers.h"
31
32#include "TMath.h"
33
34#include <cmath>
35
36using namespace std;
37
39
40////////////////////////////////////////////////////////////////////////////////
41
42RooBifurGauss::RooBifurGauss(const char *name, const char *title,
43 RooAbsReal& _x, RooAbsReal& _mean,
44 RooAbsReal& _sigmaL, RooAbsReal& _sigmaR) :
45 RooAbsPdf(name, title),
46 x ("x" , "Dependent" , this, _x),
47 mean ("mean" , "Mean" , this, _mean),
48 sigmaL("sigmaL", "Left Sigma" , this, _sigmaL),
49 sigmaR("sigmaR", "Right Sigma", this, _sigmaR)
50
51{
52}
53
54////////////////////////////////////////////////////////////////////////////////
55
57 RooAbsPdf(other,name), x("x",this,other.x), mean("mean",this,other.mean),
58 sigmaL("sigmaL",this,other.sigmaL), sigmaR("sigmaR", this, other.sigmaR)
59{
60}
61
62////////////////////////////////////////////////////////////////////////////////
63
65 Double_t arg = x - mean;
66
67 Double_t coef(0.0);
68
69 if (arg < 0.0){
70 if (TMath::Abs(sigmaL) > 1e-30) {
71 coef = -0.5/(sigmaL*sigmaL);
72 }
73 } else {
74 if (TMath::Abs(sigmaR) > 1e-30) {
75 coef = -0.5/(sigmaR*sigmaR);
76 }
77 }
78
79 return exp(coef*arg*arg);
80}
81
82////////////////////////////////////////////////////////////////////////////////
83
84namespace {
85//Author: Emmanouil Michalainas, CERN 20 AUGUST 2019
86
87template<class Tx, class Tm, class Tsl, class Tsr>
88void compute( size_t batchSize,
89 double * __restrict output,
90 Tx X, Tm M, Tsl SL, Tsr SR)
91{
92 for (size_t i=0; i<batchSize; i++) {
93 const double arg = X[i]-M[i];
94 output[i] = arg / ((arg < 0.0)*SL[i] + (arg >= 0.0)*SR[i]);
95 }
96
97 for (size_t i=0; i<batchSize; i++) {
98 if (X[i]-M[i]>1e-30 || X[i]-M[i]<-1e-30) {
99 output[i] = _rf_fast_exp(-0.5*output[i]*output[i]);
100 }
101 else {
102 output[i] = 1.0;
103 }
104 }
105}
106};
107
108RooSpan<double> RooBifurGauss::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
109 using namespace BatchHelpers;
110
111 EvaluateInfo info = getInfo( {&x, &mean, &sigmaL, &sigmaR}, begin, batchSize );
112 if (info.nBatches == 0) {
113 return {};
114 }
115 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
116 auto xData = x.getValBatch(begin, info.size);
117
118 if (info.nBatches==1 && !xData.empty()) {
119 compute(info.size, output.data(), xData.data(),
123 }
124 else {
125 compute(info.size, output.data(),
130 }
131 return output;
132}
133
134
135////////////////////////////////////////////////////////////////////////////////
136
137Int_t RooBifurGauss::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
138{
139 if (matchArgs(allVars,analVars,x)) return 1 ;
140 return 0 ;
141}
142
143////////////////////////////////////////////////////////////////////////////////
144
145Double_t RooBifurGauss::analyticalIntegral(Int_t code, const char* rangeName) const
146{
147 switch(code) {
148 case 1:
149 {
150 static Double_t root2 = sqrt(2.) ;
151 static Double_t rootPiBy2 = sqrt(atan2(0.0,-1.0)/2.0);
152
153// Double_t coefL(0.0), coefR(0.0);
154// if (TMath::Abs(sigmaL) > 1e-30) {
155// coefL = -0.5/(sigmaL*sigmaL);
156// }
157
158// if (TMath::Abs(sigmaR) > 1e-30) {
159// coefR = -0.5/(sigmaR*sigmaR);
160// }
161
162 Double_t xscaleL = root2*sigmaL;
163 Double_t xscaleR = root2*sigmaR;
164
165 Double_t integral = 0.0;
166 if(x.max(rangeName) < mean)
167 {
168 integral = sigmaL * ( RooMath::erf((x.max(rangeName) - mean)/xscaleL) - RooMath::erf((x.min(rangeName) - mean)/xscaleL) );
169 }
170 else if (x.min(rangeName) > mean)
171 {
172 integral = sigmaR * ( RooMath::erf((x.max(rangeName) - mean)/xscaleR) - RooMath::erf((x.min(rangeName) - mean)/xscaleR) );
173 }
174 else
175 {
176 integral = sigmaR*RooMath::erf((x.max(rangeName) - mean)/xscaleR) - sigmaL*RooMath::erf((x.min(rangeName) - mean)/xscaleL);
177 }
178 // return rootPiBy2*(sigmaR*RooMath::erf((x.max(rangeName) - mean)/xscaleR) -
179 // sigmaL*RooMath::erf((x.min(rangeName) - mean)/xscaleL));
180 return integral*rootPiBy2;
181 }
182 }
183
184 assert(0) ;
185 return 0 ; // to prevent compiler warnings
186}
#define e(i)
Definition: RSha256.hxx:103
double _rf_fast_exp(double x)
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
double atan2(double, double)
double sqrt(double)
double exp(double)
RooSpan< double > makeWritableBatchUnInit(std::size_t begin, std::size_t batchSize)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:58
Little adapter that gives a bracket operator to types that don't have one.
Definition: BatchHelpers.h:58
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
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.
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
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:580
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
Double_t min(const char *rname=0) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooSpan< const double > getValBatch(std::size_t begin, std::size_t batchSize) const
Double_t 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
EvaluateInfo getInfo(std::vector< const RooRealProxy * > parameters, size_t begin, size_t batchSize)
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226