Logo ROOT   6.14/05
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 
22 RooTruthModel is an implementation of RooResolution
23 model that provides a delta-function resolution model
24 The truth model supports <i>all</i> basis functions because it evaluates each basis function as
25 as a RooFormulaVar. The 6 basis functions used in B mixing and decay and 2 basis
26 functions 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>
39 using namespace std ;
40 
42 ;
43 
44 
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 /// Constructor of a truth resolution model, i.e. a delta function in observable 'xIn'
48 
49 RooTruthModel::RooTruthModel(const char *name, const char *title, RooRealVar& xIn) :
50  RooResolutionModel(name,title,xIn)
51 {
52 }
53 
54 
55 
56 ////////////////////////////////////////////////////////////////////////////////
57 /// Copy constructor
58 
59 RooTruthModel::RooTruthModel(const RooTruthModel& other, const char* name) :
60  RooResolutionModel(other,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 
81 Int_t RooTruthModel::basisCode(const char* name) const
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) {
121  removeServer(*_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 
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 
207 Int_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 
247 Double_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 
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 
354 Int_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 }
RooTruthModel is an implementation of RooResolution model that provides a delta-function resolution m...
Definition: RooTruthModel.h:21
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
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...
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.
#define R__ASSERT(e)
Definition: TError.h:96
RooRealVar & convVar() const
Return the convolution variable of the resolution model.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
STL namespace.
double cos(double)
const RooFormulaVar & basis() const
virtual ~RooTruthModel()
Destructor.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t tau
Definition: TRolke.cxx:630
double sinh(double)
virtual RooAbsGenContext * modelGenContext(const RooAbsAnaConvPdf &convPdf, const RooArgSet &vars, const RooDataSet *prototype=0, const RooArgSet *auxProto=0, Bool_t verbose=kFALSE) const
RooAbsGenContext is the abstract base class for generator contexts of RooAbsPdf objects.
double sin(double)
RooFormulaVar * _basis
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
void generateEvent(Int_t code)
Implement internal generator for observable x, x=0 for all events following definition of delta funct...
virtual void changeBasis(RooFormulaVar *basis)
Changes associated bases function to &#39;inBasis&#39;.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double cosh(double)
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:29
Class RooGenContext implement a universal generator context for all RooAbsPdf classes that do not hav...
Definition: RooGenContext.h:30
const Bool_t kFALSE
Definition: RtypesCore.h:88
#define ClassImp(name)
Definition: Rtypes.h:359
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
virtual Double_t evaluate() const
Evaluate the truth model: a delta function when used as PDF, the basis function itself, when convoluted with a basis function.
void addServer(RooAbsArg &server, Bool_t valueProp=kTRUE, Bool_t shapeProp=kFALSE)
Register another RooAbsArg as a server to us, ie, declare that we depend on it.
Definition: RooAbsArg.cxx:334
double exp(double)
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:385
const Bool_t kTRUE
Definition: RtypesCore.h:87
RooAbsArg * getParameter(const char *name) const
Definition: RooFormulaVar.h:39
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Advertise internal generator for observable x.
char name[80]
Definition: TGX11.cxx:109
virtual Int_t basisCode(const char *name) const
Return basis code for given basis definition string.
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48