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
53
54#define Debug_RooNonCPEigenDecay 1
55
56
57////////////////////////////////////////////////////////////////////////////////
58
59RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
60 RooRealVar& t,
61 RooAbsCategory& tag,
62 RooAbsReal& tau,
63 RooAbsReal& dm,
70 RooAbsReal& C,
72 RooAbsReal& S,
74 const RooResolutionModel& model,
76 : RooAbsAnaConvPdf( name, title, model, t ),
77 _acp ( "acp", "acp", this, acp ),
78 _avgC ( "C", "C", this, C ),
79 _delC ( "delC", "delC", this, delC ),
80 _avgS ( "S", "S", this, S ),
81 _delS ( "delS", "delS", this, delS ),
82 _avgW ( "avgW", "Average mistag rate",this, avgW ),
83 _delW ( "delW", "Shift mistag rate", this, delW ),
84 _t ( "t", "time", this, t ),
85 _tau ( "tau", "decay time", this, tau ),
86 _dm ( "dm", "mixing frequency", this, dm ),
87 _tag ( "tag", "CP state", this, tag ),
88 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
89 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
90 _wQ ( "wQ", "mischarge", this, wQ ),
91 _genB0Frac ( 0 ),
92 _genRhoPlusFrac( 0 ),
93 _type ( type )
94{
95 // Constructor
96 switch(type) {
97 case SingleSided:
98 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
99 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
100 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
101 break;
102 case Flipped:
103 _basisExp = declareBasis( "exp(@0)/@1)", RooArgList( tau ) );
104 _basisSin = declareBasis( "exp(@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
105 _basisCos = declareBasis( "exp(@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
106 break;
107 case DoubleSided:
108 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
109 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
110 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
111 break;
112 }
113}
114
115////////////////////////////////////////////////////////////////////////////////
116
117RooNonCPEigenDecay::RooNonCPEigenDecay( const char *name, const char *title,
118 RooRealVar& t,
119 RooAbsCategory& tag,
120 RooAbsReal& tau,
121 RooAbsReal& dm,
127 RooAbsReal& C,
129 RooAbsReal& S,
131 const RooResolutionModel& model,
133 : RooAbsAnaConvPdf( name, title, model, t ),
134 _acp ( "acp", "acp", this, acp ),
135 _avgC ( "C", "C", this, C ),
136 _delC ( "delC", "delC", this, delC ),
137 _avgS ( "S", "S", this, S ),
138 _delS ( "delS", "delS", this, delS ),
139 _avgW ( "avgW", "Average mistag rate",this, avgW ),
140 _delW ( "delW", "Shift mistag rate", this, delW ),
141 _t ( "t", "time", this, t ),
142 _tau ( "tau", "decay time", this, tau ),
143 _dm ( "dm", "mixing frequency", this, dm ),
144 _tag ( "tag", "CP state", this, tag ),
145 _rhoQ ( "rhoQ", "Charge of the rho", this, rhoQ ),
146 _correctQ ( "correctQ", "correction of rhoQ", this, correctQ ),
147 _wQ ( "wQ", "mischarge", this, std::make_unique<RooRealVar>( "wQ", "wQ", 0 ), true, false),
148 _genB0Frac ( 0 ),
149 _genRhoPlusFrac( 0 ),
150 _type ( type )
151{
152 switch(type) {
153 case SingleSided:
154 _basisExp = declareBasis( "exp(-@0/@1)", RooArgList( tau ) );
155 _basisSin = declareBasis( "exp(-@0/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
156 _basisCos = declareBasis( "exp(-@0/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
157 break;
158 case Flipped:
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 DoubleSided:
164 _basisExp = declareBasis( "exp(-abs(@0)/@1)", RooArgList( tau ) );
165 _basisSin = declareBasis( "exp(-abs(@0)/@1)*sin(@0*@2)", RooArgList( tau, dm ) );
166 _basisCos = declareBasis( "exp(-abs(@0)/@1)*cos(@0*@2)", RooArgList( tau, dm ) );
167 break;
168 }
169}
170
171////////////////////////////////////////////////////////////////////////////////
172/// Copy constructor
173
176 _acp ( "acp", this, other._acp ),
177 _avgC ( "C", this, other._avgC ),
178 _delC ( "delC", this, other._delC ),
179 _avgS ( "S", this, other._avgS ),
180 _delS ( "delS", this, other._delS ),
181 _avgW ( "avgW", this, other._avgW ),
182 _delW ( "delW", this, other._delW ),
183 _t ( "t", this, other._t ),
184 _tau ( "tau", this, other._tau ),
185 _dm ( "dm", this, other._dm ),
186 _tag ( "tag", this, other._tag ),
187 _rhoQ ( "rhoQ", this, other._rhoQ ),
188 _correctQ ( "correctQ", this, other._correctQ ),
189 _wQ ( "wQ", this, other._wQ ),
190 _genB0Frac ( other._genB0Frac ),
191 _genRhoPlusFrac( other._genRhoPlusFrac ),
192 _type ( other._type ),
193 _basisExp ( other._basisExp ),
194 _basisSin ( other._basisSin ),
195 _basisCos ( other._basisCos )
196{
197}
198
199////////////////////////////////////////////////////////////////////////////////
200/// - B0 : _tag == -1
201/// - B0bar : _tag == +1
202/// - rho+ : _rhoQ == +1
203/// - rho- : _rhoQ == -1
204/// - the charge correction factor "_correctQ" serves to implement misidentified charges
205
207{
209 assert( rhoQc == 1 || rhoQc == -1 );
210
211 double a_sin_p = _avgS + _delS;
212 double a_sin_m = _avgS - _delS;
213 double a_cos_p = _avgC + _delC;
214 double a_cos_m = _avgC - _delC;
215
216 if (basisIndex == _basisExp) {
217 if (rhoQc == -1 || rhoQc == +1) {
218 return (1 + rhoQc * _acp * (1 - 2 * _wQ)) * (1 + 0.5 * _tag * (2 * _delW));
219 } else {
220 return 1;
221 }
222 }
223
224 if (basisIndex == _basisSin) {
225
226 if (rhoQc == -1) {
227
228 return -((1 - _acp) * a_sin_m * (1 - _wQ) + (1 + _acp) * a_sin_p * _wQ) * (1 - 2 * _avgW) * _tag;
229
230 } else if (rhoQc == +1) {
231
232 return -((1 + _acp) * a_sin_p * (1 - _wQ) + (1 - _acp) * a_sin_m * _wQ) * (1 - 2 * _avgW) * _tag;
233
234 } else {
235 return -_tag * ((a_sin_p + a_sin_m) / 2) * (1 - 2 * _avgW);
236 }
237 }
238
239 if (basisIndex == _basisCos) {
240
241 if (rhoQc == -1) {
242
243 return +((1 - _acp) * a_cos_m * (1 - _wQ) + (1 + _acp) * a_cos_p * _wQ) * (1 - 2 * _avgW) * _tag;
244
245 } else if (rhoQc == +1) {
246
247 return +((1 + _acp) * a_cos_p * (1 - _wQ) + (1 - _acp) * a_cos_m * _wQ) * (1 - 2 * _avgW) * _tag;
248
249 } else {
250 return _tag * ((a_cos_p + a_cos_m) / 2) * (1 - 2 * _avgW);
251 }
252 }
253
254 return 0;
255}
256
257// advertise analytical integration
258
259////////////////////////////////////////////////////////////////////////////////
260
262 RooArgSet& analVars, const char* rangeName ) const
263{
264 if (rangeName) return 0 ;
265
266 if (matchArgs( allVars, analVars, _tag, _rhoQ )) return 3;
267 if (matchArgs( allVars, analVars, _rhoQ )) return 2;
268 if (matchArgs( allVars, analVars, _tag )) return 1;
269
270 return 0;
271}
272
273////////////////////////////////////////////////////////////////////////////////
274/// correct for the right/wrong charge...
275
277 Int_t code, const char* /*rangeName*/ ) const
278{
280
281 double a_sin_p = _avgS + _delS;
282 double a_sin_m = _avgS - _delS;
283 double a_cos_p = _avgC + _delC;
284 double a_cos_m = _avgC - _delC;
285
286 switch(code) {
287
288 // No integration
289 case 0: return coefficient(basisIndex);
290
291 // Integration over 'tag'
292 case 1:
293 if (basisIndex == _basisExp) return 2*(1 + rhoQc*_acp*(1 - 2*_wQ));
294 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
295 assert( false );
296
297 // Integration over 'rhoQ'
298 case 2:
299 if (basisIndex == _basisExp) return 2*(1 + 0.5*_tag*(2.*_delW));
300
301 if (basisIndex == _basisSin) {
302
303 return -((1 - _acp) * a_sin_m + (1 + _acp) * a_sin_p) * (1 - 2 * _avgW) * _tag;
304 }
305
306 if (basisIndex == _basisCos) {
307
308 return +((1 - _acp) * a_cos_m + (1 + _acp) * a_cos_p) * (1 - 2 * _avgW) * _tag;
309 }
310
311 assert( false );
312
313 // Integration over 'tag' and 'rhoQ'
314 case 3:
315 if (basisIndex == _basisExp) return 2*2; // for both: tag and charge
316 if (basisIndex == _basisSin || basisIndex==_basisCos) return 0;
317 assert( false );
318
319 default:
320 assert( false );
321 }
322
323 return 0;
324}
325
326////////////////////////////////////////////////////////////////////////////////
327
329 RooArgSet& generateVars, bool staticInitOK ) const
330{
331 if (staticInitOK) {
332 if (matchArgs( directVars, generateVars, _t, _tag, _rhoQ )) return 4;
333 if (matchArgs( directVars, generateVars, _t, _rhoQ )) return 3;
334 if (matchArgs( directVars, generateVars, _t, _tag )) return 2;
335 }
336 if (matchArgs( directVars, generateVars, _t )) return 1;
337 return 0;
338}
339
340////////////////////////////////////////////////////////////////////////////////
341
343{
344 if (code == 2 || code == 4) {
345 // Calculate the fraction of mixed events to generate
346 double sumInt1 = RooRealIntegral( "sumInt1", "sum integral1", *this,
347 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
348 ).getVal();
349 _tag = -1;
350 double b0Int1 = RooRealIntegral( "mixInt1", "mix integral1", *this,
351 RooArgSet( _t.arg(), _rhoQ.arg() )
352 ).getVal();
354
355 if (Debug_RooNonCPEigenDecay == 1) {
356 std::cout << " o RooNonCPEigenDecay::initgenerator: genB0Frac : " << _genB0Frac
357 << ", tag dilution: " << (1 - 2 * _avgW) << std::endl;
358 }
359 }
360
361 if (code == 3 || code == 4) {
362 // Calculate the fraction of positive rho's to generate
363 double sumInt2 = RooRealIntegral( "sumInt2", "sum integral2", *this,
364 RooArgSet( _t.arg(), _tag.arg(), _rhoQ.arg() )
365 ).getVal();
366 _rhoQ = 1;
367 double b0Int2 = RooRealIntegral( "mixInt2", "mix integral2", *this,
368 RooArgSet( _t.arg(), _tag.arg() )
369 ).getVal();
371
372 if (Debug_RooNonCPEigenDecay == 1) {
373 std::cout << " o RooNonCPEigenDecay::initgenerator: genRhoPlusFrac: " << _genRhoPlusFrac << std::endl;
374 }
375 }
376}
377
378////////////////////////////////////////////////////////////////////////////////
379
381{
382 // Generate delta-t dependent
383 while (true) {
384
385 // B flavor and rho charge (we do not use the integrated weights)
386 if (code != 1) {
387 if (code != 3) _tag = (RooRandom::uniform()<=0.5) ? -1 : +1;
388 if (code != 2) _rhoQ = (RooRandom::uniform()<=0.5) ? 1 : -1;
389 }
390
391 // opposite charge?
392 // Int_t rhoQc = _rhoQ*int(_correctQ);
393
394 double a_sin_p = _avgS + _delS;
395 double a_sin_m = _avgS - _delS;
396 double a_cos_p = _avgC + _delC;
397 double a_cos_m = _avgC - _delC;
398
399 // maximum probability density
400 double a1 = 1 + sqrt(std::pow(a_cos_m, 2) + std::pow(a_sin_m, 2));
401 double a2 = 1 + sqrt(std::pow(a_cos_p, 2) + std::pow(a_sin_p, 2));
402
403 double maxAcceptProb = (1.10 + std::abs(_acp)) * (a1 > a2 ? a1 : a2);
404 // The 1.10 in the above line is a security feature to prevent crashes close to the limit at 1.00
405
406 double rand = RooRandom::uniform();
407 double tval(0);
408
409 switch(_type) {
410
411 case SingleSided:
412 tval = -_tau*log(rand);
413 break;
414
415 case Flipped:
416 tval = +_tau*log(rand);
417 break;
418
419 case DoubleSided:
420 tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5));
421 break;
422 }
423
424 // get coefficients
425 double expC = coefficient( _basisExp );
426 double sinC = coefficient( _basisSin );
427 double cosC = coefficient( _basisCos );
428
429 // probability density
430 double acceptProb = expC + sinC*sin(_dm*tval) + cosC*cos(_dm*tval);
431
432 // sanity check...
434
435 // hit or miss...
437
438 if (accept && tval<_t.max() && tval>_t.min()) {
439 _t = tval;
440 break;
441 }
442 }
443}
#define Debug_RooNonCPEigenDecay
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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.
friend class RooRealIntegral
Definition RooAbsArg.h:572
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
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:24
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:77
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.