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