Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
rf708_bphysics.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 /// Speecial p.d.f.'s: special decay pdf for B physics with mixing and/or CP violation
5 ///
6 /// \macro_image
7 /// \macro_output
8 /// \macro_code
9 /// \author 07/2008 - Wouter Verkerke
10 
11 #include "RooRealVar.h"
12 #include "RooDataSet.h"
13 #include "RooGaussian.h"
14 #include "RooConstVar.h"
15 #include "RooCategory.h"
16 #include "RooBMixDecay.h"
17 #include "RooBCPEffDecay.h"
18 #include "RooBDecay.h"
19 #include "RooFormulaVar.h"
20 #include "RooTruthModel.h"
21 #include "TCanvas.h"
22 #include "TAxis.h"
23 #include "RooPlot.h"
24 using namespace RooFit;
25 
26 void rf708_bphysics()
27 {
28  // -------------------------------------
29  // B - D e c a y w i t h m i x i n g
30  // =====================================
31 
32  // C o n s t r u c t p d f
33  // -------------------------
34 
35  // Observable
36  RooRealVar dt("dt", "dt", -10, 10);
37  dt.setBins(40);
38 
39  // Parameters
40  RooRealVar dm("dm", "delta m(B0)", 0.472);
41  RooRealVar tau("tau", "tau (B0)", 1.547);
42  RooRealVar w("w", "flavour mistag rate", 0.1);
43  RooRealVar dw("dw", "delta mistag rate for B0/B0bar", 0.1);
44 
45  RooCategory mixState("mixState", "B0/B0bar mixing state");
46  mixState.defineType("mixed", -1);
47  mixState.defineType("unmixed", 1);
48 
49  RooCategory tagFlav("tagFlav", "Flavour of the tagged B0");
50  tagFlav.defineType("B0", 1);
51  tagFlav.defineType("B0bar", -1);
52 
53  // Use delta function resolution model
54  RooTruthModel truthModel("tm", "truth model", dt);
55 
56  // Construct Bdecay with mixing
57  RooBMixDecay bmix("bmix", "decay", dt, mixState, tagFlav, tau, dm, w, dw, truthModel, RooBMixDecay::DoubleSided);
58 
59  // P l o t p d f i n v a r i o u s s l i c e s
60  // ---------------------------------------------------
61 
62  // Generate some data
63  RooDataSet *data = bmix.generate(RooArgSet(dt, mixState, tagFlav), 10000);
64 
65  // Plot B0 and B0bar tagged data separately
66  // For all plots below B0 and B0 tagged data will look somewhat differently
67  // if the flavor tagging mistag rate for B0 and B0 is different (i.e. dw!=0)
68  RooPlot *frame1 = dt.frame(Title("B decay distribution with mixing (B0/B0bar)"));
69 
70  data->plotOn(frame1, Cut("tagFlav==tagFlav::B0"));
71  bmix.plotOn(frame1, Slice(tagFlav, "B0"));
72 
73  data->plotOn(frame1, Cut("tagFlav==tagFlav::B0bar"), MarkerColor(kCyan));
74  bmix.plotOn(frame1, Slice(tagFlav, "B0bar"), LineColor(kCyan));
75 
76  // Plot mixed slice for B0 and B0bar tagged data separately
77  RooPlot *frame2 = dt.frame(Title("B decay distribution of mixed events (B0/B0bar)"));
78 
79  data->plotOn(frame2, Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0"));
80  bmix.plotOn(frame2, Slice(tagFlav, "B0"), Slice(mixState, "mixed"));
81 
82  data->plotOn(frame2, Cut("mixState==mixState::mixed&&tagFlav==tagFlav::B0bar"), MarkerColor(kCyan));
83  bmix.plotOn(frame2, Slice(tagFlav, "B0bar"), Slice(mixState, "mixed"), LineColor(kCyan));
84 
85  // Plot unmixed slice for B0 and B0bar tagged data separately
86  RooPlot *frame3 = dt.frame(Title("B decay distribution of unmixed events (B0/B0bar)"));
87 
88  data->plotOn(frame3, Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0"));
89  bmix.plotOn(frame3, Slice(tagFlav, "B0"), Slice(mixState, "unmixed"));
90 
91  data->plotOn(frame3, Cut("mixState==mixState::unmixed&&tagFlav==tagFlav::B0bar"), MarkerColor(kCyan));
92  bmix.plotOn(frame3, Slice(tagFlav, "B0bar"), Slice(mixState, "unmixed"), LineColor(kCyan));
93 
94  // -------------------------------------------------
95  // B - D e c a y w i t h C P v i o l a t i o n
96  // =================================================
97 
98  // C o n s t r u c t p d f
99  // -------------------------
100 
101  // Additional parameters needed for B decay with CPV
102  RooRealVar CPeigen("CPeigen", "CP eigen value", -1);
103  RooRealVar absLambda("absLambda", "|lambda|", 1, 0, 2);
104  RooRealVar argLambda("absLambda", "|lambda|", 0.7, -1, 1);
105  RooRealVar effR("effR", "B0/B0bar reco efficiency ratio", 1);
106 
107  // Construct Bdecay with CP violation
108  RooBCPEffDecay bcp("bcp", "bcp", dt, tagFlav, tau, dm, w, CPeigen, absLambda, argLambda, effR, dw, truthModel,
109  RooBCPEffDecay::DoubleSided);
110 
111  // P l o t s c e n a r i o 1 - s i n ( 2 b ) = 0 . 7 , | l | = 1
112  // ---------------------------------------------------------------------------
113 
114  // Generate some data
115  RooDataSet *data2 = bcp.generate(RooArgSet(dt, tagFlav), 10000);
116 
117  // Plot B0 and B0bar tagged data separately
118  RooPlot *frame4 = dt.frame(Title("B decay distribution with CPV(|l|=1,Im(l)=0.7) (B0/B0bar)"));
119 
120  data2->plotOn(frame4, Cut("tagFlav==tagFlav::B0"));
121  bcp.plotOn(frame4, Slice(tagFlav, "B0"));
122 
123  data2->plotOn(frame4, Cut("tagFlav==tagFlav::B0bar"), MarkerColor(kCyan));
124  bcp.plotOn(frame4, Slice(tagFlav, "B0bar"), LineColor(kCyan));
125 
126  // P l o t s c e n a r i o 2 - s i n ( 2 b ) = 0 . 7 , | l | = 0 . 7
127  // -------------------------------------------------------------------------------
128 
129  absLambda = 0.7;
130 
131  // Generate some data
132  RooDataSet *data3 = bcp.generate(RooArgSet(dt, tagFlav), 10000);
133 
134  // Plot B0 and B0bar tagged data separately (sin2b = 0.7 plus direct CPV |l|=0.5)
135  RooPlot *frame5 = dt.frame(Title("B decay distribution with CPV(|l|=0.7,Im(l)=0.7) (B0/B0bar)"));
136 
137  data3->plotOn(frame5, Cut("tagFlav==tagFlav::B0"));
138  bcp.plotOn(frame5, Slice(tagFlav, "B0"));
139 
140  data3->plotOn(frame5, Cut("tagFlav==tagFlav::B0bar"), MarkerColor(kCyan));
141  bcp.plotOn(frame5, Slice(tagFlav, "B0bar"), LineColor(kCyan));
142 
143  // ---------------------------------------------------------------------------
144  // G e n e r i c B d e c a y w i t h u s e r c o e f f i c i e n t s
145  // ===========================================================================
146 
147  // C o n s t r u c t p d f
148  // -------------------------
149 
150  // Model parameters
151  RooRealVar DGbG("DGbG", "DGamma/GammaAvg", 0.5, -1, 1);
152  RooRealVar Adir("Adir", "-[1-abs(l)**2]/[1+abs(l)**2]", 0);
153  RooRealVar Amix("Amix", "2Im(l)/[1+abs(l)**2]", 0.7);
154  RooRealVar Adel("Adel", "2Re(l)/[1+abs(l)**2]", 0.7);
155 
156  // Derived input parameters for pdf
157  RooFormulaVar DG("DG", "Delta Gamma", "@1/@0", RooArgList(tau, DGbG));
158 
159  // Construct coefficient functions for sin,cos,sinh modulations of decay distribution
160  RooFormulaVar fsin("fsin", "fsin", "@0*@1*(1-2*@2)", RooArgList(Amix, tagFlav, w));
161  RooFormulaVar fcos("fcos", "fcos", "@0*@1*(1-2*@2)", RooArgList(Adir, tagFlav, w));
162  RooFormulaVar fsinh("fsinh", "fsinh", "@0", RooArgList(Adel));
163 
164  // Construct generic B decay pdf using above user coefficients
165  RooBDecay bcpg("bcpg", "bcpg", dt, tau, DG, RooConst(1), fsinh, fcos, fsin, dm, truthModel, RooBDecay::DoubleSided);
166 
167  // P l o t - I m ( l ) = 0 . 7 , R e ( l ) = 0 . 7 | l | = 1, d G / G = 0 . 5
168  // -------------------------------------------------------------------------------------
169 
170  // Generate some data
171  RooDataSet *data4 = bcpg.generate(RooArgSet(dt, tagFlav), 10000);
172 
173  // Plot B0 and B0bar tagged data separately
174  RooPlot *frame6 = dt.frame(Title("B decay distribution with CPV(Im(l)=0.7,Re(l)=0.7,|l|=1,dG/G=0.5) (B0/B0bar)"));
175 
176  data4->plotOn(frame6, Cut("tagFlav==tagFlav::B0"));
177  bcpg.plotOn(frame6, Slice(tagFlav, "B0"));
178 
179  data4->plotOn(frame6, Cut("tagFlav==tagFlav::B0bar"), MarkerColor(kCyan));
180  bcpg.plotOn(frame6, Slice(tagFlav, "B0bar"), LineColor(kCyan));
181 
182  TCanvas *c = new TCanvas("rf708_bphysics", "rf708_bphysics", 1200, 800);
183  c->Divide(3, 2);
184  c->cd(1);
185  gPad->SetLeftMargin(0.15);
186  frame1->GetYaxis()->SetTitleOffset(1.6);
187  frame1->Draw();
188  c->cd(2);
189  gPad->SetLeftMargin(0.15);
190  frame2->GetYaxis()->SetTitleOffset(1.6);
191  frame2->Draw();
192  c->cd(3);
193  gPad->SetLeftMargin(0.15);
194  frame3->GetYaxis()->SetTitleOffset(1.6);
195  frame3->Draw();
196  c->cd(4);
197  gPad->SetLeftMargin(0.15);
198  frame4->GetYaxis()->SetTitleOffset(1.6);
199  frame4->Draw();
200  c->cd(5);
201  gPad->SetLeftMargin(0.15);
202  frame5->GetYaxis()->SetTitleOffset(1.6);
203  frame5->Draw();
204  c->cd(6);
205  gPad->SetLeftMargin(0.15);
206  frame6->GetYaxis()->SetTitleOffset(1.6);
207  frame6->Draw();
208 }