Logo ROOT   6.08/07
Reference Guide
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 
35 Double_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 
48 void TSVDUnfoldExample()
49 {
50  gROOT->SetStyle("Plain");
51  gStyle->SetOptStat(0);
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 }
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title Offset is a correction factor with respect to the "s...
Definition: TAttAxis.cxx:262
SVD Approach to Data Unfolding.
Definition: TSVDUnfold.h:54
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:49
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:5936
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3125
Random number generator class based on M.
Definition: TRandom3.h:29
virtual void SetLogy(Int_t value=1)=0
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:27
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:94
TLine * line
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
return c1
Definition: legend1.C:41
R__EXTERN TStyle * gStyle
Definition: TStyle.h:418
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:373
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:659
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7124
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
Definition: TSVDUnfold.cxx:243
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:400
#define gROOT
Definition: TROOT.h:364
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...
Definition: TSVDUnfold.cxx:411
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:255
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:44
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
Definition: TSVDUnfold.cxx:610
int nbins[3]
virtual TVirtualPad * cd(Int_t subpadnumber=0)=0
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
Definition: TSVDUnfold.cxx:602
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
Definition: TSVDUnfold.cxx:582
void SetNormalize(Bool_t normalize)
Definition: TSVDUnfold.h:74
Double_t x[n]
Definition: legend1.C:17
TH1D * GetSV() const
Returns singular values vector.
Definition: TSVDUnfold.cxx:593
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition: TVirtualPad.h:59
virtual void SetBinError(Int_t bin, Double_t error)
See convention for numbering bins in TH1::GetBin.
Definition: TH1.cxx:8309
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:46
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:42
TH2D * GetAdetCovMatrix(Int_t ntoys, Int_t seed=1)
Determine covariance matrix of unfolded spectrum from finite statistics in response matrix using pseu...
Definition: TSVDUnfold.cxx:517
A simple line.
Definition: TLine.h:33
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:45
TAxis * GetYaxis()
Definition: TH1.h:325
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
The Canvas class.
Definition: TCanvas.h:41
return c2
Definition: legend2.C:14
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:280
leg
Definition: legend1.C:34
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:307
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma. ...
Definition: TRandom.cxx:186
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:48
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), errors are also recalculated.
Definition: TH1.cxx:770
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 SetRightMargin(Float_t rightmargin)
Set Pad right margin in fraction of the pad width.
Definition: TAttPad.cxx:120
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)
Automatic pad generation by division.
Definition: TPad.cxx:1089
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:1257
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:301
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:90
virtual void SetTitle(const char *title)
Change (i.e.
Definition: TH1.cxx:6027
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2475
virtual Int_t GetNbinsX() const
Definition: TH1.h:301
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:155
return c3
Definition: legend3.C:15
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
virtual void SetBorderSize(Int_t bordersize=4)
Definition: TPave.h:74
TAxis * GetXaxis()
Definition: TH1.h:324
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8173
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:296