ROOT logo
/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitModels                                                     *
 * @(#)root/roofit:$Id: RooBMixDecay.cxx 24286 2008-06-16 15:47:04Z wouter $
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// Class RooBMixDecay is a RooAbsAnaConvPdf implementation that describes
// the decay of B mesons with the effects of B0/B0bar mixing. 
// This function can be analytically convolved with any RooResolutionModel implementation
// END_HTML
//

#include "RooFit.h"

#include "Riostream.h"
#include "TMath.h"
#include "RooRealVar.h"
#include "RooBMixDecay.h"
#include "RooRandom.h"

ClassImp(RooBMixDecay) 
;



//_____________________________________________________________________________
RooBMixDecay::RooBMixDecay(const char *name, const char *title, 
			   RooRealVar& t, RooAbsCategory& mixState,
			   RooAbsCategory& tagFlav,
			   RooAbsReal& tau, RooAbsReal& dm,			   
			   RooAbsReal& mistag, RooAbsReal& delMistag,
			   const RooResolutionModel& model, 
			   DecayType type) :
  RooAbsAnaConvPdf(name,title,model,t), 
  _type(type),
  _mistag("mistag","Mistag rate",this,mistag),
  _delMistag("delMistag","Delta mistag rate",this,delMistag),
  _mixState("mixState","Mixing state",this,mixState),
  _tagFlav("tagFlav","Flavour of tagged B0",this,tagFlav),
  _tau("tau","Mixing life time",this,tau),
  _dm("dm","Mixing frequency",this,dm),
  _t("_t","time",this,t), _genMixFrac(0)
{
  // Constructor
  switch(type) {
  case SingleSided:
    _basisExp = declareBasis("exp(-@0/@1)",RooArgList(tau,dm)) ;
    _basisCos = declareBasis("exp(-@0/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
    break ;
  case Flipped:
    _basisExp = declareBasis("exp(@0/@1)",RooArgList(tau,dm)) ;
    _basisCos = declareBasis("exp(@0/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
    break ;
  case DoubleSided:
    _basisExp = declareBasis("exp(-abs(@0)/@1)",RooArgList(tau,dm)) ;
    _basisCos = declareBasis("exp(-abs(@0)/@1)*cos(@0*@2)",RooArgList(tau,dm)) ;
    break ;
  }
}



//_____________________________________________________________________________
RooBMixDecay::RooBMixDecay(const RooBMixDecay& other, const char* name) : 
  RooAbsAnaConvPdf(other,name), 
  _type(other._type),
  _mistag("mistag",this,other._mistag),
  _delMistag("delMistag",this,other._delMistag),
  _mixState("mixState",this,other._mixState),
  _tagFlav("tagFlav",this,other._tagFlav),
  _tau("tau",this,other._tau),
  _dm("dm",this,other._dm),
  _t("t",this,other._t),
  _basisExp(other._basisExp),
  _basisCos(other._basisCos),
  _genMixFrac(other._genMixFrac),
  _genFlavFrac(other._genFlavFrac),
  _genFlavFracMix(other._genFlavFracMix),
  _genFlavFracUnmix(other._genFlavFracUnmix)
{
  // Copy constructor
}



//_____________________________________________________________________________
RooBMixDecay::~RooBMixDecay()
{
  // Destructor
}



//_____________________________________________________________________________
Double_t RooBMixDecay::coefficient(Int_t basisIndex) const 
{
  // Comp with tFit MC: must be (1 - tagFlav*...)
  if (basisIndex==_basisExp) {
    return (1 - _tagFlav*_delMistag) ; 
  }

  if (basisIndex==_basisCos) {
    return _mixState*(1-2*_mistag) ;   
  }
  
  return 0 ;
}



//_____________________________________________________________________________
Int_t RooBMixDecay::getCoefAnalyticalIntegral(Int_t /*code*/, RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const 
{
//   cout << "RooBMixDecay::getCoefAI " ; allVars.Print("1") ;
  if (rangeName) {
    return 0 ;
  }

  if (matchArgs(allVars,analVars,_mixState,_tagFlav)) return 3 ;
  if (matchArgs(allVars,analVars,_mixState)) return 2 ;
  if (matchArgs(allVars,analVars,_tagFlav)) return 1 ;
  return 0 ;
}



//_____________________________________________________________________________
Double_t RooBMixDecay::coefAnalyticalIntegral(Int_t basisIndex, Int_t code, const char* /*rangeName*/) const 
{  
  switch(code) {
    // No integration
  case 0: return coefficient(basisIndex) ;

    // Integration over 'mixState' and 'tagFlav' 
  case 3:
    if (basisIndex==_basisExp) {
      return 4.0 ;
    }    
    if (basisIndex==_basisCos) {
      return 0.0 ;
    }

    // Integration over 'mixState'
  case 2:
    if (basisIndex==_basisExp) {
      return 2.0*coefficient(basisIndex) ;
    }    
    if (basisIndex==_basisCos) {
      return 0.0 ;
    }

    // Integration over 'tagFlav'
  case 1:
    if (basisIndex==_basisExp) {
      return 2.0 ;
    }    
    if (basisIndex==_basisCos) {
      return 2.0*coefficient(basisIndex) ;
    }
  default:
    assert(0) ;
  }
    
  return 0 ;
}



//_____________________________________________________________________________
Int_t RooBMixDecay::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK) const
{
  if (staticInitOK) {
    if (matchArgs(directVars,generateVars,_t,_mixState,_tagFlav)) return 4 ;  
    if (matchArgs(directVars,generateVars,_t,_mixState)) return 3 ;  
    if (matchArgs(directVars,generateVars,_t,_tagFlav)) return 2 ;  
  }

  if (matchArgs(directVars,generateVars,_t)) return 1 ;  
  return 0 ;
}



//_____________________________________________________________________________
void RooBMixDecay::initGenerator(Int_t code)
{
  switch (code) {
  case 2:
    {
      // Calculate the fraction of B0bar events to generate
      Double_t sumInt = RooRealIntegral("sumInt","sum integral",*this,RooArgSet(_t.arg(),_tagFlav.arg())).getVal() ;
      _tagFlav = 1 ; // B0 
      Double_t flavInt = RooRealIntegral("flavInt","flav integral",*this,RooArgSet(_t.arg())).getVal() ;
      _genFlavFrac = flavInt/sumInt ;
      break ;
    }  
  case 3:
    {
      // Calculate the fraction of mixed events to generate
      Double_t sumInt = RooRealIntegral("sumInt","sum integral",*this,RooArgSet(_t.arg(),_mixState.arg())).getVal() ;
      _mixState = -1 ; // mixed
      Double_t mixInt = RooRealIntegral("mixInt","mix integral",*this,RooArgSet(_t.arg())).getVal() ;
      _genMixFrac = mixInt/sumInt ;
      break ;
    }  
  case 4:
    {
      // Calculate the fraction of mixed events to generate
      Double_t sumInt = RooRealIntegral("sumInt","sum integral",*this,RooArgSet(_t.arg(),_mixState.arg(),_tagFlav.arg())).getVal() ;
      _mixState = -1 ; // mixed
      Double_t mixInt = RooRealIntegral("mixInt","mix integral",*this,RooArgSet(_t.arg(),_tagFlav.arg())).getVal() ;
      _genMixFrac = mixInt/sumInt ;
      
      // Calculate the fractio of B0bar tags for mixed and unmixed
      RooRealIntegral dtInt("mixInt","mix integral",*this,RooArgSet(_t.arg())) ;
      _mixState = -1 ; // Mixed
      _tagFlav  =  1 ; // B0
      _genFlavFracMix   = dtInt.getVal() / mixInt ;
      _mixState =  1 ; // Unmixed
      _tagFlav  =  1 ; // B0
      _genFlavFracUnmix = dtInt.getVal() / (sumInt - mixInt) ;
      break ;
    }
  }
}




//_____________________________________________________________________________
void RooBMixDecay::generateEvent(Int_t code)
{
  // Generate mix-state dependent
  switch(code) {
  case 2:
    {
      Double_t rand = RooRandom::uniform() ;
      _tagFlav = (Int_t) ((rand<=_genFlavFrac) ?  1 : -1) ;
      break ;
    }
  case 3:
    {
      Double_t rand = RooRandom::uniform() ;
      _mixState = (Int_t) ((rand<=_genMixFrac) ? -1 : 1) ;
      break ;
    }
  case 4:
    {
      Double_t rand = RooRandom::uniform() ;
      _mixState = (Int_t) ((rand<=_genMixFrac) ? -1 : 1) ;

      rand = RooRandom::uniform() ;
      Double_t genFlavFrac = (_mixState==-1) ? _genFlavFracMix : _genFlavFracUnmix ;
      _tagFlav = (Int_t) ((rand<=genFlavFrac) ?  1 : -1) ;
      break ;
    }
  }

  // Generate delta-t dependent
  while(1) {
    Double_t rand = RooRandom::uniform() ;
    Double_t tval(0) ;

    switch(_type) {
    case SingleSided:
      tval = -_tau*log(rand);
      break ;
    case Flipped:
      tval= +_tau*log(rand);
      break ;
    case DoubleSided:
      tval = (rand<=0.5) ? -_tau*log(2*rand) : +_tau*log(2*(rand-0.5)) ;
      break ;
    }

    // Accept event if T is in generated range
    Double_t dil = 1-2.*_mistag ;
    Double_t maxAcceptProb = 1 + TMath::Abs(_delMistag) + TMath::Abs(dil) ;
    Double_t acceptProb = (1-_tagFlav*_delMistag) + _mixState*dil*cos(_dm*tval);
    Bool_t mixAccept = maxAcceptProb*RooRandom::uniform() < acceptProb ? kTRUE : kFALSE ;
    
    if (tval<_t.max() && tval>_t.min() && mixAccept) {
      _t = tval ;
      break ;
    }
  }  
}
 RooBMixDecay.cxx:1
 RooBMixDecay.cxx:2
 RooBMixDecay.cxx:3
 RooBMixDecay.cxx:4
 RooBMixDecay.cxx:5
 RooBMixDecay.cxx:6
 RooBMixDecay.cxx:7
 RooBMixDecay.cxx:8
 RooBMixDecay.cxx:9
 RooBMixDecay.cxx:10
 RooBMixDecay.cxx:11
 RooBMixDecay.cxx:12
 RooBMixDecay.cxx:13
 RooBMixDecay.cxx:14
 RooBMixDecay.cxx:15
 RooBMixDecay.cxx:16
 RooBMixDecay.cxx:17
 RooBMixDecay.cxx:18
 RooBMixDecay.cxx:19
 RooBMixDecay.cxx:20
 RooBMixDecay.cxx:21
 RooBMixDecay.cxx:22
 RooBMixDecay.cxx:23
 RooBMixDecay.cxx:24
 RooBMixDecay.cxx:25
 RooBMixDecay.cxx:26
 RooBMixDecay.cxx:27
 RooBMixDecay.cxx:28
 RooBMixDecay.cxx:29
 RooBMixDecay.cxx:30
 RooBMixDecay.cxx:31
 RooBMixDecay.cxx:32
 RooBMixDecay.cxx:33
 RooBMixDecay.cxx:34
 RooBMixDecay.cxx:35
 RooBMixDecay.cxx:36
 RooBMixDecay.cxx:37
 RooBMixDecay.cxx:38
 RooBMixDecay.cxx:39
 RooBMixDecay.cxx:40
 RooBMixDecay.cxx:41
 RooBMixDecay.cxx:42
 RooBMixDecay.cxx:43
 RooBMixDecay.cxx:44
 RooBMixDecay.cxx:45
 RooBMixDecay.cxx:46
 RooBMixDecay.cxx:47
 RooBMixDecay.cxx:48
 RooBMixDecay.cxx:49
 RooBMixDecay.cxx:50
 RooBMixDecay.cxx:51
 RooBMixDecay.cxx:52
 RooBMixDecay.cxx:53
 RooBMixDecay.cxx:54
 RooBMixDecay.cxx:55
 RooBMixDecay.cxx:56
 RooBMixDecay.cxx:57
 RooBMixDecay.cxx:58
 RooBMixDecay.cxx:59
 RooBMixDecay.cxx:60
 RooBMixDecay.cxx:61
 RooBMixDecay.cxx:62
 RooBMixDecay.cxx:63
 RooBMixDecay.cxx:64
 RooBMixDecay.cxx:65
 RooBMixDecay.cxx:66
 RooBMixDecay.cxx:67
 RooBMixDecay.cxx:68
 RooBMixDecay.cxx:69
 RooBMixDecay.cxx:70
 RooBMixDecay.cxx:71
 RooBMixDecay.cxx:72
 RooBMixDecay.cxx:73
 RooBMixDecay.cxx:74
 RooBMixDecay.cxx:75
 RooBMixDecay.cxx:76
 RooBMixDecay.cxx:77
 RooBMixDecay.cxx:78
 RooBMixDecay.cxx:79
 RooBMixDecay.cxx:80
 RooBMixDecay.cxx:81
 RooBMixDecay.cxx:82
 RooBMixDecay.cxx:83
 RooBMixDecay.cxx:84
 RooBMixDecay.cxx:85
 RooBMixDecay.cxx:86
 RooBMixDecay.cxx:87
 RooBMixDecay.cxx:88
 RooBMixDecay.cxx:89
 RooBMixDecay.cxx:90
 RooBMixDecay.cxx:91
 RooBMixDecay.cxx:92
 RooBMixDecay.cxx:93
 RooBMixDecay.cxx:94
 RooBMixDecay.cxx:95
 RooBMixDecay.cxx:96
 RooBMixDecay.cxx:97
 RooBMixDecay.cxx:98
 RooBMixDecay.cxx:99
 RooBMixDecay.cxx:100
 RooBMixDecay.cxx:101
 RooBMixDecay.cxx:102
 RooBMixDecay.cxx:103
 RooBMixDecay.cxx:104
 RooBMixDecay.cxx:105
 RooBMixDecay.cxx:106
 RooBMixDecay.cxx:107
 RooBMixDecay.cxx:108
 RooBMixDecay.cxx:109
 RooBMixDecay.cxx:110
 RooBMixDecay.cxx:111
 RooBMixDecay.cxx:112
 RooBMixDecay.cxx:113
 RooBMixDecay.cxx:114
 RooBMixDecay.cxx:115
 RooBMixDecay.cxx:116
 RooBMixDecay.cxx:117
 RooBMixDecay.cxx:118
 RooBMixDecay.cxx:119
 RooBMixDecay.cxx:120
 RooBMixDecay.cxx:121
 RooBMixDecay.cxx:122
 RooBMixDecay.cxx:123
 RooBMixDecay.cxx:124
 RooBMixDecay.cxx:125
 RooBMixDecay.cxx:126
 RooBMixDecay.cxx:127
 RooBMixDecay.cxx:128
 RooBMixDecay.cxx:129
 RooBMixDecay.cxx:130
 RooBMixDecay.cxx:131
 RooBMixDecay.cxx:132
 RooBMixDecay.cxx:133
 RooBMixDecay.cxx:134
 RooBMixDecay.cxx:135
 RooBMixDecay.cxx:136
 RooBMixDecay.cxx:137
 RooBMixDecay.cxx:138
 RooBMixDecay.cxx:139
 RooBMixDecay.cxx:140
 RooBMixDecay.cxx:141
 RooBMixDecay.cxx:142
 RooBMixDecay.cxx:143
 RooBMixDecay.cxx:144
 RooBMixDecay.cxx:145
 RooBMixDecay.cxx:146
 RooBMixDecay.cxx:147
 RooBMixDecay.cxx:148
 RooBMixDecay.cxx:149
 RooBMixDecay.cxx:150
 RooBMixDecay.cxx:151
 RooBMixDecay.cxx:152
 RooBMixDecay.cxx:153
 RooBMixDecay.cxx:154
 RooBMixDecay.cxx:155
 RooBMixDecay.cxx:156
 RooBMixDecay.cxx:157
 RooBMixDecay.cxx:158
 RooBMixDecay.cxx:159
 RooBMixDecay.cxx:160
 RooBMixDecay.cxx:161
 RooBMixDecay.cxx:162
 RooBMixDecay.cxx:163
 RooBMixDecay.cxx:164
 RooBMixDecay.cxx:165
 RooBMixDecay.cxx:166
 RooBMixDecay.cxx:167
 RooBMixDecay.cxx:168
 RooBMixDecay.cxx:169
 RooBMixDecay.cxx:170
 RooBMixDecay.cxx:171
 RooBMixDecay.cxx:172
 RooBMixDecay.cxx:173
 RooBMixDecay.cxx:174
 RooBMixDecay.cxx:175
 RooBMixDecay.cxx:176
 RooBMixDecay.cxx:177
 RooBMixDecay.cxx:178
 RooBMixDecay.cxx:179
 RooBMixDecay.cxx:180
 RooBMixDecay.cxx:181
 RooBMixDecay.cxx:182
 RooBMixDecay.cxx:183
 RooBMixDecay.cxx:184
 RooBMixDecay.cxx:185
 RooBMixDecay.cxx:186
 RooBMixDecay.cxx:187
 RooBMixDecay.cxx:188
 RooBMixDecay.cxx:189
 RooBMixDecay.cxx:190
 RooBMixDecay.cxx:191
 RooBMixDecay.cxx:192
 RooBMixDecay.cxx:193
 RooBMixDecay.cxx:194
 RooBMixDecay.cxx:195
 RooBMixDecay.cxx:196
 RooBMixDecay.cxx:197
 RooBMixDecay.cxx:198
 RooBMixDecay.cxx:199
 RooBMixDecay.cxx:200
 RooBMixDecay.cxx:201
 RooBMixDecay.cxx:202
 RooBMixDecay.cxx:203
 RooBMixDecay.cxx:204
 RooBMixDecay.cxx:205
 RooBMixDecay.cxx:206
 RooBMixDecay.cxx:207
 RooBMixDecay.cxx:208
 RooBMixDecay.cxx:209
 RooBMixDecay.cxx:210
 RooBMixDecay.cxx:211
 RooBMixDecay.cxx:212
 RooBMixDecay.cxx:213
 RooBMixDecay.cxx:214
 RooBMixDecay.cxx:215
 RooBMixDecay.cxx:216
 RooBMixDecay.cxx:217
 RooBMixDecay.cxx:218
 RooBMixDecay.cxx:219
 RooBMixDecay.cxx:220
 RooBMixDecay.cxx:221
 RooBMixDecay.cxx:222
 RooBMixDecay.cxx:223
 RooBMixDecay.cxx:224
 RooBMixDecay.cxx:225
 RooBMixDecay.cxx:226
 RooBMixDecay.cxx:227
 RooBMixDecay.cxx:228
 RooBMixDecay.cxx:229
 RooBMixDecay.cxx:230
 RooBMixDecay.cxx:231
 RooBMixDecay.cxx:232
 RooBMixDecay.cxx:233
 RooBMixDecay.cxx:234
 RooBMixDecay.cxx:235
 RooBMixDecay.cxx:236
 RooBMixDecay.cxx:237
 RooBMixDecay.cxx:238
 RooBMixDecay.cxx:239
 RooBMixDecay.cxx:240
 RooBMixDecay.cxx:241
 RooBMixDecay.cxx:242
 RooBMixDecay.cxx:243
 RooBMixDecay.cxx:244
 RooBMixDecay.cxx:245
 RooBMixDecay.cxx:246
 RooBMixDecay.cxx:247
 RooBMixDecay.cxx:248
 RooBMixDecay.cxx:249
 RooBMixDecay.cxx:250
 RooBMixDecay.cxx:251
 RooBMixDecay.cxx:252
 RooBMixDecay.cxx:253
 RooBMixDecay.cxx:254
 RooBMixDecay.cxx:255
 RooBMixDecay.cxx:256
 RooBMixDecay.cxx:257
 RooBMixDecay.cxx:258
 RooBMixDecay.cxx:259
 RooBMixDecay.cxx:260
 RooBMixDecay.cxx:261
 RooBMixDecay.cxx:262
 RooBMixDecay.cxx:263
 RooBMixDecay.cxx:264
 RooBMixDecay.cxx:265
 RooBMixDecay.cxx:266
 RooBMixDecay.cxx:267
 RooBMixDecay.cxx:268
 RooBMixDecay.cxx:269
 RooBMixDecay.cxx:270
 RooBMixDecay.cxx:271
 RooBMixDecay.cxx:272
 RooBMixDecay.cxx:273
 RooBMixDecay.cxx:274
 RooBMixDecay.cxx:275
 RooBMixDecay.cxx:276
 RooBMixDecay.cxx:277
 RooBMixDecay.cxx:278
 RooBMixDecay.cxx:279
 RooBMixDecay.cxx:280
 RooBMixDecay.cxx:281
 RooBMixDecay.cxx:282
 RooBMixDecay.cxx:283
 RooBMixDecay.cxx:284
 RooBMixDecay.cxx:285
 RooBMixDecay.cxx:286
 RooBMixDecay.cxx:287
 RooBMixDecay.cxx:288
 RooBMixDecay.cxx:289
 RooBMixDecay.cxx:290
 RooBMixDecay.cxx:291
 RooBMixDecay.cxx:292
 RooBMixDecay.cxx:293
 RooBMixDecay.cxx:294
 RooBMixDecay.cxx:295
 RooBMixDecay.cxx:296
 RooBMixDecay.cxx:297
 RooBMixDecay.cxx:298
 RooBMixDecay.cxx:299
 RooBMixDecay.cxx:300