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