Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
82
83TRandom *rnd=nullptr;
84
85// generate an event
86// output:
87// negative mass: background event
88// positive mass: signal event
89Double_t GenerateEvent(Double_t bgr, // relative fraction of background
90 Double_t mass, // peak position
91 Double_t gamma) // peak width
92{
93 Double_t t;
94 if(rnd->Rndm()>bgr) {
95 // generate signal event
96 // with positive mass
97 do {
98 do {
99 t=rnd->Rndm();
100 } while(t>=1.0);
101 t=TMath::Tan((t-0.5)*TMath::Pi())*gamma+mass;
102 } while(t<=0.0);
103 return t;
104 } else {
105 // generate background event
106 // generate events following a power-law distribution
107 // f(E) = K * TMath::power((E0+E),N0)
108 static Double_t const E0=2.4;
109 static Double_t const N0=2.9;
110 do {
111 do {
112 t=rnd->Rndm();
113 } while(t>=1.0);
114 // the mass is returned negative
115 // In our example a convenient way to indicate it is a background event.
116 t= -(TMath::Power(1.-t,1./(1.-N0))-1.0)*E0;
117 } while(t>=0.0);
118 return t;
119 }
120}
121
122// smear the event to detector level
123// input:
124// mass on generator level (mTrue>0 !)
125// output:
126// mass on detector level
127Double_t DetectorEvent(Double_t mTrue) {
128 // smear by double-gaussian
129 static Double_t frac=0.1;
130 static Double_t wideBias=0.03;
131 static Double_t wideSigma=0.5;
132 static Double_t smallBias=0.0;
133 static Double_t smallSigma=0.1;
134 if(rnd->Rndm()>frac) {
135 return rnd->Gaus(mTrue+smallBias,smallSigma);
136 } else {
137 return rnd->Gaus(mTrue+wideBias,wideSigma);
138 }
139}
140
141int testUnfold2()
142{
143 // switch on histogram errors
145
146 // random generator
147 rnd=new TRandom3();
148
149 // data and MC luminosity, cross-section
150 Double_t const luminosityData=100000;
151 Double_t const luminosityMC=1000000;
152 Double_t const crossSection=1.0;
153
154 Int_t const nDet=250;
155 Int_t const nGen=100;
156 Double_t const xminDet=0.0;
157 Double_t const xmaxDet=10.0;
158 Double_t const xminGen=0.0;
159 Double_t const xmaxGen=10.0;
160
161 //============================================
162 // generate MC distribution
163 //
164 TH1D *histMgenMC=new TH1D("MgenMC",";mass(gen)",nGen,xminGen,xmaxGen);
165 TH1D *histMdetMC=new TH1D("MdetMC",";mass(det)",nDet,xminDet,xmaxDet);
166 TH2D *histMdetGenMC=new TH2D("MdetgenMC",";mass(det);mass(gen)",nDet,xminDet,xmaxDet,
167 nGen,xminGen,xmaxGen);
168 Int_t neventMC=rnd->Poisson(luminosityMC*crossSection);
169 for(Int_t i=0;i<neventMC;i++) {
170 Double_t mGen=GenerateEvent(0.3, // relative fraction of background
171 4.0, // peak position in MC
172 0.2); // peak width in MC
173 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
174 // the generated mass is negative for background
175 // and positive for signal
176 // so it will be filled in the underflow bin
177 // this is very convenient for the unfolding:
178 // the unfolded result will contain the number of background
179 // events in the underflow bin
180
181 // generated MC distribution (for comparison only)
182 histMgenMC->Fill(mGen,luminosityData/luminosityMC);
183 // reconstructed MC distribution (for comparison only)
184 histMdetMC->Fill(mDet,luminosityData/luminosityMC);
185
186 // matrix describing how the generator input migrates to the
187 // reconstructed level. Unfolding input.
188 // NOTE on underflow/overflow bins:
189 // (1) the detector level under/overflow bins are used for
190 // normalisation ("efficiency" correction)
191 // in our toy example, these bins are populated from tails
192 // of the initial MC distribution.
193 // (2) the generator level underflow/overflow bins are
194 // unfolded. In this example:
195 // underflow bin: background events reconstructed in the detector
196 // overflow bin: signal events generated at masses > xmaxDet
197 // for the unfolded result these bins will be filled
198 // -> the background normalisation will be contained in the underflow bin
199 histMdetGenMC->Fill(mDet,mGen,luminosityData/luminosityMC);
200 }
201
202 //============================================
203 // generate data distribution
204 //
205 TH1D *histMgenData=new TH1D("MgenData",";mass(gen)",nGen,xminGen,xmaxGen);
206 TH1D *histMdetData=new TH1D("MdetData",";mass(det)",nDet,xminDet,xmaxDet);
207 Int_t neventData=rnd->Poisson(luminosityData*crossSection);
208 for(Int_t i=0;i<neventData;i++) {
209 Double_t mGen=GenerateEvent(0.4, // relative fraction of background
210 3.8, // peak position
211 0.15); // peak width
212 Double_t mDet=DetectorEvent(TMath::Abs(mGen));
213 // generated data mass for comparison plots
214 // for real data, we do not have this histogram
215 histMgenData->Fill(mGen);
216
217 // reconstructed mass, unfolding input
218 histMdetData->Fill(mDet);
219 }
220
221 //=========================================================================
222 // set up the unfolding
223 TUnfold unfold(histMdetGenMC,TUnfold::kHistMapOutputVert,
225 // regularisation
226 //----------------
227 // the regularisation is done on the curvature (2nd derivative) of
228 // the output distribution
229 //
230 // One has to exclude the bins near the peak of the Breit-Wigner,
231 // because there the curvature is high
232 // (and the regularisation eventually could enforce a small
233 // curvature, thus biasing result)
234 //
235 // in real life, the parameters below would have to be optimized,
236 // depending on the data peak position and width
237 // Or maybe one finds a different regularisation scheme... this is
238 // just an example...
239 Double_t estimatedPeakPosition=3.8;
240 Int_t nPeek=3;
242 // calculate bin number corresponding to estimated peak position
243 Int_t iPeek=(Int_t)(nGen*(estimatedPeakPosition-xminGen)/(xmaxGen-xminGen)
244 // offset 1.5
245 // accounts for start bin 1
246 // and rounding errors +0.5
247 +1.5);
248 // regularize output bins 1..iPeek-nPeek
249 unfold.RegularizeBins(1,1,iPeek-nPeek,regMode);
250 // regularize output bins iPeek+nPeek..nGen
251 unfold.RegularizeBins(iPeek+nPeek,1,nGen-(iPeek+nPeek),regMode);
252
253 // unfolding
254 //-----------
255
256 // set input distribution and bias scale (=0)
257 if(unfold.SetInput(histMdetData,0.0)>=10000) {
258 std::cout<<"Unfolding result may be wrong\n";
259 }
260
261 // do the unfolding here
262 Double_t tauMin=0.0;
263 Double_t tauMax=0.0;
264 Int_t nScan=30;
265 Int_t iBest;
266 TSpline *logTauX,*logTauY;
267 TGraph *lCurve;
268 // this method scans the parameter tau and finds the kink in the L curve
269 // finally, the unfolding is done for the "best" choice of tau
270 iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
271 std::cout<<"tau="<<unfold.GetTau()<<"\n";
272 std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
273 <<" / "<<unfold.GetNdf()<<"\n";
274
275 // save point corresponding to the kink in the L curve as TGraph
276 Double_t t[1],x[1],y[1];
277 logTauX->GetKnot(iBest,t[0],x[0]);
278 logTauY->GetKnot(iBest,t[0],y[0]);
279 TGraph *bestLcurve=new TGraph(1,x,y);
280 TGraph *bestLogTauX=new TGraph(1,t,x);
281
282 //============================================================
283 // extract unfolding results into histograms
284
285 // set up a bin map, excluding underflow and overflow bins
286 // the binMap relates the output of the unfolding to the final
287 // histogram bins
288 Int_t *binMap=new Int_t[nGen+2];
289 for(Int_t i=1;i<=nGen;i++) binMap[i]=i;
290 binMap[0]=-1;
291 binMap[nGen+1]=-1;
292
293 TH1D *histMunfold=new TH1D("Unfolded",";mass(gen)",nGen,xminGen,xmaxGen);
294 unfold.GetOutput(histMunfold,binMap);
295 TH1D *histMdetFold=new TH1D("FoldedBack","mass(det)",nDet,xminDet,xmaxDet);
296 unfold.GetFoldedOutput(histMdetFold);
297
298 // store global correlation coefficients
299 TH1D *histRhoi=new TH1D("rho_I","mass",nGen,xminGen,xmaxGen);
300 unfold.GetRhoI(histRhoi,binMap);
301
302 delete[] binMap;
303 binMap=nullptr;
304
305 //=====================================================================
306 // plot some histograms
308
309 // produce some plots
310 output.Divide(3,2);
311
312 // Show the matrix which connects input and output
313 // There are overflow bins at the bottom, not shown in the plot
314 // These contain the background shape.
315 // The overflow bins to the left and right contain
316 // events which are not reconstructed. These are necessary for proper MC
317 // normalisation
318 output.cd(1);
319 histMdetGenMC->Draw("BOX");
320
321 // draw generator-level distribution:
322 // data (red) [for real data this is not available]
323 // MC input (black) [with completely wrong peak position and shape]
324 // unfolded data (blue)
325 output.cd(2);
326 histMunfold->SetLineColor(kBlue);
327 histMunfold->Draw();
328 histMgenData->SetLineColor(kRed);
329 histMgenData->Draw("SAME");
330 histMgenMC->Draw("SAME HIST");
331
332 // show detector level distributions
333 // data (red)
334 // MC (black)
335 // unfolded data (blue)
336 output.cd(3);
337 histMdetFold->SetLineColor(kBlue);
338 histMdetFold->Draw();
339 histMdetData->SetLineColor(kRed);
340 histMdetData->Draw("SAME");
341 histMdetMC->Draw("SAME HIST");
342
343 // show correlation coefficients
344 // all bins outside the peak are found to be highly correlated
345 // But they are compatible with zero anyway
346 // If the peak shape is fitted,
347 // these correlations have to be taken into account, see example
348 output.cd(4);
349 histRhoi->Draw();
350
351 // show rhoi_max(tau) distribution
352 output.cd(5);
353 logTauX->Draw();
354 bestLogTauX->SetMarkerColor(kRed);
355 bestLogTauX->Draw("*");
356
357 output.cd(6);
358 lCurve->Draw("AL");
359 bestLcurve->SetMarkerColor(kRed);
360 bestLcurve->Draw("*");
361
362 output.SaveAs("testUnfold2.ps");
363 return 0;
364}
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
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
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
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:664
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3340
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition TH1.cxx:3062
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
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:338
Int_t Fill(Double_t) override
Invalid Fill method.
Definition TH2.cxx:350
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
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:275
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
virtual ULong64_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:404
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
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)
Returns x raised to the power y.
Definition TMath.h:721
constexpr Double_t Pi()
Definition TMath.h:37
Double_t Tan(Double_t)
Returns the tangent of an angle of x radians.
Definition TMath.h:600
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
static void output()