Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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,
73 RooAbsReal& C,
74 RooAbsReal& delC,
75 RooAbsReal& S,
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/// - B0 : _tag == -1
204/// - B0bar : _tag == +1
205/// - rho+ : _rhoQ == +1
206/// - rho- : _rhoQ == -1
207/// - the charge correction factor "_correctQ" serves to implement misidentified charges
208
209double RooNonCPEigenDecay::coefficient( Int_t basisIndex ) const
210{
211 Int_t rhoQc = _rhoQ * int(_correctQ);
212 assert( rhoQc == 1 || rhoQc == -1 );
213
214 double a_sin_p = _avgS + _delS;
215 double a_sin_m = _avgS - _delS;
216 double a_cos_p = _avgC + _delC;
217 double a_cos_m = _avgC - _delC;
218
219 if (basisIndex == _basisExp) {
220 if (rhoQc == -1 || rhoQc == +1) {
221 return (1 + rhoQc * _acp * (1 - 2 * _wQ)) * (1 + 0.5 * _tag * (2 * _delW));
222 } else {
223 return 1;
224 }
225 }
226
227 if (basisIndex == _basisSin) {
228
229 if (rhoQc == -1) {
230
231 return -((1 - _acp) * a_sin_m * (1 - _wQ) + (1 + _acp) * a_sin_p * _wQ) * (1 - 2 * _avgW) * _tag;
232
233 } else if (rhoQc == +1) {
234
235 return -((1 + _acp) * a_sin_p * (1 - _wQ) + (1 - _acp) * a_sin_m * _wQ) * (1 - 2 * _avgW) * _tag;
236
237 } else {
238 return -_tag * ((a_sin_p + a_sin_m) / 2) * (1 - 2 * _avgW);
239 }
240 }
241
242 if (basisIndex == _basisCos) {
243
244 if (rhoQc == -1) {
245
246 return +((1 - _acp) * a_cos_m * (1 - _wQ) + (1 + _acp) * a_cos_p * _wQ) * (1 - 2 * _avgW) * _tag;
247
248 } else if (rhoQc == +1) {
249
250 return +((1 + _acp) * a_cos_p * (1 - _wQ) + (1 - _acp) * a_cos_m * _wQ) * (1 - 2 * _avgW) * _tag;
251
252 } else {
253 return _tag * ((a_cos_p + a_cos_m) / 2) * (1 - 2 * _avgW);
254 }
255 }
256
257 return 0;
258}
259
260// advertise analytical integration
261
262////////////////////////////////////////////////////////////////////////////////
263
265 RooArgSet& analVars, const char* rangeName ) const
266{
267 if (rangeName) return 0 ;
268
269 if (matchArgs( allVars, analVars, _tag, _rhoQ )) return 3;
270 if (matchArgs( allVars, analVars, _rhoQ )) return 2;
271 if (matchArgs( allVars, analVars, _tag )) return 1;
272
273 return 0;
274}
275
276////////////////////////////////////////////////////////////////////////////////
277/// correct for the right/wrong charge...
278
280 Int_t code, const char* /*rangeName*/ ) const
281{
282 Int_t rhoQc = _rhoQ*int(_correctQ);
283
284 double a_sin_p = _avgS + _delS;
285 double a_sin_m = _avgS - _delS;
286 double a_cos_p = _avgC + _delC;
287 double a_cos_m = _avgC - _delC;
288
289 switch(code) {
290
291 // No integration
292 case 0: return coefficient(basisIndex);
293
294 // Integration over 'tag'
295 case 1:
296 if (basisIndex == _basisExp) return 2*(1 + rhoQc*_acp*(1 - 2*_wQ));
297 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
298 assert( false );
299
300 // Integration over 'rhoQ'
301 case 2:
302 if (basisIndex == _basisExp) return 2*(1 + 0.5*_tag*(2.*_delW));
303
304 if (basisIndex == _basisSin) {
305
306 return -((1 - _acp) * a_sin_m + (1 + _acp) * a_sin_p) * (1 - 2 * _avgW) * _tag;
307 }
308
309 if (basisIndex == _basisCos) {
310
311 return +((1 - _acp) * a_cos_m + (1 + _acp) * a_cos_p) * (1 - 2 * _avgW) * _tag;
312 }
313
314 assert( false );
315
316 // Integration over 'tag' and 'rhoQ'
317 case 3:
318 if (basisIndex == _basisExp) return 2*2; // for both: tag and charge
319 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
320 assert( false );
321
322 default:
323 assert( false );
324 }
325
326 return 0;
327}
328
329////////////////////////////////////////////////////////////////////////////////
330
332 RooArgSet& generateVars, bool staticInitOK ) const
333{
334 if (staticInitOK) {
335 if (matchArgs( directVars, generateVars, _t, _tag, _rhoQ )) return 4;
336 if (matchArgs( directVars, generateVars, _t, _rhoQ )) return 3;
337 if (matchArgs( directVars, generateVars, _t, _tag )) return 2;
338 }
339 if (matchArgs( directVars, generateVars, _t )) return 1;
340 return 0;
341}
342
343////////////////////////////////////////////////////////////////////////////////
344
346{
347 if (code == 2 || code == 4) {
348 // Calculate the fraction of mixed events to generate
349 double sumInt1 = RooRealIntegral( "sumInt1", "sum integral1", *this,
350 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
351 ).getVal();
352 _tag = -1;
353 double b0Int1 = RooRealIntegral( "mixInt1", "mix integral1", *this,
354 RooArgSet( _t.arg(), _rhoQ.arg() )
355 ).getVal();
356 _genB0Frac = b0Int1/sumInt1;
357
358 if (Debug_RooNonCPEigenDecay == 1) {
359 cout << " o RooNonCPEigenDecay::initgenerator: genB0Frac : " << _genB0Frac
360 << ", tag dilution: " << (1 - 2 * _avgW) << endl;
361 }
362 }
363
364 if (code == 3 || code == 4) {
365 // Calculate the fraction of positive rho's to generate
366 double sumInt2 = RooRealIntegral( "sumInt2", "sum integral2", *this,
367 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
368 ).getVal();
369 _rhoQ = 1;
370 double b0Int2 = RooRealIntegral( "mixInt2", "mix integral2", *this,
371 RooArgSet( _t.arg(), _tag.arg() )
372 ).getVal();
373 _genRhoPlusFrac = b0Int2/sumInt2;
374
375 if (Debug_RooNonCPEigenDecay == 1) {
376 cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: " << _genRhoPlusFrac << endl;
377 }
378 }
379}
380
381////////////////////////////////////////////////////////////////////////////////
382
384{
385 // Generate delta-t dependent
386 while (true) {
387
388 // B flavor and rho charge (we do not use the integrated weights)
389 if (code != 1) {
390 if (code != 3) _tag = (RooRandom::uniform()<=0.5) ? -1 : +1;
391 if (code != 2) _rhoQ = (RooRandom::uniform()<=0.5) ? 1 : -1;
392 }
393
394 // opposite charge?
395 // Int_t rhoQc = _rhoQ*int(_correctQ);
396
397 double a_sin_p = _avgS + _delS;
398 double a_sin_m = _avgS - _delS;
399 double a_cos_p = _avgC + _delC;
400 double a_cos_m = _avgC - _delC;
401
402 // maximum probability density
403 double a1 = 1 + sqrt(TMath::Power(a_cos_m, 2) + TMath::Power(a_sin_m, 2));
404 double a2 = 1 + sqrt(TMath::Power(a_cos_p, 2) + TMath::Power(a_sin_p, 2));
405
406 double maxAcceptProb = (1.10 + TMath::Abs(_acp)) * (a1 > a2 ? a1 : a2);
407 // The 1.10 in the above line is a security feature to prevent crashes close to the limit at 1.00
408
409 double rand = RooRandom::uniform();
410 double tval(0);
411
412 switch(_type) {
413
414 case SingleSided:
415 tval = -_tau*log(rand);
416 break;
417
418 case Flipped:
419 tval = +_tau*log(rand);
420 break;
421
422 case DoubleSided:
423 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5));
424 break;
425 }
426
427 // get coefficients
428 double expC = coefficient( _basisExp );
429 double sinC = coefficient( _basisSin );
430 double cosC = coefficient( _basisCos );
431
432 // probability density
433 double acceptProb = expC + sinC*sin(_dm*tval) + cosC*cos(_dm*tval);
434
435 // sanity check...
436 assert( acceptProb <= maxAcceptProb );
437
438 // hit or miss...
439 bool accept = maxAcceptProb*RooRandom::uniform() < acceptProb ? true : false;
440
441 if (accept && tval<_t.max() && tval>_t.min()) {
442 _t = tval;
443 break;
444 }
445 }
446}
#define Debug_RooNonCPEigenDecay
#define ClassImp(name)
Definition Rtypes.h:377
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
Base class for PDFs that represent a physics model that can be analytically convolved with a resoluti...
Int_t declareBasis(const char *expression, const RooArgList &params)
Declare a basis function for use in this physics model.
A space to attach TBranches.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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:55
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.
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
Performs hybrid numerical/analytical integrals of RooAbsReal objects.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
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.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123