Logo ROOT  
Reference Guide
rf210_angularconv.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roofit
3 /// \notebook -js
4 ///
5 /// Addition and convolution: convolution in cyclical angular observables theta
6 ///
7 /// and construction of p.d.f in terms of transformed angular coordinates, e.g. cos(theta),
8 /// where the convolution is performed in theta rather than cos(theta)
9 ///
10 /// ```
11 /// pdf(theta) = T(theta) (x) gauss(theta)
12 /// pdf(cosTheta) = T(acos(cosTheta)) (x) gauss(acos(cosTheta))
13 /// ```
14 ///
15 /// This tutorial requires FFT3 to be enabled.
16 ///
17 /// \macro_image
18 /// \macro_output
19 /// \macro_code
20 ///
21 /// \date 04/2009
22 /// \author Wouter Verkerke
23 
24 #include "RooRealVar.h"
25 #include "RooDataSet.h"
26 #include "RooGaussian.h"
27 #include "RooGenericPdf.h"
28 #include "RooFormulaVar.h"
29 #include "RooFFTConvPdf.h"
30 #include "RooPlot.h"
31 #include "TCanvas.h"
32 #include "TAxis.h"
33 #include "TH1.h"
34 using namespace RooFit;
35 
36 void rf210_angularconv()
37 {
38  // S e t u p c o m p o n e n t p d f s
39  // ---------------------------------------
40 
41  // Define angle psi
42  RooRealVar psi("psi", "psi", 0, 3.14159268);
43 
44  // Define physics p.d.f T(psi)
45  RooGenericPdf Tpsi("Tpsi", "1+sin(2*@0)", psi);
46 
47  // Define resolution R(psi)
48  RooRealVar gbias("gbias", "gbias", 0.2, 0., 1);
49  RooRealVar greso("greso", "greso", 0.3, 0.1, 1.0);
50  RooGaussian Rpsi("Rpsi", "Rpsi", psi, gbias, greso);
51 
52  // Define cos(psi) and function psif that calculates psi from cos(psi)
53  RooRealVar cpsi("cpsi", "cos(psi)", -1, 1);
54  RooFormulaVar psif("psif", "acos(cpsi)", cpsi);
55 
56  // Define physics p.d.f. also as function of cos(psi): T(psif(cpsi)) = T(cpsi) ;
57  RooGenericPdf Tcpsi("T", "1+sin(2*@0)", psif);
58 
59  // C o n s t r u c t c o n v o l u t i o n p d f i n p s i
60  // --------------------------------------------------------------
61 
62  // Define convoluted p.d.f. as function of psi: M=[T(x)R](psi) = M(psi)
63  RooFFTConvPdf Mpsi("Mf", "Mf", psi, Tpsi, Rpsi);
64 
65  // Set the buffer fraction to zero to obtain a true cyclical convolution
66  Mpsi.setBufferFraction(0);
67 
68  // S a m p l e , f i t a n d p l o t c o n v o l u t e d p d f ( p s i )
69  // --------------------------------------------------------------------------------
70 
71  // Generate some events in observable psi
72  RooDataSet *data_psi = Mpsi.generate(psi, 10000);
73 
74  // Fit convoluted model as function of angle psi
75  Mpsi.fitTo(*data_psi);
76 
77  // Plot cos(psi) frame with Mf(cpsi)
78  RooPlot *frame1 = psi.frame(Title("Cyclical convolution in angle psi"));
79  data_psi->plotOn(frame1);
80  Mpsi.plotOn(frame1);
81 
82  // Overlay comparison to unsmeared physics p.d.f T(psi)
83  Tpsi.plotOn(frame1, LineColor(kRed));
84 
85  // C o n s t r u c t c o n v o l u t i o n p d f i n c o s ( p s i )
86  // --------------------------------------------------------------------------
87 
88  // Define convoluted p.d.f. as function of cos(psi): M=[T(x)R](psif(cpsi)) = M(cpsi)
89  //
90  // Need to give both observable psi here (for definition of convolution)
91  // and function psif here (for definition of observables, ultimately in cpsi)
92  RooFFTConvPdf Mcpsi("Mf", "Mf", psif, psi, Tpsi, Rpsi);
93 
94  // Set the buffer fraction to zero to obtain a true cyclical convolution
95  Mcpsi.setBufferFraction(0);
96 
97  // S a m p l e , f i t a n d p l o t c o n v o l u t e d p d f ( c o s p s i )
98  // --------------------------------------------------------------------------------
99 
100  // Generate some events
101  RooDataSet *data_cpsi = Mcpsi.generate(cpsi, 10000);
102 
103  // set psi constant to exclude to be a parameter of the fit
104  psi.setConstant(true);
105 
106  // Fit convoluted model as function of cos(psi)
107  Mcpsi.fitTo(*data_cpsi);
108 
109  // Plot cos(psi) frame with Mf(cpsi)
110  RooPlot *frame2 = cpsi.frame(Title("Same convolution in psi, expressed in cos(psi)"));
111  data_cpsi->plotOn(frame2);
112  Mcpsi.plotOn(frame2);
113 
114  // Overlay comparison to unsmeared physics p.d.f Tf(cpsi)
115  Tcpsi.plotOn(frame2, LineColor(kRed));
116 
117  // Draw frame on canvas
118  TCanvas *c = new TCanvas("rf210_angularconv", "rf210_angularconv", 800, 400);
119  c->Divide(2);
120  c->cd(1);
121  gPad->SetLeftMargin(0.15);
122  frame1->GetYaxis()->SetTitleOffset(1.4);
123  frame1->Draw();
124  c->cd(2);
125  gPad->SetLeftMargin(0.15);
126  frame2->GetYaxis()->SetTitleOffset(1.4);
127  frame2->Draw();
128 }
RooFormulaVar.h
c
#define c(i)
Definition: RSha256.hxx:119
RooPlot::Draw
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:691
RooGaussian.h
RooGaussian
Definition: RooGaussian.h:25
TCanvas.h
RooGenericPdf.h
RooDataSet.h
RooPlot::frame
static RooPlot * frame(const RooAbsRealLValue &var, Double_t xmin, Double_t xmax, Int_t nBins)
Create a new frame for a given variable in x.
Definition: RooPlot.cxx:249
RooFormulaVar
Definition: RooFormulaVar.h:30
RooFit
Definition: RooCFunction1Binding.h:29
RooAbsData::plotOn
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Definition: RooAbsData.cxx:547
RooPlot.h
RooPlot::GetYaxis
TAxis * GetYaxis() const
Definition: RooPlot.cxx:1256
RooPlot
Definition: RooPlot.h:44
RooRealVar.h
kRed
@ kRed
Definition: Rtypes.h:66
RooFit::LineColor
RooCmdArg LineColor(Color_t color)
Definition: RooGlobalFunc.cxx:56
RooFFTConvPdf
PDF for the numerical (FFT) convolution of two PDFs.
Definition: RooFFTConvPdf.h:25
TCanvas
Definition: TCanvas.h:23
TAxis.h
RooGenericPdf
Definition: RooGenericPdf.h:25
RooFFTConvPdf.h
gPad
#define gPad
Definition: TVirtualPad.h:287
RooDataSet
Definition: RooDataSet.h:33
TAttAxis::SetTitleOffset
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition: TAttAxis.cxx:293
RooRealVar
Definition: RooRealVar.h:35
TH1.h
RooFit::Title
RooCmdArg Title(const char *name)
Definition: RooGlobalFunc.cxx:173