Logo ROOT  
Reference Guide
testUnfold7a.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 example is documented in conference proceedings:
9///
10/// arXiv:1611.01927
11/// 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
12///
13/// This is an example of unfolding a two-dimensional distribution
14/// also using an auxiliary measurement to constrain some background
15///
16/// The example comprises several macros
17/// - testUnfold7a.C create root files with TTree objects for
18/// signal, background and data
19/// - write files testUnfold7_signal.root
20/// testUnfold7_background.root
21/// testUnfold7_data.root
22///
23/// - testUnfold7b.C loop over trees and fill histograms based on the
24/// TUnfoldBinning objects
25/// - read testUnfold7binning.xml
26/// testUnfold7_signal.root
27/// testUnfold7_background.root
28/// testUnfold7_data.root
29///
30/// - write testUnfold7_histograms.root
31///
32/// - testUnfold7c.C run the unfolding
33/// - read testUnfold7_histograms.root
34/// - write testUnfold7_result.root
35/// testUnfold7_result.ps
36///
37/// \macro_output
38/// \macro_code
39///
40/// **Version 17.6, in parallel to changes in TUnfold**
41///
42/// This file is part of TUnfold.
43///
44/// TUnfold is free software: you can redistribute it and/or modify
45/// it under the terms of the GNU General Public License as published by
46/// the Free Software Foundation, either version 3 of the License, or
47/// (at your option) any later version.
48///
49/// TUnfold is distributed in the hope that it will be useful,
50/// but WITHOUT ANY WARRANTY; without even the implied warranty of
51/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
52/// GNU General Public License for more details.
53///
54/// You should have received a copy of the GNU General Public License
55/// along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
56///
57/// \author Stefan Schmitt DESY, 14.10.2008
58
59#include <iostream>
60#include <map>
61#include <cmath>
62#include <TMath.h>
63#include <TRandom3.h>
64#include <TFile.h>
65#include <TTree.h>
66#include <TLorentzVector.h>
67
68#define MASS1 0.511E-3
69
70using namespace std;
71
72TRandom *g_rnd=0;
73
74class ToyEvent7 {
75public:
76 void GenerateDataEvent(TRandom *rnd);
77 void GenerateSignalEvent(TRandom *rnd);
78 void GenerateBgrEvent(TRandom *rnd);
79 // reconstructed quantities
80 inline Double_t GetMRec(int i) const { return fMRec[i]; }
81 inline Double_t GetPtRec(int i) const { return fPtRec[i]; }
82 inline Double_t GetEtaRec(int i) const { return fEtaRec[i]; }
83 inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
84 inline Double_t GetPhiRec(int i) const { return fPhiRec[i]; }
85 inline Bool_t IsTriggered(void) const { return fIsTriggered; }
86
87 // generator level quantities
88 inline Double_t GetMGen(int i) const {
89 if(IsSignal()) return fMGen[i];
90 else return -1.0;
91 }
92 inline Double_t GetPtGen(int i) const {
93 if(IsSignal()) return fPtGen[i];
94 else return -1.0;
95 }
96 inline Double_t GetEtaGen(int i) const {
97 if(IsSignal()) return fEtaGen[i];
98 else return 999.0;
99 }
100 inline Double_t GetPhiGen(int i) const {
101 if(IsSignal()) return fPhiGen[i];
102 else return 999.0;
103 }
104 inline Bool_t IsSignal(void) const { return fIsSignal; }
105protected:
106
107 void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
108 void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
109 void GenerateReco(TRandom *rnd);
110
111 // reconstructed quantities
112 Double_t fMRec[3];
113 Double_t fPtRec[3];
114 Double_t fEtaRec[3];
115 Double_t fPhiRec[3];
116 Double_t fDiscriminator;
117 Bool_t fIsTriggered;
118 // generated quantities
119 Double_t fMGen[3];
120 Double_t fPtGen[3];
121 Double_t fEtaGen[3];
122 Double_t fPhiGen[3];
123 Bool_t fIsSignal;
124public:
125 static Double_t kDataSignalFraction;
126 static Double_t kMCSignalFraction;
127
128};
129
130void testUnfold7a()
131{
132 // random generator
133 g_rnd=new TRandom3(4711);
134
135 // data and MC number of events
136 Double_t muData0=5000.;
137 // luminosity error
138 Double_t muData=muData0*g_rnd->Gaus(1.0,0.03);
139 // stat error
140 Int_t neventData = g_rnd->Poisson( muData);
141
142 // generated number of MC events
143 Int_t neventSigmc = 250000;
144 Int_t neventBgrmc = 100000;
145
146 Float_t etaRec[3],ptRec[3],phiRec[3],mRec[3],discr;
147 Float_t etaGen[3],ptGen[3],phiGen[3],mGen[3];
148 Float_t weight;
149 Int_t istriggered,issignal;
150
151 //==================================================================
152 // Step 1: generate data TTree
153
154 TFile *dataFile=new TFile("testUnfold7_data.root","recreate");
155 TTree *dataTree=new TTree("data","event");
156
157 dataTree->Branch("etarec",etaRec,"etarec[3]/F");
158 dataTree->Branch("ptrec",ptRec,"ptrec[3]/F");
159 dataTree->Branch("phirec",phiRec,"phirec[3]/F");
160 dataTree->Branch("mrec",mRec,"mrec[3]/F");
161 dataTree->Branch("discr",&discr,"discr/F");
162
163 // for real data, only the triggered events are available
164 dataTree->Branch("istriggered",&istriggered,"istriggered/I");
165 // data truth parameters
166 dataTree->Branch("etagen",etaGen,"etagen[3]/F");
167 dataTree->Branch("ptgen",ptGen,"ptgen[3]/F");
168 dataTree->Branch("phigen",phiGen,"phigen[3]/F");
169 dataTree->Branch("mgen",mGen,"mgen[3]/F");
170 dataTree->Branch("issignal",&issignal,"issignal/I");
171
172 cout<<"fill data tree\n";
173
174 //Int_t nEvent=0,nTriggered=0;
175 for(int ievent=0;ievent<neventData;ievent++) {
176 ToyEvent7 event;
177 event.GenerateDataEvent(g_rnd);
178 for(int i=0;i<3;i++) {
179 etaRec[i]=event.GetEtaRec(i);
180 ptRec[i]=event.GetPtRec(i);
181 phiRec[i]=event.GetPhiRec(i);
182 mRec[i]=event.GetMRec(i);
183 etaGen[i]=event.GetEtaGen(i);
184 ptGen[i]=event.GetPtGen(i);
185 phiGen[i]=event.GetPhiGen(i);
186 mGen[i]=event.GetMGen(i);
187 }
188 discr=event.GetDiscriminator();
189 istriggered=event.IsTriggered() ? 1 : 0;
190 issignal=event.IsSignal() ? 1 : 0;
191
192 dataTree->Fill();
193
194 if(!(ievent%100000)) cout<<" data event "<<ievent<<"\n";
195
196 //if(istriggered) nTriggered++;
197 //nEvent++;
198
199 }
200
201 dataTree->Write();
202 delete dataTree;
203 delete dataFile;
204
205 //==================================================================
206 // Step 2: generate signal TTree
207
208 TFile *signalFile=new TFile("testUnfold7_signal.root","recreate");
209 TTree *signalTree=new TTree("signal","event");
210
211 signalTree->Branch("etarec",etaRec,"etarec[3]/F");
212 signalTree->Branch("ptrec",ptRec,"ptrec[3]/F");
213 signalTree->Branch("phirec",ptRec,"phirec[3]/F");
214 signalTree->Branch("mrec",mRec,"mrec[3]/F");
215 signalTree->Branch("discr",&discr,"discr/F");
216
217 // for real data, only the triggered events are available
218 signalTree->Branch("istriggered",&istriggered,"istriggered/I");
219 // data truth parameters
220 signalTree->Branch("etagen",etaGen,"etagen[3]/F");
221 signalTree->Branch("ptgen",ptGen,"ptgen[3]/F");
222 signalTree->Branch("phigen",phiGen,"phigen[3]/F");
223 signalTree->Branch("weight",&weight,"weight/F");
224 signalTree->Branch("mgen",mGen,"mgen[3]/F");
225
226 cout<<"fill signal tree\n";
227
228 weight=ToyEvent7::kMCSignalFraction*muData0/neventSigmc;
229
230 for(int ievent=0;ievent<neventSigmc;ievent++) {
231 ToyEvent7 event;
232 event.GenerateSignalEvent(g_rnd);
233
234 for(int i=0;i<3;i++) {
235 etaRec[i]=event.GetEtaRec(i);
236 ptRec[i]=event.GetPtRec(i);
237 phiRec[i]=event.GetPhiRec(i);
238 mRec[i]=event.GetMRec(i);
239 etaGen[i]=event.GetEtaGen(i);
240 ptGen[i]=event.GetPtGen(i);
241 phiGen[i]=event.GetPhiGen(i);
242 mGen[i]=event.GetMGen(i);
243 }
244 discr=event.GetDiscriminator();
245 istriggered=event.IsTriggered() ? 1 : 0;
246
247 if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
248
249 signalTree->Fill();
250 }
251
252 signalTree->Write();
253 delete signalTree;
254 delete signalFile;
255
256 // ==============================================================
257 // Step 3: generate background MC TTree
258
259 TFile *bgrFile=new TFile("testUnfold7_background.root","recreate");
260 TTree *bgrTree=new TTree("background","event");
261
262 bgrTree->Branch("etarec",&etaRec,"etarec[3]/F");
263 bgrTree->Branch("ptrec",&ptRec,"ptrec[3]/F");
264 bgrTree->Branch("phirec",&phiRec,"phirec[3]/F");
265 bgrTree->Branch("mrec",&mRec,"mrec[3]/F");
266 bgrTree->Branch("discr",&discr,"discr/F");
267 bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
268 bgrTree->Branch("weight",&weight,"weight/F");
269
270 cout<<"fill background tree\n";
271
272 weight=(1.-ToyEvent7::kMCSignalFraction)*muData0/neventBgrmc;
273
274 for(int ievent=0;ievent<neventBgrmc;ievent++) {
275 ToyEvent7 event;
276 event.GenerateBgrEvent(g_rnd);
277 for(int i=0;i<3;i++) {
278 etaRec[i]=event.GetEtaRec(i);
279 ptRec[i]=event.GetPtRec(i);
280 phiRec[i]=event.GetPhiRec(i);
281 }
282 discr=event.GetDiscriminator();
283 istriggered=event.IsTriggered() ? 1 : 0;
284
285 if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
286
287 bgrTree->Fill();
288 }
289
290 bgrTree->Write();
291 delete bgrTree;
292 delete bgrFile;
293}
294
295Double_t ToyEvent7::kDataSignalFraction=0.75;
296Double_t ToyEvent7::kMCSignalFraction=0.75;
297
298void ToyEvent7::GenerateDataEvent(TRandom *rnd) {
299 fIsSignal=rnd->Uniform()<kDataSignalFraction;
300 if(IsSignal()) {
301 GenerateSignalKinematics(rnd,kTRUE);
302 } else {
303 GenerateBgrKinematics(rnd,kTRUE);
304 }
305 GenerateReco(rnd);
306}
307
308void ToyEvent7::GenerateSignalEvent(TRandom *rnd) {
309 fIsSignal=1;
310 GenerateSignalKinematics(rnd,kFALSE);
311 GenerateReco(rnd);
312}
313
314void ToyEvent7::GenerateBgrEvent(TRandom *rnd) {
315 fIsSignal=0;
316 GenerateBgrKinematics(rnd,kFALSE);
317 GenerateReco(rnd);
318}
319
320void ToyEvent7::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
321
322 // fake decay of Z0 to two fermions
323 double M0=91.1876;
324 double Gamma=2.4952;
325 // generated mass
326 do {
327 fMGen[2]=rnd->BreitWigner(M0,Gamma);
328 } while(fMGen[2]<=0.0);
329
330 double N_ETA=3.0;
331 double MU_PT=5.;
332 double SIGMA_PT=2.0;
333 double DECAY_A=0.2;
334 if(isData) {
335 //N_ETA=2.5;
336 MU_PT=6.;
337 SIGMA_PT=1.8;
338 //DECAY_A=0.5;
339 }
340 fEtaGen[2]=TMath::Power(rnd->Uniform(0,1.5),N_ETA);
341 if(rnd->Uniform(-1.,1.)<0.) fEtaGen[2] *= -1.;
342 fPhiGen[2]=rnd->Uniform(-M_PI,M_PI);
343 do {
344 fPtGen[2]=rnd->Landau(MU_PT,SIGMA_PT);
345 } while((fPtGen[2]<=0.0)||(fPtGen[2]>500.));
346 //========================== decay
348 sum.SetPtEtaPhiM(fPtGen[2],fEtaGen[2],fPhiGen[2],fMGen[2]);
349 // boost into lab-frame
350 TVector3 boost=sum.BoostVector();
351 // decay in rest-frame
352
353 TLorentzVector p[3];
354 double m=MASS1;
355 double costh;
356 do {
357 double r=rnd->Uniform(-1.,1.);
358 costh=r*(1.+DECAY_A*r*r);
359 } while(fabs(costh)>=1.0);
360 double phi=rnd->Uniform(-M_PI,M_PI);
361 double e=0.5*sum.M();
362 double ptot=TMath::Sqrt(e+m)*TMath::Sqrt(e-m);
363 double pz=ptot*costh;
364 double pt=TMath::Sqrt(ptot+pz)*TMath::Sqrt(ptot-pz);
365 double px=pt*cos(phi);
366 double py=pt*sin(phi);
367 p[0].SetXYZT(px,py,pz,e);
368 p[1].SetXYZT(-px,-py,-pz,e);
369 for(int i=0;i<2;i++) {
370 p[i].Boost(boost);
371 }
372 p[2]=p[0]+p[1];
373 for(int i=0;i<3;i++) {
374 fPtGen[i]=p[i].Pt();
375 fEtaGen[i]=p[i].Eta();
376 fPhiGen[i]=p[i].Phi();
377 fMGen[i]=p[i].M();
378 }
379}
380
381void ToyEvent7::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
382 for(int i=0;i<3;i++) {
383 fPtGen[i]=0.0;
384 fEtaGen[i]=0.0;
385 fPhiGen[i]=0.0;
386 }
387 TLorentzVector p[3];
388 for(int i=0;i<2;i++) {
389 p[i].SetPtEtaPhiM(rnd->Exp(15.0),rnd->Uniform(-3.,3.),
390 rnd->Uniform(-M_PI,M_PI),isData ? MASS1 : MASS1);
391 }
392 p[2]=p[0]+p[1];
393 for(int i=0;i<3;i++) {
394 fPtRec[i]=p[i].Pt();
395 fEtaRec[i]=p[i].Eta();
396 fPhiRec[i]=p[i].Phi();
397 fMRec[i]=p[i].M();
398 }
399}
400
401void ToyEvent7::GenerateReco(TRandom *rnd) {
402 if(fIsSignal) {
403 TLorentzVector p[3];
404 for(int i=0;i<2;i++) {
405 Double_t expEta=TMath::Exp(fEtaGen[i]);
406 Double_t coshEta=(expEta+1./expEta);
407 Double_t eGen=fPtGen[i]*coshEta;
408 Double_t sigmaE=
409 0.1*TMath::Sqrt(eGen)
410 +1.0*coshEta
411 +0.01*eGen;
412 Double_t eRec;
413 do {
414 eRec=rnd->Gaus(eGen,sigmaE);
415 } while(eRec<=0.0);
416 Double_t sigmaEta=0.1+0.02*TMath::Abs(fEtaGen[i]);
417 p[i].SetPtEtaPhiM(eRec/(expEta+1./expEta),
418 rnd->Gaus(fEtaGen[i],sigmaEta),
419 remainder(rnd->Gaus(fPhiGen[i],0.03),2.*M_PI),
420 MASS1);
421 }
422 p[2]=p[0]+p[1];
423 for(int i=0;i<3;i++) {
424 fPtRec[i]=p[i].Pt();
425 fEtaRec[i]=p[i].Eta();
426 fPhiRec[i]=p[i].Phi();
427 fMRec[i]=p[i].M();
428 }
429 }
430 if(fIsSignal) {
431 do {
432 Double_t tauDiscr=0.08-0.04/(1.+fPtRec[2]/10.0);
433 Double_t sigmaDiscr=0.01;
434 fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
435 } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
436 } else {
437 do {
438 Double_t tauDiscr=0.15-0.05/(1.+fPtRec[2]/5.0)+0.1*fEtaRec[2];
439 Double_t sigmaDiscr=0.02+0.01*fEtaRec[2];
440 fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
441 } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
442 }
443 fIsTriggered=false;
444 for(int i=0;i<2;i++) {
445 if(rnd->Uniform()<0.92/(TMath::Exp(-(fPtRec[i]-15.5)/2.5)+1.)) fIsTriggered=true;
446 }
447}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define e(i)
Definition: RSha256.hxx:103
#define M_PI
Definition: Rotated.cxx:105
const Bool_t kFALSE
Definition: RtypesCore.h:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
double cos(double)
double sin(double)
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:53
void SetXYZT(Double_t x, Double_t y, Double_t z, Double_t t)
Double_t Pt() const
void SetPtEtaPhiM(Double_t pt, Double_t eta, Double_t phi, Double_t m)
void Boost(Double_t, Double_t, Double_t)
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 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
virtual Double_t Landau(Double_t mean=0, Double_t sigma=1)
Generate a random number following a Landau distribution with location parameter mu and scale paramet...
Definition: TRandom.cxx:369
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma.
Definition: TRandom.cxx:214
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
TVector3 is a general three vector class, which can be used for the description of different vectors ...
Definition: TVector3.h:22
TPaveText * pt
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
Definition: VectorUtil.h:327
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
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
Double_t Gamma(Double_t z)
Computation of gamma(z) for all z.
Definition: TMath.cxx:348
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2275