Logo ROOT   6.08/07
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 /**
21 \file RooBCPGenDecay.cxx
22 \class RooBCPGenDecay
23 \ingroup Roofit
24 
25 Implement standard CP physics model with S and C (no mention of lambda)
26 Suitably stolen and modified from RooBCPEffDecay
27 **/
28 
29 #include "RooRealVar.h"
30 #include "RooRandom.h"
31 #include "RooBCPGenDecay.h"
32 #include "RooRealIntegral.h"
33 
34 using namespace std;
35 
37 ;
38 
39 
40 
41 ////////////////////////////////////////////////////////////////////////////////
42 /// Constructor
43 
44 RooBCPGenDecay::RooBCPGenDecay(const char *name, const char *title,
45  RooRealVar& t, RooAbsCategory& tag,
47  RooAbsReal& avgMistag,
49  RooAbsReal& delMistag,
50  RooAbsReal& mu,
52  RooAbsAnaConvPdf(name,title,model,t),
53  _avgC("C","Coefficient of cos term",this,a),
54  _avgS("S","Coefficient of sin term",this,b),
55  _avgMistag("avgMistag","Average mistag rate",this,avgMistag),
56  _delMistag("delMistag","Delta mistag rate",this,delMistag),
57  _mu("mu","Tagg efficiency difference",this,mu),
58  _t("t","time",this,t),
59  _tau("tau","decay time",this,tau),
60  _dm("dm","mixing frequency",this,dm),
61  _tag("tag","CP state",this,tag),
62  _genB0Frac(0),
63  _type(type)
64 {
65  switch(type) {
66  case SingleSided:
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 Flipped:
72  _basisExp = declareBasis("exp(@0)/@1)",RooArgList(tau,dm)) ;
73  _basisSin = declareBasis("exp(@0/@1)*sin(@0*@2)",RooArgList(tau,dm)) ;
74  _basisCos = declareBasis("exp(@0/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
75  break ;
76  case DoubleSided:
77  _basisExp = declareBasis("exp(-abs(@0)/@1)",RooArgList(tau,dm)) ;
78  _basisSin = declareBasis("exp(-abs(@0)/@1)*sin(@0*@2)",RooArgList(tau,dm)) ;
79  _basisCos = declareBasis("exp(-abs(@0)/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
80  break ;
81  }
82 }
83 
84 
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// Copy constructor
88 
89 RooBCPGenDecay::RooBCPGenDecay(const RooBCPGenDecay& other, const char* name) :
90  RooAbsAnaConvPdf(other,name),
91  _avgC("C",this,other._avgC),
92  _avgS("S",this,other._avgS),
93  _avgMistag("avgMistag",this,other._avgMistag),
94  _delMistag("delMistag",this,other._delMistag),
95  _mu("mu",this,other._mu),
96  _t("t",this,other._t),
97  _tau("tau",this,other._tau),
98  _dm("dm",this,other._dm),
99  _tag("tag",this,other._tag),
100  _genB0Frac(other._genB0Frac),
101  _type(other._type),
102  _basisExp(other._basisExp),
103  _basisSin(other._basisSin),
104  _basisCos(other._basisCos)
105 {
106 }
107 
108 
109 
110 ////////////////////////////////////////////////////////////////////////////////
111 /// Destructor
112 
114 {
115 }
116 
117 
118 
119 ////////////////////////////////////////////////////////////////////////////////
120 /// B0 : _tag = +1
121 /// B0bar : _tag = -1
122 
124 {
125  if (basisIndex==_basisExp) {
126  //exp term: (1 -/+ dw + mu*_tag*w)
127  return (1 - _tag*_delMistag + _mu*_tag*(1. - 2.*_avgMistag)) ;
128  // = 1 + _tag*deltaDil/2 + _mu*avgDil
129  }
130 
131  if (basisIndex==_basisSin) {
132  //sin term: (+/- (1-2w) + _mu*(1 -/+ delw))*S
133  return (_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgS ;
134  // = (_tag*avgDil + _mu*(1 + tag*deltaDil/2)) * S
135  }
136 
137  if (basisIndex==_basisCos) {
138  //cos term: -(+/- (1-2w) + _mu*(1 -/+ delw))*C
139  return -1.*(_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgC ;
140  // = -(_tag*avgDil + _mu*(1 + _tag*deltaDil/2) )* C
141  }
142 
143  return 0 ;
144 }
145 
146 
147 
148 ////////////////////////////////////////////////////////////////////////////////
149 
150 Int_t RooBCPGenDecay::getCoefAnalyticalIntegral(Int_t /*code*/, RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const
151 {
152  if (rangeName) return 0 ;
153  if (matchArgs(allVars,analVars,_tag)) return 1 ;
154  return 0 ;
155 }
156 
157 
158 
159 ////////////////////////////////////////////////////////////////////////////////
160 
161 Double_t RooBCPGenDecay::coefAnalyticalIntegral(Int_t basisIndex, Int_t code, const char* /*rangeName*/) const
162 {
163  switch(code) {
164  // No integration
165  case 0: return coefficient(basisIndex) ;
166 
167  // Integration over 'tag'
168  case 1:
169  if (basisIndex==_basisExp) {
170  return 2 ;
171  }
172 
173  if (basisIndex==_basisSin) {
174  return 2*_mu*_avgS ;
175  }
176  if (basisIndex==_basisCos) {
177  return -2*_mu*_avgC ;
178  }
179  break ;
180 
181  default:
182  assert(0) ;
183  }
184 
185  return 0 ;
186 }
187 
188 
189 
190 ////////////////////////////////////////////////////////////////////////////////
191 
192 Int_t RooBCPGenDecay::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK) const
193 {
194  if (staticInitOK) {
195  if (matchArgs(directVars,generateVars,_t,_tag)) return 2 ;
196  }
197  if (matchArgs(directVars,generateVars,_t)) return 1 ;
198  return 0 ;
199 }
200 
201 
202 
203 ////////////////////////////////////////////////////////////////////////////////
204 
206 {
207  if (code==2) {
208  // Calculate the fraction of mixed events to generate
209  Double_t sumInt = RooRealIntegral("sumInt","sum integral",*this,RooArgSet(_t.arg(),_tag.arg())).getVal() ;
210  _tag = 1 ;
211  Double_t b0Int = RooRealIntegral("mixInt","mix integral",*this,RooArgSet(_t.arg())).getVal() ;
212  _genB0Frac = b0Int/sumInt ;
213  }
214 }
215 
216 
217 
218 ////////////////////////////////////////////////////////////////////////////////
219 /// Generate mix-state dependent
220 
222 {
223  if (code==2) {
224  Double_t rand = RooRandom::uniform() ;
225  _tag = (rand<=_genB0Frac) ? 1 : -1 ;
226  }
227 
228  // Generate delta-t dependent
229  while(1) {
230  Double_t rand = RooRandom::uniform() ;
231  Double_t tval(0) ;
232 
233  switch(_type) {
234  case SingleSided:
235  tval = -_tau*log(rand);
236  break ;
237  case Flipped:
238  tval= +_tau*log(rand);
239  break ;
240  case DoubleSided:
241  tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5)) ;
242  break ;
243  }
244 
245  // Accept event if T is in generated range
246  Double_t maxDil = 1.0 ;
247 // 2 in next line is conservative and inefficient - allows for delMistag=1!
248  Double_t maxAcceptProb = 2 + fabs(maxDil*_avgS) + fabs(maxDil*_avgC);
249  Double_t acceptProb = (1-_tag*_delMistag + _mu*_tag*(1. - 2.*_avgMistag))
250  + (_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgS*sin(_dm*tval)
251  - (_tag*(1-2*_avgMistag) + _mu*(1. - _tag*_delMistag))*_avgC*cos(_dm*tval);
252 
253  Bool_t accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? kTRUE : kFALSE ;
254 
255  if (tval<_t.max() && tval>_t.min() && accept) {
256  _t = tval ;
257  break ;
258  }
259  }
260 
261 }
262 
RooRealProxy _mu
Implement standard CP physics model with S and C (no mention of lambda) Suitably stolen and modified ...
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooRealProxy _avgC
RooRealProxy _dm
virtual ~RooBCPGenDecay()
Destructor.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
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...
STL namespace.
double cos(double)
virtual Double_t coefficient(Int_t basisIndex) const
B0 : _tag = +1 B0bar : _tag = -1.
Double_t _genB0Frac
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
RooRealProxy _tau
const RooAbsCategory & arg() const
RooCategoryProxy _tag
friend class RooArgSet
Definition: RooAbsArg.h:469
double sin(double)
DecayType _type
RooRealProxy _delMistag
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void initGenerator(Int_t code)
Interface for one-time initialization to setup the generator for the specified code.
#define ClassImp(name)
Definition: Rtypes.h:279
static Double_t uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:84
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooRealProxy _avgMistag
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
void generateEvent(Int_t code)
Generate mix-state dependent.
int type
Definition: TGX11.cxx:120
Int_t declareBasis(const char *expression, const RooArgList &params)
Declare a basis function for use in this physics model.
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
friend class RooRealIntegral
Definition: RooAbsPdf.h:294
RooAbsCategory is the common abstract base class for objects that represent a discrete value with a f...
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 Double_t b
Definition: TRolke.cxx:630
virtual Double_t coefAnalyticalIntegral(Int_t coef, Int_t code, const char *rangeName=0) const
Default implementation of function implementing advertised integrals.
RooRealProxy _t
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
const Bool_t kTRUE
Definition: Rtypes.h:91
virtual Int_t getCoefAnalyticalIntegral(Int_t coef, RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Default implementation of function advertising integration capabilities.
char name[80]
Definition: TGX11.cxx:109
RooRealProxy _avgS
double log(double)