Logo ROOT  
Reference Guide
testUnfold5a.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 <map>
67#include <cmath>
68#include <TMath.h>
69#include <TRandom3.h>
70#include <TFile.h>
71#include <TTree.h>
72
73using namespace std;
74
75TRandom *g_rnd=0;
76
77class ToyEvent {
78public:
79 void GenerateDataEvent(TRandom *rnd);
80 void GenerateSignalEvent(TRandom *rnd);
81 void GenerateBgrEvent(TRandom *rnd);
82 // reconstructed quantities
83 inline Double_t GetPtRec(void) const { return fPtRec; }
84 inline Double_t GetEtaRec(void) const { return fEtaRec; }
85 inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
86 inline Bool_t IsTriggered(void) const { return fIsTriggered; }
87
88 // generator level quantities
89 inline Double_t GetPtGen(void) const {
90 if(IsSignal()) return fPtGen;
91 else return -1.0;
92 }
93 inline Double_t GetEtaGen(void) const {
94 if(IsSignal()) return fEtaGen;
95 else return 999.0;
96 }
97 inline Bool_t IsSignal(void) const { return fIsSignal; }
98protected:
99
100 void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
101 void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
102 void GenerateReco(TRandom *rnd);
103
104 // reconstructed quantities
105 Double_t fPtRec;
106 Double_t fEtaRec;
107 Double_t fDiscriminator;
108 Bool_t fIsTriggered;
109 // generated quantities
110 Double_t fPtGen;
111 Double_t fEtaGen;
112 Bool_t fIsSignal;
113
114 static Double_t kDataSignalFraction;
115
116};
117
118void testUnfold5a()
119{
120 // random generator
121 g_rnd=new TRandom3();
122
123 // data and MC number of events
124 const Int_t neventData = 20000;
125 Double_t const neventSignalMC =2000000;
126 Double_t const neventBgrMC =2000000;
127
128 Float_t etaRec,ptRec,discr,etaGen,ptGen;
129 Int_t istriggered,issignal;
130
131 //==================================================================
132 // Step 1: generate data TTree
133
134 TFile *dataFile=new TFile("testUnfold5_data.root","recreate");
135 TTree *dataTree=new TTree("data","event");
136
137 dataTree->Branch("etarec",&etaRec,"etarec/F");
138 dataTree->Branch("ptrec",&ptRec,"ptrec/F");
139 dataTree->Branch("discr",&discr,"discr/F");
140
141 // for real data, only the triggered events are available
142 dataTree->Branch("istriggered",&istriggered,"istriggered/I");
143 // data truth parameters
144 dataTree->Branch("etagen",&etaGen,"etagen/F");
145 dataTree->Branch("ptgen",&ptGen,"ptgen/F");
146 dataTree->Branch("issignal",&issignal,"issignal/I");
147
148 cout<<"fill data tree\n";
149
150 Int_t nEvent=0,nTriggered=0;
151 while(nTriggered<neventData) {
152 ToyEvent event;
153 event.GenerateDataEvent(g_rnd);
154
155 etaRec=event.GetEtaRec();
156 ptRec=event.GetPtRec();
157 discr=event.GetDiscriminator();
158 istriggered=event.IsTriggered() ? 1 : 0;
159 etaGen=event.GetEtaGen();
160 ptGen=event.GetPtGen();
161 issignal=event.IsSignal() ? 1 : 0;
162
163 dataTree->Fill();
164
165 if(!(nEvent%100000)) cout<<" data event "<<nEvent<<"\n";
166
167 if(istriggered) nTriggered++;
168 nEvent++;
169
170 }
171
172 dataTree->Write();
173 delete dataTree;
174 delete dataFile;
175
176 //==================================================================
177 // Step 2: generate signal TTree
178
179 TFile *signalFile=new TFile("testUnfold5_signal.root","recreate");
180 TTree *signalTree=new TTree("signal","event");
181
182 signalTree->Branch("etarec",&etaRec,"etarec/F");
183 signalTree->Branch("ptrec",&ptRec,"ptrec/F");
184 signalTree->Branch("discr",&discr,"discr/F");
185 signalTree->Branch("istriggered",&istriggered,"istriggered/I");
186 signalTree->Branch("etagen",&etaGen,"etagen/F");
187 signalTree->Branch("ptgen",&ptGen,"ptgen/F");
188
189 cout<<"fill signal tree\n";
190
191 for(int ievent=0;ievent<neventSignalMC;ievent++) {
192 ToyEvent event;
193 event.GenerateSignalEvent(g_rnd);
194
195 etaRec=event.GetEtaRec();
196 ptRec=event.GetPtRec();
197 discr=event.GetDiscriminator();
198 istriggered=event.IsTriggered() ? 1 : 0;
199 etaGen=event.GetEtaGen();
200 ptGen=event.GetPtGen();
201
202 if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
203
204 signalTree->Fill();
205 }
206
207 signalTree->Write();
208 delete signalTree;
209 delete signalFile;
210
211 // ==============================================================
212 // Step 3: generate background MC TTree
213
214 TFile *bgrFile=new TFile("testUnfold5_background.root","recreate");
215 TTree *bgrTree=new TTree("background","event");
216
217 bgrTree->Branch("etarec",&etaRec,"etarec/F");
218 bgrTree->Branch("ptrec",&ptRec,"ptrec/F");
219 bgrTree->Branch("discr",&discr,"discr/F");
220 bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
221
222 cout<<"fill background tree\n";
223
224 for(int ievent=0;ievent<neventBgrMC;ievent++) {
225 ToyEvent event;
226 event.GenerateBgrEvent(g_rnd);
227 etaRec=event.GetEtaRec();
228 ptRec=event.GetPtRec();
229 discr=event.GetDiscriminator();
230 istriggered=event.IsTriggered() ? 1 : 0;
231
232 if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
233
234 bgrTree->Fill();
235 }
236
237 bgrTree->Write();
238 delete bgrTree;
239 delete bgrFile;
240}
241
242Double_t ToyEvent::kDataSignalFraction=0.8;
243
244void ToyEvent::GenerateDataEvent(TRandom *rnd) {
245 fIsSignal=rnd->Uniform()<kDataSignalFraction;
246 if(IsSignal()) {
247 GenerateSignalKinematics(rnd,kTRUE);
248 } else {
249 GenerateBgrKinematics(rnd,kTRUE);
250 }
251 GenerateReco(rnd);
252}
253
254void ToyEvent::GenerateSignalEvent(TRandom *rnd) {
255 fIsSignal=1;
256 GenerateSignalKinematics(rnd,kFALSE);
257 GenerateReco(rnd);
258}
259
260void ToyEvent::GenerateBgrEvent(TRandom *rnd) {
261 fIsSignal=0;
262 GenerateBgrKinematics(rnd,kFALSE);
263 GenerateReco(rnd);
264}
265
266void ToyEvent::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
267 Double_t e_T0=0.5;
268 Double_t e_T0_eta=0.0;
269 Double_t e_n=2.0;
270 Double_t e_n_eta=0.0;
271 Double_t eta_p2=0.0;
272 Double_t etaMax=4.0;
273 if(isData) {
274 e_T0=0.6;
275 e_n=2.5;
276 e_T0_eta=0.05;
277 e_n_eta=-0.05;
278 eta_p2=1.5;
279 }
280 if(eta_p2>0.0) {
281 fEtaGen=TMath::Power(rnd->Uniform(),eta_p2)*etaMax;
282 if(rnd->Uniform()>=0.5) fEtaGen= -fEtaGen;
283 } else {
284 fEtaGen=rnd->Uniform(-etaMax,etaMax);
285 }
286 Double_t n=e_n + e_n_eta*fEtaGen;
287 Double_t T0=e_T0 + e_T0_eta*fEtaGen;
288 fPtGen=(TMath::Power(rnd->Uniform(),1./(1.-n))-1.)*T0;
289 /* static int print=100;
290 if(print) {
291 cout<<fEtaGen
292 <<" "<<fPtGen
293 <<"\n";
294 print--;
295 } */
296}
297
298void ToyEvent::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
299 fPtGen=0.0;
300 fEtaGen=0.0;
301 fPtRec=rnd->Exp(isData ? 2.5 : 2.5);
302 fEtaRec=rnd->Uniform(-3.,3.);
303}
304
305void ToyEvent::GenerateReco(TRandom *rnd) {
306 if(fIsSignal) {
307 Double_t expEta=TMath::Exp(fEtaGen);
308 Double_t eGen=fPtGen*(expEta+1./expEta);
309 Double_t sigmaE=0.1*TMath::Sqrt(eGen)+(0.01+0.002*TMath::Abs(fEtaGen))
310 *eGen;
311 Double_t eRec;
312 do {
313 eRec=rnd->Gaus(eGen,sigmaE);
314 } while(eRec<=0.0);
315 Double_t sigmaEta=0.1+0.02*fEtaGen;
316 fEtaRec=rnd->Gaus(fEtaGen,sigmaEta);
317 fPtRec=eRec/(expEta+1./expEta);
318 do {
319 Double_t tauDiscr=0.08-0.04/(1.+fPtRec/10.0);
320 Double_t sigmaDiscr=0.01;
321 fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
322 } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
323 /* static int print=100;
324 if(print) {
325 cout<<fEtaGen<<" "<<fPtGen
326 <<" -> "<<fEtaRec<<" "<<fPtRec
327 <<"\n";
328 print--;
329 } */
330 } else {
331 do {
332 Double_t tauDiscr=0.15-0.05/(1.+fPtRec/5.0)+0.1*fEtaRec;
333 Double_t sigmaDiscr=0.02+0.01*fEtaRec;
334 fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
335 } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
336 }
337 fIsTriggered=(rnd->Uniform()<1./(TMath::Exp(-fPtRec+3.5)+1.));
338}
const Bool_t kFALSE
Definition: RtypesCore.h:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
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 Double_t Exp(Double_t tau)
Returns an exponential deviate.
Definition: TRandom.cxx:240
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:635
A TTree represents a columnar dataset.
Definition: TTree.h:78
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4524
TBranch * Branch(const char *name, T *obj, Int_t bufsize=32000, Int_t splitlevel=99)
Add a new branch, and infer the data type from the type of obj being passed.
Definition: TTree.h:348
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TTree.cxx:9595
const Int_t n
Definition: legend1.C:16
T etaMax()
Function providing the maximum possible value of pseudorapidity for a non-zero rho,...
Definition: etaMax.h:50
Double_t Exp(Double_t x)
Definition: TMath.h:717
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
Short_t Abs(Short_t d)
Definition: TMathBase.h:120