Logo ROOT  
Reference Guide
testUnfold2.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_unfold
3/// \notebook
4/// Test program as an example for a user specific regularisation scheme.
5///
6/// 1. Generate Monte Carlo and Data events
7/// The events consist of:
8/// - signal
9/// - background
10///
11/// The signal is a resonance. It is generated with a Breit-Wigner,
12/// smeared by a Gaussian
13///
14/// 2. Unfold the data. The result is:
15/// - The background level
16/// - The shape of the resonance, corrected for detector effects
17///
18/// The regularisation is done on the curvature, excluding the bins
19/// near the peak.
20///
21/// 3. produce some plots
22///
23/// \macro_output
24/// \macro_code
25///
26/// **Version 17.6, in parallel to changes in TUnfold**
27///
28/// #### History:
29/// - Version 17.5, in parallel to changes in TUnfold
30/// - Version 17.4, in parallel to changes in TUnfold
31/// - Version 17.3, in parallel to changes in TUnfold
32/// - Version 17.2, in parallel to changes in TUnfold
33/// - Version 17.1, in parallel to changes in TUnfold
34/// - Version 17.0, updated for changed methods in TUnfold
35/// - Version 16.1, parallel to changes in TUnfold
36/// - Version 16.0, parallel to changes in TUnfold
37/// - Version 15, with automatic L-curve scan, simplified example
38/// - Version 14, with changes in TUnfoldSys.cxx
39/// - Version 13, with changes to TUnfold.C
40/// - Version 12, with improvements to TUnfold.cxx
41/// - Version 11, print chi**2 and number of degrees of freedom
42/// - Version 10, with bug-fix in TUnfold.cxx
43/// - Version 9, with bug-fix in TUnfold.cxx, TUnfold.h
44/// - Version 8, with bug-fix in TUnfold.cxx, TUnfold.h
45/// - Version 7, with bug-fix in TUnfold.cxx, TUnfold.h
46/// - Version 6a, fix problem with dynamic array allocation under windows
47/// - Version 6, re-include class MyUnfold in the example
48/// - Version 5, move class MyUnfold to separate files
49/// - Version 4, with bug-fix in TUnfold.C
50/// - Version 3, with bug-fix in TUnfold.C
51/// - Version 2, with changed ScanLcurve() arguments
52/// - Version 1, remove L curve analysis, use ScanLcurve() method instead
53/// - Version 0, L curve analysis included here
54///
55/// This file is part of TUnfold.
56///
57/// TUnfold is free software: you can redistribute it and/or modify
58/// it under the terms of the GNU General Public License as published by
59/// the Free Software Foundation, either version 3 of the License, or
60/// (at your option) any later version.
61///
62/// TUnfold is distributed in the hope that it will be useful,
63/// but WITHOUT ANY WARRANTY; without even the implied warranty of
64/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
65/// GNU General Public License for more details.
66///
67/// You should have received a copy of the GNU General Public License
68/// along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
69///
70/// \author Stefan Schmitt DESY, 14.10.2008
71
72#include <TMath.h>
73#include <TCanvas.h>
74#include <TRandom3.h>
75#include <TFitter.h>
76#include <TF1.h>
77#include <TStyle.h>
78#include <TVector.h>
79#include <TGraph.h>
80#include "TUnfold.h"
81
82using namespace std;
83
84
85TRandom *rnd=0;
86
87// generate an event
88// output:
89// negative mass: background event
90// positive mass: signal event
91Double_t GenerateEvent(Double_t bgr, // relative fraction of background
92 Double_t mass, // peak position
93 Double_t gamma) // peak width
94{
95 Double_t t;
96 if(rnd->Rndm()>bgr) {
97 // generate signal event
98 // with positive mass
99 do {
100 do {
101 t=rnd->Rndm();
102 } while(t>=1.0);
103 t=TMath::Tan((t-0.5)*TMath::Pi())*gamma+mass;
104 } while(t<=0.0);
105 return t;
106 } else {
107 // generate background event
108 // generate events following a power-law distribution
109 // f(E) = K * TMath::power((E0+E),N0)
110 static Double_t const E0=2.4;
111 static Double_t const N0=2.9;
112 do {
113 do {
114 t=rnd->Rndm();
115 } while(t>=1.0);
116 // the mass is returned negative
117 // In our example a convenient way to indicate it is a background event.
118 t= -(TMath::Power(1.-t,1./(1.-N0))-1.0)*E0;
119 } while(t>=0.0);
120 return t;
121 }
122}
123
124// smear the event to detector level
125// input:
126// mass on generator level (mTrue>0 !)
127// output:
128// mass on detector level
129Double_t DetectorEvent(Double_t mTrue) {
130 // smear by double-gaussian
131 static Double_t frac=0.1;
132 static Double_t wideBias=0.03;
133 static Double_t wideSigma=0.5;
134 static Double_t smallBias=0.0;
135 static Double_t smallSigma=0.1;
136 if(rnd->Rndm()>frac) {
137 return rnd->Gaus(mTrue+smallBias,smallSigma);
138 } else {
139 return rnd->Gaus(mTrue+wideBias,wideSigma);
140 }
141}
142
143int testUnfold2()
144{
145 // switch on histogram errors
147
148 // random generator
149 rnd=new TRandom3();
150
151 // data and MC luminosity, cross-section
152 Double_t const luminosityData=100000;
153 Double_t const luminosityMC=1000000;
154 Double_t const crossSection=1.0;
155
156 Int_t const nDet=250;
157 Int_t const nGen=100;
158 Double_t const xminDet=0.0;
159 Double_t const xmaxDet=10.0;
160 Double_t const xminGen=0.0;
161 Double_t const xmaxGen=10.0;
162
163 //============================================
164 // generate MC distribution
165 //
166 TH1D *histMgenMC=new TH1D("MgenMC",";mass(gen)",nGen,xminGen,xmaxGen);
167 TH1D *histMdetMC=new TH1D("MdetMC",";mass(det)",nDet,xminDet,xmaxDet);
168 TH2D *histMdetGenMC=new TH2D("MdetgenMC",";mass(det);mass(gen)",nDet,xminDet,xmaxDet,
169 nGen,xminGen,xmaxGen);
170 Int_t neventMC=rnd->Poisson(luminosityMC*crossSection);
171 for(Int_t i=0;i<neventMC;i++) {
172 Double_t mGen=GenerateEvent(0.3, // relative fraction of background
173 4.0, // peak position in MC
174 0.2); // peak width in MC
175 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
176 // the generated mass is negative for background
177 // and positive for signal
178 // so it will be filled in the underflow bin
179 // this is very convenient for the unfolding:
180 // the unfolded result will contain the number of background
181 // events in the underflow bin
182
183 // generated MC distribution (for comparison only)
184 histMgenMC->Fill(mGen,luminosityData/luminosityMC);
185 // reconstructed MC distribution (for comparison only)
186 histMdetMC->Fill(mDet,luminosityData/luminosityMC);
187
188 // matrix describing how the generator input migrates to the
189 // reconstructed level. Unfolding input.
190 // NOTE on underflow/overflow bins:
191 // (1) the detector level under/overflow bins are used for
192 // normalisation ("efficiency" correction)
193 // in our toy example, these bins are populated from tails
194 // of the initial MC distribution.
195 // (2) the generator level underflow/overflow bins are
196 // unfolded. In this example:
197 // underflow bin: background events reconstructed in the detector
198 // overflow bin: signal events generated at masses > xmaxDet
199 // for the unfolded result these bins will be filled
200 // -> the background normalisation will be contained in the underflow bin
201 histMdetGenMC->Fill(mDet,mGen,luminosityData/luminosityMC);
202 }
203
204 //============================================
205 // generate data distribution
206 //
207 TH1D *histMgenData=new TH1D("MgenData",";mass(gen)",nGen,xminGen,xmaxGen);
208 TH1D *histMdetData=new TH1D("MdetData",";mass(det)",nDet,xminDet,xmaxDet);
209 Int_t neventData=rnd->Poisson(luminosityData*crossSection);
210 for(Int_t i=0;i<neventData;i++) {
211 Double_t mGen=GenerateEvent(0.4, // relative fraction of background
212 3.8, // peak position
213 0.15); // peak width
214 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
215 // generated data mass for comparison plots
216 // for real data, we do not have this histogram
217 histMgenData->Fill(mGen);
218
219 // reconstructed mass, unfolding input
220 histMdetData->Fill(mDet);
221 }
222
223 //=========================================================================
224 // set up the unfolding
225 TUnfold unfold(histMdetGenMC,TUnfold::kHistMapOutputVert,
227 // regularisation
228 //----------------
229 // the regularisation is done on the curvature (2nd derivative) of
230 // the output distribution
231 //
232 // One has to exclude the bins near the peak of the Breit-Wigner,
233 // because there the curvature is high
234 // (and the regularisation eventually could enforce a small
235 // curvature, thus biasing result)
236 //
237 // in real life, the parameters below would have to be optimized,
238 // depending on the data peak position and width
239 // Or maybe one finds a different regularisation scheme... this is
240 // just an example...
241 Double_t estimatedPeakPosition=3.8;
242 Int_t nPeek=3;
244 // calculate bin number corresponding to estimated peak position
245 Int_t iPeek=(Int_t)(nGen*(estimatedPeakPosition-xminGen)/(xmaxGen-xminGen)
246 // offset 1.5
247 // accounts for start bin 1
248 // and rounding errors +0.5
249 +1.5);
250 // regularize output bins 1..iPeek-nPeek
251 unfold.RegularizeBins(1,1,iPeek-nPeek,regMode);
252 // regularize output bins iPeek+nPeek..nGen
253 unfold.RegularizeBins(iPeek+nPeek,1,nGen-(iPeek+nPeek),regMode);
254
255 // unfolding
256 //-----------
257
258 // set input distribution and bias scale (=0)
259 if(unfold.SetInput(histMdetData,0.0)>=10000) {
260 std::cout<<"Unfolding result may be wrong\n";
261 }
262
263 // do the unfolding here
264 Double_t tauMin=0.0;
265 Double_t tauMax=0.0;
266 Int_t nScan=30;
267 Int_t iBest;
268 TSpline *logTauX,*logTauY;
269 TGraph *lCurve;
270 // this method scans the parameter tau and finds the kink in the L curve
271 // finally, the unfolding is done for the "best" choice of tau
272 iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
273 std::cout<<"tau="<<unfold.GetTau()<<"\n";
274 std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
275 <<" / "<<unfold.GetNdf()<<"\n";
276
277 // save point corresponding to the kink in the L curve as TGraph
278 Double_t t[1],x[1],y[1];
279 logTauX->GetKnot(iBest,t[0],x[0]);
280 logTauY->GetKnot(iBest,t[0],y[0]);
281 TGraph *bestLcurve=new TGraph(1,x,y);
282 TGraph *bestLogTauX=new TGraph(1,t,x);
283
284 //============================================================
285 // extract unfolding results into histograms
286
287 // set up a bin map, excluding underflow and overflow bins
288 // the binMap relates the the output of the unfolding to the final
289 // histogram bins
290 Int_t *binMap=new Int_t[nGen+2];
291 for(Int_t i=1;i<=nGen;i++) binMap[i]=i;
292 binMap[0]=-1;
293 binMap[nGen+1]=-1;
294
295 TH1D *histMunfold=new TH1D("Unfolded",";mass(gen)",nGen,xminGen,xmaxGen);
296 unfold.GetOutput(histMunfold,binMap);
297 TH1D *histMdetFold=new TH1D("FoldedBack","mass(det)",nDet,xminDet,xmaxDet);
298 unfold.GetFoldedOutput(histMdetFold);
299
300 // store global correlation coefficients
301 TH1D *histRhoi=new TH1D("rho_I","mass",nGen,xminGen,xmaxGen);
302 unfold.GetRhoI(histRhoi,binMap);
303
304 delete[] binMap;
305 binMap=0;
306
307 //=====================================================================
308 // plot some histograms
310
311 // produce some plots
312 output.Divide(3,2);
313
314 // Show the matrix which connects input and output
315 // There are overflow bins at the bottom, not shown in the plot
316 // These contain the background shape.
317 // The overflow bins to the left and right contain
318 // events which are not reconstructed. These are necessary for proper MC
319 // normalisation
320 output.cd(1);
321 histMdetGenMC->Draw("BOX");
322
323 // draw generator-level distribution:
324 // data (red) [for real data this is not available]
325 // MC input (black) [with completely wrong peak position and shape]
326 // unfolded data (blue)
327 output.cd(2);
328 histMunfold->SetLineColor(kBlue);
329 histMunfold->Draw();
330 histMgenData->SetLineColor(kRed);
331 histMgenData->Draw("SAME");
332 histMgenMC->Draw("SAME HIST");
333
334 // show detector level distributions
335 // data (red)
336 // MC (black)
337 // unfolded data (blue)
338 output.cd(3);
339 histMdetFold->SetLineColor(kBlue);
340 histMdetFold->Draw();
341 histMdetData->SetLineColor(kRed);
342 histMdetData->Draw("SAME");
343 histMdetMC->Draw("SAME HIST");
344
345 // show correlation coefficients
346 // all bins outside the peak are found to be highly correlated
347 // But they are compatible with zero anyway
348 // If the peak shape is fitted,
349 // these correlations have to be taken into account, see example
350 output.cd(4);
351 histRhoi->Draw();
352
353 // show rhoi_max(tau) distribution
354 output.cd(5);
355 logTauX->Draw();
356 bestLogTauX->SetMarkerColor(kRed);
357 bestLogTauX->Draw("*");
358
359 output.cd(6);
360 lCurve->Draw("AL");
361 bestLcurve->SetMarkerColor(kRed);
362 bestLcurve->Draw("*");
363
364 output.SaveAs("testUnfold2.ps");
365 return 0;
366}
int Int_t
Definition: RtypesCore.h:43
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
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:27
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
virtual void Draw(Option_t *chopt="")
Draw this graph with its current attributes.
Definition: TGraph.cxx:760
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
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
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
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
Base class for spline implementation containing the Draw/Paint methods.
Definition: TSpline.h:22
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TSpline.cxx:97
virtual void GetKnot(Int_t i, Double_t &x, Double_t &y) const =0
An algorithm to unfold distributions from detector to truth level.
Definition: TUnfold.h:103
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:119
@ kRegModeNone
no regularisation, or defined later by RegularizeXXX() methods
Definition: TUnfold.h:122
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
Definition: TUnfold.h:131
@ kHistMapOutputVert
truth level on y-axis of the response matrix
Definition: TUnfold.h:145
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
double gamma(double x)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
constexpr Double_t Pi()
Definition: TMath.h:38
Double_t Tan(Double_t)
Definition: TMath.h:635
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
static void output(int code)
Definition: gifencode.c:226