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 "RooTruthModel.h"
31#include "RooGenContext.h"
32#include "RooAbsAnaConvPdf.h"
33
34#include "TError.h"
35
36#include <algorithm>
37using namespace std ;
38
40;
41
42
43
44////////////////////////////////////////////////////////////////////////////////
45/// Constructor of a truth resolution model, i.e. a delta function in observable 'xIn'
46
47RooTruthModel::RooTruthModel(const char *name, const char *title, RooAbsRealLValue& xIn) :
48 RooResolutionModel(name,title,xIn)
49{
50}
51
52
53
54////////////////////////////////////////////////////////////////////////////////
55/// Copy constructor
56
59{
60}
61
62
63
64////////////////////////////////////////////////////////////////////////////////
65/// Destructor
66
68{
69}
70
71
72
73////////////////////////////////////////////////////////////////////////////////
74/// Return basis code for given basis definition string. Return special
75/// codes for 'known' bases for which compiled definition exists. Return
76/// generic bases code if implementation relies on TFormula interpretation
77/// of basis name
78
80{
81 // Check for optimized basis functions
82 if (!TString("exp(-@0/@1)").CompareTo(name)) return expBasisPlus ;
83 if (!TString("exp(@0/@1)").CompareTo(name)) return expBasisMinus ;
84 if (!TString("exp(-abs(@0)/@1)").CompareTo(name)) return expBasisSum ;
85 if (!TString("exp(-@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisPlus ;
86 if (!TString("exp(@0/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisMinus ;
87 if (!TString("exp(-abs(@0)/@1)*sin(@0*@2)").CompareTo(name)) return sinBasisSum ;
88 if (!TString("exp(-@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisPlus ;
89 if (!TString("exp(@0/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisMinus ;
90 if (!TString("exp(-abs(@0)/@1)*cos(@0*@2)").CompareTo(name)) return cosBasisSum ;
91 if (!TString("(@0/@1)*exp(-@0/@1)").CompareTo(name)) return linBasisPlus ;
92 if (!TString("(@0/@1)*(@0/@1)*exp(-@0/@1)").CompareTo(name)) return quadBasisPlus ;
93 if (!TString("exp(-@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisPlus;
94 if (!TString("exp(@0/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisMinus;
95 if (!TString("exp(-abs(@0)/@1)*cosh(@0*@2/2)").CompareTo(name)) return coshBasisSum;
96 if (!TString("exp(-@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisPlus;
97 if (!TString("exp(@0/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisMinus;
98 if (!TString("exp(-abs(@0)/@1)*sinh(@0*@2/2)").CompareTo(name)) return sinhBasisSum;
99
100 // Truth model is delta function, i.e. convolution integral
101 // is basis function, therefore we can handle any basis function
102 return genericBasis ;
103}
104
105
106
107
108////////////////////////////////////////////////////////////////////////////////
109/// Changes associated bases function to 'inBasis'
110
112{
113 // Process change basis function. Since we actually
114 // evaluate the basis function object, we need to
115 // adjust our client-server links to the basis function here
116
117 // Remove client-server link to old basis
118 if (_basis) {
120 }
121
122 // Change basis pointer and update client-server link
123 _basis = inBasis ;
124 if (_basis) {
126 }
127
128 _basisCode = inBasis?basisCode(inBasis->GetTitle()):0 ;
129}
130
131
132
133
134
135////////////////////////////////////////////////////////////////////////////////
136/// Evaluate the truth model: a delta function when used as PDF,
137/// the basis function itself, when convoluted with a basis function.
138
140{
141 // No basis: delta function
142 if (_basisCode == noBasis) {
143 if (x==0) return 1 ;
144 return 0 ;
145 }
146
147 // Generic basis: evaluate basis function object
148 if (_basisCode == genericBasis) {
149 return basis().getVal() ;
150 }
151
152 // Precompiled basis functions
153 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
154 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
155
156 // Enforce sign compatibility
157 if ((basisSign==Minus && x>0) ||
158 (basisSign==Plus && x<0)) return 0 ;
159
160
161 Double_t tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ;
162 // Return desired basis function
163 switch(basisType) {
164 case expBasis: {
165 //cout << " RooTruthModel::eval(" << GetName() << ") expBasis mode ret = " << exp(-fabs((Double_t)x)/tau) << " tau = " << tau << endl ;
166 return exp(-fabs((Double_t)x)/tau) ;
167 }
168 case sinBasis: {
169 Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
170 return exp(-fabs((Double_t)x)/tau)*sin(x*dm) ;
171 }
172 case cosBasis: {
173 Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
174 return exp(-fabs((Double_t)x)/tau)*cos(x*dm) ;
175 }
176 case linBasis: {
177 Double_t tscaled = fabs((Double_t)x)/tau;
178 return exp(-tscaled)*tscaled ;
179 }
180 case quadBasis: {
181 Double_t tscaled = fabs((Double_t)x)/tau;
182 return exp(-tscaled)*tscaled*tscaled;
183 }
184 case sinhBasis: {
185 Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
186 return exp(-fabs((Double_t)x)/tau)*sinh(x*dg/2) ;
187 }
188 case coshBasis: {
189 Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
190 return exp(-fabs((Double_t)x)/tau)*cosh(x*dg/2) ;
191 }
192 default:
193 R__ASSERT(0) ;
194 }
195
196 return 0 ;
197}
198
199
200
201////////////////////////////////////////////////////////////////////////////////
202/// Advertise analytical integrals for compiled basis functions and when used
203/// as p.d.f without basis function.
204
205Int_t RooTruthModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
206{
207 switch(_basisCode) {
208
209 // Analytical integration capability of raw PDF
210 case noBasis:
211 if (matchArgs(allVars,analVars,convVar())) return 1 ;
212 break ;
213
214 // Analytical integration capability of convoluted PDF
215 case expBasisPlus:
216 case expBasisMinus:
217 case expBasisSum:
218 case sinBasisPlus:
219 case sinBasisMinus:
220 case sinBasisSum:
221 case cosBasisPlus:
222 case cosBasisMinus:
223 case cosBasisSum:
224 case linBasisPlus:
225 case quadBasisPlus:
226 case sinhBasisPlus:
227 case sinhBasisMinus:
228 case sinhBasisSum:
229 case coshBasisPlus:
230 case coshBasisMinus:
231 case coshBasisSum:
232 if (matchArgs(allVars,analVars,convVar())) return 1 ;
233 break ;
234 }
235
236 return 0 ;
237}
238
239
240
241////////////////////////////////////////////////////////////////////////////////
242/// Implement analytical integrals when used as p.d.f and for compiled
243/// basis functions.
244
245Double_t RooTruthModel::analyticalIntegral(Int_t code, const char* rangeName) const
246{
247
248 // Code must be 1
249 R__ASSERT(code==1) ;
250
251 // Unconvoluted PDF
252 if (_basisCode==noBasis) return 1 ;
253
254 // Precompiled basis functions
255 BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 );
256 BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ;
257 //cout << " calling RooTruthModel::analyticalIntegral with basisType " << basisType << endl;
258
259 Double_t tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ;
260 switch (basisType) {
261 case expBasis:
262 {
263 // WVE fixed for ranges
264 Double_t result(0) ;
265 if (tau==0) return 1 ;
266 if ((basisSign != Minus) && (x.max(rangeName)>0)) {
267 result += tau*(-exp(-x.max(rangeName)/tau) - -exp(-max(0.,x.min(rangeName))/tau) ) ; // plus and both
268 }
269 if ((basisSign != Plus) && (x.min(rangeName)<0)) {
270 result -= tau*(-exp(-max(0.,x.min(rangeName))/tau)) - -tau*exp(-x.max(rangeName)/tau) ; // minus and both
271 }
272
273 return result ;
274 }
275 case sinBasis:
276 {
277 Double_t result(0) ;
278 if (tau==0) return 0 ;
279 Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
280 if (basisSign != Minus) result += exp(-x.max(rangeName)/tau)*(-1/tau*sin(dm*x.max(rangeName)) - dm*cos(dm*x.max(rangeName))) + dm; // fixed FMV 08/29/03
281 if (basisSign != Plus) result -= exp( x.min(rangeName)/tau)*(-1/tau*sin(dm*(-x.min(rangeName))) - dm*cos(dm*(-x.min(rangeName)))) + dm ; // fixed FMV 08/29/03
282 return result / (1/(tau*tau) + dm*dm) ;
283 }
284 case cosBasis:
285 {
286 Double_t result(0) ;
287 if (tau==0) return 1 ;
288 Double_t dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
289 if (basisSign != Minus) result += exp(-x.max(rangeName)/tau)*(-1/tau*cos(dm*x.max(rangeName)) + dm*sin(dm*x.max(rangeName))) + 1/tau ;
290 if (basisSign != Plus) result += exp( x.min(rangeName)/tau)*(-1/tau*cos(dm*(-x.min(rangeName))) + dm*sin(dm*(-x.min(rangeName)))) + 1/tau ; // fixed FMV 08/29/03
291 return result / (1/(tau*tau) + dm*dm) ;
292 }
293 case linBasis:
294 {
295 if (tau==0) return 0 ;
296 Double_t t_max = x.max(rangeName)/tau ;
297 return tau*( 1 - (1 + t_max)*exp(-t_max) ) ;
298 }
299 case quadBasis:
300 {
301 if (tau==0) return 0 ;
302 Double_t t_max = x.max(rangeName)/tau ;
303 return tau*( 2 - (2 + (2 + t_max)*t_max)*exp(-t_max) ) ;
304 }
305 case sinhBasis:
306 {
307 Double_t result(0) ;
308 if (tau==0) return 0 ;
309 Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
310 Double_t taup = 2*tau/(2-tau*dg);
311 Double_t taum = 2*tau/(2+tau*dg);
312 if (basisSign != Minus) result += 0.5*( taup*(1-exp(-x.max(rangeName)/taup)) - taum*(1-exp(-x.max(rangeName)/taum)) ) ;
313 if (basisSign != Plus) result -= 0.5*( taup*(1-exp( x.min(rangeName)/taup)) - taum*(1-exp( x.min(rangeName)/taum)) ) ;
314 return result ;
315 }
316 case coshBasis:
317 {
318 Double_t result(0) ;
319 if (tau==0) return 1 ;
320 Double_t dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ;
321 Double_t taup = 2*tau/(2-tau*dg);
322 Double_t taum = 2*tau/(2+tau*dg);
323 if (basisSign != Minus) result += 0.5*( taup*(1-exp(-x.max(rangeName)/taup)) + taum*(1-exp(-x.max(rangeName)/taum)) ) ;
324 if (basisSign != Plus) result += 0.5*( taup*(1-exp( x.min(rangeName)/taup)) + taum*(1-exp( x.min(rangeName)/taum)) ) ;
325 return result ;
326 }
327 default:
328 R__ASSERT(0) ;
329 }
330
331 R__ASSERT(0) ;
332 return 0 ;
333}
334
335
336////////////////////////////////////////////////////////////////////////////////
337
339(const RooAbsAnaConvPdf& convPdf, const RooArgSet &vars, const RooDataSet *prototype,
340 const RooArgSet* auxProto, Bool_t verbose) const
341{
342 RooArgSet forceDirect(convVar()) ;
343 return new RooGenContext(dynamic_cast<const RooAbsPdf&>(convPdf), vars, prototype,
344 auxProto, verbose, &forceDirect) ;
345}
346
347
348
349////////////////////////////////////////////////////////////////////////////////
350/// Advertise internal generator for observable x
351
352Int_t RooTruthModel::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
353{
354 if (matchArgs(directVars,generateVars,x)) return 1 ;
355 return 0 ;
356}
357
358
359
360////////////////////////////////////////////////////////////////////////////////
361/// Implement internal generator for observable x,
362/// x=0 for all events following definition
363/// of delta function
364
366{
367 R__ASSERT(code==1) ;
368 Double_t zero(0.) ;
369 x = zero ;
370 return;
371}
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#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
RooAbsAnaConvPdf is the base class for PDFs that represent a physics model that can be analytically c...
void addServer(RooAbsArg &server, Bool_t valueProp=kTRUE, Bool_t shapeProp=kFALSE, std::size_t refCount=1)
Register another RooAbsArg as a server to us, ie, declare that we depend on it.
Definition: RooAbsArg.cxx:379
void removeServer(RooAbsArg &server, Bool_t force=kFALSE)
Unregister another RooAbsArg as a server to us, ie, declare that we no longer depend on its value and...
Definition: RooAbsArg.cxx:430
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:63
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:93
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:55
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...
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 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.
RooTruthModel is an implementation of RooResolution model that provides a delta-function resolution m...
Definition: RooTruthModel.h:21
Double_t evaluate() const override
Evaluate the truth model: a delta function when used as PDF, the basis function itself,...
void generateEvent(Int_t code) override
Implement internal generator for observable x, x=0 for all events following definition of delta funct...
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) 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_t staticInitOK=kTRUE) const override
Advertise internal generator for observable x.
RooAbsGenContext * modelGenContext(const RooAbsAnaConvPdf &convPdf, const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE) const override
~RooTruthModel() override
Destructor.
Int_t basisCode(const char *name) const override
Return basis code for given basis definition string.
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const override
Implement analytical integrals when used as p.d.f and for compiled basis functions.
void changeBasis(RooFormulaVar *basis) override
Changes associated bases function to 'inBasis'.
const char * GetTitle() const override
Returns title of object.
Definition: TNamed.h:48
Basic string class.
Definition: TString.h:136
RVec< PromoteType< T > > cosh(const RVec< T > &v)
Definition: RVec.hxx:1767
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition: RVec.hxx:1759
RVec< PromoteType< T > > sinh(const RVec< T > &v)
Definition: RVec.hxx:1766
RVec< PromoteType< T > > exp(const RVec< T > &v)
Definition: RVec.hxx:1744
RVec< PromoteType< T > > sin(const RVec< T > &v)
Definition: RVec.hxx:1758
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)