Logo ROOT  
Reference Guide
RooTruthModel.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$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/**
18\file RooTruthModel.cxx
19\class RooTruthModel
20\ingroup Roofitcore
21
22RooTruthModel is an implementation of RooResolution
23model that provides a delta-function resolution model.
24The truth model supports <i>all</i> basis functions because it evaluates each basis function as
25as a RooFormulaVar. The 6 basis functions used in B mixing and decay and 2 basis
26functions used in D mixing have been hand coded for increased execution speed.
27**/
28
29#include "Riostream.h"
30#include "RooBatchCompute.h"
31#include "RooTruthModel.h"
32#include "RooGenContext.h"
33#include "RooAbsAnaConvPdf.h"
34
35#include "TError.h"
36
37#include <algorithm>
38using namespace std ;
39
41;
42
43
44
45////////////////////////////////////////////////////////////////////////////////
46/// Constructor of a truth resolution model, i.e. a delta function in observable 'xIn'
47
48RooTruthModel::RooTruthModel(const char *name, const char *title, RooAbsRealLValue& xIn) :
49 RooResolutionModel(name,title,xIn)
50{
51}
52
53
54
55////////////////////////////////////////////////////////////////////////////////
56/// Copy constructor
57
60{
61}
62
63
64
65////////////////////////////////////////////////////////////////////////////////
66/// Destructor
67
69{
70}
71
72
73
74////////////////////////////////////////////////////////////////////////////////
75/// Return basis code for given basis definition string. Return special
76/// codes for 'known' bases for which compiled definition exists. Return
77/// generic bases code if implementation relies on TFormula interpretation
78/// of basis name
79
81{
82 std::string str = name;
83
84 // Remove whitespaces from the input string
85 str.erase(remove(str.begin(),str.end(),' '),str.end());
86
87 // Check for optimized basis functions
88 if (str == "exp(-@0/@1)") return expBasisPlus ;
89 if (str == "exp(@0/@1)") return expBasisMinus ;
90 if (str == "exp(-abs(@0)/@1)") return expBasisSum ;
91 if (str == "exp(-@0/@1)*sin(@0*@2)") return sinBasisPlus ;
92 if (str == "exp(@0/@1)*sin(@0*@2)") return sinBasisMinus ;
93 if (str == "exp(-abs(@0)/@1)*sin(@0*@2)") return sinBasisSum ;
94 if (str == "exp(-@0/@1)*cos(@0*@2)") return cosBasisPlus ;
95 if (str == "exp(@0/@1)*cos(@0*@2)") return cosBasisMinus ;
96 if (str == "exp(-abs(@0)/@1)*cos(@0*@2)") return cosBasisSum ;
97 if (str == "(@0/@1)*exp(-@0/@1)") return linBasisPlus ;
98 if (str == "(@0/@1)*(@0/@1)*exp(-@0/@1)") return quadBasisPlus ;
99 if (str == "exp(-@0/@1)*cosh(@0*@2/2)") return coshBasisPlus;
100 if (str == "exp(@0/@1)*cosh(@0*@2/2)") return coshBasisMinus;
101 if (str == "exp(-abs(@0)/@1)*cosh(@0*@2/2)") return coshBasisSum;
102 if (str == "exp(-@0/@1)*sinh(@0*@2/2)") return sinhBasisPlus;
103 if (str == "exp(@0/@1)*sinh(@0*@2/2)") return sinhBasisMinus;
104 if (str == "exp(-abs(@0)/@1)*sinh(@0*@2/2)") return sinhBasisSum;
105
106 // Truth model is delta function, i.e. convolution integral is basis
107 // function, therefore we can handle any basis function
108 return genericBasis ;
109}
110
111
112
113////////////////////////////////////////////////////////////////////////////////
114/// Changes associated bases function to 'inBasis'
115
117{
118 // Remove client-server link to old basis
119 if (_basis) {
120 if (_basisCode == genericBasis) {
121 // In the case of a generic basis, we evaluate it directly, so the
122 // basis was a direct server.
124 } else {
125 for (RooAbsArg *basisServer : _basis->servers()) {
126 removeServer(*basisServer);
127 }
128 }
129
130 if (_ownBasis) {
131 delete _basis;
132 }
133 }
134 _ownBasis = false;
135
136 _basisCode = inBasis ? basisCode(inBasis->GetTitle()) : 0;
137
138 // Change basis pointer and update client-server link
139 _basis = inBasis;
140 if (_basis) {
141 if (_basisCode == genericBasis) {
142 // Since we actually evaluate the basis function object, we need to
143 // adjust our client-server links to the basis function here
144 addServer(*_basis, true, false);
145 } else {
146 for (RooAbsArg *basisServer : _basis->servers()) {
147 addServer(*basisServer, true, false);
148 }
149 }
150 }
151}
152
153
154
155////////////////////////////////////////////////////////////////////////////////
156/// Evaluate the truth model: a delta function when used as PDF,
157/// the basis function itself, when convoluted with a basis function.
158
160{
161 // No basis: delta function
162 if (_basisCode == noBasis) {
163 if (x==0) return 1 ;
164 return 0 ;
165 }
166
167 // Generic basis: evaluate basis function object
168 if (_basisCode == genericBasis) {
169 return basis().getVal() ;
170 }
171
172 // Precompiled basis functions
173 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
174 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
175
176 // Enforce sign compatibility
177 if ((basisSign==Minus && x>0) ||
178 (basisSign==Plus && x<0)) return 0 ;
179
180
181 double tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ;
182 // Return desired basis function
183 switch(basisType) {
184 case expBasis: {
185 //cout << " RooTruthModel::eval(" << GetName() << ") expBasis mode ret = " << exp(-std::abs((double)x)/tau) << " tau = " << tau << endl ;
186 return exp(-std::abs((double)x)/tau) ;
187 }
188 case sinBasis: {
189 double dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
190 return exp(-std::abs((double)x)/tau)*sin(x*dm) ;
191 }
192 case cosBasis: {
193 double dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
194 return exp(-std::abs((double)x)/tau)*cos(x*dm) ;
195 }
196 case linBasis: {
197 double tscaled = std::abs((double)x)/tau;
198 return exp(-tscaled)*tscaled ;
199 }
200 case quadBasis: {
201 double tscaled = std::abs((double)x)/tau;
202 return exp(-tscaled)*tscaled*tscaled;
203 }
204 case sinhBasis: {
205 double dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
206 return exp(-std::abs((double)x)/tau)*sinh(x*dg/2) ;
207 }
208 case coshBasis: {
209 double dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
210 return exp(-std::abs((double)x)/tau)*cosh(x*dg/2) ;
211 }
212 default:
213 R__ASSERT(0) ;
214 }
215
216 return 0 ;
217}
218
219
220void RooTruthModel::computeBatch(cudaStream_t *stream, double *output, size_t nEvents,
221 RooFit::Detail::DataMap const &dataMap) const
222{
224
225 auto xVals = dataMap.at(x);
226
227 // No basis: delta function
228 if (_basisCode == noBasis) {
229 dispatch->compute(stream, RooBatchCompute::DeltaFunction, output, nEvents, {xVals});
230 return;
231 }
232
233 // Generic basis: evaluate basis function object
234 if (_basisCode == genericBasis) {
235 dispatch->compute(stream, RooBatchCompute::Identity, output, nEvents, {dataMap.at(&basis())});
236 return;
237 }
238
239 // Precompiled basis functions
240 const BasisType basisType = static_cast<BasisType>((_basisCode == 0) ? 0 : (_basisCode / 10) + 1);
241
242 // Cast the int from the enum to double because we can only pass doubles to
243 // RooBatchCompute at this point.
244 const double basisSign = static_cast<double>((BasisSign)(_basisCode - 10 * (basisType - 1) - 2));
245
246 auto param1 = static_cast<RooAbsReal const *>(basis().getParameter(1));
247 auto param2 = static_cast<RooAbsReal const *>(basis().getParameter(2));
248 auto param1Vals = param1 ? dataMap.at(param1) : RooSpan<const double>{};
249 auto param2Vals = param2 ? dataMap.at(param2) : RooSpan<const double>{};
250
251 // Return desired basis function
252 switch (basisType) {
253 case expBasis: {
254 dispatch->compute(stream, RooBatchCompute::TruthModelExpBasis, output, nEvents, {xVals, param1Vals}, {basisSign});
255 break;
256 }
257 case sinBasis: {
258 dispatch->compute(stream, RooBatchCompute::TruthModelSinBasis, output, nEvents, {xVals, param1Vals, param2Vals},
259 {basisSign});
260 break;
261 }
262 case cosBasis: {
263 dispatch->compute(stream, RooBatchCompute::TruthModelCosBasis, output, nEvents, {xVals, param1Vals, param2Vals},
264 {basisSign});
265 break;
266 }
267 case linBasis: {
268 dispatch->compute(stream, RooBatchCompute::TruthModelLinBasis, output, nEvents, {xVals, param1Vals}, {basisSign});
269 break;
270 }
271 case quadBasis: {
272 dispatch->compute(stream, RooBatchCompute::TruthModelQuadBasis, output, nEvents, {xVals, param1Vals},
273 {basisSign});
274 break;
275 }
276 case sinhBasis: {
277 dispatch->compute(stream, RooBatchCompute::TruthModelSinhBasis, output, nEvents, {xVals, param1Vals, param2Vals},
278 {basisSign});
279 break;
280 }
281 case coshBasis: {
282 dispatch->compute(stream, RooBatchCompute::TruthModelCoshBasis, output, nEvents, {xVals, param1Vals, param2Vals},
283 {basisSign});
284 break;
285 }
286 default: R__ASSERT(0);
287 }
288}
289
290
291////////////////////////////////////////////////////////////////////////////////
292/// Advertise analytical integrals for compiled basis functions and when used
293/// as p.d.f without basis function.
294
295Int_t RooTruthModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
296{
297 switch(_basisCode) {
298
299 // Analytical integration capability of raw PDF
300 case noBasis:
301 if (matchArgs(allVars,analVars,convVar())) return 1 ;
302 break ;
303
304 // Analytical integration capability of convoluted PDF
305 case expBasisPlus:
306 case expBasisMinus:
307 case expBasisSum:
308 case sinBasisPlus:
309 case sinBasisMinus:
310 case sinBasisSum:
311 case cosBasisPlus:
312 case cosBasisMinus:
313 case cosBasisSum:
314 case linBasisPlus:
315 case quadBasisPlus:
316 case sinhBasisPlus:
317 case sinhBasisMinus:
318 case sinhBasisSum:
319 case coshBasisPlus:
320 case coshBasisMinus:
321 case coshBasisSum:
322 if (matchArgs(allVars,analVars,convVar())) return 1 ;
323 break ;
324 }
325
326 return 0 ;
327}
328
329
330
331////////////////////////////////////////////////////////////////////////////////
332/// Implement analytical integrals when used as p.d.f and for compiled
333/// basis functions.
334
335double RooTruthModel::analyticalIntegral(Int_t code, const char* rangeName) const
336{
337
338 // Code must be 1
339 R__ASSERT(code==1) ;
340
341 // Unconvoluted PDF
342 if (_basisCode==noBasis) return 1 ;
343
344 // Precompiled basis functions
345 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
346 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
347 //cout << " calling RooTruthModel::analyticalIntegral with basisType " << basisType << endl;
348
349 double tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ;
350
351 const double xmin = x.min(rangeName);
352 const double xmax = x.max(rangeName);
353
354 switch (basisType) {
355 case expBasis:
356 {
357 // WVE fixed for ranges
358 double result(0) ;
359 if (tau==0) return 1 ;
360 if ((basisSign != Minus) && (xmax>0)) {
361 result += tau*(-exp(-xmax/tau) - -exp(-max(0.,xmin)/tau) ) ; // plus and both
362 }
363 if ((basisSign != Plus) && (xmin<0)) {
364 result -= tau*(-exp(-max(0.,xmin)/tau)) - -tau*exp(-xmax/tau) ; // minus and both
365 }
366
367 return result ;
368 }
369 case sinBasis:
370 {
371 double result(0) ;
372 if (tau==0) return 0 ;
373 double dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
374 if (basisSign != Minus) {
375 double term = exp(-xmax/tau);
376 // We only multiply with the sine term if the coefficient is non zero,
377 // i.e. if xmax was not infinity. Otherwise, we are evaluating the
378 // sine of infinity, whic is NAN! Same applies to the other terms
379 // below.
380 if(term > 0.0) term *= -1/tau*sin(dm*xmax) - dm*cos(dm*xmax);
381 term += dm;
382 result += term;
383 }
384 if (basisSign != Plus) {
385 double term = exp(xmin/tau);
386 if (term > 0.0) term *= -1/tau*sin(dm*(-xmin)) - dm*cos(dm*(-xmin));
387 term += dm;
388 result -= term;
389 }
390 return result / (1/(tau*tau) + dm*dm) ;
391 }
392 case cosBasis:
393 {
394 double result(0) ;
395 if (tau==0) return 1 ;
396 double dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
397 if (basisSign != Minus) {
398 double term = exp(-xmax/tau);
399 if(term > 0.0) term *= -1/tau*cos(dm*xmax) + dm*sin(dm*xmax);
400 term += 1/tau;
401 result += term;
402 }
403 if (basisSign != Plus) {
404 double term = exp( xmin/tau);
405 if(term > 0.0) term *= -1/tau*cos(dm*(-xmin)) + dm*sin(dm*(-xmin));
406 term += 1/tau;
407 result += term;
408 }
409 return result / (1/(tau*tau) + dm*dm) ;
410 }
411 case linBasis:
412 {
413 if (tau==0) return 0 ;
414 double t_max = xmax/tau ;
415 return tau*( 1 - (1 + t_max)*exp(-t_max) ) ;
416 }
417 case quadBasis:
418 {
419 if (tau==0) return 0 ;
420 double t_max = xmax/tau ;
421 return tau*( 2 - (2 + (2 + t_max)*t_max)*exp(-t_max) ) ;
422 }
423 case sinhBasis:
424 {
425 double result(0) ;
426 if (tau==0) return 0 ;
427 double dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
428 double taup = 2*tau/(2-tau*dg);
429 double taum = 2*tau/(2+tau*dg);
430 if (basisSign != Minus) result += 0.5*( taup*(1-exp(-xmax/taup)) - taum*(1-exp(-xmax/taum)) ) ;
431 if (basisSign != Plus) result -= 0.5*( taup*(1-exp( xmin/taup)) - taum*(1-exp( xmin/taum)) ) ;
432 return result ;
433 }
434 case coshBasis:
435 {
436 double result(0) ;
437 if (tau==0) return 1 ;
438 double dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
439 double taup = 2*tau/(2-tau*dg);
440 double taum = 2*tau/(2+tau*dg);
441 if (basisSign != Minus) result += 0.5*( taup*(1-exp(-xmax/taup)) + taum*(1-exp(-xmax/taum)) ) ;
442 if (basisSign != Plus) result += 0.5*( taup*(1-exp( xmin/taup)) + taum*(1-exp( xmin/taum)) ) ;
443 return result ;
444 }
445 default:
446 R__ASSERT(0) ;
447 }
448
449 R__ASSERT(0) ;
450 return 0 ;
451}
452
453
454////////////////////////////////////////////////////////////////////////////////
455
457(const RooAbsAnaConvPdf& convPdf, const RooArgSet &vars, const RooDataSet *prototype,
458 const RooArgSet* auxProto, bool verbose) const
459{
460 RooArgSet forceDirect(convVar()) ;
461 return new RooGenContext(convPdf, vars, prototype, auxProto, verbose, &forceDirect);
462}
463
464
465
466////////////////////////////////////////////////////////////////////////////////
467/// Advertise internal generator for observable x
468
469Int_t RooTruthModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, bool /*staticInitOK*/) const
470{
471 if (matchArgs(directVars,generateVars,x)) return 1 ;
472 return 0 ;
473}
474
475
476
477////////////////////////////////////////////////////////////////////////////////
478/// Implement internal generator for observable x,
479/// x=0 for all events following definition
480/// of delta function
481
483{
484 R__ASSERT(code==1) ;
485 double zero(0.) ;
486 x = zero ;
487 return;
488}
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char name[80]
Definition: TGX11.cxx:110
float xmin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
RooAbsAnaConvPdf is the base class for PDFs that represent a physics model that can be analytically c...
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:77
void removeServer(RooAbsArg &server, bool force=false)
Unregister another RooAbsArg as a server to us, ie, declare that we no longer depend on its value and...
Definition: RooAbsArg.cxx:401
const RefCountList_t & servers() const
List of all servers of this object.
Definition: RooAbsArg.h:204
void addServer(RooAbsArg &server, bool valueProp=true, bool shapeProp=false, std::size_t refCount=1)
Register another RooAbsArg as a server to us, ie, declare that we depend on it.
Definition: RooAbsArg.cxx:350
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:61
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:105
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:56
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:57
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Definition: RooFormulaVar.h:30
RooAbsArg * getParameter(const char *name) const
Return pointer to parameter with given name.
Definition: RooFormulaVar.h:44
Class RooGenContext implement a universal generator context for all RooAbsPdf classes that do not hav...
Definition: RooGenContext.h:30
RooResolutionModel is the base class for PDFs that represent a resolution model that can be convolute...
bool _ownBasis
Flag indicating ownership of _basis.
Int_t _basisCode
Identifier code for selected basis function.
RooAbsRealLValue & convVar() const
Return the convolution variable of the resolution model.
RooFormulaVar * _basis
Basis function convolved with this resolution model.
const RooFormulaVar & basis() const
RooTemplateProxy< RooAbsRealLValue > x
Dependent/convolution variable.
A simple container to hold a batch of data values.
Definition: RooSpan.h:34
double max(const char *rname=nullptr) const
Query upper limit of range. This requires the payload to be RooAbsRealLValue or derived.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RooTruthModel is an implementation of RooResolution model that provides a delta-function resolution m...
Definition: RooTruthModel.h:21
void generateEvent(Int_t code) override
Implement internal generator for observable x, x=0 for all events following definition of delta funct...
double evaluate() const override
Evaluate the truth model: a delta function when used as PDF, the basis function itself,...
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=nullptr) const override
Advertise analytical integrals for compiled basis functions and when used as p.d.f without basis func...
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Advertise internal generator for observable x.
void computeBatch(cudaStream_t *, double *output, size_t size, RooFit::Detail::DataMap const &) const override
Base function for computing multiple values of a RooAbsReal.
~RooTruthModel() override
Destructor.
RooAbsGenContext * modelGenContext(const RooAbsAnaConvPdf &convPdf, const RooArgSet &vars, const RooDataSet *prototype=nullptr, const RooArgSet *auxProto=nullptr, bool verbose=false) const override
double analyticalIntegral(Int_t code, const char *rangeName=nullptr) const override
Implement analytical integrals when used as p.d.f and for compiled basis functions.
Int_t basisCode(const char *name) const override
Return basis code for given basis definition string.
void changeBasis(RooFormulaVar *basis) override
Changes associated bases function to 'inBasis'.
const char * GetTitle() const override
Returns title of object.
Definition: TNamed.h:48
RVec< PromoteType< T > > cosh(const RVec< T > &v)
Definition: RVec.hxx:1808
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1780
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition: RVec.hxx:1800
RVec< PromoteType< T > > sinh(const RVec< T > &v)
Definition: RVec.hxx:1807
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1785
RVec< PromoteType< T > > sin(const RVec< T > &v)
Definition: RVec.hxx:1799
R__EXTERN RooBatchComputeInterface * dispatchCUDA
R__EXTERN RooBatchComputeInterface * dispatchCPU
This dispatch pointer points to an implementation of the compute library, provided one has been loade...
static void output()