Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
30 Time-dependent RooAbsAnaConvPdf for CP violating decays
31 to Non-CP eigenstates (eg, \f$ B_0 \rightarrow \rho^\pm \pi^\mp\f$).
32 For a description of the physics model see the
33 BaBar Physics Book, section 6.5.2.3 .
34 The set of CP parameters used in this class is equivalent to
35 the one used in the Physics Book, but it is not exactly the
36 same. Starting from the set in the BaBar Book, in order to
37 get 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]
43 where 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 
56 using namespace std;
57 
58 ClassImp(RooNonCPEigenDecay);
59 
60 #define Debug_RooNonCPEigenDecay 1
61 
62 
63 ////////////////////////////////////////////////////////////////////////////////
64 
65 RooNonCPEigenDecay::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,
76  RooAbsReal& C,
77  RooAbsReal& delC,
78  RooAbsReal& S,
79  RooAbsReal& delS,
80  const RooResolutionModel& model,
81  DecayType type )
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 
123 RooNonCPEigenDecay::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,
138  DecayType type )
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 
182 RooNonCPEigenDecay::RooNonCPEigenDecay( const RooNonCPEigenDecay& other, const char* name )
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 
210 RooNonCPEigenDecay::~RooNonCPEigenDecay( void )
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 
221 Double_t RooNonCPEigenDecay::coefficient( Int_t basisIndex ) const
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 
273 Int_t RooNonCPEigenDecay::getCoefAnalyticalIntegral( Int_t /*code*/, RooArgSet& allVars,
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 
288 Double_t RooNonCPEigenDecay::coefAnalyticalIntegral( Int_t basisIndex,
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 
338 Int_t RooNonCPEigenDecay::getGenerator( const RooArgSet& directVars,
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 
352 void RooNonCPEigenDecay::initGenerator( Int_t code )
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 
365  if (Debug_RooNonCPEigenDecay == 1)
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 
383  if (Debug_RooNonCPEigenDecay == 1)
384  cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: "
385  << _genRhoPlusFrac << endl;
386  }
387 }
388 
389 ////////////////////////////////////////////////////////////////////////////////
390 
391 void RooNonCPEigenDecay::generateEvent( Int_t code )
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 
417  Double_t rand = RooRandom::uniform();
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
436  Double_t expC = coefficient( _basisExp );
437  Double_t sinC = coefficient( _basisSin );
438  Double_t cosC = coefficient( _basisCos );
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 }