Logo ROOT  
Reference Guide
RooBCPGenDecay.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * JGS, Jim Smith , University of Colorado, jgsmith@pizero.colorado.edu *
7 * History:
8 * 15-Aug-2002 JGS Created initial version
9 * 11-Sep-2002 JGS Mods to introduce mu (Mirna van Hoek, JGS, Nick Danielson)
10 * *
11 * Copyright (c) 2000-2005, Regents of the University of California, *
12 * University of Colorado *
13 * and Stanford University. All rights reserved. *
14 * *
15 * Redistribution and use in source and binary forms, *
16 * with or without modification, are permitted according to the terms *
17 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
18 *****************************************************************************/
19
20/** \class RooBCPGenDecay
21 \ingroup Roofit
22
23Implement standard CP physics model with S and C (no mention of lambda)
24Suitably stolen and modified from RooBCPEffDecay
25**/
26
27#include "RooRealVar.h"
28#include "RooRandom.h"
29#include "RooBCPGenDecay.h"
30#include "RooRealIntegral.h"
31
32using namespace std;
33
35
36////////////////////////////////////////////////////////////////////////////////
37/// Constructor
38
39RooBCPGenDecay::RooBCPGenDecay(const char *name, const char *title,
41 RooAbsReal& tau, RooAbsReal& dm,
42 RooAbsReal& avgMistag,
44 RooAbsReal& delMistag,
45 RooAbsReal& mu,
46 const RooResolutionModel& model, DecayType type) :
47 RooAbsAnaConvPdf(name,title,model,t),
48 _avgC("C","Coefficient of cos term",this,a),
49 _avgS("S","Coefficient of sin term",this,b),
50 _avgMistag("avgMistag","Average mistag rate",this,avgMistag),
51 _delMistag("delMistag","Delta mistag rate",this,delMistag),
52 _mu("mu","Tagg efficiency difference",this,mu),
53 _t("t","time",this,t),
54 _tau("tau","decay time",this,tau),
55 _dm("dm","mixing frequency",this,dm),
56 _tag("tag","CP state",this,tag),
57 _genB0Frac(0),
58 _type(type)
59{
60 switch(type) {
61 case SingleSided:
62 _basisExp = declareBasis("exp(-@0/@1)",RooArgList(tau,dm)) ;
63 _basisSin = declareBasis("exp(-@0/@1)*sin(@0*@2)",RooArgList(tau,dm)) ;
64 _basisCos = declareBasis("exp(-@0/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
65 break ;
66 case Flipped:
67 _basisExp = declareBasis("exp(@0)/@1)",RooArgList(tau,dm)) ;
68 _basisSin = declareBasis("exp(@0/@1)*sin(@0*@2)",RooArgList(tau,dm)) ;
69 _basisCos = declareBasis("exp(@0/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
70 break ;
71 case DoubleSided:
72 _basisExp = declareBasis("exp(-abs(@0)/@1)",RooArgList(tau,dm)) ;
73 _basisSin = declareBasis("exp(-abs(@0)/@1)*sin(@0*@2)",RooArgList(tau,dm)) ;
74 _basisCos = declareBasis("exp(-abs(@0)/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
75 break ;
76 }
77}
78
79////////////////////////////////////////////////////////////////////////////////
80/// Copy constructor
81
84 _avgC("C",this,other._avgC),
85 _avgS("S",this,other._avgS),
86 _avgMistag("avgMistag",this,other._avgMistag),
87 _delMistag("delMistag",this,other._delMistag),
88 _mu("mu",this,other._mu),
89 _t("t",this,other._t),
90 _tau("tau",this,other._tau),
91 _dm("dm",this,other._dm),
92 _tag("tag",this,other._tag),
93 _genB0Frac(other._genB0Frac),
94 _type(other._type),
95 _basisExp(other._basisExp),
96 _basisSin(other._basisSin),
97 _basisCos(other._basisCos)
98{
99}
100
101////////////////////////////////////////////////////////////////////////////////
102/// Destructor
103
105{
106}
107
108////////////////////////////////////////////////////////////////////////////////
109/// B0 : _tag = +1
110/// B0bar : _tag = -1
111
113{
114 if (basisIndex==_basisExp) {
115 //exp term: (1 -/+ dw + mu*_tag*w)
116 return (1 - _tag*_delMistag + _mu*_tag*(1. - 2.*_avgMistag)) ;
117 // = 1 + _tag*deltaDil/2 + _mu*avgDil
118 }
119
120 if (basisIndex==_basisSin) {
121 //sin term: (+/- (1-2w) + _mu*(1 -/+ delw))*S
122 return (_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgS ;
123 // = (_tag*avgDil + _mu*(1 + tag*deltaDil/2)) * S
124 }
125
126 if (basisIndex==_basisCos) {
127 //cos term: -(+/- (1-2w) + _mu*(1 -/+ delw))*C
128 return -1.*(_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgC ;
129 // = -(_tag*avgDil + _mu*(1 + _tag*deltaDil/2) )* C
130 }
131
132 return 0 ;
133}
134
135////////////////////////////////////////////////////////////////////////////////
136
137Int_t RooBCPGenDecay::getCoefAnalyticalIntegral(Int_t /*code*/, RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
138{
139 if (rangeName) return 0 ;
140 if (matchArgs(allVars,analVars,_tag)) return 1 ;
141 return 0 ;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145
146Double_t RooBCPGenDecay::coefAnalyticalIntegral(Int_t basisIndex, Int_t code, const char* /*rangeName*/) const
147{
148 switch(code) {
149 // No integration
150 case 0: return coefficient(basisIndex) ;
151
152 // Integration over 'tag'
153 case 1:
154 if (basisIndex==_basisExp) {
155 return 2 ;
156 }
157
158 if (basisIndex==_basisSin) {
159 return 2*_mu*_avgS ;
160 }
161 if (basisIndex==_basisCos) {
162 return -2*_mu*_avgC ;
163 }
164 break ;
165
166 default:
167 assert(0) ;
168 }
169
170 return 0 ;
171}
172
173////////////////////////////////////////////////////////////////////////////////
174
175Int_t RooBCPGenDecay::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK) const
176{
177 if (staticInitOK) {
178 if (matchArgs(directVars,generateVars,_t,_tag)) return 2 ;
179 }
180 if (matchArgs(directVars,generateVars,_t)) return 1 ;
181 return 0 ;
182}
183
184////////////////////////////////////////////////////////////////////////////////
185
187{
188 if (code==2) {
189 // Calculate the fraction of mixed events to generate
190 Double_t sumInt = RooRealIntegral("sumInt","sum integral",*this,RooArgSet(_t.arg(),_tag.arg())).getVal() ;
191 _tag = 1 ;
192 Double_t b0Int = RooRealIntegral("mixInt","mix integral",*this,RooArgSet(_t.arg())).getVal() ;
193 _genB0Frac = b0Int/sumInt ;
194 }
195}
196
197////////////////////////////////////////////////////////////////////////////////
198/// Generate mix-state dependent
199
201{
202 if (code==2) {
204 _tag = (rand<=_genB0Frac) ? 1 : -1 ;
205 }
206
207 // Generate delta-t dependent
208 while(1) {
210 Double_t tval(0) ;
211
212 switch(_type) {
213 case SingleSided:
214 tval = -_tau*log(rand);
215 break ;
216 case Flipped:
217 tval= +_tau*log(rand);
218 break ;
219 case DoubleSided:
220 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5)) ;
221 break ;
222 }
223
224 // Accept event if T is in generated range
225 Double_t maxDil = 1.0 ;
226// 2 in next line is conservative and inefficient - allows for delMistag=1!
227 Double_t maxAcceptProb = 2 + fabs(maxDil*_avgS) + fabs(maxDil*_avgC);
228 Double_t acceptProb = (1-_tag*_delMistag + _mu*_tag*(1. - 2.*_avgMistag))
229 + (_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgS*sin(_dm*tval)
230 - (_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgC*cos(_dm*tval);
231
232 Bool_t accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? kTRUE : kFALSE ;
233
234 if (tval<_t.max() && tval>_t.min() && accept) {
235 _t = tval ;
236 break ;
237 }
238 }
239
240}
#define b(i)
Definition: RSha256.hxx:100
const Bool_t kFALSE
Definition: RtypesCore.h:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
double cos(double)
double sin(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.
friend class RooArgSet
Definition: RooAbsArg.h:572
RooAbsCategory is the base class for objects that represent a discrete value with a finite number of ...
friend class RooRealIntegral
Definition: RooAbsPdf.h:313
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().
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
Implement standard CP physics model with S and C (no mention of lambda) Suitably stolen and modified ...
RooRealProxy _mu
RooRealProxy _avgC
virtual Double_t coefAnalyticalIntegral(Int_t coef, Int_t code, const char *rangeName=0) const
Default implementation of function implementing advertised integrals.
virtual Int_t getCoefAnalyticalIntegral(Int_t coef, RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Default implementation of function advertising integration capabilities.
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...
RooRealProxy _avgS
Double_t _genB0Frac
RooRealProxy _dm
void generateEvent(Int_t code)
Generate mix-state dependent.
RooRealProxy _tau
DecayType _type
RooCategoryProxy _tag
RooRealProxy _t
RooRealProxy _avgMistag
RooRealProxy _delMistag
virtual ~RooBCPGenDecay()
Destructor.
virtual Double_t coefficient(Int_t basisIndex) const
B0 : _tag = +1 B0bar : _tag = -1.
void initGenerator(Int_t code)
Interface for one-time initialization to setup the generator for the specified code.
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.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
auto * a
Definition: textangle.C:12