Logo ROOT  
Reference Guide
RooBDecay.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * PL, Parker C Lund, UC Irvine *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
9 * *
10 * Copyright (c) 2000-2005, Regents of the University of California *
11 * and Stanford University. All rights reserved. *
12 * *
13 * Redistribution and use in source and binary forms, *
14 * with or without modification, are permitted according to the terms *
15 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
16 *****************************************************************************/
17
18
19/** \class RooBDecay
20 \ingroup Roofit
21
22Most general description of B decay time distribution with effects
23of CP violation, mixing and life time differences. This function can
24be analytically convolved with any RooResolutionModel implementation
25**/
26
27#include "RooFit.h"
28
29#include "Riostream.h"
30#include "TMath.h"
31#include "RooBDecay.h"
32#include "RooRealVar.h"
33#include "RooRandom.h"
34
35#include "TError.h"
36
37using namespace std;
38
40
41////////////////////////////////////////////////////////////////////////////////
42
43RooBDecay::RooBDecay(const char *name, const char* title,
44 RooRealVar& t, RooAbsReal& tau, RooAbsReal& dgamma,
46 RooAbsReal& dm, const RooResolutionModel& model, DecayType type) :
47 RooAbsAnaConvPdf(name, title, model, t),
48 _t("t", "time", this, t),
49 _tau("tau", "Average Decay Time", this, tau),
50 _dgamma("dgamma", "Delta Gamma", this, dgamma),
51 _f0("f0", "Cosh Coefficient", this, f0),
52 _f1("f1", "Sinh Coefficient", this, f1),
53 _f2("f2", "Cos Coefficient", this, f2),
54 _f3("f3", "Sin Coefficient", this, f3),
55 _dm("dm", "Delta Mass", this, dm),
56 _type(type)
57
58{
59 //Constructor
60 switch(type)
61 {
62 case SingleSided:
63 _basisCosh = declareBasis("exp(-@0/@1)*cosh(@0*@2/2)", RooArgList(tau,dgamma));
64 _basisSinh = declareBasis("exp(-@0/@1)*sinh(@0*@2/2)", RooArgList(tau,dgamma));
65 _basisCos = declareBasis("exp(-@0/@1)*cos(@0*@2)",RooArgList(tau, dm));
66 _basisSin = declareBasis("exp(-@0/@1)*sin(@0*@2)",RooArgList(tau, dm));
67 break;
68 case Flipped:
69 _basisCosh = declareBasis("exp(@0/@1)*cosh(@0*@2/2)", RooArgList(tau,dgamma));
70 _basisSinh = declareBasis("exp(@0/@1)*sinh(@0*@2/2)", RooArgList(tau,dgamma));
71 _basisCos = declareBasis("exp(@0/@1)*cos(@0*@2)",RooArgList(tau, dm));
72 _basisSin = declareBasis("exp(@0/@1)*sin(@0*@2)",RooArgList(tau, dm));
73 break;
74 case DoubleSided:
75 _basisCosh = declareBasis("exp(-abs(@0)/@1)*cosh(@0*@2/2)", RooArgList(tau,dgamma));
76 _basisSinh = declareBasis("exp(-abs(@0)/@1)*sinh(@0*@2/2)", RooArgList(tau,dgamma));
77 _basisCos = declareBasis("exp(-abs(@0)/@1)*cos(@0*@2)",RooArgList(tau, dm));
78 _basisSin = declareBasis("exp(-abs(@0)/@1)*sin(@0*@2)",RooArgList(tau, dm));
79 break;
80 }
81}
82
83////////////////////////////////////////////////////////////////////////////////
84///Copy constructor
85
86RooBDecay::RooBDecay(const RooBDecay& other, const char* name) :
87 RooAbsAnaConvPdf(other, name),
88 _t("t", this, other._t),
89 _tau("tau", this, other._tau),
90 _dgamma("dgamma", this, other._dgamma),
91 _f0("f0", this, other._f0),
92 _f1("f1", this, other._f1),
93 _f2("f2", this, other._f2),
94 _f3("f3", this, other._f3),
95 _dm("dm", this, other._dm),
96 _basisCosh(other._basisCosh),
97 _basisSinh(other._basisSinh),
98 _basisCos(other._basisCos),
99 _basisSin(other._basisSin),
100 _type(other._type)
101{
102}
103
104////////////////////////////////////////////////////////////////////////////////
105///Destructor
106
108{
109}
110
111////////////////////////////////////////////////////////////////////////////////
112
114{
115 if(basisIndex == _basisCosh)
116 {
117 return _f0;
118 }
119 if(basisIndex == _basisSinh)
120 {
121 return _f1;
122 }
123 if(basisIndex == _basisCos)
124 {
125 return _f2;
126 }
127 if(basisIndex == _basisSin)
128 {
129 return _f3;
130 }
131
132 return 0 ;
133}
134
135////////////////////////////////////////////////////////////////////////////////
136
138{
139 if(basisIndex == _basisCosh)
140 {
141 return _f0.arg().getVariables();
142 }
143 if(basisIndex == _basisSinh)
144 {
145 return _f1.arg().getVariables();
146 }
147 if(basisIndex == _basisCos)
148 {
149 return _f2.arg().getVariables();
150 }
151 if(basisIndex == _basisSin)
152 {
153 return _f3.arg().getVariables();
154 }
155
156 return 0 ;
157}
158
159////////////////////////////////////////////////////////////////////////////////
160
161Int_t RooBDecay::getCoefAnalyticalIntegral(Int_t coef, RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
162{
163 if(coef == _basisCosh)
164 {
165 return _f0.arg().getAnalyticalIntegral(allVars,analVars,rangeName) ;
166 }
167 if(coef == _basisSinh)
168 {
169 return _f1.arg().getAnalyticalIntegral(allVars,analVars,rangeName) ;
170 }
171 if(coef == _basisCos)
172 {
173 return _f2.arg().getAnalyticalIntegral(allVars,analVars,rangeName) ;
174 }
175 if(coef == _basisSin)
176 {
177 return _f3.arg().getAnalyticalIntegral(allVars,analVars,rangeName) ;
178 }
179
180 return 0 ;
181}
182
183////////////////////////////////////////////////////////////////////////////////
184
185Double_t RooBDecay::coefAnalyticalIntegral(Int_t coef, Int_t code, const char* rangeName) const
186{
187 if(coef == _basisCosh)
188 {
189 return _f0.arg().analyticalIntegral(code,rangeName) ;
190 }
191 if(coef == _basisSinh)
192 {
193 return _f1.arg().analyticalIntegral(code,rangeName) ;
194 }
195 if(coef == _basisCos)
196 {
197 return _f2.arg().analyticalIntegral(code,rangeName) ;
198 }
199 if(coef == _basisSin)
200 {
201 return _f3.arg().analyticalIntegral(code,rangeName) ;
202 }
203
204 return 0 ;
205}
206
207////////////////////////////////////////////////////////////////////////////////
208
209Int_t RooBDecay::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
210{
211 if (matchArgs(directVars, generateVars, _t)) return 1;
212 return 0;
213}
214
215////////////////////////////////////////////////////////////////////////////////
216
218{
219 R__ASSERT(code==1);
220 Double_t gammamin = 1/_tau-TMath::Abs(_dgamma)/2;
221 while(1) {
222 Double_t t = -log(RooRandom::uniform())/gammamin;
223 if (_type == Flipped || (_type == DoubleSided && RooRandom::uniform() <0.5) ) t *= -1;
224 if ( t<_t.min() || t>_t.max() ) continue;
225
226 Double_t dgt = _dgamma*t/2;
227 Double_t dmt = _dm*t;
228 Double_t ft = fabs(t);
229 Double_t f = exp(-ft/_tau)*(_f0*cosh(dgt)+_f1*sinh(dgt)+_f2*cos(dmt)+_f3*sin(dmt));
230 if(f < 0) {
231 cout << "RooBDecay::generateEvent(" << GetName() << ") ERROR: PDF value less than zero" << endl;
232 ::abort();
233 }
234 Double_t w = 1.001*exp(-ft*gammamin)*(TMath::Abs(_f0)+TMath::Abs(_f1)+sqrt(_f2*_f2+_f3*_f3));
235 if(w < f) {
236 cout << "RooBDecay::generateEvent(" << GetName() << ") ERROR: Envelope function less than p.d.f. " << endl;
237 cout << _f0 << endl;
238 cout << _f1 << endl;
239 cout << _f2 << endl;
240 cout << _f3 << endl;
241 ::abort();
242 }
243 if(w*RooRandom::uniform() > f) continue;
244 _t = t;
245 break;
246 }
247}
#define f(i)
Definition: RSha256.hxx:104
#define ClassImp(name)
Definition: Rtypes.h:361
#define R__ASSERT(e)
Definition: TError.h:96
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
double cosh(double)
double sinh(double)
double cos(double)
double sqrt(double)
double sin(double)
double exp(double)
double log(double)
RooAbsAnaConvPdf is the base class for PDFs that represent a physics model that can be analytically c...
Int_t declareBasis(const char *expression, const RooArgList &params)
Declare a basis function for use in this physics model.
RooArgSet * getVariables(Bool_t stripDisconnected=kTRUE) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:1909
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().
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
Definition: RooAbsReal.cxx:414
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
Definition: RooAbsReal.cxx:388
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Most general description of B decay time distribution with effects of CP violation,...
Definition: RooBDecay.h:25
Int_t getCoefAnalyticalIntegral(Int_t coef, RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Default implementation of function advertising integration capabilities.
Definition: RooBDecay.cxx:161
RooRealProxy _t
Definition: RooBDecay.h:58
RooArgSet * coefVars(Int_t coefIdx) const
Return set of parameters with are used exclusively by the coefficient functions.
Definition: RooBDecay.cxx:137
virtual ~RooBDecay()
Destructor.
Definition: RooBDecay.cxx:107
RooRealProxy _f2
Definition: RooBDecay.h:63
Int_t _basisCosh
Definition: RooBDecay.h:66
DecayType _type
Definition: RooBDecay.h:71
virtual Double_t coefficient(Int_t basisIndex) const
Definition: RooBDecay.cxx:113
Int_t _basisSin
Definition: RooBDecay.h:69
RooRealProxy _tau
Definition: RooBDecay.h:59
RooBDecay()
Definition: RooBDecay.h:32
RooRealProxy _f0
Definition: RooBDecay.h:61
@ DoubleSided
Definition: RooBDecay.h:29
@ SingleSided
Definition: RooBDecay.h:29
RooRealProxy _dgamma
Definition: RooBDecay.h:60
RooRealProxy _f3
Definition: RooBDecay.h:64
Double_t coefAnalyticalIntegral(Int_t coef, Int_t code, const char *rangeName=0) const
Default implementation of function implementing advertised integrals.
Definition: RooBDecay.cxx:185
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
Definition: RooBDecay.cxx:209
Int_t _basisSinh
Definition: RooBDecay.h:67
RooRealProxy _dm
Definition: RooBDecay.h:65
void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
Definition: RooBDecay.cxx:217
Int_t _basisCos
Definition: RooBDecay.h:68
RooRealProxy _f1
Definition: RooBDecay.h:62
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:83
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
RooResolutionModel is the base class for PDFs that represent a resolution model that can be convolute...
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.
const T & arg() const
Return reference to object held in proxy.
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
TF1 * f1
Definition: legend1.C:11
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Short_t Abs(Short_t d)
Definition: TMathBase.h:120