Logo ROOT  
Reference Guide
RooArgusBG.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
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 RooArgusBG
18 \ingroup Roofit
19
20RooArgusBG is a RooAbsPdf implementation describing the ARGUS background shape.
21\f[
22 \mathrm{Argus}(m, m_0, c, p) = \mathcal{N} \cdot m \cdot \left[ 1 - \left( \frac{m}{m_0} \right)^2 \right]^p
23 \cdot \exp\left[ c \cdot \left(1 - \left(\frac{m}{m_0}\right)^2 \right) \right]
24\f]
25\image html RooArgusBG.png
26*/
27
28#include "RooArgusBG.h"
29#include "RooRealVar.h"
30#include "RooRealConstant.h"
31#include "RooMath.h"
32#include "BatchHelpers.h"
33#include "RooVDTHeaders.h"
34
35#include "TMath.h"
36
37#include <cmath>
38using namespace std;
39
41
42////////////////////////////////////////////////////////////////////////////////
43/// Constructor.
44
45RooArgusBG::RooArgusBG(const char *name, const char *title,
46 RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _c) :
47 RooAbsPdf(name, title),
48 m("m","Mass",this,_m),
49 m0("m0","Resonance mass",this,_m0),
50 c("c","Slope parameter",this,_c),
51 p("p","Power",this,(RooRealVar&)RooRealConstant::value(0.5))
52{
53}
54
55////////////////////////////////////////////////////////////////////////////////
56/// Constructor.
57
58RooArgusBG::RooArgusBG(const char *name, const char *title,
59 RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _c, RooAbsReal& _p) :
60 RooAbsPdf(name, title),
61 m("m","Mass",this,_m),
62 m0("m0","Resonance mass",this,_m0),
63 c("c","Slope parameter",this,_c),
64 p("p","Power",this,_p)
65{
66}
67
68////////////////////////////////////////////////////////////////////////////////
69/// Constructor.
70
71RooArgusBG::RooArgusBG(const RooArgusBG& other, const char* name) :
72 RooAbsPdf(other,name),
73 m("m",this,other.m),
74 m0("m0",this,other.m0),
75 c("c",this,other.c),
76 p("p",this,other.p)
77{
78}
79
80////////////////////////////////////////////////////////////////////////////////
81
83 Double_t t= m/m0;
84 if(t >= 1) return 0;
85
86 Double_t u= 1 - t*t;
87 //cout << "c = " << c << " result = " << m*TMath::Power(u,p)*exp(c*u) << endl ;
88 return m*TMath::Power(u,p)*exp(c*u) ;
89}
90
91////////////////////////////////////////////////////////////////////////////////
92
93namespace {
94//Author: Emmanouil Michalainas, CERN 19 AUGUST 2019
95
96template<class Tm, class Tm0, class Tc, class Tp>
97void compute( size_t batchSize,
98 double * __restrict output,
99 Tm M, Tm0 M0, Tc C, Tp P)
100{
101 for (size_t i=0; i<batchSize; i++) {
102 const double t = M[i]/M0[i];
103 const double u = 1 - t*t;
104 output[i] = C[i]*u + P[i]*_rf_fast_log(u);
105 }
106 for (size_t i=0; i<batchSize; i++) {
107 if (M[i] >= M0[i]) output[i] = 0.0;
108 else output[i] = M[i]*_rf_fast_exp(output[i]);
109 }
110}
111};
112
113RooSpan<double> RooArgusBG::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
114 using namespace BatchHelpers;
115
116 EvaluateInfo info = getInfo( {&m, &m0, &c, &p}, begin, batchSize );
117 if (info.nBatches == 0) {
118 return {};
119 }
120 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
121 auto mData = m.getValBatch(begin, info.size);
122
123 if (info.nBatches==1 && !mData.empty()) {
124 compute(info.size, output.data(), mData.data(),
128 }
129 else {
130 compute(info.size, output.data(),
135 }
136 return output;
137}
138
139////////////////////////////////////////////////////////////////////////////////
140
141Int_t RooArgusBG::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
142{
143 if (p.arg().isConstant()) {
144 // We can integrate over m if power = 0.5
145 if (matchArgs(allVars,analVars,m) && p == 0.5) return 1;
146 }
147 return 0;
148
149}
150
151////////////////////////////////////////////////////////////////////////////////
152
153Double_t RooArgusBG::analyticalIntegral(Int_t code, const char* rangeName) const
154{
155 R__ASSERT(code==1);
156 // Formula for integration over m when p=0.5
157 static const Double_t pi = atan2(0.0,-1.0);
158 Double_t min = (m.min(rangeName) < m0) ? m.min(rangeName) : m0;
159 Double_t max = (m.max(rangeName) < m0) ? m.max(rangeName) : m0;
160 Double_t f1 = (1.-TMath::Power(min/m0,2));
161 Double_t f2 = (1.-TMath::Power(max/m0,2));
162 Double_t aLow, aHigh ;
163 if ( c < 0. ) {
164 aLow = -0.5*m0*m0*(exp(c*f1)*sqrt(f1)/c + 0.5/TMath::Power(-c,1.5)*sqrt(pi)*RooMath::erf(sqrt(-c*f1)));
165 aHigh = -0.5*m0*m0*(exp(c*f2)*sqrt(f2)/c + 0.5/TMath::Power(-c,1.5)*sqrt(pi)*RooMath::erf(sqrt(-c*f2)));
166 } else if ( c == 0. ) {
167 aLow = -m0*m0/3.*f1*sqrt(f1);
168 aHigh = -m0*m0/3.*f1*sqrt(f2);
169 } else {
170 aLow = 0.5*m0*m0*exp(c*f1)/(c*sqrt(c)) * (0.5*sqrt(pi)*(RooMath::faddeeva(sqrt(c*f1))).imag() - sqrt(c*f1));
171 aHigh = 0.5*m0*m0*exp(c*f2)/(c*sqrt(c)) * (0.5*sqrt(pi)*(RooMath::faddeeva(sqrt(c*f2))).imag() - sqrt(c*f2));
172 }
173 Double_t area = aHigh - aLow;
174 //cout << "c = " << c << "aHigh = " << aHigh << " aLow = " << aLow << " area = " << area << endl ;
175 return area;
176
177}
#define c(i)
Definition: RSha256.hxx:101
double _rf_fast_exp(double x)
VDT headers for RooFit.
Definition: RooVDTHeaders.h:47
double _rf_fast_log(double x)
Definition: RooVDTHeaders.h:51
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
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, const RooArgSet *const normSet=nullptr, Tag_t ownerTag=kUnspecified)
Make a batch and return a span pointing to the pdf-local memory.
Definition: BatchData.cxx:118
Little adapter that gives a bracket operator to types that don't have one.
Definition: BatchHelpers.h:58
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:342
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:450
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooArgusBG is a RooAbsPdf implementation describing the ARGUS background shape.
Definition: RooArgusBG.h:25
RooRealProxy m
Definition: RooArgusBG.h:40
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooArgusBG.cxx:82
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooArgusBG.cxx:141
RooRealProxy c
Definition: RooArgusBG.h:42
RooRealProxy p
Definition: RooArgusBG.h:43
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooArgusBG.cxx:153
RooRealProxy m0
Definition: RooArgusBG.h:41
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
Definition: RooArgusBG.cxx:113
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:580
static std::complex< double > faddeeva(std::complex< double > z)
evaluate Faddeeva function for complex argument
Definition: RooMath.cxx:542
RooRealConstant provides static functions to create and keep track of RooRealVar constants.
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
A simple container to hold a batch of data values.
Definition: RooSpan.h:32
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.
RooSpan< const typename T::value_type > getValBatch(std::size_t begin, std::size_t batchSize) const
Retrieve a batch of real or category data.
const T & arg() const
Return reference to object held in proxy.
TF1 * f1
Definition: legend1.C:11
EvaluateInfo getInfo(std::vector< const RooRealProxy * > parameters, size_t begin, size_t batchSize)
static double P[]
static double C[]
static constexpr double pi
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
auto * m
Definition: textangle.C:8
static void output(int code)
Definition: gifencode.c:226