Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
testUnfold5d.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_unfold
3/// \notebook
4/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
5///
6/// A toy test of the TUnfold package
7///
8/// This is an example of unfolding a two-dimensional distribution
9/// also using an auxiliary measurement to constrain some background
10///
11/// The example comprises several macros
12/// testUnfold5a.C create root files with TTree objects for
13/// signal, background and data
14/// -> write files testUnfold5_signal.root
15/// testUnfold5_background.root
16/// testUnfold5_data.root
17///
18/// testUnfold5b.C create a root file with the TUnfoldBinning objects
19/// -> write file testUnfold5_binning.root
20///
21/// testUnfold5c.C loop over trees and fill histograms based on the
22/// TUnfoldBinning objects
23/// -> read testUnfold5_binning.root
24/// testUnfold5_signal.root
25/// testUnfold5_background.root
26/// testUnfold5_data.root
27///
28/// -> write testUnfold5_histograms.root
29///
30/// testUnfold5d.C run the unfolding
31/// -> read testUnfold5_histograms.root
32/// -> write testUnfold5_result.root
33/// testUnfold5_result.ps
34///
35/// \macro_output
36/// \macro_code
37///
38/// **Version 17.6, in parallel to changes in TUnfold**
39///
40/// #### History:
41/// - Version 17.5, in parallel to changes in TUnfold
42/// - Version 17.4, in parallel to changes in TUnfold
43/// - Version 17.3, in parallel to changes in TUnfold
44/// - Version 17.2, in parallel to changes in TUnfold
45/// - Version 17.1, in parallel to changes in TUnfold
46/// - Version 17.0 example for multi-dimensional unfolding
47///
48/// This file is part of TUnfold.
49///
50/// TUnfold is free software: you can redistribute it and/or modify
51/// it under the terms of the GNU General Public License as published by
52/// the Free Software Foundation, either version 3 of the License, or
53/// (at your option) any later version.
54///
55/// TUnfold is distributed in the hope that it will be useful,
56/// but WITHOUT ANY WARRANTY; without even the implied warranty of
57/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
58/// GNU General Public License for more details.
59///
60/// You should have received a copy of the GNU General Public License
61/// along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
62///
63/// \author Stefan Schmitt DESY, 14.10.2008
64
65#include <iostream>
66#include <cmath>
67#include <TMath.h>
68#include <TCanvas.h>
69#include <TStyle.h>
70#include <TGraph.h>
71#include <TFile.h>
72#include <TH1.h>
73#include "TUnfoldDensity.h"
74
75using std::cout;
76
77// #define PRINT_MATRIX_L
78
79#define TEST_INPUT_COVARIANCE
80
81void testUnfold5d()
82{
83 // switch on histogram errors
85
86 //==============================================
87 // step 1 : open output file
88 TFile *outputFile=new TFile("testUnfold5_results.root","recreate");
89
90 //==============================================
91 // step 2 : read binning schemes and input histograms
92 TFile *inputFile=new TFile("testUnfold5_histograms.root");
93
94 outputFile->cd();
95
96 TUnfoldBinning *detectorBinning,*generatorBinning;
97
98 inputFile->GetObject("detector",detectorBinning);
99 inputFile->GetObject("generator",generatorBinning);
100
101 if((!detectorBinning)||(!generatorBinning)) {
102 cout<<"problem to read binning schemes\n";
103 }
104
105 // save binning schemes to output file
106 detectorBinning->Write();
107 generatorBinning->Write();
108
109 // read histograms
110 TH1 *histDataReco,*histDataTruth;
111 TH2 *histMCGenRec;
112
113 inputFile->GetObject("histDataReco",histDataReco);
114 inputFile->GetObject("histDataTruth",histDataTruth);
115 inputFile->GetObject("histMCGenRec",histMCGenRec);
116
117#ifdef TEST_ZERO_UNCORR_ERROR
118 // special test (bug in version 17.2 and below)
119 // set all errors in hisMCGenRec to zero
120 // -> program will crash
121 for(int i=0;i<=histMCGenRec->GetNbinsX()+1;i++) {
122 for(int j=0;j<=histMCGenRec->GetNbinsY()+1;j++) {
123 histMCGenRec->SetBinError(i,j,0.0);
124 }
125 }
126#endif
127
128 histDataReco->Write();
129 histDataTruth->Write();
130 histMCGenRec->Write();
131
132 if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {
133 cout<<"problem to read input histograms\n";
134 }
135
136 //========================
137 // Step 3: unfolding
138
139 // preserve the area
141
142 // basic choice of regularisation scheme:
143 // curvature (second derivative)
145
146 // density flags
147 TUnfoldDensity::EDensityMode densityFlags=
149
150 // detailed steering for regularisation
151 const char *REGULARISATION_DISTRIBUTION=nullptr;
152 const char *REGULARISATION_AXISSTEERING="*[B]";
153
154 // set up matrix of migrations
156 regMode,constraintMode,densityFlags,
157 generatorBinning,detectorBinning,
158 REGULARISATION_DISTRIBUTION,
159 REGULARISATION_AXISSTEERING);
160
161 // define the input vector (the measured data distribution)
162
163#ifdef TEST_INPUT_COVARIANCE
164 // special test to use input covariance matrix
165 TH2D *inputEmatrix=
166 detectorBinning->CreateErrorMatrixHistogram("input_covar",true);
167 for(int i=1;i<=inputEmatrix->GetNbinsX();i++) {
168 Double_t e=histDataReco->GetBinError(i);
169 inputEmatrix->SetBinContent(i,i,e*e);
170 // test: non-zero covariance where variance is zero
171 // if(e<=0.) inputEmatrix->SetBinContent(i,i+1,1.0);
172 }
173 unfold.SetInput(histDataReco,0.0,0.0,inputEmatrix);
174#else
175 unfold.SetInput(histDataReco /* ,0.0,1.0 */);
176#endif
177 // print matrix of regularisation conditions
178#ifdef PRINT_MATRIX_L
179 TH2 *histL= unfold.GetL("L");
180 for(Int_t j=1;j<=histL->GetNbinsY();j++) {
181 cout<<"L["<<unfold.GetLBinning()->GetBinName(j)<<"]";
182 for(Int_t i=1;i<=histL->GetNbinsX();i++) {
183 Double_t c=histL->GetBinContent(i,j);
184 if(c!=0.0) cout<<" ["<<i<<"]="<<c;
185 }
186 cout<<"\n";
187 }
188#endif
189 // run the unfolding
190 //
191 // here, tau is determined by scanning the global correlation coefficients
192
193 Int_t nScan=30;
194 TSpline *rhoLogTau=nullptr;
195 TGraph *lCurve=nullptr;
196
197 // for determining tau, scan the correlation coefficients
198 // correlation coefficients may be probed for all distributions
199 // or only for selected distributions
200 // underflow/overflow bins may be included/excluded
201 //
202 const char *SCAN_DISTRIBUTION="signal";
203 const char *SCAN_AXISSTEERING=nullptr;
204
205 Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
207 SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
208 &lCurve);
209
210 // create graphs with one point to visualize best choice of tau
211 Double_t t[1],rho[1],x[1],y[1];
212 rhoLogTau->GetKnot(iBest,t[0],rho[0]);
213 lCurve->GetPoint(iBest,x[0],y[0]);
214 TGraph *bestRhoLogTau=new TGraph(1,t,rho);
215 TGraph *bestLCurve=new TGraph(1,x,y);
216 Double_t *tAll=new Double_t[nScan],*rhoAll=new Double_t[nScan];
217 for(Int_t i=0;i<nScan;i++) {
218 rhoLogTau->GetKnot(i,tAll[i],rhoAll[i]);
219 }
220 TGraph *knots=new TGraph(nScan,tAll,rhoAll);
221
222 cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
223 <<" / "<<unfold.GetNdf()<<"\n";
224
225
226 //===========================
227 // Step 4: retrieve and plot unfolding results
228
229 // get unfolding output
230 TH1 *histDataUnfold=unfold.GetOutput("unfolded signal",nullptr,nullptr,nullptr,kFALSE);
231 // get Monte Carlo reconstructed data
232 TH1 *histMCReco=histMCGenRec->ProjectionY("histMCReco",0,-1,"e");
233 TH1 *histMCTruth=histMCGenRec->ProjectionX("histMCTruth",0,-1,"e");
234 Double_t scaleFactor=histDataTruth->GetSumOfWeights()/
235 histMCTruth->GetSumOfWeights();
236 histMCReco->Scale(scaleFactor);
237 histMCTruth->Scale(scaleFactor);
238 // get matrix of probabilities
239 TH2 *histProbability=unfold.GetProbabilityMatrix("histProbability");
240 // get global correlation coefficients
241 /* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal("histGlobalCorr",nullptr,nullptr,nullptr,kFALSE);
242 TH1 *histGlobalCorrScan=unfold.GetRhoItotal
243 ("histGlobalCorrScan",nullptr,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE);
244 /* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal("histCorrCoeff",nullptr,nullptr,nullptr,kFALSE);
245
246 TCanvas canvas;
247 canvas.Print("testUnfold5.ps[");
248
249 //========== page 1 ============
250 // unfolding control plots
251 // input, matrix, output
252 // tau-scan, global correlations, correlation coefficients
253 canvas.Clear();
254 canvas.Divide(3,2);
255
256 // (1) all bins, compare to original MC distribution
257 canvas.cd(1);
258 histDataReco->SetMinimum(0.0);
259 histDataReco->Draw("E");
260 histMCReco->SetLineColor(kBlue);
261 histMCReco->Draw("SAME HIST");
262 // (2) matrix of probabilities
263 canvas.cd(2);
264 histProbability->Draw("BOX");
265 // (3) unfolded data, data truth, MC truth
266 canvas.cd(3);
267 gPad->SetLogy();
268 histDataUnfold->Draw("E");
269 histDataTruth->SetLineColor(kBlue);
270 histDataTruth->Draw("SAME HIST");
271 histMCTruth->SetLineColor(kRed);
272 histMCTruth->Draw("SAME HIST");
273 // (4) scan of correlation vs tau
274 canvas.cd(4);
275 rhoLogTau->Draw();
276 knots->Draw("*");
277 bestRhoLogTau->SetMarkerColor(kRed);
278 bestRhoLogTau->Draw("*");
279 // (5) global correlation coefficients for the distributions
280 // used during the scan
281 canvas.cd(5);
282 //histCorrCoeff->Draw("BOX");
283 histGlobalCorrScan->Draw("HIST");
284 // (6) L-curve
285 canvas.cd(6);
286 lCurve->Draw("AL");
287 bestLCurve->SetMarkerColor(kRed);
288 bestLCurve->Draw("*");
289
290
291 canvas.Print("testUnfold5.ps");
292
293 canvas.Print("testUnfold5.ps]");
294
295}
#define c(i)
Definition RSha256.hxx:101
#define e(i)
Definition RSha256.hxx:103
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
#define gPad
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition TAttMarker.h:38
The Canvas class.
Definition TCanvas.h:23
void Clear(Option_t *option="") override
Remove all primitives from the canvas.
Definition TCanvas.cxx:734
TVirtualPad * cd(Int_t subpadnumber=0) override
Set current canvas & pad.
Definition TCanvas.cxx:716
Bool_t cd() override
Change current directory to "this" directory.
void GetObject(const char *namecycle, T *&ptr)
Get an object with proper type checking.
Definition TDirectory.h:212
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
Definition TFile.h:53
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
void Draw(Option_t *chopt="") override
Draw this graph with its current attributes.
Definition TGraph.cxx:809
virtual Int_t GetPoint(Int_t i, Double_t &x, Double_t &y) const
Get x and y values for point number i.
Definition TGraph.cxx:1511
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:59
virtual Int_t GetNbinsY() const
Definition TH1.h:298
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9027
virtual Int_t GetNbinsX() const
Definition TH1.h:297
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:9170
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
virtual void SetMinimum(Double_t minimum=-1111)
Definition TH1.h:404
static void SetDefaultSumw2(Bool_t sumw2=kTRUE)
When this static function is called with sumw2=kTRUE, all new histograms will automatically activate ...
Definition TH1.cxx:6667
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition TH1.cxx:6568
virtual Double_t GetSumOfWeights() const
Return the sum of weights excluding under/overflows.
Definition TH1.cxx:7881
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:338
Service class for 2-D histogram classes.
Definition TH2.h:30
TH1D * ProjectionY(const char *name="_py", Int_t firstxbin=0, Int_t lastxbin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along Y.
Definition TH2.cxx:2433
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Definition TH2.cxx:2573
TH1D * ProjectionX(const char *name="_px", Int_t firstybin=0, Int_t lastybin=-1, Option_t *option="") const
Project a 2-D histogram into a 1-D histogram along X.
Definition TH2.cxx:2393
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition TH2.h:89
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TObject.cxx:880
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Automatic pad generation by division.
Definition TPad.cxx:1153
void Print(const char *filename="") const override
This method is equivalent to SaveAs("filename"). See TPad::SaveAs for details.
Definition TPad.cxx:4700
Base class for spline implementation containing the Draw/Paint methods.
Definition TSpline.h:31
void Draw(Option_t *option="") override
Draw this function with its current attributes.
Definition TSpline.cxx:101
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
Create a TH2D histogram capable to hold a covariance matrix.
An algorithm to unfold distributions from detector to truth level.
@ kEScanTauRhoMax
maximum global correlation coefficient (from TUnfold::GetRhoI())
EDensityMode
choice of regularisation scale factors to cinstruct the matrix L
@ kDensityModeBinWidth
scale factors from multidimensional bin width
EConstraint
type of extra constraint
Definition TUnfold.h:109
@ kEConstraintArea
enforce preservation of the area
Definition TUnfold.h:115
ERegMode
choice of regularisation scheme
Definition TUnfold.h:119
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
Definition TUnfold.h:131
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition TUnfold.h:142
virtual void SetLogy(Int_t value=1)=0
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17