Logo ROOT  
Reference Guide
RooDstD0BG.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * UE, Ulrik Egede, RAL, U.Egede@rl.ac.uk *
7 * MT, Max Turri, UC Santa Cruz turri@slac.stanford.edu *
8 * CC, Chih-hsiang Cheng, Stanford chcheng@slac.stanford.edu *
9 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
10 * *
11 * Copyright (c) 2000-2005, Regents of the University of California *
12 * RAL and Stanford University. All rights reserved.*
13 * *
14 * Redistribution and use in source and binary forms, *
15 * with or without modification, are permitted according to the terms *
16 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
17 *****************************************************************************/
18
19/** \class RooDstD0BG
20 \ingroup Roofit
21
22Special p.d.f shape that can be used to model the background of
23D*-D0 mass difference distributions
24**/
25
26#include "RooDstD0BG.h"
27#include "RooFit.h"
28#include "RooAbsReal.h"
29#include "RooRealVar.h"
30#include "RooIntegrator1D.h"
31#include "RooAbsFunc.h"
32#include "RooVDTHeaders.h"
33#include "BatchHelpers.h"
34
35#include "TMath.h"
36
37#include <cmath>
38using namespace std;
39
41
42////////////////////////////////////////////////////////////////////////////////
43
44RooDstD0BG::RooDstD0BG(const char *name, const char *title,
45 RooAbsReal& _dm, RooAbsReal& _dm0,
46 RooAbsReal& _c, RooAbsReal& _a, RooAbsReal& _b) :
47 RooAbsPdf(name,title),
48 dm("dm","Dstar-D0 Mass Diff",this, _dm),
49 dm0("dm0","Threshold",this, _dm0),
50 C("C","Shape Parameter",this, _c),
51 A("A","Shape Parameter 2",this, _a),
52 B("B","Shape Parameter 3",this, _b)
53{
54}
55
56////////////////////////////////////////////////////////////////////////////////
57
58RooDstD0BG::RooDstD0BG(const RooDstD0BG& other, const char *name) :
59 RooAbsPdf(other,name), dm("dm",this,other.dm), dm0("dm0",this,other.dm0),
60 C("C",this,other.C), A("A",this,other.A), B("B",this,other.B)
61{
62}
63
64////////////////////////////////////////////////////////////////////////////////
65
67{
68 Double_t arg= dm- dm0;
69 if (arg <= 0 ) return 0;
70 Double_t ratio= dm/dm0;
71 Double_t val= (1- exp(-arg/C))* TMath::Power(ratio, A) + B*(ratio-1);
72
73 return (val > 0 ? val : 0) ;
74}
75
76namespace {
77//Author: Emmanouil Michalainas, CERN 9 SEPTEMBER 2019
78
79template<class Tdm, class Tdm0, class TC, class TA, class TB>
80void compute( size_t batchSize, double * __restrict output,
81 Tdm DM, Tdm0 DM0, TC C, TA A, TB B)
82{
83 for (size_t i=0; i<batchSize; i++) {
84 const double ratio = DM[i] / DM0[i];
85 const double arg1 = (DM0[i]-DM[i]) / C[i];
86 const double arg2 = A[i]*_rf_fast_log(ratio);
87 output[i] = (1 -_rf_fast_exp(arg1)) * _rf_fast_exp(arg2) +B[i]*(ratio-1);
88 }
89
90 for (size_t i=0; i<batchSize; i++) {
91 if (output[i]<0) output[i] = 0;
92 }
93}
94};
95
96RooSpan<double> RooDstD0BG::evaluateBatch(std::size_t begin, std::size_t batchSize) const {
97 using namespace BatchHelpers;
98
99 EvaluateInfo info = getInfo( {&dm, &dm0, &C, &A, &B}, begin, batchSize );
100 if (info.nBatches == 0) {
101 return {};
102 }
103 auto output = _batchData.makeWritableBatchUnInit(begin, batchSize);
104 auto dmData = dm.getValBatch(begin, info.size);
105
106 if (info.nBatches==1 && !dmData.empty()) {
107 compute(info.size, output.data(), dmData.data(),
112 }
113 else {
114 compute(info.size, output.data(),
120 }
121 return output;
122}
123
124////////////////////////////////////////////////////////////////////////////////
125/// if (matchArgs(allVars,analVars,dm)) return 1 ;
126
127Int_t RooDstD0BG::getAnalyticalIntegral(RooArgSet& /*allVars*/, RooArgSet& /*analVars*/, const char* /*rangeName*/) const
128{
129 return 0 ;
130}
131
132////////////////////////////////////////////////////////////////////////////////
133
134Double_t RooDstD0BG::analyticalIntegral(Int_t code, const char* rangeName) const
135{
136 switch(code) {
137 case 1:
138 {
139 Double_t min= dm.min(rangeName);
140 Double_t max= dm.max(rangeName);
141 if (max <= dm0 ) return 0;
142 else if (min < dm0) min = dm0;
143
144 Bool_t doNumerical= kFALSE;
145 if ( A != 0 ) doNumerical= kTRUE;
146 else if (B < 0) {
147 // If b<0, pdf can be negative at large dm, the integral should
148 // only up to where pdf hits zero. Better solution should be
149 // solve the zero and use it as max.
150 // Here we check this whether pdf(max) < 0. If true, let numerical
151 // integral take care of. ( kind of ugly!)
152 if ( 1- exp(-(max-dm0)/C) + B*(max/dm0 -1) < 0) doNumerical= kTRUE;
153 }
154 if ( ! doNumerical ) {
155 return (max-min)+ C* exp(dm0/C)* (exp(-max/C)- exp(-min/C)) +
156 B * (0.5* (max*max - min*min)/dm0 - (max- min));
157 } else {
158 // In principle the integral for a!=0 can be done analytically.
159 // It involves incomplete Gamma function, TMath::Gamma(a+1,m/c),
160 // which is not defined for a < -1. And the whole expression is
161 // not stable for m/c >> 1.
162 // Do numerical integral
163 RooArgSet vset(dm.arg(),"vset");
164 RooAbsFunc *func= bindVars(vset);
165 RooIntegrator1D integrator(*func,min,max);
166 return integrator.integral();
167 }
168 }
169 }
170
171 assert(0) ;
172 return 0 ;
173}
double _rf_fast_exp(double x)
VDT headers for RooFit.
Definition: RooVDTHeaders.h:47
double _rf_fast_log(double x)
Definition: RooVDTHeaders.h:51
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
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
Abstract interface for evaluating a real-valued function of one real variable and performing numerica...
Definition: RooAbsFunc.h:23
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
BatchHelpers::BatchData _batchData
Definition: RooAbsReal.h:447
RooAbsFunc * bindVars(const RooArgSet &vars, const RooArgSet *nset=0, Bool_t clipInvalid=kFALSE) const
Create an interface adaptor f(vars) that binds us to the specified variables (in arbitrary order).
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Special p.d.f shape that can be used to model the background of D*-D0 mass difference distributions.
Definition: RooDstD0BG.h:26
RooRealProxy B
Definition: RooDstD0BG.h:45
RooRealProxy dm
Definition: RooDstD0BG.h:43
RooRealProxy A
Definition: RooDstD0BG.h:45
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooDstD0BG.cxx:134
RooRealProxy dm0
Definition: RooDstD0BG.h:44
RooSpan< double > evaluateBatch(std::size_t begin, std::size_t batchSize) const
Evaluate function for a batch of input data points.
Definition: RooDstD0BG.cxx:96
RooRealProxy C
Definition: RooDstD0BG.h:45
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
if (matchArgs(allVars,analVars,dm)) return 1 ;
Definition: RooDstD0BG.cxx:127
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
Definition: RooDstD0BG.cxx:66
RooIntegrator1D implements an adaptive one-dimensional numerical integration algorithm.
virtual Double_t integral(const Double_t *yvec=0)
Calculate numeric integral at given set of function binding parameters.
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
const T & arg() const
Return reference to object held in proxy.
Double_t max(const char *rname=0) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
EvaluateInfo getInfo(std::vector< const RooRealProxy * > parameters, size_t begin, size_t batchSize)
static double B[]
static double A[]
static double C[]
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
static void output(int code)
Definition: gifencode.c:226