Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TSVDUnfoldExample.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebook
4/// Data unfolding using Singular Value Decomposition
5///
6/// TSVDUnfold example
7///
8/// Data unfolding using Singular Value Decomposition (hep-ph/9509307)
9///
10/// Example distribution and smearing model from Tim Adye (RAL)
11///
12/// \macro_image
13/// \macro_code
14///
15/// \authors Kerstin Tackmann, Andreas Hoecker, Heiko Lacker
16
17#include <iostream>
18
19#include "TROOT.h"
20#include "TSystem.h"
21#include "TStyle.h"
22#include "TRandom3.h"
23#include "TString.h"
24#include "TMath.h"
25#include "TH1D.h"
26#include "TH2D.h"
27#include "TLegend.h"
28#include "TCanvas.h"
29#include "TColor.h"
30#include "TLine.h"
31
32#include "TSVDUnfold.h"
33
34
35Double_t Reconstruct( Double_t xt, TRandom3& R )
36{
37 // apply some Gaussian smearing + bias and efficiency corrections to fake reconstruction
38 const Double_t cutdummy = -99999.0;
39 Double_t xeff = 0.3 + (1.0 - 0.3)/20.0*(xt + 10.0); // efficiency
40 Double_t x = R.Rndm();
41 if (x > xeff) return cutdummy;
42 else {
43 Double_t xsmear= R.Gaus(-2.5,0.2); // bias and smear
44 return xt+xsmear;
45 }
46}
47
48void TSVDUnfoldExample()
49{
50 gROOT->SetStyle("Plain");
52
53 TRandom3 R;
54
55 const Double_t cutdummy= -99999.0;
56
57 // --------------------------------------
58 // Data/MC toy generation
59 //
60 // The MC input
61 Int_t nbins = 40;
62 TH1D *xini = new TH1D("xini", "MC truth", nbins, -10.0, 10.0);
63 TH1D *bini = new TH1D("bini", "MC reco", nbins, -10.0, 10.0);
64 TH2D *Adet = new TH2D("Adet", "detector response", nbins, -10.0, 10.0, nbins, -10.0, 10.0);
65
66 // Data
67 TH1D *data = new TH1D("data", "data", nbins, -10.0, 10.0);
68 // Data "truth" distribution to test the unfolding
69 TH1D *datatrue = new TH1D("datatrue", "data truth", nbins, -10.0, 10.0);
70 // Statistical covariance matrix
71 TH2D *statcov = new TH2D("statcov", "covariance matrix", nbins, -10.0, 10.0, nbins, -10.0, 10.0);
72
73 // Fill the MC using a Breit-Wigner, mean 0.3 and width 2.5.
74 for (Int_t i= 0; i<100000; i++) {
75 Double_t xt = R.BreitWigner(0.3, 2.5);
76 xini->Fill(xt);
77 Double_t x = Reconstruct( xt, R );
78 if (x != cutdummy) {
79 Adet->Fill(x, xt);
80 bini->Fill(x);
81 }
82 }
83
84 // Fill the "data" with a Gaussian, mean 0 and width 2.
85 for (Int_t i=0; i<10000; i++) {
86 Double_t xt = R.Gaus(0.0, 2.0);
87 datatrue->Fill(xt);
88 Double_t x = Reconstruct( xt, R );
89 if (x != cutdummy)
90 data->Fill(x);
91 }
92
93 cout << "Created toy distributions and errors for: " << endl;
94 cout << "... \"true MC\" and \"reconstructed (smeared) MC\"" << endl;
95 cout << "... \"true data\" and \"reconstructed (smeared) data\"" << endl;
96 cout << "... the \"detector response matrix\"" << endl;
97
98 // Fill the data covariance matrix
99 for (int i=1; i<=data->GetNbinsX(); i++) {
100 statcov->SetBinContent(i,i,data->GetBinError(i)*data->GetBinError(i));
101 }
102
103 // ----------------------------
104 // Here starts the actual unfolding
105 //
106 // Create TSVDUnfold object and initialise
107 TSVDUnfold *tsvdunf = new TSVDUnfold( data, statcov, bini, xini, Adet );
108
109 // It is possible to normalise unfolded spectrum to unit area
110 tsvdunf->SetNormalize( kFALSE ); // no normalisation here
111
112 // Perform the unfolding with regularisation parameter kreg = 13
113 // - the larger kreg, the finer grained the unfolding, but the more fluctuations occur
114 // - the smaller kreg, the stronger is the regularisation and the bias
115 TH1D* unfres = tsvdunf->Unfold( 13 );
116
117 // Get the distribution of the d to cross check the regularization
118 // - choose kreg to be the point where |d_i| stop being statistically significantly >>1
119 TH1D* ddist = tsvdunf->GetD();
120
121 // Get the distribution of the singular values
122 TH1D* svdist = tsvdunf->GetSV();
123
124 // Compute the error matrix for the unfolded spectrum using toy MC
125 // using the measured covariance matrix as input to generate the toys
126 // 100 toys should usually be enough
127 // The same method can be used for different covariance matrices separately.
128 TH2D* ustatcov = tsvdunf->GetUnfoldCovMatrix( statcov, 100 );
129
130 // Now compute the error matrix on the unfolded distribution originating
131 // from the finite detector matrix statistics
132 TH2D* uadetcov = tsvdunf->GetAdetCovMatrix( 100 );
133
134 // Sum up the two (they are uncorrelated)
135 ustatcov->Add( uadetcov );
136
137 //Get the computed regularized covariance matrix (always corresponding to total uncertainty passed in constructor) and add uncertainties from finite MC statistics.
138 TH2D* utaucov = tsvdunf->GetXtau();
139 utaucov->Add( uadetcov );
140
141 //Get the computed inverse of the covariance matrix
142 TH2D* uinvcov = tsvdunf->GetXinv();
143
144
145 // ---------------------------------
146 // Only plotting stuff below
147
148 for (int i=1; i<=unfres->GetNbinsX(); i++) {
149 unfres->SetBinError(i, TMath::Sqrt(utaucov->GetBinContent(i,i)));
150 }
151
152 // Renormalize just to be able to plot on the same scale
153 xini->Scale(0.7*datatrue->Integral()/xini->Integral());
154
155 TLegend *leg = new TLegend(0.58,0.60,0.99,0.88);
156 leg->SetBorderSize(0);
157 leg->SetFillColor(0);
158 leg->SetFillStyle(0);
159 leg->AddEntry(unfres,"Unfolded Data","p");
160 leg->AddEntry(datatrue,"True Data","l");
161 leg->AddEntry(data,"Reconstructed Data","l");
162 leg->AddEntry(xini,"True MC","l");
163
164 TCanvas *c1 = new TCanvas( "c1", "Unfolding toy example with TSVDUnfold", 1000, 900 );
165
166 c1->Divide(1,2);
167 TVirtualPad * c11 = c1->cd(1);
168
169 TH1D* frame = new TH1D( *unfres );
170 frame->SetTitle( "Unfolding toy example with TSVDUnfold" );
171 frame->GetXaxis()->SetTitle( "x variable" );
172 frame->GetYaxis()->SetTitle( "Events" );
173 frame->GetXaxis()->SetTitleOffset( 1.25 );
174 frame->GetYaxis()->SetTitleOffset( 1.29 );
175 frame->Draw();
176
177 data->SetLineStyle(2);
178 data->SetLineColor(4);
179 data->SetLineWidth(2);
180 unfres->SetMarkerStyle(20);
181 datatrue->SetLineColor(2);
182 datatrue->SetLineWidth(2);
183 xini->SetLineStyle(2);
184 xini->SetLineColor(8);
185 xini->SetLineWidth(2);
186 // ------------------------------------------------------------
187
188 // add histograms
189 unfres->Draw("same");
190 datatrue->Draw("same");
191 data->Draw("same");
192 xini->Draw("same");
193
194 leg->Draw();
195
196 // covariance matrix
197 TVirtualPad * c12 = c1->cd(2);
198 c12->Divide(2,1);
199 TVirtualPad * c2 = c12->cd(1);
200 c2->SetRightMargin ( 0.15 );
201
202 TH2D* covframe = new TH2D( *ustatcov );
203 covframe->SetTitle( "TSVDUnfold covariance matrix" );
204 covframe->GetXaxis()->SetTitle( "x variable" );
205 covframe->GetYaxis()->SetTitle( "x variable" );
206 covframe->GetXaxis()->SetTitleOffset( 1.25 );
207 covframe->GetYaxis()->SetTitleOffset( 1.29 );
208 covframe->Draw();
209
210 ustatcov->SetLineWidth( 2 );
211 ustatcov->Draw( "colzsame" );
212
213 // distribution of the d quantity
214 TVirtualPad * c3 = c12->cd(2);
215 c3->SetLogy();
216
217 TLine *line = new TLine( 0.,1.,40.,1. );
218 line->SetLineStyle(2);
219
220 TH1D* dframe = new TH1D( *ddist );
221 dframe->SetTitle( "TSVDUnfold |d_{i}|" );
222 dframe->GetXaxis()->SetTitle( "i" );
223 dframe->GetYaxis()->SetTitle( "|d_{i}|" );
224 dframe->GetXaxis()->SetTitleOffset( 1.25 );
225 dframe->GetYaxis()->SetTitleOffset( 1.29 );
226 dframe->SetMinimum( 0.001 );
227 dframe->Draw();
228
229 ddist->SetLineWidth( 2 );
230 ddist->Draw( "same" );
231 line->Draw();
232}
#define R(a, b, c, d, e, f, g, h, i)
Definition RSha256.hxx:110
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
#define gROOT
Definition TROOT.h:406
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition TAttAxis.cxx:298
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition TAttLine.h:42
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:40
virtual void SetRightMargin(Float_t rightmargin)
Set Pad right margin in fraction of the pad width.
Definition TAttPad.cxx:119
The Canvas class.
Definition TCanvas.h:23
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:669
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6686
TAxis * GetXaxis()
Definition TH1.h:324
virtual Int_t GetNbinsX() const
Definition TH1.h:297
virtual Bool_t Add(TF1 *h1, Double_t c1=1, Option_t *option="")
Performs the operation: this = this + c1*f1 if errors are defined (see TH1::Sumw2),...
Definition TH1.cxx:826
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition TH1.cxx:9174
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3344
TAxis * GetYaxis()
Definition TH1.h:325
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3066
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:404
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition TH1.cxx:7909
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6572
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:357
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Definition TH2.cxx:2616
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition TH2.h:93
Int_t Fill(Double_t) override
Invalid Fill method.
Definition TH2.cxx:393
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
Use the TLine constructor to create a simple line.
Definition TLine.h:22
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition TNamed.cxx:164
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
Random number generator class based on M.
Definition TRandom3.h:27
SVD Approach to Data Unfolding.
Definition TSVDUnfold.h:46
TH1D * GetSV() const
Returns singular values vector.
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
void SetNormalize(Bool_t normalize)
Definition TSVDUnfold.h:66
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
TH2D * GetUnfoldCovMatrix(const TH2D *cov, Int_t ntoys, Int_t seed=1)
Determine for given input error matrix covariance matrix of unfolded spectrum from toy simulation giv...
TH2D * GetAdetCovMatrix(Int_t ntoys, Int_t seed=1)
Determine covariance matrix of unfolded spectrum from finite statistics in response matrix using pseu...
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition TStyle.cxx:1636
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition TVirtualPad.h:51
virtual TVirtualPad * cd(Int_t subpadnumber=0)=0
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)=0
virtual void SetLogy(Int_t value=1)=0
TLine * line
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
leg
Definition legend1.C:34
return c2
Definition legend2.C:14
return c3
Definition legend3.C:15
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662