Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
22Implements a RooResolution model that corresponds to a delta function.
23The truth model supports <i>all</i> basis functions because it evaluates each basis function as
24as a RooFormulaVar. The 6 basis functions used in B mixing and decay and 2 basis
25functions used in D mixing have been hand coded for increased execution speed.
26**/
27
28#include <RooTruthModel.h>
29
30#include <RooAbsAnaConvPdf.h>
31#include <RooBatchCompute.h>
32#include <RooGenContext.h>
33
35
36#include <Riostream.h>
37
38#include <TError.h>
39
40#include <algorithm>
41#include <array>
42#include <cmath>
43#include <limits>
44
45namespace {
46
47enum RooTruthBasis {
48 noBasis = 0,
49 expBasisMinus = 1,
50 expBasisSum = 2,
51 expBasisPlus = 3,
52 sinBasisMinus = 11,
53 sinBasisSum = 12,
54 sinBasisPlus = 13,
55 cosBasisMinus = 21,
56 cosBasisSum = 22,
57 cosBasisPlus = 23,
58 linBasisPlus = 33,
59 quadBasisPlus = 43,
60 coshBasisMinus = 51,
61 coshBasisSum = 52,
62 coshBasisPlus = 53,
63 sinhBasisMinus = 61,
64 sinhBasisSum = 62,
65 sinhBasisPlus = 63,
66 genericBasis = 100
67};
68
69enum BasisType {
70 none = 0,
71 expBasis = 1,
72 sinBasis = 2,
73 cosBasis = 3,
74 linBasis = 4,
75 quadBasis = 5,
76 coshBasis = 6,
77 sinhBasis = 7
78};
79
80enum BasisSign { Both = 0, Plus = +1, Minus = -1 };
81
82} // namespace
83
84
85////////////////////////////////////////////////////////////////////////////////
86/// Constructor of a truth resolution model, i.e. a delta function in observable 'xIn'
87
88RooTruthModel::RooTruthModel(const char *name, const char *title, RooAbsRealLValue& xIn) :
90{
91}
92
93////////////////////////////////////////////////////////////////////////////////
94/// Return basis code for given basis definition string. Return special
95/// codes for 'known' bases for which compiled definition exists. Return
96/// generic bases code if implementation relies on TFormula interpretation
97/// of basis name
98
100{
101 std::string str = name;
102
103 // Remove whitespaces from the input string
104 str.erase(remove(str.begin(),str.end(),' '),str.end());
105
106 // Check for optimized basis functions
107 if (str == "exp(-@0/@1)") return expBasisPlus ;
108 if (str == "exp(@0/@1)") return expBasisMinus ;
109 if (str == "exp(-abs(@0)/@1)") return expBasisSum ;
110 if (str == "exp(-@0/@1)*sin(@0*@2)") return sinBasisPlus ;
111 if (str == "exp(@0/@1)*sin(@0*@2)") return sinBasisMinus ;
112 if (str == "exp(-abs(@0)/@1)*sin(@0*@2)") return sinBasisSum ;
113 if (str == "exp(-@0/@1)*cos(@0*@2)") return cosBasisPlus ;
114 if (str == "exp(@0/@1)*cos(@0*@2)") return cosBasisMinus ;
115 if (str == "exp(-abs(@0)/@1)*cos(@0*@2)") return cosBasisSum ;
116 if (str == "(@0/@1)*exp(-@0/@1)") return linBasisPlus ;
117 if (str == "(@0/@1)*(@0/@1)*exp(-@0/@1)") return quadBasisPlus ;
118 if (str == "exp(-@0/@1)*cosh(@0*@2/2)") return coshBasisPlus;
119 if (str == "exp(@0/@1)*cosh(@0*@2/2)") return coshBasisMinus;
120 if (str == "exp(-abs(@0)/@1)*cosh(@0*@2/2)") return coshBasisSum;
121 if (str == "exp(-@0/@1)*sinh(@0*@2/2)") return sinhBasisPlus;
122 if (str == "exp(@0/@1)*sinh(@0*@2/2)") return sinhBasisMinus;
123 if (str == "exp(-abs(@0)/@1)*sinh(@0*@2/2)") return sinhBasisSum;
124
125 // Truth model is delta function, i.e. convolution integral is basis
126 // function, therefore we can handle any basis function
127 return genericBasis ;
128}
129
130
131
132////////////////////////////////////////////////////////////////////////////////
133/// Changes associated bases function to 'inBasis'
134
136{
137 // Remove client-server link to old basis
138 if (_basis) {
139 if (_basisCode == genericBasis) {
140 // In the case of a generic basis, we evaluate it directly, so the
141 // basis was a direct server.
143 } else {
144 for (RooAbsArg *basisServer : _basis->servers()) {
146 }
147 }
148
149 if (_ownBasis) {
150 delete _basis;
151 }
152 }
153 _ownBasis = false;
154
155 _basisCode = inBasis ? basisCode(inBasis->GetTitle()) : 0;
156
157 // Change basis pointer and update client-server link
158 _basis = inBasis;
159 if (_basis) {
160 if (_basisCode == genericBasis) {
161 // Since we actually evaluate the basis function object, we need to
162 // adjust our client-server links to the basis function here
163 addServer(*_basis, true, false);
164 } else {
165 for (RooAbsArg *basisServer : _basis->servers()) {
166 addServer(*basisServer, true, false);
167 }
168 }
169 }
170}
171
172
173
174////////////////////////////////////////////////////////////////////////////////
175/// Evaluate the truth model: a delta function when used as PDF,
176/// the basis function itself, when convoluted with a basis function.
177
179{
180 // No basis: delta function
181 if (_basisCode == noBasis) {
182 if (x==0) return 1 ;
183 return 0 ;
184 }
185
186 // Generic basis: evaluate basis function object
187 if (_basisCode == genericBasis) {
188 return basis().getVal() ;
189 }
190
191 // Precompiled basis functions
192 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
193 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
194
195 // Enforce sign compatibility
196 if ((basisSign==Minus && x>0) ||
197 (basisSign==Plus && x<0)) return 0 ;
198
199
200 double tau = (static_cast<RooAbsReal*>(basis().getParameter(1)))->getVal() ;
201 // Return desired basis function
202 switch(basisType) {
203 case expBasis: {
204 return std::exp(-std::abs((double)x)/tau) ;
205 }
206 case sinBasis: {
207 double dm = (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() ;
208 return std::exp(-std::abs((double)x)/tau)*std::sin(x*dm) ;
209 }
210 case cosBasis: {
211 double dm = (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() ;
212 return std::exp(-std::abs((double)x)/tau)*std::cos(x*dm) ;
213 }
214 case linBasis: {
215 double tscaled = std::abs((double)x)/tau;
216 return std::exp(-tscaled)*tscaled ;
217 }
218 case quadBasis: {
219 double tscaled = std::abs((double)x)/tau;
220 return std::exp(-tscaled)*tscaled*tscaled;
221 }
222 case sinhBasis: {
223 double dg = (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() ;
224 return std::exp(-std::abs((double)x)/tau)*std::sinh(x*dg/2) ;
225 }
226 case coshBasis: {
227 double dg = (static_cast<RooAbsReal*>(basis().getParameter(2)))->getVal() ;
228 return std::exp(-std::abs((double)x)/tau)*std::cosh(x*dg/2) ;
229 }
230 default:
231 R__ASSERT(0) ;
232 }
233
234 return 0 ;
235}
236
237
239{
240 auto config = ctx.config(this);
241 auto xVals = ctx.at(x);
242
243 // No basis: delta function
244 if (_basisCode == noBasis) {
246 return;
247 }
248
249 // Generic basis: evaluate basis function object
250 if (_basisCode == genericBasis) {
251 RooBatchCompute::compute(config, RooBatchCompute::Identity, ctx.output(), {ctx.at(&basis())});
252 return;
253 }
254
255 // Precompiled basis functions
256 const BasisType basisType = static_cast<BasisType>((_basisCode == 0) ? 0 : (_basisCode / 10) + 1);
257
258 // Cast the int from the enum to double because we can only pass doubles to
259 // RooBatchCompute at this point.
260 const double basisSign = static_cast<double>((BasisSign)(_basisCode - 10 * (basisType - 1) - 2));
261
262 auto param1 = static_cast<RooAbsReal const *>(basis().getParameter(1));
263 auto param2 = static_cast<RooAbsReal const *>(basis().getParameter(2));
264 auto param1Vals = param1 ? ctx.at(param1) : std::span<const double>{};
265 auto param2Vals = param2 ? ctx.at(param2) : std::span<const double>{};
266
267 // Return desired basis function
268 std::array<double, 1> extraArgs{basisSign};
269 switch (basisType) {
270 case expBasis: {
271 RooBatchCompute::compute(config, RooBatchCompute::TruthModelExpBasis, ctx.output(), {xVals, param1Vals},
272 extraArgs);
273 break;
274 }
275 case sinBasis: {
277 {xVals, param1Vals, param2Vals}, extraArgs);
278 break;
279 }
280 case cosBasis: {
282 {xVals, param1Vals, param2Vals}, extraArgs);
283 break;
284 }
285 case linBasis: {
286 RooBatchCompute::compute(config, RooBatchCompute::TruthModelLinBasis, ctx.output(), {xVals, param1Vals},
287 extraArgs);
288 break;
289 }
290 case quadBasis: {
291 RooBatchCompute::compute(config, RooBatchCompute::TruthModelQuadBasis, ctx.output(), {xVals, param1Vals},
292 extraArgs);
293 break;
294 }
295 case sinhBasis: {
297 {xVals, param1Vals, param2Vals}, extraArgs);
298 break;
299 }
300 case coshBasis: {
302 {xVals, param1Vals, param2Vals}, extraArgs);
303 break;
304 }
305 default: R__ASSERT(0);
306 }
307}
308
309
310////////////////////////////////////////////////////////////////////////////////
311/// Advertise analytical integrals for compiled basis functions and when used
312/// as p.d.f without basis function.
313
314Int_t RooTruthModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
315{
316 switch(_basisCode) {
317
318 // Analytical integration capability of raw PDF
319 case noBasis:
320 if (matchArgs(allVars,analVars,convVar())) return 1 ;
321 break ;
322
323 // Analytical integration capability of convoluted PDF
324 case expBasisPlus:
325 case expBasisMinus:
326 case expBasisSum:
327 case sinBasisPlus:
328 case sinBasisMinus:
329 case sinBasisSum:
330 case cosBasisPlus:
331 case cosBasisMinus:
332 case cosBasisSum:
333 case linBasisPlus:
334 case quadBasisPlus:
335 case sinhBasisPlus:
336 case sinhBasisMinus:
337 case sinhBasisSum:
338 case coshBasisPlus:
339 case coshBasisMinus:
340 case coshBasisSum:
341 if (matchArgs(allVars,analVars,convVar())) return 1 ;
342 break ;
343 }
344
345 return 0 ;
346}
347
348
349namespace {
350
351// From asking WolframAlpha: integrate exp(-x/tau) over x.
352inline double indefiniteIntegralExpBasisPlus(double x, double tau, double /*dm*/)
353{
354 // Restrict to positive x
355 x = std::max(x, 0.0);
356 return -tau * std::exp(-x / tau);
357}
358
359// From asking WolframAlpha: integrate exp(-x/tau)* x / tau over x.
360inline double indefiniteIntegralLinBasisPlus(double x, double tau, double /*dm*/)
361{
362 // Restrict to positive x
363 x = std::max(x, 0.0);
364 return -(tau + x) * std::exp(-x / tau);
365}
366
367// From asking WolframAlpha: integrate exp(-x/tau) * (x / tau)^2 over x.
368inline double indefiniteIntegralQuadBasisPlus(double x, double tau, double /*dm*/)
369{
370 // Restrict to positive x
371 x = std::max(x, 0.0);
372 return -(std::exp(-x / tau) * (2 * tau * tau + x * x + 2 * tau * x)) / tau;
373}
374
375// A common factor that appears in the integrals of the trigonometric
376// function bases (sin and cos).
377inline double commonFactorPlus(double x, double tau, double dm)
378{
379 const double num = tau * std::exp(-x / tau);
380 const double den = dm * dm * tau * tau + 1.0;
381 return num / den;
382}
383
384// A common factor that appears in the integrals of the hyperbolic
385// trigonometric function bases (sinh and cosh).
386inline double commonFactorHyperbolicPlus(double x, double tau, double dm)
387{
388 const double num = 2 * tau * std::exp(-x / tau);
389 const double den = dm * dm * tau * tau - 4.0;
390 return num / den;
391}
392
393// From asking WolframAlpha: integrate exp(-x/tau)*sin(x*m) over x.
394inline double indefiniteIntegralSinBasisPlus(double x, double tau, double dm)
395{
396 // Restrict to positive x
397 x = std::max(x, 0.0);
398 const double fac = commonFactorPlus(x, tau, dm);
399 // Only multiply with the sine term if the coefficient is non zero,
400 // i.e. if x was not infinity. Otherwise, we are evaluating the
401 // sine of infinity, which is NAN!
402 return fac != 0.0 ? fac * (-tau * dm * std::cos(dm * x) - std::sin(dm * x)) : 0.0;
403}
404
405// From asking WolframAlpha: integrate exp(-x/tau)*cos(x*m) over x.
406inline double indefiniteIntegralCosBasisPlus(double x, double tau, double dm)
407{
408 // Restrict to positive x
409 x = std::max(x, 0.0);
410 const double fac = commonFactorPlus(x, tau, dm);
411 return fac != 0.0 ? fac * (tau * dm * std::sin(dm * x) - std::cos(dm * x)) : 0.0;
412}
413
414// From asking WolframAlpha: integrate exp(-x/tau)*sinh(x*m/2) over x.
415inline double indefiniteIntegralSinhBasisPlus(double x, double tau, double dm)
416{
417 // Restrict to positive x
418 x = std::max(x, 0.0);
419 const double fac = commonFactorHyperbolicPlus(x, tau, dm);
420 const double arg = 0.5 * dm * x;
421 return fac != 0.0 ? fac * (tau * dm * std::cosh(arg) - 2. * std::sinh(arg)) : 0.0;
422}
423
424// From asking WolframAlpha: integrate exp(-x/tau)*cosh(x*m/2) over x.
425inline double indefiniteIntegralCoshBasisPlus(double x, double tau, double dm)
426{
427 // Restrict to positive x
428 x = std::max(x, 0.0);
429 const double fac = commonFactorHyperbolicPlus(x, tau, dm);
430 const double arg = 0.5 * dm * x;
431 return fac != 0.0 ? fac * (tau * dm * std::sinh(arg) + 2. * std::cosh(arg)) : 0.0;
432}
433
434// Integrate one of the basis functions. Takes a function that represents the
435// indefinite integral, some parameters, and a flag that indicates whether the
436// basis function is symmetric or antisymmetric. This information is used to
437// evaluate the integrals for the "Minus" and "Sum" cases.
438template <class Function>
439double definiteIntegral(Function indefiniteIntegral, double xmin, double xmax, double tau, double dm,
440 BasisSign basisSign, bool isSymmetric)
441{
442 // Note: isSymmetric == false implies antisymmetric
443 if (tau == 0.0)
444 return isSymmetric ? 1.0 : 0.0;
445 double result = 0.0;
446 if (basisSign != Minus) {
447 result += indefiniteIntegral(xmax, tau, dm) - indefiniteIntegral(xmin, tau, dm);
448 }
449 if (basisSign != Plus) {
450 const double resultMinus = indefiniteIntegral(-xmax, tau, dm) - indefiniteIntegral(-xmin, tau, dm);
452 }
453 return result;
454}
455
456} // namespace
457
458////////////////////////////////////////////////////////////////////////////////
459/// Implement analytical integrals when used as p.d.f and for compiled
460/// basis functions.
461
462double RooTruthModel::analyticalIntegral(Int_t code, const char *rangeName) const
463{
464 // Code must be 1
465 R__ASSERT(code == 1);
466
467 // Unconvoluted PDF
468 if (_basisCode == noBasis)
469 return 1;
470
471 // Precompiled basis functions
472 BasisType basisType = (BasisType)((_basisCode == 0) ? 0 : (_basisCode / 10) + 1);
473 BasisSign basisSign = (BasisSign)(_basisCode - 10 * (basisType - 1) - 2);
474
475 const bool needsDm =
477
478 const double tau = (static_cast<RooAbsReal *>(basis().getParameter(1)))->getVal();
479 const double dm =
480 needsDm ? (static_cast<RooAbsReal *>(basis().getParameter(2)))->getVal() : std::numeric_limits<Double_t>::quiet_NaN();
481
482 const double xmin = x.min(rangeName);
483 const double xmax = x.max(rangeName);
484
485 auto integrate = [&](auto indefiniteIntegral, bool isSymmetric) {
487 };
488
489 switch (basisType) {
490 case expBasis: return integrate(indefiniteIntegralExpBasisPlus, /*isSymmetric=*/true);
491 case sinBasis: return integrate(indefiniteIntegralSinBasisPlus, /*isSymmetric=*/false);
492 case cosBasis: return integrate(indefiniteIntegralCosBasisPlus, /*isSymmetric=*/true);
493 case linBasis: return integrate(indefiniteIntegralLinBasisPlus, /*isSymmetric=*/false);
494 case quadBasis: return integrate(indefiniteIntegralQuadBasisPlus, /*isSymmetric=*/true);
495 case sinhBasis: return integrate(indefiniteIntegralSinhBasisPlus, /*isSymmetric=*/false);
496 case coshBasis: return integrate(indefiniteIntegralCoshBasisPlus, /*isSymmetric=*/true);
497 default: R__ASSERT(0);
498 }
499
500 R__ASSERT(0);
501 return 0;
502}
503
504
505////////////////////////////////////////////////////////////////////////////////
506
508(const RooAbsAnaConvPdf& convPdf, const RooArgSet &vars, const RooDataSet *prototype,
509 const RooArgSet* auxProto, bool verbose) const
510{
512 return new RooGenContext(convPdf, vars, prototype, auxProto, verbose, &forceDirect);
513}
514
515
516
517////////////////////////////////////////////////////////////////////////////////
518/// Advertise internal generator for observable x
519
521{
522 if (matchArgs(directVars,generateVars,x)) return 1 ;
523 return 0 ;
524}
525
526
527
528////////////////////////////////////////////////////////////////////////////////
529/// Implement internal generator for observable x,
530/// x=0 for all events following definition
531/// of delta function
532
534{
535 R__ASSERT(code==1) ;
536 double zero(0.) ;
537 x = zero ;
538 return;
539}
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
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
float xmax
Double_t(* Function)(Double_t)
Definition Functor.C:4
Base class for PDFs that represent a physics model that can be analytically convolved with a resoluti...
Common abstract base class for objects that represent a value and a "shape" in RooFit.
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...
const RefCountList_t & servers() const
List of all servers of this object.
Definition RooAbsArg.h:149
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.
Abstract base class for generator contexts of RooAbsPdf objects.
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:63
bool matchArgs(const RooArgSet &allDeps, RooArgSet &analDeps, const RooArgProxy &a, const Proxies &... proxies) const
Definition RooAbsReal.h:428
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
Container class to hold unbinned data.
Definition RooDataSet.h:32
std::span< const double > at(RooAbsArg const *arg, RooAbsArg const *caller=nullptr)
std::span< double > output()
RooBatchCompute::Config config(RooAbsArg const *arg) const
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Implements a universal generator context for all RooAbsPdf classes that do not have or need a special...
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.
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.
void doEval(RooFit::EvalContext &) const override
Base function for computing multiple values of a RooAbsReal.
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.
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.
RooTruthModel()=default
void changeBasis(RooFormulaVar *basis) override
Changes associated bases function to 'inBasis'.
Double_t x[n]
Definition legend1.C:17
void compute(Config cfg, Computer comp, std::span< double > output, VarSpan vars, ArgSpan extraArgs={})