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 _wQ ( "wQ", "mischarge", this, *(new RooRealVar( "wQ", "wQ", 0 )), true, false, true ),
151 _genB0Frac ( 0 ),
152 _genRhoPlusFrac( 0 ),
153 _type ( type )
154{
155 switch(type) {
156 case SingleSided:
157 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
158 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
159 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
160 break;
161 case Flipped:
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 DoubleSided:
167 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
168 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
169 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
170 break;
171 }
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Copy constructor
176
178 : RooAbsAnaConvPdf( other, name ),
179 _acp ( "acp", this, other._acp ),
180 _avgC ( "C", this, other._avgC ),
181 _delC ( "delC", this, other._delC ),
182 _avgS ( "S", this, other._avgS ),
183 _delS ( "delS", this, other._delS ),
184 _avgW ( "avgW", this, other._avgW ),
185 _delW ( "delW", this, other._delW ),
186 _t ( "t", this, other._t ),
187 _tau ( "tau", this, other._tau ),
188 _dm ( "dm", this, other._dm ),
189 _tag ( "tag", this, other._tag ),
190 _rhoQ ( "rhoQ", this, other._rhoQ ),
191 _correctQ ( "correctQ", this, other._correctQ ),
192 _wQ ( "wQ", this, other._wQ ),
193 _genB0Frac ( other._genB0Frac ),
194 _genRhoPlusFrac( other._genRhoPlusFrac ),
195 _type ( other._type ),
196 _basisExp ( other._basisExp ),
197 _basisSin ( other._basisSin ),
198 _basisCos ( other._basisCos )
199{
200}
201
202////////////////////////////////////////////////////////////////////////////////
203/// Destructor
204
206{
207}
208
209////////////////////////////////////////////////////////////////////////////////
210/// - B0 : _tag == -1
211/// - B0bar : _tag == +1
212/// - rho+ : _rhoQ == +1
213/// - rho- : _rhoQ == -1
214/// - the charge correction factor "_correctQ" serves to implement mis-charges
215
216double RooNonCPEigenDecay::coefficient( Int_t basisIndex ) const
217{
218 Int_t rhoQc = _rhoQ * int(_correctQ);
219 assert( rhoQc == 1 || rhoQc == -1 );
220
221 double a_sin_p = _avgS + _delS;
222 double a_sin_m = _avgS - _delS;
223 double a_cos_p = _avgC + _delC;
224 double a_cos_m = _avgC - _delC;
225
226 if (basisIndex == _basisExp) {
227 if (rhoQc == -1 || rhoQc == +1)
228 return (1 + rhoQc*_acp*(1 - 2*_wQ))*(1 + 0.5*_tag*(2*_delW));
229 else
230 return 1;
231 }
232
233 if (basisIndex == _basisSin) {
234
235 if (rhoQc == -1)
236
237 return - ((1 - _acp)*a_sin_m*(1 - _wQ) + (1 + _acp)*a_sin_p*_wQ)*(1 - 2*_avgW)*_tag;
238
239 else if (rhoQc == +1)
240
241 return - ((1 + _acp)*a_sin_p*(1 - _wQ) + (1 - _acp)*a_sin_m*_wQ)*(1 - 2*_avgW)*_tag;
242
243 else
244 return - _tag*((a_sin_p + a_sin_m)/2)*(1 - 2*_avgW);
245 }
246
247 if (basisIndex == _basisCos) {
248
249 if ( rhoQc == -1)
250
251 return + ((1 - _acp)*a_cos_m*(1 - _wQ) + (1 + _acp)*a_cos_p*_wQ)*(1 - 2*_avgW)*_tag;
252
253 else if (rhoQc == +1)
254
255 return + ((1 + _acp)*a_cos_p*(1 - _wQ) + (1 - _acp)*a_cos_m*_wQ)*(1 - 2*_avgW)*_tag;
256
257 else
258 return _tag*((a_cos_p + a_cos_m)/2)*(1 - 2*_avgW);
259 }
260
261 return 0;
262}
263
264// advertise analytical integration
265
266////////////////////////////////////////////////////////////////////////////////
267
269 RooArgSet& analVars, const char* rangeName ) const
270{
271 if (rangeName) return 0 ;
272
273 if (matchArgs( allVars, analVars, _tag, _rhoQ )) return 3;
274 if (matchArgs( allVars, analVars, _rhoQ )) return 2;
275 if (matchArgs( allVars, analVars, _tag )) return 1;
276
277 return 0;
278}
279
280////////////////////////////////////////////////////////////////////////////////
281/// correct for the right/wrong charge...
282
284 Int_t code, const char* /*rangeName*/ ) const
285{
286 Int_t rhoQc = _rhoQ*int(_correctQ);
287
288 double a_sin_p = _avgS + _delS;
289 double a_sin_m = _avgS - _delS;
290 double a_cos_p = _avgC + _delC;
291 double a_cos_m = _avgC - _delC;
292
293 switch(code) {
294
295 // No integration
296 case 0: return coefficient(basisIndex);
297
298 // Integration over 'tag'
299 case 1:
300 if (basisIndex == _basisExp) return 2*(1 + rhoQc*_acp*(1 - 2*_wQ));
301 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
302 assert( false );
303
304 // Integration over 'rhoQ'
305 case 2:
306 if (basisIndex == _basisExp) return 2*(1 + 0.5*_tag*(2.*_delW));
307
308 if (basisIndex == _basisSin)
309
310 return - ( (1 - _acp)*a_sin_m + (1 + _acp)*a_sin_p )*(1 - 2*_avgW)*_tag;
311
312 if (basisIndex == _basisCos)
313
314 return + ( (1 - _acp)*a_cos_m + (1 + _acp)*a_cos_p )*(1 - 2*_avgW)*_tag;
315
316 assert( false );
317
318 // Integration over 'tag' and 'rhoQ'
319 case 3:
320 if (basisIndex == _basisExp) return 2*2; // for both: tag and charge
321 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
322 assert( false );
323
324 default:
325 assert( false );
326 }
327
328 return 0;
329}
330
331////////////////////////////////////////////////////////////////////////////////
332
334 RooArgSet& generateVars, bool staticInitOK ) const
335{
336 if (staticInitOK) {
337 if (matchArgs( directVars, generateVars, _t, _tag, _rhoQ )) return 4;
338 if (matchArgs( directVars, generateVars, _t, _rhoQ )) return 3;
339 if (matchArgs( directVars, generateVars, _t, _tag )) return 2;
340 }
341 if (matchArgs( directVars, generateVars, _t )) return 1;
342 return 0;
343}
344
345////////////////////////////////////////////////////////////////////////////////
346
348{
349 if (code == 2 || code == 4) {
350 // Calculate the fraction of mixed events to generate
351 double sumInt1 = RooRealIntegral( "sumInt1", "sum integral1", *this,
352 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
353 ).getVal();
354 _tag = -1;
355 double b0Int1 = RooRealIntegral( "mixInt1", "mix integral1", *this,
356 RooArgSet( _t.arg(), _rhoQ.arg() )
357 ).getVal();
358 _genB0Frac = b0Int1/sumInt1;
359
361 cout << " o RooNonCPEigenDecay::initgenerator: genB0Frac : "
362 << _genB0Frac
363 << ", tag dilution: " << (1 - 2*_avgW)
364 << endl;
365 }
366
367 if (code == 3 || code == 4) {
368 // Calculate the fraction of positive rho's to generate
369 double sumInt2 = RooRealIntegral( "sumInt2", "sum integral2", *this,
370 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
371 ).getVal();
372 _rhoQ = 1;
373 double b0Int2 = RooRealIntegral( "mixInt2", "mix integral2", *this,
374 RooArgSet( _t.arg(), _tag.arg() )
375 ).getVal();
376 _genRhoPlusFrac = b0Int2/sumInt2;
377
379 cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: "
380 << _genRhoPlusFrac << endl;
381 }
382}
383
384////////////////////////////////////////////////////////////////////////////////
385
387{
388 // Generate delta-t dependent
389 while (true) {
390
391 // B flavor and rho charge (we do not use the integrated weights)
392 if (code != 1) {
393 if (code != 3) _tag = (RooRandom::uniform()<=0.5) ? -1 : +1;
394 if (code != 2) _rhoQ = (RooRandom::uniform()<=0.5) ? 1 : -1;
395 }
396
397 // opposite charge?
398 // Int_t rhoQc = _rhoQ*int(_correctQ);
399
400 double a_sin_p = _avgS + _delS;
401 double a_sin_m = _avgS - _delS;
402 double a_cos_p = _avgC + _delC;
403 double a_cos_m = _avgC - _delC;
404
405 // maximum probability density
406 double a1 = 1 + sqrt(TMath::Power(a_cos_m, 2) + TMath::Power(a_sin_m, 2));
407 double a2 = 1 + sqrt(TMath::Power(a_cos_p, 2) + TMath::Power(a_sin_p, 2));
408
409 double maxAcceptProb = (1.10 + TMath::Abs(_acp)) * (a1 > a2 ? a1 : a2);
410 // The 1.10 in the above line is a security feature to prevent crashes close to the limit at 1.00
411
412 double rand = RooRandom::uniform();
413 double tval(0);
414
415 switch(_type) {
416
417 case SingleSided:
418 tval = -_tau*log(rand);
419 break;
420
421 case Flipped:
422 tval = +_tau*log(rand);
423 break;
424
425 case DoubleSided:
426 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5));
427 break;
428 }
429
430 // get coefficients
431 double expC = coefficient( _basisExp );
432 double sinC = coefficient( _basisSin );
433 double cosC = coefficient( _basisCos );
434
435 // probability density
436 double acceptProb = expC + sinC*sin(_dm*tval) + cosC*cos(_dm*tval);
437
438 // sanity check...
439 assert( acceptProb <= maxAcceptProb );
440
441 // hit or miss...
442 bool accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? true : false;
443
444 if (accept && tval<_t.max() && tval>_t.min()) {
445 _t = tval;
446 break;
447 }
448 }
449}
#define Debug_RooNonCPEigenDecay
#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 RooRealIntegral
Definition: RooAbsArg.h:629
A space to attach TBranches.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
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:56
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=nullptr) const override
Default implementation of function advertising integration capabilities.
~RooNonCPEigenDecay(void) override
Destructor.
double coefAnalyticalIntegral(Int_t coef, Int_t code, const char *rangeName=nullptr) const override
correct for the right/wrong charge...
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.
RooRealProxy _wQ
dummy mischarge (must be set to zero!)
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 max(const char *rname=nullptr) 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.
double min(const char *rname=nullptr) const
Query lower limit of range. This requires the payload to be RooAbsRealLValue or derived.
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition: RVec.hxx:1776
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1765
RVec< PromoteType< T > > sin(const RVec< T > &v)
Definition: RVec.hxx:1775
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:240
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123