Logo ROOT  
Reference Guide
testUnfold4.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_unfold
3/// \notebook
4/// Test program for the class TUnfoldSys.
5///
6/// Simple toy tests of the TUnfold package
7///
8/// Pseudo data (5000 events) are unfolded into three components
9/// The unfolding is performed once without and once with area constraint
10///
11/// Ideally, the pulls may show that the result is biased if no constraint
12/// is applied. This is expected because the true data errors are not known,
13/// and instead the sqrt(data) errors are used.
14///
15/// \macro_output
16/// \macro_code
17///
18/// **Version 17.6, in parallel to changes in TUnfold**
19///
20/// #### History:
21/// - Version 17.5, in parallel to changes in TUnfold
22/// - Version 17.4, in parallel to changes in TUnfold
23/// - Version 17.3, in parallel to changes in TUnfold
24/// - Version 17.2, in parallel to changes in TUnfold
25/// - Version 17.1, in parallel to changes in TUnfold
26/// - Version 16.1, parallel to changes in TUnfold
27/// - Version 16.0, parallel to changes in TUnfold
28/// - Version 15, use L-curve scan to scan the average correlation
29///
30/// This file is part of TUnfold.
31///
32/// TUnfold is free software: you can redistribute it and/or modify
33/// it under the terms of the GNU General Public License as published by
34/// the Free Software Foundation, either version 3 of the License, or
35/// (at your option) any later version.
36///
37/// TUnfold is distributed in the hope that it will be useful,
38/// but WITHOUT ANY WARRANTY; without even the implied warranty of
39/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
40/// GNU General Public License for more details.
41///
42/// You should have received a copy of the GNU General Public License
43/// along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
44///
45/// \author Stefan Schmitt DESY, 14.10.2008
46
47#include <TMath.h>
48#include <TCanvas.h>
49#include <TRandom3.h>
50#include <TFitter.h>
51#include <TF1.h>
52#include <TStyle.h>
53#include <TVector.h>
54#include <TGraph.h>
55#include <TError.h>
56
57#include "TUnfoldDensity.h"
58
59using namespace std;
60
61TRandom *rnd=0;
62
63Int_t GenerateGenEvent(Int_t nmax,const Double_t *probability) {
64 // choose an integer random number in the range [0,nmax]
65 // (the generator level bin)
66 // depending on the probabilities
67 // probability[0],probability[1],...probability[nmax-1]
68 Double_t f=rnd->Rndm();
69 Int_t r=0;
70 while((r<nmax)&&(f>=probability[r])) {
71 f -= probability[r];
72 r++;
73 }
74 return r;
75}
76
77Double_t GenerateRecEvent(const Double_t *shapeParm) {
78 // return a coordinate (the reconstructed variable)
79 // depending on shapeParm[]
80 // shapeParm[0]: fraction of events with Gaussian distribution
81 // shapeParm[1]: mean of Gaussian
82 // shapeParm[2]: width of Gaussian
83 // (1-shapeParm[0]): fraction of events with flat distribution
84 // shapeParm[3]: minimum of flat component
85 // shapeParm[4]: maximum of flat component
86 Double_t f=rnd->Rndm();
87 Double_t r;
88 if(f<shapeParm[0]) {
89 r=rnd->Gaus(shapeParm[1],shapeParm[2]);
90 } else {
91 r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3];
92 }
93 return r;
94}
95
96void testUnfold4(bool printInfo = false)
97{
98
99 // switch off printing Info messages
100 if (!printInfo) gErrorIgnoreLevel = kWarning;
101
102 // switch on histogram errors
104
105 // random generator
106 rnd=new TRandom3();
107
108 // data and MC number of events
109 Double_t const nData0= 500.0;
110 Double_t const nMC0 = 50000.0;
111
112 // Binning
113 // reconstructed variable (0-10)
114 Int_t const nDet=15;
115 Double_t const xminDet=0.0;
116 Double_t const xmaxDet=15.0;
117
118 // signal binning (three shapes: 0,1,2)
119 Int_t const nGen=3;
120 Double_t const xminGen=-0.5;
121 Double_t const xmaxGen= 2.5;
122
123 // parameters
124 // fraction of events per signal shape
125 static const Double_t genFrac[]={0.3,0.6,0.1};
126
127 // signal shapes
128 static const Double_t genShape[][5]=
129 {{1.0,2.0,1.5,0.,15.},
130 {1.0,7.0,2.5,0.,15.},
131 {0.0,0.0,0.0,0.,15.}};
132
133 // define DATA histograms
134 // observed data distribution
135 TH1D *histDetDATA=new TH1D("Yrec",";DATA(Yrec)",nDet,xminDet,xmaxDet);
136
137 // define MC histograms
138 // matrix of migrations
139 TH2D *histGenDetMC=new TH2D("Yrec%Xgen","MC(Xgen,Yrec)",
140 nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
141
142 TH1D *histUnfold=new TH1D("Xgen",";DATA(Xgen)",nGen,xminGen,xmaxGen);
143
144 TH1D **histPullNC=new TH1D* [nGen];
145 TH1D **histPullArea=new TH1D* [nGen];
146 for(int i=0;i<nGen;i++) {
147 histPullNC[i]=new TH1D(TString::Format("PullNC%d",i),"pull",15,-3.,3.);
148 histPullArea[i]=new TH1D(TString::Format("PullArea%d",i),"pull",15,-3.,3.);
149 }
150
151 // this method is new in version 16 of TUnfold
152 cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
153
154 for(int itoy=0;itoy<1000;itoy++) {
155 if(!(itoy %10)) cout<<"toy iteration: "<<itoy<<"\n";
156 histDetDATA->Reset();
157 histGenDetMC->Reset();
158
159 Int_t nData=rnd->Poisson(nData0);
160 for(Int_t i=0;i<nData;i++) {
161 Int_t iGen=GenerateGenEvent(nGen,genFrac);
162 Double_t yObs=GenerateRecEvent(genShape[iGen]);
163 histDetDATA->Fill(yObs);
164 }
165
166 Int_t nMC=rnd->Poisson(nMC0);
167 for(Int_t i=0;i<nMC;i++) {
168 Int_t iGen=GenerateGenEvent(nGen,genFrac);
169 Double_t yObs=GenerateRecEvent(genShape[iGen]);
170 histGenDetMC->Fill(iGen,yObs);
171 }
172 /* for(Int_t ix=0;ix<=histGenDetMC->GetNbinsX()+1;ix++) {
173 for(Int_t iy=0;iy<=histGenDetMC->GetNbinsY()+1;iy++) {
174 cout<<ix<<iy<<" : "<<histGenDetMC->GetBinContent(ix,iy)<<"\n";
175 }
176 } */
177 //========================
178 // unfolding
179
180 TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz,
182 // define the input vector (the measured data distribution)
183 unfold.SetInput(histDetDATA,0.0,1.0);
184
185 // run the unfolding
186 unfold.ScanLcurve(50,0.,0.,0,0,0);
187
188 // fill pull distributions without constraint
189 unfold.GetOutput(histUnfold);
190
191 for(int i=0;i<nGen;i++) {
192 histPullNC[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
193 histUnfold->GetBinError(i+1));
194 }
195
196 // repeat unfolding on the same data, now with Area constraint
197 unfold.SetConstraint(TUnfold::kEConstraintArea);
198
199 // run the unfolding
200 unfold.ScanLcurve(50,0.,0.,0,0,0);
201
202 // fill pull distributions with constraint
203 unfold.GetOutput(histUnfold);
204
205 for(int i=0;i<nGen;i++) {
206 histPullArea[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
207 histUnfold->GetBinError(i+1));
208 }
209
210 }
212 output.Divide(3,2);
213
214 gStyle->SetOptFit(1111);
215
216 for(int i=0;i<nGen;i++) {
217 output.cd(i+1);
218 histPullNC[i]->Fit("gaus");
219 histPullNC[i]->Draw();
220 }
221 for(int i=0;i<nGen;i++) {
222 output.cd(i+4);
223 histPullArea[i]->Fit("gaus");
224 histPullArea[i]->Draw();
225 }
226 output.SaveAs("testUnfold4.ps");
227}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define f(i)
Definition: RSha256.hxx:104
const Int_t kWarning
Definition: TError.h:38
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:105
R__EXTERN TStyle * gStyle
Definition: TStyle.h:410
The Canvas class.
Definition: TCanvas.h:27
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9741
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3808
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
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:6330
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH2.cxx:3816
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
Random number generator class based on M.
Definition: TRandom3.h:27
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
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:263
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:391
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition: TString.cxx:2311
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1542
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
Definition: TUnfold.cxx:3680
@ kEConstraintArea
enforce preservation of the area
Definition: TUnfold.h:115
@ kEConstraintNone
use no extra constraint
Definition: TUnfold.h:112
@ kRegModeSize
regularise the amplitude of the output distribution
Definition: TUnfold.h:125
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition: TUnfold.h:142
static void output(int code)
Definition: gifencode.c:226