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 "Riostream.h"
47#include "RooRealVar.h"
48#include "RooRandom.h"
49#include "RooNonCPEigenDecay.h"
50#include "TMath.h"
51#include "RooRealIntegral.h"
52
53using namespace std;
54
56
57#define Debug_RooNonCPEigenDecay 1
58
59
60////////////////////////////////////////////////////////////////////////////////
61
62RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
63 RooRealVar& t,
64 RooAbsCategory& tag,
65 RooAbsReal& tau,
66 RooAbsReal& dm,
67 RooAbsReal& avgW,
68 RooAbsReal& delW,
69 RooAbsCategory& rhoQ,
70 RooAbsReal& correctQ,
71 RooAbsReal& wQ,
72 RooAbsReal& acp,
74 RooAbsReal& delC,
76 RooAbsReal& delS,
77 const RooResolutionModel& model,
79 : RooAbsAnaConvPdf( name, title, model, t ),
80 _acp ( "acp", "acp", this, acp ),
81 _avgC ( "C", "C", this, C ),
82 _delC ( "delC", "delC", this, delC ),
83 _avgS ( "S", "S", this, S ),
84 _delS ( "delS", "delS", this, delS ),
85 _avgW ( "avgW", "Average mistag rate",this, avgW ),
86 _delW ( "delW", "Shift mistag rate", this, delW ),
87 _t ( "t", "time", this, t ),
88 _tau ( "tau", "decay time", this, tau ),
89 _dm ( "dm", "mixing frequency", this, dm ),
90 _tag ( "tag", "CP state", this, tag ),
91 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
92 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
93 _wQ ( "wQ", "mischarge", this, wQ ),
94 _genB0Frac ( 0 ),
95 _genRhoPlusFrac( 0 ),
96 _type ( type )
97{
98 // Constructor
99 switch(type) {
100 case SingleSided:
101 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
102 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
103 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
104 break;
105 case Flipped:
106 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
107 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
108 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
109 break;
110 case DoubleSided:
111 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
112 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
113 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
114 break;
115 }
116}
117
118////////////////////////////////////////////////////////////////////////////////
119
120RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
121 RooRealVar& t,
122 RooAbsCategory& tag,
123 RooAbsReal& tau,
124 RooAbsReal& dm,
125 RooAbsReal& avgW,
126 RooAbsReal& delW,
127 RooAbsCategory& rhoQ,
128 RooAbsReal& correctQ,
129 RooAbsReal& acp,
130 RooAbsReal& C,
131 RooAbsReal& delC,
132 RooAbsReal& S,
133 RooAbsReal& delS,
134 const RooResolutionModel& model,
136 : RooAbsAnaConvPdf( name, title, model, t ),
137 _acp ( "acp", "acp", this, acp ),
138 _avgC ( "C", "C", this, C ),
139 _delC ( "delC", "delC", this, delC ),
140 _avgS ( "S", "S", this, S ),
141 _delS ( "delS", "delS", this, delS ),
142 _avgW ( "avgW", "Average mistag rate",this, avgW ),
143 _delW ( "delW", "Shift mistag rate", this, delW ),
144 _t ( "t", "time", this, t ),
145 _tau ( "tau", "decay time", this, tau ),
146 _dm ( "dm", "mixing frequency", this, dm ),
147 _tag ( "tag", "CP state", this, tag ),
148 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
149 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
150 _genB0Frac ( 0 ),
151 _genRhoPlusFrac( 0 ),
152 _type ( type )
153{
154 // dummy mischarge (must be set to zero!)
155 _wQ = RooRealProxy( "wQ", "mischarge", this, *(new RooRealVar( "wQ", "wQ", 0 )) );
156
157 switch(type) {
158 case SingleSided:
159 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
160 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
161 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
162 break;
163 case Flipped:
164 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
165 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
166 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
167 break;
168 case DoubleSided:
169 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
170 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
171 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
172 break;
173 }
174}
175
176////////////////////////////////////////////////////////////////////////////////
177/// Copy constructor
178
180 : RooAbsAnaConvPdf( other, name ),
181 _acp ( "acp", this, other._acp ),
182 _avgC ( "C", this, other._avgC ),
183 _delC ( "delC", this, other._delC ),
184 _avgS ( "S", this, other._avgS ),
185 _delS ( "delS", this, other._delS ),
186 _avgW ( "avgW", this, other._avgW ),
187 _delW ( "delW", this, other._delW ),
188 _t ( "t", this, other._t ),
189 _tau ( "tau", this, other._tau ),
190 _dm ( "dm", this, other._dm ),
191 _tag ( "tag", this, other._tag ),
192 _rhoQ ( "rhoQ", this, other._rhoQ ),
193 _correctQ ( "correctQ", this, other._correctQ ),
194 _wQ ( "wQ", this, other._wQ ),
195 _genB0Frac ( other._genB0Frac ),
196 _genRhoPlusFrac( other._genRhoPlusFrac ),
197 _type ( other._type ),
198 _basisExp ( other._basisExp ),
199 _basisSin ( other._basisSin ),
200 _basisCos ( other._basisCos )
201{
202}
203
204////////////////////////////////////////////////////////////////////////////////
205/// Destructor
206
208{
209}
210
211////////////////////////////////////////////////////////////////////////////////
212/// - B0 : _tag == -1
213/// - B0bar : _tag == +1
214/// - rho+ : _rhoQ == +1
215/// - rho- : _rhoQ == -1
216/// - the charge correction factor "_correctQ" serves to implement mis-charges
217
218double RooNonCPEigenDecay::coefficient( Int_t basisIndex ) const
219{
220 Int_t rhoQc = _rhoQ * int(_correctQ);
221 assert( rhoQc == 1 || rhoQc == -1 );
222
223 double a_sin_p = _avgS + _delS;
224 double a_sin_m = _avgS - _delS;
225 double a_cos_p = _avgC + _delC;
226 double a_cos_m = _avgC - _delC;
227
228 if (basisIndex == _basisExp) {
229 if (rhoQc == -1 || rhoQc == +1)
230 return (1 + rhoQc*_acp*(1 - 2*_wQ))*(1 + 0.5*_tag*(2*_delW));
231 else
232 return 1;
233 }
234
235 if (basisIndex == _basisSin) {
236
237 if (rhoQc == -1)
238
239 return - ((1 - _acp)*a_sin_m*(1 - _wQ) + (1 + _acp)*a_sin_p*_wQ)*(1 - 2*_avgW)*_tag;
240
241 else if (rhoQc == +1)
242
243 return - ((1 + _acp)*a_sin_p*(1 - _wQ) + (1 - _acp)*a_sin_m*_wQ)*(1 - 2*_avgW)*_tag;
244
245 else
246 return - _tag*((a_sin_p + a_sin_m)/2)*(1 - 2*_avgW);
247 }
248
249 if (basisIndex == _basisCos) {
250
251 if ( rhoQc == -1)
252
253 return + ((1 - _acp)*a_cos_m*(1 - _wQ) + (1 + _acp)*a_cos_p*_wQ)*(1 - 2*_avgW)*_tag;
254
255 else if (rhoQc == +1)
256
257 return + ((1 + _acp)*a_cos_p*(1 - _wQ) + (1 - _acp)*a_cos_m*_wQ)*(1 - 2*_avgW)*_tag;
258
259 else
260 return _tag*((a_cos_p + a_cos_m)/2)*(1 - 2*_avgW);
261 }
262
263 return 0;
264}
265
266// advertise analytical integration
267
268////////////////////////////////////////////////////////////////////////////////
269
271 RooArgSet& analVars, const char* rangeName ) const
272{
273 if (rangeName) return 0 ;
274
275 if (matchArgs( allVars, analVars, _tag, _rhoQ )) return 3;
276 if (matchArgs( allVars, analVars, _rhoQ )) return 2;
277 if (matchArgs( allVars, analVars, _tag )) return 1;
278
279 return 0;
280}
281
282////////////////////////////////////////////////////////////////////////////////
283/// correct for the right/wrong charge...
284
286 Int_t code, const char* /*rangeName*/ ) const
287{
288 Int_t rhoQc = _rhoQ*int(_correctQ);
289
290 double a_sin_p = _avgS + _delS;
291 double a_sin_m = _avgS - _delS;
292 double a_cos_p = _avgC + _delC;
293 double a_cos_m = _avgC - _delC;
294
295 switch(code) {
296
297 // No integration
298 case 0: return coefficient(basisIndex);
299
300 // Integration over 'tag'
301 case 1:
302 if (basisIndex == _basisExp) return 2*(1 + rhoQc*_acp*(1 - 2*_wQ));
303 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
304 assert( false );
305
306 // Integration over 'rhoQ'
307 case 2:
308 if (basisIndex == _basisExp) return 2*(1 + 0.5*_tag*(2.*_delW));
309
310 if (basisIndex == _basisSin)
311
312 return - ( (1 - _acp)*a_sin_m + (1 + _acp)*a_sin_p )*(1 - 2*_avgW)*_tag;
313
314 if (basisIndex == _basisCos)
315
316 return + ( (1 - _acp)*a_cos_m + (1 + _acp)*a_cos_p )*(1 - 2*_avgW)*_tag;
317
318 assert( false );
319
320 // Integration over 'tag' and 'rhoQ'
321 case 3:
322 if (basisIndex == _basisExp) return 2*2; // for both: tag and charge
323 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
324 assert( false );
325
326 default:
327 assert( false );
328 }
329
330 return 0;
331}
332
333////////////////////////////////////////////////////////////////////////////////
334
336 RooArgSet& generateVars, bool staticInitOK ) const
337{
338 if (staticInitOK) {
339 if (matchArgs( directVars, generateVars, _t, _tag, _rhoQ )) return 4;
340 if (matchArgs( directVars, generateVars, _t, _rhoQ )) return 3;
341 if (matchArgs( directVars, generateVars, _t, _tag )) return 2;
342 }
343 if (matchArgs( directVars, generateVars, _t )) return 1;
344 return 0;
345}
346
347////////////////////////////////////////////////////////////////////////////////
348
350{
351 if (code == 2 || code == 4) {
352 // Calculate the fraction of mixed events to generate
353 double sumInt1 = RooRealIntegral( "sumInt1", "sum integral1", *this,
354 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
355 ).getVal();
356 _tag = -1;
357 double b0Int1 = RooRealIntegral( "mixInt1", "mix integral1", *this,
358 RooArgSet( _t.arg(), _rhoQ.arg() )
359 ).getVal();
360 _genB0Frac = b0Int1/sumInt1;
361
363 cout << " o RooNonCPEigenDecay::initgenerator: genB0Frac : "
364 << _genB0Frac
365 << ", tag dilution: " << (1 - 2*_avgW)
366 << endl;
367 }
368
369 if (code == 3 || code == 4) {
370 // Calculate the fraction of positive rho's to generate
371 double sumInt2 = RooRealIntegral( "sumInt2", "sum integral2", *this,
372 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
373 ).getVal();
374 _rhoQ = 1;
375 double b0Int2 = RooRealIntegral( "mixInt2", "mix integral2", *this,
376 RooArgSet( _t.arg(), _tag.arg() )
377 ).getVal();
378 _genRhoPlusFrac = b0Int2/sumInt2;
379
381 cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: "
382 << _genRhoPlusFrac << endl;
383 }
384}
385
386////////////////////////////////////////////////////////////////////////////////
387
389{
390 // Generate delta-t dependent
391 while (true) {
392
393 // B flavor and rho charge (we do not use the integrated weights)
394 if (code != 1) {
395 if (code != 3) _tag = (RooRandom::uniform()<=0.5) ? -1 : +1;
396 if (code != 2) _rhoQ = (RooRandom::uniform()<=0.5) ? 1 : -1;
397 }
398
399 // opposite charge?
400 // Int_t rhoQc = _rhoQ*int(_correctQ);
401
402 double a_sin_p = _avgS + _delS;
403 double a_sin_m = _avgS - _delS;
404 double a_cos_p = _avgC + _delC;
405 double a_cos_m = _avgC - _delC;
406
407 // maximum probability density
408 double a1 = 1 + sqrt(TMath::Power(a_cos_m, 2) + TMath::Power(a_sin_m, 2));
409 double a2 = 1 + sqrt(TMath::Power(a_cos_p, 2) + TMath::Power(a_sin_p, 2));
410
411 double maxAcceptProb = (1.10 + TMath::Abs(_acp)) * (a1 > a2 ? a1 : a2);
412 // The 1.10 in the above line is a security feature to prevent crashes close to the limit at 1.00
413
414 double rand = RooRandom::uniform();
415 double tval(0);
416
417 switch(_type) {
418
419 case SingleSided:
420 tval = -_tau*log(rand);
421 break;
422
423 case Flipped:
424 tval = +_tau*log(rand);
425 break;
426
427 case DoubleSided:
428 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5));
429 break;
430 }
431
432 // get coefficients
433 double expC = coefficient( _basisExp );
434 double sinC = coefficient( _basisSin );
435 double cosC = coefficient( _basisCos );
436
437 // probability density
438 double acceptProb = expC + sinC*sin(_dm*tval) + cosC*cos(_dm*tval);
439
440 // sanity check...
441 assert( acceptProb <= maxAcceptProb );
442
443 // hit or miss...
444 bool accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? true : false;
445
446 if (accept && tval<_t.max() && tval>_t.min()) {
447 _t = tval;
448 break;
449 }
450 }
451}
#define Debug_RooNonCPEigenDecay
RooTemplateProxy< RooAbsReal > RooRealProxy
Compatibility typedef replacing the old RooRealProxy class.
Definition: RooRealProxy.h:23
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition: TGX11.cxx:110
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:645
A space to attach TBranches.
friend class RooRealIntegral
Definition: RooAbsPdf.h:347
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
bool 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:57
Time-dependent RooAbsAnaConvPdf for CP violating decays to Non-CP eigenstates (eg,...
Int_t getCoefAnalyticalIntegral(Int_t coef, RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const override
Default implementation of function advertising integration capabilities.
double coefAnalyticalIntegral(Int_t coef, Int_t code, const char *rangeName=0) const override
correct for the right/wrong charge...
~RooNonCPEigenDecay(void) override
Destructor.
RooCategoryProxy _rhoQ
RooCategoryProxy _tag
void generateEvent(Int_t code) override
Interface for generation of an event using the algorithm corresponding to the specified code.
double coefficient(Int_t basisIndex) const override
Int_t getGenerator(const RooArgSet &directVars, RooArgSet &generateVars, bool staticInitOK=true) const override
Load generatedVars with the subset of directVars that we can generate events for, and return a code t...
void initGenerator(Int_t code) override
Interface for one-time initialization to setup the generator for the specified code.
static double uniform(TRandom *generator=randomGenerator())
Return a number uniformly distributed from (0,1)
Definition: RooRandom.cxx:81
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
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.
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition: RVec.hxx:1759
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1748
RVec< PromoteType< T > > sin(const RVec< T > &v)
Definition: RVec.hxx:1758
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:241
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:718
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123