Logo ROOT   6.18/05
Reference Guide
RooNonCPEigenDecay.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitModels *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * AH, Andreas Hoecker, Orsay, hoecker@slac.stanford.edu *
7 * SL, Sandrine Laplace, Orsay, laplace@slac.stanford.edu *
8 * JS, Jan Stark, Paris, stark@slac.stanford.edu *
9 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
10 * *
11 * Copyright (c) 2000-2005, Regents of the University of California, *
12 * IN2P3. All rights reserved. *
13 * *
14 * History *
15 * Nov-2001 WV Created initial version *
16 * Dec-2001 SL mischarge correction, direct CPV *
17 * Jan-2002 AH built dedicated generator + code cleaning *
18 * Mar-2002 JS committed debugged version to CVS *
19 * Apr-2002 AH allow prompt (ie, non-Pdf) mischarge treatment *
20 * May-2002 JS Changed the set of CP parameters (mathematically equiv.) *
21 * *
22 * Redistribution and use in source and binary forms, *
23 * with or without modification, are permitted according to the terms *
24 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
25 *****************************************************************************/
26
27/** \class RooNonCPEigenDecay
28 \ingroup Roofit
29
30Time-dependent RooAbsAnaConvPdf for CP violating decays
31to Non-CP eigenstates (eg, \f$ B_0 \rightarrow \rho^\pm \pi^\mp\f$).
32For a description of the physics model see the
33BaBar Physics Book, section 6.5.2.3 .
34The set of CP parameters used in this class is equivalent to
35the one used in the Physics Book, but it is not exactly the
36same. Starting from the set in the BaBar Book, in order to
37get the parameters used here you have to change the sign of both
38\f$a_c^+\f$ and \f$a_c^-\f$, and then substitute:
39\f[
40 a_s^Q = S + Q \cdot \delta S \\
41 a_c^Q = C + Q \cdot \delta C
42\f]
43where Q denotes the charge of the \f$\rho\f$ meson.
44**/
45
46#include "RooFit.h"
47
48#include "Riostream.h"
49#include "Riostream.h"
50#include "RooRealVar.h"
51#include "RooRandom.h"
52#include "RooNonCPEigenDecay.h"
53#include "TMath.h"
54#include "RooRealIntegral.h"
55
56using namespace std;
57
59
60#define Debug_RooNonCPEigenDecay 1
61
62
63////////////////////////////////////////////////////////////////////////////////
64
65RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
66 RooRealVar& t,
67 RooAbsCategory& tag,
68 RooAbsReal& tau,
69 RooAbsReal& dm,
70 RooAbsReal& avgW,
71 RooAbsReal& delW,
72 RooAbsCategory& rhoQ,
73 RooAbsReal& correctQ,
74 RooAbsReal& wQ,
75 RooAbsReal& acp,
77 RooAbsReal& delC,
79 RooAbsReal& delS,
80 const RooResolutionModel& model,
82 : RooAbsAnaConvPdf( name, title, model, t ),
83 _acp ( "acp", "acp", this, acp ),
84 _avgC ( "C", "C", this, C ),
85 _delC ( "delC", "delC", this, delC ),
86 _avgS ( "S", "S", this, S ),
87 _delS ( "delS", "delS", this, delS ),
88 _avgW ( "avgW", "Average mistag rate",this, avgW ),
89 _delW ( "delW", "Shift mistag rate", this, delW ),
90 _t ( "t", "time", this, t ),
91 _tau ( "tau", "decay time", this, tau ),
92 _dm ( "dm", "mixing frequency", this, dm ),
93 _tag ( "tag", "CP state", this, tag ),
94 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
95 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
96 _wQ ( "wQ", "mischarge", this, wQ ),
97 _genB0Frac ( 0 ),
98 _genRhoPlusFrac( 0 ),
99 _type ( type )
100{
101 // Constructor
102 switch(type) {
103 case SingleSided:
104 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
105 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
106 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
107 break;
108 case Flipped:
109 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
110 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
111 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
112 break;
113 case DoubleSided:
114 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
115 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
116 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
117 break;
118 }
119}
120
121////////////////////////////////////////////////////////////////////////////////
122
123RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
124 RooRealVar& t,
125 RooAbsCategory& tag,
126 RooAbsReal& tau,
127 RooAbsReal& dm,
128 RooAbsReal& avgW,
129 RooAbsReal& delW,
130 RooAbsCategory& rhoQ,
131 RooAbsReal& correctQ,
132 RooAbsReal& acp,
133 RooAbsReal& C,
134 RooAbsReal& delC,
135 RooAbsReal& S,
136 RooAbsReal& delS,
137 const RooResolutionModel& model,
139 : RooAbsAnaConvPdf( name, title, model, t ),
140 _acp ( "acp", "acp", this, acp ),
141 _avgC ( "C", "C", this, C ),
142 _delC ( "delC", "delC", this, delC ),
143 _avgS ( "S", "S", this, S ),
144 _delS ( "delS", "delS", this, delS ),
145 _avgW ( "avgW", "Average mistag rate",this, avgW ),
146 _delW ( "delW", "Shift mistag rate", this, delW ),
147 _t ( "t", "time", this, t ),
148 _tau ( "tau", "decay time", this, tau ),
149 _dm ( "dm", "mixing frequency", this, dm ),
150 _tag ( "tag", "CP state", this, tag ),
151 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
152 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
153 _genB0Frac ( 0 ),
154 _genRhoPlusFrac( 0 ),
155 _type ( type )
156{
157 // dummy mischarge (must be set to zero!)
158 _wQ = RooRealProxy( "wQ", "mischarge", this, *(new RooRealVar( "wQ", "wQ", 0 )) );
159
160 switch(type) {
161 case SingleSided:
162 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
163 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
164 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
165 break;
166 case Flipped:
167 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
168 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
169 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
170 break;
171 case DoubleSided:
172 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
173 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
174 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
175 break;
176 }
177}
178
179////////////////////////////////////////////////////////////////////////////////
180/// Copy constructor
181
183 : RooAbsAnaConvPdf( other, name ),
184 _acp ( "acp", this, other._acp ),
185 _avgC ( "C", this, other._avgC ),
186 _delC ( "delC", this, other._delC ),
187 _avgS ( "S", this, other._avgS ),
188 _delS ( "delS", this, other._delS ),
189 _avgW ( "avgW", this, other._avgW ),
190 _delW ( "delW", this, other._delW ),
191 _t ( "t", this, other._t ),
192 _tau ( "tau", this, other._tau ),
193 _dm ( "dm", this, other._dm ),
194 _tag ( "tag", this, other._tag ),
195 _rhoQ ( "rhoQ", this, other._rhoQ ),
196 _correctQ ( "correctQ", this, other._correctQ ),
197 _wQ ( "wQ", this, other._wQ ),
198 _genB0Frac ( other._genB0Frac ),
199 _genRhoPlusFrac( other._genRhoPlusFrac ),
200 _type ( other._type ),
201 _basisExp ( other._basisExp ),
202 _basisSin ( other._basisSin ),
203 _basisCos ( other._basisCos )
204{
205}
206
207////////////////////////////////////////////////////////////////////////////////
208/// Destructor
209
211{
212}
213
214////////////////////////////////////////////////////////////////////////////////
215/// - B0 : _tag == -1
216/// - B0bar : _tag == +1
217/// - rho+ : _rhoQ == +1
218/// - rho- : _rhoQ == -1
219/// - the charge correction factor "_correctQ" serves to implement mis-charges
220
222{
223 Int_t rhoQc = _rhoQ * int(_correctQ);
224 assert( rhoQc == 1 || rhoQc == -1 );
225
226 Double_t a_sin_p = _avgS + _delS;
227 Double_t a_sin_m = _avgS - _delS;
228 Double_t a_cos_p = _avgC + _delC;
229 Double_t a_cos_m = _avgC - _delC;
230
231 if (basisIndex == _basisExp) {
232 if (rhoQc == -1 || rhoQc == +1)
233 return (1 + rhoQc*_acp*(1 - 2*_wQ))*(1 + 0.5*_tag*(2*_delW));
234 else
235 return 1;
236 }
237
238 if (basisIndex == _basisSin) {
239
240 if (rhoQc == -1)
241
242 return - ((1 - _acp)*a_sin_m*(1 - _wQ) + (1 + _acp)*a_sin_p*_wQ)*(1 - 2*_avgW)*_tag;
243
244 else if (rhoQc == +1)
245
246 return - ((1 + _acp)*a_sin_p*(1 - _wQ) + (1 - _acp)*a_sin_m*_wQ)*(1 - 2*_avgW)*_tag;
247
248 else
249 return - _tag*((a_sin_p + a_sin_m)/2)*(1 - 2*_avgW);
250 }
251
252 if (basisIndex == _basisCos) {
253
254 if ( rhoQc == -1)
255
256 return + ((1 - _acp)*a_cos_m*(1 - _wQ) + (1 + _acp)*a_cos_p*_wQ)*(1 - 2*_avgW)*_tag;
257
258 else if (rhoQc == +1)
259
260 return + ((1 + _acp)*a_cos_p*(1 - _wQ) + (1 - _acp)*a_cos_m*_wQ)*(1 - 2*_avgW)*_tag;
261
262 else
263 return _tag*((a_cos_p + a_cos_m)/2)*(1 - 2*_avgW);
264 }
265
266 return 0;
267}
268
269// advertise analytical integration
270
271////////////////////////////////////////////////////////////////////////////////
272
274 RooArgSet& analVars, const char* rangeName ) const
275{
276 if (rangeName) return 0 ;
277
278 if (matchArgs( allVars, analVars, _tag, _rhoQ )) return 3;
279 if (matchArgs( allVars, analVars, _rhoQ )) return 2;
280 if (matchArgs( allVars, analVars, _tag )) return 1;
281
282 return 0;
283}
284
285////////////////////////////////////////////////////////////////////////////////
286/// correct for the right/wrong charge...
287
289 Int_t code, const char* /*rangeName*/ ) const
290{
291 Int_t rhoQc = _rhoQ*int(_correctQ);
292
293 Double_t a_sin_p = _avgS + _delS;
294 Double_t a_sin_m = _avgS - _delS;
295 Double_t a_cos_p = _avgC + _delC;
296 Double_t a_cos_m = _avgC - _delC;
297
298 switch(code) {
299
300 // No integration
301 case 0: return coefficient(basisIndex);
302
303 // Integration over 'tag'
304 case 1:
305 if (basisIndex == _basisExp) return 2*(1 + rhoQc*_acp*(1 - 2*_wQ));
306 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
307 assert( kFALSE );
308
309 // Integration over 'rhoQ'
310 case 2:
311 if (basisIndex == _basisExp) return 2*(1 + 0.5*_tag*(2.*_delW));
312
313 if (basisIndex == _basisSin)
314
315 return - ( (1 - _acp)*a_sin_m + (1 + _acp)*a_sin_p )*(1 - 2*_avgW)*_tag;
316
317 if (basisIndex == _basisCos)
318
319 return + ( (1 - _acp)*a_cos_m + (1 + _acp)*a_cos_p )*(1 - 2*_avgW)*_tag;
320
321 assert( kFALSE );
322
323 // Integration over 'tag' and 'rhoQ'
324 case 3:
325 if (basisIndex == _basisExp) return 2*2; // for both: tag and charge
326 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
327 assert( kFALSE );
328
329 default:
330 assert( kFALSE );
331 }
332
333 return 0;
334}
335
336////////////////////////////////////////////////////////////////////////////////
337
339 RooArgSet& generateVars, Bool_t staticInitOK ) const
340{
341 if (staticInitOK) {
342 if (matchArgs( directVars, generateVars, _t, _tag, _rhoQ )) return 4;
343 if (matchArgs( directVars, generateVars, _t, _rhoQ )) return 3;
344 if (matchArgs( directVars, generateVars, _t, _tag )) return 2;
345 }
346 if (matchArgs( directVars, generateVars, _t )) return 1;
347 return 0;
348}
349
350////////////////////////////////////////////////////////////////////////////////
351
353{
354 if (code == 2 || code == 4) {
355 // Calculate the fraction of mixed events to generate
356 Double_t sumInt1 = RooRealIntegral( "sumInt1", "sum integral1", *this,
357 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
358 ).getVal();
359 _tag = -1;
360 Double_t b0Int1 = RooRealIntegral( "mixInt1", "mix integral1", *this,
361 RooArgSet( _t.arg(), _rhoQ.arg() )
362 ).getVal();
363 _genB0Frac = b0Int1/sumInt1;
364
366 cout << " o RooNonCPEigenDecay::initgenerator: genB0Frac : "
367 << _genB0Frac
368 << ", tag dilution: " << (1 - 2*_avgW)
369 << endl;
370 }
371
372 if (code == 3 || code == 4) {
373 // Calculate the fraction of positive rho's to generate
374 Double_t sumInt2 = RooRealIntegral( "sumInt2", "sum integral2", *this,
375 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
376 ).getVal();
377 _rhoQ = 1;
378 Double_t b0Int2 = RooRealIntegral( "mixInt2", "mix integral2", *this,
379 RooArgSet( _t.arg(), _tag.arg() )
380 ).getVal();
381 _genRhoPlusFrac = b0Int2/sumInt2;
382
384 cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: "
385 << _genRhoPlusFrac << endl;
386 }
387}
388
389////////////////////////////////////////////////////////////////////////////////
390
392{
393 // Generate delta-t dependent
394 while (kTRUE) {
395
396 // B flavor and rho charge (we do not use the integrated weights)
397 if (code != 1) {
398 if (code != 3) _tag = (RooRandom::uniform()<=0.5) ? -1 : +1;
399 if (code != 2) _rhoQ = (RooRandom::uniform()<=0.5) ? 1 : -1;
400 }
401
402 // opposite charge?
403 // Int_t rhoQc = _rhoQ*int(_correctQ);
404
405 Double_t a_sin_p = _avgS + _delS;
406 Double_t a_sin_m = _avgS - _delS;
407 Double_t a_cos_p = _avgC + _delC;
408 Double_t a_cos_m = _avgC - _delC;
409
410 // maximum probability density
411 double a1 = 1 + sqrt(TMath::Power(a_cos_m, 2) + TMath::Power(a_sin_m, 2));
412 double a2 = 1 + sqrt(TMath::Power(a_cos_p, 2) + TMath::Power(a_sin_p, 2));
413
414 Double_t maxAcceptProb = (1.10 + TMath::Abs(_acp)) * (a1 > a2 ? a1 : a2);
415 // The 1.10 in the above line is a security feature to prevent crashes close to the limit at 1.00
416
418 Double_t tval(0);
419
420 switch(_type) {
421
422 case SingleSided:
423 tval = -_tau*log(rand);
424 break;
425
426 case Flipped:
427 tval = +_tau*log(rand);
428 break;
429
430 case DoubleSided:
431 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5));
432 break;
433 }
434
435 // get coefficients
439
440 // probability density
441 Double_t acceptProb = expC + sinC*sin(_dm*tval) + cosC*cos(_dm*tval);
442
443 // sanity check...
444 assert( acceptProb <= maxAcceptProb );
445
446 // hit or miss...
447 Bool_t accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? kTRUE : kFALSE;
448
449 if (accept && tval<_t.max() && tval>_t.min()) {
450 _t = tval;
451 break;
452 }
453 }
454}
#define Debug_RooNonCPEigenDecay
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
char name[80]
Definition: TGX11.cxx:109
int type
Definition: TGX11.cxx:120
double cos(double)
double sqrt(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:516
RooAbsCategory is the common abstract base class for objects that represent a discrete value with a f...
friend class RooRealIntegral
Definition: RooAbsPdf.h:308
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
friend class RooRealProxy
Definition: RooAbsReal.h:422
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
const RooAbsCategory & arg() const
Time-dependent RooAbsAnaConvPdf for CP violating decays to Non-CP eigenstates (eg,...
virtual Double_t coefAnalyticalIntegral(Int_t coef, Int_t code, const char *rangeName=0) const
correct for the right/wrong charge...
RooCategoryProxy _rhoQ
void generateEvent(Int_t code)
Interface for generation of an event using the algorithm corresponding to the specified code.
RooCategoryProxy _tag
void initGenerator(Int_t code)
Interface for one-time initialization to setup the generator for the specified code.
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...
virtual Double_t coefficient(Int_t basisIndex) const
virtual Int_t getCoefAnalyticalIntegral(Int_t coef, RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Default implementation of function advertising integration capabilities.
virtual ~RooNonCPEigenDecay(void)
Destructor.
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
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
static double C[]
RooArgSet S(const RooAbsArg &v1)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:723
Short_t Abs(Short_t d)
Definition: TMathBase.h:120