Loading [MathJax]/extensions/tex2jax.js
Logo ROOT   6.12/07
Reference Guide
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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 
56 #include "TUnfoldDensity.h"
57 
58 using namespace std;
59 
60 TRandom *rnd=0;
61 
62 Int_t GenerateGenEvent(Int_t nmax,const Double_t *probability) {
63  // choose an integer random number in the range [0,nmax]
64  // (the generator level bin)
65  // depending on the probabilities
66  // probability[0],probability[1],...probability[nmax-1]
67  Double_t f=rnd->Rndm();
68  Int_t r=0;
69  while((r<nmax)&&(f>=probability[r])) {
70  f -= probability[r];
71  r++;
72  }
73  return r;
74 }
75 
76 Double_t GenerateRecEvent(const Double_t *shapeParm) {
77  // return a coordinate (the reconstructed variable)
78  // depending on shapeParm[]
79  // shapeParm[0]: fraction of events with Gaussian distribution
80  // shapeParm[1]: mean of Gaussian
81  // shapeParm[2]: width of Gaussian
82  // (1-shapeParm[0]): fraction of events with flat distribution
83  // shapeParm[3]: minimum of flat component
84  // shapeParm[4]: maximum of flat component
85  Double_t f=rnd->Rndm();
86  Double_t r;
87  if(f<shapeParm[0]) {
88  r=rnd->Gaus(shapeParm[1],shapeParm[2]);
89  } else {
90  r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3];
91  }
92  return r;
93 }
94 
95 void testUnfold4()
96 {
97  // switch on histogram errors
99 
100  // random generator
101  rnd=new TRandom3();
102 
103  // data and MC number of events
104  Double_t const nData0= 500.0;
105  Double_t const nMC0 = 50000.0;
106 
107  // Binning
108  // reconstructed variable (0-10)
109  Int_t const nDet=15;
110  Double_t const xminDet=0.0;
111  Double_t const xmaxDet=15.0;
112 
113  // signal binning (three shapes: 0,1,2)
114  Int_t const nGen=3;
115  Double_t const xminGen=-0.5;
116  Double_t const xmaxGen= 2.5;
117 
118  // parameters
119  // fraction of events per signal shape
120  static const Double_t genFrac[]={0.3,0.6,0.1};
121 
122  // signal shapes
123  static const Double_t genShape[][5]=
124  {{1.0,2.0,1.5,0.,15.},
125  {1.0,7.0,2.5,0.,15.},
126  {0.0,0.0,0.0,0.,15.}};
127 
128  // define DATA histograms
129  // observed data distribution
130  TH1D *histDetDATA=new TH1D("Yrec",";DATA(Yrec)",nDet,xminDet,xmaxDet);
131 
132  // define MC histograms
133  // matrix of migrations
134  TH2D *histGenDetMC=new TH2D("Yrec%Xgen","MC(Xgen,Yrec)",
135  nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
136 
137  TH1D *histUnfold=new TH1D("Xgen",";DATA(Xgen)",nGen,xminGen,xmaxGen);
138 
139  TH1D **histPullNC=new TH1D* [nGen];
140  TH1D **histPullArea=new TH1D* [nGen];
141  for(int i=0;i<nGen;i++) {
142  histPullNC[i]=new TH1D(TString::Format("PullNC%d",i),"pull",15,-3.,3.);
143  histPullArea[i]=new TH1D(TString::Format("PullArea%d",i),"pull",15,-3.,3.);
144  }
145 
146  // this method is new in version 16 of TUnfold
147  cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
148 
149  for(int itoy=0;itoy<1000;itoy++) {
150  if(!(itoy %10)) cout<<"toy iteration: "<<itoy<<"\n";
151  histDetDATA->Reset();
152  histGenDetMC->Reset();
153 
154  Int_t nData=rnd->Poisson(nData0);
155  for(Int_t i=0;i<nData;i++) {
156  Int_t iGen=GenerateGenEvent(nGen,genFrac);
157  Double_t yObs=GenerateRecEvent(genShape[iGen]);
158  histDetDATA->Fill(yObs);
159  }
160 
161  Int_t nMC=rnd->Poisson(nMC0);
162  for(Int_t i=0;i<nMC;i++) {
163  Int_t iGen=GenerateGenEvent(nGen,genFrac);
164  Double_t yObs=GenerateRecEvent(genShape[iGen]);
165  histGenDetMC->Fill(iGen,yObs);
166  }
167  /* for(Int_t ix=0;ix<=histGenDetMC->GetNbinsX()+1;ix++) {
168  for(Int_t iy=0;iy<=histGenDetMC->GetNbinsY()+1;iy++) {
169  cout<<ix<<iy<<" : "<<histGenDetMC->GetBinContent(ix,iy)<<"\n";
170  }
171  } */
172  //========================
173  // unfolding
174 
175  TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz,
177  // define the input vector (the measured data distribution)
178  unfold.SetInput(histDetDATA,0.0,1.0);
179 
180  // run the unfolding
181  unfold.ScanLcurve(50,0.,0.,0,0,0);
182 
183  // fill pull distributions without constraint
184  unfold.GetOutput(histUnfold);
185 
186  for(int i=0;i<nGen;i++) {
187  histPullNC[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
188  histUnfold->GetBinError(i+1));
189  }
190 
191  // repeat unfolding on the same data, now with Area constraint
192  unfold.SetConstraint(TUnfold::kEConstraintArea);
193 
194  // run the unfolding
195  unfold.ScanLcurve(50,0.,0.,0,0,0);
196 
197  // fill pull distributions with constraint
198  unfold.GetOutput(histUnfold);
199 
200  for(int i=0;i<nGen;i++) {
201  histPullArea[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
202  histUnfold->GetBinError(i+1));
203  }
204 
205  }
206  TCanvas output;
207  output.Divide(3,2);
208 
209  gStyle->SetOptFit(1111);
210 
211  for(int i=0;i<nGen;i++) {
212  output.cd(i+1);
213  histPullNC[i]->Fit("gaus");
214  histPullNC[i]->Draw();
215  }
216  for(int i=0;i<nGen;i++) {
217  output.cd(i+4);
218  histPullArea[i]->Fit("gaus");
219  histPullArea[i]->Draw();
220  }
221  output.SaveAs("testUnfold4.ps");
222 }
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3244
Random number generator class based on M.
Definition: TRandom3.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:256
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4763
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
int Int_t
Definition: RtypesCore.h:41
STL namespace.
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:6139
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:2365
truth level on x-axis of the response matrix
Definition: TUnfold.h:143
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
regularise the amplitude of the output distribution
Definition: TUnfold.h:126
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9540
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:533
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
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:1218
static const char * GetTUnfoldVersion(void)
Return a string describing the TUnfold version.
Definition: TUnfold.cxx:3680
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:610
enforce preservation of the area
Definition: TUnfold.h:116
The Canvas class.
Definition: TCanvas.h:31
double Double_t
Definition: RtypesCore.h:55
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:290
use no extra constraint
Definition: TUnfold.h:113
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:1153
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:284
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH2.cxx:3807
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:383
virtual void SaveAs(const char *filename="", Option_t *option="") const
Save Pad contents in a file in one of various formats.
Definition: TPad.cxx:5486
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
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:3688
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8319
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290