Logo ROOT  
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 "RooRealVar.h"
50#include "RooRandom.h"
51#include "RooNonCPEigenDecay.h"
52#include "TMath.h"
53#include "RooRealIntegral.h"
54
55using namespace std;
56
58
59#define Debug_RooNonCPEigenDecay 1
60
61
62////////////////////////////////////////////////////////////////////////////////
63
64RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
65 RooRealVar& t,
66 RooAbsCategory& tag,
67 RooAbsReal& tau,
68 RooAbsReal& dm,
69 RooAbsReal& avgW,
70 RooAbsReal& delW,
71 RooAbsCategory& rhoQ,
72 RooAbsReal& correctQ,
73 RooAbsReal& wQ,
74 RooAbsReal& acp,
76 RooAbsReal& delC,
78 RooAbsReal& delS,
79 const RooResolutionModel& model,
81 : RooAbsAnaConvPdf( name, title, model, t ),
82 _acp ( "acp", "acp", this, acp ),
83 _avgC ( "C", "C", this, C ),
84 _delC ( "delC", "delC", this, delC ),
85 _avgS ( "S", "S", this, S ),
86 _delS ( "delS", "delS", this, delS ),
87 _avgW ( "avgW", "Average mistag rate",this, avgW ),
88 _delW ( "delW", "Shift mistag rate", this, delW ),
89 _t ( "t", "time", this, t ),
90 _tau ( "tau", "decay time", this, tau ),
91 _dm ( "dm", "mixing frequency", this, dm ),
92 _tag ( "tag", "CP state", this, tag ),
93 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
94 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
95 _wQ ( "wQ", "mischarge", this, wQ ),
96 _genB0Frac ( 0 ),
97 _genRhoPlusFrac( 0 ),
98 _type ( type )
99{
100 // Constructor
101 switch(type) {
102 case SingleSided:
103 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
104 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
105 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
106 break;
107 case Flipped:
108 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
109 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
110 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
111 break;
112 case DoubleSided:
113 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
114 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
115 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
116 break;
117 }
118}
119
120////////////////////////////////////////////////////////////////////////////////
121
122RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
123 RooRealVar& t,
124 RooAbsCategory& tag,
125 RooAbsReal& tau,
126 RooAbsReal& dm,
127 RooAbsReal& avgW,
128 RooAbsReal& delW,
129 RooAbsCategory& rhoQ,
130 RooAbsReal& correctQ,
131 RooAbsReal& acp,
132 RooAbsReal& C,
133 RooAbsReal& delC,
134 RooAbsReal& S,
135 RooAbsReal& delS,
136 const RooResolutionModel& model,
138 : RooAbsAnaConvPdf( name, title, model, t ),
139 _acp ( "acp", "acp", this, acp ),
140 _avgC ( "C", "C", this, C ),
141 _delC ( "delC", "delC", this, delC ),
142 _avgS ( "S", "S", this, S ),
143 _delS ( "delS", "delS", this, delS ),
144 _avgW ( "avgW", "Average mistag rate",this, avgW ),
145 _delW ( "delW", "Shift mistag rate", this, delW ),
146 _t ( "t", "time", this, t ),
147 _tau ( "tau", "decay time", this, tau ),
148 _dm ( "dm", "mixing frequency", this, dm ),
149 _tag ( "tag", "CP state", this, tag ),
150 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
151 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
152 _genB0Frac ( 0 ),
153 _genRhoPlusFrac( 0 ),
154 _type ( type )
155{
156 // dummy mischarge (must be set to zero!)
157 _wQ = RooRealProxy( "wQ", "mischarge", this, *(new RooRealVar( "wQ", "wQ", 0 )) );
158
159 switch(type) {
160 case SingleSided:
161 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
162 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
163 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
164 break;
165 case Flipped:
166 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
167 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
168 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
169 break;
170 case DoubleSided:
171 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
172 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
173 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
174 break;
175 }
176}
177
178////////////////////////////////////////////////////////////////////////////////
179/// Copy constructor
180
182 : RooAbsAnaConvPdf( other, name ),
183 _acp ( "acp", this, other._acp ),
184 _avgC ( "C", this, other._avgC ),
185 _delC ( "delC", this, other._delC ),
186 _avgS ( "S", this, other._avgS ),
187 _delS ( "delS", this, other._delS ),
188 _avgW ( "avgW", this, other._avgW ),
189 _delW ( "delW", this, other._delW ),
190 _t ( "t", this, other._t ),
191 _tau ( "tau", this, other._tau ),
192 _dm ( "dm", this, other._dm ),
193 _tag ( "tag", this, other._tag ),
194 _rhoQ ( "rhoQ", this, other._rhoQ ),
195 _correctQ ( "correctQ", this, other._correctQ ),
196 _wQ ( "wQ", this, other._wQ ),
197 _genB0Frac ( other._genB0Frac ),
198 _genRhoPlusFrac( other._genRhoPlusFrac ),
199 _type ( other._type ),
200 _basisExp ( other._basisExp ),
201 _basisSin ( other._basisSin ),
202 _basisCos ( other._basisCos )
203{
204}
205
206////////////////////////////////////////////////////////////////////////////////
207/// Destructor
208
210{
211}
212
213////////////////////////////////////////////////////////////////////////////////
214/// - B0 : _tag == -1
215/// - B0bar : _tag == +1
216/// - rho+ : _rhoQ == +1
217/// - rho- : _rhoQ == -1
218/// - the charge correction factor "_correctQ" serves to implement mis-charges
219
221{
222 Int_t rhoQc = _rhoQ * int(_correctQ);
223 assert( rhoQc == 1 || rhoQc == -1 );
224
225 Double_t a_sin_p = _avgS + _delS;
226 Double_t a_sin_m = _avgS - _delS;
227 Double_t a_cos_p = _avgC + _delC;
228 Double_t a_cos_m = _avgC - _delC;
229
230 if (basisIndex == _basisExp) {
231 if (rhoQc == -1 || rhoQc == +1)
232 return (1 + rhoQc*_acp*(1 - 2*_wQ))*(1 + 0.5*_tag*(2*_delW));
233 else
234 return 1;
235 }
236
237 if (basisIndex == _basisSin) {
238
239 if (rhoQc == -1)
240
241 return - ((1 - _acp)*a_sin_m*(1 - _wQ) + (1 + _acp)*a_sin_p*_wQ)*(1 - 2*_avgW)*_tag;
242
243 else if (rhoQc == +1)
244
245 return - ((1 + _acp)*a_sin_p*(1 - _wQ) + (1 - _acp)*a_sin_m*_wQ)*(1 - 2*_avgW)*_tag;
246
247 else
248 return - _tag*((a_sin_p + a_sin_m)/2)*(1 - 2*_avgW);
249 }
250
251 if (basisIndex == _basisCos) {
252
253 if ( rhoQc == -1)
254
255 return + ((1 - _acp)*a_cos_m*(1 - _wQ) + (1 + _acp)*a_cos_p*_wQ)*(1 - 2*_avgW)*_tag;
256
257 else if (rhoQc == +1)
258
259 return + ((1 + _acp)*a_cos_p*(1 - _wQ) + (1 - _acp)*a_cos_m*_wQ)*(1 - 2*_avgW)*_tag;
260
261 else
262 return _tag*((a_cos_p + a_cos_m)/2)*(1 - 2*_avgW);
263 }
264
265 return 0;
266}
267
268// advertise analytical integration
269
270////////////////////////////////////////////////////////////////////////////////
271
273 RooArgSet& analVars, const char* rangeName ) const
274{
275 if (rangeName) return 0 ;
276
277 if (matchArgs( allVars, analVars, _tag, _rhoQ )) return 3;
278 if (matchArgs( allVars, analVars, _rhoQ )) return 2;
279 if (matchArgs( allVars, analVars, _tag )) return 1;
280
281 return 0;
282}
283
284////////////////////////////////////////////////////////////////////////////////
285/// correct for the right/wrong charge...
286
288 Int_t code, const char* /*rangeName*/ ) const
289{
290 Int_t rhoQc = _rhoQ*int(_correctQ);
291
292 Double_t a_sin_p = _avgS + _delS;
293 Double_t a_sin_m = _avgS - _delS;
294 Double_t a_cos_p = _avgC + _delC;
295 Double_t a_cos_m = _avgC - _delC;
296
297 switch(code) {
298
299 // No integration
300 case 0: return coefficient(basisIndex);
301
302 // Integration over 'tag'
303 case 1:
304 if (basisIndex == _basisExp) return 2*(1 + rhoQc*_acp*(1 - 2*_wQ));
305 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
306 assert( kFALSE );
307
308 // Integration over 'rhoQ'
309 case 2:
310 if (basisIndex == _basisExp) return 2*(1 + 0.5*_tag*(2.*_delW));
311
312 if (basisIndex == _basisSin)
313
314 return - ( (1 - _acp)*a_sin_m + (1 + _acp)*a_sin_p )*(1 - 2*_avgW)*_tag;
315
316 if (basisIndex == _basisCos)
317
318 return + ( (1 - _acp)*a_cos_m + (1 + _acp)*a_cos_p )*(1 - 2*_avgW)*_tag;
319
320 assert( kFALSE );
321
322 // Integration over 'tag' and 'rhoQ'
323 case 3:
324 if (basisIndex == _basisExp) return 2*2; // for both: tag and charge
325 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
326 assert( kFALSE );
327
328 default:
329 assert( kFALSE );
330 }
331
332 return 0;
333}
334
335////////////////////////////////////////////////////////////////////////////////
336
338 RooArgSet& generateVars, Bool_t staticInitOK ) const
339{
340 if (staticInitOK) {
341 if (matchArgs( directVars, generateVars, _t, _tag, _rhoQ )) return 4;
342 if (matchArgs( directVars, generateVars, _t, _rhoQ )) return 3;
343 if (matchArgs( directVars, generateVars, _t, _tag )) return 2;
344 }
345 if (matchArgs( directVars, generateVars, _t )) return 1;
346 return 0;
347}
348
349////////////////////////////////////////////////////////////////////////////////
350
352{
353 if (code == 2 || code == 4) {
354 // Calculate the fraction of mixed events to generate
355 Double_t sumInt1 = RooRealIntegral( "sumInt1", "sum integral1", *this,
356 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
357 ).getVal();
358 _tag = -1;
359 Double_t b0Int1 = RooRealIntegral( "mixInt1", "mix integral1", *this,
360 RooArgSet( _t.arg(), _rhoQ.arg() )
361 ).getVal();
362 _genB0Frac = b0Int1/sumInt1;
363
365 cout << " o RooNonCPEigenDecay::initgenerator: genB0Frac : "
366 << _genB0Frac
367 << ", tag dilution: " << (1 - 2*_avgW)
368 << endl;
369 }
370
371 if (code == 3 || code == 4) {
372 // Calculate the fraction of positive rho's to generate
373 Double_t sumInt2 = RooRealIntegral( "sumInt2", "sum integral2", *this,
374 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
375 ).getVal();
376 _rhoQ = 1;
377 Double_t b0Int2 = RooRealIntegral( "mixInt2", "mix integral2", *this,
378 RooArgSet( _t.arg(), _tag.arg() )
379 ).getVal();
380 _genRhoPlusFrac = b0Int2/sumInt2;
381
383 cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: "
384 << _genRhoPlusFrac << endl;
385 }
386}
387
388////////////////////////////////////////////////////////////////////////////////
389
391{
392 // Generate delta-t dependent
393 while (kTRUE) {
394
395 // B flavor and rho charge (we do not use the integrated weights)
396 if (code != 1) {
397 if (code != 3) _tag = (RooRandom::uniform()<=0.5) ? -1 : +1;
398 if (code != 2) _rhoQ = (RooRandom::uniform()<=0.5) ? 1 : -1;
399 }
400
401 // opposite charge?
402 // Int_t rhoQc = _rhoQ*int(_correctQ);
403
404 Double_t a_sin_p = _avgS + _delS;
405 Double_t a_sin_m = _avgS - _delS;
406 Double_t a_cos_p = _avgC + _delC;
407 Double_t a_cos_m = _avgC - _delC;
408
409 // maximum probability density
410 double a1 = 1 + sqrt(TMath::Power(a_cos_m, 2) + TMath::Power(a_sin_m, 2));
411 double a2 = 1 + sqrt(TMath::Power(a_cos_p, 2) + TMath::Power(a_sin_p, 2));
412
413 Double_t maxAcceptProb = (1.10 + TMath::Abs(_acp)) * (a1 > a2 ? a1 : a2);
414 // The 1.10 in the above line is a security feature to prevent crashes close to the limit at 1.00
415
417 Double_t tval(0);
418
419 switch(_type) {
420
421 case SingleSided:
422 tval = -_tau*log(rand);
423 break;
424
425 case Flipped:
426 tval = +_tau*log(rand);
427 break;
428
429 case DoubleSided:
430 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5));
431 break;
432 }
433
434 // get coefficients
438
439 // probability density
440 Double_t acceptProb = expC + sinC*sin(_dm*tval) + cosC*cos(_dm*tval);
441
442 // sanity check...
443 assert( acceptProb <= maxAcceptProb );
444
445 // hit or miss...
446 Bool_t accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? kTRUE : kFALSE;
447
448 if (accept && tval<_t.max() && tval>_t.min()) {
449 _t = tval;
450 break;
451 }
452 }
453}
#define Debug_RooNonCPEigenDecay
RooTemplateProxy< RooAbsReal > RooRealProxy
Compatibility typedef replacing the old RooRealProxy class.
Definition: RooRealProxy.h:23
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
char name[80]
Definition: TGX11.cxx:110
int type
Definition: TGX11.cxx:121
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:624
RooAbsCategory is the base class for objects that represent a discrete value with a finite number of ...
friend class RooRealIntegral
Definition: RooAbsPdf.h:357
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:63
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:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:35
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:83
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
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.
static double C[]
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
RooArgSet S(Args_t &&... args)
Definition: RooArgSet.h:219
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:685
Short_t Abs(Short_t d)
Definition: TMathBase.h:120